format: text

text:
# Quantum properties of transverse pattern formation in second-harmonic generation

## I Introduction

Pattern formation has been an active area of research in many diverse systems . Numerous similarities to pattern formation in other systems have been reported in recent studies in nonlinear optics . However similar, nonlinear optics also displays properties that are wholly unique due to the relevance of quantum aspects in optical systems, one manifestation of this is the inevitable quantum fluctuations of light. In the last decade an effort has been made to study the interplay in the spatial domain between optical pattern formation, known from classical nonlinear optics, and the quantum fluctuations of light . New nonclassical effects such as quantum entanglement and squeezing in patterns were predicted . Another interesting example is the phenomenon of quantum images: below the instability threshold, information about the pattern is encoded in the way the quantum fluctuations of the fields are spatially correlated .

Nonlinear $`\chi ^{(2)}`$-materials immersed in a cavity have shown most promising quantum effects. A paradigm of spatiotemporal quantum behavior has been the optical parametric oscillator (OPO), which despite its striking simplicity is able to display highly complex behavior . In the degenerate OPO, pump photons are down-converted to signal photons at half the frequency and with a high degree of quantum correlation. This might be attributed to the fact that the signal photons are created simultaneously conserving energy and momentum, leading to the notion of twin photons. In the opposite process of second-harmonic generation (SHG) fundamental photons are up-converted to second-harmonic photons at the double frequency. On a classical level, both the OPO and intra-cavity SHG display similar spatiotemporal behavior. The essential difference between them is that in the OPO an oscillation threshold for the process exists, which simultaneously acts as the threshold for pattern formation. On the contrary, SHG always takes place no matter the strength of the pump field, but there is a threshold that marks the onset of pattern formation. This gives pronounced differences with the OPO in the linearized behavior below the threshold for pattern formation. In the OPO the pump and the signal fields effectively decouple and only the latter becomes unstable at threshold. At a microscopic level, the behavior of the OPO close to the threshold can be understood in terms of a unique process in which a pump photon decays into two signal photons with opposite wave numbers. In SHG the fundamental and second-harmonic fields are coupled and both become unstable at threshold. This complicates the picture mainly by the number of microscopic mechanisms that are relevant to describe the pattern formation process. But this complexity, on the other hand, is likely to generate interesting correlations between the fundamental and the second-harmonic field. Recently, transverse quantum properties in the singly resonant SHG setup were investigated . There, squeezing in the fundamental output was observed close to the critical wave number, but since the second-harmonic is not resonated the question of possible correlations between the two fields was not addressed. However, since the second harmonic in the singly resonant case is given directly as a function of the fundamental, correlations similar to the ones observed in the fundamental should be expected. In this paper, we will consider the case of doubly resonant SHG with the aim of investigating the spatial correlations not only within each field (fundamental field and second-harmonic field), but also between the two fields.

For this purpose we use the formalism of quasi-probability distributions . Choosing the use of the Q representation we are able to derive a set of nonlinear Langevin equations that describes the time evolution of the quantum fields in the SHG setup (Sec. II). In Sec. III, the linear stability analysis of this system will be discussed and a proper regime of parameters specified, for which the formalism adopted here is applicable. Section IV will be devoted to an analysis, on a microscopic level, of the implications of the three-wave interactions in the nonlinear crystal. These considerations allow to identify the most important spatial correlations expected in this two-field system and to define suitable quantities to be calculated. In particular, we will focus on equal time correlation functions of intensity fluctuations and we will study photon number variances when looking for nonclassical features of the intra-cavity fields. A systematic study of the spatial correlations is presented first through analytical results in the framework of a linearized theory below the threshold for pattern formation (Sec. V), and also through extensive numerical simulations of the nonlinear Langevin equations reported below (Sec. VI) and above (Sec. VII) the threshold for pattern formation. We conclude in Sec. VIII.

## II Nonlinear quantum model for intra-cavity SHG

We consider a nonlinear $`\chi ^{(2)}`$-material with type I phase matching immersed in a cavity with a high reflection input mirror $`M_1`$ and a fully reflecting mirror $`M_2`$ at the other end, cf. Fig. 1. The cavity is pumped at the frequency $`\omega _1`$ and through the nonlinear interaction in the crystal photons of frequency $`\omega _2=2\omega _1`$ are generated. This is the process of SHG. The cavity supports a discrete number of longitudinal modes, and we will consider the case where only two of these modes are relevant, namely the mode $`\omega _{1,\mathrm{c}av}`$ closest to the fundamental frequency and $`\omega _{2,\mathrm{c}av}`$ closest to the second-harmonic frequency. In the setup shown in Fig. 1 $`\omega _{2,\mathrm{c}av}=2\omega _{1,\mathrm{c}av}`$, but we will allow the cavity resonances to be independent in order to control the detunings individually. The pump beam propagates along the $`z`$-direction and using the mean field approximation, variations in the $`z`$-direction are averaged out. This approach is justified as long as the losses and detunings are small. Due to diffraction the transverse section perpendicular to the $`z`$-direction spanned by the $`xy`$-plane also comes into play. We consider the simple 1D case where only one of the transverse directions is relevant, so variations along the $`y`$-direction are neglected and only the $`x`$-direction is taken into account (this could easily be achieved experimentally by using a crystal with a small height). Let $`\widehat{A}_1(x,t)`$ and $`\widehat{A}_2(x,t)`$ denote the 1D intra-cavity boson operators of the fundamental field (FH) and second-harminic field (SH), respectively. They obey the following equal time commutation relation

$$[\widehat{A}_i(x,t),\widehat{A}_j^{}(x^{},t)]=\delta _{ij}\delta (xx^{}),i,j=1,2.$$
(1)
The Hamiltonian operator describing SHG including diffraction can be written as done in Ref.  for the OPO,

$$\widehat{H}=\widehat{H}_{\mathrm{f}ree}+\widehat{H}_{\mathrm{i}nt}+\widehat{H}_{\mathrm{e}xt},$$
(2)
where the free Hamiltonian is given by

$`\widehat{H}_{\mathrm{f}ree}`$ $`=`$ $`\mathrm{}{\displaystyle 𝑑x\widehat{A}_1^{}(x,t)\left(\delta _1\frac{c^2}{2\omega _1}\frac{^2}{x^2}\right)\widehat{A}_1(x,t)}`$ (4)
$`+\mathrm{}{\displaystyle 𝑑x\widehat{A}_2^{}(x,t)\left(\delta _2\frac{c^2}{4\omega _1}\frac{^2}{x^2}\right)\widehat{A}_2(x,t)}.`$
Here $`\delta _j=\omega _j\omega _{j,\mathrm{c}av}`$ are the detunings from the nearest cavity resonances, $`^2/x^2`$ describes the diffraction, and $`c`$ is the speed of light. The interaction Hamiltonian describes the nonlinear interaction in the material

$$\widehat{H}_{\mathrm{i}nt}=\frac{i\mathrm{}g}{2}dx(\widehat{A}_2(x,t)(\widehat{A}_1^{}(x,t))^2\mathrm{H}.c.),$$
(5)
where $`g`$ is the nonlinear coupling parameter proportional to the $`\chi ^{(2)}`$-nonlinearity of the crystal. The external Hamiltonian describes the effects of the pump injected into the cavity at the fundamental frequency, which is taken to be a classical quantity $`_{\mathrm{i}n}`$, so we have

$$\widehat{H}_{\mathrm{e}xt}=i\mathrm{}𝑑x\left(_{\mathrm{i}n}\widehat{A}_1^{}(x,t)_{\mathrm{i}n}^{}\widehat{A}_1(x,t)\right).$$
(6)
Then the master equation for the density matrix $`\widehat{\rho }`$ in the interaction picture is given by

$$\frac{\widehat{\rho }}{t}=\frac{i}{\mathrm{}}[\widehat{H},\widehat{\rho }]+(\widehat{L}_1+\widehat{L}_2)\widehat{\rho }.$$
(7)
The cavity losses are assumed to occur only through the input coupling mirror to the external continuum of modes, and are here included through the Liouvillian terms

$`\widehat{L}_j\widehat{\rho }=`$ $`{\displaystyle }dx\gamma _j(2\widehat{A}_j(x,t)\widehat{\rho }\widehat{A}_j^{}(x,t)`$ (9)
$`\widehat{\rho }\widehat{A}_j^{}(x,t)\widehat{A}_j\widehat{A}_j^{}(x,t)\widehat{A}_j(x,t)\widehat{\rho }),`$
where $`\gamma _j`$ are the cavity loss rates. Here we have assumed that thermal fluctuations in the system can be neglected.

Using the standard approach of expanding the density matrix into coherent states weighted by a quasi-probability distribution function, the master equation (7) is mapped onto a functional equation, depending on the order for creation and destruction operators . For a Hamiltonian that is quadratic in the field operators this results in a Fokker-Planck equation, implying that the dynamical evolution of the distribution function may be modelled by an equivalent set of classical stochastic Langevin equations. However, due to the contributions of higher order to the Hamiltonian (5) problems may arise. When using the Wigner representation the evolution equation of the quasi-probability functional contains third order derivatives, which means that no equivalent Langevin equations can be found. These third order terms have been shown to model quantum jump processes . In some cases these terms can be neglected leading to an approximate Fokker-Planck form. When using the P or Q representations problems of negative diffusion in the Fokker-Planck equation come into play . To avoid negative diffusion in the P representation, some techniques have been developed where the phase space is doubled , but then numerical problems due to divergent stochastic trajectories generally appear . We choose here to use the Q representation which in a restricted domain of parameters has a nonnegative diffusion matrix and has been shown to be a useful alternative in the similar problem of calculating nonlinear quantum correlations in the OPO . The Q representation has no singularity problems, is bounded and always nonnegative.

Introducing $`\alpha _i`$ and $`\alpha _i^{}`$ as the c-number equivalents of the intra-cavity boson operators $`\widehat{A}_i`$ and $`\widehat{A}_i^{}`$, the evolution equation for the quasi-probability distribution function $`Q(\underset{¯}{\alpha })`$ is

$`{\displaystyle \frac{Q(\underset{¯}{\alpha })}{t}}`$ $`=({\displaystyle \frac{}{\alpha _1}}[(\gamma _1i\delta _1)\alpha _1g\alpha _1^{}\alpha _2i{\displaystyle \frac{c^2}{2\omega _1}}{\displaystyle \frac{^2}{x^2}}_{\mathrm{i}n}]`$ (10)
$`+`$ $`{\displaystyle \frac{}{\alpha _2}}[(\gamma _2i\delta _2)\alpha _2+{\displaystyle \frac{g}{2}}\alpha _1^2i{\displaystyle \frac{c^2}{4\omega _1}}{\displaystyle \frac{^2}{x^2}}]{\displaystyle \frac{g}{2}}\alpha _2{\displaystyle \frac{^2}{\alpha _1^2}}`$ (11)
$`+`$ $`\gamma _1{\displaystyle \frac{^2}{\alpha _1\alpha _1^{}}}+\gamma _2{\displaystyle \frac{^2}{\alpha _2\alpha _2^{}}}+\mathrm{c}.c.)Q(\underset{¯}{\alpha }),`$ (12)
with $`\underset{¯}{\alpha }=\{\alpha _1,\alpha _1^{},\alpha _2,\alpha _2^{}\}`$. This is just an extension to the diffractive case of the result obtained by Savage . Equation (10) has the form of a Fokker-Planck equation, and it has positive diffusion if

$`|\alpha _2|<2{\displaystyle \frac{\gamma _1}{g}}.`$ (13)
As shown below, it is possible to fix the parameters of the system in such a way that the stable solution for the SH field is well below the value $`2\gamma _1/g`$. Fluctuations around this stable solution are small, so that the probability violating the condition (13) is almost zero. Neglecting then stochastic trajectories violating this condition, we may write a set of equivalent Langevin stochastic equations by applying the Ito formalism for the stochastic integration . We then obtain the following nonlinear Langevin equations

$`_t\alpha _1(x,t)=`$ $`(\gamma _1+i\delta _1)\alpha _1(x,t)+g\alpha _1^{}(x,t)\alpha _2(x,t)`$ (16)
$`+i{\displaystyle \frac{c^2}{2\omega _1}}{\displaystyle \frac{^2}{x^2}}\alpha _1(x,t)+_{\mathrm{i}n}+\sqrt{2\gamma _1}\xi _1(x,t),`$
$`_t\alpha _2(x,t)=`$ $`(\gamma _2+i\delta _2)\alpha _2(x,t){\displaystyle \frac{g}{2}}\alpha _1^2(x,t)`$ (18)
$`+i{\displaystyle \frac{c^2}{4\omega _1}}{\displaystyle \frac{^2}{x^2}}\alpha _2(x,t)+\sqrt{2\gamma _2}\xi _2(x,t),`$
with multiplicative Gaussian white noise sources correlated as follows

$`\xi _i^{}(x,t)\xi _j(x^{},t^{})`$ $`=`$ $`\delta _{ij}\delta (xx^{})\delta (tt^{}),`$ (20)
$`\xi _2(x,t)\xi _2(x^{},t^{})`$ $`=`$ $`0,`$ (21)
$`\xi _1(x,t)\xi _1(x^{},t^{})`$ $`=`$ $`{\displaystyle \frac{g\alpha _2(x,t)}{2\gamma _1}}\delta (xx^{})\delta (tt^{}).`$ (22)
We rescale space and time according to

$$\stackrel{~}{t}=t\gamma _1,\stackrel{~}{x}=x/l_d,$$
(23)
where $`l_d`$ is the characteristic length scale given by

$$l_d^2=\frac{c^2}{2\gamma _1\omega _1}.$$
(24)
We also normalize the fields and noise according to

$`A_j(x,t)=\alpha _j(x,t){\displaystyle \frac{g}{\gamma _1}},\stackrel{~}{\xi }_j(x,t)=\xi _j(x,t)\sqrt{{\displaystyle \frac{l_d}{\gamma _1}}},`$ (25)
$`E=_{\mathrm{i}n}{\displaystyle \frac{g}{\gamma _1^2}}.`$ (26)
This allows to rewrite the Langevin equations in dimensionless form:

$`_{\stackrel{~}{t}}A_1(\stackrel{~}{x},\stackrel{~}{t})`$ $`=`$ $`(1+i\mathrm{\Delta }_1)A_1(\stackrel{~}{x},\stackrel{~}{t})+A_1^{}(\stackrel{~}{x},\stackrel{~}{t})A_2(\stackrel{~}{x},\stackrel{~}{t})`$ (29)
$`+i{\displaystyle \frac{^2}{\stackrel{~}{x}^2}}A_1(\stackrel{~}{x},\stackrel{~}{t})+E+\sqrt{{\displaystyle \frac{2}{n_{th}}}}\stackrel{~}{\xi }_1(\stackrel{~}{x},\stackrel{~}{t}),`$
$`_{\stackrel{~}{t}}A_2(\stackrel{~}{x},\stackrel{~}{t})`$ $`=`$ $`(\gamma +i\mathrm{\Delta }_2)A_2(\stackrel{~}{x},\stackrel{~}{t}){\displaystyle \frac{1}{2}}A_1^2(\stackrel{~}{x},\stackrel{~}{t})`$ (31)
$`+{\displaystyle \frac{i}{2}}{\displaystyle \frac{^2}{\stackrel{~}{x}^2}}A_2(\stackrel{~}{x},\stackrel{~}{t})+\sqrt{{\displaystyle \frac{2\gamma }{n_{th}}}}\stackrel{~}{\xi }_2(\stackrel{~}{x},\stackrel{~}{t}),`$
where $`\gamma =\gamma _2/\gamma _1`$ and $`\mathrm{\Delta }_j=\delta _j/\gamma _1`$, and $`E`$ may be taken real. Moreover we have introduced

$$n_{th}=\frac{\gamma _1^2l_d}{g^2},$$
(32)
which in the OPO coincides with the number of photons in the characteristic “area” $`l_d`$ required to trigger the oscillation. The noise strength is seen to scale like $`n_{th}^{1/2}`$. The normalized noise sources are correlated by

$`\stackrel{~}{\xi }_i^{}(\stackrel{~}{x},\stackrel{~}{t})\stackrel{~}{\xi }_j(\stackrel{~}{x}^{},\stackrel{~}{t}^{})`$ $`=`$ $`\delta _{ij}\delta (\stackrel{~}{x}\stackrel{~}{x}^{})\delta (\stackrel{~}{t}\stackrel{~}{t}^{}),`$ (34)
$`\stackrel{~}{\xi }_2(\stackrel{~}{x},\stackrel{~}{t})\stackrel{~}{\xi }_2(\stackrel{~}{x}^{},\stackrel{~}{t}^{})`$ $`=`$ $`0,`$ (35)
$`\stackrel{~}{\xi }_1(\stackrel{~}{x},\stackrel{~}{t})\stackrel{~}{\xi }_1(\stackrel{~}{x}^{},\stackrel{~}{t}^{})`$ $`=`$ $`{\displaystyle \frac{A_2(\stackrel{~}{x},\stackrel{~}{t})}{2}}\delta (\stackrel{~}{x}\stackrel{~}{x}^{})\delta (\stackrel{~}{t}\stackrel{~}{t}^{}).`$ (36)
These noise sources turn out to only to be defined for

$`|A_2(\stackrel{~}{x},\stackrel{~}{t})|<2,`$ (37)
which coincides with the condition (13) for a positive diffusion expressed in terms of the rescaled fields.

In the following the tildes are dropped, and only normalized dimensionless equations are considered. We will also use the terminology $`\omega \omega _1`$ and $`2\omega \omega _2`$.

## III Linearized equations and bifurcation diagram

In this section we consider the linearization of the nonlinear Langevin equations in the Q representation around the homogeneous steady state solutions below the threshold for pattern formation. This approach relies on the assumption that the fluctuations are small with respect to the field mean values, and therefore we expect this approach to break down close to the instability threshold. We will come back later (Sec. VI B) to the question of the validity of the linear approximation. We write the fields as $`A_j(x,t)=𝒜_j+\beta _j(x,t)`$, where $`\beta _j(x,t)`$ represent the fluctuations around $`𝒜_j`$. The classical homogeneous values $`𝒜_j`$ of the fields given by the homogeneous steady state solutions of the deterministic limit ($`n_{th}\mathrm{}`$) of Eqs. (II), as found in Ref.. Using this in Eqs. (II) we find the following set of linearized equations

$`_t\beta _1(x,t)`$ $`=`$ $`(1+i\mathrm{\Delta }_1)\beta _1(x,t)+𝒜_2\beta _1^{}(x,t)`$ (39)
$`+`$ $`𝒜_1^{}\beta _2(x,t)+i{\displaystyle \frac{^2}{x^2}}\beta _1(x,t)+\sqrt{{\displaystyle \frac{2}{n_{th}}}}\xi _1(x,t),`$ (40)
$`_t\beta _2(x,t)`$ $`=`$ $`(\gamma +i\mathrm{\Delta }_2)\beta _2(x,t)𝒜_1\beta _1(x,t)`$ (41)
$`+`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \frac{^2}{x^2}}\beta _2(x,t)+\sqrt{{\displaystyle \frac{2\gamma }{n_{th}}}}\xi _2(x,t).`$ (42)
The correlations of the stochastic sources $`\xi _i(x,t)`$ in the linearized limit become

$`\xi _i^{}(x,t)\xi _j(x^{},t^{})`$ $`=`$ $`\delta _{ij}\delta (xx^{})\delta (tt^{}),`$ (44)
$`\xi _1(x,t)\xi _1(x^{},t^{})`$ $`=`$ $`{\displaystyle \frac{𝒜_2}{2}}\delta (xx^{})\delta (tt^{}),`$ (45)
$`\xi _2(x,t)\xi _2(x^{},t^{})`$ $`=`$ $`0.`$ (46)
With $`𝒜_2`$ being merely a constant, the noise in the linear approximation is not multiplicative any more. However, as in the nonlinear equations we have the restriction

$$|𝒜_2|<2.$$
(47)
We would like to mention that the Wigner representation, in the linear regime, would lead to equivalent results without suffering from any limitation since it satisfies a Fokker-Planck equation for any value of $`|𝒜_2|`$. However, for the sake of a consistent presentation of our results we have chosen to consider the Q representaion also in the linear case.

It is instructive to introduce the spatial Fourier transform of the fluctuations

$$\beta _j(k,t)=_{\mathrm{}}^{\mathrm{}}\frac{dx}{\sqrt{2\pi }}\beta _j(x,t)e^{ikx},$$
(48)
which physically represents the amplitude of the fluctuations in the far field. Considering Eqs. (III) and their complex conjugates, it is readily shown that these amplitudes $`\beta _j(k,t)`$ fulfill a set of equations which can be written in the following matrix form

$`_t\left(\begin{array}{c}\beta _1(k,t)\\ \beta _1^{}(k,t)\\ \beta _2(k,t)\\ \beta _2^{}(k,t)\end{array}\right)=𝐌(k)\left(\begin{array}{c}\beta _1(k,t)\\ \beta _1^{}(k,t)\\ \beta _2(k,t)\\ \beta _2^{}(k,t)\end{array}\right)`$ (58)
$`+\sqrt{{\displaystyle \frac{2}{n_{th}}}}\left(\begin{array}{c}\eta _1(k,t)\\ \eta _1^{}(k,t)\\ \sqrt{\gamma }\eta _2(k,t)\\ \sqrt{\gamma }\eta _2^{}(k,t)\end{array}\right),`$ (63)
$`𝐌(k)=\left(\begin{array}{cccc}\sigma _1(k)& 𝒜_2& 𝒜_1^{}& 0\\ 𝒜_2^{}& \sigma _1^{}(k)& 0& 𝒜_1\\ 𝒜_1& 0& \sigma _2(k)& 0\\ 0& 𝒜_1^{}& 0& \sigma _2^{}(k)\end{array}\right),`$ (68)
where $`\sigma _1(k)=1+i(\mathrm{\Delta }_1k^2)`$ and $`\sigma _2(k)=\gamma +i(\mathrm{\Delta }_2k^2/2)`$ have been introduced and each noise term $`\eta _j(k,t)`$ is the Fourier transform of the noise term appearing in the real space linearized Langevin equations (III). Their correlations are given by

$`\eta _i^{}(k,t)\eta _j(k^{},t^{})`$ $`=`$ $`\delta _{ij}\delta (kk^{})\delta (tt^{}),`$ (70)
$`\eta _1(k,t)\eta _1(k^{},t^{})`$ $`=`$ $`{\displaystyle \frac{𝒜_2}{2}}\delta (k+k^{})\delta (tt^{}),`$ (71)
$`\eta _2(k,t)\eta _2(k^{},t^{})`$ $`=`$ $`0.`$ (72)
The linear stability of the classical equations obtained as the $`n_{th}\mathrm{}`$ limit of Eqs. (III) was investigated by Etrich et al. . A rich variety of instabilities was shown to exist: A self-pulsing instability, that leads to oscillations of the homogeneous steady states without any transverse structure, was present for all parameters. The oscillatory transverse instability leading to patterns traveling in space and time was only present for certain parameters and branched out from the self-pulsing instability. Bistability was demonstrated for large detunings of same sign and for $`\gamma `$ small. Most importantly, for all parameters also stationary transverse instabilities were found to exist, i.e. instabilities at a critical transverse wave number $`k=k_c`$ and with zero imaginary eigenvalue. It was shown that stripe-type solutions exist but are always unstable, and numerical simulations showed that instead hexagons are the dominating stationary transverse instability. The 1D configuration we have chosen to consider here has the advantage that the pattern always will be a stripe and therefore leads to simpler interpretation of the correlations. We will choose a range of parameters in which the stationary transverse instability is accessible as the primary bifurcation. This bifurcation is supercritical in the 1D model.

The choice of parameters must take into account the requirement of applicability of the Q representation. One finds that Eq. (47) can only be satisfied for $`\mathrm{\Delta }_1>0`$ . Using the expressions presented in Ref.  and fixing $`\mathrm{\Delta }_1=2.0`$ and $`\gamma =0.5`$ we obtain the bifurcation diagram shown in Fig. 2 . We observe that for $`\mathrm{\Delta }_2<0`$ it is possible to obtain stationary patterns (solid line) as the primary bifurcation at a critical value of the pump, $`E_t`$; increasing the pump beyond $`E_t`$ eventually the system will also become self-pulsing unstable (dotted line). For $`\mathrm{\Delta }_2>0`$ the transverse oscillatory bifurcation (dashed line) is the primary one, and therefore travelling waves are observed in this region. The bistable area is located for $`\mathrm{\Delta }_2>8.3`$ and hence beyond the range shown here.

Expressing the onset of transverse instability, seen in Fig. 2, in terms of the intra-cavity value of the SH we have the bifurcation diagram for the transverse instability shown in Fig. 3. We see that for $`\mathrm{\Delta }_2<0`$ we are well below the limit for positive diffusion (47). Therefore, the probability of trajectories violating the condition (37) of the nonlinear equations is almost zero. For $`\mathrm{\Delta }_2>0`$, increasing $`\gamma `$ or decreasing $`\mathrm{\Delta }_1`$ towards zero, this threshold gets closer to $`|𝒜_2|=2`$.

We will therefore use the parameters $`\mathrm{\Delta }_1=2.0`$, $`\mathrm{\Delta }_2=2.0`$ and $`\gamma =0.5`$ in the rest of this paper, which gives a pattern formation threshold of $`E_t=7.481757`$ and a critical wave number $`k_c=1.833`$. The noise strength is set to $`n_{th}=10^8`$ which is a typical value for the cavity setup discussed here .

The main task of the following section is to identify the most important correlations we expect to find in the system. For this purpose it is useful to have a good knowledge of the spatial structures that emerge in the system.

Numerical simulations  of the nonlinear Eqs. (II) confirmed the instability at a finite transverse wave number $`k=k_c`$ predicted by the linear stability analysis. Above the threshold for pattern formation modulations was observed around the steady state with wavelengths corresponding to $`k_c`$. This is shown in Fig. 4 where the far field intensity shows distinct peaks at $`k=0`$, corresponding to the homogeneous background, and at $`k=\pm k_c`$ corresponding to the modulations observed in the near field, as well as higher harmonics.

Below threshold the quantum noise will excite the least damped modes and precursors of the spatial pattern are observed. This is shown in Fig. 5 where a space-time plot is presented for the FH near and far field. Clearly a stripe-type pattern is formed, but as time progresses the noise diffuses the pattern  so that averaging over time will wash out this emerging structure and a spatially homogeneous near field will remain. On the contrary, as we will show the spatial correlation functions do encode precise information about the emerging pattern, even after this time averaging has been carried out, as illustrated through the concept of quantum images .

## IV Correlations, photon interaction and pattern formation

Our general objective is the investigation of the spatial intra-cavity field correlations emerging in this system as a result of the coupling of FH and SH field through the nonlinearity of the crystal, and the implications of the spatial instability on these correlations. This study has a two-fold purpose: First, to obtain a precise picture on how pattern formation occurs in cavity SHG. In particular, we will aim at identifying the relevant mechanisms, in terms of elementary three-wave processes that are important for the understanding of the intra-cavity field dynamics. Secondly, it will be interesting to investigate whether these correlations are the manifestation of nonclassical states of the fields. Such states are identified by investigating the statistics of the intra-cavity intensities, looking in particular for possible sub-Poissonian features .

### A Photon interaction

We will start by investigating the equal time correlations between intensity fluctuations at different points in the far field. The intensity of each field being directly proportional to the number of photons in the corresponding mode, we can relate the intensity fluctuations to the creation or destruction of photons. The idea is that the way these fluctuations are correlated gives information about the microscopic mechanisms that take place in the cavity and, ultimately, that are involved in the pattern formation process. Generally speaking, a positive correlation tells us that there should exist a coherent mechanism that creates simultaneously the corresponding photons. The following normalized correlations are considered

$$C_{ij}^n(k,k^{})=\frac{\delta \widehat{N}_i(k,t)\delta \widehat{N}_j(k^{},t)}{\sqrt{\delta \widehat{N}_i(k,t)^2\delta \widehat{N}_j(k^{},t)^2}},$$
(73)
where the superscript $`n`$ denotes normalization. The intensity fluctuations are given by $`\delta \widehat{N}_j(k,t)=\widehat{N}_j(k,t)\widehat{N}_j(k,t)`$, which involves the photon number operator $`\widehat{N}_j(k,t)=\widehat{A}_j^{}(k,t)\widehat{A}_j(k,t)`$. The brackets denote quantum mechanical averages (expectations) of the operators, which in our approach based on numerical simulations of equivalent c-numbers will be translated into an average over time. The normalization of the correlations implies that $`C_{ij}^n(k,k^{})=1`$ for perfectly correlated fluctuations, whereas $`C_{ij}^n(k,k^{})=1`$ will be the signature of perfect anti-correlation between the intensity fluctuations. As usual, the absence of any correlation will translate into a vanishing correlation function $`C_{ij}^n(k,k^{})=0`$. In the following we will refer to $`C_{11}^n(k,k^{})`$ and $`C_{22}^n(k,k^{})`$ as self-correlations (between different modes of a given field) and to $`C_{12}^n(k,k^{})`$ as cross-correlations (between modes in different fields).

As a guideline for the investigation of the properties of these correlation functions, the first step consists in identifying the basic photon processes when the system is taken close to a transverse instability. These photon processes must obey the standard energy and momentum conservation laws. Whereas the former merely implies that each elementary process must connect one SH photon with two FH photons, the latter will translate into a condition on the transverse wave numbers. Keeping in mind that the cavity is pumped with an homogeneous field at the frequency $`\omega `$, the first process to consider consists in two homogeneous FH photons, $`[\omega ](k=0)[\omega ](0)`$ combining to give one homogeneous SH photon, $`[2\omega ](0)`$, what will be written as $`[\omega ](0)+[\omega ](0)[2\omega ](0)`$. This is encoded in the Hamiltonian term $`\widehat{A}_1^2\widehat{A}_2^{}`$ in Eq. (5). The inverse process, which corresponds to the degenerate OPO process, also takes place in the system, as shown by the presence of the term $`(\widehat{A}_1^{})^2\widehat{A}_2`$ in Eq. (5). Elaborating on these considerations we propose the scheme in Fig. 6 as the simplest way of obtaining a pattern in both fields.

* The first step is the basic SHG channel where two homogeneous FH photons give a SH photon and vice versa, i.e. the channel $`[\omega ](0)+[\omega ](0)[2\omega ](0)`$. It is important to realize that fluctuations around the steady state are considered, hence it is not considered how the FH photons combine to give the steady state SH photons via the channel above, but rather how the fluctuations invoke the channel beyond this.
* The second step is the down-conversion of a SH photon into two FH photons. Momentum conservation in the process implies that the two FH photons have the same value of the transverse wave number but with opposite signs. These are called twin photons since an emission of a $`[\omega ](+k^{})`$ photon must be accompanied by an emission of a $`[\omega ](k^{})`$ photon, and they therefore show a high degree of correlation. This channel written as $`[2\omega ](0)[\omega ](k^{})+[\omega ](+k^{})`$ generates off-axis FH photons.
* Off-axis SH photons are obtained by combining the created off-axis FH photon from step 2) with a photon from the homogeneous background to give a SH photon, which by momentum conservation must have the same wave number as the off-axis FH photon. This channel can be written as $`[\omega ](0)+[\omega ](+k^{})[2\omega ](+k^{})`$.

Of course, these are not the only three-wave processes which are kinematically allowed in the nonlinear crystal, since the interaction Hamiltonian (5) induces any process of the form $`[\omega ](k^{})+[\omega ](k^{\prime \prime })[2\omega ](k^{}+k^{\prime \prime })`$, with arbitrary wave numbers $`k^{}`$ and $`k^{\prime \prime }`$. In fact, the basic scheme we propose in Fig. 6 only takes into account those three-wave processes which involve at least one photon of the homogeneous background fields. Empirically, this choice is motivated by the observation that below the threshold these are the only field modes that are macroscopically populated, so that any process involving them should be stimulated in analogy to what occurs in standard stimulated emission. Formally, the selection of these particular elementary processes corresponds precisely to the approximation made by linearizing the field equations around the steady state solution. As can easily be checked, the full equations for the far field fluctuations contain additional terms quadratic in the fluctuation amplitudes, which indeed account for other three-wave processes. Linearizing we are left with Eq. (63), which only take into account the processes represented by step 2) and step 3). These processes translate into nondiagonal elements of the matrix $`𝐌(k)`$ of the linear system, and as a consequence, for any value of $`k`$, the time evolution of the four amplitudes $`\beta _1(k,t)`$, $`\beta _1(k,t)`$, $`\beta _2(k,t)`$ and $`\beta _2(k,t)`$ will be coupled. This coupling is expected to translate into correlations between the intensity fluctuations $`\delta I_1(k)`$, $`\delta I_1(k)`$, $`\delta I_2(k)`$, $`\delta I_2(k)`$.

This preliminary observation already allows to give a more explicit interpretation of the basic scheme of Fig. 6. Splitting the dynamics of the intra-cavity fields into independent elementary steps, as suggested in the discussion of Fig. 6, would not explain any correlations either between $`[\omega ](k^{})`$ and $`[2\omega ](k^{})`$ nor between $`[2\omega ](k^{})`$ and $`[2\omega ](k^{})`$. Hence the inspection of the linearized equations shows that the interpretation of Fig. 6 in terms of a cascade is too naive. Instead, we have to understand step 2) and 3) as two coherent, joint processes, which generate simultaneously correlations between the 4 modes $`[\omega ](k^{})`$, $`[\omega ](k^{})`$, $`[2\omega ](k^{})`$ and $`[2\omega ](k^{})`$. Finally, it is important to stress that the linearized analysis does not predict any correlation between intensity fluctuations in field modes with wave numbers of different modulus. Mathematically, this is due to the fact that in the linear approximation all correlation functions (73) have the structure

$$C_{ij}(k,k^{})=C_{ij}^{()}(k)\delta ^2(kk^{})+C_{ij}^{(+)}(k)\delta ^2(k+k^{}),$$
(74)
as will be shown in the next section. Close enough to threshold, however, this will not be true any more because of the emergence of additional correlations of nonlinear nature.

Let us finally briefly address the fundamental difference between OPO and SHG: Whereas in SHG, the two fields $`𝒜_1`$ and $`𝒜_2`$ are always nonzero regardless of the pump level, in the OPO case below the oscillation threshold $`𝒜_2`$ is fixed by the pump and $`𝒜_1=0`$. Considering the scheme presented in Fig. 6, the vanishing of $`𝒜_1`$ implies that there is no macroscopic population of the mode $`[\omega ](0)`$ and therefore, step 3) of Fig. 6 is not present. The route to pattern formation simply consists of step 2) in Fig. 6, generating correlations between $`\delta \widehat{N}_1(k,t)`$ and $`\delta \widehat{N}_1(k,t)`$. Mathematically, the consequence for the stability of the homogeneous solution is that the two equations (III) effectively decouple and that only the FH becomes unstable at the threshold.

### B Correlations below shot-noise

Once correlations between intensity fluctuations are identified, it is interesting to investigate if they are connected to nonclassical states of the intra-cavity fields. A coherent field obeys Poissonian photon statistics, which implies that the variance and the mean of the photon number operator $`\widehat{N}`$ are equal. Let us consider the photon number operators associated with the sum and difference of the intensities at different far-field points $`\widehat{N}_i(k)\pm \widehat{N}_j(k^{})`$, where $`\widehat{N}_i(k)=\widehat{a}_i^{}(k)\widehat{a}_i(k)`$ and $`\widehat{N}_j(k)=\widehat{a}_j^{}(k)\widehat{a}_j(k)`$ are the number operators of two states $`\widehat{a}_i(k,t)`$ and $`\widehat{a}_j(k,t)`$. Since we will consider equal time correlation functions in the steady state of the system, from now on we will drop the time argument of the field operators. Taking out the special case $`i=j`$ and $`k^{}=k`$ which will be treated separately, the variance expressed in normal order (indicated by dots) reads

$`\mathrm{V}ar(\widehat{N}_i(k)`$ $`\pm \widehat{N}_j(k^{}))=:\mathrm{V}ar(\widehat{N}_i(k)\pm \widehat{N}_j(k^{})):`$ (75)
$`+:\widehat{N}_i(k):[\widehat{a}_i(k),\widehat{a}_i^{}(k)]`$ (76)
$`+:\widehat{N}_j(k^{}):[\widehat{a}_j(k^{}),\widehat{a}_j^{}(k^{})],`$ (77)
where $`\mathrm{V}ar(X)X^2X^2`$. For a coherent state, the normal ordered variance vanishes, and the mean, given by last two terms in (77), represents the shot noise level for the considered quantity

$$S.N.=:\widehat{N}_i(k):[\widehat{a}_i(k),\widehat{a}_i^{}(k)]+:\widehat{N}_j(k^{}):[\widehat{a}_j(k^{}),\widehat{a}_j^{}(k^{})].$$
(78)
If the normal ordered variance becomes negative

$$:\mathrm{V}ar(\widehat{N}_i\pm \widehat{N}_j):<0,$$
(79)
the variance becomes less than the mean, indicating sub-Poissonian behavior. Such a nonclassical state is identified when the correlation normalized to the shot-noise level, defined as

$`C_{ij}^{(\pm )}(k,k^{})`$ (80)
$`{\displaystyle \frac{:\mathrm{V}ar[\widehat{N}_i(k)\pm \widehat{N}_j(k^{})]:}{:\widehat{N}_i(k):[\widehat{a}_i(k),\widehat{a}_i^{}(k)]+:\widehat{N}_j(k^{}):[\widehat{a}_j(k^{}),\widehat{a}_j^{}(k^{})]}}+1,`$ (81)
is such that $`C_{ij}^{(\pm )}(k,k^{})<1`$. The computation of this quantity requires to write the normal ordered quantities appearing in (81) in terms of anti-normal ordered quantities, since these are the quantities that are computed as averages in our Langevin equations associated with the Q representation. Using the identities

$`\mathrm{}\widehat{N}_i(k)\mathrm{}`$ $`=`$ $`\widehat{a}_i(k)\widehat{a}_i^{}(k)=:\widehat{N}_i(k):+[\widehat{a}_i(k),\widehat{a}_i^{}(k)],`$ (83)
$`\mathrm{}\widehat{N}_i^2(k)\mathrm{}`$ $`=`$ $`\widehat{a}_i(k)\widehat{a}_i(k)\widehat{a}_i^{}(k)\widehat{a}_i^{}(k)`$ (84)
$`=`$ $`:\widehat{N}_i(k)^2:+4:\widehat{N}_i(k):[\widehat{a}_i(k),\widehat{a}_i^{}(k)]`$ (86)
$`+2[\widehat{a}_i(k),\widehat{a}_i^{}(k)]^2`$
with three dots indicating anti-normal ordering, Eq. (81) reads, when expressed in terms of anti-normal ordered quantities

$$C_{ij}^{(\pm )}(k,k^{})=\frac{\mathrm{}\mathrm{V}ar[\widehat{N}_i(k)\pm \widehat{N}_j(k^{})]\mathrm{}\mathrm{}\widehat{N}_i(k)\mathrm{}[\widehat{a}_i(k),\widehat{a}_i^{}(k)]\mathrm{}\widehat{N}_j(k^{})\mathrm{}[\widehat{a}_j(k^{}),\widehat{a}_j^{}(k^{})]}{\mathrm{}\widehat{N}_i(k)\mathrm{}[\widehat{a}_i(k),\widehat{a}_i^{}(k)]+\mathrm{}\widehat{N}_j(k^{})\mathrm{}[\widehat{a}_j(k^{}),\widehat{a}_j^{}(k^{})][\widehat{a}_i(k),\widehat{a}_i^{}(k)]^2[\widehat{a}_j(k^{}),\widehat{a}_j^{}(k^{})]^2}.$$
(87)
Then e.g. the normalized correlation $`\mathrm{V}ar[\widehat{N}_1(k)\pm \widehat{N}_1(k)]/S.N.`$ may be found by setting $`i=j=1`$ and $`k^{}=k`$. Equation (87) is valid for $`k,k^{}0`$, while the special case $`k=0`$ will be addressed in the specific cases.

## V Linearized calculations below threshold

Below threshold, the linear approximation scheme allows to derive semi-analytical expressions for the correlation functions defined in the previous section. These may be expressed in terms of the auxiliary correlation function

$`C_{ij}^Q(k,k^{})`$ $`=`$ $`\mathrm{}\delta \widehat{N}_i(k,t)\delta \widehat{N}_j(k^{},t)\mathrm{},i,j=1,2`$ (88)
$`=`$ $`|A_i(k,t)|^2|A_j(k^{},t)|^2`$ (90)
$`|A_i(k,t)|^2|A_j(k^{},t)|^2,`$
where the superscript Q indicates that the average is done with the Q representation, corresponding to anti-normal ordered quantities, as indicated in the first line of (90).

The starting point of our analysis is the set of linearized Langevin equations (63) which have the exact solutions

$`\left(\begin{array}{c}\beta _1(k,t)\\ \beta _1^{}(k,t)\\ \beta _2(k,t)\\ \beta _2^{}(k,t)\end{array}\right)=e^{𝐌(k)t}\left(\begin{array}{c}\beta _1(k,0)\\ \beta _1^{}(k,0)\\ \beta _2(k,0)\\ \beta _2^{}(k,0)\end{array}\right)`$ (99)
$`+\sqrt{{\displaystyle \frac{2}{n_{th}}}}e^{𝐌(k)t}{\displaystyle _0^t}𝑑t^{}e^{𝐌(k)t^{}}\left(\begin{array}{c}\eta _1(k,t^{})\\ \eta _1^{}(k,t^{})\\ \sqrt{\gamma }\eta _2(k,t^{})\\ \sqrt{\gamma }\eta _2^{}(k,t^{})\end{array}\right).`$ (104)
The first term in Eq. (104) describes how the intra-cavity fields with arbitrary initial conditions relax to the steady state solution and it does not contribute to the steady state correlations. The second term in Eq. (104) gives the response of the intra-cavity fields to the vacuum fluctuations entering the cavity through the partially transparent input mirror. Starting from Eq. (104), it is possible to derive semi-analytical expressions for the correlations (90)

$`C_{ij}^Q(k,k^{})`$ $`=`$ $`|\beta _i(k,t)|^2|\beta _j(k^{},t)|^2`$ (105)
$``$ $`|\beta _i(k,t)|^2|\beta _j(k^{},t)|^2`$ (106)
$`+`$ $`2\mathrm{R}e\{𝒜_i^{}𝒜_j\beta _i(k,t)\beta _j^{}(k^{},t)`$ (108)
$`+𝒜_i^{}𝒜_j^{}\beta _i(k,t)\beta _j(k^{},t)\}\delta (k)\delta (k^{}),`$
where $`\mathrm{R}e\{\}`$ denotes the real part. Whereas the first two terms in the right-hand side (r.h.s.) of Eq. (108) measure the correlations in the intensities of the fluctuations, the last two terms can be traced back to interferences between the fluctuations and the homogeneous component of each field. Since these interferences only contribute to the equal time correlations when $`k=k^{}=0`$, we will first concentrate on $`k,k^{}0`$ and come back later to this special case. Henceforth, unless otherwise specified we consider the case $`k,k^{}0`$.

The Gaussian character of the fluctuations in this linearized Langevin model allows to factorize (108) in terms of second order moments of the field fluctuations

$`C_{ij}^Q(k,k^{})=|\beta _i(k,t)\beta _j^{}(k^{},t)|^2+|\beta _i(k,t)\beta _j(k^{},t)|^2.`$ (109)
The field correlations $`\beta _i(k,t)\beta _j^{}(k^{},t)`$ and $`\beta _i(k,t)\beta _j(k^{},t)`$ can be best evaluated for the solution Eq. (104) if we introduce the set of eigenvectors $`\{𝐯^{(l)}(k)\}_{l=1,\mathrm{},4}`$ of the matrix $`𝐌(k)`$, defined through:

$$𝐌(k)𝐯^{(l)}(k)=\lambda ^{(l)}(k)𝐯^{(l)}(k).$$
(110)
An arbitrary $`4`$-component vector $`𝐰`$ can be decomposed on this basis

$$𝐰(k)=\left(\begin{array}{c}w_1(k)\\ w_2(k)\\ w_3(k)\\ w_4(k)\end{array}\right)=\underset{l=1}{\overset{4}{}}w^{(l)}(k)𝐯^{(l)}(k),$$
(111)
and its components $`w^{(l)}`$ in the new basis are calculated via the linear transformation

$$w^{(l)}(k)=\underset{m=1}{\overset{4}{}}T_{lm}(k)w_m(k).$$
(112)
This involves a $`4\times 4`$ matrix $`T_{lm}(k)`$ calculated as $`𝐓(k)=𝐕(k)^1`$ with $`V_{lm}(k)=v_l^{(m)}(k)`$. Decomposing now the noise vector appearing on the r.h.s. of Eq. (104) on this basis

$$\left(\begin{array}{c}\eta _1(k,t^{})\\ \eta _1^{}(k,t^{})\\ \sqrt{\gamma }\eta _2(k,t^{})\\ \sqrt{\gamma }\eta _2^{}(k,t^{})\end{array}\right)=\underset{l=1}{\overset{4}{}}\eta ^{(l)}(k,t^{})𝐯^{(l)}(k),$$
(113)
allows to rewrite Eq. (104) in the large time limit as

$`\left(\begin{array}{c}\beta _1(k,t)\\ \beta _1^{}(k,t)\\ \beta _2(k,t)\\ \beta _2^{}(k,t)\end{array}\right)`$ $`=\sqrt{{\displaystyle \frac{2}{n_{th}}}}{\displaystyle _0^t}𝑑t^{}`$ (119)
$`\times {\displaystyle \underset{l=1}{\overset{4}{}}}e^{\lambda ^{(l)}(k)(tt^{})}\eta ^{(l)}(k,t^{})𝐯^{(l)}(k).`$
The needed field correlations are given as

$`\beta _i(k,t)\beta _j^{}(k^{},t)`$ $`=`$ $`{\displaystyle \frac{2}{n_{th}}}{\displaystyle _0^t}𝑑t^{}{\displaystyle _0^t}𝑑t^{\prime \prime }{\displaystyle \underset{l,m=1}{\overset{4}{}}}v_{2i1}^{(l)}(k)v_{2j1}^{(m)}{}_{}{}^{}(k^{})e^{\lambda ^{(l)}(k)(tt^{})}e^{\lambda _{}^{(m)}{}_{}{}^{}(k^{})(tt^{\prime \prime })}\eta ^{(l)}(k,t^{})\eta _{}^{(m)}{}_{}{}^{}(k^{},t^{\prime \prime }),`$ (121)
$`\beta _i(k,t)\beta _j(k^{},t)`$ $`=`$ $`{\displaystyle \frac{2}{n_{th}}}{\displaystyle _0^t}𝑑t^{}{\displaystyle _0^t}𝑑t^{\prime \prime }{\displaystyle \underset{l,m=1}{\overset{4}{}}}v_{2i1}^{(l)}(k)v_{2j1}^{(m)}(k^{})e^{\lambda ^{(l)}(k)(tt^{})}e^{\lambda ^{(m)}(k^{})(tt^{\prime \prime })}\eta ^{(l)}(k,t^{})\eta ^{(m)}(k^{},t^{\prime \prime }).`$ (122)
The noise correlations in the new basis $`\eta ^{(l)}(k,t^{})\eta _{}^{(m)}{}_{}{}^{}(k^{},t^{\prime \prime })`$ and $`\eta ^{(l)}(k,t^{})\eta ^{(m)}(k^{},t^{\prime \prime })`$ are

$`\eta ^{(l)}(k,t^{})\eta _{}^{(m)}{}_{}{}^{}(k^{},t^{\prime \prime })`$ $`=`$ $`A_{lm}(k)\delta (kk^{})\delta (t^{}t^{\prime \prime }),`$ (124)
$`\eta ^{(l)}(k,t^{})\eta ^{(m)}(k^{},t^{\prime \prime })`$ $`=`$ $`B_{lm}(k)\delta (k+k^{})\delta (t^{}t^{\prime \prime }).`$ (125)
where the matrix elements of the $`4\times 4`$ matrices $`𝐀(k)`$ and $`𝐁(k)`$ can easily be evaluated in terms of the matrix elements $`T_{lm}T_{lm}(k)`$ as

$`A_{lm}(k)=`$ $`T_{l1}T_{m1}^{}{\displaystyle \frac{𝒜_2}{2}}T_{l1}T_{m2}^{}+T_{l2}T_{m2}^{}`$ (128)
$`{\displaystyle \frac{𝒜_2^{}}{2}}T_{l2}T_{m1}^{}+\gamma T_{l3}T_{m3}^{}+\gamma T_{l4}T_{m4}^{},`$
$`B_{lm}(k)=`$ $`T_{l1}T_{m2}{\displaystyle \frac{𝒜_2}{2}}T_{l1}T_{m1}+T_{l2}T_{m1}`$ (130)
$`{\displaystyle \frac{𝒜_2^{}}{2}}T_{l2}T_{m2}+\gamma T_{l3}T_{m4}+\gamma T_{l4}T_{m3}.`$
Inserting Eqs. (V) in Eqs. (V) we can easily carry out the time integration, and neglecting transient contributions, we end up with the following expressions

$`\underset{t\mathrm{}}{lim}\beta _i(k,t)\beta _j^{}(k^{},t)`$ $`=`$ $`{\displaystyle \frac{2}{n_{th}}}G_{ij}^{()}(k)\delta (kk^{}),`$ (132)
$`\underset{t\mathrm{}}{lim}\beta _i(k,t)\beta _j(k^{},t)`$ $`=`$ $`{\displaystyle \frac{2}{n_{th}}}G_{ij}^{(+)}(k)\delta (k+k^{}).`$ (133)
with

$`G_{ij}^{()}(k)={\displaystyle \underset{l=1}{\overset{4}{}}}{\displaystyle \underset{m=1}{\overset{4}{}}}A_{lm}(k){\displaystyle \frac{v_{2i1}^{(l)}(k)v_{2j1}^{(m)}{}_{}{}^{}(k)}{(\lambda ^{(l)}(k)+\lambda _{}^{(m)}{}_{}{}^{}(k))}},`$ (135)
$`G_{ij}^{(+)}(k)={\displaystyle \underset{l=1}{\overset{4}{}}}{\displaystyle \underset{m=1}{\overset{4}{}}}B_{lm}(k){\displaystyle \frac{v_{2i1}^{(l)}(k)v_{2j1}^{(m)}(k)}{(\lambda ^{(l)}(k)+\lambda ^{(m)}(k))}}.`$ (136)
In terms of $`G_{ij}^{()}(k)`$ and $`G_{ij}^{(+)}(k)`$, Eq. (109) is given by

$`C_{ij}^Q(k,k^{})=`$ $`{\displaystyle \frac{4}{n_{th}^2}}[|G_{ij}^{()}(k)|^2\delta ^2(kk^{})`$ (138)
$`+|G_{ij}^{(+)}(k)|^2\delta ^2(k+k^{})].`$
### A Intensity fluctuation correlations

It is now easy to compute the normalized correlation function Eq. (73). This involves taking into account the commutation relation Eq. (1) which reads

$`[\widehat{A}_i(k,t),\widehat{A}_{j}^{}{}_{}{}^{}(k^{},t)]=\delta _{ij}{\displaystyle \frac{1}{n_{th}}}\delta (kk^{}),`$ (139)
after rescaling space and time according to Eq. (23) and the operators similar to the c-number fields in Eq. (25). We finally find

$`C_{ij}^n(k,k^{})=`$ $`{\displaystyle \frac{|G_{ij}^{()}(k)|^2}{\sqrt{\eta _i(k)}\sqrt{\eta _j(k)}}}{\displaystyle \frac{\delta (kk^{})^2}{\delta (0)^2}}`$ (141)
$`+{\displaystyle \frac{|G_{ij}^{(+)}(k)|^2}{\sqrt{\eta _i(k)}\sqrt{\eta _j(k)}}}{\displaystyle \frac{\delta (k+k^{})^2}{\delta (0)^2}},`$
with $`\eta _j(k)=G_{jj}^{()}(k)(G_{jj}^{()}(k)1/2)`$, the $`1/2`$ in the parenthesis reflecting the conversion from anti-normal to direct ordering. Unlike the mathematical expression (141) derived for an ideally infinite system, the correlation functions determined from the simulations will have peaks of a finite width, which will be determined by the discretization in $`k`$-space used in the numerical codes, i.e. the inverse of the total length of the system. This difference, however, will not alter the only relevant information, which is the height of each of these peaks. In fact, the quantities

$`C_{jj}^n(k,k)`$ $`=`$ $`{\displaystyle \frac{|G_{jj}^{(+)}(k)|^2}{\eta _j(k)}}`$ (143)
$`C_{12}^n(k,\pm k)`$ $`=`$ $`{\displaystyle \frac{|G_{12}^{()}(k)|^2}{\sqrt{\eta _1(k)}\sqrt{\eta _2(k)}}},`$ (144)
characterize the strength of the correlations between the modes $`[\omega ](k)`$ and $`[\omega ](k)`$, $`[2\omega ](k)`$ and $`[2\omega ](k)`$, $`[\omega ](k)`$ and $`[2\omega ](k)`$, and $`[\omega ](k)`$ and $`[2\omega ](k)`$ respectively. One easily checks that $`C_{ii}^n(k,k)=1`$, as a result of an autocorrelation.

All the expressions derived so far are only valid for nonvanishing transverse wave numbers. At $`k=k^{}=0`$, we already observed that there are extra contributions to the equal time correlation function, as expressed by Eq. (108). Furthermore, in the framework of an expansion in the small parameter $`\sqrt{2/n_{th}}`$, it is obvious that these extra terms even dominate, since they scale with $`|\beta _i(k,t)|^22/n_{th}`$, whereas the contributions on the first line of Eq. (108) scales with $`|\beta _i(k,t)|^4(2/n_{th})^2`$. Hence, in the leading order, the correlation function at $`k=k^{}=0`$ is given by

$`C_{12}^Q(k,k^{})`$ $`|_{k=k^{}=0}={\displaystyle \frac{2\delta (0)}{n_{th}}}2\mathrm{R}e(𝒜_1^{}𝒜_2G_{12}^{()}(0)`$ (146)
$`+𝒜_1^{}𝒜_2^{}G_{12}^{(+)}(0))\delta (k)\delta (k^{})|_{k=k^{}=0}.`$
Similar calculations as before allow us to derive the following expression for the value of the normalized cross-correlation at $`k=k^{}=0`$

$$C_{12}^n(0,0)=\frac{\mathrm{R}e\left(𝒜_1^{}𝒜_2G_{12}^{()}(0)+𝒜_1^{}𝒜_2^{}G_{12}^{(+)}(0)\right)}{\sqrt{\zeta _1}\sqrt{\zeta _2}},$$
(147)
where $`\zeta _j=|𝒜_j|^2(G_{jj}^{()}(0)1/4)+\mathrm{R}e\{𝒜_{j}^{}{}_{}{}^{2}G_{jj}^{(+)}(0)\}`$.

### B Nonclassical photon number variances

The photon number variances considered in Sec. IV B can be calculated in terms of the auxiliary functions $`G_{ij}^{()}(k)`$ and $`G_{ij}^{(+)}(k)`$ as well. The anti-normal ordered quantities in Eq. (87) can be directly calculated by averages in the Langevin equation, so below threshold the anti-normal ordered variance is for $`k0`$

$$\mathrm{}\mathrm{V}ar[\widehat{N}_i(k)\pm \widehat{N}_j(k)]\mathrm{}=\mathrm{V}ar[|\beta _i(k,t)|^2\pm |\beta _j(k,t)|^2].$$
(148)
Using the commutation relations (139), the commutators in Eq. (87) are $`[a_j(k),a_j^{}(k)]=\delta (0)/n_{th}`$, and the normalized self-correlations takes the form

$$C_{jj}^{(\pm )}(k,k)=\frac{2\left(|G_{jj}^{()}(k)|^2\pm |G_{jj}^{(+)}(k)|^2\right)G_{jj}^{()}(k)}{G_{jj}^{()}(k)1/2}.$$
(149)
Similarly, the cross-correlations are

$$C_{12}^{(\pm )}(k,\nu k)=\frac{2\left(_j|G_{jj}^{()}(k)|^2\pm 2|G_{12}^{(\nu )}(k)|^2\right)_jG_{jj}^{()}(k)}{_jG_{jj}^{()}(k)1},\nu =+1,1$$
(150)
When $`k=k^{}=0`$ Eq. (148) is no longer valid. Instead, following the procedure outlined for the normalized correlations we have to the leading order $`O(n_{th}^1)`$

$$C_{12}^{(\pm )}(0,0)=4\frac{\mathrm{R}e\left[_j𝒜_{j}^{}{}_{}{}^{2}G_{jj}^{(+)}(0)\pm 2𝒜_1^{}\left(𝒜_2^{}G_{12}^{(+)}(0)+𝒜_2G_{12}^{()}(0)\right)\right]+_j|𝒜_j|^2G_{jj}^{()}(0)}{_j|𝒜_j|^2}1.$$
(151)
The self-correlations become

$`C_{jj}^{()}(0,0)`$ $`=`$ $`0,`$ (153)
$`C_{jj}^{(+)}(0,0)`$ $`=`$ $`4\mathrm{R}e\left[e^{i2\varphi _{A_j}}G_{jj}^{(+)}(0)\right]+4G_{jj}^{()}(0)1.`$ (154)
where $`\varphi _{A_j}`$ is the phase of $`𝒜_j`$. Note that $`C_{jj}^{(+)}(0,0)`$ is actually $`\mathrm{V}ar[\widehat{N}_j(0)]`$ normalized to shot-noise. The result of Eq. (153) is simply because the correlation $`C_{jj}^{()}(k,k^{})`$ amounts to calculating the variance of zero for $`k=k^{}=0`$.

## VI Correlations below threshold

The linearized results of Sec. V gives an analytical insight to the behaviour below threshold for pattern formation. However, very close to the threshold this linear approximation breaks down because of critical nonlinear fluctuations, and additional contributions may emerge as for example shown in a vector Kerr model by Hoyuelos et al. . Such nonlinear correlations can be calculated through numerical simulations of the full nonlinear evolution equations.

In this section we present numerical results obtained from simulations of the nonlinear equations (II) below threshold, with the parameters discussed in Sec. III. Our numerical results are compared with the analytical results of the previous section, and therefore also serves as a cross-check of our analytical and numerical methods.

### A Linear correlations: Analytical and numerical results

We first consider the strength of the correlations between symmetric points in the far fields below the threshold for pattern formation. In Fig. 7, the four quantities defined by Eq. (73) are plotted. The data are obtained from numerical simulations and from the analytical results of Eqs. (V A) and (147). Very good agreement is found between numerics and analytical results.

There are three main features to be considered in the results of Fig. 7. First, all curves present a distinctly peaked behavior around the critical wave number $`k_c`$ for pattern formation, which means that the corresponding modes are more strongly correlated than the modes at any other wave number. Manifestly, this behavior is connected with the pattern formation mechanism and is closely related to the phenomenon of quantum images . Secondly, we also note that in all four plots the correlations show a jump at $`k=0`$. In a) and b) it is the trivial manifestation of an auto-correlation, since for $`k=0`$, $`k`$ and $`k`$ coincide, while in c) and d) the jump is due to the extra interferences with the homogeneous background fields as predicted from Eq. (147). Finally, we observe that the peaks localized around $`k_c`$ are superimposed onto smooth correlation profiles.

The strong correlations appearing between the modes associated with wave numbers around $`k_c`$ indicate a strongly synchronized emission of photons in the modes $`[\omega ](+k)`$, $`[\omega ](k)`$ and $`[2\omega ](+k)`$, $`[2\omega ](k)`$. This behavior reflects the direction of instability of the system. As a matter-of-fact, regardless that all transverse modes of both fields are equally excited by the vacuum fluctuations entering the cavity, the fluctuations of the intra-cavity field modes around the critical wave vector will be less damped than the fluctuations in the other modes. The closer to the threshold, the more the behavior of the intra-cavity fields will be dominated by the mode that becomes unstable at the threshold and gives rise to the pattern. In the 4-dimensional phase space spanned by the fluctuation amplitudes $`\{\beta _1(k,t),\beta _1^{}(k,t),\beta _2(k,t),\beta _2^{}(k,t)\}`$, this mode is characterized by a vector with a given direction. What we learn from the correlation functions is that the emerging instability results in an almost perfectly synchronized emission of photons in the modes $`[\omega ](+k)`$, $`[\omega ](k)`$ and $`[2\omega ](+k)`$, $`[2\omega ](k)`$.

The dominance of this particular mode when the threshold is approached is confirmed by the study of the strength of these correlations as a function of the pump. In Fig. 8 we follow the height of the peaks at $`k=k_c`$ of the four linear correlations displayed in Fig. 7, as a function of the pump level $`E/E_t`$. The most immediate observation is that all the correlations become perfect in the limit $`EE_t`$. This asymptotic behavior can be understood from the linearized fluctuation analysis presented in Sec. V. It is enough to observe that Eqs. (V) involve the inverse of the real part of the eigenvalues of the linear system (63). The dominance at the threshold of the undamped eigenmode of the linear system (63) emerges from the fact that here the real part of the associated eigenvalue precisely goes to zero. Thus, the decrease in the correlations as we move away from threshold can be seen as the result of the coexistence of different eigenmodes. Physically the emergence of these correlations is much less intuitive than the ones in an OPO. As a matter-of-fact, in the OPO below the threshold momentum conservation is enough to predict the existence of correlations between the fluctuations in the modes $`[\omega ](+k)`$ and $`[\omega ](k)`$. In the presence of the 4-mode interaction of SHG, the momentum conservation gives a global condition involving all four beams (at $`[\omega ](+k)`$, $`[\omega ](k)`$ and $`[2\omega ](+k)`$, $`[2\omega ](k)`$). These correlations in fact arise in connection with the emergence of an instability.

Turning now to the cross-correlation between the homogeneous components of the fields, we observe that $`C_{12}^n(k=0,k^{}=0)`$ in Fig. 7 is negative, reflecting an anticorrelation of the photons associated with the FH and SH homogeneous waves. In other words, the creation of a photon $`[2\omega ](0)`$ implies the destruction of (two) photons $`[\omega ](0)`$ and vice versa. The origin of this correlation is much simpler to understand than the previous one: The two modes $`[\omega ](0)`$ and $`[2\omega ](0)`$ being macroscopically populated, the vacuum fluctuations simply induce transitions between these two modes, according to step 1) in the scheme in Fig. 6. In Fig. 9 we plot this correlation as a function of the pump. Comparing the value of the correlations below and above threshold, we observe that very close to, but below, the threshold, the tendency of the curve is reversed and it anticipates the behavior of the correlation above threshold. These are nonlinear correlation effects that will be discussed in Sec. VI B.

Finally, we would like to discuss the smooth contributions to the correlations displayed in Fig. 7. We first note that these are not connected with the pattern instability. This was checked by considering very low pump values for which the peaks around $`k_c`$ completely vanish, while the smooth structures of the curves remain. Considering the central region of the curves, roughly for $`|k|<k_c`$, the most striking observation is the absence of correlations between the fluctuations in the modes $`[\omega ](k)`$ and $`[2\omega ](k)`$, whereas $`[\omega ](k)`$ and $`[2\omega ](k)`$ are correlated, as well as $`[2\omega ](k)`$ with $`[2\omega ](k)`$. This behavior seems to indicate the existence of a symmetry restoring principle in the dynamics of the intra-cavity fields. As a matter-of-fact, the absence of correlations between $`[\omega ](k)`$ and $`[2\omega ](k)`$ implies that the fluctuations of the numbers of pair productions through step 2) and the fluctuations of the number of conversions $`[\omega ](k)[2\omega ](k)`$ through step 3) occur independently of each other. However, while step 2) of Fig. 6 conserves the $`kk`$ symmetry of the system, step 3) does not. As a consequence, a positive fluctuation in the number of times step 3) occurs ($`[\omega ](k)+[\omega ](0)[2\omega ](k)`$), automatically implies that there will be less $`[\omega ](k)`$ than $`[\omega ](k)`$ in the system, and more $`[2\omega ](k)`$ than $`[2\omega ](k)`$. The correlations observed may indicate that the system will try to restore the $`kk`$ symmetry by down-converting $`[2\omega ](0)[\omega ](k)+[\omega ](k)`$, producing a surplus of $`[\omega ](k)`$ which again will produce more $`[2\omega ](k)`$. These mechanisms seems to fit well with the relative strengths of the correlations observed in the central region of Fig. 7. The strongest is always $`C_{11}^n(k,k)`$, in agreement with the fact that the twin photon emission is the principal source of correlations in the system. Weaker is the correlation $`C_{12}^n(k,k)`$ and even weaker $`C_{22}^n(k,k)`$. This interpretation is consistent with the way the correlations at $`k=k_c`$ depart from the value $`1`$ at threshold, when the pump is lowered, as displayed in Fig. 8.

We now turn our attention to the study of the fluctuations in the sum and difference of the photon numbers at symmetrical points of the far field. We first consider the twin beam photon variances for the FH, $`C_{11}^{(\pm )}(k,k)`$ defined in Eq. (81) and shown in Fig. 10. The results are symmetric with respect to the substitution $`kk`$, wherefore we plotted this quantity for positive $`k`$, shifting the origin for better view of the specific behavior at $`k=0`$. The linearized calculation predicts sub shot-noise statistics in the difference $`\widehat{N}_1(k)\widehat{N}_1(k)`$ for all wave numbers. For large wave numbers the analytical result for the correlation approaches the value $`1/2`$. It is interesting to keep in mind that for the OPO, the same quantity is equal to $`1/2`$ independently of the wave number . In the SHG case additional processes taking place in the cavity result in a smooth $`k`$-dependence of $`C_{11}^{()}(k,k)`$. These characteristics do not depend much on the value of pump, and are not changed significantly even when the pump level is taken beyond threshold, cf. Sec. VII. Therefore, the statistics of the intensity difference is not directly affected by the pattern formation mechanism. A radically different situation occurs for the sum-correlation $`C_{11}^{(+)}(k,k)`$, which shows a strong peak around $`k=k_c`$. For the pump value used in Fig. 10 the peaks correspond to a maximum value $`C_{11}^{(+)}(k_c,k_c)35`$. This behavior is connected with the increase of the fluctuations in the modes associated with the pattern instability when the threshold is approached, leading to a large excess noise in the statistics of the intensity of the individual modes $`[\omega ](k)`$ and $`[\omega ](k)`$. This excess noise in each intensity cancels when the difference $`\widehat{N}_1(k)\widehat{N}_1(k)`$ is considered leading to a sub-Poissonian statistics, while it is still present in the sum $`\widehat{N}_1(k)+\widehat{N}_1(k)`$. For large $`k`$ the correlation approaches 1.5, coinciding again with the corresponding value for the OPO. Finally, as before, the jumps at $`k=0`$ are due to contributions from the homogeneous steady states, cf. Eqs. (151)-(151).

The corresponding photon number variances $`C_{22}^{(\pm )}(k,k)`$ for the SH field are shown in Fig. 11. In contrast to the FH correlations there is almost no sub-shot noise behavior in the difference correlation $`C_{22}^{()}(k,k)`$. In other words, the SH beams only display very weak nonclassical correlations. As for the FH field, the emerging instability does not influence the noise level in $`C_{22}^{()}(k,k)`$, but $`C_{22}^{(+)}(k,k)`$ displays a large amount of excess noise in the vicinity of $`k_c`$. The asymptotic large $`k`$ behavior for both correlations $`C_{22}^{()}(k,k)`$ and $`C_{22}^{(+)}(k,k)`$ is analytically found to correspond to the shot-noise limit 1.0.

The cross-correlations $`C_{12}^{()}(k,k)`$ and $`C_{12}^{(+)}(k,k)`$ are shown in Fig. 12. The linearization approach predicts that these correlations are always above the shot-noise limit. Furthermore, at small wave numbers we note that the variances of the sum and difference coincide. This can only occur when the fluctuations in the individual modes $`[\omega ](k)`$ and $`[2\omega ](k)`$ are uncorrelated, what was indeed observed in Fig. 7. Moreover, both the sum and difference correlations show a large excess noise at $`k=k_c`$, which is slightly weaker for the difference, as the result of a partial noise cancellation.

The cross-correlations $`C_{12}^{()}(k,k)`$ and $`C_{12}^{(+)}(k,k)`$ are shown in Fig. 13, and here the difference correlations interestingly go below the shot-noise limit as long as $`k`$ is not too close to the critical wave number. It is worth pointing out that the difference $`\widehat{N}_1(k)\widehat{N}_2(k)`$ shows nonclassical behavior while the difference $`\widehat{N}_1(k)\widehat{N}_2(k)`$ (shown in Fig. 12) does not. This somehow paradoxical situation is related to what was observed in the normalized correlations where the cross-correlation between $`\widehat{N}_1(k)`$ and $`\widehat{N}_2(k)`$ was stronger than the almost vanishing cross-correlation between $`\widehat{N}_1(k)`$ and $`\widehat{N}_2(k)`$. At $`k=k_c`$ a large amount of excess noise dominates the behavior of both the sum and the difference correlation and the two correlations show a pronounced peak. For large $`k`$ the correlations approach the shot noise limit, as seen for the other cross-correlations in Fig. 12.

Olsen et al.  have investigated the system without spatial coupling corresponding to our results at $`k=0`$, and they find that, for certain detunings, the variance of the sum of the FH and SH intensities are more strongly quantum correlated than the variance of the individual intensities, due to the anti-correlation between them. $`\mathrm{V}ar[\widehat{N}_1(0)]/S.N.`$ and $`\mathrm{V}ar[\widehat{N}_2(0)]/S.N.`$ can be seen from Fig. 10 and 11, respectively, at $`k=0`$. Both are larger than the $`\mathrm{V}ar[\widehat{N}_1(0)+\widehat{N}_2(0)]/S.N.`$ observed in Figs. 12 and 13, so that our results confirm the ones of .

### B Nonlinear correlations: Numerical results

So far we have only considered the correlations predicted by the linearized equations. In order to go beyond this regime, we use our numerical simulations to search for nonlinear fingerprints in the correlations and in particular for the emergence of new correlations, i.e. $`C_{ij}^n(k,k^{})`$ with $`k\pm k^{}`$. Of particular interest is to look for correlations between the homogeneous steady states ($`k=0`$) and the states with $`k=\pm k_c`$, $`C_{ij}^n(0,\pm k_c)`$. From a technical point of view this task turned out to be difficult because nonlinear contributions to the correlation functions were only observable for pump values extremely close to threshold, in a region where the characteristic time of the dynamics diverges because of critical slowing down. This translates into very long transients and the need of equally long simulations.

We have observed some indication of nonlinear correlations for a pump $`E/E_t=0.99999`$, which became very clear when using $`E/E_t=0.999999`$. For this value of the pump, we show in Fig. 14 our results for $`C_{12}^n(k,k^{}=0)`$ and $`C_{12}^n(k^{}=0,k)`$: these curves put into evidence an anti-correlation between the modes $`[\omega ](\pm k_c)`$ and $`[2\omega ](0)`$, and between $`[2\omega ](\pm k_c)`$ and $`[\omega ](0)`$ . They present a very sharp peak structure, with a width determined by the distance between two adjacent points of the discretized $`k`$-space used for the simulations. This is due to the fact that we now consider the correlation functions at fixed $`k^{}`$ and let $`k`$ vary. These correlations are a result of nonlinear amplification of the diverging fluctuations as the threshold is approached. The negative nature of the correlation is connected with the fact that the fields with nonzero average values (here the homogeneous components) act as a ”reservoir” of photons for all processes occurring in the cavity. As we will show later, they are a precursor of the behavior of the correlations above the threshold. The correlations at $`k=0`$ correspond to the linear correlation shown in Fig. 9. The bottom plot in Fig. 14 shows the nonlinear correlations $`C_{ij}^n(k=+k_c,k^{}=0)`$ as the threshold is approached. The correlations are nonzero only for $`E/E_t>0.9999`$. As the nonlinear correlations set in the nonlinear channels in step 2) and 3) of Fig. 6 become stronger and this weakens the correlations induced by the channel of step 1), which is exactly what we observed in Fig. 9; $`C_{12}^n(0,0)`$ becomes less correlated very close to the threshold. Moreover, we see that the correlations $`C_{11}^n(0,+k_c)`$ and $`C_{12}^n(0,+k_c)`$ have almost identical values, and the same holds for $`C_{22}^n(+k_c,0)`$ and $`C_{12}^n(+k_c,0)`$. This interesting behavior can be traced back to the fact that close to the threshold the fluctuations $`\delta I_1(k_c)`$ and $`\delta I_2(k_c)`$ are perfectly correlated, as displayed by Fig. 8, whereas the slight anti-correlation between $`\delta I_1(0)`$ and $`\delta I_2(0)`$ is responsible for the lower values of $`C_{22}^n(+k_c,0)`$ and $`C_{12}^n(+k_c,0)`$ with respect to $`C_{11}^n(0,+k_c)`$ and $`C_{12}^n(0,+k_c)`$.

## VII Correlations above threshold

Above the threshold for pattern formation the linearized equations (III) are no longer valid. As displayed in Fig. 4, above the threshold not only the homogeneous modes, but also all modes with wave numbers $`k=\pm k_c,\pm 2k_c,\pm 3k_c,\mathrm{}`$, will present a macroscopic photon number. Linearizing around the steady state pattern solution above the threshold under the assumption of small fluctuations, one obtains new linear equations for the far field fluctuation amplitudes, which take into account three-wave processes such as $`[2\omega ](k_c)[\omega ](k)+[\omega ](k_ck)`$ or $`[2\omega ](k)[\omega ](k_c)+[\omega ](kk_c)`$. In analogy to the situation below the threshold a linear fluctuation analysis above the threshold predicts, in addition to the correlations already present below the threshold, the existence of additional correlations between the fluctuations $`\delta I_1(k)`$ and $`\delta I_1(k_ck)`$, and between $`\delta I_2(k)`$ and $`\delta I_1(kk_c)`$. We will not report here the explicit results of this cumbersome linear analysis and restore directly to the numerical analysis of the full nonlinear Langevin equations.

To investigate the implications of the new field configuration above the threshold on the intensity correlations, we first consider the correlations $`C_{ij}^n(k,k^{})`$. The same normalized correlations discussed in Fig. 7 below the threshold are plotted in Fig. 15 for a pump value above the threshold. We observe that the correlations at $`k=\pm k_c`$ decrease from their threshold value and are no longer perfect as they were at the threshold. A closer look actually reveals a dip in the correlations exactly at the pixels corresponding to $`k=\pm k_c`$. A tentative explanation for this is based on the fact that now the modes at the critical wave number have a finite average value, connected with macroscopic photon numbers in these modes, whereas the neighbouring pixels are significantly less populated, cf. the far field of Fig. 4. In comparison the normalized correlations $`C_{ij}^n(k,k^{})`$ show a much smoother behaviour around $`k_c`$. Hence, the observed reductions in the correlations above threshold at $`k=\pm k_c`$ are connected with spontaneous population exchanges between these macroscopically populated modes.

In Fig. 8 the peaks at $`k=\pm k_c`$ of Fig. 15 are followed as function of the pump. The behavior is very similar to what is seen below the threshold. Close to the threshold the correlations are perfect, and as the pump is taken further away from $`E_t`$ the correlations become weaker. Below the threshold this was explained through an eigenvalue competition, while above the threshold the explanation is that the competitions between the states become stronger.

The $`k=0`$ cross-correlation is plotted in Fig. 9, and above the threshold there is a loss of anticorrelation or there is even a small positive correlation. This might be attributed to the macroscopic and independent occurrences of the processes of step 2) and 3) in Fig. 6.

We saw in Sec. VI B nonlinear correlations just below the threshold, and in Fig. 16 the peaks corresponding to these correlations are plotted in order to follow the progress above the threshold. The strongest anti-correlation is observed just above the threshold, $`E/E_t=1.0001`$ and as the pump is increased the correlations become weaker due to increasing competition of processes involving higher harmonics. Moreover, the connection between the self-correlations and cross-correlations seen below $`E_t`$ only remains very close to the threshold, so as the pump is increased $`C_{11}^n(0,+k_c)C_{12}^n(0,+k_c)`$ and $`C_{22}^n(+k_c,0)C_{12}^n(+k_c,0)`$. This is related to the loss of perfect correlations away from the threshold.

In Fig. 17 the photon number variances $`C_{11}^{(\pm )}(k,k)`$ above the threshold are presented. Comparing these results with the corresponding ones below the threshold from Fig. 10 we observe that they are very similar. Generally, the correlation $`C_{11}^{()}(k,k)`$ does not change much with the pump level, and this fact has also been observed in the OPO . The sum correlation $`C_{11}^{(+)}(k,k)`$, however, contains peaks that are very sensitive to the pump level, both below and above the threshold. The behavior discussed here for the FH is also valid for the SH and the cross-correlations.

## VIII Conclusion and discussion

We have used the master equation approach to describe the spatiotemporal dynamics of the boson intra-cavity operators in second-harmonic generation, and we included in the model quantum noise as well as diffraction. Our study is based on the Q representation to describe the dynamics of the quantum fields in terms of a set of nonlinear stochastic Langevin equations for equivalent c-number fields. The choice of the Q representations gives some restraints on the parameter space, but we have checked that similar results are obtained using the approximated Wigner representation in other parameter regions.

A simple scheme describing the microscopic photon interaction that underlies the process of pattern formation has guided us in our analytical and numerical studies of the spatial correlations. Equal time correlations between intensity fluctuations were used to investigate the strength of the correlations between different modes. Also, possible nonclassical effects, such as twin beam correlations, were considered by calculating the photon number variances of the intensity sums and differences between spatial modes of the FH and SH fields.

We have found that at the threshold for pattern formation the Fourier modes with the critical wave number are perfectly correlated for the FH field, the SH field and also between the FH and the SH field. As the distance to the threshold is increased these correlations become weaker, which was shown analytically to be due to the competition of the eigenvalues of the linear system describing the system below the threshold. At large wave numbers, only the correlation between opposite points of the FH far field survives. This correlation is always found to be stronger than the others, which is consistent with the fact that the twin photon emission at the fundamental frequency is the primary source for correlations in the system. For far field modes around the critical wave number the self-correlations as well as the cross-correlations between FH and SH photons are linked to the pattern forming instability.

Very close to the threshold the linear analysis breaks down. The numerical simulations below the threshold showed the existence of nonlinear correlations which involve the $`k=0`$ mode and these are also seen above the threshold. The other correlations described above are also found above the threshold, but their strength decreases when moving away from the threshold. This can be understood from the fact that additional processes come into play, mainly consisting in population exchanges between the macroscopic fields at the critical wave number and its harmonics.

The intensity differences between opposite points of both the FH and SH far fields, as well as the cross-correlation between the two have been shown to exhibit nonclassical sub-shot noise behavior. These properties for the intensity difference turn out not to be sensitive to the process of pattern formation, since the corresponding correlations depend very weakly on the distance to the threshold and show no particular structure close to the critical wave number. The emerging pattern is connected with increased fluctuations in the modes with wave numbers around the critical wave number, leading to an excess noise in the corresponding individual intensities. Therefore, the sub-Poissonian statistics of the intensity differences reveal a partial noise cancellation. On the contrary, the sum of intensities clearly exhibit peaks around the critical wave number, originating from excess noise connected with the formation of a pattern.

In this work we considered equal time correlations calculated for the intra-cavity fields. This approach turned out to be very useful to understand the intra-cavity field dynamics. For the output fields we expect that the nonclassical correlations of the intra-cavity fields will remain below shot noise. The quantitative assessment of the amount of noise reduction or excess noise with respect to the shot noise level requires a specific additional calculation. For future work it would also be interesting to calculate the output fluctuation spectra at $`0`$ frequency for the difference and sum of intensities, which reflect the full amount of quantum correlations induced by the microscopic processes taking place inside the cavity, as for example considered for a vectorial Kerr model in .

## IX Acknowledgments

We acknowledge financial support from the European Commission projects QSTRUCT (FMRX-CT96-0077), QUANTIM (IST-2000-26019) and PHASE and from the Spanish MCyT project BFM2000-1108. We thank Steve Barnett and Pere Colet for helpful discussions on this topic.