Exceptional point in a trimer chain of oscillators with a quadratic driving

M Shoufie Ukhtary Research Center for Quantum Physics, National Research and Innovation Agency (BRIN), South Tangerang 15314, Indonesia [email protected]    Albert Andersen Research Center for Quantum Physics, National Research and Innovation Agency (BRIN), South Tangerang 15314, Indonesia Department of Physics, Faculty of Mathematics and Natural Sciences, Universitas Indonesia, Depok 16424, Indonesia    Donny Dwiputra    M. Jauhar Kholili Research Center for Quantum Physics, National Research and Innovation Agency (BRIN), South Tangerang 15314, Indonesia
(January 2, 2025)
Abstract

Exceptional points of a dissipative chain of three coupled oscillators (trimer), which is driven by quadratic photon, are investigated. The exceptional points emerge from the coalescence of both eigenvalues and eigenvectors of the dynamical matrix that describes the first moments of the trimer. At the exceptional point, we found that the optical spectrum is split into two peaks, instead of a conventional single peak, as in the case of a single oscillator. In particular, the positions of these peaks correspond to the natural frequency of the trimer in a closed system, which depends only on the coupling strength. Furthermore, after passing the exceptional point, the peak positions do not change, which can be used to estimate the coupling strength between oscillators.

I Introduction

Exceptional point (EP) has been the key feature of a non-Hermitian system, in which gain and loss play an important role in the realization of EP [1, 2, 3, 4]. At EP, both the eigenvalues and the eigenvectors of an operator describing the dynamics of a system coalesce into a singularity followed by a transition in dynamics [1, 2, 3, 4, 5, 6]. For example, the population dynamics in a quantum system changes from oscillatory in time to exponentially increasing after passing EP [6, 7]. The collapse of both eigenvalues and eigenvectors presents a strong response of the system to external perturbation [1]. In particular, this strong response to perturbation has been applied for developments of sensitive sensor working in the vicinity of EP [1, 8, 9, 10, 11, 12].

The dynamics of an open system is often described by the Lindbladian master equation [13, 14, 15]. The nonunitary effects are incorporated into the master equation by the so-called dissipators associated with their jump operators. From the master equation, the associated dynamical matrix of statistical moments can be derived. This matrix fully describes the evolution of the system over time, since the eigenvalues give the frequency of the system. Moreover, the coalesce of the eigenvalues and their corresponding eigenvectors of this matrix, thus, corresponds to EP of the system [6, 7, 16, 17]. To investigate the presence of EP experimentally, the optical spectrum or the power spectrum is measured [18, 19, 20, 14, 15]. Passing EP usually changes the number of peaks and the lineshape of the spectrum. Normally, the peak-to-peak separation corresponds to the splitting of the eigenvalues. Therefore, approaching EP, the separation becomes smaller since at least two eigenvalues are merging. In a simple system consisting of two objects, the spectrum changes from doublet to singlet as we pass EP [19].

Typically, EP is investigated in an interacting multipartite system, where the coupling constant, gain, or loss rate determine EP. Downing et al. considered a case of a single oscillator driven parametrically, from which they found EP determined by the driving strength and the loss rate [19]. Furthermore, they found that the optical spectrum changes from doublet to singlet when passing EP [19]. In this work, we investigate EP in a multipartite system consisting of three interacting quantum oscillators driven parametrically, as shown in Fig. 1. Our aim is to study the effect of coupling on EP and the corresponding optical spectrum. At EP, we found that the optical spectrum is split into two peaks, instead of conventional single peak as in the case of a single oscillator. In particular, the positions of these peaks correspond to the natural frequency of the trimer in a closed system, which depends only on the coupling strength. Therefore, we can apply our system to estimate the coupling strength using the optical spectrum at EP.

II Model

Refer to caption
Figure 1: A chain of interacting trimer. The system is driven quadratically with strength of Ω/2Ω2\Omega/2roman_Ω / 2 on oscillator B. We allow dissipation with rate of 2γ2𝛾2\gamma2 italic_γ on oscillator B. Oscillator B is coupled with two other oscillators with coupling constant of J𝐽Jitalic_J.

The trimer consists of three interacting oscillators as illustrated in Fig. 1, with level spacing of each oscillator ωi(i=a,b,c)subscript𝜔𝑖𝑖𝑎𝑏𝑐\omega_{i}\,(i=a,b,c)italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_i = italic_a , italic_b , italic_c ). We consider a quadratic driving on oscillator B, which is coupled with two other oscillators with a coupling constant of J𝐽Jitalic_J. In this case, we consider the system that obtains energy from external quadratic driving, but the energy dissipates to the environment incoherently. The system Hamiltonian is expressed as follows:

H𝐻\displaystyle Hitalic_H =ωaaa+ωbbb+ωccc+J(ab+ab+bc+bc)absentsubscript𝜔𝑎superscript𝑎𝑎subscript𝜔𝑏superscript𝑏𝑏subscript𝜔𝑐superscript𝑐𝑐𝐽superscript𝑎𝑏𝑎superscript𝑏superscript𝑏𝑐𝑏superscript𝑐\displaystyle=\omega_{a}a^{\dagger}a+\omega_{b}b^{\dagger}b+\omega_{c}c^{% \dagger}c+J(a^{\dagger}b+ab^{\dagger}+b^{\dagger}c+bc^{\dagger})= italic_ω start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_a + italic_ω start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b + italic_ω start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_c + italic_J ( italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b + italic_a italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT + italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_c + italic_b italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT )
+Ω2(bbei2ωDt+bbei2ωDt).Ω2superscript𝑏superscript𝑏superscript𝑒𝑖2subscript𝜔𝐷𝑡𝑏𝑏superscript𝑒𝑖2subscript𝜔𝐷𝑡\displaystyle+\frac{\Omega}{2}\left(b^{\dagger}b^{\dagger}e^{-i2\omega_{D}t}+% bbe^{i2\omega_{D}t}\right).+ divide start_ARG roman_Ω end_ARG start_ARG 2 end_ARG ( italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_e start_POSTSUPERSCRIPT - italic_i 2 italic_ω start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT italic_t end_POSTSUPERSCRIPT + italic_b italic_b italic_e start_POSTSUPERSCRIPT italic_i 2 italic_ω start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT italic_t end_POSTSUPERSCRIPT ) . (1)

To remove the explicit time dependence of the Hamiltonian, we move to a frame rotating with the driving by using the unitary operator of U=eiωDt(aa+bb+cc)𝑈superscript𝑒𝑖subscript𝜔𝐷𝑡superscript𝑎𝑎superscript𝑏𝑏superscript𝑐𝑐U=e^{-i\omega_{D}t(a^{\dagger}a+b^{\dagger}b+c^{\dagger}c)}italic_U = italic_e start_POSTSUPERSCRIPT - italic_i italic_ω start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT italic_t ( italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_a + italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b + italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_c ) end_POSTSUPERSCRIPT. The transformed Hamiltonian is given as follows:

H𝐻\displaystyle Hitalic_H =Δaaa+Δbbb+Δccc+J(ab+ab+bc+bc)absentsubscriptΔ𝑎superscript𝑎𝑎subscriptΔ𝑏superscript𝑏𝑏subscriptΔ𝑐superscript𝑐𝑐𝐽superscript𝑎𝑏𝑎superscript𝑏superscript𝑏𝑐𝑏superscript𝑐\displaystyle=\Delta_{a}a^{\dagger}a+\Delta_{b}b^{\dagger}b+\Delta_{c}c^{% \dagger}c+J(a^{\dagger}b+ab^{\dagger}+b^{\dagger}c+bc^{\dagger})= roman_Δ start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_a + roman_Δ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b + roman_Δ start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_c + italic_J ( italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b + italic_a italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT + italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_c + italic_b italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT )
+Ω2(bb+bb),Ω2superscript𝑏superscript𝑏𝑏𝑏\displaystyle+\frac{\Omega}{2}\left(b^{\dagger}b^{\dagger}+bb\right),+ divide start_ARG roman_Ω end_ARG start_ARG 2 end_ARG ( italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT + italic_b italic_b ) , (2)

where Δi=ωiωDsubscriptΔ𝑖subscript𝜔𝑖subscript𝜔𝐷\Delta_{i}=\omega_{i}-\omega_{D}roman_Δ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = italic_ω start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_ω start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT is the detuning frequency of oscillator-i𝑖iitalic_i. For simplicity, we assume that oscillators A and C are resonant with driving, while we allow detuning of the frequency of oscillator B as ΔΔbΔsubscriptΔ𝑏\Delta\equiv\Delta_{b}roman_Δ ≡ roman_Δ start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT.

The system is coupled to the environment through dissipation in oscillator B with a dissipation rate of 2γ2𝛾2\gamma2 italic_γ, as illustrated in Fig. 1. This allows us to compare it with the case of a dissipative oscillator driven parametrically studied by Downing et al. and to understand the impact of the coupling. The dynamics of the system is governed by the following master equation for density matrix ρ𝜌\rhoitalic_ρ,

tρsubscript𝑡𝜌\displaystyle\partial_{t}\rho∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_ρ =i[H,ρ]+γ(2bρbbbρρbb),absent𝑖𝐻𝜌𝛾2𝑏𝜌superscript𝑏superscript𝑏𝑏𝜌𝜌superscript𝑏𝑏\displaystyle=-i[H,\rho]+\gamma\left(2b\rho b^{\dagger}-b^{\dagger}b\rho-\rho b% ^{\dagger}b\right),= - italic_i [ italic_H , italic_ρ ] + italic_γ ( 2 italic_b italic_ρ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT - italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b italic_ρ - italic_ρ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b ) , (3)

where the first term describes the Hermitian evolution of the system, while the second term gives the dissipation of the system. Using Eq. (3) and the trace property tO=Tr[Otρ]subscript𝑡delimited-⟨⟩𝑂Trdelimited-[]𝑂subscript𝑡𝜌\partial_{t}\langle O\rangle=\textrm{Tr}[O\partial_{t}\rho]∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ⟨ italic_O ⟩ = Tr [ italic_O ∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_ρ ], we derive the following differential equation for the first moment of the system,

itΨ=Ψ,𝑖subscript𝑡ΨΨ\displaystyle i\partial_{t}\Psi=\mathcal{H}\Psi,italic_i ∂ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT roman_Ψ = caligraphic_H roman_Ψ , (4)

where Ψ=(a,b,c,a,b,c)TΨsuperscriptdelimited-⟨⟩𝑎delimited-⟨⟩𝑏delimited-⟨⟩𝑐delimited-⟨⟩superscript𝑎delimited-⟨⟩superscript𝑏delimited-⟨⟩superscript𝑐𝑇\Psi=(\langle a\rangle,\langle b\rangle,\langle c\rangle,\langle a^{\dagger}% \rangle,\langle b^{\dagger}\rangle,\langle c^{\dagger}\rangle)^{T}roman_Ψ = ( ⟨ italic_a ⟩ , ⟨ italic_b ⟩ , ⟨ italic_c ⟩ , ⟨ italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ⟩ , ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ⟩ , ⟨ italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ⟩ ) start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT and

=(0J0000JΔiγJ0Ω00J00000000J00Ω0JΔiγJ0000J0)matrix0𝐽0000𝐽Δ𝑖𝛾𝐽0Ω00𝐽00000000𝐽00Ω0𝐽Δ𝑖𝛾𝐽0000𝐽0\displaystyle\mathcal{H}=\begin{pmatrix}0&J&0&0&0&0\\ J&\Delta-i\gamma&J&0&\Omega&0\\ 0&J&0&0&0&0\\ 0&0&0&0&-J&0\\ 0&-\Omega&0&-J&-\Delta-i\gamma&-J\\ 0&0&0&0&-J&0\end{pmatrix}caligraphic_H = ( start_ARG start_ROW start_CELL 0 end_CELL start_CELL italic_J end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL italic_J end_CELL start_CELL roman_Δ - italic_i italic_γ end_CELL start_CELL italic_J end_CELL start_CELL 0 end_CELL start_CELL roman_Ω end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL italic_J end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL - italic_J end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL - roman_Ω end_CELL start_CELL 0 end_CELL start_CELL - italic_J end_CELL start_CELL - roman_Δ - italic_i italic_γ end_CELL start_CELL - italic_J end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL - italic_J end_CELL start_CELL 0 end_CELL end_ROW end_ARG ) (5)

The dynamics of the system is determined by the matrix \mathcal{H}caligraphic_H, which can be understood as the effective Hamiltonian of the system. The eigenvalues of \mathcal{H}caligraphic_H, denoted by λ𝜆\lambdaitalic_λ, give the frequencies for the evolution of the system. Due to the presence of γ𝛾\gammaitalic_γ, the eigenvalues are generally complex, which highlights the relaxation of the system to a steady state. We found analytical solutions for the eigenvalues as follows:

λ1,2subscript𝜆12\displaystyle\lambda_{1,2}italic_λ start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT =0,λ3,4=12(i(γ+Π)±ω)formulae-sequenceabsent0subscript𝜆3412plus-or-minus𝑖𝛾Πsubscript𝜔\displaystyle=0,\quad\lambda_{3,4}=\frac{1}{2}\left(-i(\gamma+\Pi)\pm\omega_{-% }\right)= 0 , italic_λ start_POSTSUBSCRIPT 3 , 4 end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( - italic_i ( italic_γ + roman_Π ) ± italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT )
λ5,6subscript𝜆56\displaystyle\lambda_{5,6}italic_λ start_POSTSUBSCRIPT 5 , 6 end_POSTSUBSCRIPT =12(i(γΠ)±ω+),absent12plus-or-minus𝑖𝛾Πsubscript𝜔\displaystyle=\frac{1}{2}\left(-i(\gamma-\Pi)\pm\omega_{+}\right),= divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( - italic_i ( italic_γ - roman_Π ) ± italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ) , (6)

where we define,

ΠΩ2Δ2,ω±ΠsuperscriptΩ2superscriptΔ2subscript𝜔plus-or-minus\displaystyle\Pi\equiv\sqrt{\Omega^{2}-\Delta^{2}},\,\omega_{\pm}roman_Π ≡ square-root start_ARG roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Δ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_ω start_POSTSUBSCRIPT ± end_POSTSUBSCRIPT Z±2γΠ,Zabsentplus-or-minus𝑍2𝛾Π𝑍\displaystyle\equiv\sqrt{Z\pm 2\gamma\Pi},\,Z≡ square-root start_ARG italic_Z ± 2 italic_γ roman_Π end_ARG , italic_Z 8J2γ2Π2.absent8superscript𝐽2superscript𝛾2superscriptΠ2\displaystyle\equiv 8J^{2}-\gamma^{2}-\Pi^{2}.≡ 8 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (7)

From Eq. (6), we found that there are three distinct driving strengths ΩΩ\Omegaroman_Ω that generate exceptional points (EPs), which are given by the zeros of the square roots argument. The EPs are given by the following driving strength,

ΩEP(1)subscriptsuperscriptΩ1EP\displaystyle\Omega^{(1)}_{\textrm{EP}}roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT =Δ,ΩEP±=Δ2+(γ±22J)2formulae-sequenceabsentΔsubscriptsuperscriptΩplus-or-minusEPsuperscriptΔ2superscriptplus-or-minus𝛾22𝐽2\displaystyle=\Delta,\quad\Omega^{\pm}_{\textrm{EP}}=\sqrt{\Delta^{2}+(\gamma% \pm 2\sqrt{2}J)^{2}}= roman_Δ , roman_Ω start_POSTSUPERSCRIPT ± end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = square-root start_ARG roman_Δ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_γ ± 2 square-root start_ARG 2 end_ARG italic_J ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG (8)

We only consider the case of steady state, which is represented by negative imaginary part of the eigenvalue. The positive imaginary part of the eigenvalue leads to spectral collapse in the case of a closed oscillator driven quadratically [19, 7]. The critical driving, above which the steady state is not present, is given by Ωcγ2+Δ2subscriptΩ𝑐superscript𝛾2superscriptΔ2\Omega_{c}\equiv\sqrt{\gamma^{2}+\Delta^{2}}roman_Ω start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT ≡ square-root start_ARG italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_Δ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG. Using this expression, we can limit our parameters to 2J2γ22superscript𝐽2superscript𝛾22J^{2}\leq\gamma^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≤ italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, so that the steady state is present. From Eq. (8), we obtain an EP at the driving strength independent of J𝐽Jitalic_J, ΩEP(1)subscriptsuperscriptΩ1EP\Omega^{(1)}_{\textrm{EP}}roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT. Interestingly, this driving strength also gives EP in the case of a single driven oscillator [19]. Therefore, at any coupling strength, we can always find the EP. In the subsequent sections, we will investigate this EP ΩEP(1)subscriptsuperscriptΩ1EP\Omega^{(1)}_{\textrm{EP}}roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT, which will give us indication whether the driven oscillator is coupled to other oscillators or isolated.

To investigate the properties of the corresponding EP, we derive the first-order correlation function gb(1)(τ)superscriptsubscript𝑔𝑏1𝜏g_{b}^{(1)}(\tau)italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ), from which we can further calculate the optical spectrum of the emitted photon from the driven oscillator B. The normalized correlation function is defined as follows [18, 19, 14, 15]:

gb(1)(τ)=limtb(t)b(t+τ)b(t)b(t),superscriptsubscript𝑔𝑏1𝜏subscript𝑡delimited-⟨⟩superscript𝑏𝑡𝑏𝑡𝜏delimited-⟨⟩superscript𝑏𝑡𝑏𝑡\displaystyle g_{b}^{(1)}(\tau)=\lim_{t\rightarrow\infty}\frac{\langle b^{% \dagger}(t)b(t+\tau)\rangle}{\langle b^{\dagger}(t)b(t)\rangle},italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ) = roman_lim start_POSTSUBSCRIPT italic_t → ∞ end_POSTSUBSCRIPT divide start_ARG ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_b ( italic_t + italic_τ ) ⟩ end_ARG start_ARG ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_b ( italic_t ) ⟩ end_ARG , (9)

where limtb(t)b(t)bbsssubscript𝑡delimited-⟨⟩superscript𝑏𝑡𝑏𝑡subscriptdelimited-⟨⟩superscript𝑏superscript𝑏ss\lim_{t\rightarrow\infty}\langle b^{\dagger}(t)b(t)\rangle\equiv\langle b^{% \dagger}b^{\dagger}\rangle_{\textrm{ss}}roman_lim start_POSTSUBSCRIPT italic_t → ∞ end_POSTSUBSCRIPT ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_b ( italic_t ) ⟩ ≡ ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ⟩ start_POSTSUBSCRIPT ss end_POSTSUBSCRIPT is the steady-state population of oscillator B. The steady-state population is solved from the differential equation for second moments, similar to Eq. (4). The resulted steady-state populations are expressed as follows:

aasssubscriptdelimited-⟨⟩superscript𝑎𝑎ss\displaystyle\langle a^{\dagger}a\rangle_{\textrm{ss}}⟨ italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_a ⟩ start_POSTSUBSCRIPT ss end_POSTSUBSCRIPT =ccss=Ω24(γ2+Δ2Ω2)absentsubscriptdelimited-⟨⟩superscript𝑐𝑐sssuperscriptΩ24superscript𝛾2superscriptΔ2superscriptΩ2\displaystyle=\langle c^{\dagger}c\rangle_{\textrm{ss}}=\frac{\Omega^{2}}{4% \left(\gamma^{2}+\Delta^{2}-\Omega^{2}\right)}= ⟨ italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_c ⟩ start_POSTSUBSCRIPT ss end_POSTSUBSCRIPT = divide start_ARG roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_Δ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG (10)
bbsssubscriptdelimited-⟨⟩superscript𝑏𝑏ss\displaystyle\langle b^{\dagger}b\rangle_{\textrm{ss}}⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b ⟩ start_POSTSUBSCRIPT ss end_POSTSUBSCRIPT =Ω22(γ2+Δ2Ω2)absentsuperscriptΩ22superscript𝛾2superscriptΔ2superscriptΩ2\displaystyle=\frac{\Omega^{2}}{2\left(\gamma^{2}+\Delta^{2}-\Omega^{2}\right)}= divide start_ARG roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 2 ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_Δ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG (11)
bbsssubscriptdelimited-⟨⟩superscript𝑏superscript𝑏ss\displaystyle\langle b^{\dagger}b^{\dagger}\rangle_{\textrm{ss}}⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ⟩ start_POSTSUBSCRIPT ss end_POSTSUBSCRIPT =i(γΩ+iΔΩ)2(γ2+Δ2Ω2).absent𝑖𝛾Ω𝑖ΔΩ2superscript𝛾2superscriptΔ2superscriptΩ2\displaystyle=\frac{i(\gamma\Omega+i\Delta\Omega)}{2\left(\gamma^{2}+\Delta^{2% }-\Omega^{2}\right)}.= divide start_ARG italic_i ( italic_γ roman_Ω + italic_i roman_Δ roman_Ω ) end_ARG start_ARG 2 ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_Δ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG . (12)

gb(1)(τ)superscriptsubscript𝑔𝑏1𝜏g_{b}^{(1)}(\tau)italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ) is derived using the quantum regression theorem [21, 14, 15]. gb(1)(τ)superscriptsubscript𝑔𝑏1𝜏g_{b}^{(1)}(\tau)italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ) is obtained by solving the following differential equation,

τv=iv,subscript𝜏𝑣𝑖𝑣\displaystyle\partial_{\tau}v=-i\mathcal{H}v,∂ start_POSTSUBSCRIPT italic_τ end_POSTSUBSCRIPT italic_v = - italic_i caligraphic_H italic_v , (13)

where \mathcal{H}caligraphic_H is given by Eq. (5) and v=(b(t)a(t+τ),b(t)b(t+τ),b(t)c(t+τ),b(t)a(t+τ),b(t)b(t+τ),b(t)c(t+τ))T𝑣superscriptdelimited-⟨⟩superscript𝑏𝑡𝑎𝑡𝜏delimited-⟨⟩superscript𝑏𝑡𝑏𝑡𝜏delimited-⟨⟩superscript𝑏𝑡𝑐𝑡𝜏delimited-⟨⟩superscript𝑏𝑡superscript𝑎𝑡𝜏delimited-⟨⟩superscript𝑏𝑡superscript𝑏𝑡𝜏delimited-⟨⟩superscript𝑏𝑡superscript𝑐𝑡𝜏𝑇v=(\langle b^{\dagger}(t)a(t+\tau)\rangle,\langle b^{\dagger}(t)b(t+\tau)% \rangle,\langle b^{\dagger}(t)c(t+\tau)\rangle,\langle b^{\dagger}(t)a^{% \dagger}(t+\tau)\rangle,\langle b^{\dagger}(t)b^{\dagger}(t+\tau)\rangle,% \langle b^{\dagger}(t)c^{\dagger}(t+\tau)\rangle)^{T}italic_v = ( ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_a ( italic_t + italic_τ ) ⟩ , ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_b ( italic_t + italic_τ ) ⟩ , ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_c ( italic_t + italic_τ ) ⟩ , ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_a start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t + italic_τ ) ⟩ , ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t + italic_τ ) ⟩ , ⟨ italic_b start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t ) italic_c start_POSTSUPERSCRIPT † end_POSTSUPERSCRIPT ( italic_t + italic_τ ) ⟩ ) start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT. In solving Eq. (13), the initial conditions (τ=0)𝜏0(\tau=0)( italic_τ = 0 ) are given by the steady-state solutions of the corresponding second moments. Therefore, the first-order correlation function is expressed as follows:

gb(1)(τ)=superscriptsubscript𝑔𝑏1𝜏absent\displaystyle g_{b}^{(1)}(\tau)=italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ) = i2ΠΩ3[e12(Γ+τ)ω(Ω3+iα+β+Ω)(βωcos(ω2τ)(α+βiΩ2)sin(ω2τ))\displaystyle\frac{i}{2\Pi\Omega^{3}}\Bigg{[}\frac{e^{-\frac{1}{2}(\Gamma_{+}% \tau)}}{\omega_{-}}\left(\Omega^{3}+i\alpha_{+}\beta_{+}\Omega\right)\left(% \beta_{-}\omega_{-}\cos\left(\frac{\omega_{-}}{2}\tau\right)-\left(\alpha_{+}% \beta_{-}-i\Omega^{2}\right)\sin\left(\frac{\omega_{-}}{2}\tau\right)\right)divide start_ARG italic_i end_ARG start_ARG 2 roman_Π roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG [ divide start_ARG italic_e start_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( roman_Γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_τ ) end_POSTSUPERSCRIPT end_ARG start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG ( roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + italic_i italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT + end_POSTSUBSCRIPT roman_Ω ) ( italic_β start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT roman_cos ( divide start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_τ ) - ( italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT - end_POSTSUBSCRIPT - italic_i roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) roman_sin ( divide start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_τ ) )
e12(Γτ)ω+(Ω3+iα+βΩ)(β+ω+cos(ω+2τ)+(α+β++iΩ2)sin(ω+2τ))],\displaystyle-\frac{e^{-\frac{1}{2}(\Gamma_{-}\tau)}}{\omega_{+}}\left(\Omega^% {3}+i\alpha_{+}\beta_{-}\Omega\right)\left(\beta_{+}\omega_{+}\cos\left(\frac{% \omega_{+}}{2}\tau\right)+\left(-\alpha_{+}\beta_{+}+i\Omega^{2}\right)\sin% \left(\frac{\omega_{+}}{2}\tau\right)\right)\Bigg{]},- divide start_ARG italic_e start_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG ( roman_Γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_τ ) end_POSTSUPERSCRIPT end_ARG start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG ( roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + italic_i italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT - end_POSTSUBSCRIPT roman_Ω ) ( italic_β start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT roman_cos ( divide start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_τ ) + ( - italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + italic_i roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) roman_sin ( divide start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG italic_τ ) ) ] , (14)

where Γ±γ±Π,α±γ±iΔ,β±Δ±iΠformulae-sequencesubscriptΓplus-or-minusplus-or-minus𝛾Πformulae-sequencesubscript𝛼plus-or-minusplus-or-minus𝛾𝑖Δsubscript𝛽plus-or-minusplus-or-minusΔ𝑖Π\Gamma_{\pm}\equiv\gamma\pm\Pi,\,\alpha_{\pm}\equiv\gamma\pm i\Delta,\,\beta_{% \pm}\equiv\Delta\pm i\Piroman_Γ start_POSTSUBSCRIPT ± end_POSTSUBSCRIPT ≡ italic_γ ± roman_Π , italic_α start_POSTSUBSCRIPT ± end_POSTSUBSCRIPT ≡ italic_γ ± italic_i roman_Δ , italic_β start_POSTSUBSCRIPT ± end_POSTSUBSCRIPT ≡ roman_Δ ± italic_i roman_Π. The optical spectrum is calculated by taking the Fourier transform of gb(1)(τ)superscriptsubscript𝑔𝑏1𝜏g_{b}^{(1)}(\tau)italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ). Since gb(1)(τ)superscriptsubscript𝑔𝑏1𝜏g_{b}^{(1)}(\tau)italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ) consists of four independent terms, the Fourier transform can be taken separately for each term. The total spectrum S(ω)𝑆𝜔S(\omega)italic_S ( italic_ω ) is expressed as follows [18, 19]:

S(ω)𝑆𝜔absent\displaystyle S(\omega)\equivitalic_S ( italic_ω ) ≡ 1πRe0gb(1)(τ)eiωτ𝑑τ=1πRefs(ω),1𝜋Resuperscriptsubscript0superscriptsubscript𝑔𝑏1𝜏superscript𝑒𝑖𝜔𝜏differential-d𝜏1𝜋Resubscript𝑓𝑠𝜔\displaystyle\frac{1}{\pi}\textrm{Re}\int\limits_{0}^{\infty}g_{b}^{(1)}(\tau)% e^{i\omega\tau}d\tau=\frac{1}{\pi}\textrm{Re}f_{s}(\omega),divide start_ARG 1 end_ARG start_ARG italic_π end_ARG Re ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_g start_POSTSUBSCRIPT italic_b end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT ( italic_τ ) italic_e start_POSTSUPERSCRIPT italic_i italic_ω italic_τ end_POSTSUPERSCRIPT italic_d italic_τ = divide start_ARG 1 end_ARG start_ARG italic_π end_ARG Re italic_f start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_ω ) , (15)

where

fs(ω)=subscript𝑓𝑠𝜔absent\displaystyle f_{s}(\omega)=italic_f start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ( italic_ω ) = i2ΠΩ3[Ω3+iα+β+Ωω[βωF1(α+βiΩ2)F2]Ω3+iα+βΩω+[β+ω+F3(α+β+iΩ2)F4]],𝑖2ΠsuperscriptΩ3delimited-[]superscriptΩ3𝑖subscript𝛼subscript𝛽Ωsubscript𝜔delimited-[]subscript𝛽subscript𝜔subscript𝐹1subscript𝛼subscript𝛽𝑖superscriptΩ2subscript𝐹2superscriptΩ3𝑖subscript𝛼subscript𝛽Ωsubscript𝜔delimited-[]subscript𝛽subscript𝜔subscript𝐹3subscript𝛼subscript𝛽𝑖superscriptΩ2subscript𝐹4\displaystyle\frac{i}{2\Pi\Omega^{3}}\left[\frac{\Omega^{3}+i\alpha_{+}\beta_{% +}\Omega}{\omega_{-}}\left[\beta_{-}\omega_{-}F_{1}-(\alpha_{+}\beta_{-}-i% \Omega^{2})F_{2}\right]-\frac{\Omega^{3}+i\alpha_{+}\beta_{-}\Omega}{\omega_{+% }}\left[\beta_{+}\omega_{+}F_{3}-(\alpha_{+}\beta_{+}-i\Omega^{2})F_{4}\right]% \right],divide start_ARG italic_i end_ARG start_ARG 2 roman_Π roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG [ divide start_ARG roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + italic_i italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT + end_POSTSUBSCRIPT roman_Ω end_ARG start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG [ italic_β start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - ( italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT - end_POSTSUBSCRIPT - italic_i roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) italic_F start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] - divide start_ARG roman_Ω start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT + italic_i italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT - end_POSTSUBSCRIPT roman_Ω end_ARG start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG [ italic_β start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_F start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT - ( italic_α start_POSTSUBSCRIPT + end_POSTSUBSCRIPT italic_β start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - italic_i roman_Ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) italic_F start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT ] ] , (16)

and

F1=subscript𝐹1absent\displaystyle F_{1}=italic_F start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = Γ+2iω(ω2)2+(Γ+2iω)2,F2=ω2(ω2)2+(Γ+2iω)2subscriptΓ2𝑖𝜔superscriptsubscript𝜔22superscriptsubscriptΓ2𝑖𝜔2subscript𝐹2subscript𝜔2superscriptsubscript𝜔22superscriptsubscriptΓ2𝑖𝜔2\displaystyle\frac{\frac{\Gamma_{+}}{2}-i\omega}{\left(\frac{\omega_{-}}{2}% \right)^{2}+\left(\frac{\Gamma_{+}}{2}-i\omega\right)^{2}},\quad F_{2}=\frac{% \frac{\omega_{-}}{2}}{\left(\frac{\omega_{-}}{2}\right)^{2}+\left(\frac{\Gamma% _{+}}{2}-i\omega\right)^{2}}divide start_ARG divide start_ARG roman_Γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - italic_i italic_ω end_ARG start_ARG ( divide start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( divide start_ARG roman_Γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - italic_i italic_ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_F start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = divide start_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG end_ARG start_ARG ( divide start_ARG italic_ω start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( divide start_ARG roman_Γ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - italic_i italic_ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG
F3=subscript𝐹3absent\displaystyle F_{3}=italic_F start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = Γ2iω(ω+2)2+(Γ2iω)2,F4=ω+2(ω+2)2+(Γ2iω)2.subscriptΓ2𝑖𝜔superscriptsubscript𝜔22superscriptsubscriptΓ2𝑖𝜔2subscript𝐹4subscript𝜔2superscriptsubscript𝜔22superscriptsubscriptΓ2𝑖𝜔2\displaystyle\frac{\frac{\Gamma_{-}}{2}-i\omega}{\left(\frac{\omega_{+}}{2}% \right)^{2}+\left(\frac{\Gamma_{-}}{2}-i\omega\right)^{2}},\quad F_{4}=\frac{% \frac{\omega_{+}}{2}}{\left(\frac{\omega_{+}}{2}\right)^{2}+\left(\frac{\Gamma% _{-}}{2}-i\omega\right)^{2}}.divide start_ARG divide start_ARG roman_Γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - italic_i italic_ω end_ARG start_ARG ( divide start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( divide start_ARG roman_Γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - italic_i italic_ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_F start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT = divide start_ARG divide start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG end_ARG start_ARG ( divide start_ARG italic_ω start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( divide start_ARG roman_Γ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG - italic_i italic_ω ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG .

As we can see from Eq. (16), the optical spectrum can be separated into four Lorentzian lineshapes. The spectrum is expressed further as follows:

Refer to caption
Figure 2: (a) and (b) The real and imaginary parts of λ𝜆\lambdaitalic_λ in the vicinity of the EP , ΩEP(1)superscriptsubscriptΩEP1\Omega_{\textrm{EP}}^{(1)}roman_Ω start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT, respectively. (c) The optical spectrum at the EP ΩEP(1)subscriptsuperscriptΩ1EP\Omega^{(1)}_{\textrm{EP}}roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT, Ω=ΔΩΔ\Omega=\Deltaroman_Ω = roman_Δ. Here we use Δ=2γΔ2𝛾\Delta=2\gammaroman_Δ = 2 italic_γ and J=0𝐽0J=0italic_J = 0. (d) - (f) The same as (a) - (c) but for J=0.25γ𝐽0.25𝛾J=0.25\gammaitalic_J = 0.25 italic_γ and (g) - (i) for J=0.5γ𝐽0.5𝛾J=0.5\gammaitalic_J = 0.5 italic_γ. The dashed lines in the optical spectra correspond to the peak positions ±2Jplus-or-minus2𝐽\pm\sqrt{2}J± square-root start_ARG 2 end_ARG italic_J. In all cases, we set γ=1𝛾1\gamma=1italic_γ = 1
S(ω)=2γω4(γ2Π2)π(ω4(γ4+2γ2(ω2Π2)+(Π2+ω2)2)+16J832J6ω2+8J4ω2(γ2+Π2+3ω2)8J2ω4(γ2+Π2+ω2))𝑆𝜔2𝛾superscript𝜔4superscript𝛾2superscriptΠ2𝜋superscript𝜔4superscript𝛾42superscript𝛾2superscript𝜔2superscriptΠ2superscriptsuperscriptΠ2superscript𝜔2216superscript𝐽832superscript𝐽6superscript𝜔28superscript𝐽4superscript𝜔2superscript𝛾2superscriptΠ23superscript𝜔28superscript𝐽2superscript𝜔4superscript𝛾2superscriptΠ2superscript𝜔2\displaystyle S(\omega)=\frac{2\gamma\omega^{4}\left(\gamma^{2}-\Pi^{2}\right)% }{\pi\left(\omega^{4}\left(\gamma^{4}+2\gamma^{2}\left(\omega^{2}-\Pi^{2}% \right)+\left(\Pi^{2}+\omega^{2}\right)^{2}\right)+16J^{8}-32J^{6}\omega^{2}+8% J^{4}\omega^{2}\left(\gamma^{2}+\Pi^{2}+3\omega^{2}\right)-8J^{2}\omega^{4}% \left(\gamma^{2}+\Pi^{2}+\omega^{2}\right)\right)}italic_S ( italic_ω ) = divide start_ARG 2 italic_γ italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG start_ARG italic_π ( italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( italic_γ start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + 2 italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - roman_Π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) + ( roman_Π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) + 16 italic_J start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT - 32 italic_J start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 8 italic_J start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_Π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 3 italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) - 8 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + roman_Π start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) ) end_ARG (17)

III Results and Discussion

Before looking at the optical spectrum, let us investigate λ𝜆\lambdaitalic_λ at the exceptional point with Ω=ΩEP(1)ΩsuperscriptsubscriptΩEP1\Omega=\Omega_{\textrm{EP}}^{(1)}roman_Ω = roman_Ω start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT. At an exceptional point, two or more eigenvalues and the corresponding eigenvectors coalesce. From Eq. (6), we have two distinct λ𝜆\lambdaitalic_λ when Ω=ΔΩΔ\Omega=\Deltaroman_Ω = roman_Δ. The λ3subscript𝜆3\lambda_{3}italic_λ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT and λ5subscript𝜆5\lambda_{5}italic_λ start_POSTSUBSCRIPT 5 end_POSTSUBSCRIPT coalesce with eigenvalues of,

λEP,3(1)=λEP,5(1)=iγ2+2J2(γ2)2subscriptsuperscript𝜆1EP,3subscriptsuperscript𝜆1EP,5𝑖𝛾22superscript𝐽2superscript𝛾22\displaystyle\lambda^{(1)}_{\textrm{EP,3}}=\lambda^{(1)}_{\textrm{EP,5}}=-i% \frac{\gamma}{2}+\sqrt{2J^{2}-(\frac{\gamma}{2})^{2}}italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,3 end_POSTSUBSCRIPT = italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,5 end_POSTSUBSCRIPT = - italic_i divide start_ARG italic_γ end_ARG start_ARG 2 end_ARG + square-root start_ARG 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - ( divide start_ARG italic_γ end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG (18)

and the corresponding coalesced eigenvectors of,

𝐮3=𝐮5=(1,λEP,3(1)J,1,1,λEP,3(1)J,1)T.subscript𝐮3subscript𝐮5superscript1subscriptsuperscript𝜆1EP,3𝐽11subscriptsuperscript𝜆1EP,3𝐽1𝑇\displaystyle\mathbf{u}_{3}=\mathbf{u}_{5}=\left(1,\frac{\lambda^{(1)}_{% \textrm{EP,3}}}{J},1,1,-\frac{\lambda^{(1)}_{\textrm{EP,3}}}{J},1\right)^{T}.bold_u start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT = bold_u start_POSTSUBSCRIPT 5 end_POSTSUBSCRIPT = ( 1 , divide start_ARG italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,3 end_POSTSUBSCRIPT end_ARG start_ARG italic_J end_ARG , 1 , 1 , - divide start_ARG italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,3 end_POSTSUBSCRIPT end_ARG start_ARG italic_J end_ARG , 1 ) start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT . (19)

On the other hand, the λ4subscript𝜆4\lambda_{4}italic_λ start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT and λ6subscript𝜆6\lambda_{6}italic_λ start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT coalesce with eigenvalues of,

λEP,4(1)=λEP,6(1)=iγ22J2(γ2)2subscriptsuperscript𝜆1EP,4subscriptsuperscript𝜆1EP,6𝑖𝛾22superscript𝐽2superscript𝛾22\displaystyle\lambda^{(1)}_{\textrm{EP,4}}=\lambda^{(1)}_{\textrm{EP,6}}=-i% \frac{\gamma}{2}-\sqrt{2J^{2}-(\frac{\gamma}{2})^{2}}italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,4 end_POSTSUBSCRIPT = italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,6 end_POSTSUBSCRIPT = - italic_i divide start_ARG italic_γ end_ARG start_ARG 2 end_ARG - square-root start_ARG 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - ( divide start_ARG italic_γ end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG (20)

and the corresponding coalesced eigenvectors of,

𝐮4=𝐮6=(1,λEP,4(1)J,1,1,λEP,4(1)J,1)T.subscript𝐮4subscript𝐮6superscript1subscriptsuperscript𝜆1EP,4𝐽11subscriptsuperscript𝜆1EP,4𝐽1𝑇\displaystyle\mathbf{u}_{4}=\mathbf{u}_{6}=\left(1,\frac{\lambda^{(1)}_{% \textrm{EP,4}}}{J},1,1,-\frac{\lambda^{(1)}_{\textrm{EP,4}}}{J},1\right)^{T}.bold_u start_POSTSUBSCRIPT 4 end_POSTSUBSCRIPT = bold_u start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT = ( 1 , divide start_ARG italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,4 end_POSTSUBSCRIPT end_ARG start_ARG italic_J end_ARG , 1 , 1 , - divide start_ARG italic_λ start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,4 end_POSTSUBSCRIPT end_ARG start_ARG italic_J end_ARG , 1 ) start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT . (21)

Therefore, there are two distinct EPs when we set Ω=ΔΩΔ\Omega=\Deltaroman_Ω = roman_Δ and J>0𝐽0J>0italic_J > 0, as shown by Fig. 2. On the other hand, in the case of J=0𝐽0J=0italic_J = 0, all eigenvalues and eigenvector coalesce into one EP creating a higher order EP as shown in Figs. 2 (a) and (b).

In the case of 2J2<(γ/2)22superscript𝐽2superscript𝛾222J^{2}<(\gamma/2)^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT < ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, λ𝜆\lambdaitalic_λ is purely imaginary as shown by Figs. 2 (d) and (e). The two EPs are clearly shown in the imaginary part of λ𝜆\lambdaitalic_λ, while the real part of λ𝜆\lambdaitalic_λ of all branches vanishes. In the case of 2J2>(γ/2)22superscript𝐽2superscript𝛾222J^{2}>(\gamma/2)^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT > ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, λ𝜆\lambdaitalic_λ is complex, with the two EPs clearly shown in the real part of λ𝜆\lambdaitalic_λ as shown in 2 (g) and (h). In all cases, J𝐽Jitalic_J determines the separation of the two EPs in λ𝜆\lambdaitalic_λ, with separation of ΔλEP=22J2(γ/2)2Δsubscript𝜆EP22superscript𝐽2superscript𝛾22\Delta\lambda_{\textrm{EP}}=2\sqrt{2J^{2}-(\gamma/2)^{2}}roman_Δ italic_λ start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = 2 square-root start_ARG 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG.

Refer to caption
Figure 3: Optical spectra for several ΩΩ\Omegaroman_Ω. (a) For J=0𝐽0J=0italic_J = 0. (b) For J=0.25γ𝐽0.25𝛾J=0.25\gammaitalic_J = 0.25 italic_γ. (c) For J=0.5γ𝐽0.5𝛾J=0.5\gammaitalic_J = 0.5 italic_γ. In all cases, we set Δ=2γΔ2𝛾\Delta=2\gammaroman_Δ = 2 italic_γ and γ=1𝛾1\gamma=1italic_γ = 1. The dashed lines in the optical spectra correspond to the peak positions ±2Jplus-or-minus2𝐽\pm\sqrt{2}J± square-root start_ARG 2 end_ARG italic_J. The EP is given by Ω=Δ=2γΩΔ2𝛾\Omega=\Delta=2\gammaroman_Ω = roman_Δ = 2 italic_γ.

Let us investigate the optical spectrum corresponding to these EPs. As we mentioned before, these EPs are generated by the same driving strength as in the case of single oscillator ΩEP(1)=ΔsubscriptsuperscriptΩ1EPΔ\Omega^{(1)}_{\textrm{EP}}=\Deltaroman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = roman_Δ, which does not depend on γ𝛾\gammaitalic_γ and J𝐽Jitalic_J. Therefore, it is intriguing to probe the optical spectrum at these EPs and look into how the coupling affects the optical spectrum.

At exactly Ω=ΔΩΔ\Omega=\Deltaroman_Ω = roman_Δ, the optical spectrum is expressed simply as

SEP,1(ω)=subscript𝑆EP,1𝜔absent\displaystyle S_{\textrm{EP,1}}(\omega)=italic_S start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ( italic_ω ) = 2γ3ω4π(((γ2)2+(ωωEP,1(1))2)((γ2)2+(ω+ωEP,1(1))2))2=2γ3ω4π(γ2ω2+(ω22J2)2)22superscript𝛾3superscript𝜔4𝜋superscriptsuperscript𝛾22superscript𝜔subscriptsuperscript𝜔1EP,12superscript𝛾22superscript𝜔subscriptsuperscript𝜔1EP,1222superscript𝛾3superscript𝜔4𝜋superscriptsuperscript𝛾2superscript𝜔2superscriptsuperscript𝜔22superscript𝐽222\displaystyle\frac{2\gamma^{3}\omega^{4}}{\pi\left(\left(\left(\frac{\gamma}{2% }\right)^{2}+(\omega-\omega^{(1)}_{\textrm{EP,1}})^{2}\right)\left(\left(\frac% {\gamma}{2}\right)^{2}+(\omega+\omega^{(1)}_{\textrm{EP,1}})^{2}\right)\right)% ^{2}}=\frac{2\gamma^{3}\omega^{4}}{\pi\left(\gamma^{2}\omega^{2}+\left(\omega^% {2}-2J^{2}\right)^{2}\right)^{2}}divide start_ARG 2 italic_γ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_π ( ( ( divide start_ARG italic_γ end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_ω - italic_ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) ( ( divide start_ARG italic_γ end_ARG start_ARG 2 end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_ω + italic_ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = divide start_ARG 2 italic_γ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_π ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG (22)

for the case of 2J2>(γ/2)22superscript𝐽2superscript𝛾222J^{2}>(\gamma/2)^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT > ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where we define ωEP,1(1)2J2(γ/2)2subscriptsuperscript𝜔1EP,12superscript𝐽2superscript𝛾22\omega^{(1)}_{\textrm{EP,1}}\equiv\sqrt{2J^{2}-(\gamma/2)^{2}}italic_ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ≡ square-root start_ARG 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG. The peaks of the spectrum are located at ωp=±(γ/2)2+(ωEP,1(1))2=±2Jsubscript𝜔𝑝plus-or-minussuperscript𝛾22superscriptsubscriptsuperscript𝜔1EP,12plus-or-minus2𝐽\omega_{p}=\pm\sqrt{(\gamma/2)^{2}+(\omega^{(1)}_{\textrm{EP,1}})^{2}}=\pm% \sqrt{2}Jitalic_ω start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT = ± square-root start_ARG ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = ± square-root start_ARG 2 end_ARG italic_J. For the case of 2J2<(γ/2)22superscript𝐽2superscript𝛾222J^{2}<(\gamma/2)^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT < ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where we define ωEP,1(2)(γ/2)22J2subscriptsuperscript𝜔2EP,1superscript𝛾222superscript𝐽2\omega^{(2)}_{\textrm{EP,1}}\equiv\sqrt{(\gamma/2)^{2}-2J^{2}}italic_ω start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ≡ square-root start_ARG ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG, the spectrum is expressed as follows:

SEP,1(ω)=subscript𝑆EP,1𝜔absent\displaystyle S_{\textrm{EP,1}}(\omega)=italic_S start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ( italic_ω ) = 2γ3ω4π(((γ2ωEP,1(2))2+4ω2)((γ+2ωEP,1(2))2+4ω2))2=2γ3ω4π(γ2ω2+(ω22J2)2)22superscript𝛾3superscript𝜔4𝜋superscriptsuperscript𝛾2subscriptsuperscript𝜔2EP,124superscript𝜔2superscript𝛾2subscriptsuperscript𝜔2EP,124superscript𝜔222superscript𝛾3superscript𝜔4𝜋superscriptsuperscript𝛾2superscript𝜔2superscriptsuperscript𝜔22superscript𝐽222\displaystyle\frac{2\gamma^{3}\omega^{4}}{\pi\left(\left((\gamma-2\omega^{(2)}% _{\textrm{EP,1}})^{2}+4\omega^{2}\right)\left((\gamma+2\omega^{(2)}_{\textrm{% EP,1}})^{2}+4\omega^{2}\right)\right)^{2}}=\frac{2\gamma^{3}\omega^{4}}{\pi% \left(\gamma^{2}\omega^{2}+\left(\omega^{2}-2J^{2}\right)^{2}\right)^{2}}divide start_ARG 2 italic_γ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_π ( ( ( italic_γ - 2 italic_ω start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 4 italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) ( ( italic_γ + 2 italic_ω start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 4 italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = divide start_ARG 2 italic_γ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG italic_π ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG (23)

The peaks are located at ωp=±(γ/2)2(ωEP,1(2))2=±2Jsubscript𝜔𝑝plus-or-minussuperscript𝛾22superscriptsubscriptsuperscript𝜔2EP,12plus-or-minus2𝐽\omega_{p}=\pm\sqrt{(\gamma/2)^{2}-(\omega^{(2)}_{\textrm{EP,1}})^{2}}=\pm% \sqrt{2}Jitalic_ω start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT = ± square-root start_ARG ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - ( italic_ω start_POSTSUPERSCRIPT ( 2 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP,1 end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = ± square-root start_ARG 2 end_ARG italic_J. Interestingly, in both cases, the final expressions of the spectrum are the same as given by Eqs. (22) and (23), giving the same peak positions, which do not depend on γ𝛾\gammaitalic_γ.

At the zero coupling J=0𝐽0J=0italic_J = 0, we recover the spectrum for a driven single oscillator [19],

SEPsingle=2γ3π(γ2+ω2)2,superscriptsubscript𝑆EPsingle2superscript𝛾3𝜋superscriptsuperscript𝛾2superscript𝜔22\displaystyle S_{\textrm{EP}}^{\textrm{single}}=\frac{2\gamma^{3}}{\pi\left(% \gamma^{2}+\omega^{2}\right)^{2}},italic_S start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT start_POSTSUPERSCRIPT single end_POSTSUPERSCRIPT = divide start_ARG 2 italic_γ start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG start_ARG italic_π ( italic_γ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , (24)

where the peak is located at ω=0𝜔0\omega=0italic_ω = 0 as shown in Fig. 2 (c). This peak position simply corresponds to λ=iγ𝜆𝑖𝛾\lambda=-i\gammaitalic_λ = - italic_i italic_γ as shown in Figs. 2 (a) and (b), where all branches coalesces at zero real part of λ𝜆\lambdaitalic_λ. The purely imaginary λ𝜆\lambdaitalic_λ gives the broadening of the spectrum centered on ω=0𝜔0\omega=0italic_ω = 0. On the other hand, when the coupling is finite, the spectrum is split into two symmetric peaks, with a peak separation of ΔωEP=22JΔsubscript𝜔EP22𝐽\Delta\omega_{\textrm{EP}}=2\sqrt{2}Jroman_Δ italic_ω start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = 2 square-root start_ARG 2 end_ARG italic_J. This peak separation is clearly shown in Figs. 2 (f) and (i), where we use J=0.25γ𝐽0.25𝛾J=0.25\gammaitalic_J = 0.25 italic_γ and J=0.5γ𝐽0.5𝛾J=0.5\gammaitalic_J = 0.5 italic_γ, respectively. From Eqs. (22) and (23), we also find that the intensity of the spectrum vanishes for ω=0𝜔0\omega=0italic_ω = 0 as soon as J𝐽Jitalic_J is finite, which distinguishes it from the spectrum of a single oscillator. More interestingly, this peak separation is different from the separation of EP eigenvalues ΔλEPΔsubscript𝜆EP\Delta\lambda_{\textrm{EP}}roman_Δ italic_λ start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT. In particular, in the case of 2J2<(γ/2)22superscript𝐽2superscript𝛾222J^{2}<(\gamma/2)^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT < ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, where the separation of EP appears in the imaginary λ𝜆\lambdaitalic_λ, the two peaks emerge, instead of one peak like the case of the single oscillator in (c).

To better understand the optical spectrum, we calculate the spectrum for several ΩΩ\Omegaroman_Ω as shown in Fig. 3. The spectrum for a single oscillator is shown in Fig. 3(a). There appear two symmetric peaks when Ω=1/2ΩEP(1)=γΩ12subscriptsuperscriptΩ1EP𝛾\Omega=1/2\Omega^{(1)}_{\textrm{EP}}=\gammaroman_Ω = 1 / 2 roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = italic_γ, which correspond to the two branches of the real part of λ𝜆\lambdaitalic_λ. Approaching EP, the two peaks merge into a single broad peak. A clear single Lorentzian spectrum starts to appear at the EP and stays centered at ω=0𝜔0\omega=0italic_ω = 0 for Ω>ΩEP(1)=ΔΩsubscriptsuperscriptΩ1EPΔ\Omega>\Omega^{(1)}_{\textrm{EP}}=\Deltaroman_Ω > roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = roman_Δ, as expected since the real part of λ𝜆\lambdaitalic_λ vanishes for all branches. Unlike a single oscillator, when J𝐽Jitalic_J is finite, we have four peaks instead of two, as shown in Figs. 3 (b) and (c), which corresponds to the four distinct real parts of λ𝜆\lambdaitalic_λ when Ω=1/2ΩEP(1)=γΩ12subscriptsuperscriptΩ1EP𝛾\Omega=1/2\Omega^{(1)}_{\textrm{EP}}=\gammaroman_Ω = 1 / 2 roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = italic_γ. Approaching the EP, each of the two peaks merges into a single peak at ω=±2J𝜔plus-or-minus2𝐽\omega=\pm\sqrt{2}Jitalic_ω = ± square-root start_ARG 2 end_ARG italic_J and stays at the same positions even for Ω>ΩEP(1)=ΔΩsubscriptsuperscriptΩ1EPΔ\Omega>\Omega^{(1)}_{\textrm{EP}}=\Deltaroman_Ω > roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT = roman_Δ. Thus, the change from quadruplet to doublet in the optical spectrum is the signature of passing EPs in this coupled oscillator. In addition to the change of the number of peaks as a signature of passing EPs, the appearance of two peaks at EP implies the presence of a coupling in the oscillator (compared with the case of a single oscillator).

This peculiar position of the peaks comes from the presence of a dispersive term ω4superscript𝜔4\omega^{4}italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT on the numerator in Eqs. (22) and (23). Without ω4superscript𝜔4\omega^{4}italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT, the expression of the optical spectrum will be pure Lorentzian, with the peak position corresponding to the real part of λ𝜆\lambdaitalic_λ. In particular, we would have a single peak at ω=0𝜔0\omega=0italic_ω = 0 for 2J2<(γ/2)22superscript𝐽2superscript𝛾222J^{2}<(\gamma/2)^{2}2 italic_J start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT < ( italic_γ / 2 ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT instead of a doublet in Fig. 2(f). In the presence of ω4superscript𝜔4\omega^{4}italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT, the spectrum at ω=0𝜔0\omega=0italic_ω = 0 vanishes, which shifts the peak positions of the pure Lorenztian spectrum to ω=±2J𝜔plus-or-minus2𝐽\omega=\pm\sqrt{2}Jitalic_ω = ± square-root start_ARG 2 end_ARG italic_J. Interestingly, these peak position corresponds to the natural frequency of the trimer system in a closed system. Therefore, operating at the EP ΩEP(1)subscriptsuperscriptΩ1EP\Omega^{(1)}_{\textrm{EP}}roman_Ω start_POSTSUPERSCRIPT ( 1 ) end_POSTSUPERSCRIPT start_POSTSUBSCRIPT EP end_POSTSUBSCRIPT, the coupling strength between oscillators can be estimated from the peaks of optical spectrum.

IV Conclusion

We investigated the optical spectrum of a trimer chain driven by a quadratic photon. We have shown the emergence of multiple exceptional points corresponding to the coalescence of eigenvalues and eigenvectors of the dynamical matrix. In a particular exceptional point Ω=ΔΩΔ\Omega=\Deltaroman_Ω = roman_Δ, the resulted optical spectra are expressed as a product of a conventional Lorenztian spectrum with a dispersive term of ω4superscript𝜔4\omega^{4}italic_ω start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT, which is absent in the case of a single oscillator. Unlike the case of a single oscillator, the presence of coupling gives rise to the emergence of two peaks, instead of a single peak. In particular, those two peaks are located at the natural frequency of the trimer in a closed system of ω=±2J𝜔plus-or-minus2𝐽\omega=\pm\sqrt{2}Jitalic_ω = ± square-root start_ARG 2 end_ARG italic_J proportional to the coupling strength J𝐽Jitalic_J, regardless of the spectra of the eigenvalues. Furthermore, we have found that these peak positions persist for Ω>ΔΩΔ\Omega>\Deltaroman_Ω > roman_Δ. Therefore, the coupling between the oscillator can be estimated from the splitting of the two peaks.

Acknowledgements.
We thank Mahameru BRIN for their HPC facilities. A. A. is supported by research assistantship from BRIN Directorate for Talent Management. We thank Dr. Charles Downing for a fruitful discussion.

References

  • Wiersig [2020] J. Wiersig, Review of exceptional point-based sensors, Photonics Research 8, 1457 (2020).
  • Miri and Alu [2019] M.-A. Miri and A. Alu, Exceptional points in optics and photonics, Science 363, eaar7709 (2019).
  • Ashida et al. [2020] Y. Ashida, Z. Gong, and M. Ueda, Non-hermitian physics, Advances in Physics 69, 249 (2020).
  • Özdemir et al. [2019] Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Parity–time symmetry and exceptional points in photonics, Nature materials 18, 783 (2019).
  • Heiss et al. [1998] W. Heiss, M. Müller, and I. Rotter, Collectivity, phase transitions, and exceptional points in open quantum systems, Physical Review E 58, 2894 (1998).
  • Downing and Saroka [2021] C. A. Downing and V. A. Saroka, Exceptional points in oligomer chains, Communications Physics 4, 254 (2021).
  • Downing and Ukhtary [2023] C. A. Downing and M. S. Ukhtary, A quantum battery with quadratic driving, Communications Physics 6, 322 (2023).
  • Hodaei et al. [2017] H. Hodaei, A. U. Hassan, S. Wittek, H. Garcia-Gracia, R. El-Ganainy, D. N. Christodoulides, and M. Khajavikhan, Enhanced sensitivity at higher-order exceptional points, Nature 548, 187 (2017).
  • Chen et al. [2017] W. Chen, Ş. Kaya Özdemir, G. Zhao, J. Wiersig, and L. Yang, Exceptional points enhance sensing in an optical microcavity, Nature 548, 192 (2017).
  • Wiersig [2014] J. Wiersig, Enhancing the sensitivity of frequency and energy splitting detection by using exceptional points: application to microcavity sensors for single-particle detection, Physical review letters 112, 203901 (2014).
  • Wiersig [2016] J. Wiersig, Sensors operating at exceptional points: General theory, Physical review A 93, 033809 (2016).
  • Kuo et al. [2020] P.-C. Kuo, N. Lambert, A. Miranowicz, H.-B. Chen, G.-Y. Chen, Y.-N. Chen, and F. Nori, Collectively induced exceptional points of quantum emitters coupled to nanoparticle surface plasmons, Physical Review A 101, 013814 (2020).
  • Manzano [2020] D. Manzano, A short introduction to the lindblad master equation, Aip advances 10 (2020).
  • Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione, The theory of open quantum systems (Oxford University Press, USA, 2002).
  • Gardiner and Zoller [2014] C. Gardiner and P. Zoller, The Quantum World of Ultra-Cold Atoms and Light, Book I: Foundations of Quantum Optics (Imperial College Pressy, 2014).
  • Arkhipov et al. [2020] I. I. Arkhipov, A. Miranowicz, F. Minganti, and F. Nori, Liouvillian exceptional points of any order in dissipative linear bosonic systems: Coherence functions and switching between pt and anti-pt symmetries, Physical Review A 102, 033715 (2020).
  • Perina Jr et al. [2022] J. Perina Jr, A. Miranowicz, G. Chimczak, and A. Kowalewska-Kudlaszyk, Quantum liouvillian exceptional and diabolical points for bosonic fields with quadratic hamiltonians: The heisenberg-langevin equation approach, Quantum 6, 883 (2022).
  • Kavokin et al. [2017] A. Kavokin, J. J. Baumberg, G. Malpuech, and F. P. Laussy, Microcavities (Oxford university press, 2017).
  • Downing and Vidiella-Barranco [2023] C. Downing and A. Vidiella-Barranco, Parametrically driving a quantum oscillator into exceptionality, Scientific reports 13, 11004 (2023).
  • Downing et al. [2023] C. Downing, E. Del Valle, and A. I. Fernández-Domínguez, Resonance fluorescence of two asymmetrically pumped and coupled two-level systems, Physical Review A 107, 023717 (2023).
  • Khan et al. [2022] S. Khan, B. K. Agarwalla, and S. Jain, Quantum regression theorem for multi-time correlators: A detailed analysis in the heisenberg picture, Physical Review A 106, 022214 (2022).