Spontaneous PT-symmetry breaking in lasing dynamics

Parity-time (PT) symmetric lasers exploit the modulation of optical gain and loss and have led to important fundamental demonstrations in non-Hermitian physics. The current theoretical analysis of PT-symmetric laser physics is performed on the basis of the adiabatic elimination of the medium polarization. This approximation doesn’t hold true for a more general optical system with strong photon-particle interactions, where the Rabi oscillation of active particles plays a non-negligible role in the lasing action. Here, we propose a model that takes into account the internal dynamics of active particles and numerically investigate the PT symmetry of macroscopicand microscopic-sized laser systems that operate in the strongcoupling regime. The distinct phase diagrams are drawn according to the features of intracavity photon numbers and emission spectra. Our work extends the PT-symmetric optics from the weakto the strong-coupling limit, potentially paving the way towards nonclassical PT-symmetric light sources for integrated photonic networks and ultrasensitive sensors. https://doi.org/10.1038/s42005-021-00575-7 OPEN

I n quantum mechanics, the reality of a physical observable and the unitary dynamics for a quantum system rely on the fundamental postulate of the Hermiticity of the associated operators. However, this postulate has been questioned by the recent discovery that a non-Hermitian Hamiltonian with PT symmetry may also possess a real energy spectrum 1 . Due to their high degree of flexibility and controllability, gain-loss-balanced optical systems provide a versatile platform to exploit the fundamental nature of the non-Hermitian PT-symmetric quantum mechanics in experiment 2 . Thus far, a number of PT-symmetryrelated phenomena, including spontaneous PT-symmetry breaking 3 , unidirectional invisibility 4 , and nonreciprocity 5 , have been observed for the propagation of light waves in a pair of coupled optical devices, such as waveguides and microcavities. Besides, the PT-symmetric optical platforms can also be employed in single dielectric nanoparticle/molecule detection, where the micrometer-sized sensors operate around the so-called exceptional points (EPs), potentially enhancing the sensitivity [6][7][8] . In all these studies, two coupled optical waves (modes) follow the Schrödinger-type equation 9 i d dt with the amplitudes A k¼1;2 ðtÞ ¼ ffiffiffiffiffiffiffiffiffiffi n k ðtÞ p exp½iθ k ðtÞ, intensities (intracavity photon numbers) n k = 1,2 (t) and phases θ k = 1,2 (t) of the light fields. The two-dimensional Hamiltonian is given by with the complex optical gain/loss parameter ξ and the real interwave/intermode coupling strength g C > 0. Obviously, H is non-Hermitian, H ≠ H † , and owns the PT symmetry, PTH = HPT with the parity operator P = σ x and the time-reversal operator T = K. Here, σ x is the x-component of Pauli matrices and K denotes the complex conjugate operator. Comparing |Im ξ| and g C , the coupled optical system exhibits two distinct phases. In the unbroken phase with |Im ξ| > g C , the energy-level spectrum of H is entirely real-valued and the powers of two optical waves/modes and their power sum present an oscillatory behavior. By contrast, the eigenvalues of H become a pair of complex conjugates in the broken phase with |Im ξ| > g C , where the amplitude of one supermode grows exponentially in time while the other decays. Generally, the PT-symmetric optical platforms are divided into two groups according to the manner of their operation. In the passive group, the optical waves in the coupled system come from the external light sources. The transmitted beams from the system are measured to analyze the energy-level spectrum of H, where the locations and linewidths of the transmission peaks in the frequency domain correspond to the real and imaginary parts of the eigenvalues of H, respectively. However, the extra probe-beam noise inevitably influences the spectral purity of the transmission response. Furthermore, the transmission spectra cannot unveil some intrinsic properties of the coupled optical system (e.g., quantum correlation and statistics of photon emission) that are determined by the physical realization itself. In the active group [10][11][12][13][14] , the coherent light waves are generated through the lasing action of the modes of interest and the intermode coupling is superimposed on the laser dynamics. The emission spectra are utilized to map the properties of the coupled system. Since no external light fields are introduced, these active platforms may reveal the nature of the PT-symmetric optics, compared to the passive ones. Thus far, the PT-symmetry breaking concept has been applied to interpret most experimental observations of various active optical schemes. A direct proof showing that the coupled laser system owns the PT symmetry from the measurement of the output light fields is still demanded. We notice that the emission spectra from a pair of PT-symmetric optical resonators were studied in theory 15 . However, the relevant optical scheme operates below or at the lasing threshold, rather than the lasing regime.
In addition, the current PT-symmetric optics is usually performed on the basis of the adiabatic elimination of the medium polarization [16][17][18][19][20] . This adiabatic approximation is valid when the relaxation constant of the medium polarization much exceeds the cavity loss rates (i.e., the good-cavity limit) and the Rabi frequency associated with the light-particle interaction (i.e., the weak-coupling limit). The resulting rate equations for the coupled-cavity fields can be mapped onto the Schrödinger Eq. (1) with the Hamiltonian (2). Indeed, the rate-equation model is commonly employed to study, for example, the semiconductor lasers 10,11 but fails in describing the laser systems that operate in the strong light-particle coupling limit, such as the atomic gas/ beam lasers [21][22][23][24][25] . Hence, a more general study is necessary. Moreover, thus far, research has been carried out mainly on the classical optical systems. Suppressing the system's size, which is measured by the numbers of particles and light quanta (photons), highlights the quantum behavior of the light-particle interaction and a full quantum treatment is needed 26 . However, rare studies have been focused on the PT-symmetric optics in this regime.
In this paper, we numerically investigate the PT symmetry of an open system composed of a macroscopic/microscopic laser interacting with an empty optical cavity. In the model, we take into account the strong light-particle interaction, which makes the adiabatic elimination of the medium polarization become invalid. Our goal is to identify the PT-symmetry breaking through the measurements of the intracavity photon numbers and the emission spectra of the cavity fields. The macroscopic-sized system presents two steady-state lasing phases, where one denotes the PTsymmetry broken phase, while the other is non-PT-symmetric. The system's phase diagram has one unstable lasing regime, which corresponds to the PT-symmetry unbroken phase, and one nonlasing regime. By contrast, the microscopic-sized system exhibits the steady-state PT-symmetry broken and unbroken phases. The phase transition is controlled via adjusting the intercavity coupling strength. Especially, the intensity auto-and cross-correlation functions of two cavity fields display the nonclassical behaviors of the photon emission of individual optical cavities and photon-pair generation. This microscopic-sized PT-symmetric lasing scheme has, to the best of our knowledge, not been studied before and may potentially operate as a nonclassical PT-symmetric light source in an integrated optical circuit. Our study explores the fundamental properties of PT-symmetric optics in the lasing regime, and our predictions can be tested by current microcavity and nanophotonic technologies.

Results and discussion
Macroscopic laser oscillation. Figure 1a, b illustrate the macroscopic and microscopic optical systems, respectively, under study. We first consider the classical laser scheme shown in Fig. 1a. Two single-mode Fabry-Pérot cavities 1 and 2 are resonantly coupled by sharing one partial reflecting mirror, through which the photons inside one cavity enter into the other cavity at an effective rate g C . Two cavities have the same resonant frequency but their photon loss rates, denoted, respectively, by κ 1 and κ 2 , may be different. An ensemble of active particles is located inside the active cavity 1, while the passive cavity 2 is empty. Each particle behaves independently of the others and is modeled by a two-state system that is composed of the upper ei j and lower gi states. The lasing ei j À gi transition (frequency ω 0 ) is resonantly coupled to the cavity 1 with a strength g A . The unavoidable ARTICLE COMMUNICATIONS PHYSICS | https://doi.org/10.1038/s42005-021-00575-7 dissipative processes (such as spontaneous emission) reduce the number of particles in ei j and jgi at the rates γ 0 e and γ g , respectively. The decay rate of the particles from ei j to gi is γ e (smaller than γ g ) and the particle polarization degrades at a rate γ eg ¼ ðγ 0 e þ γ g Þ=2. In what follows, we assume that the strong light-particle coupling regime with 4g 2 A =γ 2 eg ) 1 and 4g 2 A =κ 1 γ eg ) 1 is accessed. The particles are pumped into ei j at a rate R, producing the population inversion. The lasing action occurs once R exceeds the total loss of the open system. We emphasize that the diagram shown in Fig. 1a just presents one feasible coupling scheme. There are also other ways to achieve the intercavity coupling, for example, the spatial/spectral overlap between two whispering gallery modes 5,27 (WGMs), the photonic crystal (PhC) microcavities in one slab with asymmetric optical gains 12 , and the electric/magnetic coupling between superconducting circuits 28 . The quality factor of a glass WGM microsphere 29 can reach as high as 8 10 9 , corresponding a photon loss rate of a few 10 5 s −1 for a light wavelength of 800 nm. While the PhC microcavities 30 have the relatively low-quality factors (up to 10 7 ), the intercavity coupling strength may greatly exceed the cavity loss rates, g C ) κ k¼1;2 . By contrast, the superconducting cavities, operating at a frequency of a few GHz, have the advantages of flexibility, tunability, and scalability, but their quality factors are much lower (~10 3 ). The optical gain materials inside cavity 1 can be four-level atomic gases 21,24 , fluorescent dye molecules 31 , ion-doped crystals 32,33 , two-level atomic beams 22,25 , and superconducting qubits 34 . Recent experiments [35][36][37] have manifested the strong (artificial) particle-cavity coupling.
The open system described above may be investigated by using the Heisenberg-Langevin approach 34 (Supplementary Note 1). The equations of motion of the modified cavity fields A k¼1;2 t ð Þ ¼ ffiffiffiffiffiffiffiffiffiffi n k t ð Þ p e iθ k t ð Þ still take the form of Eq. (1) but with a different expression of the Hamiltonian Here, θ k¼1;2 t ð Þ denote the cavity-field phases relative to that of the macroscopic medium polarization. The diagonal elements G t ð Þ and L t ð Þ depend on the macroscopic variables (i.e., polarization and populations in two laser states) of the active particles and represent the net optical gain (in cavity 1) and loss (in cavity 2), respectively (see Supplementary Equations S8a and S8b). The non-diagonal elements account for the energy exchange between two cavity fields through the shared mirror. In general, G t ð Þ and L t ð Þ are complex and the common gain-loss balance G t ð Þ ¼ L * t ð Þ is not maintained at any time. One may rewrite Eq. (3) as i d dt Here, we have defined the diagonal element μ t ð Þ can be interpreted as an energy-shifted Hamiltonian.
In the laser dynamics, the pump process continuously injects the energy into the open system, while the unavoidable decay sources, such as the spontaneous emission of the active particles and the photons escaping from the optical cavities, dissipate the system's energy into the environment. The balance of these two processes brings the open system to a steady state (denoted by the subscript ss) when the time scale of interest is much longer than the lifetimes of particle states and intracavity photons. Consequently, Eq. (4) is reduced to a time-independent Schrödingerlike equation with μ ss ¼ G ss À L ss ð Þ =2 and λ ss ¼ À G ss þ L ss ð Þ =2 (see Supplementary Equations S12a and S12b  Fig. 1 Two coupled optical cavities in the lasing regime. a Conventional laser system. Two Fabry-Pérot cavities 1 and 2 share the same mirror (purple disc), resulting in an intercavity coupling strength g C . The photon loss rates of two cavities are κ 1 and κ 2 , respectively. The active particles are located inside cavity 1, which is coupled to the laser transition of the particles with a strength g A . The laser transition is formed by a pair of upper ei j and lower gi j levels of particle. The particles are pumped onto ei j at a rate R. The ei j -populated particles decay at a rate γ 0 e , where the decay rate corresponding to the ei j À gi j branch is γ e . The decay rate of the gi j -populated particles is γ g . The relaxation constant of the particle polarization is γ eg ¼ ðγ 0 e þ γ g Þ=2. In this study, we set κ 2 ¼ 2κ 1 , γ e ¼ 10κ 1 , and γ 0 e ¼ γ g ¼ g A ¼ 20κ 1 and the strong-coupling regime with g A >ðκ 1 ; γ eg Þ is accessed. b One-particle laser. An active particle is placed inside cavity 1. The particle is characterized by four states, where a pump light drives the 1i j À 4i j transition with a strength Ω and the coupling strength between the lasing 2i j À 3i j transition and cavity 1 is g A . The population inversion is achieved via the spontaneous emission from 4i j to 3i j . The decay rate of the particle from ii j to ji j with i; j ¼ 1; 2; 3; 4 and i>j is defined as γ ij . The cross-coincidence of the photon emissions from the two cavities is measured by two photodetectors (PDs) that are, respectively, placed in front of two output cavity mirrors. We set system is in the unbroken (broken) phase when the eigenvalue λ ss is real (complex). In addition, a steady-state solution of the open system may be not stable, that is, an arbitrary small perturbation can take the system away from the steady state. Thus, the linear analysis needs to be performed to examine the stability of a steady-state solution (Supplementary Note 1). Figure 2a illustrates the stable steady-state solutions of the intracavity photon numbers n ss;k¼1;2 . Two distinct lasing phases are presented, corresponding to two sets of A ss;k¼1;2 , respectively (see Supplementary Equations S9-S11). In phase 1, n ss;1 is higher than n ss;2 since the photons inside the cavity 2 are entirely from the cavity 1. Raising the pumping rate R enhances both n ss;k¼1;2 . As g C is increased, more photons inside cavity 1 can enter cavity 2 through the shared mirror before leaking out of cavity 1. Hence, n ss;1 falls while n ss;2 rises. Surprisingly, n ss;2 slightly exceeds n ss;1 when g C approaches the phase boundary that is determined by the stability of the steady-state solutions. In addition, varying g C does not change the cavity-field phase θ ss;1 and also two cavity fields maintain a constant phase difference ðθ ss;1 À θ ss;2 Þ for an arbitrary g C (Fig. 2b). This may be interpreted as follows: The intercavity coupling modulates the laser dynamics in cavity 1. For a weak g C , the relatively fast cavity losses κ 1;2 and mediumpolarization relaxation γ eg erase the modulation effects on two cavity-field phases θ ss;k¼1;2 within a time scale much shorter than g À1 C . Thus, θ ss;k¼1;2 are completely determined by the light source, that is, the medium polarization, except the extra constant phase offsets. Actually, the fixed value of θ ss;1 maximizes the interaction of the active particles with cavity 1. Moreover, it is found that in  Steady-state solutions of the conventional laser system. a Dependence of the intracavity photon numbers n ss;k¼1;2 on the pumping rate R and the intercavity coupling strength g C . Here, the subscript ss represents the steady-state solutions and k ¼ 1; 2 denote the cavity 1 and 2, respectively. The photon loss rate κ 1 of the cavity 1 is chosen as the frequency unit in the plot. Four distinct regimes are presented. In the phase 1, n ss;1 decreases, while n ss;2 increases as g C is enhanced. In the phase 2, n ss;k¼1;2 are nearly equal and slowly degrade when g C grows. In the unstable regime, both n k¼1;2 vary periodically in time. In the rest regime, both n ss;k¼1;2 stay at zero, that is, no lasing action occurs. b Steady-state cavity-field phases θ ss;1 and ðθ ss;1 À θ ss;2 Þ as a function of g C . The pumping rate is fixed at R=κ 1 ¼ 300. In the phase 1, both θ ss;1 and ðθ ss;1 À θ ss;2 Þ maintain constant. In the phase 2, the phases vary with g C . c Emission spectra of two cavity fields S 1 ω ð Þ (red) and S 2 ω ð Þ (blue) at several selected values of g C . Here, ω denotes the angular frequency in the frequency domain via the Fourier transform and ω 0 corresponds to the central frequency of the laser transition. The pumping rate is set at R=κ 1 ¼ 300. The symbols are obtained from the numerical simulation, while the solid and dashed curves correspond to the curve fitting. The full width at half-maximum of each spectral curve has been inserted. the lasing phase 1 both μ ss and λ ss are purely imaginary and one has Im μ ss À Á > 0 and Im λ ss ð Þ< 0. Thus, the shifted Hamiltonian H ss has PT symmetry and the coupled system is in the broken phase. Figure 2c displays the dependence of the emission spectra S k¼1;2 ω ð Þ of two cavity fields on g C (see "Methods"). It is seen that only one spectral peak is presented for either cavity field. Interestingly, the spectral linewidth of cavity field 2 is always narrower than that of cavity field 1, despite the fact that cavity field 2 completely comes from cavity field 1. Increasing g C broadens both spectra. In particular, the spectral peak of S 1 ω ð Þ presents a nearly flat top when the coupled system operates in the vicinity of the phase boundary, denoting the onset of the spectral splitting.
In the lasing phase 2, both n ss;k¼1;2 degrade slowly as g C is increased. The analytical expressions of the steady-state solutions (see Supplementary Equations S9-S11) illustrate that for a large g C the effect of the intercavity coupling on A ss;1 is so strong that θ ss;1 can no longer sustain the maximum light-particle interaction in cavity 1 but has to be shifted, thereby attenuating the lasing action. Also, ðθ ss;1 À θ ss;2 Þ becomes g C -dependent (see Fig. 2b) and this weakens the intercavity coupling. The boundary between phases 1 and 2 depends solely upon g C (see Fig. 2a) and the relevant phase transition is found to be of the second order. Within the regime of phase 2, the parameter μ ss in Eq. (5) is complex but not purely imaginary. One may interpret this as follows: The shifted θ ss;1 (compared to the value in the lasing phase 1) effectively introduces an extra frequency shift to cavity field 1 (see Supplementary Equations S12a and S12b), thereby lifting the degeneracy of two cavity modes and resulting in < μ ss À Á ≠0. Thus, the Hamiltonian H ss does not possess the PT symmetry, PTH ss ≠H ss PT, in the lasing phase 2. Varying g C from phase 1 to phase 2 breaks the PT symmetry of H ss . This is unlike the common PT-symmetry breaking, where the intercavity coupling affects the reality of the energy-level spectrum of a PT-symmetric Hamiltonian. By using Supplementary Equations S11d-S11e, the ratio of n ss;1 to n ss;2 is derived as n ss;1 =n ss;2 ¼ ð2γ eg À κ 2 Þ=ð2γ eg þ κ 1 Þ in the lasing phase 2. The two cavities have different photon numbers, n ss;1 <n ss;2 , for a finite γ eg . In the limit γ eg ) κ 1;2 , where the approximation of adiabatically eliminating the medium polarization becomes valid, the relative difference jn ss;1 Àn ss;2 j n ss;k¼1;2 is significantly suppressed and one approximately obtains n ss;1 ¼ n ss;2 . This result corresponds to the unbroken phase of a conventional PT-symmetric coupled-cavity system. Thus, the conventional unbroken PT-symmetric phase is actually an idealization in the adiabatic limit. As depicted in Fig. 2c, two cavity fields have the similar spectral profiles, The strong intercavity coupling gives rise to a spectral splitting whose separation is determined by 2g C .
When the intercavity coupling strength g C is further enhanced, the coupled system is no longer allowed to stay in a steady state in a stable manner. This is because none of the steady-state θ ss;k¼1;2 enable the laser dynamics in the cavity 1 to follow the modulation of the intercavity coupling. As shown in Fig. 3a, both n k¼1;2 t ð Þ oscillate at a frequency about g C and the two oscillations are nearly out of phase, indicating the two cavities exchange the energy periodically. The total photon number n 1 t ð Þ þ n 2 t ð Þ ð Þ also shows an oscillatory behavior. To understand the mechanics behind this, we define No lasing a c d Fig. 3 Unstable lasing dynamics. a Time-dependent intensities of two cavity fields n k¼1;2 t ð Þ with the intercavity coupling strength g C =κ 1 ¼ 99. Here, k ¼ 1; 2 denote the cavities 1 and 2, respectively, and the photon loss rate κ 1 of cavity 1 is chosen as the frequency unit in the plot. Both n k¼1;2 t ð Þ vary periodically with a period of about g À1 C . b Time evolution of the diagonal element μ t ð Þ of the Hamiltonian H t ð Þ with g C =κ 1 ¼ 99. The real part of μ t ð Þ maintains at zero, while its imaginary part varies periodically and diverges within a short duration around n 1 t ð Þ $ 0. c Numerical simulation of the emission spectrum S 1 ω ð Þ of cavity field 1 (red filled circles) with g C =κ 1 ¼ 99. Here, ω denotes the angular frequency in the frequency domain and ω 0 is the central frequency of the laser transition. Only the spectral section on the blue-detuned side is presented and its mirror image gives the spectral section on the reddetuned side. The solid line corresponds to the curve fitting. The spectrum S 2 ω ð Þ of cavity field 2 shows a nearly identical line shape. d Time-averaged intracavity intensities n k¼1;2 ¼ 1 ð Þdt with an arbitrary initial time t 0 and a time duration Δt $ g À1 C vs. the intercavity coupling strength g C . In the stable steady-state (denoted by the subscript ss) regime, we have n k¼1;2 are equal to the steady-state values n ss;k¼1;2 of two cavity fields. For all plots, the pumping rate R is set at R=κ 1 ¼ 300.
Figure 3b displays an example of the diagonal element μ t ð Þ of H t ð Þ. It is seen that μ t ð Þ is purely imaginary, meaning H t ð Þ possesses the PT symmetry, PTH t ð Þ ¼ H t ð ÞPT. In addition, although Im μ t ð Þ ½ exhibits a dispersion behavior around n 1 t ð Þ $ 0, μ t ð Þ j j is smaller than the off-diagonal elements g C of H t ð Þ. Thus, the two eigenvalues of H t ð Þ are both real at any time, that is, the coupled system is in the PT-symmetry unbroken phase. The transition from the lasing phase 2 to the unstable regime hardly affects the spectral profiles of two cavity fields (see Fig. 3c). Figure 3d illustrates the dependence of the time-averaged intensities n k¼1;2 ¼ 1 C on g C , where n k¼1;2 first go up and then fall to zero as g C is increased. The boundary of the unstable regime on the high-g C side is determined by the threshold of the steady-state laser solutions, beyond which the lasing action ceases because the strong mode splitting caused by the intercavity interaction pushes the modes out of the gain spectrum of the medium polarization.
The phase diagram of the macroscopic coupled system in the R À g C plane is summarized in Fig. 4, where the phase boundaries are determined by the lasing thresholds (TH1 and TH2) and the linear stability analysis (ST1 and ST2) of the steady-state solutions A ss;k¼1;2 . The lasing phase 1 is surrounded by a section of the TH1 line and the ST1 line. In this phase, θ ss;k¼1;2 are g C -independent and maximize both the light-particle interaction in cavity 1 and the intercavity coupling. The steady-state energyshifted Hamiltonian H ss is PT symmetric but has purely imaginary eigenvalues. Thus, the coupled system is in the PTsymmetry broken phase. The border of the lasing phase 2 consists of the ST1 and ST2 lines and a section of the TH2 line. In this phase, H ss does not satisfy the PT symmetry, PTH ss ≠H ss PT. The coupled system possesses an unstable regime surrounded by the ST2 line and a section of the TH2 line. This regime actually corresponds to the PT-symmetry unbroken phase, where H t ð Þ is PT symmetric and both fields exhibit an oscillatory behavior with a frequency of $ g C . The rest regime of the R À g C plane accounts for the no lasing phase, where the population inversion of the active particles cannot be achieved. The coupled system in different phases may be identified from the distinct behaviors of the photon numbers n k¼1;2 . In addition, two cavity outputs enable a direct measurement of the emission spectra from their beat note.
Microlaser. As the system's size (measured by the number of particles and photons) is significantly reduced, the quantum effects (such as light's particle nature and cavity quantum electrodynamics) become predominant. Our next study is dedicated to the PT symmetry of a microlaser (cavity 1) resonantly interacting with an empty cavity (cavity 2) as illustrated in Fig. 1b. The photon creation and annihilation operators for the cavity k ¼ 1; 2 are a y k and a k , respectively. The microlaser consists of one active particle placed between two highly reflective mirrors. The energylevel structure of the particle involves four states 1; 2; 3; 4i j , whose energies are in an ascending order. The decay rate from ii j (higher) to ji (lower) is γ ij . An external laser beam resonantly drives the 1i j À 4i j transition with a Rabi frequency Ω and cavity 1 is resonantly coupled to the lasing 2i j À 3i j transition (frequency ω 0 ) with a strength g A . Reducing the mode volume of cavity 1 may strongly enhance g A so that only a few or even less than one photon can saturate the lasing transition, that is, Again, the photon loss rates of two cavities are κ 1 and κ 2 , respectively, and g C denotes the intercavity coupling strength. In what follows, we restrict ourselves to the system operating in the strong photon-particle coupling limit, that is, the lasing regime with 4g 2 A =Γ 2 32 ) 1 and 4g 2 A =κ 1 Γ 32 ) 1. We emphasize that such a oneparticle laser is experimentally feasible by means of current optical techniques, where the active particle can be an ion 26 or a quantum dot 38 . The dynamics of the coupled system is governed by the Lindblad master equation 39 , which combines both the coherent particle-cavity and intercavity interactions and the dissipative processes of the particle and intracavity photons (Supplementary Note 2). We employ the Monte Carlo wavefunction method 40 to investigate the PT-symmetric properties of the system's unique steady state 41 .
As depicted in Fig. 5a, the dependence of the steady-state photon numbers n ss;k¼1;2 ¼ ha y k t ð Þa k t ð Þi t!1 (here a y k t ð Þ and a k t ð Þ are the photon operators in the Heisenberg picture and h i t!1 denotes the expectation value in the steady state) on the intercavity coupling strength g C can be divided into four regimes, whose boundaries are found by the slope sign change of the n ss;k¼1;2 vs. g C curves. In the weak-g C regime, n ss;1 is higher than n ss;2 and their difference Δn ss ¼ n ss;1 À n ss;2 decreases towards zero as g C is enhanced. Meanwhile, the photon number distributions of the two cavity fields tend to be the same (Fig. 5b). In the strong-g C regime, the two cavity fields have similar photon numbers that roughly maintain a constant value as g C grows. The photon distributions of the two cavity fields are almost the same and they apparently differ from the Poisson distribution. For a  Fig. 4 Phase diagram of the macroscopic-sized open system. The pumping rate R and intercavity coupling strength g C are tuned widely. The loss rate κ 1 of cavity 1 is chosen as the frequency unit in the plot. The whole R À g C plane is divided into four regimes by four curves. Two threshold lines (TH1, dotted and TH2, dash-dotted) are determined by two sets of steady-state solutions of the coupled system. Two stability boundary lines (ST1, solid and ST2, dashed) are obtained from the linear stability analysis of the steady-state solutions. In the lasing phase 1, the steady-state Hamiltonian H ss possesses the PT symmetry but the coupled system is in the broken phase. In the lasing phase 2, H ss is non-PT symmetric. In the unstable regime, H t ð Þ becomes PT symmetric and the coupled system is in the unbroken phase. In the rest regime, no lasing action occurs.
large enough g C , both n ss;k¼1;2 degrade since the cavity-mode shifts induced by the intercavity interaction weaken the lasing action in cavity 1. Eventually, n ss;k¼1;2 are suppressed to a nearzero value. Unlike the macroscopic counterpart, the microscopic system does not exhibit an unstable regime.
To map the steady-state solutions of the coupled system onto the time-independent Schrödinger-like Eq. (5), we write the cavity fields as A ss;k¼1;2 ¼ ffiffiffiffiffiffiffiffi n ss;k p e iθ ss;k with the field phases θ ss;k relative to that of the particle's polarization (Supplementary Note 2). The diagonal element and eigenvalue of the shifted Hamiltonian H ss are then given by μ ss ¼ ðg C =2ÞðA ss;1 =A ss;2 À A ss;2 =A ss;1 Þ and λ ss ¼ ðg C =2ÞðA ss;1 =A ss;2 þ A ss;2 =A ss;1 Þ, respectively. It should be noted that the expressions of μ ss and λ ss in terms of A ss;k¼1;2 are also valid for the macroscopic coupled system. The numerical simulation results manifest θ ss;1 ¼ π=2 and θ ss;2 ¼ 0 for an arbitrary g C . Consequently, μ ss is purely imaginary and H ss owns the PT symmetry, PTH ss ¼ H ss PT. In the weak-g C regime, λ ss is purely imaginary since n ss;1 differs from n ss;2 , that is, Δn ss ≠0. Hence, the coupled system is in the PT-symmetry broken phase (Fig. 5a). By contrast, n ss;1 is almost equal to n ss;2 for a large enough g C , that is, Δn ss $ 0, indicating that λ ss % 0 and the coupled system is in the unbroken phase.
The PT-symmetry breaking is accompanied by a change in the two cavity-field spectra S k¼1;2 ω ð Þ (see "Methods"). When the coupled system is in the broken phase, for a small g C both S k¼1;2 ω ð Þ present only one spectral peak located at the resonance ω ¼ ω 0 . As g C grows, the spectral peak in S 1 ω ð Þ is split into two due to the intercavity interaction while S 2 ω ð Þ is still single-peaked (Fig. 6a). The spectrum S 2 ω ð Þ starts to split when the coupled system approaches the phase transition boundary. Unlike the double-peaked S 1 ω ð Þ, the two spectral peaks in S 2 ω ð Þ are unresolvable in principle since the interpeak separation is smaller than the peak linewidth (Fig. 6b). Interestingly, in the broken phase the peak linewidth of S 1 ω ð Þ always exceeds that of S 2 ω ð Þ. In other words, the cavity field 2 has a temporal coherence superior to the cavity field 1, although the former entirely comes from the latter. As g C is further increased, the coupled system enters the unbroken phase, where both S k¼1;2 ω ð Þ exhibit a double-peaked profile. The two spectra have similar interpeak separations ($ 2g C ) and similar peak linewidths. For either S k¼1;2 ω ð Þ, the separation between two spectral peaks surpasses the peak a Dependence of n ss;k¼1;2 and the photon difference Δn ss ¼ n ss;1 À n ss;2 on the intercavity coupling strength g C . Here, the subscript ss denotes the steady-state solutions and k ¼ 1; 2 represent the cavities 1 and 2, respectively. The loss rate κ 1 of cavity 1 is chosen as the frequency unit in the plot. Four regimes are presented according to the slope sign change of the n ss;k¼1;2 vs. g C curves. In the first regime, n ss;k¼1;2 are apparently different and the coupled system is in the PT-symmetry broken phase. In the rest three regimes, Δn ss approximates zero and the coupled system is in the PT-symmetry unbroken phase. b Normalized photon distributions of two cavity fields for different g C (histograms). The Poisson distributions based on the corresponding n ss;1 (circles) and n ss;2 (squares) are shown for comparison.
linewidth, leading to two resolvable peaks. Therefore, the PT-symmetry breaking can be detected through the measurement of the lasing spectra. Moreover, we find that the crossspectrum S 12 ω ð Þ between two cavity fields (see "Methods") is almost the same as S 2 ω ð Þ as shown in Fig. 6a. That is, the cross-correlation between the two cavity fields is primarily determined by cavity field 2, whose coherence is greater than cavity field 1.
We are also interested in the second-order (intensity) autocorrela- ss;k with a time difference τ. Such detection can be performed by the well-known Hanbury Brown-Twiss setup 42 . The statistical characteristics of the photon emission from individual cavities are closely related to C 2 ð Þ k¼1;2 τ ð Þ in the short-time-delay limit τ ¼ 0. As illustrated in Fig. 7a, both C 2 ð Þ k¼1;2 0 ð Þ are equal or higher than unity. Thus, the photon distributions of two cavity fields exhibit a Poisson or a super-Poissonian statistics (Fig. 5b), corresponding to the photon bunching effect. This is ascribed to the fact that the open system operates in the good-cavity limit (i.e., κ k¼1;2 ( γ 32 ), where the ability of the cavities to accumulate the photons plays a notable role. One may expect a nonclassical behavior (i.e., photon antibunching/single-photon emission) of individual cavity fields when the cavity loss rates approximate to (see below) or much exceed the decay rates of the particle states (i.e., the bad-cavity limit, Supplementary Note 2). In this circumstance, one photon that is generated via the stimulated emission of the active particle may escape the cavities before the particle population is inverted again for the next stimulated emission event. The resulting photon distributions of individual cavity fields are characterized by a sub-Poissonian statistics. In addition, the single-photon emission can also be performed through the pulsed pump process 43 (Supplementary Note 2). When the coupled system is in the broken phase, C 2 ð Þ k¼1;2 τ ð Þ gradually go up towards unity as τ is increased. By contrast, for the coupled system in the unbroken phase, the dependence of C 2 ð Þ k¼1;2 τ ð Þ on τ displays a damped oscillatory behavior with an oscillation frequency $ 2g C (Fig. 7b).
In addition, one may examine the intensity cross-correlation functions C 2 ð Þ 12 τ ð Þ ¼ a y 2 t ð Þa y 1 t þ τ ð Þa 1 t þ τ ð Þa 2 t ð Þ t!1 =n ss;1 n ss;2 and C 2 ð Þ 21 τ ð Þ ¼ a y 1 t ð Þa y 2 t þ τ ð Þa 2 t þ τ ð Þa 1 t ð Þ t!1 =n ss;1 n ss;2 , which are measured via the coincidence counting interferometry as shown in Fig. 1b. We have C ð Þ of the coupled system go down below unity as the coupling strength g C is increased. That is, two cavities are unlikely to emit the photons at the same time, even though each cavity emits the photons in bunches (C 2 ð Þ k¼1;2 0 ð Þ ≥ 1). This is because, for instance, after cavity 2 emits a photon, a portion of the photons inside cavity 1 is consumed to compensate the loss of the cavity field 2, causing cavity 1 hardly emits a photon. A large g C facilitates such a process. Interestingly, the anticorrelation between two cavity fields happens before the occurrence of the PT-symmetry breaking. Thus, the PT symmetry does not account for all the properties of the coupled system. Both C 2 ð Þ k¼12;21 τ ð Þ undergo the damped oscillations at a frequency $ 2g C and approach unity as τ goes to infinity, but their behaviors are still slightly different because of the asymmetry of the physical system (Fig. 7b). Using the bad cavities or the pulsed pump approach,  Fig. 6 Emission spectra of two cavity fields. a Spectra S k¼1;2 ω ð Þ and cross-spectrum S 12 ω ð Þ at several different intercavity coupling strengths g C . Here, k ¼ 1; 2, respectively, represent the cavities 1 and 2, ω denotes the angular frequency in the frequency domain, and ω 0 corresponds to the central frequency of the laser transition. The loss rate κ 1 of cavity 1 is chosen as the frequency unit in the plot. b Interpeak separations of S k¼1;2 ω ð Þ and S 12 ω ð Þ as a function of g C . Inset: Linewidths of the spectral peaks in S k¼1;2 ω ð Þ and S 12 ω ð Þ vs. g C . The peak linewidths are evaluated by means of curve fitting (see "Methods"). The arrow shows the phase transition boundary, above which the interpeak separations exceed the linewidths of spectral peaks in both spectra S k¼1;2 ω ð Þ.
one may obtain both strong photon antibunching of individual cavity fields C  Fig. 7a, χ is always lower than unity in the good-cavity limit, that is, the absence of the nonclassical feature. Increasing the cavity loss rates κ k¼1;2 may suppress both C 2 ð Þ k¼1;2 0 ð Þ, giving rise to the potential of violating the inequality. Although large κ 1;2 reduce the cavity-field intensities n ss;k¼1;2 (Fig. 7c), the coupled system still operates in the lasing regime when the strong particle-cavity coupling limit (i.e., 4g 2 A =Γ 2 32 ) 1 and 4g 2 A =κ 1 Γ 32 ) 1) is reached. Figure 7c illustrates that in the weak-g C regime, the ratio C ð Þ, resulting in χ > 1, that is, the nonclassical photon-pair generation. Meanwhile, for an appropriate g C both C 2 ð Þ k¼1;2 0 ð Þ may be below unity, that is, the nonclassical photon emission of individual cavity fields. By contrast, when g C >κ 1 , C 2 ð Þ 12 0 ð Þ is lower than both C 2 ð Þ k¼1;2 0 ð Þ and the inequality χ ≤ 1 becomes satisfied again. In summary, our work explores fundamental properties of PTsymmetric optics in the lasing regime, where, to the best of our Here, γ ij with i; j ¼ 1; 2; 3; 4 and i>j corresponds to the decay rate of the particle from the ii j to ji j state, Ω is the Rabi frequency of the pump light coupling to the particle's 1i j À 4i j transition, and g A is the coupling strength between the cavity 1 and the laser 2i j À 3i j transition. The coupled system still operates in the strong photon-particle coupling limit, that is, 4g 2 A =Γ 2 32 ) 1 and 4g 2 The gray shade denotes the regime of the nonclassical photon-pair generation with χ>1, that is, the violation of the Cauchy-Schwarz inequality.
knowledge, the microscopic coupled structures have not been accessed before. We extend the PT-symmetric optics from the weak to the strong particle-cavity coupling limit, where the particle polarization cannot be adiabatically eliminated. The resulting phase diagrams differ from the usual PT-symmetric optical schemes based on the rate-equation model [16][17][18][19][20] . The phase boundaries of the macroscopic-sized system are well defined by the lasing thresholds and the stability analysis of the steady-state solutions. The emission spectra can be employed to detect the PT-symmetric properties, avoiding the extra noise sources caused by the technical issues in the transmission spectrum measurement. By contrast, the interphase boundary becomes fuzzy as the system's size is significantly reduced to the microscopic scale. The intensity auto-and cross-correlation functions of the microscopic-sized system with appropriate cavity loss rates illustrate the nonclassical photon emission of individual cavities and the nonclassical photon-pair generation.

Conclusion
Recent developments in atomic physics, microcavity technology, and nanophotonics will allow the predictions obtained in this paper to be tested. The macroscopic lasers have been demonstrated by high-Q optical cavities filled with dilute atomic gases [21][22][23][24][25] and fluorescent dye molecules 31 , ion-doped PhC cavities 32,33 , and dye-doped polymer resonators 46 . The intercavity coupling can be either in a direct manner (i.e., the spatial overlap of two cavity modes) or mediated via another optical device. The experimentally achieved g C , for example, between two PhC cavities exceeds 10 2 GHz 12 . In contrast, the microlasers can be implemented via one quantum dot embedded in a nanoscale PhC/semiconductor cavity 38,47 . The mode volume of a PhC nanocavity is extremely small,~10 −2 μm 3 , resulting in a large g A of~10 GHz 38 . Also, the microlaser structure may consist of one ion/molecule and a WGM microcavity 35 , where the localized surface plasmon resonance of metal nanoparticles is utilized as a solution to achieve the strong particle-cavity coupling 48,49 . Moreover, our study is potentially extended to coupled systems involving multi-energy-level particles 50 and particles with longrange interactions 51 .
The PT-symmetric lasing scheme may be applied for sensing and gyroscope applications [6][7][8]52 . The spectral sensitivity of a modified non-Hermitian optical system to a small change Δg C of the intercavity coupling g C scales as ffiffiffiffiffiffiffiffi Δg C p around the EP 6 , exceeding the sensitivity around the normal diabolic point. For our macroscopic coupled system, we find that the changes of the photon numbers n ss;k¼1;2 are proportional to Δg C around the phase 1-phase 2 transition point. However, due to the simulation restriction, we are unable to confirm if the dependence of the spectral splitting (see Fig. 2c) on Δg C follows the square root law around the phase transition point. In practice, our laser system facilitates this spectrum measurement directly from the beat note between two cavity outputs. In comparison, for the microscopicsized system the small numbers n ss;k¼1;2 challenge the detection of Δg C . Nevertheless, since the particle-cavity interaction depends strongly on the square root of the number N of the particles located inside the cavity 53 , this can be utilized to detect the particle number with a sensitivity that scales as 1= ffiffiffiffi N p and is maximized when N ¼ 1, that is, the single-particle level.

Methods
Lasing spectra of macroscopic-sized system. The cavity-field spectra of the classical optical system may be obtained by numerically solving the c-number Langevin equations listed in Supplementary Note 1. For simplicity, we assume the open system operates at zero temperature and the Langevin forces related to the cavity fields are eliminated. The laser spectral linewidth is primarily caused by the fluctuations in the macroscopic polarization M t ð Þ of the active particles. One may only take into account its associated Langevin force F M t ð Þ and omit the Langevin forces exerting on other particle variables. The Langevin force F M t ð Þ is characterized by the correlation functions hF M t ð Þi ¼ 0, hF M t ð ÞF M t 0 ð Þi ¼ 2D M; M ð Þ δ t À t 0 ð Þ, and hF * M ðtÞF M ðt 0 Þi ¼ 2DðM * ; MÞδðt À t 0 Þ. with an extra phase difference φ of cavity field 2 relative to cavity field 1. The corresponding cross-spectrum reads S 12 ω ð Þ / R Re½C 1 ð Þ 12 τ ð Þe Àiωτ dτ. Here we set φ ¼ π=2 due to the steady-state phase difference Δθ ss;1 À Δθ ss;2 ¼ π=2 between two cavity fields. To evaluate the linewidths and frequency shifts of the spectra S k¼1;2;12 ω ð Þ, one may fit the exponential decay function Ae Àγτ=2 cos Δτ þ θ ð Þwith an amplitude A, a decay constant γ, an oscillation frequency Δ and a phase bias θ to the correlation functions C 1 ð Þ k¼1;2;12 τ ð Þ. The linewidth of the peak(s) in a spectrum is determined by γ. Due to the nonzero θ, the separation between two peaks in a double-peaked spectrum differs from 2Δ, especially for the coupled system in the broken phase. Thus, the interpeak separation is computed directly from the double-peaked spectrum.

Data availability
All data supporting the findings of this study are available from the corresponding authors upon reasonable request.

Code availability
The computer code to simulate the dynamics is available from the corresponding authors upon reasonable request.