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

Enhancement of the radiative processes due to the localization of emitters in a subwavelength cavity (the so-called Purcell enhancement¹) 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 light-emitting diodes 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., Refs. 2–10 and references therein. The rise of two-dimensional (2D) materials with atomic monolayer thickness and excellent nonlinear optical properties, such as graphene^11–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 and 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 (QW)] 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 down-conversion. 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., Refs. 19–23. Its generalization to systems far from equilibrium, with arbitrary dissipation mechanisms and arbitrary photon density of states is nontrivial; see, e.g., Refs. 24 and 25. Recent work¹⁰ 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 this 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 quantized field (e.g., Refs. 19, 20, and 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., Eqs. (32) and (35) below.

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 a 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 (FP) cavity containing a 2D nonlinear layer. In this case, the performance trade-off 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 in the parametric amplification threshold due to Purcell enhancement in quasi-2D subwavelength cavities.

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 the operator,

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–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 Γ: |(∂$ε$_jk(ω)/∂ω)Γ| ≪ |$ε$_jk(ω)|. 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: |(∂²$ε$_jk(ω)/∂ω²)Γ| ≪ |(∂$ε$_jk(ω)/∂ω)|.

Consider a 3D cavity filled with an isotropic dielectric medium, as sketched in Fig. 1. The cavity thickness in the z-direction is much smaller than wavelength, $Lz≪c/ε¯ω$⁠, where $ε¯$ is a typical (average) value of the dielectric constant of the filling. As was shown in Ref. 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,

Consider a z-component of this equation,

Following the procedure described in Ref. 26, we integrate Eq. (6) over $∫−Lz/2Lz/2dz$ 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 Ref. 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.¹⁴ 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²⁷ take the form

The nonlinear polarization $P^NL(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 $χ(2)(ωs=ωp−ωi)=χs(2)$⁠, $χ(2)(ωi=ωp−ωs)=χi(2)$⁠, and $χ(2)(ωp=ωs+ωi)=χp(2)$⁠. 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 as

where $cs†∝eiωst$⁠, $ci†∝eiωit$⁠,

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 that 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:

where

is a modal overlap factor, and l is the thickness of the nonlinear layer in the 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₀₁₁ modes satisfying Eq. (20); in this case, $ζs(x,y)=cosπyLycosπxLx$⁠. 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₀₁₃ pumping mode (see Fig. 1), when $ζp(x,y)=cosπyLycos3πxLx$⁠, we obtain $J=LxLyπ2×1645$⁠. In this case from Eq. (20) and mode frequencies given by Eq. (11), one can get a condition for cavity sizes

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 and phonons). 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,it$ and $ĉs,i†=ĉ0s,i†(t)e+iωs,it$⁠, we obtain from Eq. (16) the following equations:

where Γ_s,i = Γ_r(s,i) + Γ_σ(s,i) and 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 the  Appendix that to preserve commutation relations $[ĉ0i,ĉ0i†]=[ĉ0s,ĉ0s†]=1$ at Γ_s,i ≠ 0 the noise operators in the right-hand side of Eq. (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 Ref. 19 to the case of two absorbers with different temperatures,

where $nTr,σ$ 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 Eq. (22) can be represented as^21,23

where λ_1,2 and $1K1$ are eigenvalues and eigenvectors of the 2 × 2 matrix,

where ĉ_s(0) and $ĉi†(0)$ are initial conditions.

From the solution (25) and (26), one can derive a standard-looking condition for the parametric instability (see, e.g., Ref. 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_xL_y) depends only on the spatial structure of the modes with frequencies ω_p,s,i. From Eqs. (28) and (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 is real and positive. Then Eqs. (25)–(27) yield

Taking into account the properties of Langevin operators in Eq. (22) and taking $ĉs†(0)ĉs†(0)=ĉs(0)ĉs(0)=0$ as an initial state, one can derive from Eq. (31) the average photon numbers for signal modes $ns=ĉs†ĉs=ĉ0s†ĉ0s$⁠,

where Γ_s = Γ_σs + Γ_rs. When the parametric amplification starts from the level of vacuum fluctuations, one should put n_s(0) = 0 in Eq. (32).

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 $nTr(ωs)$ is negligible. This gives an expression for the parametric signal power P_s ≈ 2ςℏω_sn_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 ≈ Γ_rsP_s/ς and P_σs ≈ Γ_σsP_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 T_σ can be ignored if $Γσsς⋅2exp(ℏωs/Tσ)−1≪1$⁠.

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Γ_rsℏω_sn_s becomes

The first term in the right-hand side of Eq. (35) is due to the thermal emission modified by the parametric decay of the pump photons. The second term originates from the parametric decay of the pump photons under the action of vacuum fluctuations of the intracavity field; this is a purely quantum effect. Thermal effects can be neglected if

The results in Eqs. (32)–(35) provide the dependence of the parametric signal from all relevant dissipation and fluctuation effects. In addition to dissipation and thermal fluctuations due to absorption in the cavity walls and a semiconductor heterostructure, they take into account outcoupling of a signal into the ambient space and coupling to thermal photons from the environment. Equation (35) determines the spontaneous parametric signal emitted from a cavity against the background of noise created by both thermal radiation from a cavity and reemission of thermal photons coupled into a cavity from the outside. The background noise depends from the cavity temperature and the environment temperature. In addition to the spontaneous decay process, we take into account the modification of background noise by pumping.

For a numerical example, consider a nanocavity filled with multiple quantum wells, excited with a mid-infrared pump, as reported in Ref. 8. Using Eq. (30) for their values of intersubband nonlinearity |χ⁽²⁾| ∼ 3 × 10⁻⁷ m/V, ω_s,i/Δω_s,i ∼ 20, and J/L_xL_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,³³ i.e., the intracavity signal field reaches $|Es|≃|Ep|/2$⁠. Therefore, the maximum output signal power that can be obtained per one nanocavity described in Ref. 8 is about 8 × 10⁻⁷ W for the photon leakage rate Γ_rs = 10¹² s⁻¹.

Far below the instability threshold, when ς ≪ Γ_s, the spontaneous rate of parametric down-conversion scales as $(Γrs/Γs2)ς2$⁠. For the parameters from the above numerical example, when ς = Γ_s/2, the emission rate of signal photons is around 3 × 10¹¹ s⁻¹ and the power is 6 nW.

A very high second-order nonlinear surface conductivity for graphene was reported in Ref. 14, corresponding to the effective bulk susceptibility |χ⁽²⁾| ∼ 10⁻³ m/V per monolayer in the THz range. This large susceptibility is partially offset by a small factor l/L_z in Eq. (29) where l is the thickness of the graphene layer. However, for hBN-encapsulated graphene utilized to fabricate low-disorder graphene samples,³⁴ 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.

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_xL_y, in which the waves are propagating along the nonlinear layer of thickness l much smaller than the cavity thickness L_FP transverse to the nonlinear layer, so the cavity volume is L_FPL_xL_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_xL_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|² in a subwavelength cavity to that in a quasi-optical cavity scales as $∼LzLFP2ΔωeffΔωFP2$⁠, where Δω_eff and Δω_FP ≈ Δω_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, Δω_FP < Δω_eff.

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) and (37) after replacing L_z or L_FP with a mode size transversely to the nonlinear layer.

A promising example of such a plasmonic nanocavity was reported in Ref. 35. It consists of a monolayer MoS₂, 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 Ref. 35 used their cavities to obtain a 2000-fold enhancement in the photoluminescence intensity from MoS₂ 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 β-barium borate or lithium niobate, has been reported for MoS₂ monolayer.³⁶ An even higher nonlinearity has been observed in the vicinity of exciton resonances.³⁷ Assuming conservatively that the effective second-order susceptibility for MoS₂ is |χ⁽²⁾| ∼ 10⁻¹⁰ m/V, monolayer thickness 0.6 nm, transverse mode size 5 nm, ω_s,i/Δω_s,i ∼ 20, and J/L_xL_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₂ 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 Ref. 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 Purcell enhancement.

This material is based on the work supported by the Air Force Office of Scientific Research under Award Nos. FA9550-17-1-0341 and FA9550-15-1-0153. M.T. acknowledges the support from RFBR Grant No. 17-02-00387 and from the program of the Presidium of the Russian Academy of Sciences “Nonlinear dynamics in mathematical and physical sciences” (Project No. 0035-2018-0006).

Consider first a single quantum oscillator described by the Hamiltonian Ĥ = ℏω(ĉ^†ĉ + 1/2). After substituting ĉ = ĉ₀e^−iωt and $ĉ†=ĉ0†e−iωt$⁠, the Heisenberg equations of motion take the form $ĉ̇0=0$ and $ĉ̇0†=0$⁠. The simplest model of interaction with a dissipative reservoir modifies these equations as follows: $ĉ̇0+Γĉ0=0$ and $ĉ̇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 us 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) and $ĉ0†(0)$⁠, commute with Langevin operators $L^(t)$ and $L^†(t)$ in any combination. Second, the following condition has to be satisfied:

Substituting Eq. (A2) into Eq. (A3) and using the identity $∫0tX(t′)δ(t−t′)dt′=X(t)/2$⁠, we arrive at

Now consider an ensemble of coupled oscillators [Eq. (22)]. One can find directly from the solution [Eq. (25)] that the following conditions have to be satisfied in order to preserve standard commutation relations $[ĉ0s,ĉ0s†]=[ĉ0i,ĉ0i†]=1,[ĉ0s,ĉ0i]=0$⁠, etc.:

It is easy to find out that Eq. (A5) will be satisfied if the field operators at t = 0 commute with Langevin noise operators in any combination, and the noise operators $L^s$ and $L^i$ commute with each other. In addition, substituting Eqs. (25)–(27) into Eq. (A5), one can show that in order to satisfy Eq. (A5) the following relations must hold:

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¹⁹ has been recently obtained in Ref. 39.