What are the four axes on PCA biplot? When you construct a biplot for a PCA analysis, you have principal component PC1 scores on the x-axis and PC2 scores on the y-axis. But what are the other two axes to the right and the top of the screen?
 A: Do you mean, e.g., in the plot that the following command returns?
biplot(prcomp(USArrests, scale = TRUE))


If yes, then the top and the right axes are meant to be used for interpreting the red arrows (points depicting the variables) in the plot.
If you know how the principal component analysis works, and you can read R code, the code below shows you how the results from prcomp() are initially treated by biplot.prcomp() before the final plotting by biplot.default(). These two functions are called in the background when you plot with biplot(), and the following modified code excerpt is from biplot.prcomp().
x<-prcomp(USArrests, scale=TRUE)
choices = 1L:2L
scale = 1
pc.biplot = FALSE
scores<-x$x
lam <- x$sdev[choices]
n <- NROW(scores)
lam <- lam * sqrt(n)
lam <- lam^scale
yy<-t(t(x$rotation[, choices]) * lam)
xx<-t(t(scores[, choices])/lam)
biplot(xx,yy)

Shortly, in the example above, the the matrix of variable loadings (x$rotation) is scaled by the standard deviation of the principal components (x$sdev) times square root of the number of observations. This sets the scale for the top and right axes to what is seen on the plot.
There are other methods to scale the variable loadings, also. These are offered e.g. by the R package vegan.
A: I have a better visualization for the biplot. Please check following figure. 
In the experiment, I am trying to mapping 3d points into 2d (simulated data set). 
The trick to understand biplot in 2d is finding the correct angle to see same thing in 3d. All the data points are numbered, you can see the mapping clearly.

Here is the code to reproduce the results.
require(rgl)
set.seed(0)

feature1=round(rnorm(50)*10+20)
feature2=round(rnorm(50)*10+30)
feature3=round(runif(50)*feature1)

d=data.frame(feature1,feature2,feature3)

head(d)

plot(feature1,feature2)
plot(feature2,feature3)
plot(feature1,feature3)

plot3d(d$feature1, d$feature2, d$feature3, type = 'n')
points3d(d$feature1, d$feature2, d$feature3, color = 'red', size = 10)
shift <- matrix(c(-2, 2, 0), 12, 3, byrow = TRUE)
text3d(d+shift,texts=1:50)
grid3d(c("x", "y", "z"))

pr.out=prcomp(d,scale.=T)
biplot(pr.out)
grid()

