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.
 A: The intersection of the plane with coordinate axes $x$ and $y$
(hereinafter "the $x$-$y$ plane") with a plane that is perpendicular
to the $x$-$y$ plane is a straight line $ax+b =c$ where at least
one of $a$ and $b$ is nonzero. So, what you are
are being asked to show is that when you constrain $x$ and $y$ in
$f_{X,Y}(x,y)$ to satisfy $ax+by = c$, what you get is a multiple
of a normal density function.
So, assuming that $b \neq 0$, replace $y$ by $\frac{c-ax}{b}$ in
$$\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].$$
It should be obvious (but if not, just do a lot of middle-school level
algebra) that what is left is a quadratic function of $x$,
viz. $\alpha x^2 + \beta x + \gamma$. What is not quite so obvious is that $\alpha$ is a positive number. Note that 
\begin{align}
\alpha &= \frac{1}{2(1-\rho^2)}\left[\frac{1}{\sigma_1^2}
+ 2\rho\frac{a}{b\sigma_1\sigma_2} + 
\frac{a^2}{b^2\sigma_2^2}\right]\\
&= \frac{1}{2(1-\rho^2)}\left[\left(\frac{a}{b\sigma_2}+\frac{\rho}{\sigma_1}\right)^2 + \left.\left.\frac{1}{\sigma_1^2}
\right(1-\rho^2\right)\right]\\
&> 0
\end{align}
since $1-\rho^2 > 0$. Consequently, $\alpha x^2 + \beta x + \gamma$
can be massaged (via the technique known as "completing the
square" which is illustrated above) into
$$\frac{1}{2}\left(\frac{x-\mu_3}{\sigma_3}\right)^2 + \delta$$
for $\mu_3$, $\sigma_3$ and $\delta$ which can be written in terms
of $a,b,c$ and the parameters of the bivariate distribution
(the expressions are messy). And so we are done; the bivariate
normal pdf has been shown to be reduced to a univariate normal pdf times
a constant. 
Note that it $b = 0$ (and so $a \neq 0$), we can substitute for
$x$ in terms of $y$ and get a univariate normal pdf with argument
$y$ instead of $x$.
Think of the bivariate normal pdf as defining a solid with base
the $x$-$y$ plane and upper surface $z = f_{X,Y}(x,y)$. Then
what all this is saying is that every cross-section of this
solid (by a plane perpendicular to the $x$-$y$ plane) has the shape
of a normal pdf (except that the "area under the curve" is not
necessarily equal to $1$ as all valid pdfs must have). A
fanciful way, familiar to Americans, of putting it is that 
the solid is a piece of
bologna: no matter how you slice it, it is still bologna!
A: As you can see from another answer, it is possible to carry through
your initial idea to get a proof.
The one detail that I questioned was how we establish that the
coefficient of $x^2$ in
$$
\frac{1}{(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]
$$
is positive (using $m$ and $b$ as defined in the question). In fact, that coefficient is
$$
\alpha = 
\frac{1}{(1-\rho)^2} \left(\frac{1}{\sigma_1^2}-2\rho\frac{m}{\sigma_1\sigma_2}+\frac{m^2}{\sigma_2^2}\right) .
$$
But
$$
\frac{1}{\sigma_1^2}-2\rho\frac{m}{\sigma_1\sigma_2}+\frac{m^2}{\sigma_2^2}
= \left(\frac{1}{\sigma_1}-\rho\frac{m}{\sigma_2}\right)^2 
+  (1 - \rho^2)\frac{m^2}{\sigma_2^2},
$$
and since $-1<\rho<1$ it follows that $1-\rho^2>0$ and
also that $\alpha > 0.$
Note that for any $\rho$ such that $\lvert\rho\rvert>1$ 
it is possible to choose the other parameters so that $\alpha<0,$
and of course the entire polynomial is undefined for $\lvert\rho\rvert=1,$
so the fact that $-1<\rho<1$ is a necessary condition for this
proof to go through.

An alternative is a suitable substitution of variables for both $x$ and $y$ that transform the plane so that the distribution over the new variables is standard normal with covariance zero over the transformed variables and such that the equation of the line in the transformed variables has the form $y'=k.$
Such a substitution certainly exists, but it might involve
translation (to eliminate the means), rotation (to eliminate the covariance), scaling by unequal factors along the transformed
axes (to eliminate the variances), and then a final rotation to
make the perpendicular plane in the question intersect the 
transformed coordinate plane in a "horizontal" line.
This method is conceptually simple, and avoids the concern about signs
that I had with the other approach, but you may find that actually
working through the math of the other approach is easier for you.
A: 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

