# Second Harmonic Revisited: An Analytic Quantum Approach
## 2 Results
### Quantum Regime
#### 2.2.6 Multithermal State  
The probability distribution of a multithermal state reads [16]:  
\[p_{\mu}(n)=\binom{n+\mu-1}{n}\frac{1}{(M+1)^{\mu}\left(1/M+1\right)^{n}}, \tag{47}\]  
Figure 2: Coherent input state. Solid blue line: SH autocorrelation function, from Equation (46). Dashed line: SH autocorrelation function in the absence of annihilation processes and superpositions.  
which is the discrete case of \(P(I_{\rm F})\) (see Equation (16)). The mean number \(M=\left[\exp\left(\beta h\omega\right)-1\right]^{-1}\) is the average boson number from the Bose-Einstein statistics [20]. If \(\mu=1\), the well-known thermal distribution:  
\[p(n)=\frac{M^{n}}{(M+1)^{n+1}} \tag{48}\]  
is retrieved. The computation of the mean value and the variance from Equation (47) yields:  
\[N = \mu M \tag{49}\] \[\Delta N^{2} = \mu M(M+1). \tag{50}\]  
We note that the variance can be rewritten in terms of the mean value as \(\Delta N^{2}=N(N/\mu+1)\), so that it is clear that, for large mode numbers, the statistics tends to a Poissonian distribution, as expected (see the final remarks in Section 2.1). A generic input state with multithermal statistics is the mixed state \(\hat{\rho}_{0}=\sum_{n}p_{\mu}(n)|n,0\rangle\langle n,0|\). Here, we write explicitly the mean value and the variance of the output statistics since this case is amenable to a comparison with the classical one:  
\[\langle\hat{n}_{\rm SH}\rangle = \frac{(\mu+1)}{\mu}N^{2}\hat{\gamma}^{2}-\frac{2}{3}\frac{(\mu+1 )}{\mu}N^{2}\left[1+2\frac{\mu+2}{\mu}N\right]\hat{\gamma}^{4}+o(\hat{\gamma} ^{6}) \tag{51}\] \[\langle\Delta\hat{n}_{\rm SH}^{2}\rangle = \langle\hat{n}_{\rm SH}\rangle+2\frac{(2\mu+3)(\mu+1)}{\mu^{3}}N ^{4}\hat{\gamma}^{4}+o(\hat{\gamma}^{6}). \tag{52}\]  
We remark that, at the lowest orders of \(\hat{\gamma}\), \(\langle\hat{n}_{\rm SH}\rangle\) and \(\langle\Delta\hat{n}_{\rm SH}^{2}\rangle\) in Equations (51) and (52) exhibit the same dependence on the number of modes \(\mu\) and on the input mean value \(N\), which was retrieved in Equations (18) and (20) for the classical case.  
Thus, at the first order, we find again the moments of a superthermal distribution. Furthermore, at higher orders, we still have that annihilation processes and superpositions result in a reduction of correlations, as we show in Figure 3. The computation of the SH autocorrelation function yields:  
\[g_{\rm SH}^{(2)}=\frac{(\mu+2)(\mu+3)}{\mu(\mu+1)}\left[1-\frac{4}{3}\left(1+4 \frac{N}{\mu}\right)\hat{\gamma}^{2}\right]+o(\hat{\gamma}^{4}). \tag{53}\]  
In particular, if we compare the second panel of Figure 3 with Figure 2, we find again that if the number of modes is large, we end up with a Poissonian statistics. We remark that the contribution of the only creation processes (dashed red line) is significant to point out the effect of the back-conversion throughout the whole process, but that it has no particular physical meaning itself.  
Figure 3: Chaotic input state. Solid blue line: SH autocorrelation function, from Equation (53). Dashed line: SH autocorrelation function in the absence of annihilation processes and superpositions.