Controlling qubit-oscillator systems using linear parameter sweeps

We investigate the dynamics of a qubit-oscillator system under the influence of a linear sweep of system parameters. We consider two main cases. In the first case, we consider sweeping the parameters between the regime of a weakly correlated ground state and the regime of a strongly correlated ground state, a situation that can be viewed as a finite-duration quench between two phases of matter: the normal phase and the superradiant phase. Excitations are created as a result of this quench. We investigate the dependence of the excitation probabilities on the various parameters. We find a qualitative asymmetry in the dynamics between the cases of a normal-to-superradiant and superradiant-to-normal quench. The second case of parameter sweeps that we investigate is the problem of a Landau-Zener sweep in the qubit bias term for a qubit coupled to a harmonic oscillator. We analyze a theoretical formula based on the assumption that the dynamics can be decomposed into a sequence of independent Landau-Zener transitions. In addition to establishing the conditions of validity for the theoretical formula, we find that under suitable conditions, deterministic and robust multi-photon state preparation is possible in this system.


I. INTRODUCTION
The study of cavity quantum electrodynamics (QED) has played a major role in shaping our understanding of various quantum phenomena and has led to the development of various new technologies [1][2][3]. In recent years the remarkable advances made with cavity QED systems utilizing superconducting circuits (now known as circuit QED) have accelerated the experimental realization of phenomena that had been predicted theoretically decades earlier but could not be observed in conventional cavity QED setups [4]. The development of circuit QED has also spurred the emergence of new ideas and applications of cavity QED systems.
Some of the advantages of superconducting circuits in circuit QED systems are (1) the ability to reach previously inaccessible parameter regimes (e.g. deep strong coupling), (2) the ability to set the circuit parameters by appropriate circuit design and (3) the ability to design circuits in which the parameters can be tuned in real time.
The tunability of circuit parameters is often utilized to choose optimal bias points and apply sinusoidal signals to drive transitions between different energy levels to manipulate the quantum state of the system. The tunability can also be utilized by performing parameter sweeps or quenches. While sinusoidal driving is studied more extensively, linear parameter sweeps are also important in quantum computing and related applications. For example, parameter sweeps can be used to perform quantum annealing and adiabatic quantum computing [5,6]. Understanding the quantum dynamics under such unidirectional parameter sweeps is crucial to designing good control protocols that utilize the parameter sweeps for efficient state preparation. One-way parameter variations in quantum systems can also lead to interesting phenomena such as Landau-Zener (LZ) transitions [7] and the creation of topological defects via the Kibble-Zurek mechanism [8]. Having a good grasp on the probability that the system might leave the desired quantum state or the amount of energy generated as a result of the parameter sweep is crucial for practical applications. We will therefore investigate the manifestation of these effects in the present study.
One of the important tasks in circuit QED systems is the preparation of nonclassical states, such as Fock states and Schrödinger cat states. These could be useful for various applications, including quantum communication, sensing, and error correction. We will therefore pay special attention to the possibility of obtaining practically interesting states at the end of the parameter sweep. It is worth noting in this context that there are alternative proposals in the literature for generating nonclassical states in cavity QED systems using periodic modulation of the qubit-cavity coupling strength [9,10]. One potential advantage of linear sweeps is that they can be robust against certain fluctuations, depending on the details of the experimental setup under study.
In this work, we perform a systematic analysis of some practically relevant situations in which the parameters of a cavity QED system are varied in a linear sweep. We inspect the results to gain insight into the quantum system dynamics and to search for potential applications such as controlled state preparation. Two cases arise naturally when considering circuit QED systems with variable parameters. In the case where the qubit is biased at its symmetry point, the system possesses normal and superradiant phases, depending on the relation between the various system parameters. The normal phase can be thought of as the weak-coupling phase and is associated with weak qubit-cavity correlations in the low-lying energy eigenstates, while the superradiant phase is associated with strong coupling and strong qubit-cavity correlations in the low-lying energy eigenstates. We investigate the evolution of the system following a quench between the two different phases. We find in particular that the number and energy of the excitations created in the quench is qualitatively asymmetric between the cases of a normal-to-superradiant quench and a superradiant-tonormal quench. The other case in which a linear parameter sweep arises naturally in a circuit QED system is related to the LZ problem. In addition to causing a nonadiabatic transition in the qubit, the sweep can lead to the creation of photons in the cavity. We investigate the dynamics in this case and calculate the probabilities of various final outcomes. In particular, we show that a standard LZ sweep protocol can be used to almost deterministically prepare Fock states in the cavity.
The remainder of this manuscript is organized as follows: in Sec. II, we introduce the Rabi model that describes the cavity QED systems of interest for the present study. In Sec. III, we describe our simulations and results for the problem of sweeping the system parameters across the normal-superradiant phase boundary. In Sec. IV, we consider the case of performing a linear LZ parameter sweep in a qubit coupled to a harmonic oscillator.
In Sec. V we discuss possible implementations using superconducting circuits. We summarize our conclusions in Sec. VI.

II. RABI MODEL
The quantum Rabi model (QRM) describes a qubit coupled to a harmonic oscillator [11,12]. The QRM Hamiltonian can be expressed aŝ where ∆ is the qubit gap, ϵ is the qubit bias parameter, ω is the oscillator frequency (which implies the conventionh = 1), and g is the qubit-oscillator coupling strength. The operatorŝ σ β (with β = x, y, z) are the qubit's Pauli operators, whileâ andâ † are, respectively, the oscillator's annihilation and creation operators.
Previous studies on the QRM generally focused on the symmetric case ϵ = 0, mainly because the physical systems that were modeled by the QRM in the past were typically most accurately described by the case ϵ = 0. In fact, it is common in the literature not to include the bias term in the QRM Hamiltonian. Recent technological advances, e.g. with circuit QED systems using superconducting circuits, provided setups where ϵ can be finite and is usually easily tunable [13]. We shall start by setting ϵ = 0 in the remainder of this section and in Sec. III but consider the case of varying ϵ in Sec. IV.
One reason why the case ϵ = 0 is of particular interest for the present study is that this case can exhibit the so-called superradiance phase transition. When 4g 2 ≪ ω∆ the system is in the normal phase, where the low-lying states are approximately given by the qubit states are defined byσ z |↑⟩ = |↑⟩ andσ z |↓⟩ = − |↓⟩, such thatσ x |→⟩ = |→⟩ and σ x |←⟩ = − |←⟩ with |→⟩ = (|↑⟩ + |↓⟩) / √ 2 and |←⟩ = (|↑⟩ − |↓⟩) / √ 2, and the quantum number n is the number of excitations (or photons) in the oscillator (n =â †â ). The ground state is |→, 0⟩. When 4g 2 ≫ ω∆ the system is in the superradiant phase, where the low-lying states are approximately given by the displacement operatorD(α) is defined asD(α) = exp αâ † − α * â , and α = g/ω. The ground state is |+, 0⟩. In the semiclassical limit, i.e. the limit in which treating the dynamical variables in the Hamiltonian as classical dynamical variables leads to good physical predictions, an abrupt change in behaviour (e.g. in qubit-oscillator correlations) occurs at a certain point if one of the system parameters is varied. The semiclassical limit is typically associated with the Dicke model, in which a large number of qubits are coupled to the oscillator. In the case of a single qubit, a similar semiclassical limit is realized when ω/∆ → 0 [14][15][16]. In this case, an abrupt change occurs at the point 4g 2 = ω∆, with qubit-oscillator correlations remaining negligibly small below this critical coupling strength but increasing to finite values above the critical point.
It will be useful for our analysis below to note that the energy eigenstates in Eqs. (2) and (3) are either symmetric or antisymmetric with respect to the parity operatorP = exp iπ â †â +σ +σ− , whereσ ± = (σ y ± iσ z ) /2. (Note that this definition ofσ ± is different from the standard one because the qubit term in our symmetric QRM Hamiltonian is proportional toσ x and notσ z , which is more commonly used in the cavity-QED literature.) The Hilbert space can therefore be divided into parity-symmetric and parity-antisymmertic subspaces that are completely decoupled from each other. The parity-symmetric subspace con-

III. QUENCHES BETWEEN NORMAL AND SUPERRADIANT PHASES
One interesting situation that could be realized in a circuit QED system is a finiteduration quench in which one of the system parameters is swept through the critical point that separates the normal and superradiant phases. Since we normally think of correlations as resulting from the coupling term in the Hamiltonian, the most natural scenario to think of in this context is perhaps sweeping the coupling strength g, i.e. setting g = vt where v is the sweep rate and t is the time variable with initial value t = 0 and final value t = T . The Hamiltonian can then be expressed aŝ At t = 0 the qubit and oscillator are completely decoupled (g = 0), meaning that the system starts in the normal phase. Provided that vT ≫ √ ∆ω/2, the system parameters correspond to the superradiant phase at the final time.
One complication with this scenario is that the energy eigenstates do not have Tindependent asymptotic values, since α in Eq. (3) would increase indefinitely with t if g did so. In other words, physical quantities such as the number of photons in the system diverge as t → ∞. We can in principle analyze quantities that converge to finite asymptotic values when g → ∞. However, these quantities do not arise intuitively. We therefore take a different approach.
We fix ω and g, and we vary ∆ to achieve the sweep between the normal and superradiant phases. It is worth mentioning here that tunable-∆ superconducting qubits have been demonstrated in a number of recent experiments [17,18]. We first consider a sweep from a large value (∆ i ≫ 4g 2 /ω; normal phase) at the initial time to a small value (∆ f ≪ 4g 2 /ω; superradiant phase) at the final time: with the final time given by In the QRM, the asymptotic expressions for the lowest energy eigenstates in the limits ∆ → ∞ and ∆ → 0 are, to lowest order, independent of the exact values of ∆. As a result, provided that ∆ i ≫ 4g 2 /ω and ∆ f ≪ 4g 2 /ω, the results can be expected to be almost independent of the exact values of ∆ i and ∆ f . As long as g/ω is finite, no physical divergences occur in the final state of the system.
It is worth noting here that a related problem was considered in Ref. [19]. In that study, the authors took the initial and final values of ∆ to be ±∞ and demonstrated that interesting states, e.g. a single-photon state, can be generated at the end of the parameter sweep. However, it is difficult in realistic systems, e.g. using superconducting circuits, to tune ∆ through zero between positive and negative values [17,18]. We therefore focus on the case where ∆ does not change sign. A related experiment realized a quench between dynamically engineered effective quantum Rabi models with weak and strong interactions in a suitably rotating frame [20].
It is also worth noting that the semiclassical limit, where there is a sharp boundary separating the normal and superradiant phases, is realized when ∆/ω > ∼ 10 3 (see Ref. [15]). If we vary ∆ between large and small values, the semiclassical limit will not be valid throughout the sweep. However, the key point here is that at the semiclassical phase transition point the system parameters obey the relation ∆/ω = (2g/ω) 2 . As a result, if g/ω is significantly larger than 1, then ∆/ω will also be much larger than 1 when the point ∆/ω = (2g/ω) 2 is crossed. In this case, the semiclassical picture is valid at the most crucial time during the parameter sweep, and an abrupt change is expected at this point during the sweep.
We take the initial state of the system to be the ground state, which is approximately given by |→, 0⟩ in Eq. (2). This choice is motivated mainly by its experimental relevance, as the ground state is generally the easiest state to prepare. Since the ground state is symmetric with respect to the parity operator and we are considering the symmetric QRM (ϵ = 0), the parity is a conserved quantity, and the system remains in the parity-symmetric subspace throughout its evolution. We therefore confine our analysis below to this subspace of the Hilbert space, e.g. when we refer to the first-or second-excited state, ignoring the parity-antisymmetric subspace of the Hilbert space.
In the adiabatic limit v → 0, the system remains in the ground state of the instantaneous Hamiltonian, and hence ends up in the ground state of the Hamiltonian at the final time.
In the fast-sweep limit v → ∞, the total sweep time is infinitesimally short, and the system does not have time to experience any dynamical evolution. As a result, the final state probabilities are given by the overlaps between the initial state and final energy eigenstates.
If we take the limits ∆ i → ∞ and ∆ f = 0, the initial and final energy eigenstates are given exactly by Eqs. (2) and (3), respectively. Hence, we can use the formula for the excitationnumber probability distribution in a coherent state of a harmonic oscillator to obtain exact expressions for all the probabilities in the final state: where n stands for the nth excited state (with n = 0 for the ground state), and ⟨n⟩ = (g/ω) 2 .
The function P (n) peaks at n = ⟨n⟩.
To find the probabilities at intermediate values of v, we performed numerical simulations in which we solve the time-dependent Schrödinger equation. In the simulations, we generally set ∆ i /ω = 4×10 3 and ∆ f /ω = 0. We include up to 500 states in the truncated Hilbert space.
We divide the total time into 10 5 time steps. To ensure that our results are not affected by the finiteness of these parameters, we vary ∆ i /ω, the size of the truncated Hilbert space and the number of time steps to verify that the results are essentially unchanged, as long as these parameters remain sufficiently large. The occupation probabilities P (±, n) of a few representative final states (in the paritysymmetric part of the Hilbert space) as functions of the sweep rate (in the dimensionless combination v/ω 2 ) for a few different values of g/ω are plotted in Fig. 1. As expected, the system remains in the ground state in the adiabatic limit in all cases. The probabilities change as the sweep rate increases until they reach the asymptotic values given in Eq. (6) in the fast-sweep limit. For relatively small values of g/ω (Fig. 1a), all the probabilities P (±, n) evolve monotonically between the two limits v/ω 2 → 0 and v/ω 2 → ∞. First P (−, 1) has a peak, then P (+, 2) and so on, until n reaches ⟨n⟩. The probabilities P (±, n) for n > ∼ ⟨n⟩ increase monotonically with increasing v/ω 2 and reach their asymptotic values at v/ω 2 → ∞.
One overall feature that we can see in Fig. 1 is that the maximum values of P (±, n), i.e. the heights of the P (±, n) peaks, generally decrease with increasing n. Although it is not clear in the plots shown in Fig. 1, this trend is not monotonic: after a clear decrease for small values of n, the peak heights stabilize around a certain value (with a weak oscillatory behaviour). For n > ⟨n⟩, where there are no peaks anymore, the maximum values of P (±, n) [at v/ω 2 → ∞] decrease with increasing n until they vanish in the limit n → ∞.
There are a few additional small but interesting features in the behaviour of P (±, n).
There is a weak wavy behaviour in the probabilities P (±, n) as functions of v/ω 2 . The wavy behaviour is more clearly visible at higher values of g/ω. By comparing the lines for different values of g/ω (not shown in a single plot in the figure), we find that in the adiabatic regime (around v/ω 2 ∼ 1) P (+, 0) exhibits non-monotonic behaviour as a function of g/ω.
Considering the g/ω values used in the different panels in Fig. 1, Panel (c) has ∆/ω = (2g/ω) 2 = 100 and Panel (d) has ∆/ω = (2g/ω) 2 = 1600 at the semiclassical crossing point, which are therefore expected to be at least close to the semiclassical regime [15]. However, the probabilities P (±, n) do not exhibit any drastically different behaviour that indicates the crossing of a phase transition point, except for the fact that a large number of photons are generated by a fast sweep. This result is rather surprising, since one would expect a qualitatively different behaviour associated with the crossing of a sharp phase transition boundary, compared to cases where no such sharp boundary exists.
To gain further insight into the dynamics, we plot the probabilities P (±, n) as functions of time t for the fast-sweep regime in Fig. 2. In all cases, the system remains mostly in the ground state until it reaches the semiclassical crossing point, which is given by v(t − T )/ω = −(2g/ω) 2 . Then state mixing starts and continues until the final time. In the case of relatively small g/ω (Fig. 2a), one interesting feature is that somewhat highly excited states are temporarily populated at intermediate times, even though only low-lying states are populated at the final time. Since we are dealing with a dynamic situation, and the direct nonadiabatic transitions to high energy levels can be stronger than those between neighbouring energy levels, this situation is not particularly surprising. For stronger coupling ( Fig. 2b-d) a peak in P (−, 1) appears first, followed by a peak in P (+, 2) and so on, until the final time. As g/ω increases, the drop in P (+, 0) becomes increasingly sudden. This behaviour is consistent with what is expected in the semiclassical limit: when the transition point is crossed, excitations are created in analogy with the Kibble-Zurek mechanism. As can be seen most clearly in Fig. 1d, the fast drop in P (+, 0) at v(t − T )/ω = −(2g/ω) 2 is accompanied by a fast rise in probabilities P (±, n) with small values of n. The state of the system then keeps moving towards higher values of n until the final time. In other words, photons keep being added to the system until the final time. The fact that the semiclassicalpicture dynamics occurs only for (2g/ω) 2 > ∼ 10 3 is consistent with the results of Ref. [15], where it was shown that the sharp boundary between the normal and superradiant phases requires that ∆/ω > ∼ 10 3 . One implication of the rapid succession of peaks in the case of large g/ω (Fig. 2d) is that the final state is highly sensitive to the value of ∆ f . In other words, if ∆ f is small but finite, the final state will have a large number of photons but will have small overlap with the final state for ∆ f = 0.
Since we are interested in the question of controlling the system via linear parameter sweeps, we look at the possibility of preparing specific states by a proper choice of the sweep rate. Apart from P (+, 0), which approaches 1 in the adiabatic limit, no other probability reaches a large value (all below 0.4 in Fig. 1). For example, the highest possible value of P (−, 1), which is given by 0.368, is obtained in the fast-sweep limit (v/ω 2 → ∞) at g/ω = 1.
We take the ground state in the parity-symmetric sector of the Hilbert space as the initial state. The results of the corresponding simulations are shown in Figs. 3 and 4. In these figures, the plotted quantities are the populations P (γ, n) where γ is → or ←. In the adiabatic limit, the system remains in its ground state, as expected. Away from the adiabatic limit, various excited states are populated with varying probabilities. In the fast-sweep limit v/ω 2 → ∞, the probabilities are given by Eq. (6), determined by the number of photons in the final state, independently of the qubit state. Since the qubit gap is the dominant energy scale at the final time, the system occupies highly excited states (corresponding to the qubit's excited state |←⟩) with about 50% probability. The population of these highly excited states increases with increasing g/ω. In the limit g/ω ≫ 1 and v/ω 2 → ∞, the qubit ground-state and excited-state sectors are equally populated. Another feature that is different when we sweep from the strongly-correlated to the weakly-correlated regime is that the occupation probabilities of individual energy eigenstates reach higher values (which one can see by looking at the peaks of P (γ, n) plotted as functions of v/ω 2 ), especially in the case of large g/ω (Fig. 3c,d). From the point of view of state preparation, this situation allows a more deterministic preparation of specific energy eigenstates at the end of the sweep.
Examples of the dynamics of the probabilities P (γ, n) as functions of time t for the fast-sweep regime are shown in Fig. 4. For relatively small values of g/ω (Fig. 4a,b), the probability P (→, 0) drops from its initial to its final value around the time vt/ω = (2g/ω) 2 , similarly to what we saw in Fig. 2. For larger values of g/ω (Fig. 4c,d), P (→, 0) drops to its final value (which is essentially zero) well before the semiclassical crossing point is reached.
By the time that the semiclassical crossing point is reached, the system is already in a highly excited state. For all values of g/ω, the different probabilities P (γ, n) for low-lying final states exhibit transient peaks. However, unlike the behaviour shown in Fig. 2, the probabilities P (γ, n) of higher energy levels rise monotonically and asymptotically approach their final values. As a result, the final state probabilities are relatively insensitive to the exact value of ∆ f , in contrast to the dynamics shown in Fig. 2d. As in the case of P (γ, n) as functions of v/ω 2 , the overall shapes of the curves is quite different when comparing the weak-to-strong and strong-to-weak-correlation sweeps.
One might wonder about the physical origin of the asymmetry between the two sweep directions. Since we take the initial state to be the ground state of the Hamiltonian at the initial time, in the weak-to-strong-correlation sweep, we effectively start with the ground state of a single well that evolves into a double well, while in the strong-to-weak-correlation sweep we start with the ground state of a double well that evolves into a single well. There is no reason to expect any symmetry between the two cases, except for the ground state probability. One crucial factor in the asymmetry between the two sweep directions pertains to the energy level order in the two limits ∆ → ∞ and ∆ → 0. In the limit ∆ → ∞, the lowest energy levels correspond to the states |→, 0⟩, |→, 2⟩, |→, 4⟩, ..., with the states |←, 1⟩, |←, 3⟩, |←, 5⟩, ... having much higher energies. In contrast, when ∆ → 0, the energy levels in increasing order correspond to the states |+, 0⟩, |−, 1⟩, |+, 2⟩, |−, 3⟩, ...

In relation to the normal-to-superradiant phase transition in the quantum Rabi model, it
is worth mentioning the dynamical phase transition that occurs when the qubit or oscillator is driven on or near resonance [21]. At a certain driving amplitude the quasienergy spectrum collapses, and above this threshold amplitude no normalizable steady states exist. The phase transition in our case is different, because the spectrum collapses only at the transition point but simplifies away from the critical point on both sides of the transition. As such, we cannot use our results to infer the behaviour of a driven qubit-oscillator system as the one studied in Ref. [21] with a linear sweep of parameters. It will be interesting to investigate a quench across the critical point in that situation in the future.

IV. LANDAU-ZENER SWEEP OF QUBIT COUPLED TO OSCILLATOR
We now turn to the case of finite ϵ in Eq. (1), in particular sweeping ϵ from a large negative value to a large positive value: This situation is in fact close to models studied previously in the literature. The first two terms in Eq. (7) describe the LZ problem. The last two terms can be viewed from the perspective of the spin-boson model, in which the environment of a quantum two-level system is modeled as an infinite set of harmonic oscillators. As a result, combining the LZ problem with the spin-boson model to describe the LZ problem in a two-level system that is coupled to an environment gives a generalized version of Eq. (7) in which there are an infinite number of harmonic oscillators [22][23][24][25]. The model described by Eq. (7), i.e. with a singlemode environment, was also used to study the effect of a finite-temperature environment on the LZ problem [26,27]. These past studies focused on the effect of the harmonic oscillators on the dynamics of the two-level system. Here we are equally interested in the state of the harmonic oscillator, which we view as part of the accessible and controllable quantum system, as opposed to being an unwanted environment that disrupts the dynamics of the controllable quantum system.
Instead of using the basis of bare states, i.e. |γ, n⟩ where the qubit state index γ is ↑ or ↓ and the photon number operatorn =â †â , it is physically more meaningful to use a correlated basis in which the photon numbern is defined asn ↑ = (â † + g/ω)(â + g/ω) for the qubit state ↑ andn ↓ = (â † − g/ω)(â − g/ω) for the qubit state ↓ [13,28]. In other words, the states of the oscillator are displaced to account for the effective field induced by its interaction with the qubit. We use this modified basis below. It is worth noting that the displacement of the oscillator variables has the same physical origin as the displacement in the energy eigenstates in Eq. (3).
To have an intuitive picture of the setup, we consider the case ∆ = 0. This case leads to a simple energy level structure. As t goes from −∞ to +∞, energy level crossings occur at the points vt = mω with m being any integer. At each value of m, the energy levels that correspond to the states |↓, n⟩ and |↑, n + m⟩ intersect, with n = 0, 1, 2, · · ·. In other words, the energy level structure looks like a mesh containing an infinite number of energy level crossings. Taking a nonzero value of ∆ turns all the energy level crossings into avoided crossings with gaps determined by the various system parameters. Figure 5 shows a schematic diagram of the energy level structure, keeping only a few states that have the qubit in its state ↓ and a few states that have the qubit in its state ↑.
If the initial state at t → −∞ is the ground state, i.e. |↓, 0⟩, the above intuitive picture gives a very good approximation for the final energy eigenstate probabilities, as long as ω is not much smaller than ∆ [29][30][31]. In other words, the probabilities P (γ, n) are, to a very good approximation, given by the approximation that the different state populations are determined by a sequence of independent LZ transitions that occur as the system goes through the avoided crossings one by one: . . .
This expression for ∆ n resembles a Poisson distribution with ⟨n⟩ = (2g/ω) 2 , because it is obtained by taking the overlap between the oscillator state |0⟩ displaced by g/ω in one direction with the state |n⟩ displaced by g/ω in the opposite direction [13]. It is worth noting here that the expressions for P (↓, n) are exact [24]. Note also that the relation n ∆ 2 n = ∆ 2 follows from the completeness of the photon-number basis states, and it ensures that the sum of all the probabilities in Eq. (8) is equal to one.
It can be intuitively expected that if the gaps at the different avoided crossings, i.e. Eq. (9), are small relative to the energy scale separating the energy levels, i.e. ω, the picture of independent LZ transitions will be valid. This situation is realized when ∆/ω < ∼ 1 and/or g/ω ≫ 1. Outside these regimes, the energy level structure becomes more complex, and the LZ transitions cannot be considered as occurring independently. We will show the results of numerical simulations that test the validity conditions of Eq. (8). We note here that the extreme regime where ∆/ω ≫ 1 and g/ω ≫ 1 is not of strong interest to us in this work, partly because both conditions are difficult to realize experimentally, in addition to the fact that the corresponding simulations are computationally challenging. We therefore do not investigate this regime in detail in our numerical simulations.
According to Eq. (8), the probabilities P (↑, n) depend on three independent parameters: v/∆ 2 , g/ω and n. In particular, the formulae suggest that the ratio ∆/ω does not affect the structure of the probability peaks described by Eq. (8). The qubit gap ∆ affects P (↑, n) only through the ratio ∆ 2 /v. For example, an increase in ∆ simply shifts all the P (γ, n) curves to higher values of v. The formulae in Eq. (8) for a few of the representative states are plotted in Fig. 6. As mentioned above, the probability for the qubit to remain in the initial state, i.e. |↓, 0⟩, is independent of the coupling to the oscillator. For the qubit state ↑, as one would intuitively expect, the occupation probability of the final-time ground state (i.e. |↑, 0⟩) approaches one in the adiabatic limit and approaches zero in the fast-sweep limit.
Hence the plot of P (↑, 0) as a function of v/∆ 2 looks generally like an inverted version of P (↓, 0). As g/ω increases, the curve of P (↑, 0) shifts to the left, i.e. the drop of P (↑, 0) from one to zero occurs at smaller values of v/∆ 2 .
All P (↑, n) plots with n ≥ 1 are simple peaks at intermediate values of v/∆ 2 , with asymptotic values of zero in both the adiabatic and fast-sweep limits. The locations and heights of the peaks depend on the qubit-oscillator coupling in rather nontrivial ways. In the limit of small g/ω (Fig. 6a), the shapes of the P (↑, n) peaks are extremely weakly dependent on n, except for an overall reduction in the scale with increasing n; specifically The reason behind this simple relation is that all the factors in the formula for P (↑, n) in Eq. (8) except for the first and last factors, i.e. the product e −π∆ 2 0 /(2v) 1 − e −π∆ 2 n /(2v) , are approximately equal to one at the location of the peak, and the last factor can be approximated as 1 − e −π∆ 2 n /(2v) ≈ π∆ 2 n /(2v) ∝ (2g/ω) 2n /n! at the location of the peak. This result means, among other things, that one does not have much control over the number of created photons by controlling the sweep rate v. In fact, the probability of creating any photons in the oscillator remains small for all values of v/∆ 2 .
In the case where g/ω is not much smaller than one (Fig. 6b-d), the probabilities P (↑, n) have more structure. Comparing the plots of P (↑, n) as functions of v/∆ 2 for different values of n, the peak locations shift to higher values of v/∆ 2 and decrease in height as n increases. As g/ω increases, we obtain more separation between the different peaks, in addition to having higher peaks. As a result, for large g/ω (Fig. 6d), each P (↑, n) peak reaches a height of almost one. This result means that by choosing the appropriate value of v/∆ 2 , one can deterministically prepare a Fock state with a specific value of n. This result can be understood intuitively as follows: the gaps in the different avoided crossings (Eq. 9) scale as (2g/ω) n / √ n!. For g/ω ≫ √ n, the gaps obey the relation ∆ n ≫ ∆ n−1 . If we choose v such that ∆ 2 n ≫ v ≫ ∆ 2 n−1 , the first n avoided crossings (from |↑, 0⟩ to |↑, n − 1⟩) will be traversed fast, keeping the system in the state |↓, 0⟩, while the avoided crossing with the state |↑, n⟩ is traversed adiabatically, transferring the population to the state |↑, n⟩. It should be noted that to achieve the full separation between the probability peaks up to photon number n, it is necessary to have (g/ω) 2 /n > ∼ 10, which is based on the fact that the LZ probability goes from almost zero to almost one or vice versa in a v/∆ 2 n range of about one order of magnitude (see e.g. Fig. 6a). This requirement can make the gap ∆ n extremely small because of the factor e −2(g/ω) 2 and in turn require a proportionately slow sweep. As a result, it is extremely difficult to implement this scenario in a realistic experimental setup, as we will discuss in Sec. V.
We now examine the validity conditions for Eq. (8). We performed numerical simulations of the dynamics with a variety of parameters. We varied simulation parameters such as the size of the Hilbert space and the initial and final times to ensure that we obtain wellconverged results, i.e. results that are not significantly distorted by finite-size effects in the simulations. In particular, since we do not make approximations such as the weak-coupling and/or rotating-wave approximations, the results of the numerical simulations should be valid for small or large qubit-oscillator detuning and for weak or strong coupling. The such that we can no longer think of the dynamics as a sequence of independent probability redistribution processes. We therefore obtain large deviations between the simulation results and Eq. (8). A few differences are worth highlighting here. In the case of large g/ω, Eq. (8) predicts that an extremely slow sweep is required to remain in the ground state. The case with g/ω = 3 and ∆/ω = 100 shows that this prediction fails badly when the energy level structure becomes more complex than the simple mesh assumed in deriving Eq. (8). Another difference is that some of the probabilities have peaks with side shoulders, rather than the single peaks described by Eq. (8). We suspect that this dependence arises when more than two energy levels are significantly populated during the system's evolution, which can lead to quantum interference effects and nontrivial peak shapes. We finally point out that the progression of the curves as we increase ∆/ω is not unidirectional, as can be seen in the bottom-right panel, where the center of the peak moves to the right and then to the left as we go from ∆/ω = 1 to ∆/ω = 10 to ∆/ω = 100.
So far, we have considered the case of a qubit coupled to a single oscillator. The problem can be generalized to the case of a qubit coupled to multiple oscillators, or alternatively a multi-mode resonator [32]. The Hamiltonian is then given bŷ As with the case of a single oscillator, the energy level associated with the state |↓, 0, 0, 0, · · ·⟩ encounters an infinite sequence of avoided crossings. The order at which these avoided crossings are encountered depends on the energies of the oscillators. Specifically, the crossing with the state |↑, n 1 , n 2 , n 3 , · · ·⟩ occurs at the point vt = j n j ω j . Upon making the polaron transformation for the harmonic oscillators, we obtain the effective gaps of the avoided crossings: If the relevant avoided crossings are well separated from each other, the occupation probabilities can be approximately evaluated by calculating how much of the probability is transferred to the states |↑, n 1 , n 2 , n 3 , · · ·⟩ as the avoided crossings are traversed one by one. In particular, in the case when the largest avoided-crossing gaps correspond to excited states in multiple oscillators that are deep-strongly coupled to the qubit, with g j /ω j ≫ 1, at least some of the states |↑, n 1 , n 2 , n 3 , · · ·⟩ can be prepared almost deterministically by choosing a sweep rate that is adiabatic for the relevant avoided crossing but fast for all the avoided crossings that are encountered before it. This scenario obviously requires that all the gaps that are encountered earlier in the sweep must be much smaller than the one associated with the target state.
Another interesting question in this context is related to the case where one harmonic oscillator is coupled strongly to the qubit while all other oscillators are coupled weakly. In this case, we effectively obtain a quantum system that comprises the qubit and the strongly coupled oscillator, while all the weakly coupled oscillators serve as an environment for the quantum dynamics. The physics of LZ transitions in an open multilevel system is not very well understood [33]. However, we can make some general statements about this case. In particular, since the environment couples to the operatorσ z , one effect of the environment is to suppress quantum interference between multiple LZ transitions. As such, the environment pushes the final state probabilities to the expressions given in Eq. (8). for these purposes could be engineered dynamically by modulating the system parameters to obtain radically different effective system parameters [35][36][37][38]. In particular, by applying a parametric drive to the resonator frequency, the ratio g/ω can be enhanced to an effective value that is a few orders of magnitude larger than that of the undriven system with relatively moderate drive parameters [38]. The g/ω values used in all the panels in Figs. 1-7 can be achieved in a suitably rotating frame. The sweep rates required for observing excitation generation in a normal-superradiant or superradiant-normal quench are rather moderate and should be achievable in such a dynamical realization of the effective Hamiltonian. The deterministic generation of multi-photon Fock states (Fig. 6d) would require v values well below 1 MHz per second, even using the lab frame value of ∆ ∼ 10 GHz. As such, realizing this phenomenon in superconducting circuits using a simple sweep of ϵ or resonator frequency modulation as described in Ref. [38] seems unrealistic. It is possible that this phenomenon might be observable using alternative techniques or different physical systems. It is also worth noting that if we have a highly controllable and stable quantum system, it could be possible to use a nonlinear sweep in which the sweep rate is tuned in real time to be adiabatic or fast for different avoided crossings based on the intended purpose of the experiment.

VI. CONCLUSION
We have investigated the dynamics and in particular the final states of qubit-oscillator systems following a finite-duration quench in which one of the system parameters is swept between two values that correspond to qualitatively different states of the system. We obtained a number of results, some of which add theoretical insight to our understanding of quenched quantum system dynamics, and some of which are relevant to possible experimental realizations for practical applications. In our analysis, we paid special attention to the possibility of using parameter sweeps to generate states of interest for quantum information processing tasks.
In one case we investigated a sweep from the normal to the superradiant phase or vice versa. As expected, we found that excitations are created as a result of the quench between the two different phases. While the presence or absence of a sharp phase transition did not manifest itself in the final state probabilities, it led to qualitative differences in the dynamics of the quenched system. These results improve our understanding of quantum quenches through a phase transition point, especially as they relate to a rather unusual phase transition.
In another case we investigated an LZ sweep of a qubit coupled to an oscillator. In this case we found a variety of behaviours depending on the relation between the system parameters, including a case of extreme sensitivity in the final state to the sweep rate. From the point of view of quantum state preparation, a particularly interesting result that we found is the ability to almost deterministically prepare states with a specific number of excitation quanta in the oscillator. Although the peaks are narrow, requiring a fine tuning of the sweep rate, this protocol has the advantage of being robust against static fluctuations that shift the initial and final values of the qubit bias.
Our analysis also sheds light on the topic of excitation creation in a harmonic oscillator via an LZ sweep in a coupled qubit, which has not been investigated in the literature. We analyzed the predictions of a theoretical formula based on the assumption of independent LZ transitions and established the regime of validity of the theoretical formula.
In this work we focused on a single, linear sweep of parameters. Linear sweeps can also be repeated to create a situation of periodic driving [7,[39][40][41][42], which is one of the standard approaches to generating specific target states in quantum systems. Our results can form the basis for new protocols for state preparation in cavity-QED systems for quantum computing, communication and sensing applications.