State preparation in a Jaynes-Cummings lattice with quantum optimal control

High-fidelity preparation of quantum states in an interacting many-body system is often hindered by the lack of knowledge of such states and by limited decoherence times. Here, we study a quantum optimal control (QOC) approach for fast generation of quantum ground states in a finite-sized Jaynes-Cummings lattice with unit filling. Our result shows that the QOC approach can generate quantum many-body states with high fidelity when the evolution time is above a threshold time, and it can significantly outperform the adiabatic approach. We study the dependence of the threshold time on the parameter constraints and the connection of the threshold time with the quantum speed limit. We also show that the QOC approach can be robust against control errors. Our result can lead to advances in the application of the QOC to many-body state preparation.


INTRODUCTION
Recent progresses in manipulating quantum states and dynamics in noisy intermediate-scale quantum (NISQ) devices have demonstrated the potential to solve complicated problems with various platforms [1][2][3].An important question among such problems is the preparation of many-body states with high fidelity using NISQ devices, which is crucial for quantum simulation, quantum metrology and quantum communication [4][5][6][7].In the past, a number of approaches have been developed to generate desired quantum many-body states, including adiabatic processes [8,9], quantum shortcut approach [10], quantum phase estimation [11,12], quantum eigensolvers [13][14][15], and open system approach [16,17].However, due to the intrinsic complexity of quantum many-body systems, it remains challenging to prepare such states with high accuracy.
With quantum control techniques, precisely engineered pulse sequences have been employed to manipulate quantum states with high accuracy [18].Among such techniques, the quantum optimal control (QOC) approach [19][20][21][22] provides a computational framework to generate desired quantum states or quantum processes by searching for optimal, time-dependent control parameters under given constraints.In recent years, the QOC has been widely used in a broad range of applications from the implementation of high-fidelity quantum logic gates, the suppression of environmental noise, the control of quantum transduction processes, the generation of novel entangled states, to the control of quantum many-body systems [23][24][25][26][27].The problem of preparing quantum states or processes can be formulated into an optimization problem in the QOC, where an algorithm is adopted to minimize the cost function.
Here we study the QOC approach for the preparation of many-body states in a finite-sized Jaynes-Cummings (JC) lattice.In the thermodynamic limit, a JC lattice with integer fillings (i.e., the average number of excitations per lattice site is an integer) can exhibit a quantum phase transition between the Mott-insulating (MI) and superfluid (SF) phases [28][29][30][31][32][33].At a finite size, the ground states of a JC lattice in the MI and SF regimes still exhibit distinctive behaviors [34][35][36].The preparation of the ground states in a JC lattice is non-trivial, especially in the intermediate range between the deep MI and deep SF phases.In [2], we employed an optimized adiabatic approach for the state preparation in a JC lattice.In this work with the QOC approach, we adopt the chopped random basis (CRAB) algorithm [38,39] to parametrize the time-dependent couplings of the JC lattice and use the Nelder-Mead approach to optimize these couplings.Our numerical result shows that when the total evolution time is above a threshold time T th , the QOC approach can generate the target state with a high fidelity above a designated threshold value, and it can significantly outperform the adiabatic approach.We find that the threshold time decreases and the average energy fluctuation increases with the constraints on the time-dependent couplings, which indicates the connection between the threshold time and the quantum speed limit (QSL) [40][41][42][43][44][45].Furthermore, our numerical simulation shows that the QOC approach can be robust against control errors in the time-dependent couplings.JC lattices have been explored theoretically and implemented experimentally in various systems, including circuit QED system, nanophotonic devices, atoms, and trapped ions [46][47][48][49][50][51][52][53][54][55][56][57].This work can shed light on the application of the QOC in many-body state preparation and lead to deeper understanding of the QSL in preparing quantum many-body states.

JC Lattice
A JC lattice is illustrated in Fig. 1(a), where each unit cell of the lattice contains a two-level system (qubit) coupled to a cavity mode with coupling strength g, and adjacent cavities are coupled by photon hopping with hopping rate J.The Hamiltonian of this lattice can be written as H t = H 0 + H int ( = 1).Here and all parameters are in dimensionless units [58].
is the Hamiltonian of the JC models in a finite-sized lattice of size N , with j ∈ [1, N ], ω c the cavity frequency, a j (a † j ) the annihilation (creation) operator of the cavity modes, ω z the energy splitting of the qubits, and σ jz , σ j± the Pauli operators of the qubits.Also describes photon hopping between neighboring sites in the lattice.We choose the periodic boundary condition with a N +1 = a 1 and denote ∆ = ω c −ω z as the detuning between cavity and qubit frequencies.The qubit-cavity coupling g induces a built-in nonlinearity in the energy spectrum of a single JC model [59], which can be viewed as an effective onsite interaction with strength U , as shown in Fig. 1(b).Details of the JC model spectrum can be found in the Supplementary Information.In the thermodynamic limit with N → ∞, and at integer fillings when the number of excitations is an integer multiple of N , the competition between this onsite interaction and the photon hopping can lead to a quantum phase transition between the MI and SF phases [28][29][30][31].When dominated by the qubit-cavity coupling with g ≫ J, the ground state of the JC lattice will be in a MI phase characterized by localized polariton excitations.In the limiting case of J = 0, the ground state with N excitations is the product state |G J=0 = N j=1 |1, − j with each JC model in its first excited state |1, − j .When dominated by photon hopping with J ≫ g, the ground state of the lattice will be in a SF phase with long range correlation.In the limiting case of g = 0, the ground state is the Fock state , ↓ with all excitations occupying the momentum-space mode a k=0 = 1 √ N N j=1 a j e ik•j for the quasi-momentum k = 0.For a finite-sized lattice, the ground states also exhibit features of these phases in the corresponding parameter regimes [2,35].These features can be illustrated with the single-particle density matrix ρ 1 (i, j) = G|a † i a j |G / G|a † i a i |G , which describes the spatial correlation between the cavity modes at sites i and j, with |G the ground state for given parameters.As shown in Fig. 1(c), ρ 1 (1, 3) for a N = 4 lattice, and hence, the spatial correlation of the ground state, decrease algebraically (exponentially) in the SF (MI) phase.

Couplings and Fidelity
Preparing the ground states of a JC lattice with integer fillings is a challenging task except for the limiting cases of g = 0 or J = 0.Here we will employ the QOC technique to achieve fast and high-fidelity state preparation in a JC lattice and compare our result with that of the adiabatic approach in [2].We define the fidelity of the prepared state with regards to the desired many-body ground state |ψ T for the target parameters as where |ψ(T ) is the state at the final time T of the evolution.The cost function in the QOC is chosen as the infidelity I = 1 − F. The QOC approach minimizes the cost function by optimizing the coupling constants g (t) and J (t) in the Hamiltonian H t [19][20][21][22].For simplicity of discussion, we let the detuning ∆(t) ≡ 0 during the entire evolution.The couplings are bounded by the constraints g max and J max , with |g(t)| ≤ g max and |J(t)| ≤ J max at an arbitrary time t.The numerical simulation is conducted on a JC lattice with four sites and four polariton excitations (i.e., unit filling).
The initial Hamiltonian parameters are g(0) = 0 and J(0) = 0.5, and the target parameters are g(T ) = 1 and J(T ) = 0.02.The initial state of this system is the ground state for the initial parameters, which is the SF state |G g=0 .The target state is the ground state for the target parameters, which is a MI state.During the evolution, the system is governed by the Hamiltonian H t with time-dependent couplings g(t) and J(t).We adopt the CRAB algorithm that parametrizes the couplings with truncated Fourier series [38,39], and apply the Nelder-Mead method for the optimization.In Fig. 2(a-c), we plot the optimized couplings g(t) and J(t) vs the relative evolution time t/T under the constraints J max = 2 and g max = 1, 2, 4, respectively, with total evolution time T = 3.30π.The couplings in the adiabatic approach governed by (8a) and (8b) are plotted as dashed curves.The optimized couplings are continuous curves that change smoothly over the course of the evolution.For the constraint g max = 1, g(t) includes a large plateau at the maximal strength g(t) = 1; whereas the plateau area decreases significantly for g max = 4.In contrast, J(t) has no plateau.This is because the system can already reach the deep SF phase when J = J max = 1, and it does not require a larger value of J to explore the Hilbert space.Our numerical result also shows that the fidelity for larger g max is significantly higher than that for smaller g max , as shown in Fig. 2(d).For g max = 2 and 4, the fidelity of the state at the final time exceeds the designated threshold fidelity F th = 0.99; while for g max = 1, the fidelity cannot reach 0.99.For all three g max , the fidelity at time T is much higher than that from the adiabatic approach, demonstrating that the QOC approach can greatly outperform the adiabatic approach.
The fidelity of the prepared states depends on the constraints g max , J max , and the total evolution time T .In Fig. 3(a), we plot the fidelity vs T for the constraints g max = 1, 2, 4 and J max = 2.The result shows that the fidelity exhibits an increasing trend with the total time T and the constraint g max .Meanwhile, the fidelity from the QOC approach is significantly higher than the fidelity from the adiabatic ramping.For example, the QOC fidelity is greater than 0.99 for g max = 2, J max = 2 and T = 3.30π, while the fidelity from the adiabatic ramping is only 0.42 for the same evolution time T .In the numerical simulation, we observe that when the total evolution time T is below a threshold time T th , the QOC process cannot achieve a fidelity that is higher than the designated threshold fidelity, which we choose to be F th = 0.99.In Fig. 3(a), the threshold time T th for each set of constraints is indicated by a dashed vertical line.Our result shows that the threshold time decreases as the constraint g max increases.Hence it will take less time to reach a desired fidelity when the coupling g(t) can have a larger magnitude.To analyze the dependence of the threshold time on the constraints, we plot T th vs the constraint J max for different values of g max in Fig. 3(b).It is shown that T th decreases significantly as g max increases, but only decreases slightly when J max increases.For the values of g max used in our simulation, J = J max = 1 is sufficiently large for the system to enter the deep SF phase.Thus the system does not demand larger value of J or subsequently longer evolution time in order to reach high fidelity, which leads to T th 's weak dependence on J max .This result agrees with that of Fig. 2(a-c), where J(t) does not exhibit any plateau during the evolution.The threshold time for different constraints is given in Table .I, together with the fidelity for the total evolution time T = T th from the QOC approach and from adiabatic ramping.

Threshold Time
We compare the threshold time T th from our numerical simulation with an estimation of the quantum speed limit (QSL) T QSL , which is the minimal time for a given quantum system to evolve from an initial state to a target state [40][41][42][43][44][45].We estimate the QSL with [43]: Here, arccos | ψ(0)|ψ T | describes the distance between the initial and the target states in the Hilbert space [45].
For orthogonal states with ψ(0)|ψ T = 0, the distance is π/2.For the initial and target states in our simulation, TABLE I.The threshold time T th for selected constraints and the corresponding fidelity F for total evolution time T = T th using QOC and using adiabatic ramping (adia.).
In Fig. 4(a-c), we plot the energy fluctuation ∆E(t) as a function of the relative evolution time t/T for the same constraints J max , g max , and the same evolution time T as those in Fig. 2. The result from the adiabatic approach is plotted as the dashed curve.In all three plots, the energy fluctuation is the strongest when t/T ∈ (0.3, 0.4), and it is far stronger than the energy fluctuation in the adiabatic process.For g max = 2, 4, ∆E(t) becomes very small when t approaches the final time T , indicating that the final state occupies the ground state with high probability.For g max = 1, ∆E(t) at t = T remains large and is comparable to that from the adiabatic approach, which shows that the system has a sizable probability to be in the excited states in this case.This is because the threshold time for g max = 2, 4 (for g max = 1) is shorter (longer) than the evolution time T = 3.30π, and hence the QOC process can (cannot) reach high fidelity.In Fig. 4(d), we plot the average fluctuation energy ∆E ave vs the total evolution time T for the constraints used in Fig. 4(a-c).Here ∆E ave shows a decreasing trend as T increases and is stronger than that from the adiabatic approach.Using (4) and the result of ∆E ave for the threshold time T th , we estimate the QSL.As shown in the inset of Fig. 4(d), the estimated T QSL exhibits similar behavior to the threshold time T th , decreasing with the increase of the constraint g max .Meanwhile, the estimated QSL is comparable in scale to the threshold time, but it is shorter than the threshold time.We note that this comparison is only qualitative.The estimation of the QSL presented here is a rough approximation due to the complexity of the JC lattice, and the threshold time is defined for a specific threshold fidelity chosen in our numerical simulation.

DISCUSSION Control Error and Decoherence
In superconducting quantum devices, tunable qubitcavity coupling and cavity hopping (i.e., cavity coupling) can reach a few hundreds of MHz [60][61][62][63].For example, tunable qubit-cavity coupling can be achieved via fluxtuned inductive coupling in the g-mon configuration or via a tunable coupler [64,65].Tunable cavity hopping can be achieved by connecting cavities with a tunable Josephson junction [66].We assume that the dimensionless coupling g = 1 used in our numerical simulation corresponds to g = 2π × 100 MHz in actual devices [58].A dimensionless evolution time of T = 3.30π then corresponds to T = 16.5 ns.The optimized, time-dependent couplings g(t) and J(t) need to be generated within this time scale, which can be implemented with current technology.
To explore the robustness of the QOC approach against control errors, we simulate the errors by adding a time-dependent Gaussian noise to the optimized solutions of g(t) and J(t) with where δ 1 (t) and δ 2 (t) are Gaussian noise at time t with the standard deviation σ.We obtain the fidelity of the prepared state in the presence of these errors.For a FIG. 5. Effect of control error.The fidelity vs the standard deviation σ of Gaussian control errors.The total evolution time for given constraints is the corresponding threshold time T th in Table I.
given value of σ, we conduct the simulation on 1000 samples of the time-dependent errors and calculate the average value of the fidelity.The total evolution time is chosen to be the threshold time T th for given constraints.Fig. 5 shows that the fidelity of the prepared state decreases with the standard deviation of the control errors.We observe that the decrease of the fidelity for larger g max is slower than that for smaller g max ; whereas the decrease of the fidelity remains almost the same for different values of J max .For J max = 2 and σ = 0.05, the fidelity is reduced to F = 0.9902, 0.9910, 0.9948 for g max = 1, 2, 4, respectively.Compared with the fidelity for no control errors given in Table I, the reduction of the fidelity is negligible.This result shows that the QOC approach can be robust against control error.Another factor that could affect the fidelity of the QOC process is the decoherence of the qubits and the cavity modes.In the NISQ era, the finite decoherence time of the quantum devices sets a limit on the time for coherent evolution.In a JC lattice with unit filling, the polariton excitations can decay in a time scale comparable to the decoherence times.At the current state-of-the-art, the decoherence time of superconducting qubits can reach ∼ 100µs [60][61][62].Superconducting cavities can have quality factor greater than 10 6 .For a cavity frequency of ω c = 2π × 5 GHz, such a quality factor gives a cavity decay time of ∼ 32µs.With an evolution time of 16.5 ns, the QOC can be completed in a much shorter time than the decoherence time of superconducting qubits and cavities, and hence the effect of decoherence can be neglected.This analysis has been confirmed by our numerical simulation using a mas-ter equation approach, as detailed in the Supplementary Information.We want to note that the short evolution time required in the QOC approach is one of its advantages over the adiabatic approach.

METHODS
In the CRAB algorithm used in our simulation, the timedependent parameters g(t) and J(t) are parametrized with truncated Fourier series to the 8th harmonics and can be written as [38,39] with where c i,k and d i,k (i = 1, 2 and k ∈ [1,8]) are the Fourier coefficients of the k-th harmonics in g(t) and J(t), respectively, and ω i,k = k + δω i,k is the frequency of the k-th harmonics with an adjustable offset δω i,k .
For a given set of initial (target) parameters for the JC lattice, we obtain the initial state |ψ 0 (the target state |ψ T ) by diagonalizing the corresponding Hamiltonian H t .The optimization process begins with a random set of parameters c i,k , d i,k , and δω i,k and has a maximum of 150000 iterations.The convergence of these parameters vs the iteration number n can be found in Fig. S1 of the Supplementary Information.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

SUPPLEMENTARY NOTES Energy Spectrum of JC Model
We consider a single Jaynes-Cummings (JC) model on site j ∈ [1, N ] of the lattice with the Hamiltonian with ω c the frequency of the cavity mode, a j (a † j ) the annihilation (creation) operator of the cavity mode, ω z the energy splitting of the qubit, and σ jz , σ j± the Pauli operators of the qubit.We describe the eigenstates of the JC model using the basis set {|n, s }, which are the product states of the photon number state |n (n being integer) of the cavity and the spin up or down state (s =↑ or ↓) of the qubit.
The eigenstates of the JC model include the ground state |g 0 = |0, ↓ with zero photon and the qubit in the spin down state (i.e., no excitation in the model), and the polariton doublets |n, ± with n excitations: where , and ∆ = ω c − ω z is the detuning between the cavity and the qubit [S1].
The eigenenergies for different eigenstates are E g0 = 0, and When the coupling strength g is nonzero, the energies of the eigenstates are not equally spaced with, e.g., This nonlinearity in the energy spectrum is at the root of many interesting phenomena in the JC model (and hence JC lattices), such as the photon blockade effect and the Mott insulator-superfluid phase transition.

QOC Algorithm and Optimization Parameters
In our numerical simulation, we first diagonalize the Hamiltonian of the JC lattice at the initial parameters g(0), J(0) and ∆(0) to find its ground state |ψ 0 .This state is used as the initial state of the evolution with |ψ(0) = |ψ 0 .Here, the initial parameters are chosen so that |ψ 0 is easy to prepare in realistic systems [S2].We then solve the ground state of the Hamiltonian at the target parameters g(T ), J(T ) and ∆(T ), which is the desired many-body state |ψ T to be prepared at the final time T of the evolution.
In our quantum optimal control (QOC) approach, we adopt the chopped random basis (CRAB) algorithm to parameterize the couplings g(t) and J(t) with truncated Fourier series.The parameters c i,k , d i,k and δω i,k (i = 1, 2 and k ∈ [1, 8] being integers), as defined in Eq. (7a) and (7b) in the Methods section of the main paper, include a total of 48 optimization parameters.In the beginning of the optimization process, a random set of c i,k , d i,k and δω i,k are used to generate the initial trial functions for g(t) and J(t).Note that the choice of the random parameters only affects the values of g(t) and J(t) at 0 < t < T .At t = 0 (t = T ), the couplings always remain the same as the initial (target) parameters of the JC lattice.During the optimization process, following the protocol in Eq. ( 9) in the Methods section of the main paper, the values of g(t) and J(t) are bounded by the constraints g max and J max , respectively.Given a set of time-dependent couplings g(t) and J(t), the system evolves from the initial state |ψ 0 under the Hamiltonian H t [g(t), J(t)] to reach the final state |ψ(T ) at time T .The cost function for this set of couplings can then be calculated.In this work, we use the Nelder-Mead method to minimize the cost function and find the optimal values for c i,k , d i,k and δω i,k .
In our simulation, we choose the maximal number of iterations to be 150000.For a given total time T , if the optimization process converges to a fidelity above the threshold value F th = 0.99 within the maximal number of iterations, we consider the QOC process successful.The convergence of the QOC process is determined by the difference between the optimization parameters in adjacent iterations.If the difference is smaller than the pre-defined error tolerance, then we consider that convergence has been achieved.In Fig. S1(a-f), we plot the Fourier coefficients c i,k , d i,k and the frequency offsets δω i,k vs the iteration number n under the constraints J max = 1, g max = 2 and for the total evolution time T = T th = 3.30π.From the numerical result, we find that the number of iterations for this simulation to converge is n = 43194.

Master Equation Approach for Decoherence
To quantitatively characterize the effect of decoherence, we adopt the following master equation [S3]: where ρ is the density matrix of the JC lattice, H t is the total Hamiltonian with the time-dependent, optimized

FIG. 1 .
FIG. 1. Jaynes-Cummings Lattice.a The schematic of a one-dimensional JC lattice with qubit-cavity coupling g and photon hopping rate J. b The energy spectrum of a single JC model for detuning ∆ = 0, with |0, ↓ the ground state and |n, ± (n > 1, integer) the lowest excited states.c The single particle density matrix ρ1(1, 3) vs hopping rate J and detuning ∆ for a N = 4 lattice at unit filling.Here g = 1,and all parameters are in dimensionless units[58].

FIG. 2 .
FIG. 2. Optimized couplings.a-cThe optimized couplings g(t) and J(t) vs the relative evolution time t/T for the constraints Jmax = 2 and gmax = 1, 2, 4, respectively, with the total evolution time T = 3.30π.d The fidelity F vs t/T for the couplings in a-c.The dashed curves are for the adiabatic ramping.

FIG. 3 .
FIG. 3. Fidelity and threshold time.a The fidelity F of the prepared state vs the total evolution time T .The vertical dashed lines indicate the position of the threshold time T th for each set of constraints.The constraints are Jmax = 2 and gmax = 1, 2, 4. The dashed horizontal line corresponds to the threshold fidelity F th = 0.99.b The threshold time T th vs the constraint Jmax for gmax = 1, 2, 4.

FIG. 4 .
FIG. 4. Energy fluctuation.a-c The energy fluctuation ∆E(t) vs the relative evolution time t/T for the constraints Jmax = 2 and gmax = 1, 2, 4, respectively, and T = 3.30π.The dashed curve is for the adiabatic ramping.d The average energy fluctuation ∆Eave vs the total evolution time T for the constraints in (a-c).The inset of (d) shows the threshold time T th (circle) and the estimated TQSL (triangle) vs gmax.

( a )
FIG. S1.Optimization parameters c i,k , d i,k , and δω i,k vs the iteration number n.The constraints for the couplings are Jmax = 1 and gmax = 2.The evolution time is T = T th = 3.30π.