Retrieve Z values from bivariate normal distribution ggplot2

I would like to visualise the distribution of a Z value generated with bivariate normal distribution.

I generated a plot with this code :

library(ggplot2)
library(gridExtra)
library(devtools)

bvn <- data.frame(V1=rnorm(10000, mean=5, sd=2.5), V2=rnorm(10000, mean=1, sd=3),
                  V1B=rnorm(10000, mean=15, sd=2.5), V2B = rnorm(10000, mean=0, sd=3))

htop <- ggplot(data=bvn, aes(x=V1)) + 
  stat_density(colour = "midnightblue", geom="line", size = 1.5, position="identity", show.legend=FALSE) +
  stat_density(aes(x=V2), colour = "orange", geom="line", size = 1.5, position="identity", show.legend=FALSE) +
  scale_x_continuous(" ", limits = c(-10,15))

scatter <- ggplot(data=bvn, aes(x=V1, y=V1B)) + 
  geom_point(size = 0.6) +
  geom_density_2d(colour="midnightblue")+
  geom_point(aes(x=V2, y=V2B),size = 0.6) +
  geom_density_2d(aes(x=V2, y=V2B), colour="orange")+
  scale_x_continuous("X1", limits = c(-10,15)) + 
  scale_y_continuous("X2", limits = c(-15,25))

hright <- ggplot(data=bvn, aes(x=V1B)) + 
  stat_density(colour = "red", geom="line", size = 1, position="identity", show.legend=FALSE) +
  stat_density(aes(x=V2B),colour = "purple", geom="line", size = 1, position="identity", show.legend=FALSE) +
  scale_x_continuous("X2", limits = c(-15,25)) +
  coord_flip()

gridExtra::grid.arrange(htop, scatter, ncol=1, nrow=2, heights=c(1, 3))

I now want to fix a value on one axis (X1) and to retrieve the Z value (materialised with the concentric circles) of both distributions at this X1 value, in order to plot these distributions. Is there a simple way to do this ?