{'fraction_non_alphanumeric': 0.036729751547223966, 'fraction_numerical': 0.00735491972374204, 'mean_word_length': 4.601934430347947, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact calculation of phonon effects on spin squeezing

Theoretical models of spins coupled to bosons provide a simple setting for studying a broad range of important phenomena in many-body physics, from virtually mediated interactions to decoherence and thermalization. In many atomic, molecular, and optical systems, such models also underlie the most successful attempts to engineer strong, long-ranged interactions for the purpose of entanglement generation. Especially when the coupling between the spins and bosons is strong---such that it cannot be treated perturbatively---the properties of such models are extremely challenging to calculate theoretically. Here, exact analytical expressions for nonequilibrium spin-spin correlation functions are derived for a specific model of spins coupled to bosons. The spatial structure of the coupling between spins and bosons is completely arbitrary, and thus the solution can be applied to systems in any number of dimensions. The explicit and nonperturbative inclusion of the bosons enables the study of entanglement generation (in the form of spin squeezing) even when the bosons are driven strongly and near-resonantly, and thus provides a quantitative view of the breakdown of adiabatic elimination that inevitably occurs as one pushes towards the fastest entanglement generation possible. The solution also helps elucidate the effect of finite temperature on spin squeezing. The model considered is relevant to a variety of atomic, molecular, and optical systems, such as atoms in cavities or trapped ions. As an explicit example, the results are used to quantify phonon effects in trapped ion quantum simulators, which are expected to become increasingly important as these experiments push towards larger numbers of ions.


I. INTRODUCTION
Spin-spin interactions play a crucial role in the generation of entanglement for applications in quantum information and metrology [1]. In atomic, molecular, and optical (AMO) systems, intrinsic spin-spin couplings are often extremely weak, and generating entanglement much faster than decoherence timescales remains an important and challenging task. One strategy to realize strong, long-ranged spin couplings, which is routinely employed in both trapped-ion systems and cavity QED, is to mediate them via a collection of auxiliary bosonic degrees of freedom (e. g. phonons in the case of trapped ions [2] and photons in cavity QED [3]). If these bosonic modes are far off-resonance and the temperature is sufficiently low, they are only virtually occupied and can be (perturbatively) adiabatically eliminated [4]. This procedure yields approximate spin-only models that are generally easier to treat theoretically, and often more desirable experimentally. For example, if ω is the characteristic energy input needed to create a boson, and g is the characteristic coupling strength between the spins and bosons, spinspin interactions of strength ∼ g 2 /ω can be generated. However, the limit in which this procedure is quantitatively valid (ω g) is directly at odds with the limit in which the spin-dynamics is fastest (large g 2 /ω). In order to overcome intrinsic timescale limitations, experiments are often forced to operate in parameter regimes where perturbative adiabatic elimination is not quantitatively justified, and a simple spin-only picture is questionable. When the phonons are far off-resonance, ω g, one can generally derive an approximate spin-only description of the system, which has direct spin-spin couplings of order ∼ g 2 /ω.
More generally, coupled spin-boson models play a central role in our understanding of quantum systems in contact with an environment, and have been studied extensively in both the condensed-matter and AMO communities for decades. Even in the case of a single spin coupled to many noninteracting bosons [5], or many spins coupled to a single bosonic mode [6,7], remarkably rich and complex behavior emerges. The general problem of many spins coupled to many bosons has very few analytically tractable limits, and is extremely difficult to study numerically, especially out-of-equilibrium and in more than one spatial dimension. As such, exact solutions -even of the simplest nontrivial models -can play an impor-tant role in extending our understanding of these intricate coupled quantum systems.
Here, we provide an exact solution for the far-fromequilibrium dynamics of a collection of spins (with S = 1/2) coupled uniaxially to a collection of noninteracting bosonic modes. The solution is valid for arbitrary spatial structure of the bosonic modes, and therefore applies to systems in any number of spatial dimensions. To a high degree of approximation, this model describes the dynamics of trapped-ion crystals when they are perturbed by a spin-dependent force [8]. When the coupled spinphonon system is driven far off (phonon) resonance, the phonons can be adiabatically eliminated and the dynamics is governed by an Ising Hamiltonian acting only on the spins. Several experimental groups have exploited this result to engineer spin-entangled states of trapped ions [9]. However, as system sizes increase, it is crucial to characterize and understand the discrepancies from this idealized situation that arise from finite population of the phonons, either due to their nonzero initial temperature or due to deviations from the far-off-resonance limit. As a demonstration of its utility, the solution is used to calculate the spin squeezing generated dynamically by initializing the system in a product state of the spin and phonon degrees of freedom that is far from equilibrium. The solution yields expressions that are efficient to evaluate numerically, enabling the calculation of dynamics for most experimentally achievable system sizes (N 10 3 spins and phonons).
The organization of the paper is as follows. In Sec. II we present the model, and review its realization with trapped ions in a simple context. In Sec. III A, we explain the formalism used to derive the exact results for this system. Details for how to explicitly calculate correlation functions are presented in Sec. III B, and numerical results for spin squeezing follow in Sec. IV. In Sec. V we discuss several interesting directions for future research.

II. MODEL AND ITS REALIZATION WITH IONS
The model we solve consists of a collection of N b bosonic modes coupled uniaxially to N s spins. The spinboson couplings can be time-dependent, and it is useful to break the Hamiltonian up into static and time varying parts as H(t) = H 0 + V(t) [10], with Here,â † α (â α ) creates (annihilates) a boson in a particular mode α, andσ r j are the (r = x, y, z) Pauli spin matrices for the j th spin. The boson energies ω α in H 0 are arbitrary, as are the coupling constants g α j (t) in V(t) (the overbar onḡ α j (t) denotes complex conjugation). A coupling to longitudinal fields ∼ j h jσ z j could also be included in H 0 , but such a term can be removed by working in a suitably rotating frame, and so it is ignored from the outset. Also, note that terms coupling to spin directions other than z are not included, and in general prohibit an exact solution.
For time-independent couplings (or, alternatively, at fixed t), the eigenstates of the above Hamiltonian are product states between all spins and suitably displaced vacuum states of the bosonic modes, and hence equilibrium properties of the model are essentially classical. However, we are concerned with the response of the system when driven out of equilibrium; the ensuing relaxation dynamics is highly nontrivial, generically being accompanied by entanglement growth between the spins and bosons. If the couplings g α j are independent of time, and if the interactions are weak (g α j ω α ), then V can be treated perturbatively. For a system initialized in the boson vacuum, the bosons will only be populated virtually in the dynamics, which can therefore be described by an effective time-independent spin-only Hamiltonian. For example, working to second order in V, one obtains Even at this level of approximation, the spin dynamics is nontrivial, and exact time-dependent correlation functions were only recently obtained for general coupling constants J jk [11][12][13]. If the J jk do not depend on space, then Eq. (2) reduces to the single-axis-twisting model [14], which is a special case of the more general Lipkin-Meshkov-Glick model [15]. In this case, the analysis is greatly simplified because the square of the total spin becomes a good quantum number, and the model can be solved in terms of collective spin variables [6]. Before solving for the time dependence of correlation functions induced by H(t), recall that Eq. (1) appears naturally in the description of various AMO systems. For example, in ion traps the spin is realized by some internal structure of an ionized atom, and the bosons are excitations of the vibrational modes of the crystalized ions (phonons). In cavity QED, where identical [16] or closely related [17] models can be realized, the spins are two-level neutral atoms and the bosons are photons in long-lived cavity modes. In the context of trapped ions, the model can emerge in several different ways [8], the conceptually simplest of which is through the application of a spin-dependent optical force to a crystal of ions [18][19][20] (though see Refs. [21,22] for a common alternative realization). For example, in the spirit of Ref. [23], the ions can be driven by two lasers with difference frequency µ and relative wave-vector k rel , as in Fig. 2a. Each ion is assumed to possess two long-lived hyperfine states labeled | ↑ and | ↓ ; they will represent the spin degree of freedom. If the energy splitting between these two states is ∆, then the ion Hamiltonian in the absence of driving is For simplicity, we assume that the crystal possesses a direction along which a single set of decoupled normal modes oscillate, the z-direction, and that k rel points along this direction; the index α in Eq. (3) enumerates this set of modes. We also assume that the laser couples both spin states to a single optically-excited state |e , and choose the laser frequency so that the detunings from the optical transition are equal in magnitude and opposite in sign (±∆/2) for the two spin states (see Fig. 2b). If the single-photon Rabi frequency Ω 0 (assumed to be the same for both the |↑ , |↓ ↔ |e transitions) is small compared to the single-photon detuning ∆/2, the electronic excited state can be adiabatically eliminated, leaving behind an AC-Stark shift for each spin state that oscillates in time at the difference frequency µ and in space at the difference wave vector k rel . Combined with a rotatingwave approximation (i. e. ignoring all terms with opticalfrequency time-dependences) and a frame transformation to remove the energy splitting ∆, adiabatic elimination of |e yields Here, Ω ≡ 4Ω 2 0 /∆ is the characteristic strength of the spin-dependent AC-Stark shift experienced by the states |↑ and |↓ , and k rel = |k rel |. The position operator along the z-direction for the j th ion, denotedẑ j , can be expanded in terms of creation and annihilation operators for the normal modes of the crystal as k relẑj = α η α b α j (â † α +â α ). Here, b α j is the orthogonal normalmode transformation matrix, and η α = k rel /(2mω α ) (restoring temporarily) parametrizes how small the characteristic ion displacements in the ground state of the mode α are compared to the length scale k −1 rel over which the applied spin-dependent potential changes appreciably. In the Lamb-Dicke limit, η α 1 for all α, and working to lowest order in η α , Eq. (4) becomes + Ω sin(µt) which is Eq. (1) with g α j (t) =ḡ α j (t) = Ω sin(µt)η α b α j . Having motivated the general form of the Hamiltonian in Eq. (1), we now proceed to compute correlation functions evolving under it. With the formal solution in hand, however, we will eventually return to the context of trapped ions and Eq. (5) when discussing the application of our results to computing spin squeezing in Sec. IV.

III. SOLUTION FOR CORRELATION FUNCTIONS
The following section includes technical derivations that are not essential for following most of the discussion in Sec. IV; readers wishing to skip these details can proceed directly to that section. Because V(t) in Eq. (1) is explicitly time-dependent, the time-evolution operator corresponding to H(t) must be written as a timeordered product, which complicates the calculation of observables. The first step in obtaining closed forms for correlation functions, therefore, is to obtain an explicit form of the time evolution operator that does not require time ordering. It is well known (c.f. Refs. [24,25]) that this can be accomplished via appropriate factorizations of the time-evolution operator. However, in the interest of maintaining a self-contained solution of the model, this procedure is briefly reviewed in Sec. III A. With an explicit form of the time-evolution operator in hand, we then move on to our main formal results in Sec. III B, obtaining closed-form expressions for spin-spin correlation functions.

A. Explicit form for the time-evolution operator
The time-evolution operator satisfies the equation of motion i∂U(t, t 0 )/∂t = H(t)U(t, t 0 ) with respect to the full Hamiltonian H(t) defined in Eq. (1), and can be written as a time-ordered product, with T t the time-ordering operator. The first step in rewriting the time-evolution operator without the need for time ordering is to move to the interaction picture with respect to H 0 . Defining the "perturbation" in the interaction picture, Here H 0 (t − t 0 ) is the product of the "unperturbed" Hamiltonian H 0 and the time difference t − t 0 ; the latter should not be confused as an argument of H 0 , which is manifestly time-independent. The next step follows the textbook problem of driven harmonic oscillators [26,27]. Defining the operator The benefit of this factorization becomes immediately clear upon differentiatingŨ(t, t 0 ) with respect to t. We take the derivative with the help of the equality which can be verified by Taylor-expanding both sides. Crucially, the commutator of W(t 2 , t 0 ) with V I (t 1 , t 0 ) depends only on the operatorsσ z j , and hence commutes with both W and V I at all other times, As an immediate consequence of Eqs. (12) and (13), one can make the replacement in Eq. (11), evaluate the right-hand-side, and thereby obtain Because Eqs. (12) and (13) hold for all times, Eq. (14) can be integrated without regard for time ordering, yield- The full time-evolution operator can now be written as At this point we have reduced the evaluation of a timeordered product to the evaluation of the product of three different time-evolution terms, each generated by an operator that commutes with itself at different times [but note that in general the order of the three exponential factors in Eq. (15) must be maintained]. It turns out that for the choice of Hamiltonian in Eq. (1), only the first and second factor do not commute with each other, and need to maintain their relative ordering. The operator W(t, t 0 ) can be written explicitly as Taking the commutator in the third factor of Eq. (15) and integrating yields Note that we have utilized [σ z j ,σ z k ] = 0 to write the Ising coefficients in an explicitly symmetric form, so S jk = S kj .

B. Calculating time-dependent expectation values of spin operators
To simplify the notation in what follows, we set t 0 = 0 and suppress its appearance, in which case the timedependent expectation value of an operatorÔ is given by with |ψ 0 the initial wave function at t = 0 and the evolution operator given in Eq. (18). We will only consider expectation values of spin operators; because they commute with the bosonic Hamiltonian H 0 , the boson-only part of the evolution operator always cancels out, leaving The results that follow can easily be worked out for arbitrary product states between all of the degrees of freedom (spin and boson). However, in order to simplify the discussion we present results only for initial states where all spins initially point along the +x axis, The bosons are taken to be in a product state between the different modes, but we allow the state of any particular mode, |ϕ α , to be arbitrary. While Eq. (21) may seem restrictive, it is a natural choice for the generation of spin squeezing, as will be made clear in Sec. IV.
To characterize the spin dynamics of the system, we compute a number of time-dependent expectation values of products of spin operators. Defining spin raising (+) and lowering (- with a, b = ± and m = n. Note that all nontrivial expectation values of one or two spin operators can be obtained as linear combinations of those in Eqs. (22)(23)(24) (correlation functions involving only theσ z j are independent of time, as [H(t),σ z j ] = 0). In order to calculate these expectation values, we make a few observations about Eq. (20). First, the final factor of the time-evolution operator, which involves only the Ising spin operators, can be further factorized into a product of exponentials (instead of an exponential of a sum of terms), because each term in the exponent commutes with every other term in the exponent. Second, we can factorize the exponential of W(t) into factors that have fixed α (but still sum over j), because those operators also commute with each other. Each one of these factors that commutes with the operatorÔ can be moved from the right, throughÔ, and then cancels against an inverse term on the left coming from hermitian conjugate of the time-evolution operator.
Using the identity , which follow from standard operator identities, we can simplify the expectation values of the operator averages that we are calculating. After some additional algebra, and defining the displacement operatorsD α (ϑ) = exp(ϑâ † α −θâ α ) and modified spin-spin couplings the final results are as follows: From Eq. (17), we see that . This situation is realized for the normal modes of ions in linear Paul traps, and also for the axial modes of Penning traps (though not the in-plane modes) [28]. There,ḡ α j (t) = g α j (t) = Ω sin(µt)η α b α j , since the normal-mode transformation matrix b α j can always be chosen to be real.

C. Evaluation of the boson matrix elements
A boson matrix element of the form ϕ α |D α (ϑ) |ϕ α can easily be evaluated for an arbitrary state ized Fock states of the α th mode. Writing the displacement operator asD α (ϑ) = e −|ϑ| 2 /2 e ϑâ † α e −θâ α , straightforward algebra leads to The complete boson matrix elements in Eqs. (28)(29)(30) follow by taking the product of Eq. (31) over all of the N b modes labeled by α. Expectation values such as Eq. (31) can easily be generalized to deal with finite temperature states of the bosonic modes. For example, in what follows we consider a situation where all bosonic modes are thermally populated at an inverse temperature β, such that the initial boson density matrix is (β) = α ρ α (β), with Expectation values of the form ϕ α |D α (ϑ α ) |ϕ α appearing in Eqs. (28)(29)(30) should be replaced with D α (β, ϑ α ) ≡ Tr ρ α (β)D α (ϑ α ) . Inserting Eq. (32) into this trace and utilizing Eq. (31), straightforward algebra leads to

IV. SPIN SQUEEZING
It is also possible to obtain closed-form expressions for higher-order (n-spin) versions of the 2-spin correlation functions derived above. However, already at the level of 2-spin correlation functions, we can learn a great deal about the time evolution of the system, and the nature of the entanglement that develops. For example, from the correlation functions in Eqs. (29) and (30), we can compute the variance of the spin distribution, which enables us to characterize spin squeezing [14,29]. Spin squeezing is just one of many measures of entanglement in a manybody system, which has the virtue of quantifying the potential enhancement in precision obtainable in Ramsey spectroscopy (as compared to the case of unentangled spins) [30]. Moreover, it establishes a lower bound on the depth of entanglement, i. e. the minimum number of simultaneously entangled particles in the system [31].
If the bosonic modes are initially cooled to a vacuum state, driven weakly so that they can be adiabatically eliminated, and if the resulting effective spin-spin coupling strength J is independent of the spatial distance between the spins, the dynamics is governed by the single-axis-twisting Hamiltonian H sat = 4JŜ zŜz , wherê S z ≡ (1/2) jσ z j [14]. In this model, spin squeezing is generated by first polarizing the collective spin vector along the +x axis, and then letting it evolve under H sat . Making a mean-field approximation, H sat ≈ 8J Ŝ z Ŝ z − 4J Ŝ z 2 , the dynamics can be understood as a precession about the z-axis in a direction determined by the mean z-component of the spin; the initial spin state has quantum fluctuations above and below the equator, and therefore this dynamics causes it to get sheared and elongated, as in Fig. IV. Uncertainty along one axis is reduced (squeezed), while uncertainty in an orthogonal direction is increased.
The extent of squeezing along a particular direction in the yz plane can be quantified by the parameter whereŜ θ ≡Ŝ z cos θ+Ŝ y sin θ and ∆S θ = ( Ŝ 2 θ − Ŝ θ 2 ) 1 2 . The spin-squeezing parameter is then defined by minimizing the standard deviation, ∆S min = min θ ∆S θ , such that ξ = min θ ξ(θ) = N 1/2 s ∆S min /| Ŝ x |. Straightforward algebra enables the optimal angle to be expressed explicitly in terms of spin-spin correlation functions [29], A. Connection to trapped ions Equations (28-30) are very general, enabling a complete description of spin-spin correlations in a variety of different physical systems, but they cannot be further simplified or evaluated without choosing explicit forms for the boson spectrum ω α and spin-boson couplings g α j (t). In the remainder of Sec. IV, we compute the squeezing induced by the Hamiltonian in Eq. (5) using parameters relevant for ions in a linear Paul trap, though many of the qualitative features discussed below are insensitive to the detailed trap geometry, and would be similar for ions in a Penning trap if the axial modes were being driven [23]. Specifically, we examine the dynamics of 20 ions, and assume that the wave-vector difference of the driving lasers (k rel ) points along one of the two transverse principal axes of the trap. In this configuration, the spins only couple to normal modes oscillating along a single spatial direction, and therefore the number of (coupled) phonon modes is equal to the number of spins; N ≡ N s = N b = 20. To calculate the normal-mode eigenvectors and associated frequencies, we set the ratio of longitudinal to transverse trap frequencies equal to 0.1. The mass of the ions enters the calculation only as an overall scaling of the normal-mode frequencies, and can be ignored if all energies (times) are measured in units of the center-of-mass frequency ω com (1/ω com ). The phonon modes are calculated in the standard fashion assuming a pseudopotential approximation holds [28,32,33], which neglects any effects due to micromotion. We first find the equilibrium positions of the N ions in the trap for the given trap curvatures in the different spatial directions and the mutual Coulomb repulsion of the ions. Then we determine the spring constant matrices by expanding the Coulomb interaction to quadratic order about equilibrium. Diagonalizing these dynamical matrices yields both the normal-mode oscillation frequencies ω α and the associated orthonormal eigenvectors b α j . Because the spring-constant matrices are real and symmetric, the b α j are real, and the spin-phonon coupling constants satisfy g α j (t) =ḡ α j (t) = Ω sin(µt)η α b α j . Substituting this expression for g α j (t) into Eqs. (17) and (19) and performing the integrals, we obtain Here A α j (t) is proportional to the interaction-picture phase space displacement of the j th spin as a result of periodically driving the mode α. Because this driving is periodic, the displacement amplitudes A α j (t) have a simple structure. In particular, for a single mode driven near This amplitude traces a closed circle, vanishing at times such that δ α t = 2nπ (Fig. 4a), with n an integer, at which the phonon matrix elements of the mode α in Eqs. (28)(29)(30) all become equal to unity. At these times, the mode α becomes unentangled with the spins, regardless of its initial state. Physically, the independence on initial state reflects the independence of the period of a harmonic oscillator on its displacement, and implies that the periodic disentanglement of a given phonon mode occurs even at finite temperature, which was pointed out originally in Ref. [21]. Formally, it is immediately apparent from Eq. (34) that finite-temperature phonon expectation values of the form D α β, A α j (t) will return to unity periodically for any inverse temperature β (see Fig. 4b).
However, when multiple modes participate in the dynamics, they do not all disentangle at the same time, with further off-resonant modes exhibiting phase space excursions with smaller radii but larger frequencies, as in Fig. 4c.
In Sec. IV B, we produce plots of the squeezing parameter ξ as a function of time using the exact solution and also using two useful approximations. In the first approximation, we evolve our initial state with the spin-only The phonon modes periodically disentangle at finite temperature, as can be seen from the expectation values Dα β, A α j (t) returning to unity. The overlap of a given mode with its initial value, however, decreases more sharply away from the recurrence times at higher temperature (here kBT = ωα × {0, 2, 5, 10}, from darker to lighter curves). (c) When multiple modes are driven, those further from resonance make smaller but faster excursions through phase space, and in general the different modes do not simultaneously return to their initial states.

time evolution operator
which amounts to ignoring spin-phonon entanglement. This evolution is achieved by replacing A α j (t) → 0 in Eqs. (28)(29)(30) while treating the S jm (t) [which in reality depend implicitly on the A α j (t)] as independent parameters, and then evaluating the S jk (t) using Eq. (37). In the second approximation, spin-motion entanglement is ignored and the spin-spin couplings are replaced with their time averages, yielding time evolution under a time-independent Ising spin model with time-averaged coupling constants We produce plots by varying Ω, which controls the spinphonon coupling strength, δ ≡ µ − ω com , which is the difference between the two-photon detuning of the Raman lasers and the center-of-mass mode frequency, and also the phonon temperature T , which controls the initial number distribution in the phonon modes.

B. Phonon effects at T = 0
Trapped ion experiments aimed at quantum simulations generally cool the phonon modes to temperatures on the order of, or ideally lower than, typical phonon energies. However, in general, and especially for 2D or large 1D crystals, these temperatures are decidedly not zero for all of the relevant modes. Nevertheless, the phonons are often assumed to be cooled well enough that one can approximate them as being at zero temperature; we will address the zero-temperature situation first, and later come back to analyze the consequences of having thermal phonons in the initial state. At zero temperature, the qualitative structure of the dynamics is determined by two independent considerations. First, the coupling strength to a given phonon mode, measured relative to its detuning from the drive frequency µ, determines how strongly that mode is driven. Because the vibrations of the particles in the trap are transverse, all vibrational modes have frequency ω α ≤ ω com . As a result, by choosing δ ≥ 0 (so that no modes other than the center-of-mass mode can be resonant), we ensure that Ω/δ is a suitable measure of how deeply in the "perturbative" limit the system is. When this ratio is small, all phonon modes are weakly populated in the dynamics, and we expect spinphonon entanglement to be unimportant. Conversely, when this ratio is large, spin-phonon entanglement will be important, and we do not expect the approximations described above to agree well with the exact solution. Next, we must consider the absolute size of δ and Ω relative to the mode bandwidth, denoted by D, which controls the relative extent to which different phonon modes participate in the dynamics. For example, when δ and Ω are small compared to the typical mode spacing, δ, Ω D/N , spin dynamics occurs primarily due to coupling to the center-of-mass mode. By increasing δ from much smaller than the typical mode spacing to much larger than the mode bandwidth, all the while keeping Ω δ, one can navigate between two extreme limits: (A) For δ D/N , the center-of-mass mode dominates the mediation of spin-spin interactions, which are therefore independent of the distance between two spins (since the center of mass mode is spatially homogeneous). (B) For δ D, all modes participate equally in mediating the spin-spin interactions, which fall off roughly as the cube of the distance between two spins [4]. In between these extremes, it is common to approximate the spin-spin interaction to decay as a power law, S avg jk ∝ 1/r ζ , with r the distance between ions j and k and 0 < ζ < 3. All of these considerations are summarized in Fig. 5, where a guide to the parameter space explored in the rest of this section can also be found. Figure 6 illustrates the behavior of the squeezing parameter ξ at T = 0, with the detuning chosen close to the center-of-mass mode (δ D/N ) and marginally in the perturbative limit (Ω/δ = 1/4). The squeezing parameter is normalized such that log(ξ) = 0 for a coherent state, so the region of the plot where the squeezing parameter dips below the horizontal axis denotes a period of improved uncertainty with respect to the standard quantum limit. Under these conditions, the exact solution and the two approximations exhibit fairly similar behavior. The smooth red curve produced by the timeaveraged spin-coupling approximation correctly captures Note that the time-axis scaling is different in each plot, and the optimal squeezing occurs progressively later as δ is increased. The color code here is the same as that used in Fig. 6.  8. (Color online) At zero temperature, the detuning from the center-of-mass mode is varied while keeping Ω/δ = 1/4, thereby controlling the relative contribution of the different phonon modes without greatly affecting the extent to which the individual modes are in the perturbative limit (i.e. Ω/(µ − ωα) is not changing very much). Note that as δ grows larger, the time of optimal squeezing becomes shorter and the total amount of squeezing obtained at that time is shrinking. The color code here is the same as that used in Fig. 6.
the overall trend of squeezing observed in the exact solution. As seen in the black curve, the time-dependence of the spin couplings produces small amplitude highfrequency oscillations. These oscillations are amplified by spin-phonon entanglement, as indicated by the exact solution (blue curve), showing that even at zero temperature and nominally in the perturbative limit, the creation of phonons during the dynamics can significantly affect the spin squeezing. We note that very slight improvements in squeezing over the time-averaged spin-coupling approximation do periodically occur. They appear to be due to the time dependence of the true spin couplings S jk (t), rather than the spin-phonon entanglement, as they occur in both the blue and black curves (and the latter bounds the former from below).
To better understand the role that dynamical phonon creation plays in spin squeezing, we change the detuning from the center-of-mass mode, δ, while holding Ω constant. Figure 7 illustrates a series of results in which δ is increased but kept small compared to the detuning from all other modes (i.e. δ D/N ). These plots therefore primarily reflect changes in behavior caused by variations in the strength with which the COM mode is driven relative to its detuning, i.e. the ratio Ω/δ. In the time-averaged spin-coupling approximation (red curve), the dynamics is completely insensitive to this ratio up to an overall timescale, and upon scaling the maximum time by δ/Ω, we observe nearly the same squeezing behavior in all three panels. In the first plot, δ/Ω is relatively small; the center-of-mass mode is strongly driven, resulting in strongly oscillatory spin couplings and large spin-phonon entanglement. In this limit, neither approximation accurately captures the exact dynamics, nor do they agree with each other. As δ is increased, the oscillation amplitudes of the time-dependent couplings S jk (t) diminish compared to their time-averaged values, so the two approximations begin to agree better with each other. At the same time, population in the center-of-mass mode becomes suppressed, spin-motion entanglement diminishes, and the exact result begins to converge to both approximations. As discussed in Sec. IV A, when considering coupling only to the center-of-mass mode, the phonon degree of freedom periodically becomes unentangled from the spins, even at strong driving. This behavior is reflected in the periodic agreement between the exact solution and the approximation ignoring spin-phonon entanglement.
In all panels of Fig. 7, both δ and Ω are kept small compared to the mode spacing, and therefore the variation of Ω/δ is the dominant factor contributing to the changes in behavior. However, these plots provide little insight into the other important effect of increasing the detuning δ: the increased importance of modes other than the center-of-mass mode. In order to isolate the latter effect, in Fig. 8 we again vary δ but now keep the ratio Ω/δ = 1/4 fixed. This keeps the coupling to any given mode in the (barely) perturbative limit, thus counteracting the dominant role that varying spin-phonon entanglement played in the qualitative trends observed in Fig. 7. As δ is increased, the time-averaged spin-coupling approximation correctly captures an overall trend of the exact solution: The time of optimal spin squeezing is becoming shorter, and the extent of squeezing that occurs at that time is diminishing. The former effect is simply the result of increasing the spin-phonon coupling, which increases the overall rate of dynamics. The reduction of squeezing achieved at the optimal time reflects the diminishing spatial-range of the spin-spin couplings (approaching 1/r 3 for large δ) as more phonon modes participate in mediating them. With or without the inclusion of spinphonon entanglement (i.e. in the blue or black curves), the squeezing exhibits high-frequency oscillations arising from the coupling to phonon modes other than the centerof-mass mode. As multiple phonon modes at different frequencies become entangled with the spins, it is no longer possible for all of those modes to become disentangled from the spins simultaneously. This is strongly reflected in the third panel, where the exact solution no longer agrees with either approximation at regular intervals.

C. Effects of finite temperature
In Sec. IV B, we were primarily interested in the effects of dynamical phonon generation starting from the phonon vacuum. In many experiments, especially those using large numbers of ions, this starting state is not realistic. For example, in Ref. [23], the Doppler-cooling limit of T ∼ 1mK 10 × ( ω com /k B ) leads to 10 phonons per transverse mode. As explained in Sec. III C, spinspin correlation functions can also be computed at finite temperatures, and the above analysis can therefore be extended to capture the consequences of non-zero initial motional temperature on spin squeezing.
With the addition of temperature, there is a large parameter space to be explored; here we focus our attention on the barely perturbative limit at fixed Ω/δ = 1/4 (the same ratio used in Fig. 8), and consider both variations of δ and T (see Fig. 9 for a guide to the parameter space explored). We first examine the case of near-detuning from the center-of-mass mode (δ, Ω D/N ), but now taking the phonon modes to be at a temperature T at time t = 0 (Fig. 10). Here we plot just the exact solution and the time-averaged spin-coupling approximation (both the time-averaged spin-coupling approximation and the approximation of ignoring spin-phonon entanglement are insensitive to the phonon distribution, so neither varies with T ). The primary effect of increasing temperature on log log ⌦ log T  FIG. 9. (Color online) Schematic illustration of a twodimensional cross-section of the full three-dimensional parameter space spanned by Ω, δ, and T . Working just barely in the perturbative limit (and at a fixed value of Ω/δ), we explore the effects of temperature for multiple values of δ, which controls the relative participation of the various phonon modes in mediating spin-spin interactions. The parameter space explored in Fig. 10 and Fig. 11 are indicated as vertical arrows. the squeezing is that the amplitude of oscillations above the curve obtained from the time-averaged spin-coupling approximation increases. Nevertheless, the exact solution continues to agree with the approximation at regular intervals. As discussed in Sec. IV A, this behavior can be understood as the insensitivity of a harmonic oscillator's period to its state of excitation; at finite temperature many Fock states of the center-of-mass mode are occupied, but as they are driven periodically they all return to their initial point in phase space simultaneously, at which point they are unentangled from the spins. As T becomes larger, the system spends most of its time in states with large uncertainty (poor squeez- , finite temperature has a much more detrimental effect on the spin squeezing, eventually preventing any squeezing from occurring beyond a critical temperature (here T ≈ 10). Inset: Note that, due to the relatively large population of many modes, there is significant high frequency structure that is hidden on the time scale over which squeezing occurs.
ing), but precisely timed measurements of the system could nevertheless yield a significantly improved resolution over the standard quantum limit. Interestingly, at times short compared to 2π/δ, the spin distribution is always antisqueezed, and the degree of antisqueezing could in principle be used to perform in-situ temperature measurements of the phonon modes. In Fig. 11, we explore the behavior of the system detuned away from the center-of-mass frequency and at finite temperature. The behavior exhibited in Fig. 11 reflects the general trends observed in the previous plots: increasing δ induces additional high-frequency structure in the exact solution, and increasing T produces an overall growth in the amplitude of oscillations in the spin squeezing. Unlike in Fig. 10, however, the exact solution not only shows increased oscillation amplitude at higher temperature, but in addition the local squeezing minima are also displaced increasingly far from the curve calculated in the time-averaged spin-coupling approximation. Even though at T = 0 the squeezing nearly agrees with this approximation throughout the entire time region plotted, spin squeezing is completely destroyed when the temperature reaches T ≈ 10.
Indeed, for a given detuning there will always be some temperature threshold above which no squeezing takes place at all (ξ ≥ 1 for all t). By varying T at evenly spaced intervals of log δ, in Fig. 12 we produce a "phase diagram" demarcating this boundary in the parameter space. The qualitative downward trend of the boundary can be understood as an increased sensitivity of the spin dynamics to the initial phonon temperature as more phonon modes participate in the dynamics. The additional modes prevent the spins from becoming periodically unentangled with the phonons; the consequences of this residual spin-motion entanglement on spin squeezing showing that no squeezing occurs above a critical temperature (note that Ω/δ is being held constant as δ is changed, see Fig. 9).
are exacerbated at finite temperature because the occupation of a given phonon mode affects the extent to which that mode remains entangled with the spins at a given time (see Fig. 4b). At small δ, squeezing persists even for temperatures on the order of 100 times the centerof-mass-mode energy. However, note that at these high temperatures the approximation that the phonons are non-interacting (equivalently that the potential seen by an individual ion is harmonic) is likely to break down, and the system may well be outside of the Lamb-Dicke limit used to justify a description in terms of Eq. (1).

V. DISCUSSIONS AND FUTURE DIRECTIONS
While Eq. (1) is quite general, and can describe a diverse array of practically relevant experiments, there is certainly a sense in which it is extremely constrained-it has a large number of locally conserved quantities. Indeed, a careful inspection of our calculations reveals that the local conservation ofσ z j is, at several different levels, directly responsible for the solvability of the model. Nevertheless, our solutions form a useful benchmark for powerful (but computationally expensive) tools capable of numerically solving the non-equilibrium dynamics of more general models, such as time-dependent versions of matrix-product-state algorithms [34]. Especially in 2D-but generally even in 1D given the long-range nature of the spin-couplings induced by the (often) delocalized bosonic modes-it remains unclear to what extent there are any general purpose and controlled numerical tools for studying this dynamics. We also point out that questions about equilibration and thermalization, which require studies of long-time dynamics, are difficult even in 1D for all but the smallest systems.
The exact solution will be useful in testing a variety of experimental idealizations that are frequently made but generally not quantitatively justified. For example, many trapped-ion experiments employ a spin-echo pulse in or-der to obtain a coherence time on the order of the spinspin interaction time. If the spins are unentangled with the phonons, then a spin echo completely removes the effects of an inhomogeneous magnetic field (∼ j h jσ z j ) added to the Hamiltonian in Eq. (1). However, the existence of spin-phonon entanglement at the time of the echo pulse invalidates this picture [35], and the consequences on spin-spin correlation functions could be explored using the solution developed in this paper.
There are also many interesting purely theoretical questions about the present model, many of which we believe the tools developed here are well suited to begin answering. For example, while spin-squeezing calculations reveal the impact of spin-boson entanglement on attempts to generate strictly spin-spin entanglement, it should also be possible to analyze spin-boson entanglement more directly by calculating correlation functions involving both spin and boson operators. It would also be interesting to explore to what extent the "solvability" of this model can be extended to calculating more general quantities than low-order correlation functions. While our solution enables efficient calculation of arbitrary-order correlation functions of the form given in Eqs. (22)(23)(24), calculating the full counting statistics along any spin direction orthogonal to z is still very difficult. In particular, computation of the full counting-statistics is equivalent to the computation of ∼ N th -order correlation functions of Pauli (x, y, z) matrices, which involves the summation of an exponentially large (in N ) number of high-order correlation functions of the form in Eqs. (22)(23)(24). It seems very plausible that this difficulty is related to computational hardness results for classically sampling spin distributions following dynamics under commuting spin Hamiltonians [36]. The present model adds an interesting twist, in that it builds bosonic degrees of freedom into a commuting spin Hamiltonian in a way that preserves its solvability (in the sense of calculating loworder correlation functions); the bosons alone, despite being non-interacting and therefore "solvable" in the same sense as the model studied here, are thought to be hard to simulate clasically in a precise sense [37].