{'fraction_non_alphanumeric': 0.04527260715812381, 'fraction_numerical': 0.01830334589826293, 'mean_word_length': 4.51994851994852, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Photon and phonon statistics in a qubit-plasmon-phonon ultrastrong coupling system

We study photon/phonon statistics of a qubit-plasmon-phonon hybrid system in the ultrastrong coupling regime. The introduced qubit coupling causes parity conserving and non-conserving situations. We employ an analytic approximation approach for the parity conserving case to reveal the statistical behaviors of photons and phonons. It indicates that both photons and phonons show strong antibunching at the same frequency. Even though the bunching properties of photons/phonons occupy the dominant regions of the considered frequencies, phonons tend to weakly antibunching within the photonic strong-bunching area. In contrast, one can find that the configurations of correlation functions for both photons and phonons in the parity conserving case are squeezed towards the central frequency by parity breaking, which directly triggers the reverse statistical behaviors for the different parties at the low-frequency regions and the strong bunching properties at other frequency regions. The photon-phonon cross-correlation function also demonstrates similar parity-induced differences, indicating that the non-conserving parity induces the photon-phonon bunching behavior. We finally analyze the delayed second-order correlation function with different driving frequencies, which illustrates striking oscillations revealing the occurrence of simultaneous multiple excitations.


I. INTRODUCTION
Photon blockade has received much attention in recent years due to its nonclassical behavior of a quantum emitter, which leads to the realization of single-photon sources and single-photon detectors [1][2][3][4][5][6][7]. Photon blockade is the nonlinear excitation of the first photon with a high probability of blocking the transmission of the second photon [8,9] so that the emitted photons have a strong antibunching trend. With nanotechnologies, phonon lasers and single-photon generations attract wide interest. In particular, photons and phonons exhibit mutual antibunching in their correlated behavior in some photon-phonon hybrid systems. The output photon, phonon or their correlated statistics can be characterized by the equal-time second-order correlation function g (2) (0) and two-time correlation function g (2) (τ ), which can reveal nonclassical characters of the fields [10,11].
Cavity quantum electrodynamics is a powerful platform for studying the interactions between light and matter and between photons. The typical candidate is the Jaynes-Cummings (JC) model. In such a qubit-cavity system, the ratio of the qubit-cavity coupling constant g and the cavity-mode frequency ω 0 characterize the coupling strength η, which is closely related to the nonlinear effects in the system [12,13]. In recent years, the ultrastrong coupling regime with η > 0.1 has been achieved [12][13][14][15][16][17][18][19][20][21]. In this regime, the rotating wave approximation is no longer valid, which can give many intriguing physical effects [22][23][24][25][26][27], and much extensive research has been * Electronic address: ycs@dlut.edu.cn carried out regarding the modification of weak-coupling case quantum phenomena [28][29][30][31][32]. In particular, it is found that parametric processes induced by strong coupling can greatly influence photon blockade [3]. Recently, it has also been shown that one photon could simultaneously excite two or more atoms in the ultrastrong coupling regime [34]. These effects in the ultrastrong coupling regime can open up applications in quantum information processing [27,31,33].
The ultrastrong plasmon-phonon coupling has been realized via epsilon-near-zero (ENZ) nanocavities, which can drastically reduce the size of the system and thus the amount of material involved in the realization of midinfrared ultrastrong coupling [35]. It is shown that the coupling strength between the ENZ mode of the cavity and SiO 2 phonon with a normalized coupling strength can be larger than 0.25. In particular, it provides a new form of photon and phonon coupling in experiments, where the coaxial ENZ mode is coupled to the lattice vibrations of SiO 2 in a way different from the traditional optomechanical coupling mechanism. In addition, the ENZ nanocavity has also been designed and studied in many aspects [36][37][38][39][40][41][42][43], which can provide a new way to explore the quantum nonlinear optical processes. This paper theoretically studies photon (phonon) statistics in the ENZ nanocavity, especially concerning the new optomechanical coupling mechanism. To investigate the nonlinear optical effects, we introduce a single two-level atom into the ultrastrong photon-phonon coupling regime and show the photon and phonon statistics in the parity conserving and non-conserving cases induced by the atom-photon coupling. The nonlinearity of the cavity is mediated by the two-level atom. In such an ultrastrong coupling regime, all counter-rotating wave contributions shouldn't be neglected [12,13], and the equal-time correlation function g (2) (0) and two-time correlation function g (2) (τ ) will also be changed compared with the usual forms [3,5]. We find that in the parity conserving case, photons and phonons demonstrate strong antibunching behavior at the same driving frequency ω d and approach to coherent state at the high-frequency region ω d ∈ [1.3ω 0 , 1.4ω 0 ]. In the middle-frequency region ω d ∈ [0.9ω 0 , 1.15ω 0 ], photons tend to be strong bunching, but phonons are inclined to weakly anti bunch. Interestingly, non-conserving parity leads to the shift of eigenenergies. As a result, the configurations of correlation functions g (2) (0) and intensities (average particle numbers) are squeezed towards the middle frequency region, which leads to quite different and even reverse statistical behaviors compared with the parity conserving case. Within the range of relatively high frequencies ω d ∈ [1.2ω 0 , 1.4ω 0 ], the photons/ phonons are more inclined to bunch in the parity non-conserving case. In addition, the delayed second-order correlation functions through different driving frequencies illustrate striking oscillations, which reveals the simultaneous multiple excitations. The remainder of this paper is organized as follows. Sec. II introduces our considered model and derives the corresponding dynamical equation. In Sec. III, we employ an analytic approximation to study the photon/phonon statistics in the parity conserving case. In Sec. IV, we mainly consider the parity non-conserving case. We finally give the conclusions and discussion in Sec. V. The introduced qubit coupling can lead to parity conserving and parity non-conserving situations. In each case, we study the photon statistics of plasmon mode, the phonon statistics of the SiO 2 vibrations, and their crosscorrelation characterized by the cross-correlation function g (2) ab (0).

II. MODEL AND DYNAMICS
We consider a hybrid system consisting of an embedded qubit and an ENZ nanocavity with the ultrastrong plasmon-phonon coupling. The schematic diagram is sketched in Fig. 1, where the ENZ nanocavity is fabricated by the vibrational ultrastrong coupling of the plasmon mode and the SiO 2 phonon [35,44]. Especially we suppose that the plasmon modes are linearly coupled to a qubit [3,[34][35][36]. The Hamiltonian of the whole ultrastrong coupling qubit-plasmon-phonon hybrid system reads corresponds to the free Hamiltonian ( = k B = 1) of the photon, phonon, and qubit with ω 0 , ω T O , ω representing their corresponding frequency, respectively, FIG. 1: Schematic diagram of the coherently-driving qubitplasmon-phonon hybrid nanocavity. The driving strength is Ω, the driving frequency is ω d . The coupling of the plasmon mode and SiO2 phonon mode is ultrastrong. γσ − and γ b are the spontaneous emission rate of the qubit and phonon respectively, and γa is the decay rate of the cavity.
represents the photon-phonon interaction and the interaction between the photon and the qubit [35,[45][46][47] with denoting the coupling constants of the plasmon-phonon interaction, and a (a † ) and b (b † ) being the annihilation(creation) operators of the photons and phonons. Nevertheless, the term related to g D contains only photon operators due to the squared electromagnetic vector potential part of the light-matter interaction [46]. Meanwhile, g describes the strength of the coupling between the qubit and photon or phonon mode [48]. θ is a crucial physical quantity that has a significant impact on the spectra and the transitions of the qubit-plasmon-phonon system [14]. Note that there is no direct coupling between the qubit and phonons, which is quite different from the compared to the fully coupled hybrid cavity optomechanics [49]. The system is parity conserved for θ = π/2, which can be seen from that the parity of the system [22]. The parity is determined by the total excitation numbers. The even excitations are of odd parity corresponding to the eigenvalue −1 of Π, and the odd excitations are of even parity corresponding to +1. Obviously, the ground state is of even parity. The transitions between the states in the same parity space are forbidden because the transitions induced by the operators a−a † , b−b † , and σ−σ † will alter the parity of the state [30]. The energy levels versus the coupling strength are plotted in Fig. 2 for both the parity conserving (θ = π/2) and nonconverving (θ = π/4) cases. We are mainly interested in the coupling strength g = 0.25ω 0 which is explicitly marked in Fig. 2. It can be found from Fig. 2 that the cross points among the spectra occur with the coupling strength increase in the parity-conserving case. Therefore, the coupling strength g is crucial to the spectra of this hybrid system. It can affect the transition of the system and the parity of each state [3,30]. Compared with the parity non-conservation case, the energy levels of the same photon(phonon) state are more sparse when θ = π/2. Considering a coherent driving on the plasmon cavity mode as with the driving frequency ω d and driving strength Ω, the total Hamiltonian of this driven hybrid system is given by H total = H s + H d . Since the hybrid system is within the ultrastrong coupling and weakly driving regime, local treatments on each subsystem like individual dissipations couldn't be a good approximation anymore. They could lead to wrong results [50][51][52][53][54], which means that the local Lindblad master equation can not well describe the dynamics of the system. So we have to employ the global (nonlocal) master equation to study the dynamics [50,55]. Now we follow the standard process to derive the global master equation [55] and first turn to the H s representation. Based on the eigen-decomposition H s = E i |i i| with E i and |i corresponding to the eigen-energies and eigenstates, respectively, one can find H d in the H s representation as Thus the gobal master equation can be derived as [28,55] where is the dissipator, and Γ jk note the relaxation coefficients corresponding to the transition induced by c [3]. In addition, the parameters used in the paper show that all transition frequencies are different. Here we'd like to emphasize that all the transition operators are given by the dressed-state basis, which indicates the collective roles of the cavity, the atom and the phonon. This should be distinguished from the local master equation which is usually valid only for ω 0 , ω T O , ω ≫ g C , g D , g ∼ γ c , i.e., the weak intenal coupling.

III. PARITY-CONSERVING CASE
In order to reveal the photon/phonon statistics, we will consider the equal-time correlation function g (2) (0). In the traditional treatment of the field emitted by a cavity with the photon annihilation operator a and decay rate γ, its positive-frequency part is represented by the operator a out = a in + √ γa, where a in is the vacuum field impinging on the cavity mirror from the outside [56]. It means that a click of an external detector corresponds to annihilation of an intracavity photon. In the case of several interfering bosonic modes, the operator a should be replaced by a weighted sum of the boson annihilation operators of these modes [57]. In the case of ultrastrong coupling, though the operator a in the above expression should be replaced by the positive-frequency part of the intracavity fieldẊ + , having a more complicated form [3]. We make such replacements for the photon and phonon fields and consider the modified g (2) c (τ ) as [3,58,59] g (2) c whereẊ + c = j,k>j Y jk |j k| is the c related operator given in the H s representation with It is obvious that τ = 0 in Eq. (8) defines the equaltime correlation function g c (0). In addition, one can getẊ + c |0 = 0, which corresponds to the annihilation operator in the weak and strong coupling regime [16]. The intensity of the output photon or phonon flux emitted by a resonator can be represented by n c = Ẋ − cẊ + c [34,60].
To proceed, let' s first focus on the special case with θ = π/2, in which the parity of the system is conserved [3]. Since we are only interested in the weak driving limit Ω ≪ ω 0 , g, the multiple-photon (phonon) processes correspond to the high order of Ω, which allows us to consider our model by truncating different photon (phonon) numbers. For example, let's consider a simple case of N photons and N phonons. According to the spectra and the allowable transitions of the hybrid system in Fig. 2, one can find that the transitions between the ground state |0 and the excited states |j , j = 1, 2, 3 correspond to the single-photon process, and the transitions between the ground state |0 and the excited states with j = 4, 5, 6, 7, 8 correspond to the double-photon process, etc. Let the bare states |n a , n b , e(g) correspond to photonic, phononic modes and qubit, then the state |j can be spanned as |j = nan b ncC nan b nc |n a , n b , n c with n c = 1(0) denoting qubit's excited or ground state |e(g) . One can find that the state |j with j ∈ {k 2 , · · · , k 2 + 2k} for k ≤ N and j ∈ {f (k), excitations from the ground state |0 . We'd like to emphasize that the ground state |0 is not a vacuum state but a dressed state superposed by the bare states. Essentially, this is attributed to the counter-rotating wave terms as well as the squared vector potential which break the conservation of the excitation number in this ultrastrong coupling regime [12,27].
Due to the weakly driving with Ω ≪ ω 0 , ω, g, ∆ nm , we assume that the driving field does not change the eigenstate of the hybrid system, so the RWA can be safely taken over H d . Thus the total Hamiltonian in the H s representation reads where the second term corresponds to H d , K mn = Ω 2 m| a + a † |n and σ mn = |m n|.

Perform a unitary transformation on the Hamiltonian
with K = 0, 1, 2, · · · corresponding to m ∈ {K 2 , · · · , (K + 1) 2 − 1}, the time-dependence of the total Hamiltonian of the system will be eliminated [5] as Now we can employ a similar approach as Ref. [61] to deal with the dynamics with dissipations. When no quantum jumps occur, the evolution of the system is governed by the non-Hermitian effective Hamiltonian [62] In this sense, the dynamics of the qubit-plasmon-phonon hybrid system can be described by the Schrödinger equation i ∂Ψ(t) ∂t = H ef f |Ψ(t) . With the N 2 -state truncation, |Ψ(t) can be expanded as which leads to the following equations of the probability amplitudes: where k takes odd or even numbers for even or odd K ∈ [0, N − 1], respectively; Since the driving field is weak enough, the hybrid system will stay at the ground state with almost 1 probability, which can be explicitly illustrated by the populations of the lower 9 eigenstates plotted in Fig. 3, where all the parameters are selected within the current experimental conditions [35,63]. Note that the population |C j | 2 at ω d ≈ 1.05ω 0 is about two order of the magnitude less than the population |C 1 | 2 at ω d ≈ 0.74ω 0 in Fig. 3 (b). This is also the reason why the population |C 0 | 2 in Fig. 3 (a) has no apparent decreasing at ω d ≈ 1.05ω 0 . In order to give an explicit solution, we take N = 3 for example. Thus the above equations for steady states can be completely solved and one can obtain where with Substituting C n into Eq. (14), one can obtain the state |Ψ(+∞) . Thus the equal-time second-order correlation function of Eq. (8) can be directly calculated for photons as with the mean photon number n a given by Note that g (2) b (0) can also be obtained similarly, and for a precise solution, one can solve the equations with large N . But the solution is tedious, so we don't explicitly provide them here.
To provide an illustration of the statistical behaviors, we plot the equal-time second-order correlation functions g (2) a (0) and g (2) b (0) for photons and phonons, respectively, in Fig. 4. As a comparison, we plot the results by numerically solving the master equation (6) [64] and the analytic results by solving the similar equations as Eq. (15). Here we consider the state space with 5 photons and 5 phonons as an example. It is shown that our analytic method and the numerical simulation are in good agreement. One can find that there is a deep trough at ω d ≈ 0.74ω 0 in Fig. 4 (a) and (c), respectively, which indicates the evident photon/phonon antibunching. This corresponds to the resonantly driving transition from |0 to |1 via a single excitation process, which blocks the second excitation with a high probability and hence demonstrates antibunching in both figures so that the photons/phonons tend to appear in the nanocavity one by one. At ω d ≈ 1.3ω 0 , both photons and phonons approach the Poisson distribution. Meanwhile, these phenomena are accompanied by the relatively large intensity of the output photon/phonon flux, as shown in Fig. 4 (b) and (d). In Fig. 4 (a) there are 4 apparent peaks at ω d ≈ 0.8ω 0 , 0.94ω 0 , 1.05ω 0 , 1.13ω 0 corresponding to the double-photon resonantly driving transitions from |0 to |4 , · · · , respectively, which shows photon bunching tendency. Similarly, at the 4 frequencies, one can find from Fig. 4 (c) that phonons have similar bunching behaviors. The reason is that the bare basis spanning the 5 eigenstates |4 , · · · , |5 are symmetric on photons and phonons and especially lead to the comparable transition rates for photons and phonons. One significant difference is that phonons demonstrate the weakly-antibunching behavior around ω d ≈ 1.0ω 0 , accompanied by a slightly enhanced phonon flux. It can be understood that the driving field resonant with the frequency ω 0 of the cavity mode excites the photons in the cavity with a large probability; meanwhile, the strong optomechanical coupling induces a high probability of one-to-one photon-phonon conversion, which mainly corresponds to a single-phonon excitation. In this sense, the statistical properties of photons and phonons in the considered frequency region are roughly similar in the parity conserving case, but their details are different and even reverse around ω d ≈ 1.0ω 0 .

IV. PARITY NON-CONSERVING CASE.
In the previous section, we have employed an analytic approximation method to investigate the photon/phonon statistics in the parity conserving case. However, in the parity non-conserving case, the transitions will become a little bit complicated [14], which cannot be easily dealt with by a universal analytic treatment. So we will take a numerical process to study the correlation functions. In the following, one will find more striking physical effects with θ changing. In Fig. 5, we plot the secondorder correlation functions for both θ = π/4 and θ = π/2 and the intensities of the output flux (n c , c = a, b) for θ = π/4 versus the driving frequency ω d . Compared with the conserving parity, the obvious phenomenon in Fig. 5 is that the configurations of both the correlation functions and the intensities are squeezed from two sides to the central frequency ω d ≈ ω 0 , which further leads to quite different and even reverse statistical behaviors. For example, the most typical is that around ω d ≈ 0.74ω 0 (0.8ω 0 ), the strong antibunching (bunching) for both photons and phonons become the opposite statistical behaviors. The intensities of output flux are ab (0). The blue dashed line and the red solid line correspond to θ = π/2, and θ = π/4, respectively.
modestly increased in the middle range of frequencies and slightly decreased in the other range. This reason is that the non-conserving parity induces the eigenenergy translation to change the transition frequencies. For example, the transition |0 → |1 is resonant with ω d ≈ 0.74ω 0 in the parity conserving case, but resonant with ω d ≈ 0.8ω 0 if the parity is broken. In addition, within the range of relatively high frequencies, the bunching behaviors are greatly enhanced, which shows that the photons/ phonons are more inclined to be bunching in the parity non-conserving case, which is caused by the change of the transitions. For example, when we set an excited state in |3 via resonantly driving the |0 → |3 . The system is able to emit by |3 → |2 → |1 → |0 instead of only |3 → |0 as θ = π/2, which tends to cause photon bunching.
function g (2) ab (0) [3,65], which is defined as Similar to the usual g (2) (0), g ab (0) < 1 indicates that one photon and one phonon don't tend to exist simultaneously, on the contrary, g (2) ab (0) > 1 means that photon and phonon tend to bunch together. Roughly speaking, the cross statistical behaviors of photons and phonons, as shown in Fig. 6 are very similar to the phonons shown in Fig. 5 (c). The reason can also be basically understood as g (2) b (0), which won't be repeated. But photon and phonon for parity conserving case are inclined to much stronger antibunching at ω d ≈ 1.0ω 0 , which corresponds to the driving resonant with the cavity mode. As explained previously, the strong optomechanical coupling can induce the photon-phonon conversion, but the conserving parity greatly suppresses their simultaneous double excitation, which shows the antibunching behavior of photon and phonon. However, in the parity nonconserving case, all the bare basis |n a , n b , e(g) covered in the ground state |0 have nonzero probability amplitude. Therefore, photons excited by the resonant driving will trigger photon and phonon simultaneous excitation with a relatively large probability. Thus photon and phonon tend to the stronger bunching at ω d ≈ 1.0ω 0 . In addition, it is worth noting that, when ω d ≈ 0.74ω 0 (0.8ω 0 ) in the parity conserving (nonconserving) case, g ab (0) are all at the most obvious frequency dips of antibunching, which indicates that at ω d ≈ 0.74ω 0 photons and phonons in the nanocavity tend to appear one by one with large probability and disentangle with each other [66]. This result has important implications for the detection of phonons and photons [67].
We also study the delayed second-order correlation function g (2) a (τ ) in three different situations with the driving frequency resonating with the three lower excitations. As can be seen in Fig. 7, when we drive the transitions |0 → |2 and |0 → |3 , they both give rise to noteworthy oscillations of g (2) a (τ ). Obviously when we drive |0 → |2 , g (2) a (τ ) oscillates between bunching and antibunching at the frequency ∆ 21 . This oscillation is originated from that the driving on the transition |0 → |2 leads to the excitations of both the states |1 and |2 separated by ∆ 21 , which is similar to Ref. [3]. In particular, when driving |0 → |3 , three excited states can be excited. |3 and |1 states can be excited with relatively large probability, which leads to the oscillation of g (2) a (τ ) at the frequency ∆ 31 in a short interval, and then further triggers the excitation |2 , which leads to the oscillation at the frequency ∆ 21 . The photons are gradually inclined to be photon antibunching in this case. In addition, the oscillations also occur at the beginning in ω d = ∆ 10 with a much smaller amplitude.

V. CONCLUSIONS AND DISCUSSION
We have investigated photon and phonon statistics of a coherently driven qubit-plasmon-phonon hybrid system in the ultrastrong coupling regime. We have considered the parity-conserving and non-conserving regimes. Except for the regions where photons show strong bunching behaviors, but phonons tend to antibunching weakly, the statistical behaviors for both phonons and photons are pretty similar in the parity conserving case. The broken parity essentially leads to the translation of energy levels which squeezes the correlation function towards the central frequency. This squeezing triggers the reverse statistical behaviors in the low-frequency region and enhances the bunching properties in the medium and highfrequency regions. A similar phenomenon is also found in the photon-phonon statistics characterized by their crosscorrelation functions. The delayed second-order correlation function with different driving frequencies illustrates the striking oscillations, revealing simultaneous multiple excitations. Finally, we would like to mention that the strong photon-phonon coupling in Eq. (2) has been reported for a wafer-scale resonant ENZ nanocavity (cf. [35] and its references), and the strong coupling between a qubit and optical fields is also controllable by the bias flux of the qubit loop (cf. [3,34] and references therein). Their combination could still be a challenge. However, these strong couplings provide new insights into quantum nonlinear optical processes, especially the nonlinear optical effects of this nanocavity requires further addressing on the exotic implications.