Schrodinger's Equation in three dimensions

The graphical representation of the probability density distributed over the three-dimensional space would be a four-dimensional plot--just like the plot of a probability density distribution over one dimension is two-dimensional and that of a probability density distribution over two dimensions is three-dimensional. There is no direct way to visualize a four-dimensional plot except via its projections onto lower dimensional spaces.

Now, for the specific theory of decoupled harmonic oscillators, the ground state would be the multiplication of Gaussians in each of the dimensions. Thus, the full probability density distribution over the three-dimensional space, when projected onto one (or two) dimension(s), would simply look like Gaussians in those lower dimensions. But, this simplification is owing to the decoupling of oscillators, such a simplification would not be generically possible.

It is very easy to show that the ground state will be again a Gaussian by using the separation of variables in the Schrödinger equation. I.e. you assume that the wave function can be represented as $$Psi(x,y,z) = psi_1(x)psi_2(y)psi_3(z),$$ and show that you have three independent equations for each of the components in the product. It will work in the same way for any higher number of dimensions.

The easiest way to visualize it is by using two-dimensional contour plots for every pair of dimensions.

By separation of variables the ground state will be the product of $n$ Gaussians in $n$ dimensions. For $n=3$ you’d get begin{align} psi_{000}(x,y,z)sim e^{-(lambda_xx^2+lambda_yy^2+lambda_z z^2)/2} tag{1} end{align} where $lambda_x=sqrt{m_1omega_1/hbar}$ etc. This is again a Gaussian but it’s not a surface in 3D so it’s difficult to visualize. You would need to slice it and plot the cross-section of the slices as surfaces.

In the special case of the isotopic oscillator, where $lambda_x=lambda_y=lambda_z=lambda$, then begin{align} psi_{000}sim e^{-lambda r^2/2}, ,qquad r^2=x^2+y^2+z^2 tag{2} end{align} which is a Gaussian in the radius $r$.

Eqs(1) and (2) automatically generalize in $n$ dimensions, although the n-dimensional Gaussian solid is even more difficult to imagine.

In $N$-dimensions start with the Hamiltonian $$ H = sum_{j=1}^{N} left( frac{p_j^2}{2m_j} + frac{1}{2} m omega_{j} x_j^2 right) , $$ you can then define the ladder operators (setting $hbar to 1$) $$ a_j := sqrt{ frac{m_{j}omega_j}{2} } left( x_{j} + frac{i}{m_{j}omega_j} p_j right) a_j^{dagger} = sqrt{ frac{m_{j}omega_j}{2} } left( x_{j} - frac{i}{m_{j}omega_j} p_j right) $$ and you can directly verify that (using $[x_j, p_{k}] = i delta_{jk}$) $$ H = sum_{j} left( omega_j a_{j}^{dagger} a_{j} + frac{1}{2} right) , $$ which is just a collection of 1D oscillators. You are interested in the ground state $| 0 rangle$ in the position basis, where $$ psi_0(mathbf{x}) := langle mathbf{x} | 0 rangle = langle x_1, x_2, ldots , x_{N} | 0 rangle . $$ As others have suggested, assume that the above function is separable such that $$ psi_{0}(mathbf{x}) = psi^1_0(x_1) psi^2_0(x_2) cdots psi^{N}_0(x_N) . $$ Note that each annihilation operator kills the ground state, with $a_{j} | 0 rangle = 0 $ for each $j$. We can use this to derive a differential equation for each $psi_0^j(x_j)$. Use: $$ langle mathbf{x} | a_{j} | 0 rangle = 0 quad implies quad langle mathbf{x} | left( x_j + frac{i}{m_{j}omega_j} p_j right) | 0 rangle = 0 $$ when you use $langle mathbf{x} | x_j | 0 rangle = x_{j} psi_{0}(mathbf{x})$ and $langle mathbf{x} | p_j | 0 rangle = - i frac{partial}{partial x_j} psi_{0}(mathbf{x})$, and then assume separability of the solution, you eventually derive the DE for each $j$ $$ frac{d}{d x_j} left( logbig[ psi^j_0(x_j) big] right) = - m_{j} omega_{j} x_j , $$ which you can easily integrate to give $$ psi^j_0(x_j) = C_{j} e^{- tfrac{ m_{j} omega_{j}}{2} x_j^2 } , $$ where $C_{j}$ is some constant of integration. Combining the overall solution for the ground state we get $$ psi_0(mathbf{x}) = C expleft( - frac{1}{2} sum_{j=1}^{N} m_{j} omega_{j} x_j^2 right) $$ where we have absorbed all the integration constants into one constant $C = C_1 C_2 cdots C_{N}$. Once you normalize the solution properly the constant is $$ C = left( pi^{-N} prod_{j=1}^{N} m_{j} omega_j right)^{1/4} . $$