Dynamical $N$-photon bundle emission

Engineering multiphoton resources is of importance in quantum metrology, quantum lithography, and biological sensing. Here we propose a concept of dynamical emission of $N$ strongly-correlated photons. This is realized in a circuit quantum electrodynamical system driven by two Gaussian-pulse sequences. The underlying physical mechanism relies on the stimulated Raman adiabatic passage that allows efficient and selective preparation of target multiphoton states. Assisted by the photon decay, a highly pure $N$-photon bundle emission takes place in this system. In particular, the dynamical $N$-photon bundle emission can be tuned by controlling the time interval between consecutive pulses so that the device behaves as an $N$-photon gun, which can be triggered on demand. Our work opens up a route to achieve multiphoton source devices, which have wide potential applications in quantum information processing and quantum metrology.

Physically, to realize an efficient N -photon emission, it is desired to generate a perfect N -photon number states in advance. From an application perspective, it is expected to create the states in a dynamical and deterministic manner. Here, the dynamical requirement confirms that the photon states can be generated under control and that the preparation of the photon states can be completed on demand. Meanwhile, the deterministic way ensures the high efficiency of the state generation.
Here, we propose a concept of dynamical N -photon bundle emission and present a feasible scheme to generate dynamical N -photon (N = 2, 3) bundle emission in a circuit-QED system, in which a microwave resonator is longitudinally coupled to a qubit driven by two Gaussian-pulse sequences. Based on the feature of stimulated Raman adiabatic passage (STIRAP) [47][48][49], the population transfer between zero-and N -photon states can be realized. Combined with photon decay of the resonator, the N -photon states can be emitted as photon bundles. This means that the dynamical N -photon bundle emission is realized via the dissipative channel. This dynamical Nphoton bundle emission can be tuned on demand by controlling the time interval between consecutive pulses, then the device behaves as an N -photon gun, with wide applications for quantum science and technology.

Model
We consider a circuit-QED system, which is composed of a qubit longitudinally coupled to a microwave resonator [ Fig. 1(a)]. The Hamiltonian of the system reads (h = 1) where b † (b) is the creation (annihilation) operator of the microwave resonator with resonant frequency ω b . The operators σ + = |e g| and σ − = |g e| are, respectively, the raising and lowering operators of the qubit, with transition frequency ω 0 between the excited state |e and the ground state |g . The last term in Eq. (1) denotes the interaction between the qubit and the resonator with the coupling strength λ. Note that in Eq. (1) we have displaced the photonic field by introducing a driving λ(b † + b)/2 to the resonator, as an assistance of the longitude qubit-resonator coupling λσ z (b † + b)/2, such that the generated photon bundles are not displaced by the qubit. This point can be seen based on the relation Note that the physical effect induced by the longitudinal qubit-resonator interaction has been studied both theoretically [50][51][52][53][54][55] and experimentally [56][57][58][59].
In the case of without the driving λ(b † + b)/2 of the resonator, the eigenstates of the Hamiltonian being the displacement operator. In this case, the generated photon bundles will be coherently displaced.
To study the dynamical bundle emission of n strongly-correlated photons in this system, the qubit is driven by two Gaussian-pulse sequences with corresponding carrier frequencies ω 1 and ω 2 . Each Gaussian wave packet of the two Gaussian-pulse sequences has the same amplitude and width. The driving Hamiltonian is with Here, Ω 0 and √ 2σ are the amplitude and width of each Gaussian wave packet, respectively. The t i + mT denotes the time when the pulse Ω i (t) reaches its maximum value, where m is an integer and T is the time interval between consecutive pulses. N 0 + 1 is the number of the pulses. In the rotating frame with respect to ω 0 σ + σ − , the total Hamiltonian of the system becomes (4) where ∆ i = ω i − ω 0 (i = 1, 2) is the detuning of the driving carrier frequency ω i with respect to the transition frequency ω 0 of the qubit.

Approximate Hamiltonian
In the resolved-sideband regime (i.e., the photon frequency ω b is much larger than the decay rate γ of the qubit, ω b /γ 1), we choose the proper driving carrier frequencies ω 1 and ω 2 such that the detunings ∆ 1 and ∆ 2 satisfy the resonant conditions: In this case, the Hamiltonian H I can be decomposed into two parts: where H I corresponds to the resonant-transition part [Ω 1 (t) ñ|n |e,ñ g, n| + Ω 2 (t) ñ|n + 1 |e,ñ g, n + 1| + H.c.], (10) and H I denotes the off-resonant-transition part Here the primed summation in H I eliminates these terms of H I , and the off-resonance detunings are δ i (n, m) = δE n,m − ∆ i for i = 1, 2. When the off-resonance detunings |δ i (n, m)| (i = 1, 2) are much greater than Ω 0 | ñ|m |, then the high-frequency oscillating Hamiltonian H I can be neglected by the rotating-wave approximation, namely, H I ≈ H I . Since | ñ|m | ≤ 1 [61], the condition of neglecting H I is reduced to Ω 0 ω b . In addition, we choose the proper coupling strength λ = λ N , where λ N is the minimal positive value for satisfying the equation [61] Ñ with N being a positive integer. Equation (12) denotes that the transition matrix element Ω 1 (t) Ñ |N between the states |g, N and |e,Ñ is zero. Hence, the dimension of the photon Hilbert space can be truncated up to n = N . To investigate how the coupling strength λ N depends on N , we plot in Fig. 1(b) the coupling strength λ N /ω b as a function of the index N . It can be seen that the coupling strength λ N /ω b decreases as the index N increases. Under the above mentioned conditions, the system can be well described by the approximate Hamiltonian [Ω 1 (t) ñ|n |e,ñ g, n|+Ω 2 (t) ñ|n+1 |e,ñ g, n+1|+H.c.], (13) which describes the resonant transition chain: |g, 0 ↔ |e,0 ↔ |g, 1 ↔ |e,1 ↔ · · · |e, N − 1 ↔ |g, N , as shown in Fig. 1(c). In Fig. 1(d), we show the two pulsed driving fields Ω i=1,2 (t) (m = 0) and the state populations P |g,n (t) = | g, n|ψ(t) | 2 as functions of ω b t at λ/ω b = λ 2 /ω b ≈ 0.765, i.e., corresponding to 2 |2 = 0. Here we consider that the initial state of the system is |g, 0 . We see a perfect population transfer from |g, 0 to |g, 2 in the absence of the system dissipation when the two Gaussian pulses Ω i=1,2 (t) satisfies an adiabatic evolution. This state transfer is determined by the physical mechanism of STIRAP [47][48][49]. To include the dissipation of the system, we consider the case where the resonator and the qubit are connected with two individual vacuum baths. Then the system is governed by the dressed-state master equation working in the ultrastrong-coupling regime [62][63][64][65][66][67][68] where H is given by Eq. (4) and κ (γ) is the decay rate of the resonator (qubit). The Lindblad superoperators are defined by Fig. 2, we show the state populations P |g,n (t) as functions of ω b t at various values of κ/ω b . Here we only show one period of the STIRAP pulses. It can be seen from Figs. 2(a)-2(c) that the maximum value of the population P |g,2 (t) is less than 1 in the presence of dissipation. Owing to the photon dissipation, the two photons are emitted out of the cavity, the system is then brought back to the initial state |g, 0 . The state |g, 2 is again generated for the next Gaussian pulse. Hence, the two Gaussian-pulse sequences lead to a transition |g, 0 ↔ |g, 2 under the assistant of the photon dissipation. In addition, we find that the maximum value of the population P |g,2 (t) decreases as κ/ω b increases. Meanwhile, for a smaller photon dissipation rate, a longer relaxation time is needed to reach the steady state. It should be noted that, in principle, the lowing operator associated with the coupled system in the ultrastrong-coupling regime should be b + (λ/ω b )σ + σ − rather than the annihilation operator b of the resonator. However, in the STIRAP scheme, the system is mainly in the state |g, n , and the population of the state |e,ñ is approximately 0. For example, we can see from Fig. 1(d) that the population is P |g,0 + P |g,1 + P |g,2 ≈ 1 during the STIRAP, which indicates that P |e,ñ ≈ 0, i.e., σ + σ − ≈ 0. In this case, it is reasonable to approximately calculate the photon statistics with the annihilation operator b instead of the displaced operator b + (λ/ω b )σ + σ − .

Dynamical N -photon bundle emission
To confirm the dynamical N -photon bundle emission, we employ a quantum Monte-Carlo approach to follow individual trajectory of the system [11,15,40]. Figure 3(a) shows a quantum trajectory of dynamical two-photon bundle emission under the driving of two Gaussian-pulse sequences. Here P |g,n (t) are the populations of the states |g, n  for n = 0, 1, 2. Initially, we consider that the system is in state |g, 0 . At time t ≈ t s , the system is transferred via the STIRAP from the initial state |g, 0 to the two-photon state |g, 2 with a prefect probability. The dissipation of the resonator triggers a quantum collapse of the system, from the two-photon state |g, 2 to the one-photon state |g, 1 with an almost unit probability, which causes the emission of the first photon. Subsequently, the system transits from the one-photon state |g, 1 to the zero-photon state |g, 0 , completing the two-photon emission within the resonator lifetime. After the two-photon emission, the system goes back to the initial state |g, 0 . Next, the system is again prepared in |g, 2 by the next Gaussian pulses, as the starting of the next emission of two photons. Hence, the dynamical two-photon bundle emission can be realized based on the sequential STIRAP. Here we choose the time T between coterminous pulses is much greater than the lifetime 1/κ of the resonator such that the system can go back to the initial state |g, 0 before the arrival of the next Gaussian pulses. In principle, the time duration between two STIRAPs can be controlled such that the emission of two photons can be triggered on demand.
We also prove the two-photon emission by checking both the density matrix and the Wigner function of the photon state. Figures 3(b) and 3(c) show the full density matrix elements |ρ ss | of the system and the Wigner functions W b of the reduced density matrix of the photon mode at three different moments. It can be seen that the cascadephoton-emission process |g, 2 → |g, 1 → |g, 0 occurs in a very short time window. Figure 4 shows a quantum trajectory of the average photon number b † b (t) at different values of λ/ω b and Ω 0 /ω b . In Figs. 4(a) and 4(b), we observe that the photon mode evolves from |N to |0 (N = 2, 3) during each cycle of bundle emission. Owing to the dissipation of the resonator, the system undergoes a rapid cascade emission through the series of the Fock states |m where 0 ≤ m ≤ N in very short temporal windows. These results indicate that the dynamical N -photon (N = 2, 3) bundle emission is realized in this system under the combination of the sequential STIRAP and the photon dissipation.
We point out that during the N -photon bundle emission process, there is no wavepacket driving. Therefore, the bundle emission of the N photons generated by the STIRAP scheme can be approximately understood as the N -photon emission physical process in a free cavity. Namely, for each bundle emission, there are two steps. The first step is the generation of N -photon Fock state in the cavity by the STIRAP, and the second step is the emitting of N photons by the decay process of the free cavity. According to the quantum optics theory, we know that the emissions of the N photons are independent processes and the N emitted photons in the outside fields are in N independent wave packets in frequency space [Its wave function in position space can be expressed as N j=1 where κ is the decay rate of the free cavity, x j is the coordinate of the jth photon in the outside fields, v g is the group velocity of the photon in the outside fields, and θ(t) is the Heaviside step function]. This result is understandable in physics because the cavity is free. Therefore, we can define a single temporal mode (namely N independent wave packets in frequency space) in the output fields to describe the state of the N emitted photons.
To characterize the statistical properties of the dynamical N -photon bundle emission, we consider the generalized nth-order correlation functions defined by [11]: where T ± represents the time-ordering operators and we neglect the time interval within photon bundle. Note that the standard correlation function g (n) (t 1 , . . . , t n ) corresponds to the case of N = 1 in the generalized correlation function (16). In Fig. 5(a), we show the evolution of the equal-time second-order correlation functions g 1 (t, t) and g (2) 2 (t, t) at λ/ω b ≈ 0.765. We can see that the correlation functions periodically change under the action of two Gaussian-pulse sequences. Figure 5(b) shows one period of the correlation functions g 1 (t, t) and g (2) 2 (t, t). It can be found that g (2) 1 (t, t) reaches maximum at time t = t s1 and g (2) 1 (t s1 , t s1 ) > 1, which indicates the super-Poisson distribution of a single photon. Moreover, g 2 (t, t) reaches its minimum at time t = t s2 , and g  1 (t s1 , t s1 + τ ) and g (2) 2 (t s2 , t s2 +τ ) with delayed time τ . It can be seen that g 1 (t s1 , t s1 ) > g (2) 1 (t s1 , t s1 +τ ) and g (2) 2 (t s2 , t s2 ) < g (2) 2 (t s2 , t s2 + τ ), which shows bunching of a single photon, but antibunching between two-photon bundle. In addition, we also study the statistical properties of the dynamical three-photon bundle emission in Figs. 5(d) and 5(e) when λ/ω b ≈ 0.645. It can be seen from Fig. 5(d) that g (2) 1 (t q1 , t q1 ) > 1 [g (2) 3 (t q3 , t q3 ) < 1], which indicates the super-Poisson (sub-Poisson) distribution of a single photon (threephoton bundle). In Fig. 5(e), we observe that g (2) 1 (t q1 , t q1 ) > g (2) 1 (t q1 , t q1 + τ ) and g (2) 3 (t q3 , t q3 ) < g (2) 3 (t q3 , t q3 +τ ). This means bunching of a single photon and antibunching between three-photon bundle.         It should be pointed out that the photon emission in this system occurs due to the intrinsic decay of the resonator, which means that the time of each photon emission has an inherent uncertainty [69]. However, in this work we consider the case where the time interval between any two consecutive STIRAP pulses is much larger than the time duration of N -photon emission in a bundle. Therefore, the time interval within the N photons can be ignored on the timescales of two Gaussian-pulse sequences. In this sense, the N photons can be approximately considered as a photon bundle. This point can also be confirmed from the correlation function (16), which indicates that the N photons are emitted simultaneously.
In conclusion, we have proposed the dynamical N -photon bundle emission and designed a feasible scheme to implement the physical process in a circuit-QED system based on the feature of the STIRAP. We have demonstrated that the dynamical twophoton (three-photon) bundle emission behaves as the bunching of a single photon and antibunching between two-photon (three-photon) bundle when the time interval between the consecutive pulses is much larger than the photon lifetime T 1/κ. Our work opens up a route to achieve multiphoton source emitter, which are very useful for quantum information processing and medical applications.
with two independent heat baths, then the Hamiltonian of the whole system including the circuit-QED system and its baths reads where H sd = H s + H d and the Hamiltonian related to the baths are given by Here the creation and annihilation operators b † j (c † k ) and b j (c k ) describe the jth (kth) mode with resonance frequency ω j (ω k ) of the resonator (qubit) bath, and λ b j (λ σ k ) is the coupling strength between the resonator (qubit) and the jth (kth) mode of the corresponding heat bath. In the interaction picture with respect to H 0 = H sd +H B b +H B σ , the equation of motion for the density operator ρ SB (t) of the whole system is given bẏ , where the interaction Hamiltonian H I (t) is given by In Eq. (A.4), we introduce the operators Note that in the derivation of Eq. (A.5), we have ignored the driving term H d when Ω 0 {ω 0 , ω b } [15,66], namely, we take in the derivation of the quantum master equation.
Under the Born-Markov approximation [67], the quantum master equation of the reduced density matrix ρ I (t) = Tr B [ρ SB (t)] in the interaction picture can be derived aṡ where Tr B denotes the trace operation over the bath mode and ρ B is the initial density operator of the bath. By substituting H I (t) and H I (t−s) into Eq. (A.6) and applying the rotating-wave approximation to neglect the fast-oscillating terms, the quantum master equation of the reduced density matrix ρ I (t) in the interaction picture can be obtained asρ where the two dissipation parts are given by In Eqs. (A.8) and (A.9), the correlation functions of the resonator and qubit baths are, respectively, defined as and wheren(ω j , respectively, the average occupation number associated with the resonator and the qubit. Below we will derive the contributions from the resonator and qubit baths in detail.

Resonator bath contribution-For the Hamiltonian
, its eigenvalues are given by ε g,n = nω b and ε e,n = nω b + ω 0 − λ 2 /ω b , with the corresponding eigenstates |g, n and |e,ñ . By using the eigenstates of the Hamiltonian H s , the time-dependent operatorb(t) can be expressed as |e,m e,m|, (A.12) where we use the relation b|m = √ m| m − 1 − (λ/ω b )|m . Therefore, we obtain the relationsb Assume that the spectral density of the resonator bath is Ohmic, i.e., J b (ω) = j |λ b j | 2 δ(ω − ω j ). It appears as J b (ω) = κω/(2πω b ) in the continuum limit of bath frequency, where κ = 2πJ b (ω b ) is the decay rate of the resonator. Using the formula [68] and omitting the fast-oscillating terms, the contribution of the resonator bath can then be obtained as Under the condition of ω 0 ω b , we further assume that the spectral density of the qubit bath, defined as J 0 (ω) = k |λ σ k | 2 δ(ω − ω k ), is a slow varying function in the vicinity of ω = ω 0 . Then the spectral density can be approximated as J 0 (ω) ≡ γ/2π in the entire range of the photon sidebands. In this case, we can obtain the following relations [68]: In the rotating frame with respect to ω 0 σ + σ − , the dressed-state master equation for the density operator ρ in the ultrastrong-coupling regime can be expressed as where H is given by Eq. (4).

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).