Yes, since it can be shown that these products define a basis on the space of functions, so any function can be expressed as a linear combination of them. In particular it's possible to show that these products are eigenfunctions of the spatial differential operator and, for the boundary conditions in question the spectrum is discrete. The slightly tricky part is to prove that the infinite series converge.

That's what the ``sound mathematical proof'' amounts to.

The book "Conduction of Heat in Solids", from Carslaw and Jaeger, shows a simplified proof to why this works, in section 1.15.

Basically, if you assume that the solution for the 1D partial differential equation is known, it can be shown that a function f(x,y,z,t) defined as f(x,y,z,t)=f_x(x,t)*f_y(y,t)*f_z(z,t) will satisfy the corresponding 3D partial differential equation.

This proof is not as "sound" as the explanation given by Stam Nicolis but perhaps will be easier to understand

There's no need to assume that the solution factorizes; what is sufficient to prove is that the eigenfunctions of the spatial part of the differential operator, that appears in the differential equation and satisfy the boundary conditions can be written as products of functions of one variable, which can be immediately shown by direct substitution.

There's nothing ``unsound'' about that.

The more complicated part of the proof involves showing that the eigenfunctions of the differential operator, that satisfy the boundary conditions, define a basis.

One way of showing this is presented here: https://people.math.osu.edu/costin.10/5102/Completeness.pdf

The last part involves showing that the series defined this way converges.

For the diffusion equation this is, in fact, very easy to show, since the terms of the series are of the form Cn*exp(-t*λn)*fn(x,y,z), where Cn are constants, that depend on the boundary conditions, the functions fn are bounded, for finite domains, and the eigenvalues tend to infinity with n, and zero isn't an eigenvalue. (Here n=(nx,ny,nz) ). So the only ``non-trivial'' part is to show that the initial condition can be represented by the series.

These statements define what's called ``the spectral method'', since it relies on the knowledge of the eigenfunctions and eigenvalues of the differential operator that appears in the differential equation.

What's interesting is that, in fact, it's not necessary to use the eigenfunctions of the spatial part of the differential operator itself, that appears in the equation; any set of functions, that satisfy the boundary conditions, and define a basis, can be used in expressing the solution as a linear combination of this basis.

This is known as (the family of) ``the pseudospectral methods'', since one isn't interested in the differential operator whose eigenfunctions these might be. What's more relevant, in practice is that these functions can be evaluated using fast transforms.

Cf. for instance, https://depts.washington.edu/ph506/Boyd.pdf