Scully-Lamb quantum laser model for parity-time-symmetric whispering-gallery microcavities: Gain saturation effects and non-reciprocity

TL;DR

This paper models light propagation in parity-time symmetric whispering-gallery microcavities using a Scully-Lamb laser approach, revealing how gain saturation influences non-reciprocity and exceptional points in non-Hermitian systems.

Contribution

It introduces a non-Lindbladian Scully-Lamb laser model for analyzing gain saturation effects and non-reciprocity in parity-time symmetric microcavities, extending understanding of exceptional points.

Findings

Gain saturation causes non-reciprocal light propagation.

Gain saturation modifies the position of exceptional points.

The model explains experimentally observed non-reciprocity.

Abstract

We use a non-Lindbladian master equation of the Scully-Lamb laser model for the analysis of light propagation in a parity-time symmetric photonic system composed of coupled active and passive whispering-gallery microresonators. Performing the semiclassical approximation, we obtain a set of two nonlinear coupled differential equations describing the time evolution of intracavity fields. These coupled equations are able to explain the experimentally-observed light non-reciprocity [Peng {\em et al.}, Nature Physics {\bf 10}, 394 (2014), Chang {\em et al.}, Nature Photonics {\bf 8}, 524 (2014)]. We show that this effect arises from the interplay between gain saturation in the active microcavity, intercavity coupling, and losses in the cavities. Additionally, using this approach, we study the effect of the gain saturation on exceptional points, i.e., exotic degeneracies in non-Hermitian systems. Namely, we demonstrate that the inclusion of gain saturation leads to a modification of the exceptional points in the presence of intense intracavity fields. The Scully-Lamb master equation for systems of coupled optical structures, as proposed and applied here, constitutes a promising tool for the study of quantum optical effects in coupled systems with losses, gain, and gain saturation.

Equations

$$\hat{H}=\sum _{k=1}^2\hslash {\omega }_k{\hat{a}}_k^{†}{\hat{a}}_k+i\hslash \left[\kappa {\hat{a}}_1{\hat{a}}_2^{†}+ϵ{\hat{a}}_1{e}^{i{\omega }_{l}t}-h.c.\right],$$

$$\frac{d}{dt}\hat{\rho }$$

$$A=\frac{2g^{2}r}{Y^2},\text{and}B=\frac{4g^2}{Y^2}A,$$

$${\Gamma }_i=C_i+{\gamma }_i,\mathrm{where}{C}_i=\frac{{\omega }_i}{Q_i}$$

$$\frac{d}{dt}\hat{\rho }=\frac{1}{i\hslash }\left[\hat{H},\hat{\rho }\right]-\frac{1}{2}\sum _{i=1}^4\left({\hat{L}}_i^{†}{\hat{L}}_i\hat{\rho }+\hat{\rho }{\hat{L}}_i^{†}{\hat{L}}_i-2{\hat{L}}_i\hat{\rho }{\hat{L}}_i^{†}\right),$$

$${\hat{L}}_1=\sqrt{A}{\hat{a}}_1^{†}\left(1-\frac{B}{2A}{\hat{a}}_1{\hat{a}}_1^{†}\right),{\hat{L}}_2=\frac{1}{2}\sqrt{3B}{\hat{a}}_1{\hat{a}}_1^{†},$$

$${\hat{L}}_3=\sqrt{{\Gamma }_1}{\hat{a}}_1,{\hat{L}}_4=\sqrt{{\Gamma }_2}{\hat{a}}_2.$$

$$\frac{d}{dt}⟨{\hat{a}}_j⟩=\mathrm{Tr}\left[{\hat{a}}_j\frac{d}{dt}\hat{\rho }\right],j=1,2,$$

$$\frac{d}{dt}⟨{\hat{a}}_1⟩$$

$$\frac{d}{dt}⟨{\hat{a}}_2⟩$$

$$\frac{d}{dt}{a}_1$$

$$\frac{d}{dt}{a}_2$$

$$\frac{d}{dt}{A}_1$$

$$\frac{d}{dt}{A}_2$$

$$\left(i\Delta +\frac{G_1^{\prime }}{2}\right)A_1-\kappa A_2$$

$$\left(i\Delta -\frac{{\Gamma }_2}{2}\right)A_2+\kappa A_1$$

$$i\frac{\mathrm{d}\alpha }{\mathrm{d}t}=M\alpha ,$$

$$M=\left(\begin{array}{cc}-\Delta +i\frac{G_1^{\prime }}{2}& -i\kappa \\ i\kappa & -\Delta -i\frac{{\Gamma }_2}{2}\end{array}\right).$$

$$\left(\Delta -i\frac{G_1^{\prime }}{2}\right)\left(\Delta +i\frac{{\Gamma }_2}{2}\right)-{\kappa }^2=0,$$

$${\Delta }_{\pm }=i\frac{G_1^{\prime }-{\Gamma }_2}{4}\pm \sqrt{4{\kappa }^2-\frac{1}{4}{\left(G_1^{\prime }+{\Gamma }_2\right)}^2}.$$

$${\omega }_{\pm }$$

$$A-C_1-C_2-B{I}_1=0.$$

$${\omega }_{\pm }^{\mathrm{PT}}={\omega }_c-i\frac{\gamma }{2}\pm \frac{1}{2}\sqrt{4{\kappa }^2-C_2^2}.$$

$${\kappa }_{\mathrm{EP}}=\frac{C_2}{2}.$$

$$A-C_1-C_2=0.$$

$${\kappa }_{\mathrm{EP}}=\frac{1}{4}\left|A-C_1+C_2-B{I}_1\right|=\frac{1}{4}\left|2C_2-B{I}_1\right|.$$

$$B{I}_1\approx C_2\ll A\approx C_1.$$

$${\lambda }_1{I}_1^3+{\lambda }_2{I}_1^2+{\lambda }_3{I}_1+{\lambda }_4=0,$$

$${\lambda }_1=\frac{B^2}{4},{\lambda }_2=BF,{\lambda }_4=-{ϵ}^2,$$

$$T(\omega )={\left|\frac{A_{\mathrm{out}}}{A_{\mathrm{in}}}\right|}^2.$$

$$T_{1\to 4}(\Delta )=\frac{4{\kappa }^2{\gamma }_1{\gamma }_2}{{ϵ}^2({\Gamma }_2^2+4{\Delta }^2)}{I}_1.$$

$$T_{4\to 1}(\omega )=\frac{{\gamma }_1{\gamma }_2}{{ϵ}^2}{I}_1.$$

$${\hat{H}}_I=ig\left({\hat{\sigma }}_{+}{\hat{a}}_1-{\hat{\sigma }}_{-}{\hat{a}}_1^{†}\right)+i\kappa \left({\hat{a}}_1{\hat{a}}_2^{†}-{\hat{a}}_1^{†}{\hat{a}}_2\right).$$

Full text

I Introduction

Recent years has witnessed an increasingly intense research effort to explore a class of non-Hermitian systems described by parity-time (${\cal PT}$) symmetric Hamiltonians Bender and Boettcher (1998), (for reviews see Bender (2007); Bender et al. (2003)). A system described by the Hamiltonian $H$ is ${\cal PT}$-symmetric if it is invariant under the combined action of the parity ${\cal P}$ and the time-reversal ${\cal T}$ operators (i.e., $H$ commutes with the ${\cal PT}$ operator: $[H,{\cal PT}]=0$) but not necessarily with the $\cal P$ or $\cal T$ operator alone. An important consequence of this is the necessary, but not sufficient, condition for ${\cal PT}$-symmetry: The complex potential $V(x)=V_{\rm r}(x)+iV_{\rm i}(x)$ of the Hamiltonian should satisfy $V(x)=V^{\ast}(-x)$, where the superscript $\ast$ denotes complex conjugation. In other words, the real part of the potential should be an even function of $x$, while its imaginary part should be an odd function of $x$, i.e., $V_{\rm r}(x)=V_{\rm r}(-x)$ and $V_{\rm i}(x)=-V_{\rm i}(-x)$. A ${\cal PT}$-symmetric system exhibits two very distinct phases: an unbroken ${\cal PT}$ phase (also known as the exact-${\cal PT}$ regime), where the Hamiltonian supports real eigenvalues despite being non-Hermitian, and a broken ${\cal PT}$ phase where some eigenvalues form complex conjugate pairs. The transition between these two phases takes place spontaneously as a result of parametric variation of the Hamiltonian. This real-to-complex spectral phase transition or the ${\cal PT}$-phase transition point exhibits all properties of an exceptional point (EP), which is defined as a singularity in the parameter space of a non-Hermitian system at which two or more eigenvalues and their associated eigenvectors coalesce.

A decade after the ground-breaking work of Bender and Boettcher in 1998 which initiated the mathematical framework and fundamental understanding of ${\cal PT}$-symmetric systems Bender and Boettcher (1998), it was realized that ${\cal PT}$-symmetry and its breaking at an EP can be observed in photonics by imposing the necessary condition for ${\cal PT}$-symmetry on a complex optical potential, that is on the complex refractive index, $n(x)=n_{\rm r}(x)+in_{\rm i}(x)$, which leads to $n_{\rm r}(x)=n_{\rm r}(-x)$ and $n_{\rm i}(x)=-n_{\rm i}(-x)$ El-Ganainy et al. (2007); Makris et al. (2008). Thus, an optical system with the ${\cal PT}$-symmetric potential has a symmetric index profile but an asymmetric gain/loss profile. Such a refractive index profile can be obtained in two coupled optical structures, such as waveguides or resonators: one having loss and the other having gain compensating the loss of the other. This discovery opened a very fertile research direction, where the interplay between gain, loss, and the strength of the coupling between them provides entirely new features and device functionalities Feng et al. (2017); El-Ganainy et al. (2018); Miri and Alu (2019); Ş. K. Özdemir et al. (2019). In these non-Hermitian systems, with coupled loss and gain components, EPs can be observed by controlling (or tuning) the coupling strength to balance amplification (gain) and dissipation (loss). Thus, EPs can drastically alter the overall response of the system. This leads to a plethora of nontrivial phenomena Feng et al. (2017); El-Ganainy et al. (2018); Miri and Alu (2019), such as enhanced light-matter interactions Liu et al. (2016); Chen et al. (2017); Hodaei et al. (2017), unidirectional invisibility Lin et al. (2011); Regensburger et al. (2012), lasers with enhanced mode selectivity Feng et al. (2014); Hodaei et al. (2014), low-power nonreciprocal light transmission Peng et al. (2014a); Chang et al. (2014), loss-induced lasing Brandstetter et al. (2014); Peng et al. (2014b), thresholdless phonon lasers Jing et al. (2014); Lü et al. (2017) to name a few. In parallel to these efforts in photonics, the concepts of ${\cal PT}$-symmetry have been put into use in electronics Schindler et al. (2011), optomechanics Jing et al. (2014); Xu et al. (2016); Jing et al. (2017), acoustics Zhu et al. (2014); Fleury et al. (2015), plasmonics Benisty et al. (2011), and metamaterials Kang et al. (2013). More recently, there is a trend in investigating the features of ${\cal PT}$-symmetric quantum systems and the effect of ${\cal PT}$-symmetry and its breaking on the quantum states of light and the properties of quantum information Quijandria et al. (2018); Kawabata et al. (2017); Ashida et al. (2017); Liu et al. (2019).

Although ${\cal PT}$-symmetry and related concepts have their roots in quantum field theories, the experimental demonstrations and the majority of theoretical works are focused on classical systems, as in the example of two coupled optical microresonators: where the energy loss in one of them is compensated by the gain in the other, and the system is probed with light from a laser. In experiments, gain in such systems can be provided optically via parametric gain or Raman gain of the amplifying material; while the resonator is made from emitters or rare-earth ions embedded in it Peng et al. (2014a); Chang et al. (2014). The theoretical framework to analyze such systems relies on linearly-coupled rate equations of classical fields, where the gain and loss correspond to a different sign of the imaginary part of the complex frequencies (e.g., minus for gain and plus for loss, or vice versa) without a reference to how the gain and loss are generated. For example, the theoretical model, developed in Ref. Peng et al. (2014a), was linear without any nonlinearity, although the experimental results reported in that paper showed the presence of nonlinearity leading to nonreciprocal light transmission. Thus, the theoretical framework, applied there, failed to describe the observed non-reciprocity. On the other hand, in Ref. Chang et al. (2014), closely following Ref. Peng et al. (2014a), a term was added phenomenologically (i.e., by hand) to the linear rate equations to include the effect of gain saturation, which provides the required nonlinearity for the nonreciprocal light transmission, without explicitly describing where this term comes from. One is also able to incorporate a nonlinear term in the rate equations to explain light non-reciprocity, but in the ${\cal PT}$-symmetric system of coupled waveguides Ramezani et al. (2010), by resorting to a semi-classical Maxwell-Bloch approximation, where the dynamics is governed by a spatial variable, not time. On a different note, as the optical field is shifting from the classical to the quantum realm, it is important that the rate equations for the field operators are studied and the quantum-mechanical origins of gain and loss are properly described and incorporated into the models. Thus, a theoretical framework that addresses these concerns is highly desirable. Moreover, since the EP depends sensitively on the nonlinear coefficients, we can argue that having an ab-initio model is important for the question that this paper addresses. For this, we revisit the system, studied in Ref. Peng et al. (2014a), and later on in Ref.Chang et al. (2014), of two coupled optical structures with loss, gain, and gain saturation by using a non-Lindbladian master equation originally derived for the Scully-Lamb laser model III et al. (1974); Yamamoto and Imamoǧlu (1999); Orszag (2008). We apply this master equation for the density operators of the optical fields in optical structures. Our approach explains (at a fundamental level) the results obtained in Refs. Peng et al. (2014a); Chang et al. (2014) and explicitly describes the interplay of loss, gain, and gain saturation.

The paper is organized as follows. Section II contains the Hamiltonian of the system and the Scully-Lamb laser quantum master equation for the density operators of the optical fields. Section III presents the results, which include the rate equations of the field operators and their classical limit showing explicitly the presence of gain saturation nonlinearity and its effect on the eigenfrequencies and exceptional points of the system, as well as on the transmission spectra. Sections IV and V include a discussion of our results, future prospects, and a summary of the findings of this study.

II Laser quantum master equation for optical fields

in coupled active-passive microresonators

Let us start from a description of the physical system that we would like to investigate in the Scully-Lamb laser model III et al. (1974); Yamamoto and Imamoǧlu (1999); Orszag (2008). Our interest is focussed on the system which was experimentally studied first in Ref. Peng et al. (2014a) and later on in Chang et al. (2014). This system consists of two coupled whispering-gallery microresonators $R_{1}$ and $R_{2}$, where $R_{1}$ is an active microresonator used for the amplification of optical fields, and $R_{2}$ is a passive microcavity, which only damps the propagating fields. We denote the $Q$-factors of the first and second microresonator as $Q_{1}$ and $Q_{2}$, respectively. The schematic diagram of the described system is shown in Fig. 1. The coupling strength between the two microresonators is characterized by the real-valued parameter $\kappa$. Moreover, each microresonator is coupled to a different waveguide denoted as WG1 and WG2 with coupling constants $\gamma_{1}$ and $\gamma_{2}$, respectively. The system shown in Fig. 1 is $\cal PT$-symmetric for balanced gain and loss. This is possible because the microresonators $R_{1}$ and $R_{2}$ become interchanged under the parity reflection $\cal P$, while loss and gain are interchanged under the time-reversal operation $\cal T$.

The active microresonator $R_{1}$ can be considered as a laser system with a laser gain medium. Naturally, in order to describe the dynamics of the electromagnetic field in the active laser microresonator one would resort to the quantum laser theory and its master equation III et al. (1974); Yamamoto and Imamoǧlu (1999); Orszag (2008).

The Hamiltonian of the coupled microresonators with one driving coherent field, which drives the active cavity, can be written as

[TABLE]

where $\hat{a}_{k}$ ($\hat{a}_{k}^{\dagger}$) is the boson annihilation (creation) operator of the mode $k=1,2$, with frequency $\omega_{k}$; and h.c. denotes Hermitian conjugate. Moreover, $\kappa$ is the coupling strength between the microresonators, and $\epsilon$ is the coupling strength between the active cavity and the input signal driving field, with power ${P}$ and frequency $\omega_{l}$, which are related by $\epsilon\equiv\sqrt{\gamma_{1}{P}/(\hbar\omega_{l})}$. For the case when the passive cavity is driven, it is enough to swap the boson operator $\hat{a}_{1}$ with $\hat{a}_{2}$ in the last term of the Hamiltonian in Eq. (1).

By following the analysis of Yamamoto and Imamoǧlu Yamamoto and Imamoǧlu (1999), we consider our system as an ideal laser system in the Scully-Lamb laser model, which can be described by a non-Lindbladian master equation. This equation for the density operator $\hat{\rho}$ of the optical fields and for the Hamiltonian $\hat{H}$ reads as follows (see Appendix A, for details):

[TABLE]

where the gain $A$ and gain saturation $B$ coefficients for the field in the active cavity are expressed as

[TABLE]

respectively. The parameter $g$ stands for the coupling strength between the atoms of the gain medium and the optical field in the active cavity, $Y$ is a decay rate of the atoms, and $r$ accounts for the pump rate of the gain medium. In Eq. (2), the decaying rates for both cavities are denoted by ($i=1,2$)

[TABLE]

is the intrinsic loss of the $i$th cavity, and $\gamma_{i}$ stands for the loss due to the coupling of the $i$th cavity to the $i$th waveguide.

We note that the derivation of the master equation in Eq. (2), as carried out in Appendix A, is based on two main assumptions: (i) the adiabatic elimination of the population in the gain medium of the active resonator, and (ii) the weak-gain-saturation regime, i.e., when the laser in the active cavity operates not far from the lasing threshold. The non-Linbladian form of the master equation in Eq. (2), which is obtained within the fourth-order field approximation Yamamoto and Imamoǧlu (1999); Orszag (2008), is due to quantum jump-operator terms, which account for gain saturation.

We also note that the master equation, given in Eq. (2), can be recast to the Lindbladian form as Gea-Banacloche (1997):

[TABLE]

where the Lindblad operators $\hat{L}_{i}$ (for $i=1,\dots,4$) are defined as:

[TABLE]

The Lindblad form in Eq. (5) is equivalent to the master equation (2) if the terms of second order in $B\hat{a}_{1}\hat{a}_{1}^{\dagger}/(2A)$ are neglected in Eq. (5), which holds true for the weak gain saturation regime.

III Results

III.1 Rate equations

The master equation for the field operators, given in Eq. (2), and the Hamiltonian in Eq. (1), yield the rate equations for the averaged boson operators $\hat{a}_{1}$ and $\hat{a}_{2}$. Namely, by using the formula

[TABLE]

and utilizing the cyclic property of the trace operation, after substituting Eq. (2) in Eq. (7), we obtain

[TABLE]

where $G_{1}=A-\Gamma_{1}-\frac{7}{4}B$.

As one can see, the rate equations, given in Eq. (III.1) for the averaged quantum amplitudes, are nonlinear due to the presence of the gain saturation in the active cavity, i.e., the term $(B/2)\langle\hat{a}_{1}^{\dagger}\hat{a}_{1}\hat{a}_{1}\rangle$ represents the nonlinearity.

III.2 Classical Limit

In the classical limit, i.e., in the case of large intensities of the fields, the quantum field operators can be represented by $c$-number amplitudes as $\hat{a}_{i}\rightarrow\langle\hat{a}_{i}\rangle\equiv a_{i}$. Then the rate equations in Eq. (III.1) can be rewritten in the classical limit as

[TABLE]

where the term $\frac{B}{2}|a_{1}|^{2}a_{1}$ represents the nonlinearity due to gain saturation. By substituting $a_{k}=A_{k}(t)\exp({-i\omega t})$ into Eq. (III.2), one arrives at

[TABLE]

where $\Delta=\omega-\omega_{c}$ is the frequency-detuning parameter, and we assumed that the frequency of the driving coherent field is $\omega_{l}=\omega$, and the frequencies of the cavities are the same, $\omega_{1}=\omega_{2}=\omega_{c}$. We note that when the input driving field enters the passive cavity from port 4 in Fig. 1, then the term $(-\epsilon)$ would appear in the second equation (instead of the first) in Eq. (III.2). Importantly, in the semiclassical approximation, the coupling constant $\epsilon$ can be expressed via the complex amplitude $A_{\rm in}$ of the input driving field as $\epsilon\equiv\sqrt{\gamma_{j}}A_{\rm in}$, for $j=1,2$. Nevertheless, for convenience, we will keep the current notation for the input driving field via $\epsilon$, because it shares the same dimensionality with the other parameters describing the system. Moreover, by keeping this notation it will be more transparent and easier for us to compare the semiclassical and fully quantum approaches in related future research. In what follows, we will always recall the connection between the driving coupling constant $\epsilon$ and the input driving field amplitude $A_{\rm in}$ to avoid any confusion.

It is clearly seen from Eq. (III.2) that the rate equations simplify to a linear form if the gain saturation in the active laser cavity are neglected.

We note that in the semiclassical Maxwell-Bloch picture, which is widely used for describing the characteristic $\cal PT$-symmetric properties of classical optical fields, one can introduce gain saturation via a modified gain coefficient $A$ for the field $A_{1}$ in the active cavity as $A\rightarrow A/(1+|A_{1}|^{2}/|A_{\rm s}|^{2})$ in Eq. (III.2) for $B=0$, where $|A_{\rm s}|^{2}$ is the gain saturation threshold Chang et al. (2014). In this case, by decomposing this modified gain coefficient $A$ in a Taylor series (up to the first order in $|A_{1}|^{2}/|A_{\rm s}|^{2}\ll 1$, which is justified in the weakly saturated regime), one obtains a one-to-one correspondence between the semiclassical Maxwell-Bloch picture and the semiclassical Scully-Lamb laser model for the optical fields in Eq. (III.2) by setting $B=A/|A_{\rm s}|^{2}$. Nevertheless, compared to the Maxwell-Bloch picture, where the gain saturation term is introduced phenomenologically, the semiclassical Scully-Lamb laser theory explains gain saturation not only qualitatively, but also quantitatively, via the laser system parameters given in Eq. (3).

We also note that a set of nonlinear equations, similar to Eq. (III.2), can be obtained within the semiclassical Maxwell-Bloch approximation for $\cal PT$-symmetric coupled waveguides Ramezani et al. (2010). Nonetheless, there are two main differences with respect to Ref. Ramezani et al. (2010). The first is that we consider the case of coupled resonators instead of coupled waveguides. In our case, feedback effects and two different kinds of losses (i.e., input-output and intrinsic losses) are considered. Second, even if we apply the semiclassical approximation, our theory is based on a fully quantum approach, which can be useful to investigate quantum effects and quantum noise.

III.3 Eigenfrequencies in the steady-state and exceptional point

In the steady state, by considering $\gamma_{1}=\gamma_{2}=\gamma$ for symmetry reasons, we find from Eq. (III.2) that:

[TABLE]

where $G^{\prime}_{1}=G_{1}-B|A_{1}|^{2}$. Because the system of equations in Eq. (III.3) is written for the complex fields $A_{j}$, we can incorporate the real steady-state intensity $|A_{1}|^{2}$ into the prefactor $G_{1}^{\prime}$. One can obtain the eigenvalues of the system, given in Eq. (III.3), by setting the driving field to zero ($\epsilon=0$). Defining the vector $\alpha=(A_{1},A_{2})^{T}$, we can rewrite Eq. (III.3) as follows

[TABLE]

with the evolution matrix

[TABLE]

For the case when gain and losses are balanced, i.e., when $G_{1}^{\prime}=\Gamma_{2}$ holds, the evolution matrix $M$, given in Eq. (13), becomes $\cal PT$-symmetric, i.e., $[M,{\cal PT}]=0$, with the parity operator $\cal P=\begin{pmatrix}0&1\\ 1&0\end{pmatrix}$ and $\cal T$ performing complex conjugation.

The characteristic equation of the matrix $M$ is

[TABLE]

from which we find the formal solution for the eigenfrequencies, as follows

[TABLE]

By recalling that $\Delta_{\pm}=\omega_{\pm}-\omega_{c}$, we obtain

[TABLE]

where $I_{1}=|A_{1}|^{2}$ is the dimensionless intensity of the field in the active cavity in the steady-state, and which is by itslef a function of the frequency $\omega$ (for the solutions of the steady-state intensity $I_{1}$, see Section III E). The solutions in Eq. (16) give the energy eigenspectra of the system, i.e., the eigenvalues of the system Hamiltonian, which determine the evolution matrix $M$.

We note again that, in reality, the loss rates $\gamma_{i}$, arising from the coupling with the input-output channels, are not true losses, because these describe the energy transfer from the system to the output (or from the input to the system). Hence, the concept of the effective ${\cal PT}$-symmetry in our system with balanced gain and loss, can be expressed as

[TABLE]

This is valid because the system is ${\cal PT}$-symmetric regardless of how it is probed (i.e., waveguides in the coupled resonator systems are used only to probe the system). When the condition, given in Eq. (17), is satisfied we find from Eq. (16) that:

[TABLE]

Note that, when changing the input signal $\epsilon$ (and hence the resulting steady-state intensity $I_{1}$), one has to adjust the losses correspondingly in order to satisfy the ${\cal PT}$-symmetry condition.

The analysis of the frequency spectrum, given in Eqs. (16) and (18), provides two different regimes, or so-called unbroken and broken $\cal PT$-symmetry phases, depending on the sign of the expression under the square-root sign. In the unbroken $\cal PT$-symmetry phase, the expression under the root is positive (that is, $\kappa>\frac{C_{2}}{2}$), and there are always two supermodes with non-degenerate real frequencies $\omega_{\pm}$ that propagate in the system. Note that this is true for the system itself, where the coupling loss $\gamma$ (which is not inherent to the system of coupled resonators) is zero (i.e., $\gamma=0$). In the broken $\cal PT$-symmetry phase, that expression is negative (that is, $\kappa<C_{2}/2$), and the real spectrum becomes degenerate, indicating that the system displays two modes with the same resonance frequency but with different decay rates. The transition between these two regimes takes place at an EP given by

[TABLE]

In Fig. 2 we show the real and imaginary parts of the eigenfrequencies of the supermodes. As expected by the inspection of Eq. (19), the EP at $\kappa=\kappa_{\rm EP}$ does not depend on the field intensity. For $\kappa>\kappa_{\rm EP}$, we observe that the imaginary part is different from zero. This is due to the contributions of the loss rates $\gamma_{i}$ arising from the coupling with the input-output channels. The gray dash-dotted curve in Fig. 2 describes the imaginary parts of the complex eigenfrequencies, when the input-output coupling losses are neglected (i.e., $\gamma=0$).

We now consider the case, where the ${\cal PT}$-symmetry condition is achieved at low input rates, so the gain saturation effects are negligible ($BI_{1}\simeq 0$). In this case, as for the ${\cal PT}$-symmetry condition, we can use Eq. (17) calculated for $I_{1}=0$:

[TABLE]

If the input drive is increased without adjusting the other parameters ($\kappa$ and $C_{i}$), the ${\cal PT}$-symmetry condition is not satisfied any more. We will investigate the effect of gain saturation on the spectral properties of the considered ${\cal PT}$-symmetric system, which is usually studied without the inclusion of gain saturation. Thus, hereafter, we assume that the ${\cal PT}$-symmetric condition is given in Eq. (20). The condition for the exceptional point becomes [see Eq. (16)]:

[TABLE]

From inspection of Eq. (21), we observe that, in this case, the EP changes when the steady-state field intensity increases in the active cavity (see Fig. 3). In particular, as an example, in Fig. 3 we plot the real and imaginary parts of the eigenfrequencies of the supermodes for the linear [Fig. 3(a)] and nonlinear [Figs. 3(b)–(d)] regimes for different values of the driving-field coupling $\epsilon$ when driving the system at a resonant frequency $\omega_{c}$.

In the linear regime, the gain saturation term $BI_{1}$ is either exactly zero or negligible compared to the system parameters $A$ and $C_{i}$; thus, the behavior of the EP is identical to that presented in Fig. 2 [see also Figs. 3(a) and 4(a)]. The same conclusion applies when the driving field is far away from the resonance $\omega_{c}$, as in that case, the steady-state intensity also tends to zero.

An interesting situation arises when the gain saturation term $BI_{1}$ becomes comparable with the passive cavity loss $C_{2}$. However, at the same time, it is much less than the gain coefficient $A$ in the active cavity, i.e., $BI_{1}\ll A$, so the weakly saturated regime still holds and the validity of the master equation in Eq. (2) remains. This also implies that $C_{1}\gg C_{2}$. In this case, the system starts exhibiting some nonlinear features in its eigenspectrum. In short, the described condition can be written as

[TABLE]

In what follows, we always call the system to be in the nonlinear ${\cal PT}$-symmetric regime, whenever both conditions, given in Eqs. (20) and (22), are satisfied. In that case, one can observe that the critical value of $\kappa_{\rm EP}$ significantly changes depending on the gain saturation term $BI_{1}$. Moreover, the steady-state intensity $I_{1}$ in the active cavity by itself becomes dependent on the direction of the propagation of the driving field, i.e., on whether the driving field is coupled to the active (from port 1 to port 4) or the passive cavity (from port 4 to port 1) [see Fig. 1]. Consequently, the gain saturation term $BI_{1}$ also depends on the driving-field direction.

When the input driving field is coupled to the active cavity, it can experience a significant gain saturation for large input intensities, and, as a result, its losses increase. This especially happens when one decreases the intercavity coupling $\kappa$. In the latter case, the strong signal field becomes localized in the active resonator, before being transferred to the lossy passive cavity. This leads to the gain decrease of the intense driving field due to gain saturation, and, therefore, to the observed signal-field losses (see blue dashed curve in Fig. 3). Moreover, when the losses induced by the gain saturation become comparable to the losses in the passive cavity (i.e., when $BI_{1}\approx C_{2}$), then the critical value of $\kappa_{\rm EP}$ first decreases and then increases when increasing the intensity of the input field (see Fig. 3). As Eq. (21) implies, this shift of $\kappa_{\rm EP}$ can be explained by the interplay between losses in the passive cavity and losses induced by gain saturation in the active resonator. Namely, by increasing the gain saturation term $BI_{1}$ (i.e., by increasing the input signal) in the proximity of $C_{2}$ in Eq. (21), the critical value $\kappa_{\rm EP}$ first decreases, when $BI_{1}\leq C_{2}$, and then increases, when $BI_{1}\geq C_{2}$.

III.4 Non-reciprocity of light propagation

Another important result of our work is the theoretical microscopic prediction of non-reciprocity of the propagating light. Specifically, in broken $\cal PT$-symmetry phase, i.e., when $\kappa<\kappa_{\rm EP}$, the signal field can experience smaller losses for smaller values of $\kappa$, when the driving field propagates in the direction $4\rightarrow 1$, in comparison to the case when it propagates in the opposite direction $1\rightarrow 4$ (see Fig. 4). This stems from the fact, that when the strength of the input signal field is increasing, by decreasing the intercavity coupling $\kappa$, the driving field experiences large losses in the passive cavity before passing into the active resonator. Now, because the initially strong input signal field is strongly damped by the passive resonator, it enters the active cavity having an intensity which is not sufficient to induce gain saturation. As such, the propagating field in the active cavity can even undergo notable amplification before being detected at port 1. As a result, the critical values of $\kappa_{\rm EP}$ only decrease when increasing the intensity of the input signal field, since the propagating field cannot reach high intensity in the active cavity to make the term $BI_{1}$ larger than $C_{2}$ in Eq. (21) (see also Fig. 4).

We plotted Figs. 5–7 to demonstrate the aforementioned asymmetry of the steady-state intensities in both cavities depending on the propagation direction of the resonant signal field. As it follows from Fig. 7, for very small values of $\kappa<\kappa_{\rm EP}$, the field intensity in the active cavity for the direction $4\rightarrow 1$ can be two or three orders of magnitude larger than the intensity in the passive cavity for the opposite propagation direction $1\rightarrow 4$. As a consequence, this asymmetric property can lead to the observation of non-reciprocal light behavior for small $\kappa$. Therefore, the non-reciprocity of light propagation, as experimentally demonstrated in, e.g., Ref. Peng et al. (2014a), arises from the combination of loss, gain, and gain saturation. Specifically, owing to the saturation effects, the gain experienced by the signal entering the amplifying cavity strongly depends on its intensity. We would like to stress again that the input signal with a large intensity undergoes a much lower amplification with respect to a weaker signal due to gain saturation. Hence, feeding the gain cavity with a quite strong field gives rise to a modest amplification (see also Fig. 8). This signal is finally strongly absorbed by the lossy cavity, before being detected. By contrast to this, if the same signal is first sent to the lossy cavity, then it is strongly absorbed before entering the gain cavity. Such a small signal does not saturate the gain medium and is strongly amplified before detection.

As Figs. 3 and 4 also indicate, in unbroken $\cal PT$-symmetry phase, for large values of the intercavity coupling $\kappa\gg\kappa_{\rm EP}$, the system exhibits a linear character, regardless of both strengths of the input signal field and the propagation direction ($1\rightarrow 4$ and $4\rightarrow 1$). Indeed, in this case, the coupling $\kappa$ between two microcavities becomes large enough, enabling the input signal fields to freely propagate in either direction. Hence, there is no localization of the fields in the system, and, thus, there is no observed nonlinearity due to gain saturation.

III.5 Transmission spectra

Here we focus on the spectral properties of the driving fields that propagate through the system.

By rewriting the complex amplitudes $A_{k}$ of the fields [as in Eq. (III.2)] as $A_{k}=|A_{k}|\exp({i\phi_{k}})$, one arrives at a cubic equation for the field intensity $I_{1}$ in the active cavity in the steady state (see Appendix B, for details):

[TABLE]

with coefficients $\lambda_{i}$ defined as:

[TABLE]

where $f=4\kappa^{2}/(\Gamma_{2}^{2}+4\Delta^{2})$. Equation (23) has only one real solution (see Appendix B, for details), when its discriminant is negative, which is always the case when, e.g., $A\approx\Gamma_{1}$ regardless of $\Gamma_{2}$ and $\kappa$.

The transmission spectrum can be calculated as follows

[TABLE]

To obtain the transmission spectrum $T_{1\rightarrow 4}(\omega)$ at port 4 when sending the signal from port 1, one needs to know the expressions for the corresponding input-output fields. The output field at port 4 is found as $A_{\rm out}=\sqrt{\gamma_{2}}A_{2}$. The input driving field sent from port 1 can be expressed, as $A_{\rm in}=\epsilon/\sqrt{\gamma_{1}}$, where, again, $\epsilon$ is the coupling strength between input driving field and microresonator (Fig. 1). By rewriting the field $A_{2}$ via $A_{1}$ and using Eq. (III.2), one finally obtains

[TABLE]

The same analysis can be carried out for the case when the driving coherent field is sent from port 4, and the signal is detected at port 1. In that case, one obtains the same cubic equation for the field intensity $I_{1}$, as in Eq. (23), but with different $\lambda_{k}$; $k=1,\dots,4$ (see Appendix B). The expression for the transmission spectra $T_{4\rightarrow 1}(\omega)$ then attains the following form

[TABLE]

Similarly, the transmission spectrum $T_{1\rightarrow 2}$ can be found using the input-output relation: $A_{\rm out}=A_{\rm in}+\sqrt{\gamma_{1}}A_{1}$ (see Appendix B, for details).

As an example, in Fig. 9 we plot the normalized transmission spectrum $\tilde{T}_{1\rightarrow 4}=T_{1\rightarrow 4}(\Delta)/{\rm max}[T_{1\rightarrow 4}(\Delta)]$ for the ${\cal PT}$-symmetric condition $A-C_{1}-C_{2}=0$, for different values of the intercavity coupling $\kappa$, and for different intensities of the driving field (i.e., by varying $\epsilon$) in both linear ($BI_{1}/C_{2}\ll 1$) and nonlinear ($BI_{1}/C_{2}\approx 1$) regimes [for details regarding the nonlinearity condition, see Eq. (22)]. One can observe an increasing spectral line broadening of the transmitted light for stronger driving fields $\epsilon$ in the nonlinear regime, indicating rising losses due to the gain saturation for small values of $\kappa$ [green dashed and blue dash-dotted curves in Fig. 9(a)]. Moreover, the transmission spectra in Figs. 9(b) and 9(c) display a squared Lorentzian line shape, which is characteristic for EPs Yoo et al. (2011); Sweeney et al. (2018); Pick et al. (2017). Also, in the nonlinear case, the splitting of the supermodes occurs for larger $\kappa>\kappa_{\rm EP}$ and more intensive fields [Figs. 9(b), 9(c)]. A similar behavior is observed in the transmission spectrum $T_{4\rightarrow 1}$.

IV A comparison with experimental results of

Ref. Peng et al. (2014a)

For simplicity, when plotting the graphs in this section, we assume that the losses in both cavities are comprised mainly by the losses due to the coupling of the cavities to the waveguides $\Gamma_{i}=C_{i}+\gamma_{i}=\gamma_{i}$, $i=1,2$, i.e., we set the intrinsic losses $C_{i}$ to zero. The latter assumption also implies that the system considered has a broken ${\cal PT}$-symmetry, i.e., $A-C_{1}-C_{2}\neq 0$, according to Eq. (20). At the same time, the total losses in the system are expected to be larger than the gain, i.e., $(A-\gamma_{1}-\gamma_{2})<0$. Also, the active microcavity is assumed to operate near the threshold $A\approx\gamma_{1}>\gamma_{2}$.

In this section, we discuss possible applications of the semiclassical Scully-Lamb laser theory in the prediction of some nontrivial light behavior that was experimentally observed in Refs. Peng et al. (2014a); Chang et al. (2014). In those papers, the authors experimentally studied a system of coupled ${\cal PT}$-symmetric whispering-gallery microcavities, i.e., a system that is identical to that presented in Fig. 1 and which is the focus of the theoretical study of this work. Below, we theoretically reproduce some of the experimental graphs of Ref. Peng et al. (2014a) in a qualitative rather than quantitative way (meaning that we make some additional assumptions regarding the system parameters, used in constructing the graphs here). Nevertheless, as was just mentioned, a qualitative comparison can be made, and positive conclusions can be inferred regarding the applicability of the Scully-Lamb laser model, in its semiclassical limit, to explain some of the results of Ref. Peng et al. (2014a). We note that in the construction of the graphs shown here, which are presented in Figs. 10–13, we do not invoke ${\cal PT}$-symmetry in the system (see also the text below).

For clarity, in the captions of some of our figures we indicate the corresponding experimental plots of Ref. Peng et al. (2014a) that we try to theoretically reproduce. Also, in order to stress the similarity between the figures reproduced here and the original experimental graphs of Ref. Peng et al. (2014a), we keep the axis scales of the plots to be the same as those given in Ref. Peng et al. (2014a).

For example, in Fig. 10, which corresponds to the experimentally-obtained Fig. 1(f) in Ref. Peng et al. (2014a), the transmission spectrum $T_{1\rightarrow 2}$ is shown when only the WG1 and the microresonator $R_{1}$ are coupled in the system and $\kappa=0$, i.e., only the first resonator is considered. Thus, by coupling the signal sent from port 1 to the active cavity, one obtains a substantial signal amplification, when detecting the output signal at port 2. In this case, the active cavity just enhances the incoming field. The inset of Fig. 10 demonstrates that amplification does not occur if there is no input signal.

In Fig. 11, which corresponds to the experimentally-obtained Fig. 3 in Ref. Peng et al. (2014a), we show the transmission spectra for the light propagating in the directions $1\rightarrow 4$ and $4\rightarrow 1$ in the linear regime, i.e., the laser cavity is assumed to be below the lasing threshold $A<\gamma_{1}$. One can see the linear response of the propagating signal when the inter-cavity coupling coefficient $\kappa$ is lower or larger than $\kappa_{\rm cr}$, where $\kappa_{\rm cr}$ denotes the critical point, when the supermodes start splitting. We note that, in Fig. 11, a linear behavior of the transmitted light is already observed in the unbalanced gain-loss regime, i.e., when ${\cal PT}$-symmetry is broken in the system.

Figure 12, which corresponds to the experimentally-obtained Fig. S6 in Ref. Peng et al. (2014a), displays the transmitted spectrum $T_{1\rightarrow 2}$ versus the detuning $\Delta=\omega-\omega_{c}$, for various values of the gain $A$ and the intercavity coupling $\kappa$. One can see the appearance of the supermodes splitting in the system with the increasing values of $\kappa$. Conversely, for smaller values of $\kappa$, the two supermodes coalesce resulting in only one peak in the spectrum.

Non-reciprocity of light propagation

Here, to complement Sec. III D, we further discuss the theoretical microscopic prediction of non-reciprocity of the propagating light in the considered coupled active-passive microresonators system as shown Fig. 13, which corresponds to the experimentally-obtained Fig. 4 in Ref. Peng et al. (2014a). It is seen that there is an enhancement in the transmitted light from port 4 to port 1, and tending to zero transmission $T_{1\rightarrow 4}$ in the opposite direction for small values of $\kappa$. Meaning that the system starts behaving nonreciprocally. The latter nonlinear effect was observed in Ref. Peng et al. (2014a), but, naturally, could not be explained based on the linear rate equations used there. Utilizing the semiclassical laser theory, one can attain the needed nonlinear term arising from the laser gain saturation in the active microcavity. Moreover, this non-reciprocity is observed without invoking ${\cal PT}$-symmetry, because it can be already observed when the gain and loss are unbalanced in the system (see Fig. 13). The inset of Fig. 13(f) indicates that the observed nonlinearity is not caused by the lasing initiated by spontaneous emission in the active cavity.

V Conclusions

We have applied the quantum Scully-Lamb laser theory to a pair of ${\cal PT}$-symmetric coupled whispering-gallery microcavities, i.e., a system, which consists of both active and passive microring cavities, such that gain and losses are balanced in the system. It has been shown that, in the nonlinear regime, or more precisely, under the condition in which the gain saturation in the active cavity is comparable to the losses in the passive cavity, the intense intracavity fields of the steady state lead to the modification of the eigenmodes and of the EPs of the ${\cal PT}$-symmetric system. Namely, the imaginary part of the eigenspectrum acquires an extra negative term due to the gain saturation effects. This effect leads to the shift of the EP either to lower or larger values depending on the gain saturation $B$ and the propagation direction of the driving fields. Starting from the master equation for this coupled system, including dissipation, gain, and gain saturation, and applying the semiclassical approximation, we are able to describe the experimental results obtained in Refs. Peng et al. (2014a); Chang et al. (2014). In particular, this approach is able to reproduce the observed non-reciprocal light propagation in the coupled system of whispering-gallery microcavities. We have also shown that the gain saturation mechanism in the active cavity is crucial for the observation of light non-reciprocity. Moreover, we have found that the unidirectional light propagation can be observed even when the ${\cal PT}$-symmetry condition is not fulfilled. It should be stressed that neither ${\cal PT}$-symmetry nor its breaking is required for nonreciprocity. The nonreciprocity observed in our system is a result of a nonlinearity and, in the broken ${\cal PT}$-regime, such nonlinearity can be observed at much lower input intensity.

In summary, we proposed, applied, and validated the Scully-Lamb laser model, with the non-Lindbladian master equation, for coupled resonators with losses, gain, and gain saturation. Although we studied non-reciprocity and exceptional points applying the semiclassical approximation, this master equation allows for a quantum description of the cavity fields. This approach constitutes a promising tool for the study of quantum optical effects in coupled resonators with balanced (or unbalanced) gain and losses.

Acknowledgements.

I.I.A. was supported by Grant Agency of the Czech Republic (Project No. 17-23005Y), the Project CZ.02.1.010.00.016_0190000754, and the Project LO1305 of the Ministry of Education, Youth and Sports of the Czech Republic. A.M. and F.N. acknowledge the support of a grant from the John Templeton Foundation. F.N. is supported in part by the: MURI Center for Dynamic Magneto-Optics via the Air Force Office of Scientific Research (AFOSR) (FA9550-14-1-0040), Army Research Office (ARO) (Grant No. W911NF-18-1-0358), Asian Office of Aerospace Research and Development (AOARD) (Grant No. FA2386-18-1-4045), Japan Science and Technology Agency (JST) (the Q-LEAP program, and CREST Grant No. JPMJCR1676), Japan Society for the Promotion of Science (JSPS) (JSPS-RFBR Grant No. 17-52-50023), and the RIKEN-AIST Challenge Research Fund. S.S. is supported in part by the Army Research Office (ARO) (Grant No. W911NF1910065). S.K.O. is supported by Army Research Office (ARO) grant No. W911NF-18-1-0043, Air Force Office of Scientific Research (AFOSR) award No. FA9550-18-1-0235, National Science Foundation (NSF) (1807485), and by The Pennsylvania State University, Materials Research Institute (MRI).

Appendix A Derivation of the master equation in Eq. (2)

To make this article self-consistent, in this Appendix, we present a derivation of the quantum laser master equation, given in Eq. (2), for the Scully-Lamb laser model based on the derivation of Yamamoto and Imamoǧlu in Ref. Yamamoto and Imamoǧlu (1999). This derivation bears a phenomenological character and, as such, naturally allows to incorporate all the terms of the interaction Hamiltonian of the fields without the need to solve the Schrödinger equation directly. Another derivation of the master equation (2) can be found in Ref. Orszag (2008).

The active cavity is represented by a general four-level laser system in which the two intermediate energy levels are coupled by the laser mode (see Fig. 14). In this limit, the uppermost level of the atom may be adiabatically eliminated to give an effective three-level system. The latter assumption is valid as long as the decay rate from the uppermost state $|1\rangle$ to the upper laser level $|e\rangle$ is much faster than all other rates in the atom-field system. In this limit, one has an effective incoherent pumping rate $r$ from the atomic ground state $|0\rangle$ into $|e\rangle$. Additionally, the laser mode in the active cavity is coupled to the passive microresonator.

The interaction Hamiltonian $\hat{H}_{I}$, describing such two coupled active-passive microresonators system, is given by

[TABLE]

Specifically, the second term in Eq. (28) describes the linear interaction between the modes in the active ($a_{1}$) and passive ($a_{2}$) resonators; while the first term describes the interaction between the mode $a_{1}$ with only two levels ($|g\rangle$ and $|e\rangle$) within the standard Jaynes-Cummings model. The operators $\hat{\sigma}_{k}$, $k=+,-$, are spin-raising and spin-lowering operators of the atom in the active medium, respectively. The constants $g$, $\kappa$ denote the coupling strength between the atom and the field in the laser cavity, and between the two fields propagating in two different cavities, respectively. We also assumed that the atomic and active cavity field resonances coincide.

The quantum Liouville equation for the density operator $\hat{\rho}$ in the interaction picture, which describes the atom-field-field dynamics is:

[TABLE]

In the active cavity, the optical gain is provided by the excited atoms, which are pumped by an external field. The interaction between the atom and an external pumping field, that provides an inverse population in our effective three-level atom laser (depicted in Fig. 14), can be described by the following equation

[TABLE]

where $\hat{\sigma}_{00}=\hat{\sigma}_{0e}\hat{\sigma}_{e0}$, and $\hat{\sigma}_{0e}$ ($\hat{\sigma}_{e0}$) is the spin operator for the atomic transition from $|0\rangle$ ($|e\rangle$) to $|e\rangle$ ($|0\rangle$). The coefficient $r$ accounts for the pumping rate of the atom.

To include spontaneous emission and the emission caused by external dephasing processes into an external reservoir field, we need to add in Eq. (29) the following terms:

[TABLE]

where $\gamma_{\rm sp}$ and $\gamma_{d}$, are, respectively, the rates of spontaneous emission, and the emission imposed by additional dephasing processes.

Now, collecting together Eqs. (29)–(31), one arrives at the master equation for the density operator $\hat{\rho}$

[TABLE]

where we introduced a Lindbladian damping super operator as $\hat{L}^{d}_{i}(\hat{O})=\hat{a}_{i}^{\dagger}\hat{a}_{i}\hat{O}+\hat{O}\hat{a}_{i}^{\dagger}\hat{a}_{i}-2\hat{a}_{i}\hat{O}\hat{a}_{i}^{\dagger}$.

To obtain the master equation for the reduced field density operator $\hat{\rho}_{f}$, which describes the dynamics of the optical fields in the cavities, one has to trace out the density operator $\hat{\rho}$ in Eq. (32) over the atom states. Namely,

[TABLE]

[TABLE]

where the Hamiltonian $H_{f}$ accounts for the field interaction between active and passive cavities, i.e., it is the second term in Eq. (28). The operator \hat{\rho}_{ge}=\left\langle g\Big{|}\frac{d}{dt}\hat{\rho}\Big{|}e\right\rangle, from Eq. (32), obeys the following equation

[TABLE]

where $\gamma_{T}=\gamma_{\rm sp}+\gamma_{d}$ is the total decay rate of the atom. The same relation, given in Eq. (35), also holds true for the operator $\hat{\rho}_{eg}=\langle e|\hat{\rho}|g\rangle$.

Assuming $\gamma_{T}\gg C_{1},C_{2},\kappa$, one eliminates $\hat{\rho}_{ge}$ and $\hat{\rho}_{eg}$ in the adiabatic approximation, i.e., $\frac{d}{dt}\hat{\rho}_{ge}=0$. Thus, one obtains

[TABLE]

Substituting $\hat{\rho}_{ge}$ and $\hat{\rho}_{eg}$ in Eq. (36) into Eq. (34) one obtains

[TABLE]

where $L_{i}^{a}(\hat{O})=\hat{a}_{i}\hat{a}_{i}^{\dagger}\hat{O}+\hat{O}\hat{a}_{i}\hat{a}_{i}^{\dagger}-2\hat{a}_{i}^{\dagger}\hat{O}\hat{a}_{i}$ is the Lindbladian amplification super operator.

Moreover, $\hat{\rho}_{ee}$ and $\hat{\rho}_{gg}$ satisfy the following equations

[TABLE]

If the population of the lower energy level $|g\rangle$ is very low, i.e., the energy quickly decays into the ground state $|0\rangle$, and if the gain saturation is weak, then one may write

[TABLE]

Applying standard perturbation techniques to Eq. (A), one arrives at

[TABLE]

Combining now Eqs. (39), (40), and Eq. (37), one attains the master equation for the two field operator in the interaction picture

[TABLE]

with

[TABLE]

where the coefficient $A$ stands for the linear gain, and $B$ is the gain saturation coefficient.

For the case of the ideal laser system with $\gamma_{d}=\gamma_{sp}=0$, one can attain the quantum laser master equation within the Scully-Lamb laser theory. In the weakly saturated regime, the spontaneous emission on the laser transition can be discarded, and one obtains Yamamoto and Imamoǧlu (1999)

[TABLE]

where $\Gamma$ is a total decay rate of the states $|g\rangle$ and $|e\rangle$ into the ground state $|0\rangle$ of the atom. Substituting equations in Eqs. (A) into Eq. (37), we obtain

[TABLE]

where the gain and gain saturation coefficients are now expressed as

[TABLE]

Note a typo in the prefactor of the gain saturation coefficient $B$ in Eq. (45) in Ref. Yamamoto and Imamoǧlu (1999). Namely, there is a prefactor 1, instead of 4.

We may also add the Hamiltonian term related to the coupling between the cavity field $\hat{a}_{1}$ and the external driving coherent classical field into the master equation. Such a Hamiltonian can be given by

[TABLE]

where the coupling constant $\epsilon\equiv\sqrt{\gamma_{1}{P}/(\hbar\omega_{l})}$ accounts for the coupling between the driving external coherent field with power $\cal P$ and the cavity field $\hat{a}_{1}$.

By rewriting the field Hamiltonian $\hat{H}_{f}$ in the Schrödinger picture, we finally arrive at Eq. (2).

Appendix B Derivation of Eq. (23), and some exact formulas from Sec. III E

Working in the reference frame where the phase of the driving field is zero, and by expressing the complex amplitudes in the rate equations in Eq. (III.2) as $A_{k}=|A_{k}|e^{i\phi_{k}}$, one obtains the following equations for the steady-state:

[TABLE]

Replacing now $|A_{2}|\exp[i(\phi_{2}-\phi_{1})]$ by $|A_{1}|$ in the second equation of Eq. (47) and inserting it into the first equation, we attain

[TABLE]

Separating the real and imaginary parts in Eq. (48) and equalizing them to zero, one arrives at

[TABLE]

where $f$ is given in Eq. (III.5). Utilizing now the standard trigonometric relation $\cos^{2}\phi_{1}+\sin^{2}\phi_{1}=1$, and collecting together the coefficients at each order of the real amplitude, we finally obtain the cubic equation for the field intensity $I_{1}=|A_{1}|^{2}$, given in Eq. (23).

The unique real solution of the steady-state intensity $I_{1}$, given in Eq. (23), is

[TABLE]

where

[TABLE]

with the coefficients $\lambda_{k}$ introduced in Eq. (III.5). For the transmission spectra $T_{4\rightarrow 1}$, given in Eq. (27), the steady-state intensity $I_{1}$ has the same solution as in Eq. (50), and with the same coefficients $\lambda_{k}$ in Eq. (III.5), except the coefficient $\lambda_{4}$, which reads now as

[TABLE]

By combining together the input-output relation $A_{\rm out}=A_{\rm in}+\sqrt{\gamma_{1}}A_{1}$ with the solution for the phase of the complex amplitude $A_{1}$ in Eq. (B), one can straightforwardly find the transmission spectrum $T_{1\rightarrow 2}$, after applying Eq. (25), as follows

[TABLE]

where $F$ is given in Eq. (III.5).