I did not find a "nice" proof, but the result doesn't seem really nice either, so...
I'd first notice that, using independence and linearity, where is a standard Gaussian r.v..
The direct computation of is not obvious to justify. I think a bit of complex analysis is needed any way.
- One way: first find the Laplace transform (or moment generating function) . Since and , we deduce, for ,
Now we can deduce using analytic continuation. I explain in details. Indeed, for , we know on , and we know a priori that it is analytic on the whole complex plane (using the "Analycity under the theorem"). Therefore: since the analytic function coincides with on , it coincides on and we can take to deduce
We can turn to the variable . Let be fixed. We know for . And we know a priori that it is analytic on . We deduce, like above, the value for where : for any ,
where the complex square root is chosen with a positive real part.
- Another way: through integration in the complex plane. We have hence
In the integral we have a kind of density of a Gaussian r.v. with complex mean and variance and if we apply the above formula we get the result (with a complex square root to be precised...). To make this argument rigorous, you can procede (with some care) to the change of variable that turns the integral on to an integral on a complex line, and finally using Cauchy formula (...), we can change to contour to integrate on again and use the above formula for the Gaussian...
Did you have something else in mind?
It comes up to my mind that one could probably generalize the differential proof of the value of : differentiate with respect to , integrate by parts, to find a first order ODE...
After further investigation : indeed it gives a "a bit" nicer solution. Fix , and consider . We get after integration by part. We deduce . What is ? Taking , we see that . If we know that is and the characteristic function of gamma distribution, then we conclude. Or we do the computation, which again implies some complex analysis.