Purcell enhancement of the parametric down-conversion in two-dimensional nonlinear materials

Ultracompact nonlinear optical devices utilizing two-dimensional (2D) materials and nanostructures are emerging as important elements of photonic circuits. Integration of the nonlinear material into a subwavelength cavity or waveguide leads to a strong Purcell enhancement of the nonlinear processes and compensates for a small interaction volume. The generic feature of such devices which makes them especially challenging for analysis is strong dissipation of both the nonlinear polarization and highly confined modes of a subwavelength cavity. Here we solve a quantum-electrodynamic problem of the spontaneous and stimulated parametric down-conversion in a nonlinear quasi-2D waveguide or cavity. We develop a rigorous Heisenberg-Langevin approach which includes dissipation and fluctuations in the electron ensemble and in the electromagnetic field of a cavity on equal footing. Within a relatively simple model, we take into account the nonlinear coupling of the quantized cavity modes, their interaction with a dissipative reservoir and the outside world, amplification of thermal noise and zero-point fluctuations of the electromagnetic field, and other relevant effects. We derive closed-form analytic results for relevant quantities such as the spontaneous parametric signal power and the threshold for parametric instability. We find a strong reduction in the parametric instability threshold for 2D nonlinear materials in a subwavelength cavity and provide a comparison with conventional nonlinear photonic devices.


I. INTRODUCTION
Enhancement of the radiative processes due to the localization of emitters in a subwavelength cavity (so-called Purcell enhancement [1]) is a fundamental cavity quantum electrodynamics (QED) effect of great importance for numerous applications. The bulk of the research has been focused on exploring the enhancement of spontaneous emission in various compact radiation sources from single quantum emitters to LEDs and nanolasers.
The nonlinear optics has received relatively less attention; however, recent advancements in strong light localization using subwavelength cavities, photonic crystals, metamaterials, and metasurfaces enabled the nonlinear optics in ultrasmall volumes and at relatively low power levels; see e.g. [2][3][4][5][6][7][8][9][10] and references therein. The rise of 2D materials with atomic monolayer thickness and excellent nonlinear optical properties, such as graphene [11][12][13][14] and transition metal dichalcogenide monolayers [15,16] has enabled quasi-2D cavities and waveguides only a few nm thick [17,18]. These advances create new exciting opportunities for ultracompact nonlinear optical devices, but also raise important issues of the correct description of quantum fields in systems with strong dissipation both in a macroscopic ensemble of fermionic emitters (e.g. a 2D semiconductor or monolayer graphene, or a 2D electron gas in a quantum well) and for the electromagnetic (EM) field in a cavity.
One important application for Purcell-enhanced nonlinear optics is compact systems for generation of squeezed and entangled photon states as a result of parametric downconversion. Such systems are inevitably lossy. The general approach to introducing dissipation and corresponding fluctuations has been known for a long time and is based on the Heisenberg-Langevin formalism; e.g. [19][20][21][22][23]. Its generalization to systems far from equilibrium, with arbitrary dissipation mechanisms and arbitrary photon density of states is nontrivial; see e.g. [24,25]. Recent work [10] considered the process of spontaneous parametric down-conversion in hyperbolic metamaterials, in which the EM field dissipation and fluctuations are due to an equilibrium thermal reservoir. In the present paper we consider both spontaneous and stimulated parametric down-conversion in a generic quasi-2D subwavelength cavity, taking into account dissipation and fluctuations both due to absorption in the intracavity material and due to in/outcoupling of the intracavity EM field with the outside world.
We generalize the properties of Langevin noise sources known for a single mode of a quan-tized field (e.g. [19,20,26]) to an ensemble of coupled field oscillators. We derive the properties of the Langevin sources needed to conserve their commutation relations and show that they are not affected by a more complicated dynamics of coupled Heisenberg field operators; moreover, this statement does not depend on any specific microscopic model of a dissipative reservoir. We are able to derive closed-form analytic results for the spontaneous parametric signal, the parametric gain, and the threshold for parametric amplification. These expressions include the contributions from all relevant dissipation and fluctuation effects such as absorption and radiation losses, interaction with thermal and zero-point fluctuations, parametric amplification of thermal noise and seed photons at the signal frequency, etc.; see e.g.
Our approach has obvious limitations of a Heisenberg-Langevin formalism, namely it assumes that the coupling of a dynamic subsystem to a dissipative reservoir is sufficiently weak. If this is not the case and the coupling to other EM modes, photons, etc. is strong, one would have to include it as part of an "exact" Hamiltonian dynamics, in which case there would be no need in adding the corresponding Langevin sources and the commutation relations would be satisfied automatically. We also do not investigate the nonlinear stage of parametric oscillations accompanied by the pump depletion, nonlinear evolution of phasematching conditions, nonlinear modification of refraction and diffraction losses, and other nonlinear effects that are essentially classical and depend on a particular experimental setup.
Section II describes the spatial structure of the EM field in a subwavelength quasi-2D electrodynamic structure, develops the quantization procedure in a dissipationless system, and discusses three-wave mode matching conditions. Section III introduces the Heisenberg-Langevin approach for the parametric down-conversion in a dissipative cavity. It derives convenient analytic expressions for the spontaneous parametric signal, the parametric amplification threshold in plane-parallel cavities, and the signal evolution at the linear stage.
We discuss several numerical examples for the parametric down-conversion in quasi-2D systems studied by other groups. Section IV compares parametric amplification threshold in a subwavelength cavity with the one in a standard Fabry-Perot cavity containing a 2D nonlinear layer. In this case the performance tradeoff is between the cavity losses and the modal overlap with a nonlinear layer. Larger cavities tend to have a higher Q-factor but lower coupling to a nonlinear 2D layer. Our results show that it is possible to achieve a significant reduction of the parametric amplification threshold due to Purcell enhancement in quasi-2D subwavelength cavities.

II. PARAMETRIC DOWN-CONVERSION IN A CONSERVATIVE SYSTEM
Consider three cavity modes with frequencies related by the energy conservation in the parametric down-conversion process: Here the pumping at frequency ω p will be considered a classical coherent field, The field at signal and idler frequencies, ω s and ω i , will be the quantum field described by whereĉ ν andĉ ν † are boson annihilation and creation operators. The functions E p,s,i (r) in Eqs. (2) and (3) are determined by the spatial structure of the cavity modes. The normalization of functions E ν (r) needs to be chosen in such a way that the commutation relation for boson operatorsĉ ν andĉ ν † have a standard form [ĉ ν ,ĉ ν † ] = δ νν . Following [27][28][29], one can obtain where ω ν is the eigenfrequency of a cavity mode, E νj (r) is a Cartesian component of the vector field E ν (r), ε jk (ω, r) is the dielectric tensor, and V is a cavity volume (a quantization volume).
Equation (4) is valid when the dissipation is weak enough. Specifically, the following three conditions have to be satisfied for a dissipation rate Γ of a given cavity mode. The first condition is obvious: Γ ω has to be true for the frequencies of all modes involved in the parametric process. The second condition implies that the change of the Hermitian dielectric function ε jk (ω) has to be small over the frequency interval of the order of Γ : The third condition states that the change in the derivative of ε jk (ω) which enters the expression for the EM energy density in Eq. (4) must also be small: Consider a 3D cavity filled with an isotropic dielectric medium, as sketched in Fig. 1.
The cavity thickness in z-direction is much smaller than wavelength: L z c/ √ε ω , wherē ε is a typical (average) value of the dielectric constant of the filling. As was shown in [26], if the dielectric filling consists of plane-parallel layers, i.e. ε = ε(z), the structure of the cavity eigenmodes is quasi-electrostatic along the z-axis, i.e. E z ε(z) ≈ const, E x,y E z . In this case the field of a cavity mode can be written as where the constants D p,s,i are coordinate-independent amplitudes of the electric induction.
To find the functions ζ p,s,i (x, y) we solve the wave equation Following the procedure described in [26], we integrate Eq. (6) over −Lz/2 dz taking into account the boundary conditions on the metal planes of the cavity, E x,y (±L z /2) = 0. Then we substitute Eq. (5) into the result of integration, which gives The solution to Eq. (7) with zero boundary conditions at the edges of the cavity determines eigenfrequencies and the structure of the eigenmodes for a quasi-2D cavity with an arbitrary shape in the (x, y)-plane.
Similar equations can be derived if one simply utilizes jumps of the dielectric constants on the sides instead of metal coating. Even without any jump in the dielectric constants, an open end of a thin waveguide with vertical size much smaller than wavelength is a good reflector and therefore any radiation losses through the facets are small and are not affecting the mode spatial structure significantly.
The expression for the constants D p,s,i for quantized fields can be obtained from the general equation Eq. (4) (see also [26]), where ν = p, s, i. For a simple case of a rectangular-shaped cavity, when S = L x × L y , where L x and L y are the cavity dimensions along x and y directions, it is easy to obtain useful analytic expressions for the modal spatial structure and frequencies. For eigenmodes with one half-wavelength along the y-axis and N half-wavelengths along the x-axis, we obtain the modal profile The eigenfrequencies are given by For a particular case of a uniform dielectric constant, Eqs. (9) and (11) are exact, i.e. they do not require an assumption of a quasielectrostatic field structure.
We assume that a 2D electron system with the second-order nonlinear susceptibility is positioned in the cavity. The material can be a quantum well (QW), a 2D semiconductor, and even graphene, which has a strong second-order nonlinearity beyond electric-dipole approximation, despite being centrosymmetric [14]. The second-order nonlinearity gives rise to the nonlinear polarization at signal and idler frequencies. The excitation equations for the cavity modes derived from the operator-valued Maxwell's equations [27] take the forṁ The nonlinear polarizationP N L (r, t) should be determined for a given electron system; in general, it has a nonuniform distribution over the cavity cross-section. However, it is obvious from Eq. (12) that only the integral over the nonlinear polarization matters. Therefore it is convenient to consider a nonlinear layer with uniform dielectric constant ε QW (ω p,s,i ) = ε p,s,i and uniform second-order nonlinear susceptibilities p . For a general case of a nonuniform layer the above quantities can be considered as parameters obtained as a result of integration in Eq. (12).
For a uniform layer the nonlinear polarization can be expressed aŝ and E p = const. In Eq. (14) we assumed a rectangular cavity shape for simplicity.
If the nonlinearity is non-dissipative, the nonlinear susceptibilities satisfy the symmetry properties which ensure than Manley-Rowe relationships are satisfied in a conservative system [30,31].
Using the rotating wave approximation, Eqs. (12) and (13) give the following parametrically coupled equations for a given classical pumping: is a modal overlap factor and l is the thickness of the nonlinear layer in z direction.
Equations (16) are Heisenberg equations for the effective Hamiltonian For a parametric down-conversion of a pump photon into two identical photons, when we arrive at the Hamiltonian The condition J = 0 is similar to the three-wave phase matching condition for the wave vectors of modes participating in the parametric down-conversion. The phase matching needs to be satisfied together with frequency matching in Eq. (1), which could be highly nontrivial in a 3D geometry and in a dispersive medium. An important advantage of a subwavelength cavity is that these conditions are straightforward to satisfy by adjusting the cavity geometry. Indeed, consider the decay of the pump into two lowest-order TE 011 modes satisfying Eq. (20); in this case ζ s (x, y) = cos πy L y cos πx L x . For a pumping mode of TE 01N type with N even, J = 0; however for N odd we get J = 0. For example, for a TE 013 pumping mode (see Fig. 1), when ζ p (x, y) = cos πy L y cos 3πx L x , we obtain: . In this case from Eq. (20) and mode frequencies given by Eq. (11) one can get a condition for cavity sizes:

III. EQUATIONS FOR PARAMETRIC DOWN-CONVERSION IN A DISSIPA-TIVE SYSTEM. HEISENBERG-LANGEVIN APPROACH
Here we take into account absorption and radiative losses within the Heisenberg-Langevin formalism. We remind the reader that this approach assumes that the coupling of a dynamic subsystem to a dissipative reservoir is sufficiently weak. If this is not the case and the coupling is strong, the process of energy loss by a given EM mode should be described within a closed Hamiltonian system (e.g. as a coupling to other EM modes, phonons etc.).
In this case one does not need any Langevin sources, because in a Hamiltonian system proper commutation relations are satisfied automatically. Whenever the energy exchange of a dynamical subsystem with a reservoir is relatively weak and can be considered within a phenomenological approach, the "dissipation + the Langevin noise" model should be valid for any mechanism of dissipation.
For example, we assume here that the spatial mode structure corresponds to the one in an ideal cavity, whereas diffraction losses of the field out of a cavity can be described by an effective loss rate. This assumption obviously works as long as losses do not affect the mode structure significantly. If the latter is not true, one would have to solve a rigorous diffraction problem which couples the field modes in the cavity and outside. For such a rigorous problem all commutation relations would be satisfied automatically.
Introducing operators of slowly varying field amplitudes, namelyĉ s,i =ĉ 0s,i (t)e −iω s,i t , c † s,i =ĉ † 0s,i (t)e +iω s,i t , we obtain from Eqs. (16) the following equations: where Γ s,i = Γ r(s,i) + Γ σ(s,i) , the coefficients Γ r(s,i) and Γ σ(s,i) denote, respectively, radiative losses due to the outcoupling of radiation from the cavity and absorption losses due to intracavity absorption.L s,i are the Langevin noise operators. We show in Appendix A that to preserve commutation relations [ĉ 0i ,ĉ † 0i ] = [ĉ 0s ,ĉ † 0s ] = 1 at Γ s,i = 0 the noise operators in the right-hand side of Eqs. (22) should satisfy the same commutation relations as in the case of one quantum oscillator [19,20,26] and they should also commute with each other: The fact that commutation relations are the same for one quantum oscillator and for two (or more) interacting oscillators is expected, since the processes within the Hamiltonian dynamics do not affect the commutators; this can be easily checked, for example for the system described by the Hamiltonian Eq. (19). Noise correlators can be defined by generalizing the simplest expression in [19] to the case of two absorbers with different temperatures: where n Tr,σ is the average number of thermal photons at temperature T r,σ ; T r and T σ denote the temperature outside and inside the cavity, respectively. Expressions (24) imply that dissipative and radiative noises are not correlated.
The solution to Eqs. (22) can be represented as [21,23]  ĉ 0ŝ c † where λ 1,2 and   1 K 1,2   are eigenvalues and eigenvectors of the 2 × 2 matrix: c s (0) andĉ † i (0) are initial conditions. From the solution (25)-(26) one can derive a standard-looking condition for the parametric instability(see e.g. [32]): Consider the inequality (28) in more detail, neglecting for clarity the frequency dispersion of the dielectric filling of a cavity. Taking into account Eq. (14), one can derive from Eq. (17) that where the dimensionless factor J/(L x L y ) depends only on the spatial structure of the modes with frequencies ω p,s,i . From Eqs. (28), (29) the instability condition takes the form where ∆ω s,i = 2Γ s,i are the linewidths for the signal and idler modes.
To avoid cumbersome expressions, consider the decay of a pump photon into identical quanta as in Eq. (20). In this case the instability condition is |ς| > Γ s . It is convenient to choose the phase of the pump mode so that the value of is real and positive. Then
In the limit of zero pumping intensity, Eq. (32) gives an expression which describes how the initial perturbation of a photon number approaches equilibrium: Above the instability threshold, when ς Γ s , it is enough to keep only exponentially growing terms in Eq. (32). We can also assume that an average number of thermal photons in an ambient space n Tr (ω s ) is negligible. This gives an expression for the parametric signal power P s ≈ 2ςhω s n s : Obviously this expression is valid only at the initial linear stage. The subsequent evolution is governed by the nonlinear pump depletion and nonlinear modification of phasematching conditions and losses. An order-of magnitude estimate of the maximum steady-state power can be obtained from Manley-Rowe relations as shown below for a specific example.
The fractions of the power emitted outside and absorbed inside a cavity are P rs ≈ Γ rs P s /ς and P σs ≈ Γ σs P s /ς respectively; most of the power is accumulated in a cavity. From Eq. (34) one can see that the amplification of intrinsic thermal noise of a QW layer with temperature When the parametric growth rate is lower than losses, ς < Γ s , the general solution Eq. (32) describes the regime of spontaneous parametric down-conversion. In the stationary limit, when (Γ s − ς)t → ∞, the radiated signal power P s ≈ 2Γ rsh ω s n s becomes The first term in the right-hand side of Eq. (35)  For a numerical example, consider a nanocavity filled with multiple quantum wells, excited with a mid-infrared pump, as reported in [8]. Using Eq. (30) for their values of intersubband nonlinearity |χ (2) | ∼ 3 × 10 −7 m/V, ω s,i /∆ω s,i ∼ 20 and J/L x L y ∼ 0.3 we obtain the intracavity pump field at the instability threshold to be E p 8 kV/cm, which is achievable and is lower than the saturation field for the intersubband nonlinearity. Above the threshold, the signal and idler fields start growing inside the cavity until they get limited by the Manley-Rowe relations [33], i.e. the intracavity signal field reaches |E s | |E p |/ √ 2.
Therefore, the maximum output signal power that can be obtained per one nanocavity described in [8] is about 8 × 10 −7 W for the photon leakage rate Γ rs = 10 12 s −1 .
Far below the instability threshold, when ς Γ s , the spontaneous rate of parametric down-conversion scales as (Γ rs /Γ 2 s )ς 2 . For the parameters from the above numerical example, when ς = Γ s /2 the emission rate of signal photons is around 3 × 10 11 s −1 and the power is 6 nW.
A very high second-order nonlinear surface conductivity for graphene was reported in [14], corresponding to the effective bulk susceptibility |χ (2) | ∼ 10 −3 m/V per monolayer in the THz range. This large susceptibility if partially offset by a small factor l/L z in Eq. (29) where l is a thickness of the graphene layer. However, for hBN-encapsulated graphene utilized to fabricate low-disorder graphene samples [34] the total cavity thickness L z can be as small as several nm, so the factor l/L z can be as large as 0.1. Even factoring in enhanced cavity losses, this can yield a lower parametric instability threshold as compared to semiconductor quantum well samples.

IV. COMPARING PARAMETRIC INSTABILITY IN A SUBWAVELENGTH CAVITY AND IN A FABRY-PEROT CAVITY
Compare the parametric instability in a subwavelength cavity with similar instability of modes in a Fabry-Perot (FP) cavity with all three dimensions larger than wavelength, which we will call a quasi-optical cavity. Consider a planar quasi-2D cavity of the surface area L x L y , in which the waves are propagating along the nonlinear layer of thickness l much smaller than the cavity thickness L F P transverse to the nonlinear layer, so the cavity volume is L F P L x L y .
The dielectric constant of a cavity filling is ε. In this case the parametric down-conversion is still described by Eqs. (22), (17), and (18), in which the relaxation constants Γ s,i and the overlap integral are determined by the FP cavity Q-factor and the corresponding spatial structure of the modes. For the normalization constants of the quantum fields entering Eq. (17) we use standard expressions for a two-mirror FP cavity: The resulting parametric instability threshold is where ∆ω s,i ≈ 2Γ s,i .
As we already pointed out, in a cavity with all three dimensions larger than the wavelength the phase matching condition for a three-wave mixing may be more difficult to satisfy.
Even if we assume that phase matching is somehow arranged and the geometric factor J/L x L y is of the same order as in a subwavelength cavity, the latter is expected to have a lower parametric threshold. Indeed the ratio of the threshold pump intensity |E p | 2 in a subwavelength cavity to that in a quasi-optical cavity scales as ∼ L z L F P 2 ∆ω ef f ∆ω F P

2
, where ∆ω ef f and ∆ω F P ≈ ∆ω s,i are typical linewidths of the subwavelength cavity and FP cavity modes, respectively. This ratio can be much smaller than 1, which indicates that a much lower pumping is needed to reach the parametric instability threshold in a subwavelength cavity, even if the FP cavity has a higher Q-factor as compared to the subwavelength cavity, A plane-parallel quasi-2D subwavelength cavity geometry considered in this paper is the most natural choice for integration with 2D nonlinear materials. However, other geometries are also possible, for example plasmonic or grating structures supporting surface modes. To get an order of magnitude estimate of the parametric threshold, one can use our results in Eqs. (30), (37) after replacing L z or L F P with a mode size transversely to the nonlinear layer.
A promising example of such a plasmonic nanocavity was reported in [35]. It consists of a monolayer MoS 2 , which is a 2D semiconductor, sandwiched between a gold substrate and a patch silver nanoantenna. Such a cavity has high radiative and absorption losses but a very small transverse mode size of less than 10 nm and an ultrasmall effective mode volume of ∼ 10 −3 (λ/ √ ) 3 . The authors of [35] used their cavities to obtain a 2000-fold enhancement in the photoluminescence intensity from MoS 2 monolayer. However, a cavity of similar design can also be used for parametric down-conversion from visible to the near-IR range. A high second-order nonlinearity, about an order of magnitude higher than in BBO or lithium niobate, has been reported for MoS 2 monolayer [36]. An even higher nonlinearity has been observed in the vicinity of exciton resonances [37]. Assuming conservatively that the effective second-order susceptibility for MoS 2 is |χ (2) | ∼ 10 −10 m/V, monolayer thickness 0.6 nm, transverse mode size 5 nm, ω s,i /∆ω s,i ∼ 20 and J/L x L y ∼ 0.3 we obtain from Eq. (30) the intracavity pump field at the parametric amplification threshold to be E p 30 MV/cm, which is much higher than the estimate above for a nonlinear cavity based on mid-infrared resonant intersubband nonlinearity of quantum wells, but is below the saturation threshold for MoS 2 and easily achievable with pulsed lasers.
Ultracompact subwavelength electrodynamic structures utilizing 2D materials are promising for applications in integrated photonic circuits, whenever one needs a compact planar architecture. At the same time, due to strong dissipation they are unlikely to outperform conventional nonlinear devices made of bulk transparent nonlinear materials when it comes to the nonlinear conversion efficiency and power. For example, in [38] the authors realized low-threshold mode-matched parametric generation in whispering gallery mode resonators made entirely of bulk lithium niobate. In this case the bulk nonlinear material occupies all modal volume. The lower nonlinearity and the loss of Purcell enhancement are compensated by lower dissipation and increased interaction volume.
In conclusion, we applied a consistent Heisenberg-Langevin formalism to the process of nonlinear parametric down-conversion of cavity modes in planar subwavelength cavities containing 2D nonlinear materials. We derived general analytic formulas for the spontaneous parametric signal and threshold of stimulated parametric down-conversion of a pump cavity mode into the signal and idler modes. We found that a significant reduction in the parametric instability threshold can be achieved for realistic materials and cavity parameters due to Consider first a single quantum oscillator described by the HamiltonianĤ =hω(ĉ †ĉ + 1/2). After substitutingĉ =ĉ 0 e −iωt andĉ † =ĉ † 0 e −iωt the Heisenberg equations of motion take the formċ 0 = 0,ċ † 0 = 0. The simplest model of interaction with a dissipative reservoir modifies these equations as follows:ċ 0 + Γĉ 0 = 0,ċ † 0 + Γĉ † 0 = 0. However, this modification leads to violation of boson commutation relation [ĉ 0 ,ĉ † 0 ] = 1. To resolve this issue and preserve the commutator one has to add the Langevin sources to the right-hand side of Heisenberg equations [19,20,26]: Langevin noise operators in Eq. (A1) describe fluctuations in a dissipative system. Note that L = 0 ; the notation · · · means averaging over the statistics of the dissipative reservoir and over the initial quantum state |Ψ within the Heisenberg picture.
The commutation relations for a noise operator can be obtained directly from the given form of the relaxation operator if we require that standard commutation relations [ĉ 0 ,ĉ † 0 ] = 1, [ĉ 0 ,ĉ 0 ] = 0, be satisfied at any moment of time [20,26]. Indeed, let's substitute the solution of the operator-valued equations (A1) into the commutators. It is easy to see that the standard commutation relations will be satisfied if, first of all, the field operators at an initial moment of time,ĉ 0 (0) andĉ † 0 (0), commute with Langevin operatorsL(t) andL † (t) in any combination. Second, the following condition has to be satisfied: K 2 e λ 1 (t −t) − K 1 e λ 2 (t −t) L s (t ),L † s (t) dt t 0 K 2 e λ 2 (t −t) − K 1 e λ 1 (t −t) L i (t ),L † i (t) dt From Eqs. (A6) and (A7) one can obtain the requirement Eq. (23) which preserves correct commutators of the field operators. Therefore, the commutation properties of correct noise operators for coupled oscillators have to be exactly the same as for uncoupled isolated oscillators.
Here we presented a general proof which does not rely on any specific microscopic model of a dissipative subsystem coupled to the field oscillators. The proof for a particular case of two identical coupled oscillators interacting with a standard dissipative reservoir of equilibrium harmonic oscillators [19] has been recently obtained in [39].