Matthieu Vanhoecke

Orazio Scarlatella

Marco Schirò

Dynamics Near a Photonic Band-Edge: Strong Coupling Effects Beyond Rotating-Wave Approximation

We study the dynamics of a quantum emitter coupled to a two-dimensional photonic crystal featuring a finite bandwidth with sharp edges and a Van-Hove singularity. We study the effect of strong system-bath coupling and non-Markovianity of the photonic environment using a nonperturbative approach based on the recently introduced NCA dynamical map for open quantum systems. We show that several characteristic features of the dynamics near a photonic band-edge such as the freezing of spontaneous emission and the maximum light-matter entanglement, get strongly modified in presence of counter-rotating terms in the system-bath coupling, beyond the rotating-wave approximation. Furthermore, by computing the spectral function of the quantum emitter we comment on the role played by atom-photon bound-state and show that this acquires a much larger lifetime once the rotating-wave approximation is relaxed.

I. INTRODUCTION

An atom coupled to its electromagnetic environment is one of the most fundamental example of open quantum system and the understanding of spontaneous emission has triggered fundamental advances in this field [1,[START_REF] Breuer | The Theory of Open Quantum Systems[END_REF]. Since then, the idea of controlling the electromagnetic environment of real or artificial atoms in order to modify and affect their dissipative dynamics has been a much explored and successful one. It has lead for example to the development of cavity [START_REF] Raimond | [END_REF] and circuit QED [4], where the atom is coupled to a resonant mode of a cavity or superconductor resonator. In more recent years the progress in controlling confinement of light in different types of nanostructures and platforms has made the design of photonic environment a practical experimental tool to modify and shape light-matter interactions between quantum emitters and electromagnetic modes, and to increase their coupling in the quantum regime. Examples include arrays of coupled cavities realizing photonic crystals [5], waveguide coupled to a Quantum Emitter (QE) in Waveguide-QED [START_REF] Sheremet | Waveguide quantum electrodynamics: collective radiance and photon-photon correlations[END_REF] or cold atoms arrays [START_REF] Goban | [END_REF][8][9][10].

Earlier investigations on QEs coupled to structured photonic environments such as photonic crystals displaying a photonic band gap revealed the emergence of atomradiation bound states, freezing of spontaneous emission and entanglement between matter and light [11][12][13][14][15]. Recently this topics has received renovated attention [16][17][18]. Much of these investigations, however, have focused on the weak light-matter coupling regime, where the Rotating-Wave Approximation (RWA), a conventional approximation in quantum optics, is valid. In the regime of strong or ultrastrong light-matter coupling [19,20] this approximation can break down and novel phenomena are expected to emerge [21][22][23][24][START_REF] Terradas-Briansó | Ultrastrong waveguide qed with giant atoms[END_REF]: in this work we focus on such regimes.

The theoretical description of small quantum systems strongly coupled to structured photonic environments in the non-perturbative light-matter coupling regime poses a number of challenges. One cannot make standard Markovian approximations leading to well known master equations, as these rely on the bath spectral function being smooth, which is not the case in presence of bandgaps with sharp edges or Van-Hove types of singularities. Also, one needs approaches able to capture the physics in non-perturbative regimes of the light-matter coupling. For these reasons, in this work we use a recently developed approach for open-quantum systems based on a selfconsistent dynamical map [START_REF] Scarlatella | Non-crossing dynamical maps for open quantum systems[END_REF]. This approach, based on the self-consistent resummation of infinite class of diagrams in the system-bath coupling corresponding to the so called Non-Crossing Approximation (NCA), allows to treat both the non-perturbative regime of light-matter coupling and to consider a bath density of states which is not necessarily smooth.

In this work we consider a model for a QE, described as a two-level system, coupled to a two-dimensional array of photonic cavities realizing a photonic crystal with finite bandwidth and sharp band-edges. Previous works have already studied this system within the RWA: when the emitter frequency lies outside the band [11] the dynamics is mainly dominated by bound states that are coherent states with a long relaxation time. When the emitter frequency lies inside the band the singularity in the center of the band leads to a non-perturbative regime where the Fermi golden rule fails [START_REF] González-Tudela | [END_REF]28]. We investigate the system beyond the RWA, computing its dynamics with the NCA dynamical map both without and within the RWA. We consider the cases of an emitter frequency well within the photonic band, outside of the band and at the band edge and we show that in the latter case deviations from the RWA are most pronounced. We show that several characteristic features of the dynamics near a photonic band-edge, such as the slow-down of the spontaneous decay or the light-matter entanglement are strongly modified in presence of counter-rotating terms. Furthermore, by computing the spectral function of the QE we comment on the role played by atom-photon bound-state and show that this acquires a much larger lifetime once the RWA is relaxed. This paper is organized as follows. In Sec. II we introduce our model for a quantum emitter coupled to a photonic bath, write down the Hamiltonian for the full lightmatter coupling and introduce the RWA to which we will compare our results. In Sec. III we discuss the method we use to solve for the dynamics of the QE, namely the NCA dynamical map. In Sec. IV we present our results for the dynamics of spontaneous emission, entanglement entropy and for the emitter spectral function. In Sec. V, we discuss a qualitative picture to understand the results obtained in the previous section. Finally, In Sec. VI we draw our conclusions. Appendixes contain more details about technical aspects of our work.

II. THE MODEL

We consider a two-level system (TLS) {|↑ , |↓ } coupled to a photonic bath, describing a Quantum Emitter (QE) or a defect embedded in a photonic crystal. The full Hamiltonian H of the system can be written as

H = H S + H B + H SB
where H S = ω q σ z is the TLS Hamiltonian with transition frequency ω q , H SB is the light-matter coupling while H B describes the photonic bath that we take to be a twodimensional photonic lattice described by a tight-binding model with nearest neighbor coupling J and on-site energy ω * . The Hamiltonian for the bath is given by

H B = r ω * a † r a r -J <r,r > a † r a r + H.c ,
where r = (x, y) denotes the position of the photonic mode in the lattice. By going to Fourier space we can diagonalize the Hamiltonian and bring it to the following form

H B = k ωk a † k a k (1)
where we have introduced the dispersion ωk = ω * + ω k = ω * -2J [cos (k x ) + cos (k y )] and the Fourier component of the photon operators

a k = 1 √ N r e -ik•r a r with k x,y = π N (-N, -N + 1, • • • , N -1)
. The main interest of this bath consists in its density of states (DoS)

D(ω) = 1 (2π) 2 dkδ (ω -ωk ) (2)
illustrated in the Fig. 1 when N → ∞. First, the bath has a finite bandwidth, i.e. its frequencies extend between ωk ∈ [ω * -W, ω * + W ] where we define W = 4J and are centered around ω * . This allows us to discuss different situations depending on the value of the TLS frequency ω q , as we discuss below. Furthermore, the DoS has strong non-Markovian features, including sharp band-edges and a Van-Hove singularity in the middle of the band.

The QE is coupled to the photonic bath locally at a given site through the dipole operator of the TLS which is linearly coupled to the electric field, giving rise to a dipole gauge type of Hamiltonian which reads

H SB = σ x k g k a k + a † k (3)
In the following, we consider only one mode of polarization of light and a k-independent coupling constant g k = g. Relaxing these constraints also gives rise to interesting phenomena [29]. Combining all the terms together we obtain an Hamiltonian for the light-matter system of the form

H = ω q σ z + k ωk a † k a k + σ x k g k a k + a † k (4)
In the following we will be interested in comparing the dynamics generated by Eq. ( 4) with the one obtained under the RWA. This is a widely used approximation in quantum optics which is valid when the bath and QE are close to resonance and weakly coupled: in this case, their dynamics is dominated by their bare frequencies and thus the "counter-rotating" terms in the light-matter coupling are rapidly oscillating and can be neglected. This yields

H RW A = ω q σ z + k g k a k σ + + a † k σ -+ k ωk a † k a k (5)
In this form the Hamiltonian conserves the total number of excitations N exc = σ z + k a † k a k and allows for a simple solution in the subspace at fixed N exc . In particular, for N exc = 1 several results have been known in the literature. In the single excitation sector, the total state of the system described by H RW A can be written explicitly, if we suppose a empty bath initially |0 .

|ψ(t) = C eg (t)σ + + k C k (t)a † k |↓ ⊗ |0 (6)
using this expression and the propagator we can diagonalize the hamiltonian H RW A and also obtain the dynamics of spontaneous emission [START_REF] González-Tudela | [END_REF]. Note that when the band admits discontinuities or singularities, it is necessary to carry out analytical continuations or approximations in order to obtain the dynamics. Here we will discuss the effects of counter-rotating terms on the physics of this model, as a function of the light-matter coupling g and the frequency of the emitter ω q that we parametrize in terms of ∆ = ω q -ω * , the detuning with respect to the middle of the band ω * (See Fig. 1). This formally corresponds to go to a frame rotating at the central frequency of the photonic band ω * , by applying the unitary transformation Û = exp itω * σ z + k a † k a k . We note that in the limit ω * , ω q → ∞ we expect the effect of counter-rotating terms to disappear and the RWA to become exact.

III. NCA DYNAMICAL MAP

In this section we briefly discuss the theoretical approach we use to compute the dynamics of our model. Since we are interested to go beyond the weak-coupling regime in which RWA is valid, we resort to the recently developed Non-Crossing dynamical map [START_REF] Scarlatella | Non-crossing dynamical maps for open quantum systems[END_REF]. This approach is similar in spirit to a master equation, describing the reduced dynamics of the QE after the bath is integrated out, but it yields an equation for the dynamical map of the system, that is the superoperator V(t) evolving the reduced density matrix, rather than for the density matrix itself. The dynamical map is defined by ρ(t) = T r B ρ tot (t) = V(t)ρ(0); expanding in the systembath coupling, one finds that it obeys the Dyson equation

∂ t V(t)• = -i H S , V(t)• + t 0 Σ(t -t 1 ) V(t 1 )• (7)
where Σ is defined as an infinite series involving the 2times correlations functions of the bath and where the bullets indicate the arguments of super-operators, when necessary. Within a non-crossing approximation, described in Fig. 2, and without making the RWA, the self-energy takes the simple analytical expression

Σ NCA (τ ) = Γ(τ ) V(τ ) [σ x •] σ x -σ x V(τ )σ x • + h.c. (8) Γ(τ ) is the 2-times correlation function of the bath operator Γ(τ ) = Tr B [B(τ )B(0)ρ B (0)] with B = k g k a k + a † k and a ( †)
k (τ ) the creation (annihilation) bosonic operator in the interaction picture. We note that

መ 𝑆 𝑁𝐶𝐴 = + + +… = FIG. 2.
Self-energy in the non-crossing approximation (NCA). Solid lines correspond to bare time evolution superoperators (without the bath), dashed lines correspond to bath 2-times correlation functions, bold solid lines correspond to the dynamical map V. Only diagrams where the dashed lines do not cross are kept in the non-crossing approximation.

the self-energy has a similar structure to the "dissipator" of standard master equations, and in fact the Born Master equation can be recovered by replacing V(τ ) by e H S τ in the expression (8) of the self energy. The expression for the self-energy within the RWA is reported in the Appendix A.

IV. RESULTS

In this section, we discuss our results for the dynamics of the system obtained with the NCA dynamical map. We consider the emitter initially in the excited state ρ (0) = |↑ ↑| coupled to an initially empty bath, described by ρ B (0) = |0 0|. Then we let the system evolve through the Hamiltonian with (5) and without RWA (4), and study the emitter relaxation.

We fix ω * W and we voluntarily choose a large bandwidth, in order to minimize the influence of the edge (singularity in the middle) of the band when the emitter frequency lies in the middle of the band (edge).

We start discussing the dynamics of spontaneous emission, described by the time-dependent population of the TLS

P (t) ≡ 1 2 (1 + σ z (t)) (9)
In Fig. 3, we show the results obtained numerically for the dynamics of the emitter, for different values of detuning between the emitter and the center of the photonic band ∆/W = -1.2, -1.0, -0.25, 0, 0.5 and for a fixed coupling constant g/W = 0.05. Solid and dotted lines correspond respectively to the dynamics without and within the RWA. As expected, we find that the decay of the spontaneous emission is faster when the frequency of the emitter lies within the photonic band, |∆/W | < 1, while it slows down as its frequency is moved outside of the band since the emitter cannot hybridize with the bath. The decay is exponential in time and can be estimated using perturbation theory and Fermi Golden's Rule (FGR) to give Γ F GR (∆) = 2πg 2 D(∆). This is true except if the emitter frequency is close enough to the middle of the band, where the density of states features a Van-Hove singularity: in this case the decay is non-exponential and perturbation theory breaks down as discussed in the literature [START_REF] González-Tudela | [END_REF]28,30].

We see from Fig. 3 that for (almost) all detunings, including when the emitter is resonant with the singularity ∆/W = 0, the dynamics within and without the RWA agrees for the (weak) light-matter coupling considered. We also remark (not shown) that the NCA dynamical map is able to predict a deviation from exponential relaxation when the emitter frequency is sufficiently close to the singularity. Particularly interesting for this work is the regime in which the emitter frequency is at the lower edge of the photonic band, corresponding to ∆/W = -1 in Fig. 3: in this case we see that significant deviations arise depending on whether one makes the RWA or not. This observation suggests that going beyond the RWA is necessary to describe the physics in this regime, due to the sharp discontinuity of the bath DoS at the edge of the band.

In Fig. 4 we focus on the regime in which the emitter frequency is at the edge of the photonic band ∆/W = -1, where we expect the deviations from RWA to be most important, and vary the light-matter coupling g from weak to strong. The top panel corresponds to the evolution obtained within the RWA and the bottom panel to the dynamics generated by the full Hamiltonian.

In this latter case we see that upon increasing g the exponential relaxation of the spontaneous decay observed before is modified. In particular we observe the expected freezing of spontaneous emission which appears in the form of intermediate time plateau in the time evolution. For g/W = 0.05 for example we clearly see a rapid decay of P (t) followed by an almost constant evolution up to some longer time scales at which the dynamics relaxes again towards zero. As the coupling is further increased the metastable plateau is reduced, however we see that at stronger couplings multiple plateau re-emerge. This is a phenomenology that is strongly reminiscent of prethermalization in weakly non-integrable many-body systems which possess several well separated energy scales corresponding to the unlocking of almost conserved quantities. This phenomenology is modified by the presence of counter-rotating terms (top panel) which give rises to important differences in the dynamics of the spontaneous emission. In particular we note a much faster decay towards the steady state, which quite importantly does not need to be the ground-state of the Qubit since counterrotating terms can support a non trivial steady state. Quite interestingly, we see that the plateau structure reminiscent of metastability is not completely washed away, as we see at intermediate and strong coupling g/W = 0.1 -0.2.

Within the RWA the freezing of spontaneous emission is usually understood with the appearance of coherent state in the dynamics, which live outside the band (bound state) [13,21,23,31]. At the lower edge of the band, the frequency of the emitter splits between coherent states in the gap and modes in the band [11], is it in this that we find a dynamic which for a long time relaxes towards the ground state of the QE. Here, for the full dynamics this picture will need to be modified as we will discuss in Sec. V.

B. Dynamics of Entanglement Entropy

To better understand the dynamics at strong coupling and the role of counter-rotating terms we now discuss the evolution of the entanglement entropy between the QE and its photonic environment. We emphasize that the full system plus environment evolve unitarily and therefore we can consider a bipartition containing the emitter and define its reduced density matrix as ρ(t) = Tr B ρ tot (t), where the trace is taken over the photonic bath degrees of freedom. Then the entanglement entropy between emitter and its environment reads

S(t) = -Tr (ρ(t) log ρ(t)) . ( 10
)
We note incidentally that the emitter reduced density matrix ρ(t) is a quantity we can naturally access with the NCA dynamical map. For spin-boson types of models the study of entanglement entropies in equilibrium and at the quantum phase transition has received major attention in the past [32][33][34].

We first consider the case in which the QE lies within the photonic band, ∆/W = -0.75, and different values of system-bath coupling g/W and compare in Fig. (5) the full dynamics obtained within NCA with the RWA. We see that in both cases the entanglement grows with time, starting from zero since the initial state we consider is a product state between system and bath, reaches a maximum value and then decreases. The value of this maximum entanglement does not depend on the coupling g and concides with the maximum entanglement for a twolevel system, i.e. the system is maximally entangled at short times due to the coupling to the photonic bath. On the other hand, the time scales to reach the maximum entanglement and then to leave it towards the stationary value depend strongly on g: we see that quite generically the dynamics is much faster as g is increased. The comparison between full dynamics and RWA reveals that the dynamics of entanglement is similar, while the long-time limit is different. In particular the counter-rotating terms are able to sustain a finite entanglement entropy at long times, while within RWA the entanglement always goes to zero at long times since the QE relaxes back to the ground state.

We now discuss the more relevant case for us of a QE sitting at the edge of the photonic band, i.e. ∆/W = -1. In Fig. 6 we plot the resulting dynamics of the entanglement entropy for the full Hamiltonian (top panel) and the RWA case (bottom panel). Focusing on the latter case first we note that the main difference that appears with respect to Fig. 5 is the fact that the maximum entanglement is stabilized and mantained for a finite timeinterval, before decaying to zero. This maximally entangled time interval depends on the system-bath coupling, i.e. for weaker coupling the plateau is longer-lived while it is destroyed rapidly for strong g. This phenomenon, which corresponds to a back action of the reservoir between the emitter and the mode of the photonic band [12,[35][36][37] is linked to bounds state in the dynamics. It can be understood as the entanglement entropy counterpart of the prethermalization plateau observed in the dynamics of spontaneous emission discussed before. In particular we can use the same kinematic argument and the conservation of total excitation number to understand the resilience of entanglement to decay in the case of RWA dynamics. This picture is further confirmed by looking at the full dynamics beyond RWA, shown in the top panel. There we see that the maximum entanglement is preserved at short times: the rapid increase of entanglement entropy seems not affected by the counterrotating terms which instead play a more relevant role at longer times, destabilizing the maximally entangled state on much shorter time scales and leading to a stationary state which is not the trivial ground-state of the isolated emitter but contains non-zero system-bath entanglement.

C. Photonic Bound-State and Quantum Emitter Green's Function

A typical approach to understanding the link between dynamics and the nature of excitations is to look at the spectral function of the system encoded in the retarded Green's function of the emitter, which is defined as

χ(t) = -iθ(t) [σ x (t), σ x (0)] . (11)
Here, the average is taken over the stationary density matrix of the coupled system and bath and the timeevolution is performed with respect to the full Hamiltonian of system and bath. This quantity can be computed within our NCA dynamical map, by expressing Eq. ( 11) in terms of the evolution superoperator V(t) defined in Sec. III [START_REF] Scarlatella | Non-crossing dynamical maps for open quantum systems[END_REF]. From this quantity, we can extract the Fourier transform and take the imaginary part which contains information about the spectrum of the system, i.e. χ(ω) = Im dte -iωt χ(t).

We start considering the case in which the emitter frequency lies within the photonic band, i.e. ∆/W = -0.5 and plot in Fig. 7 the spectral function for different values of light-matter coupling g, both in the case of RWA (right panel) and for the full Hamiltonian (left panel). In both cases we see some common features emerging in the spectral function, including a peak within the band (dotted lines represent the upper and lower edges FIG. 8. Spectral function of the QE, as defined in the text, for differents values of coupling g/W and for a QE frequency lying at the edge of the photonic band, ∆/W = -1. We plot only positive frequencies and represent with the dotted lines the upper and lower edges of the photonic bath. Left panel shows the result for the full Hamiltonian, while the right panel within the RWA. We see that the deviations appear already for moderately weak couplings, where the RWA results show a narrow resonance within the gap, corresponding to a photonic bound-state, while the inclusion of counter-rotating terms leads to a broadening of the in-gap peak and a transfer of spectral weight within the band.

of the band) corresponding to the frequency of the emitter ω q /W = 0.625 which shifts and becomes broader as the coupling g is increased. At stronger values of the coupling we see also the appearance (in both panels) of a spectral feature in the middle of the band corresponding to the Van-Hove singularity. On the other hand, we see that the behavior of the spectral function at frequencies near the band-edge is rather different and that the presence of counter-rotating terms in the full Hamiltonian has direct consequences on the spectral features of the system. In particular as we increase the coupling g we see that states in between the gap appears in the spectrum, corresponding to processes involving virtual photons that cannot appear within RWA due to the conservation of total number of excitations. Furthermore, we see the emergence of a sharp peak at the edge of the band which is completely absent in the RWA data and which gets stronger as g increases. We will comment on the origin of this spectral feature in Sec. V.

We now discuss the spectral function for a case in which the emitter frequency lies at the edge of the band, which for the dynamics corresponds to the stronger differences between RWA and full dynamics. In Fig. 8 we plot the spectral function for ∆/W = -1 and different values of the coupling g respectively for the full Hamiltonian (left panel) and the RWA case (right panel). We see that in this case major differences appear in the spectrum already for relatively small values of light-matter coupling. In particular within the RWA there is a narrow peak within the photonic band gap whose position depends on g and which corresponds to a photonic bound-state.

Once counter-rotating terms are included they lead to two major effects, namely a broadening of the in-gap peak which becomes a resonance and acquires a finite-lifetime and a transfer of spectral weight into the photonic bath. We also notice the presence of a features right at the band-edge. In both cases, we observe a splitting of the emitter excitation into a coherent excitations outside the band and others modes living in the band. This splitting corresponds to the freezing of the spontaneous emission, indeed the band modes have a faster decay time than the coherent excitations outside the band. Moreover, this splitting is all the more important and brings into play all the more modes that the coupling with the bath is strong.

V. DISCUSSION

The results shown in the previous section highlight the fact that counter-rotating terms in the system-bath Hamiltonian become relevant for the dynamics of the system not only at strong coupling g, but also at weak and intermediate values provided the frequency of the quantum emitter is resonant with the sharp band edge. In order to qualitative understand these results it is useful to perform a time-dependent unitary transformation and to rewrite the full Hamiltonian of H in the interaction picture with respect to the free evolution of emitter and bath, H 0 = ω q σ z + k ωk a † k a k . This gives

H(t) = g k e i(ωq-ω k )t σ + a k + h.c. + + g k e i(ωq+ω k )t σ + a † k + h.c. ( 12
)
where the system-bath coupling is now explicitly timedependent and contains two types of terms, those conserving the total number of excitations and oscillating at frequency ω q -ωk and the counter-rotating terms oscillating at frequency ω q + ωk . The RWA amounts to disregard the latter terms which are rapidly oscillating as compared to the number-conserving couplings. This approximation is usually valid provided that

|ω q + ωk | |ω q -ωk |
This is true in particular whenever the emitter frequency is resonant with the photonic band, ω q = ωk or equivalently ∆ = ω k , which is the case (for some value of k) whenever |∆|/W < 1 (see Fig. 9). We expect therefore the RWA to work well whenever the emitter is resonant with many modes of the photonic bath. When this is the case the Hamiltonian in Eq. ( 12) can be split into a sub-set of k modes which are resonant with the emitter and static, their counterrotating terms oscillating at frequency 2ω q acting as a drive and a set of complementary bath modes which are off-resonant and oscillate with multiple frequencies and include both rotating and counter-rotating terms. The condition ∆ = ω k is met by an extensive number of k = 0 points if, for example, ∆ = 0 in the middle of the band, so the subset of coherently coupled resonant modes act itself as a bath and provides dissipation. In this regime one expect a fast dynamics of the emitter and a good agreement between RWA and the full dynamics, except on long time scales where the counter-rotating terms can stabilize a non-trivial stationary state. However, as ∆ approaches the lower edge of the photonic band the number of resonant k modes decreases. Right for ∆ = -W , corresponding to the emitter at the edge of the photonic band, this condition is only met at k x = k y = 0. This suggests an effective Jaynes-Cumming or Rabi type of model in which the emitter is coherently coupled to the resonant k = 0 mode, with time-dependent coherent drive of counterrotating terms and off-resonantly coupled to the rest of the bath modes which provide dissipation.

H eff (t) = g σ + a + h.c. + g σ + a † e 2iωqt + h.c. + + g k e i(ωq-ω k )t σ + a k + h.c. + + g k e i(ωq+ω k )t σ + a † k + h.c. ( 13
)
Furthermore all the couplings in the RWA have to conserve the total number of excitations. In this picture the freezing of spontaneous emission and the entanglement plateaux arise because of the interplay between off resonant system-bath couplings, which are responsible for bringing the system back to its ground-state but are kinematically blocked at small g because of the conservation of excitation number, and resonant coherent coupling with the k = 0 mode which instead protects entanglement and a finite value of the emitter polarization. On the other hand when the counter-rotating terms are included the picture changes substantially: from one side new coherent coupling terms arise between the emitter and the k = 0 mode and new off resonant dissipative channels emerge, both of which break the conservation of excitations and lead to a non-trivial stationary state but also spoils the separation of energy scales associated to the metastable states observed in the RWA.

The argument above can also be used to understand qualitatively the origin of the strong spectral response of the system at frequency corresponding to the lower edge of the band, as shown in Fig. 7. In fact we note that in the rotating frame picture discussed above the full Hamiltonian can be obtained from the RWA one by including a time-dependent perturbation given by the counter-rotating terms, which oscillate at ω q + ωk . In the high-frequency and weak coupling limit, corresponding to the regime of validity of RWA as discussed before , we expect this perturbation to be irrelevant for the physics of the system. However we note that the modes at the edge of the band correspond to the frequencies k = 0 and are always the slowest ones ω q + ωk . For this reason we expect that in the response to a weak drive at that frequency, which corresponds essentially to the spectral function discussed in Sec. IV, the deviations between RWA and full dynamics to be more visible at frequencies near the lower band edge.

Finally, we note that while the validity of RWA is usually assumed for resonant couplings, we can provide an extended criterion of validity of RWA, which reads, in the case the initial density matrix contains a single excitation (See Appendix B),

g 2 max ωk 1 ,ω k 2 1 (ω k1 + ωk2 ) (ω q + ωk1 ) 1 ( 14
)
and that gives the condition

g 2 (2ω * -8J) (2ω * + ∆ -4J) (15)
explaining why the RWA fails more strongly when the emitter is at the lowest edge of the photonic band.

VI. CONCLUSION

In this work we have studied the dynamics of a quantum emitter strongly coupled to a photonic environment featuring a finite band-widht and a band gap. We have discussed in particular the role of the counter-rotating terms in the full Hamiltonian and the consequences they bring with respect to the dynamics within the Rotating Wave Approximation (RWA). To solve the resulting spinboson model we have used the recently developed NCA dynamical map, which gives direct access to the reduced density matrix of the emitter, from which local properties as well as spontaneous emission, entanglement entropy can be readily obtained as well as the frequency resolved spectral function. We have shown that the major deviations from RWA due to the counter-rotating terms arise when the emitter frequency is at the lower edge of the photonic band where the freezing of spontaneous emission and the maximum entanglement, usually interpretred in terms of photonic bound-states, are destroyed and the system is driven towards a non-trivial and entangled stationary state. We have shown that spectral features corresponding to this physics arise in the emitter spectral function. Specifically we have shown that the narrow in-gap peak appearing in the RWA case becomes a much broader resonance due to counter-rotating terms and also hybridizes with modes within the photonic band. We have provided a qualitative picture to understand these results in terms of an effective Rabi-like model, featuring a resonant coherent coupling between emitter and the k = 0 mode and to a bath of off-resonant excitations which provide dissipation. Our results show the importance of counter-rotating terms in the dynamics of wave-guide QED systems not only at ultrastrong coupling but also at intermediate light-matter coupling regimes, depending on the relative detuning of the emitter. Directions for future work could include for example the interplay of two or more quantum emitters, which can be still treated within the NCA dynamical map, or the calculation of photonic bath properties to inquire the detail structure of the system-bath wave function at strong coupling.

the order 1 in g does not contribute because initially we have the density matrix ρ tot (0) = |↑ ⊗ |0 . In the same way for the second order, we have 5 possible contributions for the dynamics of the Qubit:

CR (τ ) (B7)
where in its expressions, we denote H

(n) CR = O (g n ),
that is, we perform an expansion of the counter rotating term into the coupling constant g. Due to the initial state corresponding to a single excitation for the Qubit, one can easily convince oneself that the contributions (1-3) do not participate in the dynamics, because we have H e it(ω k 1 +ω k 2 ) -1 (ω k1 + ωk2 ) (ω k2 -ω q ) -e it(ω k 1 +ωq) -1 (ω k1 + ω q ) (ω k2 -ω q ) a † k1 a † k2 σ + σ -ρ tot (0) (B8)

So, the Rotating wave approximation is then valid when the following criterion is met:

g 2 max ωk 1 ,ω k 2 1 (ω k1 + ωk2 ) (ω q + ωk1 ) 1 (B9)
using the definition of detuning ∆ = ω q -ω * and remembering that we choose ω * such that ωk > 0, for all k, the expression can be reduced to g 2 (2ω * -8J) (∆ + 2ω * -4J) (B10)

FIG. 1 .

 1 FIG.1. Density of states for the photonic bath realised by a two-dimensional lattice of coupled cavities with nearestneighbor hopping J, as described in the main text. The spectrum is centered around the frequency ω * and has a finite half-bandwidth W = 4J. We note the Van Hove singularity at ω * and the sharp band edges at ω = ω * ±W . The quantum emitter frequency is ωq and we parametrize it in terms of its distance from the central frequency ω * , through the parameter ∆ = ωq -ω * . In this respect, an emitter resonant with the bath has |∆|/W < 1, while for ∆ = ±W the emitter is at the edge of the band.

FIG. 3 .

 3 FIG.3. Spontaneous emission of a single quantum emitter, initially prepared in its excited state, and then suddenly coupled to a photonic bath, with a coupling constant g/W = 0.05 and for different detuning values ∆/W = -1.2, -1, -0.25, 0, 0.5, corresponding to a QE outiside the band (∆/W < -1) at the edge (∆/W = -1) or within the band (∆/W > -1). Solid lines correspond to the dynamics of H and that dotted at HRW A.

FIG. 4 .

 4 FIG. 4. Spontaneous emission of a single quantum emitter with frequency at the lower edge of the photonic band, corresponding to ∆/W = -1, initially prepared in its excited state and then suddenly coupled to the bath. The top panel correspond to the dynamic of the full Hamiltonian while the bottom panel is the evolution of the entanglement entropy for the Hamiltonian under RWA.

FIG. 5 .

 5 FIG. 5. Dynamics of the entanglement entropy for a QE with detuning within the photon band, ∆/W = -0.75 and for different values of coupling g/W with the bath. Solid lines correspond to the dynamics of H, and the dotted one the dynamics of HRW A.

FIG. 6 .

 6 FIG. 6. Dynamics of the entanglement entropy for a QE with detuning at the edge of the photonic band, ∆/W = -1 and for different values of coupling g/W with the bath. The top panel correspond to the dynamic of the full Hamiltonian while the bottom panel is the evolution of the entanglement entropy for the Hamiltonian under RWA.

FIG. 7 .

 7 FIG. 7. Spectral function of the QE, as defined in the text, for differents values of coupling g/W and for a QE frequency lying within the photonic band, ∆/W = -0.5, corresponding to a frequency ωq/W = 0.625 We plot only positive frequencies and represent with the dotted lines the upper and lower edges of the photonic bath. Left panel shows the result for the full Hamiltonian, while the right panel within the RWA. For weak couplings we see in both cases a a narrow resonance centered around ω = ωq and a small feature in correspondence of the singularity in the middle of the band. For stronger couplings the deviations between full Hamiltonian and RWA become more evident, both in terms of width of the resonance and behavior at the edge of the band.

FIG. 9 .

 9 FIG.9. Graphical representation of the modes resonant with the frequency of the Quantum Emitter, ω k = ∆. We notice a continuum of resonant modes when the quantum emitter frequency is located in the band, |∆|/W < 1, which shrink to a point as ∆ approaches the lower edge of the photonic band, ∆ = -W . On the other hand, at the upper edge ∆ = W the resonant condition is satisfied by four state.

  CR = e i(H S +H B )τ H CR e -i(H S +H B )τ .With regard to contribution 4) (in the same way for the contribution 5) we can rewrite it in the following formCR (τ )ρ tot (0) = g 2 k1,k2

ACKNOWLEDGMENTS

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 101002955 -CONQUER).

Appendix A: Dyson equation for RWA dynamics

In the same spirit as in the main text, we can write the dyson equation and more particularly the self energy for the Hamiltonian H RW A :

where the operator X is X = σ + = σ x + iσ y and X γ are superoperators labeled by indices γ ∈ {+, -} such that X + = X• and X -= •X. The bath correlation functions are given by

Appendix B: Validity of the rotating wave approximation

In this section, we will determine the validity criterion for the RWA described in the main text. For this, we will consider the evolution of the density matrix described by the following Hamiltonian:

we recall that the evolution of the density matrix is given by ρ tot (t) = e -iHt ρ tot (0)e iHt , where e iHt is the time evolution operator for the total Hamiltonian. By moving to the interaction picture according to the unperturbed Hamiltonian H RW A , the following identites can be found

here H CR (τ ) = e iH RW A τ H CR e -iH RW A τ and T t is the real-time ordering operator, which for any operator product changes the order such that each operator has only later operators to the left and earlier operators to the right.

From this expression, we see that the rotating wave approximation is valid when the initial density matrix remains unchanged by the evolution according to the Hamiltonian H CR (t), i.e when:

Since this condition is in general not satisfied, we perform an expansion according to the orders in g, and we seek the condition such that the lowest orders are negligible. In order to manage this expansion, it is convenient to rewrite the evolution operator

where we used the interaction picture on the real-time evolution operator of the Hamiltonian H RW A = H S +H B +H SB . By using the definitions of H B and H S from the main text, the following identities can be found for the part without coupling with the bath

this expression simply corresponds to the rotating frame of the counter rotating terms. Moreover, since the couplings are linear in g, we can expand the exponentials for each order. For the first order we get ρ CR,( 1) tot

Appendix C: Comparison with Born Master equation

In this appendix we show that many features of the emitter dynamics discussed in the main text are missed if one solves the system-bath dynamics within the Born Master Equation. As we mentioned in the main text, this corresponds to perform a non-self-consistent NCA calculation in which the propagator entering the self-energy is the bare emitter propagator rather than the fully dressed one. Therefore one expect that this approach is not able to capture the physics at strong system-bath couplings. To see this we consider in Fig. C the dynamics of the spontaneous emission as in the main text and compare the NCA and Born Markov dynamics. In the left panel we consider the dynamics within RWA and for the emitter frequency on resonance with the lower edge of the photonic band. We show that the emergence of plateaux in the spontaneous emission, corresponding to freezing, are absent within Born-Markov. In fact this only features an exponential decay in time with a rate that grows with system-bath coupling. In the right panel we plot the dynamics at fixed value of the coupling g and for different detuning ∆. We note a large deviation for ∆ = 0, an emitter frequency at the center of the band in correspondence of the van Hove singularity.

Here the relxation dynamics is non-exponential and clearly cannot be captured by Born-Markov. Similar devations appear also for ∆ = -W , i.e. in correspondence of the lower edge of the photonic band. These result highlight the importance of the self-consistent dressing in our dynamical map, particularly at strong coupling. We conclude by noting that a substantial difference between NCA and Born-Markov appears also in the spectral properties of the quantum emitter. In particular the large broadening of the resonance observed in Sec. IV cannot be reproduced with a Born-Markov calculation, whose spectral function always features a narrow in gap peak which does not broaden and only slightly shifts in frequency.


A R Khojasteh

D Heitz

Generic energy function to predict particle positions in 4D-PTV

Recent developments in time-resolved Particle Tracking Velocimetry (4D-PTV) consistently improved tracking accuracy and robustness by introducing predictive algorithms such as Shake-The-Box (STB). We propose a generic nondimensional energy function to predict particle positions by leveraging additional information from coherent neighbour motions. The proposed energy function is a combination of position history, neighbour velocity and acceleration terms. Coherent neighbours are quantified by locally utilising Finite Time Lyapunov Exponent (FTLE) as an objective Lagrangian Coherent Structure (LCS) diagnostic method. Synthetic analysis of the wake behind a smooth cylinder at Reynolds number equal to 3900 showed that the optimal solution of the minimised energy function could be modelled as a function of the measurement uncertainties. The model was assessed with the 2D homogeneous isotropic turbulent flow (HIT) as well and was found to be case and Reynolds independent. Results of the 4D-PTV experimental study of the same wake flow configurations are re ported. We compared predicted positions with the optimised final positions of STB. It was found that the generic energy function succeeded in estimating particle positions with minimum deviation to the optimised positions in comparison with other techniques.

Generic energy function minimisation

This paper discusses a parametric approach in time-resolved Particle Tracking Velocimetry (4D-PTV) to obtain a generic particle position predictor function. In the classic 4D-PTV algorithms such as Shake-The-Box (STB) [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF], predicted particle positions are given to the optimisation process for further corrections. The optimisation can deal with slight deviations between the predicted position and the true position. However, the optimisation fails to find the true position if the deviation is large enough to have multi-candidates for a single particle at t n+1 . This shows the importance of having an appropriate prediction applicable in different flow complexities. The proposed idea is initiated by arguing that predictions in PTV techniques such as polynomial and Wiener filter [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF] focus on a single particle individually, while this single particle is not acting alone (Khojasteh, Yang, et al., 2021). The present study was designed as a complementary function for 4D-PTV algorithms such as STB [START_REF] Schanz | Shake-The-Box: Lagrangian particle tracking at high particle image densities[END_REF] and KLPT (Yang & Heitz, 2021). Khojasteh, Heitz, et al. (2021) showed that an energy function obtained from three terms of position, velocity, and acceleration could significantly reduce the bias error caused by the flow complexities compared with polynomial and Wiener filter predictors. However, an estimated position in a real experiment needs to be optimised for both bias error and measurement uncertainties. The derived energy function suggested by Khojasteh, Heitz, et al. (2021) has dimensional terms (position, velocity, and acceleration) with equal weights to predict the particle positions. In the present study, we defined weighting parameters to the energy function to reach the minimum estimation bias error. The weighted energy function can be written as,

J = History ( 1 n n i=1 (X i -y i ) 2 ) + α 1 ( Ẋn -ẏn ) 2 Coherent velocity + α 2 ( Ẍn -ÿn ) 2 Coherent acceleration , ( 1
)
where X is the estimated positions based on known y observations in time. α 1 and α 2 are dimensional velocity and acceleration weights, respectively. The energy function contains three terms.

The first is the least mean square minimisation problem of the polynomial predictor based on the history of particles. The second and the third terms are neighbouring coherent velocity and accelerations. Khojasteh, Heitz, et al. (2021) showed that coherent velocity and acceleration terms could be computed by using Lagrangian Coherent Structures (LCS) to determine. Equal weights of α 12 mean that the impact of the neighbouring coherent velocity is at the same level as the history of the target particle. However, the particle history should be the most significant signal in the estimation process. These terms can be non-dimensionalised based on turbulent integral scales as,

x ′ = x D , y ′ = y D (2) ẋ′ = ẋ U , ẏ′ = ẏ U (3) ẍ′ = ẍ D U 2 , ÿ′ = ÿ D U 2 , ( 4
)
where D is an integral length scale, and U a velocity reference. Therefore, the energy function will become

J = D 2 ( 1 n n i=1 (X ′ i -y ′ i ) 2 ) + α 1 U 2 ( Ẋ′ n -ẏ′ n ) 2 + α 2 U 4 D 2 ( Ẍ′ n -ÿ′ n ) 2 , (5)
which can be simplified as,

J ′ = 1 n n i=1 (X ′ i -y ′ i ) 2 + α ′ 1 ( Ẋ′ n -ẏ′ n ) 2 + α ′ 2 ( Ẍ′ n -ÿ′ n ) 2 . (6) α ′ 1 = α 1 ( U D ) 2 , α ′ 2 = α 2 ( U D ) 4 (7)
Two non-dimensional α ′ 1 and α ′ 2 weights determine how much velocity and acceleration signals can constrain the overall energy function minimisation process. In equation 6, weight of the target position history is directly linked with two α ′ 1 and α ′ 2 values. If both terms are set for any range below one, the history is the most significant signal. We need to minimise the energy function for finding the optimal solution for the predictor function. The minimisation process of J ′ starts by solving ∂J ′ ∂a = 0. The optimal solution of the minimised energy function leads to a minimum estimation error of the predictor function. We plot the RMS error in a range of α ′ weights to find the optimal solution.

Synthetic evaluation

Since a real experiment is associated with uncertainties and inaccuracies, it is crucial to have an appropriate uncertainty estimation of each term in the energy function (see equation ( 6)). To this end, we employed data reported by the 1st LPT challenge [START_REF] Sciacchitano | Main results of the first Lagrangian Particle Tracking Challenge[END_REF] as a reference starting point to estimate inaccuracies that might be introduced into a predictor function. The 1st LPT challenge assessed position estimation accuracy of six time-resolved tracking algorithms, including the coherency based tracking (Khojasteh, Yang, et al., 2021) for particle densities from 0.05 to 0.2ppp. The non-dimensionalized position ϵ X /D, velocity ϵ Ẋ /U ∞ , and acceleration ϵ Ẍ (U ∞ 2 /D) errors were employed to create the synthetic data at each particle density. As a result of the LPT challenge, the averaged RMS position error was 0.005mm where the integral scale was D = 10 mm at the density of 0.12ppp. The averaged velocity and acceleration error at the same density were found to be 0.01 U ∞ and 0.3U ∞ 2 /D, respectively. This shows that the acceleration estimation has at least an order of magnitude higher estimation error than other terms. We introduced these errors into the Lagrangian trajectories of 2D homogeneous isotropic turbulent and 3D cylinder wake flow and then computed the correlation of the final estimation error and α ′ weights. For such a scenario, we can expect the optimal solution should have less α ′ 2 weight than α ′ 1 . Figure 1 shows the RMS prediction error based on two weighing parameters. Since the acceleration error is higher than other terms in both synthetic data, the minimum error happens with relatively lower α ′ 2 values (i.e., α ′ 1 /α ′ 2 >> 1). Both weighting parameters are found to be significantly smaller than << 1. This means that the history of the target particle with unit weight is the most valuable signal in such a scenario. Comparing the error behaviour in 2D HIT and 3D cylinder wake flow cases shows that both cases need similar weighting parameters and require α ′ 1 /α ′ 2 >> 1 with small weighting magnitudes α ′ << 1.

The result of the minimum energy function from two synthetic test cases suggests that the optimal solution can be directly linked with uncertainty levels. Therefore, we designed three further parametric test cases to determine the mentioned statement. As shown in Figure 2, we increased the uncertainty level of each parameter solely while other terms are fixed to assess how the optimal solution behaves with respect to three position, velocity, and acceleration uncertainties. The optimal solution tends to linearly move away toward higher weighting magnitude with increased position uncertainty (see Figure 2 a). If we fix the position uncertainty, the optimal solution rotates with the same distance around the coordinate centre, depending on which parameter is increased (see Figure 2 bc). As a result of this assessment, an appropriate estimation of three uncertainty levels would provide enough information to set weighting parameters. On the other hand, we can model the optimal solution by roughly estimating position, velocity, and acceleration uncertainties by fitting a surface over the minimum α ′ 12 values. Therefore, we can set the weighting system depending on the uncertainty level of each term.

20th LISBON Laser Symposium 2022

Experimental evaluation

An experimental study of the cylinder wake flow at Reynolds number equal to 3900 (same value as the synthetic data) was performed. We designed an experimental setup with four cameras. Four CMOS SpeedSense DANTEC cameras with a resolution of 1280 × 800 pixels and the maximum frequency of 3 kHz are empowered. Cameras are equipped with Nikon 105 mm lenses. The first two cameras are positioned in backward light scattering, while the second two cameras receive maximum intensity signal in forward scattering. The calibration error was lower than 0.06 pixel and reduced to 0.04 after the volume self calibration. The volume of interest was 280 mm × 160 mm × 46 mm starting from roughly 4D downstream of the cylinder, knowing that the vortex formation zone ends at 4D. The aperture was set at 11 to achieve 46 mm depth of focus. We used an LED system to illuminate this large volume. The seeding particles were Helium Filled Soap Bubbles (HFSB) [START_REF] Scarano | On the use of helium-filled soap bubbles for large-scale tomographic PIV in wind tunnel experiments[END_REF] resulting desired intensity signal with appropriate particle size. However, bubbles are limited by three main factors in the wind tunnel experiments, including generation rate, lifetime, and image glare points. We placed 50 bubble generator nozzles with airfoil-shaped structures inside the wind tunnel chamber. The nozzles were far upstream of the measurement section to ensure a sufficient number of bubbles were created, and the main flow field was not disturbed by the existence of nozzles. The bubble lifetime is very short (less than 2 -3 minutes) inside the wind tunnel, mainly because they explode by passing through honeycomb layers. To overcome this issue, we injected bubbles inside the chamber for up to 5 minutes when the wind tunnel is off before starting the acquisition. We found that particles larger than three pixels create two glare points on two sides of the bubble. This requires more image treatments before running the 4D-PTV algorithm to avoid false particle reconstruction. However, the intensity of two glare points can diffuse and merge if the particle size is around two pixels. Therefore, we adjusted the camera magnification to reach two particle pixel sizes on average to surpass the glare point issue. One snapshot of the experiment is shown in Figure 3 To quantify the results of different schemes, we compared predictions with optimised positions obtained from STB Davis. As a result of the experiment, STB managed to successfully build nearly 12000 particles as shown in Figure 3.b. Noisy particle reconstruction of four time steps was used as an input of the prediction functions. We compared three techniques, Polynomial, Wiener filter, and coherent predictors, with final optimised positions. The deviation of position estimated of each technique is shown in Figure 3. The distribution shows that the coherent predictor has more accurate estimations within 1 pixel deviation from the optimised positions. Position estimations of Wiener filter and coherent predictors stay below 2.5 pixels deviation for nearly all particles. On the contrary, the Polynomial predictor has maximum deviation with STB Davis.

Conclusion

In the present study, we obtained a generic position predictor as a function of the measurement uncertainties. We performed DNS synthetic analyses of the 3D wake behind a cylinder at Reynolds number of 3900, and 2D homogeneous isotropic turbulent flow at a Reynolds number 3000. It was found that the optimal solution of the predictor function is a case and Reynolds independent model. We also used 4D-PTV experiment data to demonstrate that the generic energy function has better estimation errors than recently available techniques.

Figure 1 .

 1 Figure 1. Contour plots of the estimation error with respect to α ′ 1 and α ′ 2 weighting terms. The synthetic trajectories were generated based on the 1st LPT challenge estimation errors. a) the 3D wake behind a cylinder at a Reynolds number equal to 3900, b) 2D homogeneous isotropic turbulent flow at a Reynolds number equal to 3000

Figure 2 .

 2 Figure 2. Contour plots of the estimation error as a function of uncertainty levels. White arrows show the changing direction of the optimal solution (i.e. minimum estimation error) by increasing each uncertainty term.

  .b. Trajectory results of the current experiment with superimposed vorticity iso-surfaces are shown in Figure 3.b.

Figure 3 .

 3 Figure 3. Experiment study of the cylinder wake flow at Reynolds 3900 with four cameras: a) Experiment normal pdf results of particle position error in x direction of three predictors. Each predictor is compared with final optimised positions of STB Davis; b) Side view of particle trajectories superimposed by vorticity iso-surfaces.


Ali Rahimi Khojasteh
email: ali.rahimi-khojasteh@inrae.fr

Dominique Heitz

A R Khojasteh

Diffusion properties of time-resolved trajectories in the wake behind a cylinder

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Recent advances in time-resolved three-dimensional particle tracking velocimetry (4D-PTV) extend the possibility of exploring the Lagrangian viewpoint of turbulent properties along particle trajectories. The recent work of [START_REF] Viggiano | Lagrangian diffusion properties of a free shear turbulent jet[END_REF] showed how we could obtain fundamental Lagrangian turbulence quantities, particularly the Lagrangian structure function scaling constant C 0 , from anisotropic and inhomogeneous dynamics of a jet flow by normalising the trajectories based on local Eulerian scales. The C 0 constant in the Lagrangian framework is in a similar role as the Kolmogorov constant in the Eulerian framework [START_REF] Viggiano | Lagrangian diffusion properties of a free shear turbulent jet[END_REF]. In the present study, we are interested in examining the stationarisation process proposed by [START_REF] Viggiano | Lagrangian diffusion properties of a free shear turbulent jet[END_REF] in the wake behind a smooth cylinder and studying Lagrangian statistics from Spatio-temporal trajectories.

Lagrangian diffusion properties

Taylor's turbulent diffusion theory [START_REF] Taylor | Diffusion by continuous movements[END_REF] has been used widely to study homogeneous isotropic turbulent (HIT) flows. In a given time τ , Taylor's theory computes the Lagrangian twopoint correlation function R L uu (τ ) for an ensemble of particle trajectories based on the mean square displacements of particles σ 2 (τ ) that can be written as,