Proof that any plane perpendicular to x,y plane intersection with the bivariate normal distribution has the shape of a normal distribution. The bivariate normal distribution is given by the equation:
$$f(x,y)=\frac{\exp\left(-\frac{1}{2(1-\rho)^2}\left[\left(\frac{x-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x-\mu_1}{\sigma_1}\right)\left(\frac{y-\mu_2}{\sigma_2}\right)+\left(\frac{y-\mu_2}{\sigma_2}\right)^2\right]\right)}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}$$
for $-\infty<x<\infty$ and $-\infty<y<\infty$, where $\sigma_1>0$, $\sigma_2>0$, and $-1<p<1$.
The normal distribution is given by the equation:
$$f(x)=\frac{e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}}{\sigma\sqrt{2\pi}}$$
for $-\infty<x<\infty$ where $\sigma>0$.



I am trying to show that if I take any plane that is perpendicular to the $x,y$ plane, it's intersection with the bivariate normal distribution  is equivlent to $c\cdot f(x)$, where $f(x)$ is the normal distribution and $c$ is a real constant.
This can also be written as showing the following statement to be true.
$$\frac{\exp\left(-\frac{1}{2(1-\rho)^2}\left[\left(\frac{x-\mu_1}{\sigma_1}\right)^2-2\rho\left(\frac{x-\mu_1}{\sigma_1}\right)\left(\frac{(mx+b)-\mu_2}{\sigma_2}\right)+\left(\frac{(mx+b)-\mu_2}{\sigma_2}\right)^2\right]\right)}{2\pi\sigma_1\sigma_2\sqrt{1-\rho^2}}=c\cdot \frac{e^{-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2}}{\sigma\sqrt{2\pi}}$$
or
$$f(x,mx+b)=c\cdot f(x)$$
Where $m$ and $b$ are any real number.
This does leave out all the planes parallel to the $x$ axis, although those will be trivial to show separately, after showing that $f(y,x)$ is also a bivariate normal distribution.
I tried simplifying this but have made no progress so far, It's quite likely that I am going about this in completely the wrong way.

This is one way of doing it with a geometric flavour, whatever your original Gaussian is you can carry out a linear transformation to a new basis given by the eigenvectors of the covariance matrix so that the elliptical contours of our Gaussian become circles, lets call the new coordinate system $(x,y)$, are of the form
$$
x^2 + y^2 = c,
$$
your original line will also be transformed to a new line, which again to keep the notation simple I will just say is given by $y = mx + d$, now using Lagrangian multiplier or whatever you can find the circle tangent to this line, say the point $\hat{x},\hat{y}$, so we then represent the line in vector form in terms of the single variable
$$
\begin{align}
\ell(t) &= \hat{p} + td \\ 
&=\begin{bmatrix}
\hat{x} \\ \hat{y}
\end{bmatrix} + t
\begin{bmatrix}
1 \\ m
\end{bmatrix},
\end{align}
$$
so defining the function $f(t) = F(x(t),y(t))$ we have
$$
\begin{align}
f^{\prime}(t) &= \nabla F \cdot d \\
&= -\left[ x(t) + m\cdot y(t)\right]f(t)
\end{align}
$$
now I think with a little bit of tidying up you could show this satisfies the same ODE as the univariate Gaussian does, i.e. $\sigma^2 f^{\prime}(z) + (z-\mu)f(z) = 0.$ However going down a slightly different track we see that 
$$
\begin{align}
\frac{d^2}{dt^2} \log f(t) = -1-m^2,
\end{align}
$$
telling us that the log of this function is quadratic in $t$ and in fact that
$$
\begin{align}
f(t) &= \sqrt{2\pi}F(\hat{x},\hat{y})\frac{1}{\sqrt{2\pi}}e^{-\frac{(1+m^2)t^2}{2} } \\
&= \hat{c}\sqrt{2\pi (1+m^2)} \cdot \varphi(t),
\end{align}
$$
where $\hat{c} = F(\hat{x},\hat{y})$, i.e. the value of contour of the circle the line is tangent to, and $\varphi$ is the density of a standard univariate Gaussian. 
To make this fully general one could then transform back to the original plane, but this was the basic strategy


*

*Find the ellipse that the line is tangent to.

*Reparameterise the line in a variable $t$ such that $t=0$ corresponds to the tangent point

*You can then easily perform calculus using the directional derivative to either show $f(t)$ satisfies a certain ODE or probably easier show that the log function is quadratic in $t$ and then form the Laplace approximation

