Why averaging in Spectral magnitude domain not in complex domain to estimate spectrum of a process Consider we need the magnitude spectrum of the signal. Signal is recorded in $N$ trials to a certain stimuli. Signal varies from each trial because of noise. So, to estimate the actual magnitude spectrum, there are two approaches:


Avg. in temporal domain: Average out trials in temporal domain to obtain signal of less noise and then take its fourier transform to obtain signal's complex frequency content. Visualize magnitude spectrum. The output would be the same, even if we average the fourier transforms(complex), as it is a linear operator.
Avg. in magnitude spectrum domain: Compute fourier transforms of each trial, get magnitude spectrum and average this across trials.


According to me, approach (1) is correct as we are not neglecting any of it. Averaging has to be done in linear domain, which is real and imaginary parts of the complex output (Fourier transform). Whereas the second approach neglects phase spectrum of each trial, i.e. averaging in Polar representation of complex number individually. Many of my friends, even a professor approve the second approach (2). Their argument is that, because the requirement is the magnitude spectrum which is random, because of noise in each trial, so averaging the required output (Fourier transform) over many trials to obtain its estimate is a reasonable step. I am not convinced with it, as spectrum consists of both phase and magnitude and just neglecting one of it to get an estimate individually, by linear approaches, (here Arithmetic mean) is not proper, as the true underlying magnitude spectrum, getting modulated by noise, may be some form of non-linearity (function mapping temporal signal to magnitude spectrum). Thus, averaging magnitude spectrum out won't null out the noise. So, please help by providing opinions of what may be correct and why ?
You seem to understand the issues at hand.  In your situation, if the trial runs would be identical, except for the noise, averaging the signals (or the complex values) will definitely work, and I recommend that it how you do it.

The notion of averaging magnitudes (or magnitudes squared) across DFTs comes from when you are chunking a signal into multiple DFTs and want to aggregate the results.  In that case, if you average the raw DFTs you will get all sorts of cancellation/reinforcement effects depending on the frequency values of the signal compared to your DFT frame length.  When you average the magnitudes, the phase values are ignored and you will get more desirable results.

As is often the case wiht policy, it sometimes gets followed blindly without the understanding of the rationale behind the policy.  That may be what is behind some of the recommendations you have heard.

If you are trying to decompose your signal (assuming it is periodic) into its constituent harmonics, retaining the phase information is essential.

Hope this helps,

Ced