Robust Preparation of Many-body Ground States in Jaynes-Cummings Lattices

Strongly-correlated polaritons in Jaynes-Cummings (JC) lattices can exhibit quantum phase transitions between the Mott-insulating and superfluid phases at integer fillings. The prerequisite to observe such phase transitions is to pump polariton excitations into a JC lattice and prepare them into appropriate ground states. Despite previous efforts, it is still challenging to generate many-body states with high accuracy. Here we present an approach for the robust preparation of many-body ground states of polaritons in finite-sized JC lattices by optimized nonlinear ramping. We apply a Landau-Zener type of estimation to this finite-sized system and derive the optimal ramping index for selected ramping trajectories, which can greatly improve the fidelity of the prepared states. With numerical simulation, we show that by choosing an appropriate ramping trajectory, the fidelity in this approach can remain close to unity in almost the entire parameter space. This approach can shed light on high-fidelity state preparation in quantum simulators and advance the implementation of quantum simulation with practical devices.


INTRODUCTION
The Jaynes-Cummings (JC) model is a prototype for studying light-matter interaction, where a quantum two-level system is coupled to a cavity mode [1]. This model has been utilized to study cavity or circuit quantum electrodynamics (QED) in a wide range of systems, from individual particles in the atomic scale to collective modes in mesoscopic devices [2][3][4]. More recently, advances in device fabrication and quantum technology enabled the exploration of many-body physics in arrays of JC models, i.e., JC lattices, which can be realized with optical cavities coupled to defects in semiconductors [5][6][7][8][9] and superconducting circuit QED systems [10][11][12][13][14][15]. The light-matter coupling in a JC model induces intrinsic nonlinearity in the energy spectrum, which can be mapped to an onsite repulsive interaction between polariton excitations. The competition between this onsite interaction and polariton hopping between neighboring sites gives rise to rich manybody physics for strongly-correlated polaritons in JC lattices, such as quantum or dissipative phase transitions and photon blockade effects [16][17][18]. Moreover, when the counterrotating terms in the qubit-cavity interaction cannot be neglected, the system becomes a quantum Rabi lattice, where distinctively different many-body phase transitions have been studied [19,20]. One effect of particular interest is the quantum phase transitions between the Mott-insulating (MI) and superfluid (SF) phases for polaritons in JC lattices at integer fillings, featured by the occurrence of off-diagonal longrange order in the correlation functions. It was shown that such phase transitions can be observed in coupled cavity arrays [5][6][7] and multi-connected JC lattices [13][14][15].
The prerequisite to observe the MI-SF phase transitions * email:ltian@ucmerced.edu is to pump polariton excitations into a JC lattice and prepare them into appropriate ground states. However, preparing many-body ground states is a challenging task in engineered systems such as quantum simulators [21][22][23] and adiabatic quantum computers [24,25]. A number of approaches have been studied to tackle this problem, including adiabatic quantum evolution [26][27][28], quantum shortcut method by applying counter-diabatic interactions [29,30], quantum phase estimation via quantum Fourier transformation [31,32], variational quantum eigensolver [33,34], full quantum eigensolver [35], and engineered dissipative environment for the preparation and stabilization of entangled states [36][37][38][39]. Despite these efforts, it is still hard to generate desired manybody states with high fidelity in the noisy intermediate-scale quantum (NISQ) era [40], in particular, for systems working with excitations such as the JC lattices. The barriers to generating desired many-body states efficiently and accurately include the lack of a priori knowledge of the energy spectrum, the difficulty in engineering complicated counter-diabatic interactions, the rapid decrease of the energy gap and quick increase of the dimension of the Hilbert space with the size of the quantum simulators, and the finite decoherence times in NISQ devices. Furthermore, many-body states in stronglycorrelated systems can be highly entangled, unknown, and hence, often impossible to be generated with pre-programed quantum logic gates.
Here we study the robust generation of many-body ground states in finite-sized JC lattices at unit filling using optimized nonlinear ramping. In previous works [41][42][43], it was shown that nonlinear ramping can reduce diabatic transitions to excited states or the production of domain walls when a manybody system in the thermodynamic limit evolves across a quantum critical point due to the scaling of the phase transition. We apply nonlinear ramping to a finite-sized system, where the energy gap between the ground and the excited states remains finite. By exploiting a Landau-Zener type of estimation [44,45] and the spectral feature along a selected ramping trajectory, we derive the optimal ramping index for the trajectory, which can significantly improve the fidelity of the prepared state. Our estimation agrees well with the result from our numerical simulation of the ramping process. Moreover, we show that by selecting an appropriate trajectory for a given set of target parameters in combination with the optimal ramping index, the fidelity can remain close to unity in almost the entire parameter space. The ramping trajectory can be adjusted by varying the initial parameters or the ratios between the ramping indices for different parameters. The initial states of this nonlinear ramping process can be prepared with high accuracy by applying engineered pulse sequences [46] when tuning the system parameters to either the deep MI regime with no hopping between adjacent unit cells or the deep SF regime with diminishing light-matter coupling.
JC lattices have been implemented with superconducting quantum devices in recent experiments [18,47,48], and the ramping process studied here is within reach of current technology [49][50][51]. Using practical parameters from the experiments, we show that high fidelity can be achieved for the prepared states on a time scale much shorter than the observed decoherence times of these devices. Meanwhile, the approach of optimized nonlinear ramping for finite-sized systems is general and can be applied to many other models. The study of state generation in finite-sized systems with this approach can provide insights to the problem of preparing complex quantum states in quantum computers. Our result can hence shed light on the high-fidelity preparation of manybody states in engineered quantum systems such as quantum simulators, and advance the implementation of quantum simulation with NISQ devices.

Model and quantum phase transition
Consider the JC lattice depicted in Fig. 1a. Here each unit cell contains a qubit coupled to a cavity mode with coupling strength g, and adjacent unit cells are connected via photon hopping with hopping rate J. The total Hamiltonian of this model is is the Hamiltonian of uncoupled JC models with ω c the frequency of cavity modes, ω z the level splitting of the qubits, a j (a † j ) the annihilation (creation) operator of the j-th cavity mode, and σ j± , σ jz the Pauli operators of the j-th qubit, and is the photon hopping between neighboring unit cells. with light-matter coupling g and hopping rate J. b Single-particle density matrix ρ 1 (1, 4) vs hopping rate J and detuning ∆ for a finitesized lattice at unit filling with N = L = 6. Here we let g be the energy unit with g ≡ 1.
In the limit of J = 0, the JC lattice is composed of isolated JC models. The ground state at unit filling, where the number of polaritons N is equal to the number of lattice sites L, is with one polariton excitation occupying the state |1, − per site, which is in the deep MI regime. States with more than one excitation at the same site are energetically unfa-vorable due to the effective onsite interaction. In the opposite limit of g = 0 at finite hopping rate J, the cavity modes are decoupled from the qubits. The hopping Hamiltonian (2), now the dominant term, can be transformed to the momentum space under the periodic boundary condition with H int = −2J ∑ k cos (k) a † k a k , where a k (a † k ) is the annihilation (creation) operator of a collective cavity mode at the quasimomentum k = πm/N with integer m ∈ [−(N − 1), N] and a k = ∑ j a j e ik· j / √ N. At ∆ < 0 with the cavity frequency below the qubit energy splitting, the ground state at unit filling is with all polaritons occupying the k = 0 mode, which is a nonlocal state in the deep SF regime.
With the mean-field approximation [5][6][7] and numerical methods [8,9,[13][14][15], it was shown that quantum phase transitions between the MI and SF phases due to the competition between the onsite interaction and the photon hopping can occur in the intermediate regimes of the parameter space in JC lattices. For a finite-sized lattice with N = L = 6, we numerically calculate the many-body ground states using the exact diagonalization method. In this finite-sized system, the energy separation between the ground and the excited states decreases as the parameters approach the intermediate regimes, but maintains a finite energy gap. The spatial correlation in the many-body ground state |G can be characterized by the normalized single-particle density matrix defined as [53,54] which reveals the off-diagonal long-range order of the state. The single-particle density matrix decreases algebraically with the spatial separation |i − j| in the SF phase and decreases exponentially in the MI phase [13][14][15]. For a finite |i − j|, ρ 1 (i, j) of the SF phase is much larger than that of the MI phase. In Fig. 1b, we plot our numerical result of ρ 1 (1, 4) as functions of the hopping rate J and the detuning ∆, with the coupling g as the energy unit (g ≡ 1). It can be seen that ρ 1 (1, 4) increases with J at arbitrary detuning. In the deep MI regime with J = 0, ρ 1 (1, 4) = 0 with the polaritons localized in the lattice. In the deep SF regime, ρ 1 (1, 4) can approach unity. This result clearly indicates the occurrence of the MI-SF phase transition in the thermodynamic limit.

State initialization
We present methods to pump N = L polaritons to the JC lattice in the limiting cases of J = 0 and g = 0, respectively, by applying engineered pulses. The polaritons are pumped into the many-body ground states at the corresponding parameters. These states will be used as the initial state of the nonlinear ramping approach.
In the deep MI limit of J = 0 and finite g, the ground state is given by (5) with each unit cell in the polariton state |1, − . Because the unit cells are decoupled, we can perform a Rabi rotation between the states |g 0 and |1, − on each The corresponding Rabi frequency can be derived as Ω d1 = |ε cos(θ /2)| following Eq. (4). The duration of the Rabi flip from the initial state |g 0 to the final state |1, − is τ d1 = π/2Ω d1 . To prevent the driving pulse from inducing unwanted transitions to higher states such as |1, + , it requires that |ε| ≪ g.
In the deep SF limit of g = 0 and finite J, the ground state is given by (6) with all polaritons occupying the collective (nonlocal) mode a k=0 . To generate this state, we introduce an auxiliary qubit with Pauli operators σ 0± , σ 0z . This qubit has the Hamiltonian which includes the qubit energy splitting ω 0 , a driving on the qubit with amplitude ε (t) and frequency ω L , and a tunable coupling between the qubit and the cavity modes with coupling strength g d (t). By choosing ω L , ω 0 to be both in resonance with the mode a k=0 , we have H d2, in the rotating frame. The first term of H d2,r generates a Rabi rotation on the auxiliary qubit, and the second term is the coupling between the auxiliary qubit and the mode a k=0 . Both terms can be turned on and off within nanoseconds, as has been demonstrated in recent experiments on superconducting transmon qubits [47,48]. As the qubits in the JC lattice are decoupled from the cavities, the state of mode a k=0 is only affected by its coupling to the auxiliary qubit. Let the initial state of the coupled system of mode a k=0 and the auxiliary qubit be |0, ↓ . To generate the state (6), we utilize the approach in [46] to design a pulse sequence, which can be implemented by switching on ε (t) and g d (t) alternately. The unitary operator for this pulse sequence is where the unitary operator C l (l ∈ [1, N]) incurs a Rabi flip on the auxiliary qubit by applying a driving pulse with amplitude ε for a duration of τ cl = π/2|ε|, and the unitary operator Q l enables the exchange of excitations between the auxiliary qubit and the mode a k=0 by turning on the coupling g d for a duration of τ ql = π/2 √ Nl|g d |. Following this pulse sequence, the state evolves as |0, ↓ → |0, ↑ → |1, ↓ · · · → |N, ↓ , as shown in Fig. 2b. The total duration of this pulse sequence is τ d2 = ∑ l τ cl + τ ql . Assuming that the magnitudes of ε and g d are the fixed for all l's, we find τ d2 = Nπ/2|ε| + ∑ l π/2 √ Nl|g d |, which increases with the total number of polaritons as τ d2 = O(N). Meanwhile, it requires that |ε|, √ N|g d | ≪ ω L to achieve high fidelity for the generated state. Note that the other collective modes a k =0 of the cavities are not coupled to the auxiliary qubit due to the symmetry of the Hamiltonian H d2 . The excitation of these modes will not occur during this pulse sequence.

Optimized nonlinear ramping
Many-body ground states in the intermediate regimes of the parameter space cannot be calculated analytically, and we cannot design quantum logic operations to generate such states, in contrast to the ground states in the deep MI or SF regimes. We employ optimized nonlinear ramping to reach such states via adiabatic evolution. In this approach, a parameter p has the time dependence: where p = g, J, ∆ is a tunable parameter of the JC lattice, p(0) is the initial value of the parameter at time t = 0, p(T ) is the target value at the final time T , and r p is the ramping index of parameter p. For r p = 1, it is the linear ramping studied in [24,25]; and r p = 1 corresponds to nonlinear ramping [41][42][43]. It can be shown that for any parameter p at an arbitrary time t, Hence when the initial and target parameters are given, the ramping trajectory in the parameter space is only determined by the ratios of the ramping indices r g /r J and r ∆ /r J and is independent of the specific value of an individual ramping index (see Supplementary Notes) By varying the ramping index r p , p ′ (p), and hence dH/dt , can be tuned in a large range. We also find that , which reveals that p ′ (p)/J ′ (J) at a given position only depends on the trajectory, i.e., it only depends on the ratios r g /r J and r ∆ /r J that define the trajectory.
Let  [44,45], where the energy gap E gp is defined as the minimal energy separation between the ground and the excited states along the evolution trajectory, and H ′ gp = dH/dt gp denotes the sweeping rate of the Hamiltonian at the position of the energy gap. To reach the desired state with high fidelity, the adiabatic criterion, commonly expressed as H ′ gp ≪ E 2 gp , needs to be satisfied so that diabatic transitions are negligible. For a given trajectory, we can optimize the ramping indices r p to minimize H ′ gp so as to suppress diabatic transitions in the most vulnerable region of the evolution and improve the fidelity of the final state. With the relation between different p ′ (p)'s, we find that H ′ gp = c p p ′ (p gp ) for any parameter p, with p gp the value of the parameter p at the position of the energy gap. Here c p only depends on the ratios r g /r J and r ∆ /r J , not the individual ramping indices; whereas p ′ (p gp ) depends on the ramping index r p as shown in (12). At the optimal ramping index r which only depends on the position of the energy gap for a given trajectory.

Numerical simulation
Below we conduct numerical simulation to calculate the fidelity of the final states via nonlinear ramping with selected trajectories and compare the numerical result with the above estimation. In the simulation, we employ an algorithm developed in our previous work [13], which only involves basis states with the total excitation number N = L. This algorithm greatly speeds up the calculation of the eigenstates and dynamics of JC lattices. We first consider a trajectory following (10) with g(t) ≡ 1, J(0) = 0, J(T ) = 0.5, and ∆(t) ≡ 0, where the photon hopping rate is continuously increased from zero to a finite value. The initial state is the deep MI phase in (5).
Using the exact diagonalization method, we can calculate the eigenstates and eigenenergies of the JC lattice along this trajectory. The energy spectrum of several lowest excited states is plotted as a function of the hopping rate J in Fig. 3a. The solid curve corresponds to the energy of the lowest state that is symmetric with regard to all lattice sites, and the dotted curves are for asymmetric states. As both the initial state and the Hamiltonian H(t) are symmetric with regard to lattice sites, the wave function |ψ(t) at any time t during the evolution must remain symmetric. Hence diabatic transitions can only happen between the ground state and symmetric states, and the energy gap related to the adiabatic criterion is determined by the energy separation between the ground state and the lowest symmetric state. From our numerical result, we find that the gap position is at J gp = 0.122 with the energy gap E gp = 0.31.
The sweeping rate of the Hamiltonian can be written as H ′ gp = c J J ′ gp with J ′ gp the time derivative of the hopping rate J at the gap position. Using (12), we plot J ′ gp as a function of r J in Fig. 3b, where J ′ gp has a local minimum at the optimal ramping index r (min) J = log J(T )/J gp . For the selected trajectory, r (min) J = 1.41, which indicates that the best fidelity for the final state can be achieved with a ramping index in-between the linear and the quadratic forms. At a total evolution time T = 15π/g, the optimal ramping index gives J ′ gp = 0.01. With E gp = 0.31, the adiabatic criterion is well satisfied. We numerically simulate this ramping process and calculate the fidelity of the final state. In Fig. 3c, the fidelity vs J(T ) = J is plotted for several values of r J at T = 15π/g. The fidelity decreases quickly with J(T ), as J ′ gp increases with J(T ). It can also be seen that for J(T ) sufficiently far away from J gp , where the Landau-Zener estima- tion becomes valid, the fidelity is much higher for r J = 1, 2 than that for r J = 1/3, 1/2. As shown in Fig. 3d, the best fidelity for J(T ) = 0.5 is achieved when r J ∈ (1, 2). These results agree well with our derivation of the optimal ramping index.
We also obtain the fidelity of the final states for a wide range of target parameters J(T ), ∆(T ) following the trajectory (10) with g(t) ≡ 1, J(0) = ∆(0) = 0, r ∆ /r J = 1, and the MI initial state (5). The fidelity is presented in Fig. 4a-d for r J = 1/3, 1/2, 1, 2, respectively. It can be seen that the fidelity decreases as the target parameters move further towards the SF phase. In particular, the fidelity exhibits a sharp decrease when the parameters cross the gap positions into the SF phase. Meanwhile, the fidelity demonstrates strong dependence on the ramping index in the intermediate regimes of the parameter space, which also agrees with our analytical prediction.
Next we consider trajectories that starts from the deep SF phase with g(0) = 0, g(T ) = 1, J(0) = 0.5, J(T ) = 0, ∆(t) ≡ 0, and the initial state (6). With both J and g being time-dependent, we can choose different ramping indices for them. As shown in Fig. 5a, the ramping trajectory in the parameter space of J and g depends on the ratio r g /r J , which affects the energy spectrum and the value of the energy gap. In Fig. 5c, we plot the energy spectrum of the lowest excited states vs the hopping rate J for r g /r J = 1, where the solid curve is the energy of the lowest symmetric state. The energy gap occurs at J gp = 0.104 with E gp = 0.25. The energy spectrum for r g /r J = 1 can be found in Supplementary Figure 1a Ramping from SF to MI phase. a Ramping trajectory in the parameter space of J and g for r g /r J = 1/3, 1/2, 1, 2, 3 from top to bottom. b Time derivative J ′ gp vs ramping index r J at r g /r J = 1, J(T ) = 0, and T = 5π/g, 10π/g, 15π/g from top to bottom. c Energy spectrum of the lowest excited states vs hopping rate J for r g /r J = 1. Solid (dotted) curve is for the symmetric (asymmetric) state with the ground-state energy set to zero. d Fidelity F vs r J at r g /r J = 1, J(T ) = 0, and T = 5π/g, 10π/g, 15π/g from bottom to top. e Fidelity F vs J(T ) = J for r J = 2 (solid), 1 (dashed), 1/2 (dot-dashed), and 1/3 (dotted) at r g /r J = 1 and T = 15π/g. f Fidelity F vs r J at r g /r J = 1/3 (square), 1/2 (circle), 1 (diamond), 2 (triangle), and 3 (inverted triangle), J(T ) = 0, and T = 15π/g. In all plots, g(0) = 0, g(T ) = 1, J(0) = 0.5, and ∆(t) ≡ 0. and (11), it can be shown that For a given ratio r g /r J , g ′ gp /J ′ gp is a constant that does not depend on the specific value of r J or r g . The dependence of H ′ gp on the ramping indices can hence be charactered by the dependence of J ′ gp on r J , which is shown in Fig. 5b. For r g /r J = 1, J ′ gp has a minimum at r (min) J = 0.234. In Fig. 5d, we plot the fidelity of the final state vs r J from our numerical simulation, which indicates that the best fidelity can be achieved when r J ∈ (1/2, 1) at T = 10π/g, 15π/g and when r J ∈ (1/3, 1/2) at T = 5π/g. This result confirms our analysis that the optimal ramping index for this trajectory will shift to a smaller value with r  could be owing to the small difference |J gp − J(T )| between the gap position and the target parameter, which affects the accuracy of the Landau-Zener formula in adiabatic processes [44,45]. We also numerically simulate the ramping process for the target hopping rate J(T ) ∈ [0, 0.5] and obtain the fidelity of the final state vs J(T ) for several values of r J , as plotted in Fig. 5e. The fidelity decreases as J(T ) becomes smaller, as |J ′ gp | increases with the difference |J(T ) − J(0)|. The fidelity vs r J for r g /r J = 1 is given in Fig. 5f. It can be seen that the optimal ramping index r (min) J for different r g /r J can be quite different. This is due to the change of the ramping trajectory and the energy spectrum as r g /r J is varied. Detailed results on E gp , J gp , and r (min) J for different values of r g /r J can be found in Supplementary Figure 1b.
We present the fidelity of the final state vs the target parameters following the trajectory (10) with g(0) = 0, g(T ) = 1, J(0) = 0.5, ∆(0) = 0, r g /r J = r ∆ /r J = 1, and the SF initial state (6) for r J = 1/3, 1/2, 1, 2, respectively, in Fig. 6a-d. Our numerical result shows that the fidelity decreases quickly as the target parameters enter the MI phase and strongly depends on the ramping index in the intermediate regimes of the parameter space.

DISCUSSION
We have shown that the fidelity of the prepared state in the intermediate regimes of the parameter space can be improved by choosing the optimal ramping index for a given trajectory and by increasing the total ramping time T . Another approach to increase the fidelity is by choosing a favorable trajectory for a given set of target parameters. When the tar-  Fig. 6a-d, and b between the linear ramping data in Fig. 4c and Fig. 6c. get parameters are in the MI phase, it is better to start from an initial state in the deep MI regime such as (5) so that the adiabatic evolution does not need to cross a region with narrow energy gap to reach the target parameters so that diabatic transitions can be negligible. Similarly, when the target parameters are in the SF phase, we can choose the initial state to be in the deep SF regime such as (6). Combing the selection of the initial state with optimized nonlinear ramping can have dramatic impact on the fidelity of the prepared states. For illustration, in Fig. 7a, we plot the maximal fidelity among all eight sets of data in Fig. 4a-d and Fig. 6a-d for two initial states and various values of linear or nonlinear ramping index r J . It can be seen that the maximal fidelity remains close to unity in almost the entire parameter space. For comparison, in Fig. 7b, we plot the maximal fidelity between the data for linear ramping (r J = 1) in Fig. 4c and Fig. 6c. The result in Fig. 7a outperforms that of Fig. 7b, and both results are much better than the individual plots in Fig. 4 and Fig. 6. We expect that further improvement can be achieved by optimizing the trajectory, e.g., using optimized r g /r J , r ∆ /r J or the optimal control technique.
An obvious approach to improve the fidelity of adiabatic processes is to increase the ramping time T , which can reduce the time derivatives of the parameters and the sweeping rate of the Hamiltonian. This can be seen from the numerical result in Fig. 3d and Fig. 5d. For quantum devices in the NISQ era, however, the decoherence times of qubits and cavity modes set a limitation on the evolution time. The many-body ground states studied here involve finite number of polariton excitations. In the presence of decoherence, excitations can decay in a timescale comparable to the decoherence times. The ramping time needs to be much shorter than the decoherence times. In experiments, superconducting resonator cavities with frequency ω c /2π = 10 GHz and quality factor Q = 10 5 can be readily realized, which corresponds to a decay time of 1.6 µs. Superconducting qubits can have a decoherence time of ∼ 100 µs [51]. With a typical coupling strength of g/2π, J/2π = 200 MHz, the evolution time T = 15π/g ≈ 37.5 ns. The state initialization pulses can be completed within a few 10's of ns. These time scales are much shorter than the decoherence times.
To quantitatively characterize the effect of dissipation, we utilize a phenomenological approach with a non-Hermitian Hamiltonian [55]: H t = H t − i(κ/2) ∑ j a † j a j − i(γ/2) ∑ j σ jz , where κ is the cavity damping rate and γ is the qubit decay rate. We numerically simulate the ramping process under the Hamiltonian H t and calculate the fidelity of the final state. For r J = 1 and T = 15π/g studied in Fig. 3d, with Q = 5 × 10 4 (κ/2π = 200 kHz) and γ/2π = 2 kHz, the fidelity F = 0.9737. For comparison, F = 0.9738 when dissipation is not included. This result shows that the effect of dissipation with practical parameters is negligible at time scales of interest. Details of this result can be found in Supplementary Discussion and Supplementary Figure 2. Note that dissipation can be used for robust preparation and stabilization of entangled states, as studied in [38,39]. In addition, the qubitcavity interaction in (2) has the form of JC coupling with the counter-rotating terms omitted. This is because we study the system in the strong-coupling regime with g ≪ ω c , ω z , where the effect of the counter-rotating terms can be neglected.

METHODS
We employ a method developed in our previous work [13] to conduct numerical simulation on the JC lattice studied in this work. This method allows us to efficiently solve the ground states and the dynamics of a finite-sized JC lattice with given number of polariton excitations. In a JC lattice, each unit cell contains a qubit and a cavity mode. If we choose the photon cutoff on each lattice site to be equal to the total number of excitations in the lattice N, then the total number of basis states for this lattice is (2N + 1) L , which depends exponentially on the size of the lattice L. For N = L = 6, the number of basis states is 4826809. This dependence sets a serious limitation on the computable size of JC lattices. In our method, we only consider basis states that have exactly N excitations. We developed a code to find out such basis states. For N = L = 6, we find that the total number of basis states is 5336, which shows that this method can greatly reduce the demand on computing power. We also developed codes to derive the matrix elements for the Hamiltonian and other operators, such as the creation and annihilation operators of cavity modes, the Pauli operators of qubits, and the number operator of cavity modes, under these basis states. With these matrices, we can calculate the eigenstates and eigenenergies of a JC lattice using the exact diagonalization method. We also simulate the dynamical evolution of this system under time-dependent Hamiltonians. The matrices for the creation operators and the Pauli operators enable us to generate the initial state in the MI and the SF phases numerically. Our method can be applied to different types of JC lattices as it can be used regardless of the specific form of interaction between neighboring lattice sites.

DATA AVAILABILITY
The data that support the findings of this study are available from the authors upon reasonable request. Fig. 1 Quantum phase transition in JC lattice. a Schematic of a 1D JC lattice. Circles (rectangles) represent qubits (cavity modes) with light-matter coupling g and hopping rate J. b Single-particle density matrix ρ 1 (1, 4) vs hopping rate J and detuning ∆ for a finite-sized lattice at unit filling with N = L = 6. Here we let g be the energy unit with g ≡ 1. Fig. 2 Pulse sequence for state initialization. a Pulses for MI initial state at J = 0 and finite g. The vertical arrows are Rabi flips between the states |g 0 and |1, − in each JC model. b Pulses for SF initial state at g = 0 and finite J. The vertical (slanted) arrows are the operations C l (Q l ) with l ∈ [1, N] on the coupled system of the auxiliary qubit and mode a k=0 . Fig. 3 Ramping from MI to SF phase. a Energy spectrum of the lowest excited states vs hopping rate J. Solid (dotted) curve is for the symmetric (asymmetric) state with the ground-state energy set to zero. b Time derivative J ′ gp vs ramping index r J at J(T ) = 0.5 and T = 5π/g, 10π/g, 15π/g from top to bottom. c Fidelity F vs J(T ) = J for r J = 2 (solid), 1 (dashed), 1/2 (dot-dashed), and 1/3 (dotted) at T = 15π/g. d Fidelity F vs r J at J(T ) = 0.5 and T = 5π/g, 10π/g, 15π/g from bottom to top. In all plots, g(t) ≡ 1, J(0) = 0, and ∆(t) ≡ 0.  Ramping from SF to MI phase. a Ramping trajectory in the parameter space of J and g for r g /r J = 1/3, 1/2, 1, 2, 3 from top to bottom. b Time derivative J ′ gp vs ramping index r J at r g /r J = 1, J(T ) = 0, and T = 5π/g, 10π/g, 15π/g from top to bottom. c Energy spectrum of the lowest excited states vs hopping rate J for r g /r J = 1. Solid (dotted) curve is for the symmetric (asymmetric) state with the ground-state energy set to zero. d Fidelity F vs r J at r g /r J = 1, J(T ) = 0, and T = 5π/g, 10π/g, 15π/g from bottom to top. e Fidelity F vs J(T ) = J for r J = 2 (solid), 1 (dashed), 1/2 (dot-dashed), and 1/3 (dotted) at r g /r J = 1 and T = 15π/g. f Fidelity F vs r J at r g /r J = 1/3 (square), 1/2 (circle), 1 (diamond), 2 (triangle), and 3 (inverted triangle), J(T ) = 0, and T = 15π/g. In all plots, g(0) = 0, g(T ) = 1, J(0) = 0.5, and ∆(t) ≡ 0.  6. Fidelity with SF initial state. a-d Fidelity of the prepared state F vs target hopping rate J(T ) and target detuning ∆(T ) for ramping index r J given in the panel. Here g(0) = 0, g(T ) = 1, J(0) = 0.5, ∆(0) = 0, r g /r J = r ∆ /r J = 1, and T = 15π/g.  Fig. 4a-d and Fig. 6a-d, and b between the linear ramping data in Fig. 4c and Fig. 6c.