Dynamical dimerization phase in Jaynes–Cummings lattices

We report on an emergent dynamical phase of a strongly-correlated light–matter system, which is governed by dimerization processes due to short-range and long-range two-body interactions. The dynamical phase is characterized by the spontaneous symmetry breaking of the translational invariance and appears in an intermediate regime of light–matter interaction between the resonant and dispersive cases. We describe the quench dynamics from an initial state with integer filling factor of a finite-sized array of coupled resonators, each doped with a two-level system, in a closed and open scenario. The closed system dynamics has an effective Hilbert space description that allows us to demonstrate and characterize the emergent dynamical phase via time-averaged quantities, such as fluctuations in the number of polaritons per site and linear entropy. We prove that the dynamical phase is governed by intrinsic two-body interactions and the lattice topological structure. In the open system dynamics, we show evidence about the robustness of dynamical dimerization processes under loss mechanisms. Our findings can be used to determine the light–matter detuning range, where the dimerized phase emerges.


Introduction
The development of technology encompasses a broad range of opportunities to harness quantum phenomena. For instance, it is possible to manipulate light-matter quasiparticles or polaritons which have new properties such as stimulated scattering [1,2], lasing [3][4][5], parametric amplification [6][7][8], and superfluidity [9,10]. These characteristics can be used to enhance the experimental realization of polaritonic devices, such as semiconductor microcavities, where the coupling between quantum-well excitons and cavity photons gives rise to hybrid light-matter quasiparticles [11]. In the microwave regime, superconducting circuits based on Josephson junctions also allow us to harness light-matter interaction for simulating strongly correlated phenomena with light [12][13][14][15][16][17][18][19][20][21][22]. The underlying physics of light-matter based quantum simulators is governed by the Jaynes-Cummings-Hubbard (JCH) model [23][24][25] which describes the dynamics of coupled-resonator arrays (CRAs), each doped with a two-level system (TLS). In this case, the manipulation of polaritonic excitations locally depends on the detuning between the light and matter frequencies but also is largely influenced by the lattice structure. As the detuning increases from the resonant to the dispersive regime, the system transits from the Mott-insulating state characterized by the hybridization of light and matter states to a superfluid phase of photons [23][24][25][26][27][28][29][30].
In this work, we demonstrate by using numerical calculation and an analytic model that during the phase transition from the Mott-insulating to superfluid state, a dynamical dimerization phase (DDP) emerges, which is characterized by the spontaneous symmetry breaking of the translational invariance. As dynamical dimerization, we refer to the dynamics of a finite-sized Jaynes-Cummings (JC) lattice that exhibits resonances related to the two-sites JC lattice. In order to identify the new regime of DDP, we analyze the purely unitary quench dynamics of few-body JC lattices. In particular, we consider a quantum quench from an initial state with integer filling factor in a JC dimer, which has been proven useful to simulating second-order like phase

The model
The JCH model [23][24][25] describes a lattice of L interacting coupled QED resonators, each supporting a single mode of the electromagnetic field which interacts with a TLS. This situation is schematically shown in figure 1 In the resonant regime, Δ=ω 0 −ω=0, the hybridization of the light-matter yields to localized polariton excitations (Mott-insulator phase), while in the dispersive regime, Δ?g, the system is dominated by photonic excitation behavior (superfluid phase) [23,24]. This phase transition can also be described as a transition driven by the photon blockade effect from the Mott-insulator phase, where the intersite polariton exchange is forbidden, so effectively J ij =0, to a superfluid phase dominated by a uniform photon hopping J ij =J, in both cases, there is no cavity-embedded effect involved. As we demonstrate in section 4, the intermediate regime of light-matter interaction, which we define in the range 1<Δ/g<10, can be identified by the parameter with L the number of nonlinear coupled resonators, J the hopping parameter, and ν j the connectivity of node j. Notice that k i =0 for the resonant case and k i =J for the dispersive case. This way, the origin of translational symmetry breaking can be explained by introducing a local order parameter of μth phase Y m i , with μ=(I, II, III) represents the resonant, intermediate, and dispersive case, respectively. Indeed, , which means that in the resonant and dispersive regimes there is translational invariance. Since the order parameter shows a spatial dependence in an intermediate regime , then translational symmetry is broken. As a consequence, a DDP will happen due to intrinsic two-body interactions and the connectivity of each lattice site.
Using the above defined polaritonic basis, the Hamiltonian (1) can be rewritten as [23,27] where the matrix elements . The first term in equation (2) stands for the local polaritonic energy with an anharmonic spectrum and gives rise to an effective on-site polaritonic repulsion. This is analog to the on-site photon repulsion in the Bose-Hubbard model [36]. The last term in equation (2) represents the polariton hopping between resonators. This interaction also allows the interchange of polaritonic species of one or both sites involved [23,27,32], leading to a quite involve quantum dynamics.
A detailed understanding of the equilibrium properties of the JCH model (2) resorts on approximated analytical solutions [37] or numerical approaches such as density matrix renormalization group [38][39][40][41]. In nonequilibrium situations, one can understand the underlying physics using the time-evolving block decimation algorithm [42][43][44], or simplifying the description using effective Hilbert spaces [32,[45][46][47][48][49]. The latter is particularly appropriate for studying the quench protocol presented in this article, as we consider the closed system scenario.

Quench dynamics in a JCH dimer
In this section, we introduce a sudden quench protocol [31] and its effects on the quantum dynamics of the JCH dimer. Also, we provide quantitative explanations of time-averaged quantities such as fluctuations in the number of polaritons per site and linear entropy using an effective Hilbert space. In section 4, we will show that the underlying physics of the quench dynamics allows us to understand the emergent DDP. We stress that DDP occurring in the JCH lattice happens in the frequency regime  w Jn g n n [23], where the rotating wave approximation holds.
For each detuning Δ, we set the initial condition as the lowest energy state with integer filling factor of one excitation per site, that is, (L is the number of lattice sites) which corresponds to a Mottinsulating state at hopping rate J=0. Then, at time t=0, the parameter J is suddenly quench to a new value ¹ J 0 f such that the Hamiltonian has changed to H=H JCH (J f ). Hence, the JCH lattice dynamics is described by , and ρ i the reduced density matrix of the leftmost or rightmost site of the JC dimer.
The quench protocol described above allows us to simulating second-order like phase transition captured via the Var(n i ), see [31], which is analog to the adiabatic dynamics studied in [23]. Also, the simulated phase transitions in a JCH lattice can be characterized via linear entropy, as shown in this article, which implies that the observation of the emergent DDP is independent of the choice of the order parameter.
Effective Hilbert space for a closed system.-In order to introduce an effective description of the system dynamics, it is useful to consider the JCH Hamiltonian written in the polaritonic basis, see equation (2).
, the JCH Hamiltonian (2) may lead to processes such as the exchange of polaritonic species, or the interchange of polaritonic species of one or both sites involved (i, j). In this case, the full quantum dynamics should involve all states within the two-excitations subspace, which we define as , which results in fast oscillating contributions, and we can apply the rotating-wave approximation [23,27,32]. Figure 2 shows that the interchange of polaritonic species is suppressed over the evolution. Here we plot the populations of above defined states as a function of time. We identify populations as Due to the symmetry of the JC dimer, it is cleat that = figure 2). In this work, we carry out numerical calculations with the quantum toolbox in Python QuTiP [50,51].
Since the interchange of polaritonic species is suppressed over the evolution, we can introduce an effective Hilbert space involving states of the lower polaritonic branch for describing the quench dynamics. In this case, the effective Hamiltonian reads . Thus, the full dynamics can be solved analytically by diagonalizing the above 3×3 matrix. Starting from the initial condition | where the probability amplitudes are

Time
, and ρ i the reduced density matrix of the leftmost or rightmost site of the JC dimer. Thus, the time-averaged variance reads ⎡ Figure 3(a) shows the behavior of Var(n i ) as a function of log 10 (Δ/g) calculated from the full numerics (red diamonds) and the analytical prediction (continuous blue line) in equation (6). We see a good agreement between both predictions as the relative error shows in figure 3(b). Also, equation (6)  (a=c) (see figure 2(b)). Also, | | = b J 2 and Ω 0 =4J, so the asymptotic value of Var(n i ) reads It is worth mentioning that the analytical result (6) represents the hallmark for the dimer dynamics. In section 4, we will prove that as one increases the number of lattice sites, the time-averaged variance (6) allows us to identify resonances due to intrinsic short-and long-range two-body interactions, which govern the DDP.
On the other hand, one can also characterize the dimer dynamics via the linear entropy of the reduced density matrix as mixedness measure [52]. First, notice that the quantum state (4) is already in its Schmidt decomposition, which leads to a diagonal reduced density matrix   Figure 4(a) shows the behavior of E as a function oflog 10 (Δ/g) calculated from the full numerics (red diamonds) and the analytical prediction (continuous blue line) in equation (9). We see a good agreement between analytical and numerical predictions as shown in figure 4(b). We can also estimate the asymptotic value of the time-averaged entropy as one increases the ratio Δ/g, it reads In what follows, we will use the time-averaged variance (6) for demonstrating the existence of DDP in the intermediate regime of light-matter interaction. This choice establishes the physical framework for the subsequent discussion, but a similar analysis with the linear entropy leads to the same conclusion about DDP.

Dynamical dimerization phase
4.1. Closed system Emergent dynamical critical phenomena following a quantum quench have experienced much interest in recent years . Here, we demonstrate an emergent DDP as we extend the JC lattice to three (trimer) and four (tetramer) sites (see figure 1), and in an intermediate regime of light-matter coupling strength, that is, 1< Δ/g<10. We define the concept of dimerization as the dynamical process where short-range and long-range  (6)). We use parameters g=10 −2 ω, J=10 −4 ω, where ω is the resonator frequency.
two-body interactions govern the quantum dynamics. Here, short-range (long-range) two-body interaction is due to the direct (mediated) exchange of polaritons. We will prove that the connectivity associated with each node of the lattice plays a crucial role to define DDP.
Let us start the discussion with a three-sites JC lattice initialized in the state where the subindexes i, j, and k refer to the leftmost, center, and rightmost lattice site, respectively. If we let the system evolves according to the quantum quench described is section 3, one should impose the conditions for neglecting interchange of polaritonic species between nearest-neighbor sites, that is, {| | - . As for the dimer, only the lower polaritonic branch will be activated and the dimension of the effective Hilbert space () is given by (N + d−1)!/N!(d−1)!, where N is the number of excitations that should be distributed into d lattice sites. In our case, N=3 and d=3 results in ( ) =  dim 10. At time t, the wave function may be written as a linear combination of states belonging to the three-excitations subspace, that is where we define states |y ñ Notice that some probability amplitudes are equal due to symmetry of the trimer with respect to the lattice center ( j). Here, the state (11) will be computed via full numerics. Figure 5 shows the ratio between the absolute value of the time-averaged nearest-neighbor correlation function C ij (next nearest-neighbor correlation function C ik ) of the trimer, and the dimer variance (6). Twopoint correlation functions are ( ) , where τ=J −1 . Here, we identify two critical values of detuning, vertical dashed lines, Δ/g=(2.43, 2.73) within the intermediate regime of lightmatter interaction, 1<Δ/g<10. At these critical points the trimer experiences dynamical dimerization processes, where short-range correlations rule the dynamics at Δ/g=2.43, while at Δ/g=2.73 a combination of both short-and long-range correlations govern the dynamics. The resonances shown in figure 5 demonstrate that the intrinsic dimer dynamics, characterized by the time-averaged variance (6), governs the quench dynamics of the trimer. Furthermore, the dimer variance allows us to identify short-range and long-range twobody interactions. The former is a consequence of direct cavity-cavity coupling of sites (i, k), while the latter results from an indirect interaction between sites (i, k) mediated by the center lattice site j. Notice that figure 5 also exhibits an anti-resonance at Δ=1.82g (continuous vertical line). At this point, no dynamical dimerization happens, and the JC lattice remains approximately in the Mott insulating state.  (9)). We use parameters g=10 −2 ω, J=10 −4 ω, where ω is the resonator frequency.
Let us discuss the results for the tetramer. Figure 6 shows the ratio between the absolute value of timeaveraged two-point correlation functions C ij , C ik , C il of the tetramer, and the dimer variance in equation (6). We order the lattice sites from left to the right according to indexes (i, j, k, l). It is noticeable that a larger JC lattice also exhibits resonances at critical values of detuning Δ/g=(2.43, 2.73) (vertical dashed lines), and the antiresonance at Δ=1.82g (continuous vertical line). Moreover, the two-point correlation function C il , associated with edges of the lattice, has the same resonance (Δ=2.43g) as compared with nearest-neighbor correlation function C ij . The latter suggests that for a finite one-dimensional JC lattice of L sites, the number of resonances associated with dimerization processes corresponds to a universal number of different connectivities of the lattice, that is, connectivity ν=1 for lattice edges, and connectivity ν=2 for bulk lattice sites. These results are a consequence of the broken translational symmetry.

Open system
A realistic implementation of a strongly-correlated light-matter system should consider the system-bath interaction which leads to loss mechanisms in the initial state preparation and along the dynamics, e.g. if we Figure 5. Ratio between the absolute value of time-averaged nearest-neighbor correlation function C ij (next nearest-neighbor correlation function C ik ) of the trimer, and the dimer variance in equation (6). Vertical dashed lines, from left to right, indicate critical values of detuning Δ/g=(2.43, 2.73) where dynamical dimerization processes happen. Continuous vertical line stands for the antiresonance. We use parameters g=10 −2 ω, J=10 −4 ω, where ω is the resonator frequency, and we consider up to 5 Fock states per resonator. Figure 6. Ratio between the absolute value of time-averaged two-point correlation functions C ij , C ik , C il of the tetramer (four-sites JC lattice), and the dimer variance in equation (6). Vertical dashed lines, from left to right, indicate critical values of detuning Δ/g = (2.43, 2.73) where dynamical dimerization processes happen. Continuous vertical line stands for the anti-resonance. We use parameters g=10 −2 ω, J=10 −4 ω, where ω is the resonator frequency, and we consider up to 5 Fock states per resonator. consider a experimental realization based on superconducting circuits [15,17]. In these experiments, the dissipative dynamics is described by a Markovian Lindblad master equation ( . We consider the same loss mechanisms for each lattice site including energy relaxation, dephasing; and photon losses at rates γ, γ f , and κ, respectively.
In order to prepare the initial state, we propose to include an ancillary TLS on each lattice site, which interacts with the cavity mode. In this case, the Hamiltonian describing a single lattice site reads where ω A is the ancilla frequency, g A the ancilla-cavity coupling strength, and H i JC the JC Hamiltonian of site i. The initialization protocol makes use of Gaussian and Stark pulses as described in [76]. First, we let the system to cold down to its ground state | ⨂ | | y ñ = -ñ ñ = 0, i . Second, we apply individual Gaussian π pulses acting upon each ancilla TLS in order to prepare the state ⨂ | | -ñ ñ = 0, i . Third, a Stark pulse is applied to each ancilla TLS bringing it into resonance with its respective lattice site, i.e. w = -E A 1 , during a time interval Δτ=π/(2g A t −− 1 ). In this way, the strong lattice site-ancilla interaction governed by equation (13) yields the desired initial state | ⨂ | | y ñ = -ñ ñ = 1, i . Notice that one should satisfy the condition | |  - . The latter avoids unwanted population of the state ⨂ | | +ñ ñ = 1, i . Then, the ancilla-site interaction is suppress by applying a Stark pulse tuning ω A below the frequency -E 1 . We stress that Stark pulses can be implemented by means of external magnetic fluxes applied upon superconducting quantum interference devices that form a transmon qubit [77,78].
We have numerically calculated the initial state preparation and the sudden quench using the equation (12), using physical parameters taken from state-of-the-art circuit QED setups [15,17,78,79]. Figures 7(a), (b) shows the ratio ( ) ( ) C n Var ij j Trimer Dimer as a function of Δ/g, where (C ij ) Trimer stands for the absolute value of the timeaveraged correlation function, and Var(n j ) Dimer corresponds to the variance of the dimer. As seen in figure 7(a), we identify resonances corresponding to the critical values Δ/g=(2.57, 3.08), vertical dashed lines, within the intermediate regime of light-matter interaction, 1<Δ/g<10. In the same way, as in the closed system dynamics 4.1, the trimer experiences dynamical dimerization processes, where short-and long-range correlations dominate over dissipation. As we increase the energy relaxation and dephasing rates of TLSs, resonance peaks are spreading throughout to a wider detuning range, see figure 7(b). These results show evidence about the stability of dynamical dimerization processes that happen in a finite-sized JC lattice, and allow us to establish a parameter threshold for the appearance of DDP in the dissipative case.

Conclusions
In summary, we have reported on the emergence of a DDP in a finite-sized JC lattice as a result of a quantum quench from an initial state with integer filling factor. We have thoroughly analyzed the quench dynamics in a close two-sites JC lattice, which allows us to obtain analytical results for time-averaged order parameters such as the local variance of the number of polaritons, and the linear entropy. Further, these order parameters can be used to analyze and predict the resulting quenched dynamics for more complex architectures. When comparing the dimer variance with two-point correlation functions of the trimer and tetramer, it allows us to determine critical values for the detuning where dynamical dimerization processes happen. Recognizing resonances and anti-resonance for detuning values, in turn, allow controlling what kind of correlation dominates over the dynamics be short-range or a combination of both short-and long-range, and may also allow controlling polariton propagation along the lattice. We stress that the intrinsic dimer dynamics, characterized by the timeaveraged variance (6), governs the quench dynamics of closed finite-sized JC lattices, and we expect similar results as one increases the number of lattice sites.
In a realistic situation, it is necessary to include dissipative mechanisms in the state preparation and over the quench dynamics. Considering parameters of state-of-the-art circuit QED technology, permit numerical results to show that as one increases the coherence times of TLSs and cavities, two sharp resonance peaks become more evident. These results demonstrate that DDP remains in the dissipative case. We conjecture that for a finite onedimensional JC lattice of L sites, the number of resonances associated with dimerization processes corresponds to the number of different connectivities of the lattice, that is, connectivity ν=1 for lattice edges, and connectivity ν=2 for bulk lattice sites. Our findings could be tested with state-of-the-art quantum technologies. For instance, in trapped ions technology, the JCH model has been theoretically proposed in [33] and physically implemented in [34]. In superconducting circuits, the JC dimer has been implemented in [15]. In this case, homodyne signal detection may allow measuring the local variance of photon number, which can also be used as an order parameter.