Photoinduced polarization enhancement in biased bilayer graphene in the Landau level regime

We investigate the charge carrier dynamics in bilayer graphene subject to monochromatic laser irradiation within the Landau level quantization regime. Even though the radiation field does not lift the energy degeneracy of the lowest Landau levels (LLs), it nevertheless has a strong effect on the photoinduced pseudospin polarization response for higher LLs ($n\ge2$). Our results show that the photoinduced bandgaps lead to a finite response of the averaged pseudospin polarization with nontrivial oscillating behavior. It is shown that the contribution from these higher LL transitions turns out to be crucial to achieve an enhanced photoinduced polarization in radiated bilayer graphene. The experimental feasibility of our findings is also discussed.


I. INTRODUCTION
The proposals of driven phase transitions 1 in graphene [2][3][4][5] , so called Floquet topological insulators 6 , opened the road to explore many interesting phenomena in the family of two-dimensional materials that support Dirac fermions when they are subject to periodically modulated time-dependent interactions [7][8][9][10][11][12] . These photoinduced topological phases extended the static results of topological insulators to the dynamical regime showing that these topological phases could be dynamically generated even if the material showed a trivial topological phase in the static scenario. Indeed, within the static regime, both experimental and theoretical works have shown that the transport properties of topological insulating materials present very distinct properties that contrast those of conventional two dimensional electron gases 2DEG. An interesting example of those distinct features occurs in the Landau level structure of monolayer samples of graphene which, in contrast to the semiconductor 2DEG, shows a nonequidistant energy spectrum that, in turn, could allow the realization of a tunable laser in the Terahertz domain [13][14][15] . Moreover, the lowest LL in monolayer graphene can only be occupied in one of its sublattice degrees of freedom. This special behavior of the n = 0 Landau level (LL) in graphene monolayer renders the associated quantum Hall effect particularly interesting since at the charge neutral Dirac point it splits into four sublevels 16 at high magnetic fields (one for each valle K, K and one per spin state. When considering the bilayer graphene scenario [17][18][19][20][21][22][23] at low energies, the effective charge carriers behave as massive chiral particles in absence of a quantizing magnetic field. In addition, it has been experimentally shown 24 that it is possible to realize a gate induced insulating phase in bilayer graphene. Upon addition of a perpendicular magnetic field, the low energy excitations show an energy degeneracy at n = 0, 1 which is manifested in the transport measurements as additional van Hove singularities in the density of states 17 . Another experimental work shows that for the bilayer case in presence of a quantizing magnetic field, a 2π Berry phase is observed 25 which contrasts the π Berry phase acquired by Dirac fermions in monolayer graphene. One could expect that upon introduction of electromagnetic radiation, novel features should be feasible to be realized in the dynamical evolution of the effective charge carriers at low energies. Indeed these novel features have already been explored in irradiated bilayer samples with and without trigonal warping effects 17,[26][27][28] . However, the interplay of an applied bias voltage and circularly polarized monochromatic radiation on the pseudospin polarization of bilayer graphene and to what extend enhanced polarization inversion capabilities can be achieved is still a physics to be explored in bilayer graphene and is one of the main focus of the present work. The standard theoretical approach to describe the electronic properties of these materials rely on either first principles (DFT) numerical calculations or the use of a tight binding description, either of these gives detailed account of a number of the electronic features of both monolayer as well as bilayer samples 29 . However, it is well known that many phenomena of interest emerge already from the low energy physics 30 . In such regime, one can find analytically tractable models which, in turn, can shed light on the underlying physical processeses allowing for instance, the realization of novel transport features. In particular, within this low energy regime, most of the relevant physical features of bilayer graphene can be captured via an effective two-band model. In this work we theoretically analyze the dynamical manipulation of the LL structure of spinless charge carriers in biased bilayer graphene (in the Bernal stacking configuration), subject to a periodically driving radiation field applied perpendicular to the sample. In our approach we make use of Floquet's theorem [31][32][33][34][35][36] to recast the dynamics in an explicitly time-independent fashion but without the need to resort to the brute force numerical solution of the infinite-dimensional Fourier-mode expansion technique. Our approach has the advantage of providing an analytical description of the driven evolution of relevant physical quantities such as the pseudospin polarization. Since we are mostly interested in the leading order dynamical effects induced by the radiation field, we shall neglect trigonal warping effects that render the energy spectrum anisotropic at very low energies (< 1meV) 26 . We will also discard any spin-orbit effects (see, for instance, reference 23 for the interplay of spin-orbit effects and quantizing magnetic fields). In doing so, our analysis allows us to explicitly address each LL in an independent fashion and we show that the photoinduced bandgap depends on the Landau level index such that the n ≥ 2 LL quasienergy spectrum gives rise to a level-dependent bandgap. We then use this in order to infer the physical consequences in the dynamical evolution of physical quantities. Although bilayer graphene can also be experimentally realized in the AA stacking, where the two layers are laid on top of each other with the corresponding A2 (B2) atoms laying on top of the A1 (B1) atoms, this configuration shows a static unbiased energy spectrum consisting of two shifted copies of the monolayer graphene spectrum 37,38 , with shifted Dirac cones separated by an energy of the order of the interlayer coupling. Upon introduction of a biased potential U among the layers, a static bandgap develops. Therefore, although the interlayer coupling is larger than the typical values of the bias gate voltage, the dominant physical mechanism for static bandgap generation is the bias voltage term. Indeed, in the absence of gate voltage the Landau level physics under electromagnetic radiation would be that of two decoupled copies of monolayer graphene which we have already addressed in reference 39 . This is the reason why we are choosing the AB configuration in order to describe a distinct physical scenario as that of monolayer graphene. Thus, it is shown that in contrast to the single layer scenario, in the AB stacking configuration, the additional layer degree of freedom in bilayer graphene offers a richer physical structure for the pseudospin polarization. We find that its amplitude and decaying time can be enhanced via the radiation field, within experimentally accessible parameter regimes. In addition, the quasienergy spectrum LL anti-crossings emerging under radiation offer means to explicitly address an effective two-level system dynamics that is independent of the intensity of the quantizing magnetic field. We show that by properly tuning the laser parameters the photoinduced bandgap of the different Landau levels can lead to regimes from semi-conducting to metallic transitions with a finite to vanishing effective bandgap transition, respectively. We find that the photoinduced bandgap opening behavior enables a larger effective photoinduced polarization contribution for the higher order (n ≥ 2) LL states. As we will show below, our results on the photoinduced enhancement of the polarization effects could also lead to potential applications in quantum optics. In this realm, we could suggest using irradiated bilayer graphene in quantum optics for realizing a tunable laser taking advantage of the tunable effective bandgap that we have obtained. The paper is organized as follows. In section II we present the model and summarize the results for the quasienergy spectrum and the dynamics of the pseudospin polarization. Next, in section III we discuss the main results and we give concluding remarks, arguing on the possible experimental implementation of our proposed theoretical setup. Finally, in the appendix we summarized some mathematical calculations used during the perturbative analysis.

II. MODEL
Let us begin by considering the Landau level structure of biased bilayer graphene subject to intense circularly polarized monochromatic radiation. For definitness, we consider the planes of the layers to be the x − y plane, and the perpendicular direction to be that of the quantizing static magnetic field. We measure the momentum of the particles from the K point. Thus, the presence of the magnetic field modifies the momentum operator as p → p + e A ≡ π, where the vector potential satisfies ∇ × A = B and we are considering B = Bẑ. Our model for the bilayer graphene is such that the layers are arranged according to Bernal A 2 − B 1 stacking (i.e. atom A 2 from the upper layer lies directly on top of atom B 1 in the lower layer, see figure FIG. 1). In this case, near the K Dirac point at low energies and taking into account the time-dependent radiation field, we have a 4 × 4 Hamiltonian H(t) = H 0 + V(t), where H 0 is the static part and V(t) describes the light-matter interaction. The static contribution is explicitly given by ( = 1) where we have introduced the operatorsâ = (π x − iπ y )/ω c andâ † = (π x + iπ y )/ω c where the quantity ω c = √ 2v/ B describes the cyclotron frequency for massless Dirac fermions in monolayer graphene, γ is the interlayer coupling and U > 0 a biasing gate voltage strength. In addition, v ≈ 10 6 m/s is the Fermi velocity of charged particles in monolayer graphene, whereas −2 B = eB is the magnetic length defined in terms of the strength of the quantizing magnetic field B and e is the electric charge. Notice that the static Hamiltonian H 0 is written in the basis such that its upper left and lower right 2 × 2 sub-block matrices describe the lower and upper layers with inequivalent atoms labelled as A 1 , B 1 and A 2 , B 2 , respectively. In addition, the upper right and lower left 2 × 2 sub-block matrices describe the interlayer coupling of strength γ (between A 1 − B 2 sites). Here we are neglecting warping effects caused by weaker couplings γ 3 among sites A 1 − B 2 as well as a γ 4 coupling among A 1 − A 2 and B 1 − B 2 sites which give rise to electron-hole assymetry. Additional hoppings can be neglected in a zero order bilayer Hamiltonian. It is known that they can play a role at very low energies < 1 meV. The low energy physics is described via an effective two-band Hamiltonian to leading order in the interlayer coupling. In what follows, we shall concentrate in the range of energies |ε| γ, so in such case the reduced Hamiltonian takes the form, where Ω c = ω 2 c /γ = 2v 2 F eB/γ. In this manner, the largest energy scale is determined by the interlayer coupling parameter γ. Within this regime, the associated Landau spectrum in the absence of radiation reads E ns = sE n with E n = U 2 + n(n − 1)Ω 2 c , for n = 0, 1. The corresponding eigenstates are explicitly given by where n = 2, 3 . . . label the unperturbed Landau levels and s = ±1 label the conduction and valence band. We have also introduced the coefficients The degenerate LL corresponding to n = 0, 1 have energy E 0,1 = U and read We notice that in the presence of the bias voltage the electron density of the states can be changed and the associated bandgap can be tuned accordingly 27,28 . Yet, we will assume a small bias such that the electron density is essentially fixed. Thus, our results will be valid within the low light-matter coupling strength. Using this approximation, one can for instance evaluate the effective effective inter Landau level polarization polarization τ z , where τ z is a 2 × 2 Pauli matrix. Thus, one neads to evaluate τ z ns = φ ns |τ z |φ ns . After some algebraic manipulation we get which certainly implies that biasing the system does indeed introduces a finite pseudospin polarization via the energy bandgap and this manifest as a finite value of the effective inter Landau level polarization τ z ns .
We now deal with the photoinduced effects. In order to take into account the light-matter interaction in the model, we start from the standard minimal coupling interaction term −ev · A(t), which can be introduced, via the Peierls substitution, in the full 4 × 4 Hamiltonian. Hence, at the Dirac point, the effects of the driving field is described in the basis of ec. (1) by with 1 being the 2 × 2 unit matrix, σ = (σ x , σ y ) a vector of Pauli matrices and A(t) = A(cos ωt, sin ωt) the inplane associated radiation field which is related to the electic field via E(t) = −∂ t A(t), where A = E/ω, with E and ω = 2π/T being respectively, the amplitude and frequency of the radiation field, with T being its period. This periodic interaction makes the total Hamiltonian periodic in time H(t + T ) = H(t), with T = 2π/ω the period of oscillation of the driving field. Then, the timedependent contribution reads where we have introduced the effective light-matter coupling strength ξ = eEv F /ω, given in terms of E and ω the amplitude and frequency of the electric field, respectively. We assume that the beam radiation spot is large enough compared to the lattice spacing so we can ignore any spatial variation. On the other hand, we notice that the static Hamiltonian (1) commutes with the operator where the operator N a is defined as We can then perform a time-dependent unitary transformation given explicitly as which yields the time-independent Floquet Hamiltonian 31,32 where H F is given by We focus our analysis in the far infrared frequency domain 8 , where the laser energy is of the order of ω ≈ 10 meV and the effective radiation field intensity has the value E ∼ 1 kV/m, for which ξ ≈ 10µeV. Yet, we will show that our results could apply for larger electric field intensities E ∼ 0.15 MV/m for which one gets for the coupling constant ξ ≈ 10 meV. For frequencies ω in the Terahertz (ω = 3 THz) one gets ξ ≈ ω. This is an order of magnitude smaller than the Landau level separation ω c ≈ 116 meV, for B = 10 T, which is a typical experimental value at such fields. For larger frequencies and stronger magnetic field intensities, the ratio ξ/ω c tends to be smaller and our approximation scheme should provide values for the physical quantities that could be closer to those experimentally achievable. Thus, we can perform a perturbative treatment in the effective coupling parameter λ = ξ/ω c 1. We also notice that for intermediate values of the quantizing magnetic field satisfying ξ ω c , γ we can write down an effective two-band Hamiltonian among non-dimer sites; this is justified by recalling that quasienergies can be defined within the first Brillouin zone −ω/2 < ε < ω/2 and thus the effective two-band Hamiltonian approximation can be justified whenever ω γ (see discussion below). Let us now write down a perturbative Hamiltonian using the perturbation parameter of interest λ = ξ/ω c . For this purpose we use the antihermitian operator I − = a † σ − − aσ + , and build the 4 × 4 unitary matrix which transforms the Floquet Hamiltonian given in equation (13) as Since λ is small we restrict our analysis up to first order in the effective perturbation parameter λ. Introducing the shifted harmonic oscillator operators b = a + λ, we get (to leading order in λ) the effective Floquet Hamiltoniañ with H given by that takes into account corrections of order ξ. To get the result (16) we have neglected the off-diagonal contributions that can be treated by nondegenerate perturbation theory. They give corrections of order O(λ 2 γ 2 ) ≈ ξ 2 and O(λ 4 γ 2 ) ≈ ξ 4 , respectively. Then, they turn out to be less important than the last coupling term given in (17) that give corrections of order ξ. In addition, we have also neglected the higher order diagonal terms that could also be dealt by higher order perturbation theory. The contributions (18) and (19) can be relevant in the regime of small quantizing static magnetic fields B ≈ 1T , but can be discarded in our following discussion since, a posteriori shows that these higher order contributions do not qualitatively change our main results.
For bilayer graphene some works have reported values for γ ≈ 400meV (see 30 and references therein). Then, for a quantizing magnetic field used in experimental setups 8 B ∼ 10 T, such that Ω c ≈ 0.25γ we can safely use this effective low energy two-band approximation. Thus, for our purposes we can work within the effective two-band reduced Floquet Hamiltonian (16). First, we find useful to perform a unitary transformation H B =RHR −1 where R explictly reads that leads to, wheren b =b †b . The corresponding Floquet eigenstate is R|Φ = |Ψ , which has the two-component bi-spinor form |Ψ = (|ψ l |ψ h ) T , where the upperscript T denotes transpose and we have separated the lower energy |ψ l and |ψ h spinors corresponding to non-dimer and dimer coupling among the two layers. Thus, after eliminating the high energy spinor component we get the effective low energy two-band quasienergy problem H 2F |ψ = ε|ψ , where the effective two-band Floquet Hamiltonian reads now with Ω c = 2v 2 F eB/γ. The effective Hamiltonian (21) is valid whenever the condition γ Ω c , U , ξ is fulfilled. The Hamiltonian given in equation (21) has quasienergies ns = s [U − (n − 1)ξ] 2 + Ω 2 c n(n − 1) = s n , mod ω, with s = ±1, whereas the corresponding eigenstates read as where n = 1, 2, . . . label the shifted Landau levels. We also have defined the coefficients which, as expected, reduce to the unperturbed expressions as ξ → 0. The normalized quasienergies for n = 2 → 5 are given in FIG.2 for the interesting set of parameters U/γ = 2ξ/γ = 0.1. The upper left (right) panel shows the quasi-energy spectrum as a function of normalized cyclotron frequency for the n = 2 (n = 3) LL, whereas the lower left (right) panel corresponds to n = 4 (n = 5). In all panels, the dotted black line gives the static unbiased spectrum, the continuous black (red) curve corresponds to static (driven) biased regimes. Interestingly, we notice that the leveldependend bandgap is such that one can realize configurations where the driven regime is gapless (upper right panel), or the driven regime mimics the static biased scenario (lower right panel) which indeed shows the tunability of the photodinduced bandgap ∆ n . However, we emphasize that although the photoinduced bandgap might seem to lead to identical physical behaviour of the pseudospin polarization, we will show below that this actually not the case since the interference among the driven eigenstates mixes the static eigentstates with different weights. The latter, gives rise to a time-dependent term that is directly proportional to the driving strength. Notice however that, for higher LL the photoinduced bandgap continuously grows until the restriction ∆ n = ω is reached which is a consequence of the periodicity of the quasienergy spectrum.
Having dealt with the photoinduced bangdap spectrum we can study the radiation field effects on the layer-dependent pseudospin polarization τ z (t) = Ψ(t)|τ z |Ψ(t) , and we can interpret its fluctuations as an indirect measure of the angular momentum exchange among the graphene Dirac fermions and the electromagnetic field. That is, it provides information about the photodinduced dynamical hopping between the upper and lower layer of BLG. In order to gain further physical insight, we first show the effects of the radiation field by considering the initial state as an eigenstate (3) of the static effective two-band Hamiltonian H 2 |φ ns = E ns |φ ns , with n = 0, 1. After somewhat lengthy calculations (presented in the appendix) we get where ∆ n = U − (n − 1)ξ. As expected, in the limit ξ → 0, for which ∆ n → U and n → E n , one recovers the result (6). The time average of the polarization in one period τ z = 1/T T 0 τ z (t, ξ) dt gives with sincx = sin x/x. The pseudospin polarization is plotted in FIG.3 for the n = 2 → 5 LL states. We can observe that at low magnetic fields (Ω c → 0) the effective driven pseudospin polarization effects can be enhanced within the driven scenario when the bias voltage is present. As mentioned before, this can be explained by the interplay of the driving field and this bias which provides a LL-dependent bandgap favoring the interlayer hopping and hence the fluctuation in this physical quantity. Indeed, for any finite value of the quantizing magnetic field, the unbiased driven polarization shows a LL independent behavior which indicates that addressing the LLs requires the presence of the driving field for pseudospin inversion.
The most general scenario can be considered by writing the initial state as a linear superposition of the static eigenstates where the expansion coefficients satisfy the normalization condition |c 0 | 2 + |c 1 | 2 + s=±1 n≥2 |c ns | 2 = 1, and we have explicitly separated the n = 0, 1 eigenstates since they are degenerate in the pseudospin degree of freedom, as was discussed previously. Thus, the calculation of the pseudospin polarization gives now where the time-dependent contribution reads τ z (t) = s=±1 n≥2 1 2 n E n s|c ns | 2 ∆ n [U ∆ n + Ω 2 c n(n − 1)] + Ω 2 c n(n − 1) 2 ξ cos 2 n t + c * n,−s c ns Ω c n(n − 1) [Ω 2 c n(n − 1) + U ∆ n ] cos 2 n t − (n − 1)ξ∆ n − is n E n sin 2 n t .
As the n = 0, 1 solutions are still eigenstates of H 2F , there are not pseudospin oscillations in this subspace. Thus, this quantity cannot detect any relative phase among these two lower LL states but it does so for the higher LL. Then, we focus our attention on the oscillatory partτ z (t). For clarity of the analysis we have assumed for simplicity, but without loss of generality, that the expansion coefficients are independent of the pseudospin degree of freedom c ns = c n / √ 2. We notice that under this assumption, the first and last terms in equation (29) vanish which further simplifies the resulting expression for the polarization. Then the n ≥ 2 LL contributions to the pseudospin polarization now have the formτ Within this approximated scenario we notice that the selection of real expansion coefficients allow us a more transparent theoretical description of the pseudospin polarization effects. Clearly by selecting c 0 = c 1 = 0 implies that only pseudospin oscillations are considered. Upon averaging we get the result We have previously considered a coherent state superposition for which |c n | 2 = e −|α| 2 |α| 2n /n!, and show that interesting dynamical effects might arise in monolayer graphene LL 39 . Within the formulation of quantum optics, the parameter n = |α| 2 gives a measure of the average occupation of the coherent state. In order to determine the radiation field effects in a coherent state superposition, we select again as initial configuration a coherent state and in FIG.4 we show τ z as a function of Ω c /γ for different cases. The general outcome is that for both biased and unbiased regimes, the pseudospin polarization amplitude can be enhanced by means of the radiation field, and it takes typically longer for the driven pseudospin oscillation to decay. This shows that driving the system by monochromatic radiation affords a better control mechanism to address the pseudospin degree of freedom as compared to the role of the bias gate voltage allowing for longer and more pronounced pseudospin polarization effects.

III. DISCUSSION AND CONCLUDING REMARKS
We have analyzed the dynamical modulation of the Landau level structure of biased bilayer graphene subject to circularly polarized terahertz radiation. By means of a perturbative semi-analytical treatment we found that the most salient feature of the photoinduced modulation is to introduce a level dependent bandgap that provides an additional control parameter to modulate the electronic properties of the low energy particles present in the bilayer configuration. Nontrivial behavior of the pseudospin degree of freedom can be observed via the oscillations in the associated polarization dynamics. We show that the n ≤ 2 LL transitions are crucial to obtain a finite polarization for an initially prepared coherent state. In this manner, driving the system by monochromatic radiation could afford a better control mechanism to manipulate the pseudospin degree of freedom as compared to the sole application of the static bias gate voltage. In addition, we have also shown that in the driven scenario longer and more pronounced pseudospin polarization effects can be realized. Indeed, this is a novel feature of the driven scenario since the main new physical features are absent within the subspace spanned by the degenerate n = 0, 1 LL states. We would expect that the reported photoinduced gap modulation and pseudospin oscillations could be detected through the reemitted dipolar radiation from the oscillating charge carriers as it was proposed in reference 13 or with an appropriate modification of the experimental result reported recently 40 using a pump-probe femtosecond time-and angle-resolved photoemission spectroscopy (tr-ARPES) technique. In such experiment, a laser source is used to both map the energies of the excited states as well as follow the associated momentum-resolved population dynamics and serves to capture the transient population decay processes. The setup allow them to reach the femtosecond time scales associated to In/Si(111). We emphasize that the frequencies of interest in our model would lie in the near infrared region. Therefore, we could expect the detection scheme of our proposal to be in a much lower frequency scale associated to picosecond processes as has already been experimentally implemented recently in reference 41 . Since in this work we have considered pristine graphene samples, we would like to make a final remark on the role of defects in our results. Experimental evidence shows that defects might appear during the synthesis process of the sample which, at the nanoscale, might lead to interesting new phenomena since they could be exploited to generate novel, innovative and useful materials and devices 42 . For instance, these defects have been observed in situ via transmission electron microscopy 43 . The authors of reference 44 , have reported that point defects lead to notable paramagnetism but no magnetic ordering could be detected down to liquid helium temperatures, whereas the authors of reference 41 address the role of topological defects in photoinduced phase transitions. They show that long-range order is inhibited and is only restored when the defects annihilate. They also argue that their results would provide a framework for understanding other photoinduced phase transitions by identifying the generation of defects as a governing mechanism. Thus, considering defects in our setup should lead to further interesting results that would be addressed in future research. A very interesting review on the role of structural defects is given in reference 45 . In addition, we would expect that our results could also pave the road to discussing photoinduced Landau levels in graphene heterostructures with other materials such as black phosphorus 46-48 which would be addressed elsewhere. In summary, we have shown that photoinduced enhancement of the pseudospin polarization in AB bilayer graphene can be achieved within experimentally accessible parameter regimes. We expect that our results could lead to further interesting physical scenarios in other two-dimensional materials as black phosphorus or transition metal dichalcogenides, among others.
Acknowledgments-AL acknowledges useful discussions with Ernesto Medina and Leonardo Basile. This work has been supported by CEDIA via the project CEPRA-XII-2018-06 "Espectroscopía Mecánica: Transporte interacción materia radiación". AL thanks the University of Lorraine for partial financial support through research visits via CNRS-(PICS) project. AL and FM acknowledge funding from the project DGPA-PAPIIT IN111317.

IV. APPENDIX: DERIVATION OF THE PSEUDOSPIN POLARIZATION
Here we present some algebraic steps leading to the expressions of the photoinduced polarization effects. First, we give the results for an initially prepared eigenstate of the static effective two-band Hamiltonian where we have defined the coefficients whereas, the approximate Floquet states are given as with corresponding coefficients Then the polarization calculation is as follows = ss mm φ nτ |ψ m s ψ m s |τ z |ψ ms ψ ms |φ nτ e −(εms−ε m s )t , Using the definitions (33) and (35), we get τ z (t) nτ = τ ∆ n ε n E n U ∆ n ε n + Ω 2 c n(n − 1) E n + τ Ω 2 c n(n − 1) 2 ξ E n ε 2 n cos 2ε n t.
Upon substitution of the matrix elements φ n,τ |ψ n ,τ , and using equations (33) and (35) we arrive at the result given in equation (29). It is important to mention that in order to obtain the reported results we are assuming that, to leading order in the parameter λ = ξ/ω c , the relation m|n = δ nm among the original number operator eigenstates a † a |n = n |n and the shifted ones b † b |m = m |m is valid. Indeed, we expect that the scenario described in this work should hold in typical experimental setups in which the higher order corrections are negligible whenever ξ Ω c γ.