{'fraction_non_alphanumeric': 0.03916722588670475, 'fraction_numerical': 0.012488338216279031, 'mean_word_length': 4.542271407837446, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bichromatic four-wave mixing and quadrature-squeezing from biexcitons in atomically thin semiconductor microcavities

Nonlinear optical effects such as four-wave mixing and generation of squeezed light are ubiquitous in optical devices and light sources. For new devices operating at low optical power, the resonant nonlinearity arising from the two-photon sensitive bound biexciton in a semiconductor microcavity is an interesting prospective platform. Due to the particularly strong Coulomb interaction in atomically thin semiconductors, these materials have strongly bound biexcitons and operate in the visible frequency range of the electromagnetic spectrum. To remove the strong pump laser from the generated light in optical devices or to simultaneously excite non-degenerate polaritons, a bichromatic-pump configuration with two spectrally separated pump lasers is desirable. In this paper, we theoretically investigate spontanous four-wave mixing and quadrature-squeezing in a bichromatically pumped atomically thin semiconductor microcavity. We explore two different configurations that support degenerate and non-degenerate scattering from polaritons into bound biexcitons, respectively. We find that these configurations lead to the generation strongly single- and two-mode quadrature-squeezed light.

Nonlinear optical effects such as four-wave mixing and generation of squeezed light are ubiquitous in optical devices and light sources. For new devices operating at low optical power, the resonant nonlinearity arising from the two-photon sensitive bound biexciton in a semiconductor microcavity is an interesting prospective platform. Due to the particularly strong Coulomb interaction in atomically thin semiconductors, these materials have strongly bound biexcitons and operate in the visible frequency range of the electromagnetic spectrum. To remove the strong pump laser from the generated light in optical devices or to simultaneously excite nondegenerate polaritons, a bichromatic-pump configuration with two spectrally separated pump lasers is desirable. In this paper, we theoretically investigate spontanous four-wave mixing and quadrature-squeezing in a bichromatically pumped atomically thin semiconductor microcavity. We explore two different configurations that support degenerate and non-degenerate scattering from polaritons into bound biexcitons, respectively. We find that these configurations lead to the generation strongly single-and two-mode quadrature-squeezed light.

I. INTRODUCTION
Strongly bound biexcitons in atomically thin semiconductor microcavities provide an avenue for low-power nonlinear optical devices, because the resonant scattering of unbound exciton-polaritons into a bound biexciton yields a powerful enhancement of the nonlinear optical response [1][2][3]. Previous analyses of degenerate four-wave mixing have shown that strong parametric gain from biexcitons [4,5] can provide quadrature-squeezing with significantly lower power than when using conventional third-order nonlinear materials [6]. Thereby, atomically-thin semiconductors have a significant potential as a platform for nonlinear optics. In the interest of spectrally separating the pump photons from the generated signal, it is highly convenient to use a bichromaticpump scheme, where two non-degenerate pumps create a degenerate signal at the mean of the two pump photon energies [7][8][9][10][11][12]. In addition, such a bichromatic-pump setup allows to excite non-degenerate polaritons, which opens up a non-degenerate nonlinear scattering channel [13].
Four-wave mixing in pump-probe experiments with semiconductors [14][15][16][17][18][19][20] and semiconductor microcavities [21,22] is a heavily studied theoretical and experimental topic, and it was early recognized that the bound biexciton plays a paramount role in four-wave mixing. However, spontaneous four-wave mixing, also known as parametric flourescence or hyper-raman scattering, in semiconductors is far less studied [23,24], despite being a central topic for quantum light sources with conventional nonlinear materials [25][26][27][28]. Fourwave mixing is a third-order nonlinear process that transforms two pump photons with frequencies ω 1 , ω 2 to a pair of signal/idler photons with frequencies ω 3 , ω 4 . Stimulated fourwave mixing is facilitated by applying an idler laser field with frequency ω 3 , which will produce a signal at ω 4 = ω 1 + ω 2 − ω 3 due to energy conservation. In contrast, in the absence of a stimulating idler field, spontaneous four-wave mixing produce a broad and continuous spectrum of photon pairs. This is the regime that is studied in this paper. Due to the pairwise creation of photon pairs, the field properties of the generated light can exhibit strong non-classical correlation signatures in the form of quadrature squeezing, two-mode squeezing and photon-number correlations. In semiconductors, particularly in semiconductor microcavities with a strong light-matter interaction, spontaneous four-wave mixing can be mediated by the biexciton, which can break into a correlated pair of polaritons [see Fig. 1 In this paper, we theoretically investigate spontaneous fourwave mixing in bichromatically pumped semiconductor microcavities with atomically-thin semiconductors. We investigate the spectrum and quadrature squeezing of the generated light and discuss how the cavity resonance and pump frequencies can be tuned in two configurations to facilitate efficient light generation and squeezing. Specifically, we investigate two configurations close to resonances, where four-wave mixing and squeezing is efficient. In the first configuration [Configuration A, cf. Fig. 1(c)], the cavity frequency is tuned such that the lower polariton (LP) energy E − 0 matches half the bound biexciton energy 1 2 E xx b,− . Here, the pump lasers should be nearly degenerate around the lower-polariton energy in order to efficiently excite LPs, which scatter via the Coulomb interaction into bound biexcitons through a degenerate polaritonic Feshbach resonance [1,29]. A Feshbach resonance occurs when the energy of an open scattering channel (in this case a LP pair) matches the energy of a bound multiparticle complex (here a bound biexciton) [30]. The biexciton decays spontaneously into two degenerate LPs, which leads to single-mode squeezed light emission. In the second configuration [Configuration B, cf. Fig. 1(d)], the cavity frequency is tuned such that the sum of the LP and upper polariton (UP) energies E − 0 + E + 0 matches the bound biexciton energy E xx b,− , thereby giving rise to a non-degenerate polaritonic Feshbach resonance, also known as a polaritonic cross Feshbach resonance [13]. The driving field is resonant, when the two pump laser energies match the two polariton energies. An uncorrelated UP-LP pair can then resonantly scatter into a bound biexciton, which decays into a non-degenerate UP-LP pair . The biexciton spontaneously breaks into a pair of polaritons via the Coulomb interaction W − b , which are outcoupled from the cavity as photons (blue and green wiggly arrows). (c) Polariton energy bands (E ± q , black lines) in Configuration A, where the cavity frequency is tuned in order to bring the lower polariton branch into Feshbach resonance with the bound biexciton (blue dotted line), 2E − 0 = E xx b,− . The driving frequencies are symmetrically centered around the lower polariton (orange and red arrows). The uncoupled cavity and exciton energy bands are shown with grey lines. (d) Similar to (c) for Configuration B, where the cavity frequency is tuned such that an upper-lower polariton pair matches the energy of a bound biexciton, E − 0 + E + 0 = E xx b,− , forming a non-degenerate Feshbach resonance. The driving frequencies are symmetrically centered around 1 2 E xx b,− with a resonance occurring when that leads to the emission of two-mode squeezed light. In this context, two-mode quadrature-squeezing refers to a reduction of the quadrature noise in the cross-correlation between spectrally distinct frequency bands [31]. For the investigation, we employ a rigorous perturbative expansion of the electronic and photonic correlations up to third order in the driving field through the dynamics-controlled truncation (DCT) scheme [32,33]. In order to calculate not only intracavity dynamics, but also the properties of the outcoupled and thus detectable field, we combine DCT with a Heisenberg-Langevin approach [6,34]. Due to the presence of two non-degenerate pumps, the steady state of the system is not constant [35], but can be expressed as a discrete Fourier series. This series expansion complicates the solution of the Heisenberg-Langevin equations of motion in frequency space compared to the case of a constant steady-state, which we solve through a discrete Fourier-index formalism.
The paper is organized as follows: In Sec. II, we describe the semiconductor model that is used for the electronic states in the semiconductor, the photonic states in the cavity and the external driving field. In Sec. III, the DCT and Heisenberg-Langevin methods are applied to calculate the dynamics, spontaneous four-wave mixing spectrum and squeezing. In Sec. IV, we present and discuss the results for Configurations A and B. Finally, we conclude in Sec. V.

II. MODEL
Here, we describe the used two-band semiconductor model for electrons and holes, as well as the coupling to cavity photons and the external laser drive.

A. Hamiltonian
We consider an atomically thin semiconductor placed in a planar cavity, which is driven with two coherent pump beams with frequencies ω 1 and ω 2 at normal incidence. The electromagnetic field in the cavity is quantized through the bosonic annihilation and creation operators a σ ,k and a † σ ,k , which describe cavity photons with polarizaton σ and in-plane momentum k. We assume the cavity to be rotationally symmetric in the plane, such that the cavity mode is polarizationdegenerate.
The semiconductor is described by electronic states in the conduction and valence bands with fermionic annihilation and creation operators c ζ k , c † ζ k for the conduction band and v ζ k , v † ζ k for the valence band. The index ζ = (ξ , s) labels valley (ξ ) and spin (s). We shall restrict our analysis to the lowest-energy optical transitions in transition-metal dichalcogenide monolayers, which are located at the valleys ξ = K and ξ = K , respectively. Due to spin-orbit coupling, the lowestenergy transitions allow photons with right-hand (σ = R) or left-hand (σ = L) circular polarization to excite electron-hole pairs with spin-valley combinations ζ = (K, ↑) and ζ = (K , ↓ ). Due to this unambiguous relation between photon polar-ization and electron spin and valley, we can absorb photon polarization into the index ζ and use the shorthand notation ζ ∈ {K, K } to denote all three degrees of freedom.
The total Hamiltonian H is divided into a noninteracting term H 0 and a Coulomb interaction term H C , such that H = H 0 + H C . The noninteracting Hamiltonian describing electrons, photons and their coupling is given by where the first two terms account for the energy of electrons in the conduction and valence band, the quasiparticle bandgap and m e , m h the effective electron and hole masses. The third term accounts for the energy of cavity photons, E p k =h[ω 2 p,0 + (ck/n) 2 ] 1/2 , where ω p,0 is the resonance frequency at k = 0, c is the speed of light andn is the effective refractive index of the cavity mode [36,37]. The last two terms describe electron-hole-photon coupling with strength which depends on the out-of-plane confinement of the cavity mode and the Bloch momentum matrix element of the semiconductor [38,39].
The Coulomb interaction is given by [40] Here, V q = e 2 0 [2Sε 0 ε q q] −1 is the screened 2D Coulomb potential, where e 0 is the elementary charge, S is the quantization surface area, ε 0 is the vacuum permittivity, and ε q is the dielectric function for 2D semiconductors which is described in Appendix E. The inter-and intra-valley Coulomb exchange interaction has been neglected here, because in transitionmetal dichalcogenides it is significantly weaker [20,40] than the direct interaction in Eq. (2). We note that the exchange interaction leads to effects such as biexciton fine structure [41] and a splitting of the exciton into branches with linear and quadratic dispersion [42] that coincide at zero center-of-mass momentum. While these effects are rich and interesting, we consider here only the physics of the dominating nonlinear response from the direct Coulomb interaction and leave the additional inclusion of exchange effects to more detailed future analyses.

B. Driving field
The two pump laser drives are introduced via the inputoutput formalism [43,44], which is based on the microscopic interaction between the internal cavity mode with the quantized continuum of external modes. By formally solving the equation of motion of the external field operators, the external input field is linked to the equation of motion for the cavity field operator as where the last two terms describe outcoupling from the cavity mode and incoupling of the driving field, respectively. We take the expectation value of the driving field to be in a coherent state and to have the bichromatic form where ω 1 and ω 2 are the frequencies of the two driving lasers and a in ζ 1 and a in ζ 2 are the corresponding amplitudes. We shall take the total power in each of the driving lasers to be equal and denote by λ in i , i = 1, 2 the polarization vector of the two drives in the circular basis. We then express the input field components as polarization vectors as where the superscript T denotes transposition. The input-field vectors are related to the total driving power P in as [43]

III. METHODS
Here, we introduce the DCT scheme for time evolution of the expectation values and the Heisenberg-Langevin equations for the fluctuation operators. The model and techniques are a generalization of the work presented in Ref. 6 to bichromatic driving fields.

A. Time evolution
To calculate the time evolution of the system, we apply the DCT scheme [32,33] to perturbatively expand the equations of motion for the coherent expectation values to third order in the driving field a in ζ . The first step in this procedure is the Heisenberg equation of motion for a general operator Q, −ih∂ t Q = [H, Q] (see Ref. 6 for the explicit derivation).
The Coulomb interaction and the fermionic commutation relations generate expectation values with an unequal number of conduction and valence band operators c ζ ,k and v ζ ,k , such Here, conduction-band electron densities c † ζ ,k −q c ζ ,k are expressed perturbatively in terms of electron-hole pair operators using a unit-operator expansion [40,45] 4 ], and similarly for valence-band hole densities. In this expansion, the Hilbert space has been restricted to pairs of conduction band electrons and valence band holes [40], which means that the effects of unpaired electrons or holes have not been accounted for. In this paper, this assumption is valid, because the only source of carriers is optical excitation through the driving field, which creates electrons and holes in pairs. However, in the presence of unpaired electrons or holes through e.g. n-or p-doping, which we do not consider in this paper, the pair expansion is no longer valid. Within the DCT scheme, the third-order perturbative expansion corresponds to only keeping terms with up to three normal-ordered electron-hole pair or cavity photon operators [46], since the number of pair or photon operators corresponds to the order in the driving field. Such a perturbative expansion is possible because all correlations have their origin in the action of the external driving field [32,33]. Furthermore, in the coherent-response limit, which is considered here, the third-order expectation values are systematically factorized as This third-order factorization is exact in the coherent regime [46], i.e. when the only source of electrons and holes is excitation with coherent light near the resonances of the system, and when incoherent scattering processes e.g. via phonons can be neglected. Phonon scattering is later included phenomenologically through an exciton dephasing rate, which is obtained from a separate self-consistent microscopic calculation [47]. This means that the validity of our approach is limited to the regime where cavity outcoupling (γ p ) dominates over exciton dephasing (γ x ), such that polaritons will be outcoupled before significant dephasing and scattering takes place. In practise, this limits the theory to the low-temperature limit, where phonon dephasing is slow. For the calculations presented in this paper, the temperature is 30 K, which yields a phonon-induced dephasing rate below 1 meV, while the cavity outcoupling is 9 meV, consistent with fabricated dielectric microcavities with transition-metal dichalcogenide monolayers [48]. Thus, there is a separation between these two time scales by an order of magnitude.
The electron-hole are expanded on the complete set of excitonic eigenstates as k is the ith exciton wavefunction (solution to the Wannier equation) in momentum space with zero-momentum energy E x i,0 and P i † ζ ,q is the creation operator of the corresponding exciton with center-of-mass momentum q, and α = m e /(m e + m h ) and β = m h /(m e + m h ). We restrict the analysis to the lowestenergy (i = 1s) exciton, which is energetically separated from the next excitonic state by hundreds of meV for atomically thin transition-metal dichalcogenides [49,50] and we shall thus drop the index i from here onwards, meaning that the electronic states are projected onto the lowest exciton wavefunction. This is valid, when the excitation energy is far away from any of the neighbouring exciton states. In practise, the projection of electronic states onto the lowest-energy exciton wavefunction ensures an efficient numerical evaluation of the dynamics and exciton-exciton scattering coefficients.
Two-photon expectation values are partitioned into factorized parts and correlations, defined as D ζ ζ Three-particle photon-electron-hole correlations are defined as i.e. with the contribution that can be factorized in electronhole pairs and photons subtracted.
These correlations are projected onto the exciton wavefunction as C ζ ζ q := ∑ k φ k a † ζ ,q c † ζ ,k−αq v ζ ,k+β q c . The four-particle correlations of two electron-hole pairs have a more complicated structure due to the two possible electron-hole pairings: These correlations are projected on the 1s exciton wavefunction in terms of the biexcitonic correlationsB ζ ζ q,± in the triplet (+) and singlet (−) linear combinations through the relation [51] Due to the Coulomb interaction, the equation of motion for the biexcitonic correlationB ζ ζ q,± is coupled to correlations with different momenta,B ζ ζ q ,± . To alleviate this momentum offdiagonal coupling and thereby simplify the structure of the equations of motion, the biexcitonic correlations are expanded on the biexcitonic wavefunctions Φ ± µ,q , which are the solutions to an effective two-exciton Schrödinger equation [20] with corresponding energies E xx µ,± (see Appendix A for details). Bound (µ = b, E xx b,− < 2E x 0 ) and unbound (E xx µ,− > 2E x 0 ) solutions exist in the singlet channel, where the effective exciton-exciton Coulomb interaction is attractive. The triplet channel supports only unbound solutions [52], because the effective interaction is repulsive. The unbound solutions constitute a two-exciton scattering continuum. In this biexcitonic eigenbasis, the correlations are expressed asB At the level of the equations of motion, we phenomenologically include phonon-induced broadening of the exciton and biexciton by introducing complex-valued energies in the equations of motion: The broadening γ x is calculated through a microscopic, selfconsistent approach for the phonon interaction as in Refs. 47, 53-56. We approximate the biexcitonic broadening as twice the exciton broadening [57][58][59].
Similarly, photon outcoupling from the cavity is included in the complex cavity frequencyẼ The equations of motion of the photon and exciton amplitudes and the three types of correlations as described above form a closed set within the third-order DCT scheme [46]. In a rotating reference frame with respect to the mean drive fre-quency ω r = (ω 1 + ω 2 )/2 the equations of motion are The full details of the derivation can be found in Ref. 6, but for completeness, all quantities are defined in Appendix A. In each equation, the first term on the right-hand side describes free evolution, and the remaining terms describe couplings or driving. For the photon amplitude a † ζ ,0 , the second term is the linear vacuum-Rabi coupling to the exciton with coupling strength Ω 0 (where 2Ω 0 is the vacuum Rabi splitting). The last term is input-field driving, where in the rotating frame, the two pumps rotate with ±ω 12 , with the difference frequency For P † ζ ,0 , the second term describes vacuum Rabi coupling. The third term arises from the fermionic substructure of excitons and generates nonlinear saturation of the light-matter interaction due to Pauli blockingΩ q . The last two terms describe Coulomb exciton-exciton interactions at the mean-field level (W 0 ) and Coulomb-induced interactions with the biexcitonic correlations (W ± µ ) beyond mean-field. For the biexcitonic correlations B ζ ζ µ,± , the second term contains Coulomb-scattering of uncorrelated excitons (W ± µ ) and coupling to exciton-photon correlations through the lightmatter interaction (Ω ± µ,q ). For the exciton-photon correlations C ζ ζ q , the second term describes linear coupling to two-photon correlations by exchanging an exciton with a photon (Ω 0 ). The third term describes a nonlinear scattering of two uncorrelated excitons (Ω q ), and the last term describes coupling to biexcitonic correlations via optical fields (Ω ± µ,q ). The second and third terms in the equation of motion for the two-photon correlations D ζ ζ q describe coupling to exciton-photon correlations by exchanging a photon with an exciton through the light-matter coupling Ω 0 .
B. Steady-state discrete Fourier decomposition As noted, the driving terms oscillate with the frequency ±ω 12 in the rotating frame. In contrast to the single-pump case [6] where the driving term is constant in the rotating frame, this means that the steady state evolution of the expectation values is not constant, but contains terms oscillating with integer multiples of ω 12 . This means that the steady-state dynamical variables Y (t), where Y is any of the expectation values in Eq. (5), are periodic with period T = 2π/|ω 12 | and can thus be represented as a discrete Fourier series In order to find the steady-state discrete Fourier series, Eq. (5) is propagated numerically in time until the coefficients Y m have converged sufficiently.

C. Polaritons
When photons and excitons are strongly coupled, as considered in this paper, the relevant linear-response eigenstates are polaritonic states. These states emerge naturally out of the equation of motion by diagonalization in the linear limit. Taking the first two lines of Eq. (5), removing the nonlinear terms and the driving, and casting them in the general nonzero momentum form in the laboratory (i.e. non-rotating) reference frame, we have In the limit where photon outcoupling and exciton dephasing are neglected, these equations of motion can be diagonalized through a unitary Hopfield transformation [60] by introducing the polariton operators Ξ ± † ζ ,q = u p± q a † ζ ,q + u x± q P † ζ ,q where the expansion coefficients are given by In this basis, the linear evolution of Eq. (7) becomes with the polariton energies The commutation relations of the polariton operators are not bosonic, but have correction terms due to the fermionic substructure of the exciton operator P † ζ ,q [40]. In principle, the full equations of motion Eq. (5) can be expressed in the polaritonic basis. However, this is not necessary and will not change the dynamics itself, only the basis that it is expressed in. In this work, we only use the polaritonic energies to identify the linear resonances of the system and to understand the nonlinear scattering between zero-momentum polaritons and the bound biexciton as depicted in Fig. 1.

D. Fluctuations
With the equations of motion of the single-time expectation values as in Eq. (5), one can access the instantaneous quantum statistics in the steady state or in the transient evolution. This can in principle be used to calculate the intracavity squeezing or photon number. However, the relevant detectable quantity is not the intracavity field, but the outcoupled field. Since the cavity field at zero-momentum couples out of the cavity into a continuum of external radiation modes with different frequencies, outcoupled quantities are always described and often measured by a spectrum rather than a single number, in contrast to intracavity quantities, which can be expressed in terms of a single mode [44,61].
To calculate the spontaneous four-wave mixing spectrum and quadrature squeezing of the generated light, we need not only the single-time expectation values as in Sec. III A, but also multitime averages. To calculate these, we employ a Heisenberg-Langevin approach for the time evolution of the fluctuation operator of the cavity field, δ a † ζ , and similarly for the exciton fluctuation operator δ P † ζ ,0 . Importantly, the cavity in/out-coupling and the phonon-induced dephasing give rise to Langevin noise sources [62]. For the cavity field, the Langevin noise enters directly from microscopic theory through the input-output formalism, as seen in Eq.
is a Langevin noise term with the properties [43,44] For the phenomenologically introduced phonon-broadening, we introduce a similar Langevin noise term to accompany the decay process, where the noise correlation properties are derived from the rules in Ref. 62 . Calculating the commutator in Eq. (10) in a third-order expansion in the input field, it is found that the nonlinear terms couple δ P † ζ ,0 to exciton-photon pair fluctuation operators δ C δ Pδ B, δ P in † δ P † , δ a in † δ a † , δ P in † δ a † and δ a in † δ P † can be neglected and a set of equations that are linear in the fluctuation operators is obtained, Being interested in the steady-state emission properties, we consider the limit of t → ∞ and write all expectation values using their discrete Fourier series (see Eq. (6)). We then transform the fluctuation equations to Fourier space as where δ Q is any of the fluctuation operators. Terms on the right-hand side of Eq. (11) of the form δ Q(t)Y (t) transform as ∑ m δ Q(ω − nω 12 )Y m , where Y is a steady-state expectation value and Y m is its discrete Fourier decomposition from Eq. (6). Thus, the fluctuation equations are not diagonal in frequency space, since δ Q(ω) is coupled to δ Q(ω + nω 12 ) for n ∈ Z.
To handle this challenge, we partition the frequency axis into zones of size ω 12 and define the frequency ν to be in the interval [−ω 12 /2, ω 12 /2]. We then introduce the discrete . Instead, the frequency offdiagonal coupling in the fluctuation equations appears as a coupling between different Fourier indices [see Fig. 2(c)]. At the same time, we suppress the momentum-index 0 on δ P and δ a and use only the Fourier index, such that δ P † ζ ,m (ν) and δ a † ζ ,m (ν) implicitly refer to the zero-momentum exciton and cavity fluctuation operators.
By solving the frequency-space equations for δ B  For a compact notation, we introduce the combined fluctuation vector The frequency-space Heisenberg-Langevin equation for δ ψ reads where the inverse Green's function is given by where I 4 is the 4 × 4 identity matrix, ν m := ν − mω 12 , Σ mm ζ ζ (ν) is the exciton self-energy, andΩ mm ζ ζ ,0 is a renor-malised exciton-photon coupling strength. The latter two stem from renormalisations due to the formal elimination of the multiparticle fluctuations and are given in Appendix B.
These renormalisations give small quantitative corrections to the spontaneous four-wave mixing spectra. Most importantly, ∆ m ζ ζ is the parametric gain connecting different fourier components m 1 and m 1 + m, which represents the core mechanism behind spontaneous four-wave mixing. It is given by The parametric gain describes a coherent pairwise driving of the excitons, and is a purely a result of the nonlinearities in the system, specifically the Coulomb interaction and the Pauliblocking from the fermionic substructure of excitons. Moreover, without parametric gain, the generated field would have the same spectral and coherence properties as the classical input field. The parametric gain contains terms from the instantaneous Coulomb interaction between excitons (W 0 ), biexcitonic correlations (W ± µ ) from the bound biexciton (µ = b) and two-exciton continuum (µ = b), and Pauli-blockingΩ q . All matrix elements and coefficients are given in Appendices A and B.
The matrix T mm ζ ζ (ν) describes the coupling to the input field, δ ψ in der of the input field is given by T mm Due to elimination of the multiparticle fluctuations, additional second-order contributions are also present, which are given in Appendix B.
The formal solution of Eq. (13) is given by where G mm ζ ζ (ν) = ∑ m ζ G mm ζ ζ (ν)T m m ζ ζ (ν). We note that in practise, the renormalizations from formal elimination of the multiparticle fluctuations δ B, δ C and δ D influence the numerical calculations presented in this paper only with small quantitative corrections. This means that in many cases, one can neglect δ B and δ C in Eq. (11), leading to the simplifications This is demonstrated explicitly in Appendix D, where we compare the full calculation of fluctuation spectra with the calculation under the simplifications in Eq. (17). All calculations presented in the main text have been performed using the full contributions of Eq. (14). However, the negligible influence of multiparticle fluctuations supports the linearization of the Heisenberg-Langevin equations, since products such as δ Cδ P can be expected to have an even smaller effect than the terms linear in δ C. This fact does not mean that many-body correlation effects are unimportant: indeed, the dominating contribution to the parametric gain stems from the bound biexciton,

E. Emission spectrum and squeezing
The detected field is the output from the cavity perpendicular to the surface, i.e. at zero in-plane momentum k = 0. We take the field to be polarisation-filtered before detection and denote the polarisation-projected cavity operator by a = λ T out a 0 [44], where a 0 = (a K,0 , a K ,0 ) T and λ out is the polarisation of the detection channel in the circular basis. The corresponding fluctuation operator is δ a(t) = lim t→∞ [a(t) − a(t) ].

Spectral correlation functions
The emission properties of the generated field from the cavity are characterized by two spectral correlation functions where j ∈ [1,4] is the vector index of δ ψ as defined in Eq. (12). As shown in Sections III E 2 and III E 3, the emission spectrum is described through the correlation function in Eq. (18), whereas the squeezing spectrum is described through both correlation functions in Eqs. (18) and (19). The spectral correlation of the input field is derived from the temporal correlations from Sec. III D as where η j j = δ j ,1 δ j,3 + δ j ,2 δ j,4 . Although no correlations are present between photon and exciton input fields, such correlations will generally be nonzero in the internal field δ ψ. With this, we obtain where S mm For convenience, we shall use the shorthand notation in order to label the fluctuation spectra on the continuous frequency axis rather than in the Fourier-zone partitioning. This relation can alternatively be cast as .

Emission spectrum
The emission spectrum of the outcoupled field S tot (ω) = S coh (ω) + S FWM (ω) contains a coherent contribution S coh = S coh,1 δ (ω − ω 1 ) + S coh,2 δ (ω − ω 2 ), which shares the spectral distribution and temporal coherence of the driving field [43,44], and a spontaneously generated field S FWM , which is created by spontaneous four-wave mixing. The former is of minor interest in the present investigation, because it simply shares the quantum statistics with the classical input field. The primary emission spectrum of interest is the spontaneous fourwave mixing spectrum, which is given by where S 1 (ω) is defined in Eq. (24).

Homodyne noise spectrum and quadrature squeezing
Quadrature squeezing is described through homodyne detection of the output field from the cavity (see Fig. 3). While the quadrature squeezing of the internal cavity mode can be evaluated simply as the variance of the quadrature operator δ X(θ ) := e iθ δ a † + e −iθ δ a with respect to the quadrature phase θ , the outcoupled and thus measurable field quadrature fluctuations are more complicated [44,63]. In homodyne detection, the source field from the cavity E s (t) = √ 2γ p a(t) is mixed with a strong local oscillator E lo (t) on a beamsplitter, such that the fields leaving the beamsplitter E 1 and E 2 are given by

BS
Det. 1 Det. 2   FIG. 3. Homodyne detection. The source field E s is mixed with a strong, coherent local oscillator field E lo on a 50/50 beamsplitter (BS). The outgoing field modes E 1 and E 2 are detected on detectors 1 and 2, which produces the photocurrents I 1 and I 2 . The homodyne signal is the current difference I 1 − I 2 .
The local oscillator has a frequency equal to the rotating frame frequency 1 2 (ω 1 + ω 2 ). Thus, when working in the rotating frame, the local oscillator expectation values are timeindependent: E lo (t) = e iϕ √ F lo , where F lo and ϕ are the photon flux and reference phase of the local oscillator, respectively, and the expectation value is evaluated in the rotating frame. Taking the local oscillator to be in a coherent state, we have E † lo (t)E lo (t) = F lo . In the homodyne measurement process, the two output fields E 1 and E 2 are measured on separate photodetectors (which are labeled with corresponding indices 1 and 2) and the difference of these photocurrents are recorded as the signal.
The statistical properties of this photocurrent difference are derived using normal-ordered detection theory [63][64][65][66][67]. Here, we extend the derivation of Ref. 63 to the case of balanced homodyne detection and to bichromatic driving, where no time-independent stationary source field exists. The detection model is based on the assumption that a single photoelectric detection event produces a current pulse of duration τ d and amplitude ge/τ d , where e and g denote the electronic charge and photodetector gain, respectively. The photocurrents of detector µ = 1, 2 is then given by [63,65] where n µ is a classical stochastic variable representing the number of overlapping pulses in the detection electronics of detector µ, i.e. the number of pulses initiated in the time interval t − τ d to t. For practical purposes, we shall work with a scaled current, I µ (t) := I µ (t)/(ge), which has units of inverse time and represents the flux of detected photons. The joint probability of detecting n 1 photons in detector 1 and n 2 in detector 2 within the time interval [t − τ d ,t] is given by [63,65] p (1) (n 1 ,t − τ d ,t; n 2 ,t − τ d ,t) where the generator G of the distribution is given by and η is the detection efficiency. The symbol :: denotes normal-and time-ordering at the level of the electric field operators (E † µ arranged to the left of E µ and time arguments increasing to the right (left) in products of E † µ (E µ )). The probability distribution Eq. (28) is essentially Mandel's counting formula, as presented in e.g. Ref. [44]. From p (1) , the mean photocurrent is derived as [63] I µ (t) = 1 where we note that we use overlines to denote the mean value of the classical stochastic variable I µ and angular brackets to denote quantum-mechanical expectation values. Similarly, the two-time probability of detecting n counts in detector µ in the time interval [t − τ d ,t] and m counts in detector µ in the time interval [t + τ − τ d ,t + τ] is given by [63] p From this probability, the detector current correlation function is found to be (for τ > 0) where Θ is the Heaviside function. The measured spectral noise function N(ω) of the homodyne signal I − := I 1 − I 2 is then the Fourier transformation of the photocurrent fluctuation correlation function, averaged over one period of the signal oscillation T = 2π/|ω 12 | = 4π/|ω 2 − ω 1 |: From Eqs. (26), (30) and (32), the following expression for the photocurrent fluctuation correlation function is derived (see Appendix C) where X s (θ ,t) = e iθ E † s (t) + e −iθ E s (t) is the source-field quadrature operator with angle θ = ϕ + π/2 and its fluctuation operator is δ X s (θ ,t) := X s (θ ,t) − X s (θ ,t) . The measured quadrature angle θ stems from the local oscillator phase ϕ and a phase displacement π/2 from the beamsplitter as seen in Eq. (26). The first term in Eq. (34) is the shot noise correlation function, which in the limit where the local oscillator is much stronger than the source field is given by In order to characterize the intrinsic noise properties of the source field and to disentangle these properties from the detector properties, we take the limit of infinite detection bandwidth (τ d → 0) and unity detection efficiency (η = 1). In this limit, inserting Eq. (34) into Eq. (33) yields the noise spectrum relative to the shot-noise level (see Appendix C) where N 0 (ω) = ∞ 0 dτ cos(ωτ)N 0 (τ) is the shot noise spectrum and where Λ i (ω) = 4 ∞ 0 dτ cos(ωτ)C i (τ) and C i (τ) = 1 2π ∞ −∞ dω e iωτ 2γ p S i (ω). All details of the derivation are given in Appendix C.
The optimal quadrature angle, where the noise is minimized, is given by e −2iθ = −Λ 2 (ω)/|Λ 2 (ω)|. At this angle, the spectrum of squeezing is where the homodyne phase θ has been suppressed in Λ(ω, θ ) for notational brevity. In the analysis of squeezing, we shall use the normalized optimal-angle spectrum 1 + Λ(ω) to characterize the amount of squeezing in the generated light. This spectrum is normalized such that a value of 1 + Λ(ω) = 1 corresponds to the shot noise level, i.e. no squeezing at all, whereas a value of 1 + Λ(ω) = 0 corresponds to complete elimination of the shot noise and thus perfect squeezing. We note that the spontaneous four-wave mixing spectrum S FWM (ω) and the squeezing spectrum Λ(ω) are of fundamentally different nature and are measured with very different techniques. The spontaneous four-wave mixing spectrum S FWM (ω) can be detected by filtering out the bichromatic pump from the generated light and sending the filtered signal into an optical spectrum analyser. The frequency argument ω is connected to the emitted photon energy on the order of 2 eV.
In contrast, the squeezing spectrum Λ(ω) is measured through homodyne detection, i.e. by beating the generated signal with a local oscillator at frequency ω r := (ω 1 + ω 2 )/2. Here, the frequency argument ω is the beat frequency between the signal and the local oscillator, and is thus much smaller, on the order of a few meV.

IV. RESULTS
By numerically solving the equations of motion Eq. (5) explicitly to reach the periodic steady state (see Appendix E), the corresponding Fourier components are calculated from the resulting time series over a steady-state period from Eq. (6). The steady-state Fourier components are then used to construct the fluctuation Green's function from Eq. (14), and calculate the spontaneous four-wave mixing spectrum S FWM (ω) and the squeezing spectrum Λ(ω) as described in Sec. III E. In this section, we present these results for Configurations A and B, respectively, as shown in Fig 1(c)-(d). The numerical calculations become increasingly challenging ashω 12 := (ω 2 − ω 1 )/2 approaches 0, because the oscillation period becomes longer and longer. However, in the case wherehω 12 = 0, the spectra S FWM (ω) and Λ(ω) can be calculated using the simpler calculation methods for a monochromatic driving field described in Ref. 6. For completeness, we include the results forhω 12 = 0 in this way. The calculations are performed for atomically thin MoS 2 encapsulated by hexagonal BN on both sides. The coupling strength of the cavity is taken to be Ω 0 = 20 meV and the cavity outcoupling rate is taken to bē hγ p = 9 meV, consistent with fabricated devices [48,68]. The temperature is taken to be 30 K, leading to a phonon-induced exciton dephasing ofhγ x = 0.8 meV (see Appendix E). The pump polarization is taken to be linear, with the two drives having the same linear polarization, and the pump power is 5 mW with a pump spot size of 9 µm 2 . The detected polarization is taken to be co-linear with the drive, and we note that very similar results are seen for cross-polarized detection. Further details about the numerical calculations and all parameters are given in Appendix E.

A. Configuration A
Configuration A is defined by setting the cavity frequency such that the lower polariton energy E − 0 matches half the bound biexciton energy 1 2 E xx b,− , and setting the mean value of the driving fields to 1 2 Fig. 1(c)]. When the drive energy differencehω 12 is zero, the degenerate driving field resonantly creates lower polaritons, which resonantly scatter into bound biexcitons via the Coulomb interaction. The latter process is also known as a polaritonic Feshbach resonance [1,2,29] and gives a strong enhancement of the nonlinear response, which in this case leads to strong spontaneous four-wave mixing and singlemode squeezing. Since the in-scattering particles in this configuration (lower polaritons) are identical, this is a degenerate Feshbach resonance. In Fig. 4(a), the spontaneous four-wave mixing spectrum is shown for four values of the drive frequency differencē hω 12 . The spectra feature multiple peaks, which we can identify by considering the resonant four-wave mixing channels. Due to energy conservation the energy of a generated photon pair with frequencies ω 3 and ω 4 must match the sum of two pump photon energies. With two different pump frequencies ω 1 and ω 2 , there are three possible combinations of pump photon pairs, meaning that the output photon pair frequencies must fulfill one of the relations The resonant output channels appear where one of the output photon energies equals a polariton energy. Thus, by settinḡ hω 3 = E ± 0 , the resonance conditions Eq. (39) become These four-wave mixing resonance conditions give rise to a total of eight possible output photon frequencies (two values of ω 3 and three values of ω 4 for each of these). In Fig. 4(a), the resonances E ± 0 andhω 1 +hω 2 −E ± 0 are indicated with vertical grey lines. These resonance frequencies are independent ofhω 12 . The strongest (middle) resonance peak is the lower polariton energy E − 0 , which in Configuration A matches the mean drive energy, 1 2 (hω 1 +hω 2 ). The corresponding photon pair partner from Eq. (40a) has the same energy. The upper polariton resonance E + 0 and its partnerhω 1 +hω 2 − E + 0 are positioned symmetrically around the lower polariton. Since the excitation of the upper polariton is much further from resonance, these peaks appear with a much weaker amplitude than the central lower polariton peak.
In addition, the resonance condition in Eq. (40c) for the lower polariton, 2hω 2 − E − 0 is also shown with vertical dotted lines. This resonance depends onhω 12 and becomes degenerate with E − 0 in the limithω 12 = 0. For the largest value ofhω 12 , the resonance in Eq. (40c) is also shown for the upper polariton, i.e. 2hω 2 − E + 0 , with a dashed vertical line. This resonance cannot be seen for the other values ofhω 12 : whenhω 12 is too small, the drive energieshω 1 andhω 2 are close to resonance with E − 0 , but very far away from the upper polariton E + 0 . Whenhω 12 matches the polariton splitting E + 0 − E − 0 , as is approximately the case in the panel withhω 12 = 35 meV, the emission linehω 2 − E + 0 is degenerate with the dominating lower polariton E − 0 . When hω 12 is larger, the drives are far away from resonance with the lower polariton, meaning that the dominating polaritonic features (as indicated with vertical gray lines) are diminished, thereby revealing the much weaker peak athω 2 − E + 0 . Fig. 4(b) shows the homodyne squeezing spectrum as 1 + Λ(ω). As discussed in Section III E 3, this spectrum is normalized such that a value of 1 + Λ(ω) = 1 corresponds to no squeezing at all, whereas a value of 1 + Λ(ω) = 0 corresponds perfect squeezing. Since the dominant emission channel in Configuration A is degenerate photon pairs at the lower polariton energy E − 0 , the output field is single-mode squeezed, which is observed as a squeezing spectrum that is minimal (providing the strongest squeezing) at zero frequency. The squeezing bandwidth is several meV, which stems from the typical scale of the resonance linewidths of the cavity photons, excitons and biexcitons in the system [69]. Fig. 4(c) shows the dependence of the squeezing at zero homodyne frequency as a function of the drive frequency differencehω 12 . As expected for Configuration A, the squeezing is strongest in the degenerate limithω 12 = 0, where both drives are resonant with the lower polariton. In addition, there is a weak resonance appearing athω 12 where the drive frequency ω 2 is resonant with the upper polariton. The resonance is weak, because the energy of two upper polaritons overshoot the bound biexciton energy, 2E + 0 E xx b,− . We note that the bichromatic-pump scheme can be utilised to generate strong single-mode squeezing. Since the width of the squeezing dip aroundhω 12 = 0 in Fig. 4(c) is several meV, one can choose a pump frequency difference of e.g.hω 12 = 100 µeV, which allows to filter out the pump lasers spectrally and still have strong single-mode squeezing over a bandwidth that exceeds the resolution of any standard photodetector.
In Fig. 4(d), the leading Fourier components m = 0 and m = 2 of the parametric gain ∆ m is shown as the sum of the absolute values of the spin-diagonal and spin off-diagonal components. The m = −2 component is simply the complex conjugate of m = 2, and the remaining components are negligible in comparison to the ones shown. The reason behind this is twofold: first, the contributions to the parametric gain are products of a † and P † or quantities that are driven by such products, and therefore only contain even Fourier orders. The Fourier components with m larger than 2 are only driven by higher-order processes and are therefore strongly suppressed in comparison with the leading components. This observation corroborates the validity of the DCT perturbative expansion up to third order in the pump field. However, we can not be entirely certain that the contributions to next (5th) order are of a different nature and can introduce new effects. In addition to the total parameric gain, the contribution from the bound biexciton alone is also shown with dashed lines. As can be seen, this constitutes the dominating contribution to the parametric gain, and thus the other contributions in Eq. (15) from Pauli blocking and from the two-exciton scattering continuum are small corrections.

B. Configuration B
Configuration B is defined by setting the cavity frequency such that ) and the mean driving frequency such thathω 1 Fig. 1(d)]. Thus, we expect the driving to be resonant, when the drive frequency differencehω 12 In this case, the drive resonantly excites upper and lower polaritons, which scatter resonantly into bound biexcitons via the Coulomb interaction. The latter process is a non-degenerate polaritonic Feshbach resonance, because the in-scattering particles (upper and lower polaritons) are different with non-degenerate energies. Fig. 5(a) shows the spontaneous four-wave mixing spectra for Configuration B at different drive-frequency differences hω 12 . In contrast to Configuration A, the emission spectrum features two equally bright peaks from the upper and lower polariton, respectively. In terms of the four-wave mixing resonances introduced in Sec. IV A, these peaks correspond to Eq. (40a) withhω 3 = E + 0 and thushω 4 =hω 1 +hω 2 − E + 0 = E − 0 . These two main emission energies are indicated with grey vertical lines in Fig. 5(a). In addition, the weaker four-wave mixing resonance from Eq. (40c) 2hω 2 − E − 0 is indicated with vertical dotted lines. When the drive energy difference becomes larger, the power of the main emission peaks is weakened, whereby the additional four-wave mixing resonances 2hω 2 − E + 0 and 2hω 1 − E + 0 become visible as well (indicated with vertical dashed and dash-dotted lines, respectively).
In Fig. 5(b), the homodyne squeezing spectrum is shown. Here, the most remarkable difference from Configuration A is the appearance of maximal squeezing at the sideband frequency 1 2 (E + 0 − E − 0 ) rather than at zero homodyne frequency. This is because the non-degenerate Feshbach resonance creates two-mode squeezing, i.e. squeezing from the strong correlations between two different frequency bands. This is strongly related to the two equally strong peaks in the spontaneous four-wave mixing spectrum from Fig. 5(a) in contrast to Fig. 4(a), which features a single strong central peak. Although the detection of such high-frequency twomode squeezing is challenging in standard homodyne detection due to finite detection bandwidth, succesful detection can be carried out with bichromatic heterodyne detection, where the high-frequency sideband squeezing is mixed down to a low-frequency beat signal [70,71]. Fig. 5(c) shows the squeezing at zero homodyne frequency (solid line) and at the optimal homodyne sideband frequency (dotted line) as a function of the drive energy differencehω 12 . The optimal sideband frequency is close to 1 2 (E + 0 − E − 0 ), but due to small nonlinear shifts, we have taken the numerically optimal homodyne frequency. The behaviour shows a resonance aroundhω 12 = 1 2 (E + 0 − E − 0 ) (indicated with vertical line), wherehω 1 = E − 0 andhω 2 = E + 0 and polaritons are efficiently excited by the drive fields. This resonance is also seen in Fig. 5(d), where the parametric gain is shown as a function of the drive energy differencehω 12 . Here, it is seen that the parametric gain is strongest aroundhω 12 . Furthermore, as is the case for Configuration A, the bound biexciton (dashed lines) dominates the parametric gain in this configuration.

V. CONCLUSION
In conclusion, we have presented a theoretical investigation of spontaneous four-wave mixing and squeezing in bichromatically pumped atomically-thin semiconductor cavity polaritonic systems, in particular focusing on the strong nonlinear response from the bound biexciton. By applying a rigorous truncation scheme of the many-body state, we have derived a tractable set of equations of motions for exciton and photon fields, as well as the correlated multiparticle fields. In addi- tion, we have employed a Heisenberg-Langevin approach to calculate the fluctuation spectra in the presence of two nondegenerate pump laser fields. The combination of these two methods gives access to the spontaneous four-wave mixing spectra and the squeezing properties of the outcoupled field from the cavity in the presence of strong Coulomb-generated correlations. We have focused on two resonant configurations, corresponding to a degenerate and non-degenerate polaritonic Feshbach resonance, respectively. In the degenerate configuration, a pair of lower polaritons is resonant with the bound biexciton, thereby giving rise to a single dominating peak in the spontaneous four-wave mixing spectrum and strong single-mode squeezing. In the non-degenerate configuration, a upper and lower polariton pair is resonant with the bound biexciton, thereby giving rise to two balanced peaks in the spontaneous four-wave mixing spectrum and strong two-mode squeezing. We believe that these results will open new opportunities in the cross-field between semiconductor physics and nonlinear optics.
Although the numerical calculations in this paper have been performed for an atomically-thin semiconductor, the overall features of the presented phenomena can be expected to be seen in other semiconductor materials with spectrally resolved biexcitons as well. We note that a previous experimental investigation [4] measured parametric gain from biexcitons in a bulk CuCl microcavity in the near-UV spectral range, although no homodyne detection was performed in the experiment. Furthermore, ZnO quantum wells with biexciton binding energies around 15 meV [72] are another interesting platform to potentially observe the predicted squeezing mecha-nism in the near-UV spectrum. Polariton Feshbach resonance has been observed in pump-probe experiments with InGaAs quantum wells [13,29,73].
The matrix elements in Eq. (5) are given bỹ with The biexcitonic wavefunctions Φ ± µ,q are the solutions to the eigenvalue equation where M = m e + m h is the total exciton mass. The numerical details of the solution of this eigenvalue equation is presented in Appendix E.

Appendix B: Formal solutions of multiparticle fluctuations
The first thing to notice is that the coupling coefficients between the correlated fluctuations δ B, δ C, and δ D are linear and do not involve any steady-state expectation values. This means that the formal solution of these fluctuations follows the procedure in Ref. 6 without any changes to the many-body quantities Π and K. The only difference is the source terms involving δ P † , δ a † , δ a in † , δ P in † , which occur together with time-varying amplitudes. To incorporate this into the derivation, we first consider a general equation of motion involving two arbitrary fluctuation operators δ Q and δ R as where A(t) is periodic with the period T = 2π ω 12 . We can then express A(t) in terms of its discrete Fourier series as A(t) = ∑ n A n e −inω 12 t and write δ Q(t) = 1 2π ∞ −∞ dω e −iωt δ Q(ω) and similarly for δ R, such that Multiplying by e iωt and integrating over t, we find Similarly, if we have a product of two amplitudes occurring on the right-hand side, we can use the properties of the Fourier series of a product, Thus a time-domain equation of motion of the form where is the Green's function for δ C(ω) with self-energy The Fourier-transformation of the equation of motion for δ P † , Eq. (11), becomes where ∆ m ζ ζ is defined in Eq. (15), and where the Q-matrix is defined slightly different than in Ref. 6, Inserting the formal solution for δ C into Eq. (B8) and writing all frequency-dependent quantities using the Fourier-index form, we find where the self-energy and renormalized coupling matrices are given by and the renormalized incoupling matrices are given by with K ζ 1 ζ 2 q ζ 1 ζ 0,m (ν) := K ζ 1 ζ 2 q ζ 1 ζ 0 (ν − mω 12 ) and similarly for Q ζ ζ 2 ζ 1 ζ 2 q,m (ν).

Appendix C: Homodyne noise spectrum
In this appendix, we provide additional details of the derivation of the homodyne noise spectrum.
In terms of the fields E lo and E s given in Eq. (26), the mean photocurrent difference from Eq. (30) takes the form where X s (θ ,t) = e iθ E † s (t) + e −iθ E s (t) is the source-field quadrature operator with angle θ = ϕ + π/2. The two terms in θ stem from the local oscillator phase ϕ and a phase displacement π/2 from the beamsplitter.
Similarly, the correlation function of the photocurrent difference from Eq. (32) becomes The last term is the shot noise correlation function, N 0 (t, . In typical homodyne detection setups, the local oscillator is significantly stronger than the signal, which means that the shot noise will be dominated by the local oscillator. In this limit, we can neglect the signal contribution to N 0 (t, τ) which then becomes over t, such that To resolve the integrals over t and t , we now make the assumption that the detector response time τ d is much faster than any time scale in the dynamics of the source field. This corresponds to taking the limit of infinite detection bandwidth. Naturally, the detection bandwidth will restrict the observable homodyne spectrum, but in the interest of characterizing the intrinsic properties of the source field, τ d → 0 is the correct limit to consider. Since the integral over τ in Eq. (C6) only runs over positive values, this means that only the second Heaviside function in Eq. (C9) can be nonzero. Thus, we end up with (for τ > 0) The term in Eq. (C6) proportional to e 2iθ can be connected to S * 2 (ω) by following an analogous derivation. The final expression for the quadrature noise spectrum becomes N(ω) = 1 2 ηF lo sin 2 (ωτ d /2) where the functions Λ i (ω) are defined below Eq. (37) in the main text. In order to characterize the intrinsic noise properties of the source field, we set the detection efficiency η to unity. Thereby, when normalizing the quadrature noise spectrum to the shot-noise level, we obtain Appendix D: Effect of multiparticle fluctuation renormalization As mentioned in the main text, the multiparticle renormalizations of the Heisenberg-Langevin equations of motion only give small quantitative corrections to the spontaneous fourwave mixing and squeezing spectra. This is explicitly demonstrated in Fig. 6, which compares the full and simplified calculations for the same parameters as in Figs. 4 and 5. Specifically, Fig. 6(a)-(b) show the spontaneous four-wave mixing spectrum and squeezing spectrum calculated from the full theory (solid lines) and using the simplifications in Eq. (17) (dots). In Fig. 6(c)-(d), the relative deviation between the two calculations are shown. This relative deviation δ x is calculated as δ x (ω) = x full (ω) − x simplified (ω) /x full (ω), where x(ω) is either S FWM (ω) or 1 + Λ(ω). Evidently, the relative deviation is below 1% for the spontaneous four-wave mixing spectrum, and below 10 −5 for the squeezing spectrum. environment [76] described by the dielectric constant ε e .
The parameters used for the calculations in this paper are summarized in Table I. The phonon-induced dephasing rate γ x is calculated according to the methods given in Ref. [47], without self-consistent inclusion of radiative broadening, because this is contained in the interaction with the quantized electromagnetic field.
In the numerical representation of the equations of motion, the dynamical variables such as a † and P † are expressed in surface-density units a † / √ S and P † / √ S. This ensures that all coupling coefficients and scattering matrices in the equations of motion are independent of the quantization surface area S. Only the input-field driving term contains explicit reference to S, when converting the driving power P in to surface-density units P in → P in /S. For the numerical calculations, S is taken as the laser spot area (see Table I).
In the numerical representation of matrix elements and dynamical variables that depend on wavevectors, we exploit the rotational symmetry and express these quantities in their angle-averaged form as described in detail in Ref. 52. The complete set of exciton wavefunctions is calculated as the solutions to the angle-averaged Wannier equation (see Ref. 6 for details), which is solved in momentum space using a polar coordinate system with 4000 points in the radial direction and 4000 points in the angular direction. As discussed in Sec. III A, only the lowest-energy exciton state is kept, and the higher-lying wave functions are discarded. To construct the matrix elements for the biexcitonic eigenvalue equation (the left-hand-side of of Eq. (A3)) the 2D wavevector-summations in Eq. (A2) must be evaluated for two 2D wave vectors, which presents a significant numerical challenge. Therefore integration in polar coordinates with Gaussian quadratures is used and the matrix elements are calculated for an adjusted, nonequidistant grid of 2D wave vectors, such that the point density is higher for smaller k-values and lower for larger k-values which are required for convergence. This leads to a matrix with a resolution of 100 k-points. Finally, a bicubic spline interpolation is applied to the resulting matrix to form a regular k-grid with 100 × 100 k-space points and the eigenvalues and eigenvectors are calculated.
In the evaluation of some of the matrix elements in Eq. (A1), GPU-based parallelization is used to accelerate the momentum summations using the pytorch library in Python. This makes the evaluation several orders of magnitude faster than a corresponding single-CPU evaluation.
The equations of motion Eq. (5) are solved using the DOP853 adaptive 8th order Runge-Kutta algorithm as implemented in the scipy.integrate library in Python. The matrix inversion required to solve the frequency-space Heisenberg-Langevin equation Eq. (13)