Phase-matching-free parametric oscillators based on two dimensional semiconductors
A. Ciattoni A. Marini C. Rizza C. Conti
physics.optics



alessandro.ciattoni@spin.cnr.it

andrea.marini@icfo.es


^1Consiglio Nazionale delle Ricerche, CNR-SPIN, Via Vetoio 10, 67100 L'Aquila, Italy
^2ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
^3Department of Industrial and Information Engineering and Economics, Via G. Gronchi 18, University of L'Aquila, I-67100 L'Aquila, Italy
^4Institute for Complex Systems (ISC-CNR), Via dei Taurini 19, 00185, Rome, Italy
^5Department of Physics, University Sapienza, Piazzale Aldo Moro 5, 00185, Rome, Italy



Optical parametric oscillators are widely-used pulsed and continuous-wave tunable sources for innumerable applications, as in quantum technologies, imaging and biophysics. A key drawback is material dispersion imposing the phase-matching condition that generally entails a complex setup design, thus hindering tunability and miniaturization. Here we show that the burden of phase-matching is surprisingly absent in parametric micro-resonators adopting monolayer transition-metal dichalcogenides as quadratic nonlinear materials. By the exact solution of nonlinear Maxwell equations and first-principle calculation of the semiconductor nonlinear response, we devise a novel kind of phase-matching-free miniaturized parametric oscillator operating at conventional pump intensities. We find that different two-dimensional semiconductors yield degenerate and non-degenerate emission at various spectral regions thanks to  doubly-resonant mode excitation, which can be tuned through the incidence angle of the external pump laser. In addition we show that high-frequency electrical modulation can be achieved by doping through electrical gating that efficiently shifts the parametric oscillation threshold. Our results pave the way for new ultra-fast tunable micron-sized sources of entangled photons, a key device underpinning any quantum protocol. Highly-miniaturized optical parametric oscillators may also be employed in lab-on-chip technologies for biophysics, environmental pollution detection and security.


Phase-matching-free parametric oscillators based on two dimensional semiconductors
    C. Conti^4,5
    December 30, 2023
==================================================================================

KEYWORDS: two-dimensional materials, nonlinear response, parametric oscillators, sensors, quantum sources, entangled photons, micro-resonators.





§ INTRODUCTION


Optical nonlinearity in photonic materials enables an enormous amount of applications such as frequency conversion <cit.>, all-optical signal processing <cit.>, and non-classical sources <cit.>. Parametric down-conversion (PDC) furnishes tunable sources of coherent radiation <cit.> and generators of entangled photons and squeezed states of light <cit.>. In traditional configurations, a nonlinear crystal with broken centrosymmetry and second-order nonlinearity sustains PDC <cit.>; more recently, effective PDC was reported in centrosymmetric crystals with third-order nonlinearity <cit.> and semiconductor microcavities <cit.>.









Since three-wave parametric coupling is intrinsically weak, one can achieve low oscillation thresholds only by doubly or triply resonant optical cavities. In addition, parametric effects are severely hampered by the destructive interference among the three waves propagating with different wavenumbers k_1,2,3 in the dispersive nonlinear medium because the momentum mismatch Δ k = k_3-k_2-k_1 does not generally vanish (see Fig.1a). To avoid this highly detrimental effect, the use of phase-matching (PM) strategies is imperative. The commonly-adopted birefringence-PM method <cit.> is critically sensible to the nonlinear medium orientation. Quasi-PM <cit.> exploits the momentum due to a manufactured long-scale periodic reversal of the sign of the nonlinear susceptibility and cannot be easily applied in miniaturized system. In semiconductors, PM is achieved by the S-shaped energy-momentum polariton dispersion in the strong coupling of excitons and photons <cit.>, only accessible at low temperatures and large pump angles. Cavity PM <cit.>, also denoted “relaxed” PM <cit.>, occurs in Fabry-Perot microcavities with cavity length ℓ shorter than the coherence length π / Δ k; this technique drastically reduces the effective quadratic susceptibility χ_ eff^(2) (see Fig.1a).
Any of the above mentioned PM techniques entails a non-trivial setup design that is further constrained by the need of resonance operation.

In this manuscript, we show that emerging two-dimensional (2D) materials with high quadratic nonlinearity open unprecedented possibilities for tunable parametric micro-sources. Very remarkably, when illuminated with different visible and infrared waves, these novel 2D materials provide a negligible dispersive dephasing owing to their atomic-scale thickness (i.e ℓ≃ 0, see Fig.1b). Due to the lack of destructive interference, 2D materials support PDC without any need of satisfying a PM condition. Furthermore, these “phase-matching-free” devices turn out to be very versatile and compact, with the additional tunability offered by electrical gating of 2D materials, which provides ultrafast electrical-modulation functionality.



The most famous 2D material, graphene, is not the best candidate for PDC owing to the centrosymmetric structure. In principle, a static external field may break centrosymmetry and induce a χ_ eff^(2), but
the spectrally-flat absorption of graphene is severely detrimental for PDC. Recent years have witnessed the rise of transition metal dichalcogenides (TMDs) as promising photonic 2D materials. TMDs possess several unusual optical properties dependent on the number of layers. Bulk TMDs are semiconductors with an indirect bandgap, but the optical properties of their monolayer (ML) counterpart are characterized by a direct bandgap ranging from ∼ 1.55 eV to ∼ 1.9 eV <cit.> that is beneficial for several optoelectronic applications <cit.>. In addition, ML-TMDs have broken centrosymmetry and thus undergo second-order nonlinear processes <cit.>. Here we study PDC in micro-cavities embedding ML-TMDs; we find that the cavity design is extremely flexible if compared to standard parametric oscillators thanks to their phase-matching-free operation (see Figs.1a,1b). We demonstrate that, at conventional infrared pump intensity, parametric oscillation occurs in wavelength-sized micro-cavities with ML-TMDs. We show that the output signal and idler frequencies can be engineers thanks to the mode selectivity of doubly-resonant cavities; these frequencies are tuned by the pump incidence angle and modulated electrically by an external gate voltage.




§ RESULTS

Two hexagonal lattices of chalcogen atoms embedding a plane of metal atoms arranged at trigonal prismatic sites between the chalcogen neighbors form the structure of ML-TDMs. <cit.> Figure 1c shows the lattice structure of MX_2 ML-TMDs (M = Mo, W, and X = S, Se), and Figs.2a and 2b report the valence and conduction bands of MoS_2 obtained from tight-binding calculations <cit.>. The electronic band structure of other MX_2 materials considered is qualitatively similar. The direct bandgap is about 1.5 eV and implies transparency for infrared radiation; the linear surface conductivity has very small real part (corresponding to absorption) and higher imaginary part
at infrared wavelengths. Figure 2c shows the wavelength dependence of the linear surface conductivities of MX_2. In the presence of an external pump field with angular frequency ω_3, the ML-TMD second-order nonlinear processes lead to down-converted signal and idler waves with angular frequencies ω_1 and ω_2, such that ω_3 =ω_1+ω_2. Figure 2e illustrates the PDC mixing surface conductivities for MoS_2. Both linear and nonlinear conductivities are calculated by a perturbative expansion of the tight-binding Hamiltonian of MX_2 [see Methods and Supplementary Information (SI)]. For infrared photons with energy smaller than the bandgap, extrinsic doping by an externally applied gate voltage (see Fig.1c) modifies the optical properties and leads to an increase of absorption due to free-carrier collisions and to smaller PDC mixing conductivities. Figures 2d and 2f show the dependence of linear and nonlinear surface conductivities on the Fermi level E_ F. As detailed below, extrinsic doping generally leads to a decrease of PDC efficiency.

Figure 1b shows the parametric oscillator design with ML-TMDs. The cavity consists of a dielectric slab (thickness L) surrounded by two Bragg grating mirrors (BGs); the ML-TMD is placed on the left BG inside the cavity. The cavity is illuminated from the left by an incident (i) pump field (frequency ω_3) and the oscillator produces both reflected (r) and transmitted (t) signal and idler fields with frequencies ω_1 = (ω_3 + Δω)/2 and ω_2 = (ω_3 - Δω)/2, where Δω is the beat-note frequency of the parametric oscillation (PO).

As detailed in Methods, the cavity equations for the fields do not contain the momentum mismatch Δ k. Indeed, due to their atomic thickness, ML-TMDs are not optically characterized by a refractive index but rather by a surface conductivity. Hence, the parametric coupling produced by the quadratic surface current of ML-TMDs is not hampered by dispersion and no PM condition is required accordingly.
In order to observe signal and idler generation, only the PO condition is required along with the signal resonance (SR) and idler resonance (IR) conditions leading to a dramatic reduction of the intensity threshold (see Methods). Since there is no PM requirement, such requirements can be met by adjusting either the cavity length L or the pump incidence angle θ as tuning parameters.
For SR and IR, one needs highly reflective mirrors for both signal and idler (see Methods), as obtained by locating stop band of the micron-sized BGs at the half of the pump frequency ω_3/2 <cit.>.
Figure 3 shows the PO analysis for a cavity composed of two BGs with polymethyl methacrylate (PMMA) and MoS_2 deposited on the left mirror. The infrared pump has wavelength λ_3 = 780 nm in the spectral region where the nonlinear properties of MoS_2 are very pronounced (see Fig.2e). The BGs are tuned with their stop bands centered at 1560 nm (=2 λ_3) <cit.>. In Fig.3a we consider the case of normal incidence θ = 0 and we plot the PO (black), SR (red), and IR (green) curves in the (L/λ_3,Δω / ω_3)  plane. Doubly resonant POs (DRPOs) corresponding to the intersection points of these three curves <cit.> are labeled by dashed circles. Therefore, at pump normal incidence, degenerate (Δω =0) and non-degenerate (Δω≠ 0) DRPOs exist at specific cavity lengths. Note that such oscillations also occur for sub-wavelength cavity lengths (L < λ_3). Each oscillation starts when the incident pump intensity I_3^(i) is larger than a threshold I_3Th^(i) (see Methods) <cit.>. Figures 3a2 and 3a3 show the threshold for two specific degenerate and non-denegenerate DRPOs.Panels a2.1 and a3.1 of Fig.3 report the thresholds (black curves on the shadowed vertical planes) corresponding to the PO (black) curves; one can observe that the minimum thresholds occur at SR and IR (identified by the intersection between red and green curves). The minimum intensity thresholds are of the order of GW/cm^2 and the non-degenerate DRPO threshold is greater than the degenerate DRPO one because the reflectivity of the Bragg mirror is maximum at Δω =0 (i.e. at half the pump frequency, as discussed above). In panels a2.2, a2.3 and a2.4 of Fig. 3 (and, seemingly, panels a3.2, a3.3 and a3.4) we report the basic DRPO features by plotting the intensities I_1^(t),I_2^(t),I_3^(t) of the transmitted signal, idler and pump fields as functions of the scaled cavity length L/λ_3 and the incident pump intensity.
Note that, in the considered example, the range of L/λ_3 where the oscillation actually occurs is rather narrow owing to the adopted BGs high reflectivity.

We emphasize that tuning of the PO may be realized by the pump incidence angle θ, which negligibly affects the oscillation thresholds. In Fig.3b and 3c, we analyze the DRPOs by using θ as tuning parameter for a given cavity length. In particular, in Fig.3b we consider a cavity with fixed length as in Fig.3a2.
The PO, SR and IR curves of Fig.3b1 intersect at a degenerate DRPO point at θ≃ 6 deg.

In Fig.3b2 we plot the transmitted signal intensity I_1^(t) as a function of the pump incidence angle and intensity I_3^(i);
one can observe that the intensity threshold is comparable to the case in Fig.3a2 and the range of angles θ where PO occurs is of the order of a hundredth of degree and experimentally feasible.
We show similar results in Figs.3c1 and 3c2, where the non-degenerate DRPO of Fig.3a3 is investigated in a cavity with slightly different length, and achieved at a finite incident angle with unchanged note-beat frequency Δω. A more accurate analysis of Fig.3c1 also reveals that, for a given L, the cavity sustains multiple DRPOs (both degenerate and non-degenerate) at different incidence angles θ.
In Fig.3c3 we plot the transmitted intensity of a degenerate DRPO that grows with the pump intensity above the ignition threshold.

The novel PO with ML-TMDs as nonlinear media are PM-free because of the atomic size of ML-TMDs. The reported several examples of POs with MoS_2 can be also designed by other families of ML-TMDs leading to qualitatively similar results. In the Supplementary Material, we compare the calculated dependence of the pump intensity threshold as function of wavelength λ_3 for parametric oscillators embedding MoS_2, WS_2 and MoSe_2, WSe_2; we find that the chosen material affects the minimal threshold intensity in a given spectral range. One can optimize the choice of the material for a desired spectral content and threshold level.

A further degree of freedom offered by ML-TMDs lies in the electrical tunability through an external gate voltage, as depicted in Fig.1c. The gate voltage increases the Fermi level, and hence affects nonlinearity and absorption because of the electron-electron collision in the conduction band (see Figs.2d,f). Although electrical tunability of MX_2 has not been hitherto experimentally demonstrated, to the best of our knowledge, we emphasize that such a further degree of freedom is absent in traditional parametric oscillators. In the Supplementary Material, we report the pump intensity threshold as a function of the Fermi level of MoS_2, and we show that the threshold may increase by one order of magnitude. The external gate voltage can switch-off PO at fixed optical pump, and fast electrical modulation of the output signal and idler fields can be achieved.







§ CONCLUSIONS


POs can be excited in micron-sized cavities embedding ML-TMDs as nonlinear media at conventional pump intensities in a PM-free regime.
The cavity design remains inherently free of the complexity imposed by the need for PM and may result into doubly resonant PDC of signal and idler waves.
The flexibility offered by such a novel oscillator design enables the engineering of selective degenerate or non-degenerate down-converted excitations by simply modifying the incident angle of the pump field.
Furthermore, the electrical tunability of ML-TMDs can modulate fast output signal and idler waves by bringing POs below threshold.
Based on our calculations, we envisage that novel parametric oscillators embedding ML-TMDs are a new technology for all the applications in which highly-miniaturized tunable source are relevant, including enviromental detection, security, biophyics, imaging and spectroscopy. PM-free ML-TMD microresonators may also potentially boost the realization of micrometric sources of entangled photons when pumped slightly below threshold, thus paving the way for the development of integrated quantum processors.


AM acknowledges useful discussions with F. Javier García de Abajo, Emanuele Distante and Ugo Marzolino.
AC and CR thank the U.S. Army International Technology Center Atlantic for financial support (Grant No. W911NF-14-1-0315).
CC acknowledges funding from the Templeton foundation (Grant number 58277 ) and  PRIN NEMO (reference 2015KEZNYM).


Author Contributions
A.C. and A.M. conceived the idea and worked out the theory. All the authors discussed the results and wrote the paper.




§ METHODS


Parametric down-conversion of MX_2. We calculate the linear and PDC mixing surface conductivities of MX_2
starting from the tight-binding (TB) Hamiltonian of the electronic band structure <cit.>. Since the properties of infrared photons with energies smaller than the bandgap
are determined by small electron momenta around the K and K' valleys, we approximate the full TB Hamiltonian
as a sum of k· p Hamiltonians of first and second order H_0( k,τ,s) <cit.>,
where k is the electron wavenumber and τ and s are the valley and spin indexes, respectively.
We then derive the light-driven electron dynamics through the minimal coupling prescription leading to the
time-dependent Hamiltonian H_0[ k+(e/ħ) A(t),τ,s],
where -e is the electron charge, ħ is the reduced Planck constant, and A(t) is the radiation potential vector, and we obtain Bloch equations
for the interband coherence and the population inversion. Finally, we solve perturbatively the Bloch equations
of ML-TMDs in the weak excitation limit, obtaining the surface current density K(t) after integration over the reciprocal
space

    K(t)     =     Re{∑_j = 1^3  [ σ̂^ L (ω_j)  E_j  e^-iω_j t] + σ̂^(1,2) E_1  E_2  e^-iω_3 t + . 
          .+ σ̂^(1,3) E_1^*  E_3  e^-iω_2 t + σ̂^(2,3) E_2^*  E_3
     e^-iω_1 t},

where σ̂^ L(ω_j) (j=1,2,3) and σ̂^(l,m) (l,m=1,2,3) are the linear and
PDC surface conductivity tensors, respectively. Note that our approach is based on the independent-electron approximation
and is fully justified only for infrared photons far from exciton resonances occurring at photon energies higher than 1.5 eV <cit.>.


Parametric oscillations. The signal, idler and pump fields, labelled with subscripts 1,2,3 respectively, have frequencies ω_n satisfying ω_1+ω_2 = ω_3. By the Transfer Matrix approach, the full electromagnetic analysis of the cavity (see Supplementary Material) yields the equations


    Δ _1 Q_1  + σ̃_23 Q_2^* Q_3     =    0, 
    Δ _2 Q_2  + σ̃_13 Q_1^* Q_3     =    0, 
    Δ _3 Q_3  + σ̃_12 Q_1 Q_2     =    P_3,



where Q_1,Q_2,Q_3 are complex amplitudes proportional to the output fields produced by the pump field which is proportional to the amplitude P_3. Here σ̃_nm are scaled quadratic conductivities of the MX_2 ML-TMD and


    Δ _n     =   σ̃_n  - c /ω_n q_n  ( r_n^( R )  - 1/r_n^( R )  + 1 + r_n^( R ) e^iq_n L  - e^ - iq_n L/r_n^( R ) e^iq_n L  + e^ - iq_n L).


are parameters characterizing the linear cavity where σ̃_n are scaled linear surface conductivities, q_n  = ω_n/c√(ε( ω _n ) - sin ^2 θ) are the longitudinal wavenumbers inside the dielectric slab, ϵ(ω) is the relative permittivity of the dielectric slab, θ is the pump incidence angle whereas r_n^( R ) are the complex reflectivities for right illumination of the left Bragg mirror (with vacuum and the dielectric slab on its left and right sides, respectively). It is worth stressing that the phase-mismatch Δ k = k_3 - k_1 - k_2 does not appear in the basic cavity equations (<ref>). Hence parametric coupling is here not affected by the fields destructive interference and the phase-matching constraint is strictly avoided. Paramateric oscillations (POs) are solutions of Eqs.(<ref>) with Q_1 ≠ 0 and Q_2 ≠ 0 and in this case the compatibility of the first two equations yields (see Supplementary Material)


    | P_3 | ^2 ≥Δ _1 Δ _2^* /σ̃_23σ̃_13^* | Δ _3 |^2,


which is the leading PO condition. As the right hand side of Eq. (<ref>) is generally a complex number, for the PO we have the condition


    arg(Δ _1/σ̃_23) =  arg(Δ _2/σ̃_13).


Eq. (<ref>) can be physically interpreted as a locking of the phase difference Q_1 -  Q_2^* allowing the signal and idler to oscillate. Once Eq.(<ref>) is satisfied, Eq.(<ref>) provides the pump threshold for the onset of PO. Due to the absolute smallness of the nonlinear surface conductivities, in order to have a feasible threshold, the cavity parameters |Δ_n| must be minimized. This can be obtained by choosing the doubly resonant condition for signal and idler corresponding to the minima of |Δ_1| and |Δ_2|, respectively. In order for this minima to be very small, we need that  |r_1^( R )| and |r_2^( R )| are very close to one. One can satisfy such a constraint by a suitable Bragg mirror design to have the stop-band centered at the half of the pump frequency ω_3/2 since, in this case, signal and idler fields experience large mirror reflectance.




Here we provide additional information on technical aspects of the theoretical methods used to model parametric down-conversion and the resulting phase-matching free resonant oscillations within the micro-cavities described in the main paper.



§ PARAMETRIC DOWN-CONVERSION OF MX_2


We calculate the linear and parametric down-conversion (PDC) mixing surface conductivities of monolayer (ML) transition metal dichalcogenides (TMDs) MX_2
(M = Mo, S, and X = S, Se) starting from a tight-binding (TB) description of the electronic band structure of these materials <cit.> and
studying the light-driven electron dynamics by means of Bloch equations for the valence and conduction bands. Note that our approach is based on the
independent-electron approximation and is fully justified at the frequencies considered in the main paper since they are far from exciton resonances
happening at photon energies of ≃ 1.5 eV or higher <cit.>. In addition, for the infrared photon energies
considered, the relevant valence and conduction band regions affecting the infrared response are the ones closer to the band gap around the
K and K' band edges, for which the full TB Hamiltonian can be approximated by a two-band k· p Hamiltonian
H_0( k,τ,s) =  H_1( k,τ,s) +  H_2( k,τ) <cit.>, where

    H_1( k,τ,s) = [ [                   Δ/2 t_0 a (τ k_x - i k_y); t_0 a (τ k_x + i k_y)           τ s Λ - Δ/2 ]], 
           H_2( k,τ,s) = [ [               γ_1a^2k^2 γ_3a^2(τ k_x + i k_y)^2; γ_3a^2(τ k_x - i k_y)^2               γ_2a^2k^2 ]],

and τ,s=± 1 label non-degenerate valleys and spins, while k = (k_x k_y) indicates the electron wave-vector.

The physical parameters of H_0( k,τ,s) are obtained by fitting the k· p valence and conduction energy bands with the ones
obtained from first-principles GW simulations <cit.> accounting for both non-degenerate valleys and spin-orbit coupling, and are listed in Table S1.
We calculate the linear and PDC mixing conductivities of ML-TMDs by introducing the time-dependent Hamiltonian H_0((t),τ,s),
where we have replaced the electron wave-vector with the minimum coupling prescription for the electron quasi-momentum
ħ(t) = ħ k + e  A(t), where e is the electron charge and A(t) is the electromagnetic potential vector
accounting for pump, signal, and idler waves. With this prescription, we define unperturbed and interacting Hamiltonians H_0( k,τ,s) and H_ I( k,τ,s,t), respectively, and write the total Hamiltonian as H_ T( k,τ,s,t) =  H_0[(t),τ,s] =  H_0( k,τ,s) +  H_ I( k,τ,s,t), where

    H_ I( k,τ,s,t)     =    e/ħ [D_x A_x(t) + D_y A_y(t)] + e^2/ħ^2 [D_xx A_x^2(t) + D_xy A_x(t)A_y(t) + D_yy A_y^2(t)],

and the interaction operators are explicitly given by

    D_x     = t_0aτ[ |ψ_ V⟩⟨ψ_ C| + |ψ_ C⟩⟨ψ_ V| ], 
           D_y     = i t_0 a  [ |ψ_ V⟩⟨ψ_ C| - |ψ_ C⟩⟨ψ_ V| ], 
           D_xx  = γ_1a^2 |ψ_ C⟩⟨ψ_ C| + γ_2a^2 |ψ_ V⟩⟨ψ_ V|
    + γ_3a^2 [ |ψ_ V⟩⟨ψ_ C| + |ψ_ C⟩⟨ψ_ V| ], 
           D_yy  = γ_1a^2 |ψ_ C⟩⟨ψ_ C| + γ_2a^2 |ψ_ V⟩⟨ψ_ V|
    - γ_3a^2 [ |ψ_ V⟩⟨ψ_ C| + |ψ_ C⟩⟨ψ_ V| ], 
           D_xy  = 2iγ_3a^2τ[ |ψ_ C⟩⟨ψ_ V| - |ψ_ V⟩⟨ψ_ C| ].

In the expressions above we use the Dirac notation for the conduction |ψ_ C⟩ and valence |ψ_ V⟩ band eigenstates,
and we approximate the matrix elements by their values at the band edges (k = 0). Inserting the Ansatz
|ψ⟩ = c_-|ψ_ V⟩ + c_+ |ψ_ C⟩ in the time-dependent Schrödinger equation
iħ∂_t |ψ⟩ =  H_ T |ψ⟩, and defining the inversion population n_ k=|c_+|^2-|c_-|^2 and the interband
coherence ρ_ k=c_+c_-^*, one gets

    ρ̇_ k    =     -i/ħ (E_ C - E_ V)ρ_ k - γρ_ k + ie/ħ^2 n_ k{ D_x^ CV A_x(t) + D_y^ CV A_y(t) + e/ħ [ D_xx^ CV A_x^2(t) + D_xy^ CV A_x(t)A_y(t) + D_yy^ CV A_y^2(t) ] } + 
            +  ie^2/ħ^3[ (D_xx^ VV - D_xx^ CC)A_x^2(t) + (D_yy^ VV - D_yy^ CC)A_y^2(t) ] ρ_ k, 
    ṅ_ k    =     - 4e/ħ^2 Im{ρ_ k[ D_x^ VCA_x(t) + D_y^ VCA_y(t) + e/ħ [ D_xx^ VCA_x^2(t) + D_xy^ VCA_x(t)A_y(t) + D_yy^ VCA_y^2(t)] ]},

where E_ C ( k) and E_ V ( k) are the conduction and valence energy bands of the unperturbed Hamiltonian H_0,
D^ CV_j=⟨ψ_ C | D_j | ψ_ V⟩ are the interaction matrix elements, and we have introduced a phenomenological relaxation
rate γ = 10 ps^-1 accounting for coherence dephasing <cit.>. In order to obtain the PDC surface conductivities, we consider a coherent superposition
of three monochromatic fields
E(t) =  Re{ E_1 e^-iω_1 t+ E_2 e^-iω_2 t+ E_3 e^-iω_3 t}
with amplitudes E_1, E_2, and E_3, and with angular frequencies ω_1, ω_2, and ω_3, respectively.
In our notation, the field E_3 indicates the external pump field, while E_1, E_2 label the down-converted signal and idler
fields, respectively. The down-converted angular frequencies ω_1 and ω_2 are not independent, but are such that ω_1+ω_2=ω_3
owing to energy conservation. The electromagnetic potential vector related to the coherent superposition of pump, signal, and idler waves is thus given by

    A(t) =  Re{( E_1/iω_1) e^-iω_1 t+( E_2/iω_2) e^-iω_2 t+( E_3/iω_3) e^-iω_3 t}.

We then solve perturbatively the equations above in the vanishing temperature T→ 0 and weak excitation limits such that
n_ k≈ - Θ[ E_ C( k) - E_ F], where Θ(x) indicates the Heaviside step function and E_ F is the Fermi energy. Taking the Ansatz
ρ_ k = ∑_j = ±1,±2,±3ρ_|j|^(j/|j|) e^i(j/|j|)ω_|j| t and disregarding generation of higher harmonics we obtain analytical expressions for the coefficients ρ_|j|^(j/|j|), finding that the macroscopic surface current density given by

    K(t)     =     - e/4π^2ħ∑_τ,s=-1,1∫_-∞^+∞ dk_x ∫_-∞^+∞ dk_y [ ⟨ψ(t)| ∇_ k H_ T (t) | ψ (t) ⟩ - ⟨ψ_ V| ∇_ k H_ T (t) | ψ_ V⟩] = 
        =     - e/2π^2ħ∑_τ,s=-1,1∫_-∞^+∞ dk_x ∫_-∞^+∞ dk_y  Re{ρ_ k (t)[ ∇_ k H_0^ VC + ∇_ k D_x^ VCe/ħ A_x(t) + ∇_ k D_y^ VCe/ħ A_y(t) + . . 
          . . + ∇_ k D_xx^ VCe^2/ħ^2 A_x^2(t) + ∇_ k D_xy^ VCe^2/ħ^2 A_x(t)A_y(t) + ∇_ k D_yy^ VCe^2/ħ^2 A_y^2(t) ] + . 
          . + Θ[ E_ F - E_ C( k) ] [D_xx^ CCe/ħA_x(t)+D_yy^ CCe/ħA_y(t)]},

can be recast into

    K(t) =  Re{∑_j = 1^3  [ σ̂^ L (ω_j)  E_j  e^-iω_j t] + σ̂^(1,2) E_1
     E_2  e^-iω_3 t + σ̂^(1,3) E_1^*  E_3  e^-iω_2 t + σ̂^(2,3) E_2^*  E_3
     e^-iω_1 t},

where σ̂^ L(ω_j) (j=1,2,3) and σ̂^(l,m) (l,m=1,2,3) are the linear and PDC surface conductivity tensors, respectively, and we have neglected again generation of higher harmonics. Since centrosymmetry is broken along the y-direction
in the notation used, the relevant components of the surface conductivity tensors for PDC are the ones such that pump, signal, and idler fields are polarized
along the y-direction, for which

    σ^ L_yy (ω) = ie^2D_yy^ CC/2π^2ħ^2(ω+iγ)∑_τ,s = ± 1∫_-∞^+∞ d k_x ∫_-∞^+∞ d k_y Θ[ E_ F - E_ C( k) ] + 
           + e^2/4iπ^2ħ^2ω∑_τ,s = ± 1∫_-∞^+∞ d k_x ∫_-∞^+∞ d k_y Θ[ E_ C( k) - E_ F] {|D_y^ CV|^2/[(E_ C - E_ V) - ħ (ω + iγ)] + |D_y^ CV|^2/[(E_ C - E_ V) + ħ (ω + iγ)]}, 
          σ^ (1,2)_yyy (ω_1,ω_2,ω_3) = -e^3/4π^2ħ^3ω_1ω_2∑_τ,s = ± 1∑_j = 1^3 ∫_-∞^+∞ d k_x ∫_-∞^+∞ d k_y Θ[ E_ C( k) - E_ F] {D_y^ CVD_yy^ VC/[(E_ C - E_ V) - ħ (ω_j + iγ)] + . 
          . + D_y^ VCD_yy^ CV/[(E_ C - E_ V) + ħ (ω_j + iγ)]}, 
          σ^ (1,3)_yyy (ω_1,ω_2,ω_3) = e^3/4π^2ħ^3ω_1ω_3∑_τ,s = ± 1∫_-∞^+∞ d k_x ∫_-∞^+∞ d k_y Θ[ E_ C( k) - E_ F] {D_y^ CVD_yy^ VC/[(E_ C - E_ V) + ħ (ω_1 - iγ)] + . 
          . + D_y^ VCD_yy^ CV/[(E_ C - E_ V) - ħ (ω_1 - iγ)] + D_y^ CVD_yy^ VC/[(E_ C - E_ V) + ħ (ω_2 + iγ)] + D_y^ VCD_yy^ CV/[(E_ C - E_ V) - ħ (ω_2 + iγ)] + . 
          . + D_y^ CVD_yy^ VC/[(E_ C - E_ V) - ħ (ω_3 + iγ)] + D_y^ VCD_yy^ CV/[(E_ C - E_ V) + ħ (ω_3 + iγ)]}, 
          σ^ (2,3)_yyy (ω_1,ω_2,ω_3) = e^3/4π^2ħ^3ω_2ω_3∑_τ,s = ± 1∫_-∞^+∞ d k_x ∫_-∞^+∞ d k_y Θ[ E_ C( k) - E_ F] {D_y^ CVD_yy^ VC/[(E_ C - E_ V) + ħ (ω_1 + iγ)] + . 
          . + D_y^ VCD_yy^ CV/[(E_ C - E_ V) - ħ (ω_1 + iγ)] + D_y^ CVD_yy^ VC/[(E_ C - E_ V) + ħ (ω_2 - iγ)] + D_y^ VCD_yy^ CV/[(E_ C - E_ V) - ħ (ω_2 - iγ)] + D_y^ CVD_yy^ VC/[(E_ C - E_ V) - ħ (ω_3 + iγ)] + . 
          . + D_y^ VCD_yy^ CV/[(E_ C - E_ V) + ħ (ω_3 + iγ)]}.

Data reported in the main paper are obtained through the expressions above. In what follows, for convenience we will assume the simplified notation σ_n = σ^ L_yy(ω_n) and σ_nm = σ^(n,m)_yyy(ω_1,ω_2,ω_3) since the pump, signal, and idler electric fields are polarized in the y-direction for maximizing PDC within the micro-cavity.





§ EQUATIONS FOR THE OUTPUT FIELDS


In Fig.<ref> we sketch the geometry of the parametric oscillator (PO) considered in our calculations. A dielectric (PMMA) slab of thickness L with a MX_2 monolayer lying on its left side (at z=0) is placed between two Bragg mirrors of thickness d (for convenience we choose the right mirror to be the reflected z → -z copy of the left one). The left side of the cavity is illuminated with an incident (i) pump field which is a monochromatic Transverse Electric (TE) plane wave of frequency ω_3 with incidence angle θ. In addition to the reflected (r) and transmitted (t) pump fields, due to PDC, the cavity also produces (r) and (t) TE plane waves at the frequencies ω_1 (signal) and ω_2 (idler) such that ω_3 = ω_1 + ω_2. It is convenient to set


    ω_1    =   1/2( ω_3 +  Δω), 
    ω_2    =   1/2( ω_3 -  Δω),







since the note-beat frequency Δω = ω_1 - ω_2 is sufficient to label the signal and idler frequencies produced by a pump field of frequency ω_3. Conservation of transverse momentum of the three fields implies that the their complex amplitudes (∼ e^-iω_n t, n=1,2,3) are


    E_n     =    e^iω_n/c xsinθ[ A_ny( z )ê_y ], 
    H_n     =    e^iω_n/c xsinθ√(ε _0 /μ _0 )[ A_nx( z )ê_x  + A_ny( z ) sinθê_z ].


Accordingly the three two-component column vectors (A_nx(z)   A_ny(z))^T fully describe the field and in vacuum, i.e. outside the cavity, they are


    [ A_nx; A_ny;      ] = E_n^( i )[ - cosθ;      1;        ]e^iω_n/c( z + d)cosθ  + E_n^( r )[ cosθ;    1;      ]e^ - iω_n/c( z + d)cosθ,     z <  - d,  
    [ A_nx; A_ny;      ] = E_n^( t )[ - cosθ;      1;        ]e^iω_n/c( z - L - d)cosθ,      z> L+d,


where E_n^( i ), E_n^( r ) and E_n^( t ) are field amplitudes with E_1^( i ) = E_2^( i ) =0. Resorting to the transfer matrix approach, the fields at the left (z=0^-) and right (z=0^+) sides of the MX_2 monolayer are


    [ A_nx; A_ny;      ]_z = 0^ -     =    F_n [ A_nx; A_ny;      ]_z =  - d , 
    [ A_nx; A_ny;      ]_z = 0^ +     =    B_n B'_n [ A_nx; A_ny;      ]_z = L + d,


where F_n, B_n and B'_n are the transfer matrix describing the forward, backward and backward propagations through the left Bragg mirror, the dielectric slab and the right Bragg mirror, respectively. The transfer matrix F_n of the left Bragg mirror is the (ordered) product of the transfer matrices of the slabs composing the mirror. For later convenience it is useful to represent this matrix as <cit.>


    F_n     =   [ c q_n /ω_n cosθ[ t_n^( L ) t_n^( R )  + ( 1 + r_n^( L ))( 1 - r_n^( R ))/2t_n^( R )]  c q_n /ω_n [  - t_n^( L ) t_n^( R )  + ( 1 - r_n^( L ))( 1 - r_n^( R ))/2t_n^( R )];       1/cosθ[  - t_n^( R ) t_n^( L )  + ( 1 + r_n^( L ))( 1 + r_n^( R ))/2t_n^( R )]                [ t_n^( L ) t_n^( R )  + ( 1 - r_n^( R ))( 1 + r_n^( R ))/2t_n^( R )];                                                                                      ],


where q_n  = ω_n/c√(ε( ω _n ) - sin ^2 θ) are the longitudinal wavenumbers inside the dielectric slab, ϵ(ω) is the relative permittivity of the dielectric slab whereas r_n^( L ),t_n^( L ),r_n^( R ),t_n^( R ) are the complex reflectivities r and transmittivities t for left (L) and right (R) illumination of the left Bragg mirror (with vacuum and the dielectric on its left and right sides, respectively). The other relevant transfer matrices are <cit.>


    B_n     =   [             cos( q_n L) ic q_n /ω_n sin( q_n L);   iω_n/c q_nsin( q_n L)             cos( q_n L);                         ], 
    
    B'_n    =   [   1   0;   0 - 1;     ] F_n [   1   0;   0 - 1;     ],


where the last of Eqs. (<ref>) is a consequence of the fact that the right Bragg mirror is the reflected image of the left one. Using Eqs. (<ref>), Eqs. (<ref>) yield


    [ A_nx; A_ny;      ]_z = 0^ -     =   [ V_nx^( i ); V_ny^( i );            ]E_n^( i )  +[ V_nx^( r ); V_ny^( r );            ]E_n^( r ) , 
    
    [ A_nx; A_ny;      ]_z = 0^ +     =   [ V_nx^( t ); V_ny^( t );            ]E_n^( t ),


where


    [ V_nx^( i ); V_ny^( i );            ]   =   1/t_n^( R )[ c q_n /ω_n [  - t_n^( L ) t_n^( R )  + r_n^( L )( r_n^( R )  - 1)];                  t_n^( L ) t_n^( R )  - r_n^( L )( r_n^( R )  + 1);                                                                    ], 
    
    [ V_nx^( r ); V_ny^( r );            ]   =   1/t_n^( R )[ - c q_n /ω_n ( r_n^( R )  - 1);                 r_n^( R )  + 1;                                ], 
    [ V_nx^( t ); V_ny^( t );            ]   =   1/t_n^( R )[ c q_n /ω_n ( r_n^( R ) e^iq_n L  - e^ - iq_n L);            ( r_n^( R ) e^iq_n L  + e^ - iq_n L) ].


The monolayer of MX_2 in the presence of the above TE electromagnetic field hosts a surface current whose harmonic complex amplitudes are K_n = K_n ê_y  where


    K_1   =   [ σ _1 E_1y  + σ _23 E_2y^* E_3y]_z=0  , 
    
     K_2   =   [ σ _2 E_2y  + σ _13 E_1y^* E_3y]_z=0, 
    
     K_3   =   [ σ _3 E_3y  + σ _12 E_1y E_2y]_z=0,


showing both a linear and a quadratic response to the electric field. The effect of such surface current on the field is provided by the electromagnetic boundary conditions at z=0, namely ê_z  ×{[ E_n ]_z = 0^ +    - [ E_n ]_z = 0^ -  } = 0 and ê_z  ×{[ H_n ]_z = 0^ +    - [ H_n ]_z = 0^ -  } = K_n which, using Eqs. (<ref>), can be casted within the two-component column vector description as


    [ A_1x; A_1y;      ]_z = 0^ +    - [ A_1x; A_1y;      ]_z = 0^ -    = [ σ̃_1 A_1y  + σ̃_23 A_2y^* A_3y;                              0;                                ]_z = 0^ +  , 
    [ A_2x; A_2y;      ]_z = 0^ +    - [ A_2x; A_2y;      ]_z = 0^ -    = [ σ̃_2 A_2y  + σ̃_13 A_1y^* A_3y;                              0;                                ]_z = 0^ +  , 
    [ A_3x; A_3y;      ]_z = 0^ +    - [ A_3x; A_3y;      ]_z = 0^ -    = [ σ̃_3 A_3y  + σ̃_12 A_1y A_2y;                            0;                              ]_z = 0^ +  ,


where, for each component, we have set σ̃= √(μ _0 /ε _0 )σ, where ε _0 and μ _0 indicate the dielectric permittivity and magnetic permeability of vacuum, respectively. After inserting Eqs. (<ref>) along with E_1^( i ) = E_2^( i ) =0 into Eqs. (<ref>) we obtain


    V_1x^( t ) E_1^( t )  - V_1x^( r ) E_1^( r )   =   σ̃_1 V_1y^( t ) E_1^( t )  + σ̃_23 V_2y^( t )* V_3y^( t ) E_2^( t )* E_3^( t ) , 
    
     V_1y^( t ) E_1^( t )  - V_1y^( r ) E_1^( r )   =    0, 
    
     V_2x^( t ) E_2^( t )  - V_2x^( r ) E_2^( r )   =   σ̃_2 V_2y^( t ) E_2^( t )  + σ̃_13 V_1y^( t )* V_3y^( t ) E_1^( t )* E_3^( t ) , 
    
     V_2y^( t ) E_2^( t )  - V_2y^( r ) E_2^( r )   =    0, 
    
     V_3x^( t ) E_3^( t )  - V_3x^( i ) E_3^( i )  - V_3x^( r ) E_3^( r )   =   σ̃_3 V_3y^( t ) E_3^( t )  + σ̃_12 V_1y^( t ) V_2y^( t ) E_1^( t ) E_2^( t ) , 
    
     V_3y^( t ) E_3^( t )  - V_3y^( i ) E_3^( i )  - V_3y^( r ) E_3^( r )   =    0,


which are six equations for the six unknown amplitudes E_n^( t ), E_n^( r ) of the signal, idler and pump fields transmitted and reflected by the cavity illuminated by incident pump field of amplitude E_3^( i ). The second, fourth and sixth of Eqs. (<ref>) can be written as


    E_1^( r )   =   V_1y^( t )/V_1y^( r )E_1^( t ), 
    
     E_2^( r )   =   V_2y^( t )/V_2y^( r )E_2^( t ), 
    
     E_3^( r )   =   V_3y^( t )/V_3y^( r )E_3^( t )  - V_3y^( i )/V_3y^( r )E_3^( i ),


showing that the reflected fields can be evaluated once the transmitted fields are known. Substituting the reflected fields from Eqs. (<ref>) into Eqs. (<ref>) and using Eqs. (<ref>) we eventually get


    Δ _1 Q_1  + σ̃_23 Q_2^* Q_3     =    0, 
    Δ _2 Q_2  + σ̃_13 Q_1^* Q_3     =    0, 
    Δ _3 Q_3  + σ̃_12 Q_1 Q_2     =    P_3,


where


    Q_n     =   ( r_n^( R ) e^iq_n L  + e^ - iq_n L/t_n^( R ))E_n^( t ), 
    
     P_3     =   ( 2t_3^( R )cosθ/r_3^( R )  + 1)E_3^( i ), 
    Δ _n     =   σ̃_n  - c q_n /ω_n( r_n^( R )  - 1/r_n^( R )  + 1 + r_n^( R ) e^iq_n L  - e^ - iq_n L/r_n^( R ) e^iq_n L  + e^ - iq_n L).


Equations (<ref>) are the basic equations for the output fields. Note that Eqs. (<ref>) have been derived without resorting to any electromagnetic approximation commonly used in cavity nonlinear optics (e.g. slowly varying amplitude approximation, decoupled approximation for counter-propagating waves, etc.) and this is a consequence of the fact the overall nonlinear response of the MX_2 monolayer is confined to a single plane.





§ DOUBLY RESONANT PARAMETRIC OSCILLATION CONDITIONS


Equations (<ref>) provide the amplitudes Q_n (proportional to the amplitudes of the transmitted signal, idler and pump fields) for a given amplitude P_3 (proportional to the amplitude of the incident pump field). Note that they always admit the solution


    Q_1 = Q_2 = 0,    Q_3  = P_3/Δ_3,


which describes the linear response of the cavity to the pump field without parametric oscillations (POs) in turn characterized by Q_1 ≠ 0 and Q_2 ≠ 0. On the other hand, the first and the complex conjugate of the second of Eqs. (<ref>) are a linear system for Q_1 and Q_2^* and it admits nontrivial solutions only if its determinant vanishes (see Section IV below) or


    | Q_3 |^2  = Δ _1 Δ _2^* /σ̃_23σ̃_13^* .


This condition entails the occurrence of POs since if it is fulfilled, Eqs. (<ref>) have solutions with Q_1 ≠ 0 and Q_2 ≠ 0 which are stable whereas the linear one (Q_1 = Q_2 = 0) becomes unstable. Since the right hand side of Eq. (<ref>) is generally a complex number, it is evident that POs can occur only if such complex number is real and positive or


    Δ _1 Δ _2^* /σ̃_23σ̃_13^* =
    | Δ _1 Δ _2^* /σ̃_23σ̃_13^* |,


which, due to Eqs. (<ref>) and (<ref>), is a constraint joining the cavity length L, the note-beat frequency Δω, and the incident angle θ. Geometrically, Eq. (<ref>) represents a surface Σ of the three-dimensional cavity state space (L,Δω,θ), and its typical slices (θ = 0 or L = L_0) are illustrated in Figs.3a1, 3b1 and 3c1 of the main paper (green curves). At each point of the surface Σ, PO ignites if |Q_3| is sufficiently large to fulfill Eq. (<ref>). Therefore, considering an experiment where the incident pump |P_3| is gradually increased starting from the linear regime where Q_1=Q_2=0, the PO threshold  is obtained by inserting the linear solution Q_3  = P_3/Δ _3 into Eq. (<ref>), thus obtaining


    ( | P_3 |^2 )_th = Δ _1 Δ _2^* /σ̃_23σ̃_13^* | Δ _3 |^2,


which, through the second of Eqs. (<ref>) and the relation I_3^(i) = 1/2√(ϵ_0/μ 0)| E_3^( i )|^2, entails the intensity threshold for the incident pump. Note that the denominator of the right hand side of Eq. (<ref>) contains the nonlinear conductivities σ̃_23 and σ̃_13 whose moduli are so small to generally yield exceedingly large and unfeasible intensity thresholds. A viable way for observing POs thus necessitates the identification of the points (L,Δω,θ) of the surface Σ for which |Δ_1| and |Δ _2| are very close to zero. An inspection of the third of Eqs. (<ref>) reveals that |Δ _n| can not be small if |r_n^( R )| is not close to 1. Therefore, by choosing Bragg mirrors with high reflectivities, the third of Eqs. (<ref>) can be expanded up to the first order in the parameter 1 - | r_n^( R )| ≪ 1 thus getting


    Δ _n  = [ σ̃_n  - i( 2 c q_n /ω_n )sinΓ _n /cosΓ _n  + cos( q_n L)] + {( 2 c q_n /ω_n )1 + ( cosΓ _n  + 2isinΓ _n )cos( q_n L)/[ cosΓ _n  + cos( q_n L)]^2 }( 1 - | r_n^( R )|),


where Γ_n  = q_n L + r_n^( R ). The minima of |Δ_n| are easily seen to occur for Γ_n = m π (where m is any integer) which is exactly the cavity resonance condition for the frequency ω_n <cit.> and where, up to the first order of 1 - | r_n^( R )|,


    Δ _n  = σ̃_n  + [ ( 2 c q_n /ω_n )1/1 + cos(  r_n^( R ))] ( 1 - | r_n^( R )|).


Therefore, as for standard POs based on bulk nonlinear media, the pump intensity threshold is here minimum when one or more of the three fields meet the cavity resonant condition. Designing a Bragg mirror with high reflectivity for both signal and idler fields is relatively simple (see the Section III below) and therefore in this paper we consider only doubly resonant (DR) states where the signal resonance (SR) and idler resonance (IR) conditions


    q_1 L +  r_1^( R )   =    m_1 π, 
    
     q_2 L +  r_2^( R )   =    m_2 π,


are both achieved whereas the pump is non-resonant. Such two equations represents two surfaces Σ_1 and Σ_2 of the space   (L,Δω,θ) whose typical slices are reported in Figs.3a1, 3b1 and 3c1 of the main text (black and red lines respectively) and whose intersection describes the DR cavity states where |Δ_1 Δ_2| is minimum. The (nontrivial) intersection among the three surfaces Σ, Σ_1 and Σ_2 is the set of the DRPO cavity states with feasible intensity threshold. Note that if Σ_1 and Σ_2 intersect at a specific (L,0,θ) point this point also belongs to the surface Σ since for Δω = 0 Eqs. (<ref>) is trivially satisfied since evidently ω_1 = ω_2, Δ_1 = Δ_2 and σ_23 = σ_13. In other words a degenerate (ω_1 = ω_2) DR state always supports a PO which we refer to as a degenerate DRPO. In addition, note that if | Re σ̃_n| ≪ | Im σ̃_n| and | Re σ̃_nm| ≪ | Im σ̃_nm|, if both Eqs. (<ref>) are satisfied with Δω≠ 0, Eq. (<ref>) implies that |  Im( Δ _1 Δ _2^* /σ̃_23σ̃_13^* ) | ≪| Re( Δ _1 Δ _2^* /σ̃_23σ̃_13^* ) | so that, remarkably, if both linear and nonlinear absorption are small the non-degenerate DR states are always very close to PO states. As a consequence the non-degenerate DRPOs with feasible intensity threshold are associated to those points of the surface Σ which are as close as possible to points of the intersection between the surface Σ_1 and Σ_2. Both degenerate and non-degenerate DRPO states are labelled with a dashed disk in Figs.3a1, 3b1 and 3c1 of the main text.









§ BRAGG MIRROR DESIGN


The Bragg mirror is a periodic structure composed of N bi-layers whose dielectric materials have refractive indexes n^(a) and n^(b) and thicknesses a and b. If the layers' thicknesses are chosen to satisfy the Bragg interference condition


    a n^(a) = b n^(b) = π c/2 ω̅,


the mirror has (for normal incidence θ = 0), a spectral stop-band centered at ω̅ whose width is proportional to the refractive index contrast |n^(a) - n^(b)| <cit.>. Within the stop-band the mirror reflectivity is very large, the larger N the closer |r^(R) (ω)| to 1. As explained above in Section II, in order for one of the three fields (pump, signal and idler) to be resonant, it is necessary a very large reflectivity of the Bragg mirror at the field angular frequency, or in other words the frequency ω_n has to lie within the mirror stop-band. As noted above, the degenerate DR states with Δω = 0 (where, from Eqs. (<ref>), ω_1= ω_2 =ω_3/2) rigorously supports POs so that it is convenient to set the center of the mirror stop-bad at ω̅= ω_3/2. Due to the refractive index contrast |n^(a) -n^(b)|, this condition assures that both signal and idler fields experience very large mirror reflectivity in a range of Δω and can accordingly be resonant at the same time. On the other hand the pump frequency ω_3 is twice the central mirror frequency ω̅ and requiring also the pump to resonate would require very large refractive contrast. To avoid this difficulty we have chosen to leave the pump out of resonance.

In the analysis reported in Fig.3 of the main text, we have set as pump wavelength λ_3 = 780   nm. For the Bragg Mirror we have chosen the refractive indexes n_a = 1.2 and n_b = 2.5 so that, in order to have the center of the stop-band at ω̅= ω_3/2 we have chosen the thicknesses a = λ_3/(2 n^(a)) = 325   nm and b = λ_3/(2 n^(b)) = 156   nm. We have also set N=8 for dealing with an efficient, feasible and compact Bragg mirror of length d = N(a+b) =  3848   nm. Using the transfer matrix approach, the complex reflectivity r^(R) (for normal incidence θ=0) of the Bragg mirror which has vacuum and the dielectric at its left and right sides, respectively, is easily evaluated and we plot its absolute value and argument in panel (a) and (b), respectively, of Fig.<ref>. Accordingly, the Bragg mirror stop-band is centered at ω_3/2 and its spectral width is ≃ 0.22 ω_3. As a consequence, if ω_1 and ω_2 lie within this stop-band, signal and idler waves can resonate simultaneously since r_1^(B) = r^(B) (ω_1) and r_2^(B) = r^(B) (ω_2) have moduli very close to 1. The mirror stop-band width therefore yields the note-beat frequency range 0 ≤Δω < 0.22 ω, which is the one considered in the analysis reported in Fig.3 of the main text. Note that ω_3 lies outside the mirror stop-band and thus the pump field does not resonate.





§ OUTPUT FIELDS


In order to evaluate the PO output fields E_n^( t ), Eqs. (<ref>) have to be solved for a given incident pump field E_3^( i). POs are characterized by Q_1 ≠ 0 and Q_2 ≠ 0 for a given Q_3 ≠ 0. First note that Eqs. (<ref>) are left invariant by the gauge transformation


    Q_1     →    Q_1 e^i θ, 
    
     Q_2     →    Q_2 e^ - i θ,


and this implies that for a given P_3 there are infinite pairs (Q_1,Q_2), all with the same Ψ =  Q_1 +  Q_2. In other words, the phase difference Φ =  Q_1 -  Q_2 is not set by the pump field P_3. Evidently, such symmetry is spontaneously broken in actual experiments where a single pair (Q_1,Q_2) (i.e. a single value of Φ) is selected by the specific way chosen to trigger POs.

In the case of POs, the first and the complex conjugate of the second of Eqs. (<ref>) can be casted as


    Q_1 /Q_2^*  = - σ̃_23 Q_3/Δ _1 , 
    Q_1 /Q_2^*  =  - Δ _2^* /σ̃_13^* Q_3^* ,


whose consistency requires their right and left hand sides to coincide or


    | Q_3 |^2  = Δ _1 Δ _2^* /σ̃_23σ̃_13^*


which is the PO oscillation condition of Section II. [see Eq. (<ref>)]. The right hand side of Eq. (<ref>) is a positive real number if and only if the complex numbers Δ _1 /σ̃_23 and Δ _2 /σ̃_13 have the same argument φ or


    Δ _1     =   σ̃_23| Δ _1 /σ̃_23|e^iφ, 
    Δ _2     =   σ̃_13| Δ _2 /σ̃_13|e^iφ,


which are equivalent to Eq. (<ref>) of Section II so that, considering only those states for which Eqs. (<ref>) are satisfied, Eqs. (<ref>) yield


    | Δ _1 /σ̃_23|e^iφ Q_1  + Q_2^* Q_3     =    0, 
    | Δ _2 /σ̃_13|e^iφ Q_2  + Q_1^* Q_3     =    0, 
    Δ _3 Q_3  + σ̃_12 Q_1 Q_2     =    P_3.


To solve this equation, after noting that Eqs. (<ref>) require that |Q_1 /Q_2 |^2  = |Δ _2σ̃_23/Δ _1 σ̃_13| and exploiting the above discussed gauge symmetry, we set


    Q_1     =   √(| Δ _2 /σ̃_13|)| Q |e^i1/2( Ψ  + Φ), 
    
    Q_2     =   √(| Δ _1 /σ̃_23|)| Q |e^i1/2( Ψ  - Φ),


where we have used the symbol |Q| to stress than this quantity is real. Inserting Eqs. (<ref>) into Eqs. (<ref>) yield


    Q_3     =    - √(| Δ _1 Δ _2 /σ̃_23σ̃_13|) e^i( Ψ  + φ) , 
    | Q |^2  +  Δ _3 /σ̃_12√(| σ̃_23σ̃_13/Δ _1 Δ _2 |) e^ - iΨ Q_3     =   √(| σ̃_23σ̃_13/Δ _1 Δ _2 |)P_3 /σ̃_12e^ - iΨ.


Note that the first two of Eqs. (<ref>) both reduces to the first of Eqs. (<ref>) through the change of variables given by Eqs. (<ref>) and this is a consequence of the necessary Eq. (<ref>). Substituting Q_3 from the first of Eqs. (<ref>) into the second we get


    | Q |^2  - Δ _3 /σ̃_12e^iφ  =
    √(| σ̃_23σ̃_13/Δ _1 Δ _2 |)P_3 /σ̃_12
    e^ - iΨ,


which is a single complex equation for |Q| and Ψ. After equating the square-moduli of the left and right sides of this equation we get


    | Q |^4  - 2 Re( Δ _3 /σ̃_12e^iφ)| Q |^2  + ( | Δ _3 /σ̃_12e^iφ|^2  -
    | σ̃_23σ̃_13/Δ _1 Δ _2 |
    | P_3 /σ̃_12|^2 ) = 0,


which is a biquadratic equation for | Q | whose solutions are


    | Q | = √( Re( Δ _3 /σ̃_12e^iφ) +ξ√(| σ̃_23σ̃_13/Δ _1 Δ _2 |
    | P_3 /σ̃_12|^2
     -  Im^2 ( Δ _3 /σ̃_12e^iφ) )),


where ξ = ± 1. Note that, due to the ξ factor, generally there are two |Q| corresponding to a |P_3| and hence bistable POs can in principle occur. In addition, it is fundamental stressing that |Q| is real and hence Eq. (<ref>) provides its value only if the arguments of the square roots are positive. Before discussing the range of |P_3| where this is the case (see below), we assume |Q| real and we deduce the output field amplitudes. Equation (<ref>) yields


    e^iΨ  = √(| σ̃_23σ̃_13/Δ _1 Δ _2 |)P_3 /σ̃_12( | Q |^2  - Δ _3 /σ̃_12e^iφ),


which is consistent since the modulus of its right hand side, due to Eq. (<ref>), is equal to 1. Hence Eqs. (<ref>) and the first of Eqs. (<ref>) eventually yield


    Q_1     =   ζ| Δ _2 σ̃_23/Δ _1 σ̃_13|^1/4√(P_3 /σ̃_12( | Q |^2  - Δ _3 /σ̃_12e^iφ)) e^iΦ/2| Q |, 
    
     Q_2     =   ζ| Δ _1 σ̃_13/Δ _2 σ̃_23|^1/4√(P_3 /σ̃_12( | Q |^2  - Δ _3 /σ̃_12e^iφ)) e^ - iΦ/2| Q |, 
    
     Q_3     =   P_3 /σ̃_12( | Q |^2  - Δ _3 /σ̃_12e^iφ)e^i( φ  + π),


where ζ = ± 1 and the principal branch is assumed for all the complex square-roots.

Note that the output fields of Eqs. (<ref>) satisfy Eq. (<ref>) so that they describe all the possible cavity POs whenever they exist or, in other words, whenever |Q| of Eq. (<ref>) is a positive real number. Such requirement evidently sets a range for the input pump intensity |P_3|^2 and there are four different cases corresponding to the two values of ξ and of the two signs of
Re( Δ _3 /σ̃_12e^iφ). The results of this analysis are reported in the following table.


    [                                                                          Re( Δ _3 /σ̃_12e^iφ)>0                         Re( Δ _3 /σ̃_12e^iφ) <0;                                                                                                                                                ;                                             ξ=1                    | P_3 |^2  > ( | P_3 |^2 )_-                   | P_3 |^2  > ( | P_3 |^2 )_th;                                                                                                                                                ;                                                                                                                                                ;                                            ξ=-1 ( | P_3 |^2 )_- < | P_3 |^2  < ( | P_3 |^2 )_th                                     no  |P_3|^2;                                                                                                                                                ;                                                 ]


Here we have set


    ( | P_3 |^2 )_- = | Δ _1 Δ _2 /σ̃_23σ̃_13| Im ^2 ( Δ _3 | σ̃_12|/σ̃_12e^iφ) <
    | Δ _1 Δ _2 /σ̃_23σ̃_13|| Δ _3 |^2  =
    ( | P_3 |^2 )_th.


Therefore at each state where PO can occur (i.e. at a each point of the surface Σ) the scenario is the following one. If Re( Δ _3 /σ̃_12e^iφ) < 0, there is only one PO (with ξ =1) that effectively starts when Eq. (<ref>) is satisfied, thus confirming the analysis of Section III. On the other hand, if Re( Δ _3 /σ̃_12e^iφ) > 0 the scenario changes qualitatively since in this case there are two allowed POs (with ξ=1 and ξ=-1) when  ( | P_3 |^2 )_- < | P_3 |^2  < ( | P_3 |^2 )_th  and a single PO (with ξ=1) when Eq. (<ref>) is satisfied. As a consequence in this case POs also exist below the threshold.

In order to grasp the reason why the sub-threshold POs have not been entailed in Section III, note that in the case Re( Δ _3 /σ̃_12e^iφ) > 0, if | P_3 |^2 = ( | P_3 |^2 )_-, Eq. (<ref>) implies that | Q | = √( Re( Δ _3 /σ̃_12e^iφ)) 0 so that Q_1 ≠ 0 and Q_2 ≠ 0. In other words in this situation the signal and idler fields do not vanish at the threshold and accordingly this case is ruled out from the reasoning of Section III where the threshold has been obtained for PO oscillation starting from the linear regime. In a realistic experiment PO is switched on starting from the linear regime, with the intensity threshold given by Eq. (<ref>). However, once PO ignites, by changing the pump intensity, incidence angle θ, or the cavity length, we argue that one can in principle access sub-threshold PO states. In every design considered in the main paper we have focused on the case Re( Δ _3 /σ̃_12e^iφ) < 0, where sub-threshold PO does not occur.





§ PUMP INTENSITY THRESHOLDS

In Fig. <ref> we compare the calculated pump intensity thresholds versus the pump wavelength λ_3 for parametric oscillators embedding MoS_2, WS_2 (Fig. <ref>a) and MoSe_2, WSe_2 (Fig. <ref>b). Note that, while the minimal pump intensity threshold occurs at λ_3 ≈ 780 nm for MoS_2 and WS_2, it shifts to λ_3 ≈ 940 nm for MoSe_2 and WSe_2. None of the ML-TMDs examined enables feasible PO with low pump intensity threshold at optical frequencies owing to the enlarged absorption in this frequency range, which is the main responsible for oscillation quenching. In addition, at optical frequencies such materials exhibit exciton resonances <cit.> (not taken into account in our theoretical approach) that are also detrimental for POs owing to the enhanced absorption they are accompanied with. In Fig. <ref>c we plot the pump intensity threshold as a function of the Fermi level of MoS_2, showing that it can be increased efficiently. Thus, the external gate voltage quenches POs when the optical pump is fixed and fast modulation of the output signal and idler fields can be achieved with novel parametric oscillators embedding ML-TMDs.



10


Franken1961 P. A. Franken, A. E. Hill, C. W. Peters, and G. Weinreich, “Generation of Optical Harmonics,´´ Phys. Rev. Lett. 7, 118 - 119 (1961).

Birks2000 T. A. Birks, W. J. Wadsworth, P. St.J. Russell, “Supercontinuum generation in tapered fibers,´´ Opt. Lett. 25, 1415 - 1417 (2000).

Stegeman1996 G. I. Stegeman, D. J. Hagan, and L. Torner, “χ^(2) cascading phenomena and their applications to all-optical signal processing, mode-locking, pulse compression and solitons,´´ Opt. Quant. Electron. 28, 1691 - 1740 (1996).

Koos2009 C. Koos, P. Vorreau, T. Vallaitis, P. Dumon, W. Bogaerts, R. Baets, B. Esembeson, I. Biaggio, T. Michinobu, F. Diederich, W. Freude, and J. Leuthold, “All-optical high-speed signal processing with silicon–organic hybrid slot waveguides,´´ Nat. Photon. 3, 216 - 219 (2009).

Kwiat1995 P. G. Kwiat, K. Mattle, H. Weinfurter, A. Zeilinger, A. V. Sergienko, and Y. Shih, “New High-Intensity Source of Polarization-Entangled Photon Pairs,´´ Phys. Rev. Lett. 75, 4337 - 4341  (1995).

Stevenson2006 R. M. Stevenson, R. J. Young, P. Atkinson, K. Cooper, D. A. Ritchie, and A. J. Shields, “A semiconductor source of triggered entangled photon pairs,´´ Nature 439, 179 - 182 (2006).

Giordmaine1965  J. A. Giordmaine and R. C. Miller, “Tunable Coherent Parametric Oscillation in LiNbO3 at Optical Frequencies,´´ Phys. Rev. Lett. 14, 973 - 976 (1965).

Yariv1966 A. Yariv and W. H, Louisell, “Theory of the Optical Parametric Oscillator,´´ IEEE J. Quant. Electron. 9, 418 - 424 (1966).

Brosnan1979 S. J. Brosnan and R. L. Byer, “Optical Parametric Oscillator Threshold and Linewidth Studies,´´ IEEE J. Quant. Electron. 15, 415 - 431 (1979).

Eckardt1991 R. C. Eckardt, C. D. Nabors, W. J. Kozlovsky, and R. L. Byer, “Optical parametric oscillator frequency tuning and control,´´ J. Opt. Soc. Am. B 8, 646 - 667 (1991).

Debuisschert1993 T. Debuisschert, A. Sizmann, E. Giacobino, and C. Fabre, “Type-II continuous-wave optical parametric oscillators: oscillation and frequency-tuning characteristics,´´ J. Opt. Soc. Am. B 10, 1668 - 1680 (1993).

Fabre1997 C. Fabre, P. F. Cohadon, and C. Schwob, “CW optical parametric oscillators: single mode operation
and frequency tuning properties,´´ Quantum Semiclass. Opt. 9, 165 - 172 (1997).

Levy2010 J. S. Levy, A. Gondarenko, M. A. Foster, A. C. Turner-Foster, A. L. Gaeta, and M. Lipson, “CMOS-compatible multiple-wavelength oscillator for on-chip optical interconnects,´´ Nat. Photon. 4, 37 - 40 (2010).

Razzari2010 L. Razzari, D. Duchesne, M. Ferrera, R. Morandotti, S. Chu, B. E. Little, and D. J. Moss, “CMOS-compatible integrated optical hyperparametric oscillator,´´ Nat. Photon. 4, 41 - 45 (2010).

Wu1987 L.-A. Wu, M. Xiao, and H. J. Kimble, “Squeezed states of light from an optical parametric oscillator,´´ J. Opt. Soc. Am. B 4, 1465 - 1475 (1987).

Lu2000 Y. J. Lu and Z. Y. Ou, “Optical parametric oscillator far below threshold: Experiment versus theory,´´ Phys. Rev. A 62, 033804 (2000).

Ciuti2003 C. Ciuti, P. Schwendimann and A. Quattropani, “Theory of polariton parametric interactions in semiconductor microcavities,´´ Semicond. Sci. Technol. 18, S279 – S293 (2003).

Diederichs2006 C. Diederichs, J. Tignon, G. Dasbach, C. Ciuti, A. Lemaitre, J. Bloch, Ph. Roussignol and C. Delalande, “Parametric oscillation in vertical triple microcavities´´ Nature 440, 904 - 907 (2006).

Abbarchi2011 M. Abbarchi, V. Ardizzone, T. Lecomte, A. Lemaıtre, I. Sagnes, P. Senellart, J. Bloch, P. Roussignol, and J. Tignon, “One-dimensional microcavity-based optical parametric oscillator: Generation of balanced twin beams in strong and weak coupling regime,´´ Phys. Rev. B 83, 201310(R) (2011).

Giordmaine1962 J. A. Giordmaine, “Mixing of Light Beams in Crystals,´´ Phys. Rev. Lett. 8, 19 – 20 (1962).

Vodopyanov2004 K. L. Vodopyanov, O. Levi, P. S. Kuo, T. J. Pinguet, J. S. Harris, M. M. Fejer, B. Gerard, L. Becouarn, and E. Lallier, “Optical parametric oscillation in quasi-phase-matched GaAs,´´ Opt. Lett. 29, 1912 – 1914 (2004).

Canalias2007 C. Canalias, and V. Pasiskevicius, “Mirrorless optical parametric oscillator,´´ Nature Photonics 1 459 - 462 (2007).

Savvidis2000 P. G. Savvidis, J. J. Baumberg, R. M. Stevenson, M. S. Skolnick, D. M. Whittaker, and J. S. Roberts, “Angle-Resonant Stimulated Polariton Amplifier,´´ Phys. Rev. Lett. 84, 1547 – 1550 (2000).

Ciuti2000 C. Ciuti, P. Schwendimann, B. Deveaud, and A. Quattropani, ´´Theory of the angle-resonant polariton amplifier,´´ Phys. Rev. B 62, R4825–-R4828 (2000).

Xie2011 Z. D. Xie, X. J. Lv, Y.H. Liu, W. Ling, Z. L. Wang, Y. X. Fan, and S. N. Zhu, ´´Cavity Phase Matching via an Optical Parametric Oscillator Consisting of a Dielectric Nonlinear Crystal Sheet,´´ Phys. Rev. Lett. 106, 083901 (2011).

Clément2015 Q. Clément, J.-M. Melkonian, M. Raybaut, J.-B. Dherbecourt, A. Godard, B. Boulanger, and M. Lefebvre, ´´Ultrawidely tunable optical parametric oscillators based on relaxed phase matching: theoretical analysis,´´ J. Opt. Soc. Am. B 32, 52 - 68 (2015).

Mak2010 K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz,  “Atomically Thin MoS_2: A New Direct-Gap Semiconductor,´´ Phys. Rev. Lett. 105, 136805 (2010).

Splendiani2010 A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C. Y. Chim, G. Galli, and F. Wang, “Emerging photoluminescence in monolayer MoS_2,´´Nano Lett. 10 1271 - 1275 (2010).

Wang2012 Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman, and M. S. Strano, “Electronics and optoelectronics of two-dimensional transition metal dichalcogenides,´´ Nat. Nanotechnol. 7, 699 - 712 (2012).

Sun2016 Z. Sun, A. Martinez, and F. Wang, “Optical modulators with 2D layered materials,´´ Nat. Photon. 10, 227 - 238 (2016).

Kumar2013 N. Kumar, S. Najmaei, Q. Cui, F. Ceballos, P. M. Ajayan, J. Lou, and H. Zhao, “Second harmonic microscopy of monolayer MoS_2,´´Phys. Rev. B 87, 161403(R) (2013).

Li2013 Y. Li, Y. Rao, K. F. Mak, Y. You, S. Wang, C. R. Dean, T. F. Heinz, “Probing symmetry properties of few-layer MoS_2 and h-BN by
optical second-harmonic generation,´´ Nano Lett. 13, 3329 - 3333 (2013).

Malard2013 L. M. Malard, T. V. Alencar, A. P. M. Barboza, K. F. Mak, and A. M. de Paula, “Observation of intense second harmonic generation
from MoS_2 atomic crystals,´´ Phys. Rev. B 87, 201401(R) (2013).

Janisch2014 C. Janisch, Y. Wang, D. Ma, N. Mehta, A. L. Elias, N. Perea-Lopez, M. Terrones, V. Crespi, and Z. Liu,
Sci. Rep. 4, 5530 (2014).

Le2016 C. T. Le, D. J. Clark, F. Ullah, V. Senthilkumar, J. I. Jang, Y. Sim, M.-J. Seong, K.-H. Chung, H. Park, and Y. S. Kim, “Nonlinear optical characteristics of monolayer MoSe_2,´´ Ann. Phys. 528, 551 - 559 (2016).

LWY2013  G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao, “Three-band tight-binding model for monolayers of group-VIB transition
metal dichalcogenides,´´ Physical Review B 88, 085433 (2013).

EPAPS See Supplemental Material for more details on the technical aspects of the theory.

Ugeda2014 M. M. Ugeda, A. J. Bradley, S.-F. Shi, F. H. da Jornada, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, F. Wang, S. G. Louie, and M. F. Crommie, “Giant bandgap renormalization and excitonic effects in a monolayer transition metal dichalcogenide semiconductor´´, Nat. Mater. 13, 1091-1095 (2014).

Selig2016 M. Selig, G. Berghäuser, A. Raja, P. Nagler, C. Schüller, T. F. Heinz, T. Korn, A. Chernikov, E. Malic, and A. Knorr, “Excitonic linewidth and coherence lifetime in monolayer transition metal dichalcogenides´´, Nat. Commun 7, 13279 (2016).
LWY2013  G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao, “Three-band tight-binding model for monolayers of group-VIB transition
metal dichalcogenides,´´ Physical Review B 88, 085433 (2013).
UBS2014 M. M. Ugeda, A. J. Bradley, S.-F. Shi, F. H. da Jornada, Y. Zhang, D. Y. Qiu, W. Ruan, S.-K. Mo, Z. Hussain, Z.-X. Shen, F. Wang, S. G. Louie,
and M. F. Crommie, “Giant bandgap renormalization and excitonic effects in a monolayer transition metal dichalcogenide semiconductor,´´
Nat. Mater. 13, 1091 - 1095 (2014).
CCB2016 F. Ceballos,  Q. Cui, M. Z. Bellusa, and H. Zhao, “Exciton formation in monolayer transition metal dichalcogenides,´´ Nanoscale 8,
11681 - 11688 (2016).
YYC2014 X. Yin, Z. Ye, D. A. Chenet, Y. Ye, K. O'Brien, J. C. Hone, X. Zhang, “Edge Nonlinear Optics on a MoS2 Atomic Monolayer´´, Science 344, 488–490 (2014).
Kavok A. V. Kavokin, J. J. Baumberg, G. Malpuech and F. P. Laussy, Microcavities (Oxford University Press, USA 2008)


