The quantum transverse-field Ising chain in circuit QED: effects of disorder on the nonequilibrium dynamics

We study several dynamical properties of a recently proposed implementation of the quantum transverse-field Ising chain in the framework of circuit QED. Particular emphasis is placed on the effects of disorder on the nonequilibrium behavior of the system. We show that small amounts of fabrication-induced disorder in the system parameters do not jeopardize the observation of previously-predicted phenomena. Based on a numerical extraction of the mean free path of the system, we also provide a simple quantitative estimate for certain disorder effects on the nonequilibrium dynamics of the circuit QED quantum simulator. We discuss the transition from weak to strong disorder, characterized by the onset of Anderson localization of the system's wave functions, and the qualitatively different dynamics it leads to.

In [20], we have proposed and analyzed a circuit QED design that implements the quantum transverse-field Ising chain (TFIC) coupled to a microwave resonator for readout. The TFIC is an elementary example of an integrable quantum many-body system. Despite its simplicity, it still exhibits interesting features, e.g. a quantum phase transition (QPT), and therefore serves as a model example system in the theory of quantum criticality [21] and nonequilibrium thermodynamics [22]. Our circuit QED quantum simulator can be used to study quench dynamics, the propagation of localized excitations and other nonequilibrium phenomena in the TFIC, based on a design that could easily be extended to break the integrability of the system. While in [20] we have focused on an idealized implementation of the TFIC with perfectly uniform parameters, the main purpose of the present paper is to investigate the effects of disorder in the system parameters on the dynamical behavior of our quantum simulator.
The study of disorder effects on quantum simulators is relevant for two reasons. Firstly, on the more practical level, any real experimental system will come with a degree of unwanted disorder (especially in condensed matter settings). In the case of circuit QED systems, inhomogeneities of the system parameters are caused by fabrication issues as well as by static noise fields (e.g. produced by defects). It is important to verify that the basic behavior of a quantum simulator survives the amounts of disorder which are present in realistic systems or even to estimate the amount of disorder that can be tolerated. Secondly, on a more fundamental level, simulating quantum many-body systems with built-in (potentially tunable) disorder is interesting in its own right. Many physical phenomena, from free propagation of wave packets to quench dynamics to (quantum) phase transitions, can be affected in significant ways by disorder, and this leads to phenomena such as Anderson localization or disorder-induced phases.
To prepare for our study, we briefly review the system (section 2.1), discuss sources of disorder and how disorder scales with the tunable system parameters (section 2.2) and explain the mathematical approach to and some properties of the quantum Ising chain (section 2.3). We start our main discussion by considering the time-dependent correlations of the order parameter of the chain, where the finite-size effects and the long-time behavior will be analyzed in the absence of disorder (section 3.1). Based on this, we will move on to the spectrum of the resonator coupled to the quantum Ising chain in our system, which is closely related to the aforementioned time-dependent correlations. To that end, we employ a very useful approximation which we have introduced in [20] and which will presumably become important also for future studies of quantum many-body systems coupled to resonators. In this approximation, the full quantum many-body system is replaced by a bath of harmonic oscillators with an identical spectrum. We show here that this approximation actually works very well under appropriate circumstances (section 3.2). We then calculate the spectrum of the resonator coupled to a slightly disordered Ising chain and find that the effects of disorder on the spectrum are small (section 3.3). The Ising chain in our circuit-QED quantum simulator can be driven out of equilibrium in several ways. This allows one to perform various types of nonequilibrium experiments, a particularly appealing application of our setup. In our previous work, we have suggested to observe the propagation of a localized excitation through the chain or the nonequilibrium dynamics of the system after a quantum quench. Here, we show that the predicted phenomena are insensitive to a small amount of disorder in the system parameters (sections 4.1 and 4.2, respectively). Moreover, we provide a simple estimate of the amount of disorder that will qualitatively change the wave functions and, thus, strongly affect the dynamics even of small systems (that is, on the scale of neighboring artificial atoms). However, as argued above, it would be highly desirable to possess also a quantitative theory of disorder effects. Since the nonequilibrium dynamics of the uniform TFIC is determined by the ballistic propagation of quasiparticles (QPs; wave packets), we formulate and numerically verify for the weakly disordered case a relation between the mean free path of the latter and the parameters of the system and the disorder potential. By means of this relation we are able to predict the dynamical behavior of our quantum simulator given a certain disorder strength, and to estimate the amount of disorder that a particular experiment can tolerate (section 4.1).

Setup
We consider a circuit QED quantum simulator of the TFIC as proposed in [20]. It consists of a chain of N capacitively coupled charge-based superconducting artificial atoms [23], such as transmons or Cooper-pair boxes (the latter have to be biased to their charge degeneracy point [23] to properly simulate the TFIC). For a review on superconducting artificial atoms, see [23]. The first artificial atom is capacitively coupled to a microwave resonator (see figure 1). This resonator A is required for initialization and readout of the first artificial atom. For  [20]). Charge-based artificial atoms are capacitively coupled to their nearest neighbors. Coupling the first (N th) artificial atom to resonator A (B) allows one to use standard circuit QED techniques for initialization and read-out of the first (N th) artificial atom.
certain types of experiments, e.g. for measuring end-to-end correlators, one also needs a second resonator B, coupled to the N th artificial atom. For details of the implementation and the theoretical description of the system, see [20]. The system (at first, only with resonator A) can be approximately described by the Hamiltonian and H I is the Hamiltonian of the TFIC, Here, σ j x/z is a Pauli matrix. That is, the artificial atoms are considered as two-level systems (qubits), and their two states are described as spin states. The operators a † and a generate and annihilate a photon of energy ω 0 . The transition frequency j > 0 of the jth qubit corresponds to a local magnetic field acting on the jth spin in the usual interpretation of the TFIC. As such, it would be transverse to the direction of the qubit-qubit coupling J j . The latter can be either ferromagnetic (J j > 0, as in the geometry of figure 1) or anti-ferromagnetic (J j < 0, if the qubits in figure 1 are rotated by 90 • ). While in our previous work we have focused on the uniform case J j = J and j = for all j, here we are often interested in the case where these system parameters are explicitly nonuniform. This is because, on the one hand, a slight nonuniformity of the j and J j has to be expected from imperfections of the fabrication process. On the other hand, one can also intentionally detune one or several qubits (by threading the SQUID-like loops of the qubits with different fluxes) and observe how the system's properties change depending on the detuning.

Disorder and tunability of the system parameters
Let us discuss the flux tunability and the undesired disorder of the system parameters in some more detail. We will argue that the qubit transition frequencies j and the qubit-qubit couplings J j , when normalized to their respective mean values, may be assumed to be flux independent. This will be relevant for our theoretical description of the disorder in the system.
In reality, it should be possible to engineer the geometry of the qubits essentially uniform. That is, the areas of the qubits' SQUID loops, their charging energies and the coupling capacitances between the qubits will only vary weakly in the chain. However, the (flux-tunable) total Josephson energies E J ( ) of the artificial atoms should be experimentally harder to control since these depend exponentially on the properties of the Josephson junctions. For a flux-tunable (i.e. SQUID-type) artificial atom with two Josephson junctions [24], [24]. Thus, under the assumption of identical geometry, both for Cooper-pair boxes and for transmons the transition frequencies j ( ) of all qubits j scale with a j-independent function α( ) of the (global) flux , j ( ) = α( ) j (0). Here, α( ) = cos( π/ 0 ) for Cooper-pair boxes and α( ) = [cos( π/ 0 )] 1/2 for transmons. This result implies that the qubit transition frequencies, when normalized to their flux-dependent mean value, do not depend on and, hence, have the same statistical properties for all . Explicitly, the mean value of the j is given by j ( ) = α( ) j (0). Thus, the mean value of the j is flux tunable. However, the normalized qubit transition frequencies j ( )/ j ( ) are independent of , which must also be the case, for instance, for their standard deviation. This will become important for our numerical implementation of disorder in the j when we consider changes of the external magnetic flux .
The qubit-qubit couplings J j can also depend on the E J j ( ) and, thus, on the j . This is the case for transmons, where approximately J j ∝ ( j j+1 ) 1/2 ∝ (E J j E J j+1 ) 1/4 [20,25]. That is, the disorder in the j and the J j will not be independent for transmons. Moreover, j , J j and their mean values j and J j change with the external flux approximately in the same proportion (∝ [cos( π/ 0 )] 1/2 ). For Cooper-pair boxes, on the other hand, the J j depend only on charging energies and not on the E J j ( ) [20]. This means that the J j are not affected by changes of the external flux. Furthermore, the disorder in the J j should be less pronounced than and hardly correlated with the disorder in the j . Concerning the relative strength and the correlation of the disorder in the J j and the j , we remark that also static noise fields can play a role, producing some disorder also in the various charging energies of the system (in particular for Cooper-pair boxes, which have small electrostatic capacitances). Apart from that, disorder in the J j will turn out to have a much weaker effect than disorder in the j . These deliberations allow us to assume for simplicity that, both for Cooper-pair boxes and for transmons, disorder in the j and J j can be present to a comparable degree and that disorder in the j (J j ) would be uncorrelated with the disorder possibly present in the J j ( j ). We finally remark that many properties of the TFIC are determined by the ratio j /J j , since this ratio essentially (in the limit of weak disorder) determines the eigenstates of the system (see below). For standard transmons, the ratio j /J j is not straightforwardly flux tunable. One of the experiments we suggest to perform with our quantum simulator relies on the possibility to change the eigenfunctions of the system (cf section 4.2), which can be done only by changing the ratio j /J j . All other possible experiments discussed in this paper can be performed, in principle, with Cooper-pair boxes and transmons equally well, irrespective of the J j being fluxdependent or not [20]. Therefore, when plotting our results as a function of a flux-tunable system parameter, we will assume for definiteness that our circuit-QED quantum simulator of the TFIC is implemented with Cooper-pair boxes and that the J j do not change with the external magnetic flux.

The transverse-field Ising chain
The Hamiltonian (2) can be exactly diagonalized by means of a Jordan-Wigner transformation, which was first used in this context in [26,27]. This transformation maps the spin degrees of freedom to fermionic operators c j , c † j via σ + j = c † j exp(iπ j−1 k=1 c † k c k ) and yields Up to a constant − j j /2, this Hamiltonian is of the form Note that the conditions H = H † and {c j , c † j } = 1 require that A = A † and B = −B T . By introducing new fermions η k = N j=1 g k, j c j + h k, j c † j , such Hamiltonians can be transformed into the diagonal form H = k k (η † k η k − 1/2) + j A j, j /2 [26]. The components g k, j and h k, j of the vectors g k and h k and the excitation energies k of H are determined by defining normalized vectors φ k = g k + h k and ψ k = g k − h k and by solving the equations In our case and B is obtained by substituting A j, j = j → 0 and A j+1, j = −J j → J j into A. For uniform j and J j , the φ k , ψ k and k can be analytically calculated from equations (6) (see, e.g., [20]). For nonuniform system parameters, these quantities have to be determined numerically. In both cases, the Hamiltonian H I of the TFIC can be written in the form and knowledge of the φ k and ψ k allows one to express spin observables in terms of the η k -fermions, which is the basis of many of our calculations. For instance, We collect some important facts about the TFIC. In the uniform case, Here, J = |J | and ξ = /2J is the normalized transverse field. The possible values of k are solutions of sin k N = ξ sin k(N + 1). For N → ∞, the uniform TFIC undergoes a second-order QPT at = ±1 from a ferromagnetic [ ∈ (0, 1)] or an anti-ferromagnetic [ ∈ (−1, 0)] ordered phase with doubly degenerate eigenstates (one k → 0) to a paramagnetic disordered phase with k > 0 for all k. The QPT is signaled by the disappearance of long-range correlations in σ x . This QPT will also occur in a nonuniform system (at some mean transverse field strength j ) [21]. However, there can be weakly (dis)ordered Griffith-McCoy 'phases' in the vicinity of the critical point [28][29][30][31]. Finally, we introduce a convenient notation for nonuniform j and J j . In this case, we will frequently write j = τ j and J j = J τ j , where τ j and τ j usually have mean 1, or, if j and J j follow probability distributions, expectation value 1. We will refer to as the 'mean' qubit transition frequency, even if = j is the expectation value of a probability distribution and the actual mean value j is (for finite N ) in general different from . We use the same convention for the qubit-qubit coupling J . Furthermore, we define the local and the 'mean' normalized transverse magnetic field, ξ j = j /2J j and ξ = /2J . Note that in general both ξ = ξ j and ξ = ξ j (but for the probability distributions we will consider, (ξ − ξ j )/ ξ j < 1%). We usually characterize H I by the parameters ξ , J , τ j and τ j . Under the assumptions formulated in section 2.1, and thus ξ are fluxtunable without changing the τ j in the proposed circuit-QED quantum simulator of the TFIC.

Spectrum of the system
In order to provide a guideline for the initial experimental characterization of our setup, we have calculated in [20] the transmission spectrum S of the resonator as a function of the probe frequency ω and the flux-tunable qubit transition frequency (see below equation (28)). To that end, we have first calculated the spectrum of the bare TFIC for coupling to the first qubit via σ 1 x ,ρ which is the Fourier transform of the qubit autocorrelator ρ(t) = σ 1 . We have argued that for sufficiently large (but finite) N , qubit decay processes will render the measured spectrum continuous and akin to the spectrum one would obtain by taking the limit N → ∞ in the calculation of ρ. Assuming small coupling g/ω 0 1 of the first qubit and the resonator, we have then considered the TFIC as a linear bath for the resonator, and this approximation allowed us to calculate the resonator spectrum S in the coupled system. In this section, we add some remarks on the interpretation of the autocorrelator, the transition N → ∞ and the linear approximation. Moreover, we discuss how a small amount of disorder in the qubit parameters due to imperfections in the fabrication process affects the resonator spectrum S.

Time-dependent correlations in the transverse-field Ising chain
By means of the spin-free-fermion mapping described in section 2.3, one readily finds that Here and in the following, expectation values are calculated under the assumption of zero temperature. This is justified because the band gap of the Ising chain is of the same order of magnitude as the qubit transition frequencies ∼ 5 GHz (except near the critical point) and, thus, much bigger than the usual mK temperatures of a cryogenic environment. In the uniform case j = and J j = J , where explicit expressions for φ k and k can be found, the limit N → ∞ can be taken analytically and yields [20] Here, (x) is the Heaviside step function and (k) stands for k with continuous k (equation (10)). The first term on the rhs of (13) causes a nonzero mean value of Re ρ(t) in the ordered phase. Figure 2 shows Im ρ(t) for ξ = 8 in the cases N = 20 (equation (12)) and N → ∞ (equation (13)) (the time evolutions of Re ρ and Im ρ are qualitatively similar and agree for |ξ | 1 up to a phase). For small times, the curves coincide (the second covers the first). However, the finite size of the TFIC with N = 20 causes a revival of ρ at . This can be understood in the following way. The autocorrelator ρ is related to the linear response σ 1 x (t) of the TFIC to a perturbation ∝ δ(t)σ 1 x relative to the equilibrium value σ 1 x = 0. Indeed, Kubo's formula predicts that σ 1 x (t) ∝ Im ρ(t). The δ-pulse at t = 0 forces the first spin in the −x-direction. This local excitation in position space is composed of many excitations in k-space. Since most of them have velocity v [20], the local excitation propagates with velocity v through the system, is reflected at the far end of the chain and causes revivals of ρ at multiples of T r = 2N /v. To further clarify the transition N → ∞, we note that for large t, ρ has a standard deviation from its mean ∝ 1/ √ N . This can be expected from (12) since |ρ(t)| 2 ∼ 1/N 2 k,k e it ( k − k ) and, for t → ∞, all terms in the sum except for those with k = k will cancel. In general, the t → ∞ fluctuations that we find for all time-dependent observables considered in this work are due to the finite system size and decrease with N (but not all of them behave like ∝1/ √ N ).

Spectrum of the resonator-the linear approximation
Taking the Fourier transform of equations (12) and (13) yields the spectrumρ(ω) of the TFIC for a force that couples to σ 1 x for finite N and N → ∞, respectively. In order to calculate the spectrum S of the resonator, whose coordinate (a † + a) couples to σ 1 x (cf equation (1)), we have suggested [20] a useful approximation: we consider the TFIC as a linear bath for the resonator. That is, we replace the TFIC by a set of harmonic oscillators having the spectrumρ of the TFIC. This approximation can be straightforwardly generalized to other contexts, where a different many-body system couples to a resonator. It is justified in the limit of small qubit-resonator coupling g/ω 0 1, as we discuss in the following.
The linear approximation for the TFIC-bath fails as soon as probing the resonator sufficiently excites the TFIC so that its nonlinearity becomes important. Thus, the linear approximation requires a small coupling g and is worst if the TFIC is on resonance with the resonator (ω 0 within the band k of the TFIC). The 'most nonlinear' bath possible for the resonator, that is, the bath whose nonlinearity becomes important for the smallest value of g, is a bath consisting of only a single qubit on resonance with the resonator. If the linear approximation is adequate for such a system in the limit g/ω 0 1, it will be also sufficient for our purposes. Therefore, we now consider the case N = 1 and = ω 0 of equation (1) and calculate the spectrum of the resonator by linearizing the single-qubit bath. Since the atomic Hilbert space is small for N = 1, we can then numerically check the accuracy of our approximation. We also compare our approximation with the resonator spectrum calculated analytically within the rotating wave approximation (RWA), which is the standard approximation of H in this specific situation.
For N = 1 and = ω 0 , the Hamiltonian H (equation (1)) becomes That is, the resonator coordinate (a † + a) couples to a single-qubit bath with the Hamiltonian H q = ω 0 2 σ z via a force gσ x . The spectrum of this force is where the time evolution of σ x and the expectation value q . q are to be calculated with respect to (the ground state of) H q . Now we linearize the system and replace with bosonic b, b † and parameters g and w to be determined. In (16), the resonator couples to a force g (b † + b) exerted by a bath that consists of a single harmonic oscillator with the Hamiltonian H ho = wb † b. The spectrum of this force reads Thus, we choose g = g and w = ω 0 such thatF ho =F q . With this substitution, we now calculate the autocorrelator of the resonator coordinate and its Fourier transform, the resonator spectrum according to equation (16). To that end, we express the resonator coordinate (a † + a) in terms of Using (20), one readily finds that Before we go on and compare these approximate analytical results with numerical finitesize calculations for H N =1 (equation (14)), we calculate the same quantities on the basis of the standard approximation to H N =1 for g/ 0 1, the RWA (see, e.g., [32]). This will be a helpful benchmark for estimating the quality of the linear approximation. In the RWA, the Hamiltonian H N =1 reduces to the Jaynes-Cummings Hamiltonian This Hamiltonian can be straightforwardly diagonalized, and one can therefore analytically calculate the autocorrelator ρ RWA (t) and the spectrumρ RWA (ω) of the resonator in the approximation provided by H RWA , On the basis of (23), the results (24) and (25) are exact. The autocorrelator and the spectrum of the resonator can also be calculated numerically after truncating the photonic Hilbert space. This is achieved by expanding H N =1 and the resonator coordinate (a † + a) in the product basis {|s z , ν }, where s z =↑, ↓ and ν ∈ N 0 , and dropping all matrix elements with ν > ν max . In this finite-size approximation, the eigenvalues E n and eigenvectors |n of H N =1 can be numerically calculated (n = 0, . . . , n max = 2ν max + 1) and give ρ(t) andρ(ω) according to  27)). In both plots, we choose g/ω 0 = 0.12, which is the largest ratio of g/ω 0 used in this work and in [20]. The autocorrelator ρ(t) of the resonator (red) is well approximated both by the RWA (ρ RWA , blue) and the linear approximation (ρ lin , green), and the quality of these approximations is essentially equal. For small t, the linear approximation might be even more accurate than the RWA, but becomes worse at large t. This can be understood in the frequency domain. In figure 3(b), we plot the spectral weights of the delta peaks in the spectraρ,ρ RWA andρ lin (red, blue and green) at the corresponding peak positions. The spectrumρ contains also delta peaks at higher frequencies than the ones plotted, but their weight is virtually zero (2π | 0|(a † + a)|n | 2 < 10 −7 for all n = 1, 2). Both approximations yield good predictions for the positions and the spectral weights of the peaks inρ. The RWA is more precise in predicting the peak positions and the linear approximation in predicting the spectral weights (note, however, that the peak positions inρ lin andρ RWA agree up to first order in g/ω 0 ). Thus, the linear approximation is more precise for small t, in particular at t ≈ 0 and where the envelope of ρ(t) has a minimum, but becomes worse for large t. In summary, we conclude that even for the situation N = 1 and = ω 0 , the linear approximation yields good results for the autocorrelator and the spectrum of the resonator in the limit g/ω 0 1 that are qualitatively comparable to the usual RWA in this context. This implies that the linear approximation is well justified in our calculation of the spectrum of a resonator coupled to a TFIC.

Spectrum of the resonator-disorder effects
The linear approximation for the TFIC allows one to express the spectrum S(ω) of the (coupled) resonator as a function of the spectrumρ(ω) of the TFIC [20], Here, κ is the full-linewidth at half-maximum of the Lorentzian spectrum of the uncoupled (g = 0) resonator and χ(ω 2 ) denotes the principal-value integral χ( . This result is actually general and holds for any linear bath coupled to a resonator, with an arbitrary spectrumρ. Plots of S, withρ(ω) being the Fourier transform of (13), are presented in [20]. However, in an actual implementation of the proposed setup, the qubit parameters J j and j will not be perfectly uniform, due to imperfections in the fabrication process. We now investigate how this modifies the characteristic features of the spectrum S of the uniform system. It is known in the field of random-matrix theory that disorder would have to be very strong in order to have a dominant effect on (average) spectra. We will observe the same here, in this concrete model system.
For a nonuniform TFIC, no closed analytical expressions forρ(ω) are available. Thus, we have to consider finite system sizes and calculate numerically the relevant quantities, specifically, the spectrum of a finite-size nonuniform TFIC, which is the Fourier transform of equation (12). To take the effect of qubit decay processes into account, we phenomenologically broaden the delta peaks in (29) and replace them by Lorentzians of width γ around the k . We model the nonuniformity of the qubit parameters by writing j = τ j and J j = J τ j and choosing τ j and τ j to be random variables, which follow Gaussian distributions with means 1 and standard deviations σ τ = σ τ = 0.02. Uniformity of the qubit parameters j and J j of this degree will turn out to be sufficient for all proposed experiments. Much stronger disorder is not generally tolerable, as we will see below. However, from the experimental data for a sample with three (even spatially separated) qubits presented in [33], we calculate a standard deviation of the qubit transition frequencies from their mean of 0.8% (for zero flux bias). Thus, the requirements on the uniformity of j and J j appear to be attainable. With a typical set of system parameters that was also used in [20], we numerically calculateρ(ω) according to (29) and the corresponding resonator spectrum S according to (28). In order to judge the effects of disorder, we also reproduce our calculation of S for the corresponding uniform system [20] (figure S6). Figure 4 shows S as a function of ω and the (mean) normalized transverse field ξ = /2J for the uniform system (figures 4(a) and (b)) and for a typical disorder configuration (figures 4(c) and (d)). In the uniform case, the signatures of the QPT at ξ = 1, the dispersive shift of the resonator frequency and, on resonance, the double peak with a separation of 4J (rather than 2g as in the case N = 1) that we have discussed in detail for N → ∞ in [20] are clearly visible also for N = 20. These characteristic features are insensitive with respect to a small amount of disorder in the system parameters, as figures 4(c) and (d) demonstrate. We remark that in several recent circuit QED experiments the qubits have been found to be unexpectedly hot [34][35][36]. A corresponding non-negligible equilibrium population of the excited many-body eigenstates of the Ising chain in our setup would lead to additional lines in the described spectroscopy experiment, at frequencies smaller than the bandwidth of the Ising chain. In the experimentally realistic case that the Ising chain is deeply in the paramagnetic phase ( 2J ), these resonances at ω 4J (the bandwidth of the chain in the paramagnetic phase) would be well below the lower band edge − 2J . Thus, they would be distinguishable from the band of the Ising chain as plotted in figure 4, and their intensity might allow one to estimate the spurious population of the excited states. However, for the proposed time-domain experiments with our circuit QED quantum simulator that we discuss in the following sections, a non-negligible equilibrium excitation of the Ising chain might necessitate post-selection or initialization techniques.

Disorder effects on the system dynamics
A particularly interesting application of the proposed system would be to simulate the nonequilibrium dynamics of the TFIC. In [20], we have suggested to experimentally track the propagation of a localized excitation in the (uniform) TFIC that can be easily created in our system and to measure the system dynamics after quenching the transition frequencies of all qubits. In this section, we show that none of the predicted experimental results changes qualitatively if the parameters of the TFIC are slightly disordered, as has to be expected in reality. Stronger disorder, accessible, e.g., by deliberately detuning individual qubits, is shown to produce qualitatively different physics in the previously proposed experiments, like Anderson localization of the propagating excitation. For the realistic case J , we give an estimate of the corresponding disorder strength. Finally, we develop a quantitative theory of the effects of weak disorder on the system's nonequilibrium dynamics that explains the results of numerical experiments with the disordered TFIC. This theory might be helpful for experimentalists to estimate system and disorder parameters for successfully performing nonequilibrium experiments with the TFIC (e.g. for a given measurement resolution) without having to do numerical simulations.

Propagation of localized excitations
For the first type of experiments we have suggested in [20], it is assumed that the TFIC is deeply in the paramagnetic phase (ξ 1) and detuned from the resonator. In this situation, the TFIC is essentially decoupled from the resonator and its ground state is characterized by σ j z ≈ −1. Applying a fast π-pulse to the first qubit thus creates a localized excitation in the system that subsequently propagates through the chain due to the qubit-qubit coupling J . The time evolution of the observable σ j z after the π-pulse can be approximately described by [20] σ j We plot this result in figure 5(a) for all j in a chain of length N = 20 and for a mean normalized transverse field ξ = /2J = 8 (the same system parameters as in [20]) and, again, we randomly choose j = τ j and J j = J τ j according to Gaussian distributions with standard deviations of 2% from the mean values and J as before (right panel). The experimentally measurable observable σ 1 z (t) is singled out in the left panel. The propagation of a localized excitation through the chain, and its reflection at the far end of the chain that leads to a distinct revival of σ 1 z (t) at t ≈ N /J , are still clearly visible in this slightly nonuniform system. If the transition frequencies j of the qubits can be tuned individually, the effective length of the TFIC has been shown to be adjustable by strongly detuning one qubit from the others [20]. This holds true also for a slightly nonuniform system: figure 5(b) shows the typical result for a system with the same parameters and disorder strength as in (a), but with qubit 11 strongly detuned by setting τ 11 = 1.3. This result is qualitatively identical with the result for the corresponding nondisordered system [20]. The strong nonuniformity at j = 11 acts as a barrier for the propagating excitation and leads to its reflection. Thus, it effectively changes the length of the TFIC. Having shown that the experiments with propagating localized excitations proposed in [20] yield qualitatively the same results for ordered and slightly disordered systems, we now proceed and study disorder effects on this type of experiment quantitatively. Parts of the following analysis also apply to other nonequilibrium experiments with the TFIC, as will be discussed in the context of quantum quenches (section 4.2).
Since it is assumed that the system is deeply in the paramagnetic phase, the mean qubit transition frequency is larger than the modulus of the mean qubit-qubit coupling J , /J 1. As before, we further assume uncorrelated disorder of the system parameters via j = τ j and J j = J τ j , where τ j and τ j follow Gaussian distributions with standard deviations σ τ and σ τ from 1. That is, for σ τ = σ τ , the absolute variation of the j will be larger than the absolute variation of the J j . Therefore, the dynamics of the system may be expected to be much more sensitive to increasing σ τ than σ τ . Moreover, one may expect that disorder effects start to qualitatively affect the system dynamics even of small systems (that is, on the scale of neighboring qubits j and j + 1) when the disorder in the qubit transition frequencies becomes comparable to the modulus of the mean qubit-qubit coupling, σ τ = J . These deliberations are confirmed by numerical experiments: we first consider the wave functions g k, j and h k, j in position space (η k = N j=1 g k, j c j + h k, j c † j ). For zero disorder, they are extended over the whole chain (except for the mode with k → 0 in the ordered phase [37]). Increasing σ τ localizes the wave functions much more strongly than increasing σ τ , and the localization length of the wave functions indeed reduces from many ( 1) sites to a few ( 1) sites at σ τ ≈ J . Correspondingly, the propagation of an excitation initially localized at site 1 is only weakly affected by disorder in J . However, if σ τ J/ , it propagates only a few sites before becoming completely trapped due to the disorder. This manifestation of Anderson localization [38] is illustrated in figure 6(a), where we have used the same system parameters as in figure 5, but we have randomly chosen τ j and τ j according to Gaussian distributions around 1 with standard deviations σ τ = J/ = 0.0625 = σ τ . For definiteness, we always choose σ τ = σ τ in the following. We have seen that for |ξ | 1 (paramagnetic phase) the effective disorder strength in the quantum Ising chain is set by σ τ /J ∝ σ τ | |. Now we try to determine how the relevant observables in the currently considered type of experiment depend on this quantity. The observable we focus on in the following is the maximum excitation probability (maximized over time) of the jth qubit caused by the propagation of the localized excitation through the disordered chain. In an experiment, one would for instance create an excitation of the first qubit and measure the excitation probability of some other (e.g. the N th) qubit as a function of time. The maximum excitation probability of the jth qubit is an important quantity since it will determine whether the effect of the propagating excitation can be measured at site j, given a certain measurement resolution. In a single disordered system, the maximum excitation probability of qubit j will depend on the specific (random) disorder configuration of this system. Therefore, a study of the effect of disorder as characterized by the statistical quantity σ τ |ξ | can only refer to the statistical average of the maximum excitation probability of qubit j in one disordered system over an ensemble of many disordered systems (disorder configurations), all chosen according to the same probability distribution. Stated as a formula, this ensemble average of the maximum excitation probability of qubit j is given by Here, the double overbar · denotes the ensemble average over many disordered systems (disorder configurations) with the same system and disorder parameters ξ , J and σ τ = σ τ . This average is taken after one has maximized σ j z (t) for a specific disordered system over time. Our goal is to find the explicit functional dependence of p j σ τ ,|ξ | on σ τ and |ξ | (in fact, we expect dependence only on the product σ τ |ξ |). Note that we assume that p j σ τ ,|ξ | depends neither on the sign of ξ nor explicitly on the mean qubit-qubit coupling J , but only on the ratio of and J (via |ξ |). This is strictly true for σ τ = σ τ = 0. By explicitly solving equations (6) for this case [20], one can show that after substituting ξ → −ξ , the new allowed wave vectors are q = π − k with q = k , φ q, j = (−1) N − j φ k, j , and ψ q, j = (−1) N − j ψ k, j . With that one can easily see that equation (30) does not depend on the sign of ξ . Moreover, φ k and ψ k are independent of J (which also follows from equations (6)), and k ∝ J such that changing J corresponds only to a rescaling of time. The influence of disorder, however, is essentially set by σ τ |ξ | (for |ξ | 1), as we have argued above. Consequently, we may take p j to be independent of J and of the sign of ξ . Nevertheless, to keep notation short, we write ξ instead of |ξ | for the remainder of this section. For simplicity, we first focus on a semi-infinite system (N → ∞) and discuss later the increase of p j at the end of the chain (due to the refocusing of the dispersed wave packet of the propagating excitation).
As usual for disordered systems (e.g. [39]), we will try to characterize the disorder effects on the ensemble-averaged maximum qubit excitation p j σ τ ,ξ via a mean free path. To that end, it pays to first discuss in more detail the uniform case, p j 0,ξ . Even there, analyzing the propagation of the dispersive wave packet that determines the maximum excitation probability of a qubit requires some care. For ξ 1, this excitation probability does not depend on ξ . This is because the dispersion relation of the TFIC becomes that of the tight binding model, cos k). Thus, only sets the band gap but does not influence the shape of the dispersion relation. Except for the aforementioned boundary effects, p j 0,ξ also does not depend on N . This is evident from figure 6(c), where we plot p j 0,ξ for several ξ and N . The curves for different ξ but the same N lie almost on top of each other (henceforth, we drop the index ξ from p j 0,ξ ), and curves for different N can be distinguished only by the boundary effects, that is, by the strong increase of p j 0 at j = N (which will be discussed later). The decay of the p j 0 with j is relatively slow (slower than 1/j), which should considerably simplify the experiments proposed in [20]. This slow decay of p j 0 can be understood from the dispersion relation k of the system which, for ξ 1, is quadratic in k at k ≈ 0, π, and linear at k ≈ π/2: if an initially localized wave packet with width s and momentum q, ψ(x, 0) = α e −x 2 /2 s 2 +iq x , α = (s 2 π) −1/4 , is evolved in time by the Hamiltonians H 1 = h 1 k and H 2 = h 2 k 2 , respectively, one finds that That is, for H 1 , the maximum and width of the probability distribution for finding the particle at a position x are constant, while for H 2 and strong initial localization (or large times) the width is ∝ t and the maximum is ∝ 1/t. As the dispersion relation of the TFIC interpolates between these two cases, one may expect a decay of p j 0 slower than 1/j. Coming back now to the disordered case, one might suspect that the ensemble-averaged maximum qubit excitation p j σ τ ,ξ is related to the corresponding quantity for a nondisordered system p j 0 via an exponential decay, governed by a finite mean free path l σ τ ,ξ for the propagation of the localized excitation, If (34) holds, should be independent of j. This observation can be used to check our assumption (34). We numerically calculate p j σ τ ,ξ for all combinations of ξ = 3, . . . , 8 and 100 × σ τ ∈ {1, 1.5, 2, 2.5, 3, 4, 5, 6, 7, 8} in a chain of length N = 20, and we average over 100 disorder configurations. This turns out to be a good compromise between calculation time and ensemble and system size as long as the effective disorder σ τ ξ is not too small (see below). With these p j σ τ ,ξ , we calculated the rhs of (35) for j = 5, . . . , 16. Other j are not considered, in order to minimize boundary effects. Our results for j = 5, . . . , 16 are approximately equal, with the ratio of standard deviation to mean value being < 0.1 for given σ τ and ξ . We note that for very weak effective disorder σ τ ξ 0.1 we have to average over 500 disorder configurations such that this ratio is < 0.1, because with decreasing ratio p j 0 / p j σ τ ,ξ the slope of the logarithm on the rhs of (35) increases. Thus, the numerical data seem to confirm our assumption (34), and the influence of disorder on the considered experiment is captured by a mean free path l σ τ ,ξ . In our subsequent analysis, we try to find simple expressions for this quantity.
The propagation of the localized excitation in the Gaussian disordered TFIC is akin to the propagation of a particle in an uncorrelated random potential V (r ) with V (r )V (r ) = V 2 0 δ(r − r ). To lowest order in perturbation theory (Fermi's golden rule, e.g. [39]), the mean free path of the latter decreases as the inverse square of the disorder strength, ∝ 1/V 2 0 . In our case, the effective disorder strength is determined by the dimensionless quantity σ τ ξ . Therefore, we expect that To check this, we calculate l σ τ ,ξ for the same combinations of σ τ and ξ as before by averaging the rhs of (35) over j = 5, . . . , 16. Then we fit the function l(σ τ , ξ ) = 1/σ a τ ξ b to our data for l σ τ ,ξ . We find the exponents a ≈ 2.002 and b ≈ 2.071, which comes close to our expectation of a = b = 2. Numerical data and fit are plotted on the log-log scale in figure 6(b). As long as the effective disorder strength is not too big (σ τ ξ 0.2), l(σ τ , ξ ) with the fit values of a and b reproduces the numerically (by ensemble-averaging) extracted mean free path l σ τ ,ξ . Here, one may attribute the deviations of a and b from 2 to the finite ensemble sizes. For stronger disorder, however, the fit of l(σ τ , ξ ) = 1/σ a τ ξ b begins to deviate from l σ τ ,ξ . Thus, higher-order effects (beyond Fermi's golden rule) and/or the disorder in J seem to be no longer negligible.
Finally, in setups with a second readout resonator (cf figure 1), the maximum population p N σ τ ,ξ of the N th qubit will be an experimentally relevant quantity. Since the dispersed wave packet of the propagating excitation is refocused at the end of the chain, the maximum excitation probability of the N th qubit is considerably enhanced compared to nearby bulk qubits (see figure 6(c)). It turns out that p N σ τ ,ξ can also be estimated by means of (34), the mean free path (36), and the value of p N 0 for the corresponding nondisordered system, which we plot for N = 1, . . . , 50 in figure 6(d).
Summing up, equations (34) and (36), together with figures 6(c) and (d), allow one to easily estimate suitable system and disorder parameters for successfully implementing the presently considered type of experiment. For instance, if in a system with N = 30 and ξ = 4 the N th qubit should get a population of p 30 σ τ , = 0.3 (which corresponds to max σ 30 z (t) = −0.4), then one can roughly (i.e. averaged over many systems) afford a standard deviation of the qubit transition frequencies from their mean of σ τ = (N ξ 2 ) −1 log p 30 0 /0.3

1/2
≈ 0.034, where we have extracted p 30 0 ≈ 0.53 from figure 6(d). A typical result for these parameters is plotted in figure 7(a). Here, the maximum excitation probability is found to be p 30 0.034,4 ≈ 0.29 (since max t σ 30 z (t) ≈ −0.43). We remark that the foregoing deliberations only hold for uncorrelated disorder of the system parameters and do not take into account qubit decay. Correlated disorder can yield qualitatively different results and has to be studied explicitly via equation (30). We also remark that if the j are individually tunable, it becomes possible to study the propagation of localized excitations in arbitrary potentials. For instance, it might be interesting to choose j = [1 + √ 2σ τ cos(2π j/N )] and to compare the system dynamics with the Gaussian disordered case.
For large N , both distributions of j have the same mean and the same standard deviation, but in the former case the system is not disordered and the localization of the propagating excitation is much weaker than in the genuinely disordered case. Figure 7(b) shows the propagating excitation in such a system with N , ξ and σ τ as in figure 7(a) (with uniform J j ).

Quench dynamics
The second type of nonequilibrium experiment we have proposed in [20] relies on the possibility to rapidly change the transition frequency of a superconducting qubit in a circuit QED system by tuning the magnetic flux through its SQUID loop. This has been shown to be possible virtually instantaneously on the dynamical time scale of a circuit QED system [4,6,7], without changing the system's wave function. Let us now assume that the circuit QED quantum simulator of the (uniform) TFIC proposed in [20] is implemented with Cooper-pair boxes. For this system, such a sudden change of all j = corresponds to a global quantum quench of the normalized transverse magnetic field ξ = /2J . We remark that one can also produce quenches of ξ by using transmons in a non-standard parameter regime instead of Cooper-pair boxes, or by using usual transmons with tunable coupling capacitances [20,40]. We also remark that the observation of the phenomena described in the following will set higher requirements on the energy relaxation and phase coherence times of the collective many-body quantum states of the Ising chain than the experiments proposed in sections 3.3 and 4.1. The global quantum quench brings the Ising chain in a globally excited state whose time evolution has to be coherent on the time scale N /J of these phenomena (see below). Nevertheless, meeting this constraint seems feasible, since even for N = 30 and a moderate coupling strength J/2π = 100 MHz, we find that N /J ∼ 50 ns, which is far below the energy relaxation times T 1 ∼ 7.3 µs and coherence times T 2 ∼ 500 ns achieved for individual Cooper-pair boxes [41].
The nonequilibrium dynamics of the TFIC following a quantum quench is currently subject to much theoretical research, e.g. [22,[42][43][44][45][46][47][48][49][50][51][52][53], and should be experimentally observable with our circuit QED quantum simulator. In this context it is usually assumed that for t < 0 the system is in the ground state |0 a of a Hamiltonian H I,a (characterized by ξ a ). At t = 0, the overall transverse field is changed, a → b , and the nonequilibrium time evolution of some observable O under H I,b is investigated, Here, and G and H are matrices that, respectively, contain the g k and h k as columns. In these equations, a quantity carrying the index a or b is to be calculated from equations (6) with parameter set a or b.
To implement disorder of the system parameters before the quantum quench, we write again a j = a τ j and J j = J τ j , and we randomly choose τ j and τ j according to Gaussian distributions with standard deviations σ τ and σ τ from 1. As we have argued in section 2.1, tuning the flux through the SQUID loops of the qubits only changes the mean qubit transition frequency a → b (and, thus, the mean transverse field ξ a = a /2J → ξ b = b /2J ), but leaves τ j , J and τ j unaffected. Hence, by fixing ξ a/b and σ τ/τ , the system is fully specified before and after the quench (as in section 4.1, the absolute values of a/b and J can be absorbed in the time scale J t of the dynamics), and we are ready to evaluate equations (38) and (39). Figure 8(a) shows the local magnetization σ j z (t) for all j and for the same system parameters as in figure 4 of [20], but with a/b j and J j having standard deviations σ τ = σ τ = 2% around their mean values (right panel). The experimentally easily measurable trace of σ 1 z is singled out in the left panel (black). For comparison, we also plot (green) the local magnetization of the first qubit σ 1 z of the uniform system (as plotted in the left panel of figure 4 of [20]). Correspondingly, figure 8(b) shows (39) for a uniform system as in figure S6 of [20] (green), and with 2% disorder in a/b j and J j (black). The plots demonstrate that the quench dynamics of the considered observables is not qualitatively affected by the presence of a small amount of disorder.
For a more systematic analysis of the disorder effects on the quench experiments considered here, we make use of our findings for the mean free path of a propagating localized excitation from the previous section. This is possible because the quench dynamics of the TFIC is governed by the propagation of QPs through the system [20,43,44,46,48,54]. These correspond to flipped spins, essentially like the localized excitation of the previous section. Indeed, if the system is initially in the paramagnetic phase, the time evolution immediately after the quantum quench e −itH b |0 a ∝ j e −itJ (ξ b −ξ a )/ξ a σ j x σ j+1 x |0 a flips pairs of adjacent spins so that they point in the +z-direction. Due to the qubit-qubit coupling, these local excitations propagate as QPs with velocity v ≈ 2J through the chain. For an interpretation of the quench dynamics and the time scales indicated in the plots (all of which scale as N /J ) in terms of these QPs, see [20]. If ξ b is in the paramagnetic phase, the mean free path l of the QPs in a disordered TFIC can be estimated by l = 1/(σ τ ξ b ) 2 according to the previous section. The characteristic quasi-Tperiodic behavior (T = N /v) of the local magnetization after the quench in the nondisordered TFIC can be understood as a revival of coherence each time QPs initially generated at the same spot meet again [20,48]. This happens when the QPs have traveled multiples of the chain length N . If there should be a significant probability that two contiguously generated QPs meet again at least once before being scattered and thus decrease the local magnetization at t = T , the mean free path has to be sufficiently large, l > 2N . The appearance of significant end-toend correlations after the quench (that are stronger than those for t → ∞) requires that QPs in a disordered TFIC of length N = 30 after a quench of the mean normalized transverse field ξ = 8 → 1.5 (black), along with the corresponding trace for a uniform system (green). In both plots the qubit transition frequencies j and qubit-qubit couplings J j are randomly chosen according to Gaussian distributions with standard deviations of 2% from the mean values and J . generated in the middle of the chain reach the edges of the chain without being scattered, hence l > N . We have performed numerical experiments which indeed suggest that the corresponding values of σ τ mark the transition to a degree of disorder where the described phenomena are no longer present. In that sense, the distinctive features of the quench dynamics of the end-toend correlator are less sensitive to disorder than those of the local magnetization (and, due to the shorter time scale, less sensitive to decoherence or decay). We finally note that also here the effective chain length can be adjusted by strongly detuning individual qubits (this can also be used to create local quantum quenches by 'joining' two initially independent chains) and arbitrary effective potentials j can be chosen.

Conclusion
In the quest for controllable large-scale quantum systems, the framework of circuit QED offers several advantages, such as fast, high-fidelity readout, a great flexibility in design and steadily increasing coherence times. However, a potentially significant disadvantage arises from the hardly avoidable static noise and disorder sources in these man-made devices. The central result of this work is that also in this respect, there is reason to be optimistic: the requirements on the homogeneity of the system parameters for observing interesting (and predictable) many-body physics in a circuit QED system are not too high to be achievable with present-day or nearfuture technology. This underlines the prospects of circuit QED as a promising platform for implementing quantum simulations of complex quantum many-body Hamiltonians. In addition, we have shown that circuit QED quantum simulators could be used to study deliberately the effects of tunable disorder on quantum many-body dynamics.