Hmm, no good deed goes unpunished. ;)
The PD is square law, but it isn't E**2, it's |E**2|, so to be completely safe you have to stay in real-valued functions (sin and cos) and not complex exponentials. (At least at first.)
Assuming that the relative phase of the two beams depends only on time, i.e. that they're both really really single transverse mode, then shining E1(x, y, t) and E2(x, y, t) on a sufficiently large photodiode will get you
i_photo(t) = int(-inf < x < inf) int(-inf < y < inf)[ R|(E1(x,y,t) + E2(x,y,t)|**2]dxdy
= R int int[ |E1|**2 + |E2|**2 + 2 Re{E1 E2*}]dxdy
where R is some constant (R is the responsivity if E**2 has units of watts per square metre, but it doesn't matter here).
With a large enough beat frequency, we can regard E1**2 and E2**2 as DC and filter them out. (Since the negative frequencies are twice as far away as that, we can use the complex exponential notation with no worries.) So at AC, all we get is
i_AC = 2R Re {int int[ E1(x,y,t)E2*(x,y,t)]dxdy }
If both lasers are single transverse mode, this integral is proportional to Re{E1(0,0,t) E2*(0,0,t)}. Since the beat frequency is large, we sample all relative phases of E1 and E2 much more rapidly than their fluctuations, so we pretty much get an honest multiplication. Multiplying them in the time domain is convolving them in the frequency domain, so yes, the frequency variances add, *provided they're computed around the nominal beat frequency*. Of course nobody in his right mind would do anything else, but a computer might. ;)
The only thing that modifies this significantly is the frequency response of the photodiode and the electrical measurement system.
Cheers
Phil Hobbs