# Simulating generic spin-boson models with matrix product states

###### Abstract

The global coupling of few-level quantum systems (“spins”) to a discrete set of bosonic modes is a key ingredient for many applications in quantum science, including large-scale entanglement generation, quantum simulation of the dynamics of long-range interacting spin models, and hybrid platforms for force and spin sensing. We present a general numerical framework for treating the out-of-equilibrium dynamics of such models based on matrix product states. Our approach applies for generic spin-boson systems: it treats any spatial and operator dependence of the two-body spin-boson coupling and places no restrictions on relative energy scales. We show that the full counting statistics of collective spin measurements and infidelity of quantum simulation due to spin-boson entanglement, both of which are difficult to obtain by other techniques, are readily calculable in our approach. We benchmark our method using a recently developed exact solution for a particular spin-boson coupling relevant to trapped ion quantum simulators. Finally, we show how decoherence can be incorporated within our framework using the method of quantum trajectories, and study the dynamics of an open-system spin-boson model with spatially non-uniform spin-boson coupling relevant for trapped atomic ion crystals in the presence of molecular ion impurities.

## I Introduction

The coupling of spins to bosonic modes provides a paradigmatic setting for various important phenomena in many-body physics, such as decoherence, thermalization, virtually mediated long-range spin-spin interactions, and large-scale entanglement generation Leroux et al. (2010); Bohnet et al. (2014); Hosten et al. (2016). The study of spin-boson models has historically belonged to the domain of condensed matter theory Caldeira and Leggett (1983). However, experimental setups available to atomic, molecular and optical (AMO) physics provide versatile, clean, and controllable realizations of such models. In particular, spin-boson models can be engineered in trapped ion systems by coupling the discrete phonon modes of the ion crystal to the ion spins Sørensen and Mølmer (1999); Porras and Cirac (2004); Kim et al. (2009); Britton et al. (2012); Porras et al. (2008) and in cavities via cold atoms Ritsch et al. (2013); Leroux et al. (2010) or artificial atoms built from superconducting qubits Niemczyk et al. (2010); Hur et al. (2015) coupled to quantized cavity modes; similar models can also be achieved in optomechanical setups Aspelmeyer et al. (2014). While these AMO platforms differ from the most common condensed matter spin-boson analog–a two-state version of the Caldeira-Leggett model of a particle interacting with a continuum of harmonic oscillators–due to their discrete boson mode spectra, they offer clean systems in which one can observe out-of-equilibrium dynamics.

The long history of spin-boson models has led to the development of a range of theoretical approaches. The most prevalent one is to adiabatically eliminate the bosons, keeping only their virtual effect on the spins in the form of long-range spin-spin couplings Porras and Cirac (2004). While this approach is useful in contexts where there is a very large separation of energy scales between the bosons and spins, in many experimental realizations and applications such an approach is invalid. An additional powerful approach, valid in the case that only uniform, collective spin operators appear in the Hamiltonian, uses the permutational symmetry of the density matrix to reduce the computational scaling of direct numerical diagonalization from exponential to polynomial in the number of spins Carmichael (1999); Hartmann (2012); Richter et al. (2015). Outside of the uniform coupling regime, analytical results have recently been obtained in a special case in which the coupling Hamiltonian is linear in boson operators and commutes with the spin Hamiltonian Dylewsky et al. (2016). However, in spite of the importance of the many-spin, many boson problem to myriad fields of physics, a general, systematic approach to study their quantum dynamics has yet to be developed.

In what follows, we show that spin-boson models with global coupling can be systematically treated with well-characterized error using the framework of matrix product states (MPSs). Our methodology consists of two complementary methods for treating the out-of-equilibrium dynamics. The first, using the technology of swap gates Shi et al. (2006) to deform the MPS topology dynamically, is most efficient in the case of few bosons or at weak coupling, but has a rigorously bounded, systematically correctable error. The second uses an exact matrix product operator (MPO) representation of the spin-boson Hamiltonian, and time-evolves under this Hamiltonian using recently developed techniques for time evolution with long-range couplings Zaletel et al. (2015), together with safeguards against variational metastability. This latter method is most efficient when the number of bosonic modes is equal to the number of spins, relevant for trapped ion quantum simulators in certain cases.

This paper is organized as follows: Sec. II presents our numerical formalism, which includes a brief review of matrix product states (MPSs) and operators as well as our new MPS-based methodologies for spin-boson systems; Sec. III discusses benchmark calculations of our methods for 1D and 2D systems of trapped ions in the parameter regimes of recent experiments Richerme et al. (2014); Bohnet et al. (2015); Sec. IV extends our formalism to open quantum systems through the technique of quantum trajectories; and Sec. V presents simulations of trapped ion systems in the presence of impurity ions where decoherence and spatially non-uniform spin-boson couplings invalidate many other approaches; Finally, in Sec. VI we conclude and give an outlook on broad applications of our methods. Some technical results are given as appendices.

## Ii Formalism

We consider a collection of spins localized on a lattice and coupled globally to bosonic modes, shown schematically in Fig. 1. We take the Hamiltonian to have the generic form

(1) |

where acts only on the spins, only on the bosons, and the spin-boson coupling has the form

(2) | ||||

(3) |

Here and throughout, roman letters refer to spin indices and greek letters to boson indices. The operators and can be taken to be Hermitian, and act on bosons and spins, respectively. The index counts the number of such pairs of Hermitian operators, e.g., has while has . While we will continually refer to our degrees of freedom as being spins, our approach can be straightforwardly generalized from spins to mobile particles with any quantum statistics. A key feature of this general model is that the coupling between bosons and spins is arbitrary, i.e. each boson mode can be coupled to each spin with an arbitrary amplitude and operator pairs . In contrast to approaches based on permutational symmetry of the density matrix Carmichael (1999); Hartmann (2012); Richter et al. (2015), we do not require that this coupling is uniform; our approach applies for any amplitudes and operators and . Many iconic models can be cast in the form Eq. (1), including the Rabi Gardas (2011); Braak (2011), Jaynes-Cummings Jaynes and Cummings (1963); Eberly et al. (1980), and Tavis-Cummings Agarwal et al. (2012) models of quantum optics and optomechanics Aspelmeyer et al. (2014) and the Holstein-Hubbard Jeckelmann and White (1998) model of condensed matter physics. Also, we point out that our setting is more general than previous spin-boson approaches using related density-matrix renormalization group (DMRG) methods, in which each spin was coupled locally to a single boson Jeckelmann and White (1998); Bursill et al. (1998); Brockt et al. (2015), a single spin was coupled to many bosons Guo et al. (2012); Frenzel and Plenio (2013); Schröder et al. (2015), or a single collective boson mode was coupled to many spins Gammelmark and Mølmer (2012).

### ii.1 Matrix product states and operators

Before proceeding to the details of how we simulate the out-of-equilibrium dynamics of Eq. (1), we remind the reader a few basic facts about matrix product states (MPSs) and matrix product operators (MPOs); a more detailed discussion can be found in, e.g., Schollwöck (2011); Orús (2014). The natural setting of MPSs is a 1D chain with sites, where lattice site can be in one of the states . An MPS is a representation of a many-body quantum state on such a lattice, defined as

(4) |

Here, is a dimensional matrix, and the maximum value of any of the is called the bond dimension of the MPS. MPSs are exact ground states of 1D gapped models with short-range interactions Hastings (2007). In addition, MPSs have proven to be a useful variational ansatz for long-range interacting systems Crosswhite et al. (2008); Schachenmayer et al. (2013); Wall and Carr (2012), gapless systems, e.g., quantum critical points Pollmann et al. (2009); Rams et al. (2015), and quasi-higher-dimensional systems Mishmash et al. (2011); Stoudenmire and White (2012); Depenbrock et al. (2012) by “snaking” a 1D line across the sites of the higher-dimensional lattice.

A matrix product operator (MPO) is the natural operator-valued extension of the MPS definition McCulloch (2007); Verstraete et al. (2004). It is a representation of a many-body operator acting on the same space defined for the MPS above, and has the form

(5) |

MPOs are of interest because (1) they take MPSs to MPSs (with, in general, larger bond dimension), (2) most many-body operators of interest, e.g. Hamiltonians, can be represented as MPOs with constant bond dimension by using a small set of MPO construction “rules” Crosswhite and Bacon (2008); Fröwis et al. (2010); Pirvu et al. (2010); Wall and Carr (2012), and (3) many MPS algorithms, such as variational search for the ground state and time evolution, can be formulated generically given that all many-body operators of interest have a known MPO form Wall and Carr (2012). It is useful to re-write Eq. (5) as

(6) |

where each one of the is an operator-valued matrix of the form, e. g.,

(7) |

The matrix indices are the fictitious “bond” indices contracted in the trace of Eq. (6), and the indices of the individual operators, e.g., , run over the states of the local physical Hilbert space. For simplicity, we have given a bond dimension two MPO matrix as an example, but any dimensionality is possible (see Eq. (10)). A useful means of displaying MPSs, MPOs, and their operations is the Penrose tensor network diagram notation exemplified in Fig. 2. In this notation, a tensor is an object with lines extending from it, with the number of lines equal to the rank of the tensor. A line connecting two tensors denotes that index is summed, or “contracted,” between the two tensors.

### ii.2 Swap operator decomposition of the spin-boson state and propagator

In this section, we propose a method for simulating the out-of-equilibrium dynamics of Eq. (1) which is based on Trotter-Suzuki decompositions. First, we use a Trotter-Suzuki decomposition to isolate the action of from and . The simplest non-trivial such decomposition is ( unless otherwise specified throughout)

(8) |

which is valid to order . While this is written in a form appropriate for time-independent , a time-ordered version of this formula also exists Huyghebaert and Raedt (1990), as well as higher order decompositions in Omelyan et al. (2002). Since and commute by construction, the evolutions associated with these operators can be performed in any order. In the examples given in this paper both and are sums of single-site operators and so their exponentials are product operators which are trivial to apply to an MPS. Alternatively, in the most commonly encountered case that they are short-ranged, this propagation can be performed with ordinary t-DMRG Vidal (2004); Daley et al. (2004).

We now turn to the application of to an MPS. We do so by using another Trotter-Suzuki decomposition to write this exponential as a product of two-site unitaries, each acting on a single spin and boson pair. The simplest such decomposition is

(9) |

and, as before, higher-order decompositions can be devised Sornborger and Stewart (1999). In this expression, means the product . For a fixed ordering of MPS sites, the application of Eq. (9) requires coupling of sites which are not contiguous. While this can be done with a variety of available long-range time evolution methods for MPSs Manmana et al. (2005); García-Ripoll (2006); Wall and Carr (2012); Zaletel et al. (2015); Haegeman et al. (2014), such an approach can suffer from numerical issues, e.g. getting stuck far from the variational optimum. In contrast, applying an operator which couples neighboring sites in the 1D MPS representation can be performed efficiently, without issues of variational metastability, and with a precisely controlled error Vidal (2004).

We can apply the decomposed operator Eq. (9) to an MPS using only operations on neighboring sites by inserting swap gates which interchange the positions of boson and spin in the MPS representation between different terms appearing in the product Shi et al. (2006). Swap gates have been previously used to implement MPS algorithms for periodic boundary conditions Danshita and Naidon (2009), as well as for generic long-range time evolution Stoudenmire and White (2010). In the latter case, long-range time evolution requires swaps per time step, where is the total number of lattice sites; this can be prohibitively costly. In our case, only swaps are required, which is more favorable in the often realized case that as well as the case .

To see how the swapping works in practice, we show the tensor network diagram for a full application of for and in Fig. 3. Here, the tensors in the initial MPS are circles (squares) for boson (spin) sites, , and crossed lines indicate the application of a swap gate. Motivated by the form of the decomposition Eq. (9), a natural 1D chain ordering of the spins and bosons is to begin with all of the bosons left of the spins. In order to keep track of the position of the boson sites in the MPS representation, we draw a circle containing the mode index next to index lines corresponding to bosons, showing that each boson is moved past all of the spins using the swap gates. From the diagram it is clear that the combination of Trotter-Suzuki decomposition and swap gates leads to a propagator which involves only operations on nearest-neighbor sites.

### ii.3 Exact MPO representation of the spin-boson Hamiltonian

As the number of bosons approaches the number of spins, the scaling of the swapping method described in the last section approaches , which can be prohibitively expensive. However, as we will discuss below, in this situation a direct simulation of the dynamics of the spin-boson system using a representation of the Hamiltonian with long-range couplings between spins and bosons can be efficient. Here, the ordering of the spins and bosons in the MPS representation can be chosen to minimize issues of variational metastability that can pose problems for long-range time evolution with MPSs. Before we discuss how to choose this ordering, we first show that Eq. (1) has an exact, compact MPO representation for spin-boson couplings of the form Eq. (3), irrespective of how the spins and bosons are ordered in the 1D array. For simplicity, we will write down the MPO in the case that and do not couple multiple sites and the sum over in Eq. (3) contains a single term; this is the form of the Hamiltonian relevant for trapped ion experiments, as shown in Sec. III. More complex spin and boson Hamiltonians can be formulated using known MPO constructions in the literature Crosswhite and Bacon (2008); Fröwis et al. (2010); Pirvu et al. (2010); Wall and Carr (2012), and more terms in the summation can be included by matrix direct sums McCulloch (2007). Additionally, for compactness, we consider the case of two boson modes; the MPO for many modes is also a straightforward generalization.

To begin, we write the collection of spins and bosons as a 1D chain with sites, and let the boson modes and 2 be placed at sites and , respectively, with . Then, the MPO matrices representing Eq. (1) take the form

(10) | ||||

Here, () is the identity operator in the spin (boson) space, is the spin index for lattice site , and the first (last) MPO matrix for open boundary conditions is the last row (first column) of the associated MPO matrix . Remarkably, in spite of the generality and complexity of the model, Eq. (1), the Hamiltonian admits a very simple and compact representation. This MPO is useful not only for the out-of-equilibrium dynamics we discuss in the next section, but can also be used to find eigenstates variationally Schollwöck (2011).

### ii.4 Long-range time evolution with

Many methods for time evolution of MPSs exist which are applicable to general Hamiltonians in MPO form, such as the Krylov Manmana et al. (2005); García-Ripoll (2006); Wall and Carr (2012) and the local Runge-Kutta Zaletel et al. (2015) methods, which require some form of variational optimization of an MPS, and the time-dependent variational principle Haegeman et al. (2014), which has an additional error associated with projection of the Hamiltonian onto a local subspace. While such methods are useful in many cases, the error analysis and convergence behavior of these algorithms can be complex. Of particular concern when applying these algorithms to our Hamiltonian (1) via the exact MPO representation Eq. (10) is that the effectiveness will depend strongly on the ordering of the spin and boson modes in the 1D MPS representation. This is especially true when many strongly coupled boson modes are present, as this situation can lead to glassy spin physics Gopalakrishnan et al. (2009); Graß et al. (2015) which may pose problems for variational optimization.

The main issue of variational metastability for systems with long-range interactions is that efficient MPS algorithms use an iterative local variational optimization to search for the global optimum Schollwöck (2011), where local means that two neighboring sites are optimized at a time. Hence, MPS algorithms build up entanglement and correlations locally, two sites at a time, and then build up longer-ranged correlations and entanglement by repeated “sweeping,” in which each pair of neighboring sites are optimized in a round-robin fashion. If there are long-ranged (farther than nearest-neighbor) terms in the Hamiltonian, entanglement resulting from these terms is not directly built in by the optimization, and so we must either build in the proper entanglement structure by hand or it must be incorporated in by other short-range terms in the Hamiltonian upon repeated sweeping.

For the system at hand and in the case , we can ensure that all spins and bosons have a least some short-range entanglement built in during the optimization by alternating the spins and bosons in the 1D representation, see Fig. 4. In this representation, each spin neighbors boson and , and hence entanglement between this spin and these bosons is directly generated as part of a two-site variational optimization procedure. The fluctuations of these boson modes can then couple to the fluctuations of further away spins, improving convergence to the global optimum. A further safeguard against variational metastability, in which a state with long-range coherence between spins and bosons far-separated in the MPS representation is explicitly mixed in during optimization to speed convergence, is described in Appendix A. Our method of choice for long-range time evolution is the second-order local Runge-Kutta method of Ref. Zaletel et al. (2015).

### ii.5 Observables and fidelities

While the computation of observables with MPSs, e.g. single-spin expectations or two-point spin-spin correlation functions, is standard and discussed in review articles Schollwöck (2011); Orús (2014), here we point out the measurement of two quantities which are non-standard. The first is the measurement of the full probability distribution of outcomes for a collective spin measurement along direction , i.e. the full counting statistics (FCS) for spins to be aligned along . To compute the FCS directly involves the computation of all correlation functions of any order, which in general scales poorly with the number of spins Dylewsky et al. (2016). However, the FCS can be computed efficiently given an MPS representation for the state of the spins by measuring the Fourier transform of this distribution, known as the characteristic function

(11) |

, and then inverse Fourier transforming. Noting that

(12) |

the characteristic function for each is the expectation of an MPO with bond dimension 1, which can be performed at the same cost as taking the overlap of two MPSs, see Fig. 5. The FCS is useful in trapped ion quantum simulators because it is readily accessible thanks to near single-ion resolution, and it provides detailed information about the structure of the state beyond low-order correlation functions Bohnet et al. (2015).

Often, rather than simulating a complex model of spins coupled to bosons, one would like to integrate out the bosons and keep only their virtual effect in the form of an effective model acting only on the spins. The most stringent test of the faithfulness of this spin-only description is the fidelity between the reduced density matrix of the spin-boson system obtained by tracing out the bosons with that of a pure spin state evolving under an effective spin Hamiltonian. This fidelity is defined as

(13) |

where is the density matrix obtained from tracing the bosons out from the pure spin-boson state and is the density matrix of pure state containing only spins. As is pure, we find that

(14) |

As we intend to trace out the bosons, it is useful to separate the bosons and spins spatially in the MPS representation. Such a representation naturally appears in our swapping dynamics protocol, see Fig. 3, and can be obtained straightforwardly from any representation using swap gates. With this ordering, the fidelity has a simple tensor network representation as the norm of a vector , provided both and have known MPS representations and the boson tensors are in left-canonical form Schollwöck (2011), as shown in Fig. 6. The restriction on canonical form, which involves no loss of generality, simply uses the freedom inherent in the MPS ansatz to transform the contraction over boson states into the identity, making the trace over bosonic degrees of freedom particularly simple.

## Iii Benchmark calculations for 1D and 2D trapped ion systems

As an application of the above methods, we will consider crystals of trapped ions subject to a spin-dependent force. Such systems have been a topic of great recent interest due to their capability of performing high-fidelity quantum simulation of long-range interacting quantum spin models Porras and Cirac (2004). In these systems, a competition between the Coulomb repulsion of the ions and external trapping from electromagnetic potentials leads to a crystalline structure for the ions in equilibrium. The normal modes of these equilibrium crystal structures comprise a set of global boson (phonon) modes, which are then coupled to the spin of the ions using lasers in a Raman scheme with characteristic wavevector to impart a spin-dependent force along the direction perpendicular to the crystal structure Leibfried et al. (2003); Milburn et al. (2000); Sørensen and Mølmer (2000, 1999); Islam et al. (2013). After performing a frame transformation on the spins, the Hamiltonian is expanded to lowest order in the Lamb-Dicke parameters , where is the mass of an ion and are the phonon normal mode frequencies. In addition, we will transform the phonons to a frame rotating with , the Raman beatnote frequency of the beams creating the spin-dependent force, and enact a rotating wave approximation. While this approximation will fail for very large detunings, we have verified that it is an excellent approximation for the parameter regimes explored in this work, and it leads to time-independent Hamiltonians which are more amenable to numerical computation. Following these steps, the Hamiltonian takes the form of Eq. (1) with , , and

(15) |

Here, is the detuning of mode from the Raman beatnote frequency, , with the magnitude of the spin-dependent force, quantifies the strength of the spin-dependent force for mode , and is the normal mode amplitude vector of phonon .

The effective Hamiltonian describing the time-dependent dynamics of Eq. (15) can be found by moving to an interaction picture which rotates with and using the Magnus series Magnus (1954); Blanes et al. (2009). Due to the fact that all spin operators in commute and the boson operators form copies of the Heisenberg algebra, the Magnus series truncates exactly at second order, and the propagator may be written exactly within the rotating wave approximation as

(16) |

with the spin-phonon coupling propagator

(17) |

and the spin-spin coupling propagator

(18) |

In these expressions, we have defined

(19) |

and

(20) |

Noting that the last term in the braces of Eq. (20) is bounded while the first grows without bound, is commonly approximated as , where

(21) |

In this approximation, is the propagator of a long-range Ising model .

Trapped ion quantum spin simulator experiments can be divided into two regimes parameterized by their hierarchy of energy scales. This categorization also currently coincides with the effective dimensionality of the ion crystal in present experiments, with 1D crystals corresponding to linear Paul traps Kim et al. (2009), and 2D systems forming in Penning traps Britton et al. (2012). Experiments in 1D Paul traps typically operate in the regime , in which all modes of the ion crystal participate in the dynamics. In contrast, current experiments in the 2D Penning trap operate in the regime where for all modes except for the center of mass mode , for which . Hence, the center of mass mode dominates the dynamics in these experiments. We stress that both types of traps can operate in both parameters regimes in principle, and our grouping of parameter regimes with dimensionality refers only to current experiments.

We numerically find the equilibrium positions of the ions following Ref. James (1998) for the Paul trap and Ref. Wang et al. (2013) for the Penning trap; typical equilibrium structures and effective spin-spin couplings are shown in Fig. 7. The spin-spin interactions in the 2D array with detuning close to the center of mass, shown in Fig. 7(d), are well-approximated by a uniform, all-to-all interaction. In contrast, the effective range of the interactions for typical parameters in linear Paul traps, modeled as a power law decay Britton et al. (2012); Khan et al. (2015), depends on the detuning, see Fig. 7(b)-(c).

Remarkably, an exact solution for the dynamics of Eq. (15) exists when starting from an uncorrelated product state of phonons and spins Dylewsky et al. (2016). We will use this exact solution to numerically benchmark our methodologies in both the 2D case, where we keep only a single boson mode and use the swapping method described in Sec. II.2, and in the 1D case, where we keep all modes and use the methods of Sec. II.4. In particular, we will compare the one- and two-point correlation functions

(22) |

(23) |

(24) |

where and we have taken the initial state in which all phonons start in the vacuum state and all spins initial point along the direction. In general, the phonon temperature in trapped ion experiments is nonzero. Our methods as presented here, as well as the exact solution, can be readily extended to any phonon pure state, and hence be used to generate thermal expectations by summing over Fock state expectations. The effects of finite phonon temperature have been analyzed using the exact solution in Ref. Dylewsky et al. (2016). A more efficient method for large numbers of phonons may be to use MPS methods specifically designed for finite-temperature systems Feiguin and White (2005); White (2009).

In what follows, we also use an MPS ansatz which explicitly conserves the symmetry of Eq. (15) associated with the parity operators and for spin-boson and spin systems, respectively, as conservation of symmetries leads to significant computational gains in MPS algorithms Singh et al. (2010, 2011). For the MPS calculations on pure spin models, we calculate the dynamics by fitting a long-range spin model with time-independent spin-spin couplings determined by Eq. (21) to an MPO, and then time-evolve using the local Runge-Kutta method of Ref. Zaletel et al. (2015). As these spin-spin couplings are not translationally invariant, we use the non-uniform exponential MPO fitting method developed in Refs. Koller et al. (2016); Wall (2016).

We begin our benchmarking by looking at the 1D case, in which the drive of the spin-dependent force is far detuned above the center of mass mode and so all phonon modes contribute to the dynamics. In particular, we consider the detunings above the center of mass mode kHz and kHz, employed in Ref. Richerme et al. (2014), and which results in the effective spin-spin couplings shown in Fig. 7(c). When the spin-dependent force is turned on, the effective Ising interactions will cause a coherent depolarization of the collective spin length concurrent with the development of collective spin-spin correlations . In addition, spin-phonon entanglement resulting from spin-dependent displacements of the phonon modes will cause oscillations in collective spin observables. Both effects are seen in Fig. 8, which compares the numerically computed collective spin expectations (red symbols) with the exact, analytic expressions Eq. (22)-(24) (blue lines). Panel (a) shows the depolarization of the collective spin for the 80kHz detuning, with the inset showing the small amplitude oscillations due to dynamical phonons which are perfectly captured within our numerical approach. Panel (b) shows the development of collective spin-spin correlations for the 150kHz detuned situation. Here, the numerical points are indistinguishable from the analytical curves on the scale of the plot, again demonstrating excellent agreement.

As mentioned above, the most stringent test of the fidelity of quantum simulation of a pure spin model can be obtained by measuring the fidelity of the spin reduced density matrix obtained from the wavefunction of the full spin-boson dynamics by tracing out the bosons, , with the pure spin density matrix obtained from evolution under the time-independent Ising model, Eq. (21). We compute this quantity using the methods of Sec. II.5; the results are shown in Fig. 9 for the two detunings we consider. We see rapid oscillations at many frequencies, corresponding to the generation of spin-phonon entanglement with all phonon modes. However, throughout the course of reasonable experimental timescales, the fidelity remains quite high, . We stress that the computation of this fidelity in the multimode case using the exact solution of the spin-boson dynamics Dylewsky et al. (2016) scales poorly with the number of spins; in contrast, this quantity is easily calculated within our numerical framework

We now turn to the 2D case, in which only the center of mass mode participates in the dynamics and we use the swapping method (Sec. II.2) to compute the dynamics. We use the parameters of Fig. 7, which are similar to recent experiments Bohnet et al. (2015). A key observable in these experiments is the Ramsey squeezing parameter, defined as Kitagawa and Ueda (1993)

(25) |

with and . The squeezing parameter is of interest because the condition witnesses entanglement in the system Sørensen and Mølmer (2001), and also quantifies the signal to noise enhancement of the spin state for performing Ramsey spectroscopy compared to an uncorrelated spin state Wineland et al. (1992). Fig. 10 displays the squeezing as a function of time starting from all spins aligned along , and compares this with the results of simulating the time-independent Ising model Eq. (21). While the pure spin approach predicts a smooth buildup of squeezing to a maximal point followed by monotonic destruction of squeezing, the spin-boson dynamics displays oscillations resulting from the development of spin-phonon entanglement, which tends to antisqueeze the spin distribution, as well as from the time dependence of the spin-spin couplings ( see Eq. (20)). However, the two results agree at the decoupling times , an integer, where the spin-phonon entanglement vanishes. While we do not show it here, we again find excellent agreement between our computed values and the exact solution, validating the swapping approach to dynamics.

The modification of spin observables due to spin-phonon entanglement is strikingly seen by comparing the full counting statistics (FCS) of collective spin observables between the spin-boson and pure spin approaches. The computation of these quantities from the analytical solution scale polynomially in the number of spins in the case that only the center of mass mode participates in the dynamics Bohnet et al. (2015), but also requires high-precision arithmetic to avoid catastrophic cancellation. In the general case, the computation of these quantities scales exponentially in . However, as discussed in Sec. II.5, the computation of the FCS is particularly straightforward within our MPS-based approach. Fig. 11 displays the FCS along the antisqueezed quadrature as well as the squeezed quadrature such that minimizes the spin variance . The upper panels (a)-(b) display the squeezed quadrature, which show the development of a narrow feature which is useful for sub-shot-noise spectroscopy, with panel (a) corresponding to the coupled spin-boson dynamics and panel (b) to the pure spin dynamics. The pure spin dynamics shows a monotonic development of this narrow feature, while this feature is periodically broadened and its height modulated in the spin-boson dynamics due to spin-phonon entanglement. The antisqueezed quadratures, shown in panels (c)-(d), initially broaden with respect to the initial uncorrelated state before developing fine-scale features whose usefulness for metrology are not witnessed by spin squeezing but can be captured via other measures, such as the quantum Fisher information Pezzé and Smerzi (2009); Bohnet et al. (2015). A comparison between the pure spin and spin-boson calculations again shows periodic differences driven by spin-phonon entanglement, but the differences are more slight compared to the squeezed quadrature.

## Iv Including decoherence: Quantum Trajectories

An important component of many AMO experiments is losses and other forms of dissipation, which are usually thought to degrade coherence and destroy quantum correlations, but can also lead to the production of entangled steady states Krauter et al. (2011); Foss-Feig et al. (2012); Lin et al. (2013). A major source of decoherence in trapped ion systems arises from scattering of the light used to create the spin-dependent force Uys et al. (2010). A full description of the system including losses is found by solving the master equation

(26) |

where is the system density matrix and the action of the Lindbladian superoperator accounting for spontaneous emission from the Raman beams is given as

(27) |

Here, the terms proportional to and correspond to spontaneous excitation and de-excitation, respectively, from Raman scattering, and the term proportional to represents elastic dephasing.

It is possible to apply MPS-based approaches to directly simulate master equations by replacing states by operators (“superkets”) and operators by superoperators, in which the local dimension of superket is , the square of the physical on-site Hilbert space Zwolak and Vidal (2004); Verstraete et al. (2004). Here, we instead use the method of quantum trajectories, in its first-order incarnation Plenio and Knight (1998); Daley (2014). To do so, it is useful to re-write the master equation as

(28) |

where the non-Hermitian effective Hamiltonian is

(29) |

and the recycling term is

(30) |

In these expressions, the “jump operators” are a basis in which the Lindbladian may be written as

(31) |

In the present case, we take the basis of jump operators to be the operators for . Now, at each time step, we evolve under to obtain . Because the effective Hamiltonian is non-Hermitian, the norm of this state, call it , will be less than 1. At first order in , is comprised of a sum of terms , where

(32) |

We now stochastically choose whether to use this state at the next time step, or undergo a “quantum jump” by applying the recycling term. Namely, we choose a random number , and select the state at the next time step to be

and then renormalizing the state. The equally weighted average of the expectation of an operator taken with such stochastically generated trajectories converges to the dynamics of this same observable generated by the master equation with an error that is asymptotically statistical, . In closing, we note that the quantum trajectories method may be trivially parallelized.

## V Application: Trapped ion spin-boson simulators in the presence of molecular impurities

For a pure spin Ising model (i.e., without dynamical phonons), the summation over all quantum trajectories can be performed analytically, leading to an exact solution of arbitrary Ising models in the presence of decoherence Foss-Feig et al. (2013a, b). When the dynamical motion of the phonons is included, however, no closed-form summation over quantum trajectories exists, and the dynamics must be obtained numerically. In what follows, we will consider the dynamics in the presence of decoherence as parameterized above in the 2D Penning trap where only the center of mass mode is relevant. In principle, these dynamics can be calculated efficiently, even in the presence of decoherence, for a spatially uniform initial state by restricting the Hilbert space to the space of permutationally symmetric density matrices Carmichael (1999); Hartmann (2012); Richter et al. (2015). In order to exemplify the power and generality of our approach compared to these other methods, we thus consider a case in which the center of mass mode is no longer spatially uniform. In particular, we simulate the case of an ion crystal in which some atomic ions have converted to heavier mass molecular ions through collisions with background gas, focusing on a crystal of Be ions with BeH impurities as a particular example Sawyer et al. (2015).

Impurity molecular ions have a heavier mass than the atomic ions, and so move to the edge of the rotating crystal due to centrifugal forces. We compute the equilibrium positions and axial normal modes of the ion crystal in the presence of these impurities using the methods of Ref. McAneny et al. (2013); results for a crystal of 37 Be ions and 44 BeH impurities are shown in Fig. 12. Panel (a) shows the equilibrium crystal structure for trapping parameters similar to current experiments Bohnet et al. (2015), with filled (empty) circles denoting molecular impurities (atomic ions). The highest frequency axial mode, analogous to the center of mass mode in a crystal with no impurities, is shown in in panel (b). The impurity ions have a greatly reduced amplitude of oscillation due to their heavier mass, and, more importantly for our purposes, the amplitudes of the center of mass mode on the atomic ions are non-uniform. The spread in amplitudes is roughly 10% of the mean for the particular parameters we consider. Only the atomic ions respond to the light-induced spin-dependent force that generate effective Ising spin-spin interactions, and so the effective number of spins in the quantum spin simulator is the number of atomic ions: 37 in our example case. Still, the molecular ions influence the spin-phonon and spin-spin dynamics through their modification of the center of mass mode amplitudes on the atomic ions.

Our calculations of the spin dynamics for a spatially non-uniform center of mass mode with decoherence are shown in Fig. 13. We take decoherence rates of s, s, and s, comparable to current experiments, and average over trajectories. Panel (a) shows collective spin correlations such as the normalized total magnetization and its square . The black dashed line is the mean-field prediction of the decay of magnetization, which occurs solely due to decoherence at a rate . Strong deviations from this curve indicate the decay of magnetization due to the build-up of higher-order spin correlations. The blue dashed lines for the spin correlations use the actual inhomogeneous analog of the center of mass mode shown in Fig. 12(b). In contrast, the solid red lines use its nearest uniform approximation, in which the amplitude of the mode at each position is replaced by the average of the amplitude vector. Almost no difference between the exact and approximate predictions can be seen in the decay of the magnetization, whereas a slight difference can be seen in the squared collective magnetization at late times. In contrast, the effect of the mode inhomogeneity is clearly seen in the spatially resolved magnetization shown in Fig. 13(b). Here, the size of the circles are proportional to and the red solid (blue dashed) circles are for the uniform approximation and the true inhomogeneous mode, respectively. Our results demonstrate that, even in the presence of a relatively large number of impurity ions, the inhomogeneities in the analog of the center of mass mode have only a slight impact on low-order collective spin correlations over experimentally relevant timescales, whereas inhomogeneities in the spatial distribution of the magnetization can be quite sizable.

## Vi Conclusions and outlook

We have developed a generic approach to the out-of equilibrium dynamics of spins globally coupled to bosonic modes, which encompasses many paradigmatic models from all areas of quantum science, based on matrix product states (MPSs). In contrast to existing methods, our framework applies for any spatial and operator dependence of the spin-boson coupling, and places no restrictions on relative energy scales. In the regime of many fewer bosons than spins, , an efficient approach is to use a Trotter-Suzuki decomposition of the spin-boson coupling propagator and dynamically deform the MPS topology using swap gates, resulting in a method with rigorously controlled error. Alternatively, for equal numbers of spins and bosons, , we propose to perform direct time evolution using an exact matrix product operator (MPO) representation of the spin-boson Hamiltonian with long-ranged couplings, and provide techniques to safeguard against variational metastability. This exact MPO representation of the Hamiltonian can also be employed to find eigenstates and static properties of spin-boson systems.

Using our newly developed approach, we calculated the dynamics of 1D and 2D trapped ion quantum simulators in experimentally relevant parameter regimes, fully accounting for the dynamics of the phonons that mediate effective spin dynamics between ions. We first benchmarked our numerical results against a recent exact analytical solution for a particular spin-phonon coupling, finding excellent agreement. Beyond comparing low-order correlation functions which are straightforward and efficient to compute with the exact solution, we also demonstrated the power of our approach for computing detailed properties of the quantum state which are difficult to obtain through other means, such as the full counting statistics of collective spin observables and the fidelity of the spin-phonon system as a pure spin quantum simulator.

We showed that Markovian decoherence, which is an important component of trapped ions and many other AMO platforms, can be incorporated within our approach by solving the associated master equation using the method of quantum trajectories. Within this framework, we simulated the dynamics of a 2D ion crystal in the presence of realistic decoherence. In order to evince the power of our approach beyond those that require spatial uniformity of the spin-boson coupling, we treated the case of a non-uniform spin-boson coupling, physically motivated by the presence of impurity molecular ions in an atomic ion quantum spin simulator.

While our results show that MPSs are a powerful framework for studying spin-boson models, there are certain regimes which remain difficult or impossible to simulate efficiently with MPSs. In particular, 2D systems with many boson modes and finite-range interactions are challenging, and it is known that the scaling of MPSs for finding equilibrium states in 2D is exponential in the system size. Even in one dimension, simulations in which spins are strongly entangled with many bosons in a spatially inhomogeneous fashion are difficult, and particularly subject to variational metastability. An additional source of difficulty from a numerical point of view is that there is often a large separation of timescales between the dynamics of direct spin-boson coupling and boson-mediated spin-spin interactions, with the former limiting the size of the time step. Here, applying multiscale methods Kevorkian and Cole (May 15, 1996) is a promising avenue for improving numerical performance.

Myriad future research directions are available based on the framework developed in this paper. For one, our approach is able to treat several competing interaction scales, as well as decoherence, on the same footing, and so provide a rigorous means for justifying effective models of driven spin-boson systems in various parameter regimes. This is especially important for trapped ion quantum spin simulators in effective transverse fields, where the presence of non-commuting interactions with comparable energy scales can fundamentally change the effective spin physics in ways not captured by a perturbative approach Wall et al. (2016). Additionally, while we focused on applications to trapped ion systems, our techniques can be applied to many other platforms, including quantum optics and optomechanics, and provide a scalable means for studying strongly correlated physics in these diverse systems. Our framework also applies to higher-dimensional spin representations, such as spin-1, which can also be realized with trapped ion systems Cohen and Retzker (2014); Senko et al. (2015) and can shed light on fundamental questions such as the nature of topological phases with long-range interactions Gong et al. (2016). A generalization from discrete boson spectra to continuous boson spectra is possible, given a single global coupling operator between spins and bosons, using techniques based on orthogonal polynomials Prior et al. (2010); Chin et al. (2010) while keeping well-defined error Woods et al. (2015). This can lead to new insights into quantum phases of correlated spins which are coupled to an external bath.

## Vii Acknowledgements

We would like to acknowledge useful discussions with Justin Bohnet, John Bollinger, Johannes Schachenmayer, Zhexuan Gong, and Phil Richerme, and support from NSF-PHY 1521080, JILA-NSF-PFC-1125844, ARO, MURI-AFOSR and AFOSR. MLW thanks the NRC postdoctoral program for support.

## Appendix A Improving the convergence of long-range time evolution

In this appendix, we discuss additional safeguards against variational metastability when using the long-range time evolution methods of Sec. II.4. Time evolution proceeds by optimization of the functional , with an MPO representation of the propagator over a small time step Zaletel et al. (2015) and an MPS representation of the state at time , over all MPSs with fixed resources. The optimal MPS then becomes the new state at time . The global optimization over the entire MPS cannot be performed efficiently, and so optimization is performed over the tensors of two neighboring sites in the MPS representation (see Eq. (4)) with all other tensors held fixed Schollwöck (2011). A complete optimization cycle over all neighboring pairs of sites is referred to as a sweep. In this local optimization scheme, the main cause of metastability is that entanglement is only built up between neighboring sites in a single step, and longer-ranged entanglement must be built up by repeated sweeping.

In some cases, the algorithm is not able to properly build up the appropriate entanglement structure for a long-ranged Hamiltonian, especially when starting from an unentangled state of spins and bosons. Our idea to improve the convergence is to also endow our initial variational guess with long-range coherence between spins and bosons which are far-separated in the MPS representation. We do so by generating a state which contains long-range coherence between all spins and all bosons and mixing in a small amount of this state to the optimal state obtained from local variational minimization, with the mixing parameter taken to zero towards the end of the sweeping optimization. More precisely, the mixing is done by projecting the state onto the orthonormal basis formed by holding all tensors except two in the MPS representation fixed, and then adding this projected two-site wavefunction directly to the wavefunction obtained from the optimization procedure.

This approach is similar in spirit to the modified single-site DMRG algorithm by White White (2005), which also builds fluctuations on top of the local variational optimum to accelerate convergence. While a detailed study of the optimal states and mixing parameters is outside of the scope of this work, we find that choosing to be the ground state of and to perform well. Note that the ground state can be efficiently found using ordinary two-site DMRG with the MPO representation of Fig. 4, as entanglement which is directly built up between neighboring spins and bosons propagates to farther separated spins and bosons by repeated sweeping. We also note that an accurate representation of the ground state is not essential–only non-vanishing long-range coherence is required–and so relatively coarse tolerances can be used in obtaining the ground state.

## References

- Leroux et al. (2010) Ian D. Leroux, Monika H. Schleier-Smith, and Vladan Vuletić, “Implementation of cavity squeezing of a collective atomic spin,” Phys. Rev. Lett. 104, 073602 (2010).
- Bohnet et al. (2014) Justin G Bohnet, Kevin C Cox, Matthew A Norcia, Joshua M Weiner, Zilong Chen, and James K Thompson, “Reduced spin measurement back-action for a phase sensitivity ten times beyond the standard quantum limit,” Nature Photonics 8, 731–736 (2014).
- Hosten et al. (2016) Onur Hosten, Nils J. Engelsen, Rajiv Krishnakumar, and Mark A. Kasevich, “Measurement noise 100 times lower than the quantum-projection limit using entangled atoms,” Nature 529 (2016).
- Caldeira and Leggett (1983) AO Caldeira and A J Leggett, “Quantum tunnelling in a dissipative system,” Annals of Physics 149, 374–456 (1983).
- Sørensen and Mølmer (1999) Anders Sørensen and Klaus Mølmer, “Quantum computation with ions in thermal motion,” Phys. Rev. Lett. 82, 1971–1974 (1999).
- Porras and Cirac (2004) D. Porras and J. I. Cirac, “Effective quantum spin systems with trapped ions,” Phys. Rev. Lett. 92, 207901 (2004).
- Kim et al. (2009) K. Kim, M.-S. Chang, R. Islam, S. Korenblit, L.-M. Duan, and C. Monroe, “Entanglement and tunable spin-spin couplings between trapped ions using multiple transverse modes,” Phys. Rev. Lett. 103, 120502 (2009).
- Britton et al. (2012) Joseph W Britton, Brian C Sawyer, Adam C Keith, C-C Joseph Wang, James K Freericks, Hermann Uys, Michael J Biercuk, and John J Bollinger, “Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of spins,” Nature 484, 489–492 (2012).
- Porras et al. (2008) D. Porras, F. Marquardt, J. von Delft, and J. I. Cirac, “Mesoscopic spin-boson models of trapped ions,” Phys. Rev. A 78, 010101 (2008).
- Ritsch et al. (2013) Helmut Ritsch, Peter Domokos, Ferdinand Brennecke, and Tilman Esslinger, “Cold atoms in cavity-generated dynamical optical potentials,” Rev. Mod. Phys. 85, 553–601 (2013).
- Niemczyk et al. (2010) Thomas Niemczyk, F Deppe, H Huebl, EP Menzel, F Hocke, MJ Schwarz, JJ Garcia-Ripoll, D Zueco, T Hümmer, E Solano, et al., “Circuit quantum electrodynamics in the ultrastrong-coupling regime,” Nature Physics 6, 772–776 (2010).
- Hur et al. (2015) Karyn Le Hur, Loïc Henriet, Alexandru Petrescu, Kirill Plekhanov, Guillaume Roux, and Marco Schiró, “Many-body quantum electrodynamics networks: non-equilibrium condensed matter physics with light,” arXiv preprint arXiv:1505.00167 (2015).
- Aspelmeyer et al. (2014) Markus Aspelmeyer, Tobias J. Kippenberg, and Florian Marquardt, “Cavity optomechanics,” Rev. Mod. Phys. 86, 1391–1452 (2014).
- Carmichael (1999) Howard J. Carmichael, Statistical Methods in Quantum Optics 1 (Springer, Berlin, 1999).
- Hartmann (2012) Stephan Hartmann, “Generalized Dicke states,” arXiv preprint arXiv:1201.1732 (2012).
- Richter et al. (2015) Marten Richter, Michael Gegg, T. Sverre Theuerholz, and Andreas Knorr, “Numerically exact solution of the many emitter-cavity laser problem: Application to the fully quantized spaser emission,” Phys. Rev. B 91, 035306 (2015).
- Dylewsky et al. (2016) D. Dylewsky, J. K. Freericks, M. L. Wall, A. M. Rey, and M. Foss-Feig, “Nonperturbative calculation of phonon effects on spin squeezing,” Phys. Rev. A 93, 013415 (2016).
- Shi et al. (2006) Y.-Y. Shi, L.-M. Duan, and G. Vidal, “Classical simulation of quantum many-body systems with a tree tensor network,” Phys. Rev. A 74, 022320 (2006).
- Zaletel et al. (2015) Michael P. Zaletel, Roger S. K. Mong, Christoph Karrasch, Joel E. Moore, and Frank Pollmann, “Time-evolving a matrix product state with long-ranged interactions,” Phys. Rev. B 91, 165112 (2015).
- Richerme et al. (2014) Philip Richerme, Zhe-Xuan Gong, Aaron Lee, Crystal Senko, Jacob Smith, Michael Foss-Feig, Spyridon Michalakis, Alexey V Gorshkov, and Christopher Monroe, “Non-local propagation of correlations in quantum systems with long-range interactions,” Nature 511, 198–201 (2014).
- Bohnet et al. (2015) Justin G Bohnet, Brian C Sawyer, Joseph W Britton, Michael L Wall, Ana Maria Rey, Michael Foss-Feig, and John J Bollinger, “Quantum spin dynamics and entanglement generation with hundreds of trapped ions,” arXiv preprint arXiv:1512.03756 (2015).
- Gardas (2011) Bartłomiej Gardas, “Exact solution of the schrödinger equation with the spin-boson hamiltonian,” Journal of Physics A: Mathematical and Theoretical 44, 195301 (2011).
- Braak (2011) D. Braak, “Integrability of the rabi model,” Phys. Rev. Lett. 107, 100401 (2011).
- Jaynes and Cummings (1963) Edwin T Jaynes and Frederick W Cummings, “Comparison of quantum and semiclassical radiation theories with application to the beam maser,” Proceedings of the IEEE 51, 89–109 (1963).
- Eberly et al. (1980) J. H. Eberly, N. B. Narozhny, and J. J. Sanchez-Mondragon, “Periodic spontaneous collapse and revival in a simple quantum model,” Phys. Rev. Lett. 44, 1323–1326 (1980).
- Agarwal et al. (2012) S. Agarwal, S. M. Hashemi Rafsanjani, and J. H. Eberly, ‘‘Tavis-cummings model beyond the rotating wave approximation: Quasidegenerate qubits,” Phys. Rev. A 85, 043815 (2012).
- Jeckelmann and White (1998) Eric Jeckelmann and Steven R. White, “Density-matrix renormalization-group study of the polaron problem in the holstein model,” Phys. Rev. B 57, 6376–6385 (1998).
- Guo et al. (2012) Cheng Guo, Andreas Weichselbaum, Jan von Delft, and Matthias Vojta, “Critical and strong-coupling phases in one- and two-bath spin-boson models,” Phys. Rev. Lett. 108, 160401 (2012).
- Frenzel and Plenio (2013) Max F Frenzel and Martin B Plenio, “Matrix product state representation without explicit local hilbert space truncation with applications to the sub-ohmic spin-boson model,” New Journal of Physics 15, 073046 (2013).
- Bursill et al. (1998) Robert J. Bursill, Ross H. McKenzie, and Chris J. Hamer, “Phase diagram of the one-dimensional holstein model of spinless fermions,” Phys. Rev. Lett. 80, 5607–5610 (1998).
- Brockt et al. (2015) C. Brockt, F. Dorfner, L. Vidmar, F. Heidrich-Meisner, and E. Jeckelmann, “Matrix-product-state method with a dynamical local basis optimization for bosonic systems out of equilibrium,” Phys. Rev. B 92, 241106 (2015).
- Schröder et al. (2015) Florian A. Y. N. Schröder, Alex W Chin, and Richard H Friend, “Simulating open quantum dynamics with time-dependent variational matrix product states: Towards microscopic correlation of environment dynamics and reduced system evolution,” arXiv preprint arXiv:1507.02202 (2015).
- Gammelmark and Mølmer (2012) Søren Gammelmark and Klaus Mølmer, “Interacting spins in a cavity: Finite-size effects and symmetry-breaking dynamics,” Phys. Rev. A 85, 042114 (2012).
- Schollwöck (2011) Ulrich Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Annals of Physics 326, 96–192 (2011).
- Orús (2014) Román Orús, ‘‘A practical introduction to tensor networks: Matrix product states and projected entangled pair states,” Annals of Physics 349, 117–158 (2014).
- Hastings (2007) Matthew B Hastings, “An area law for one-dimensional quantum systems,” Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007).
- Crosswhite et al. (2008) Gregory M. Crosswhite, A. C. Doherty, and Guifré Vidal, “Applying matrix product operators to model systems with long-range interactions,” Phys. Rev. B 78, 035116 (2008).
- Schachenmayer et al. (2013) J. Schachenmayer, B. P. Lanyon, C. F. Roos, and A. J. Daley, “Entanglement growth in quench dynamics with variable range interactions,” Phys. Rev. X 3, 031015 (2013).
- Wall and Carr (2012) M L Wall and Lincoln D Carr, “Out-of-equilibrium dynamics with matrix product states,” New Journal of Physics 14, 125015 (2012).
- Pollmann et al. (2009) Frank Pollmann, Subroto Mukerjee, Ari M. Turner, and Joel E. Moore, “Theory of finite-entanglement scaling at one-dimensional quantum critical points,” Phys. Rev. Lett. 102, 255701 (2009).
- Rams et al. (2015) Marek M. Rams, Valentin Zauner, Matthias Bal, Jutho Haegeman, and Frank Verstraete, “Truncating an exact matrix product state for the xy model: Transfer matrix and its renormalization,” Phys. Rev. B 92, 235150 (2015).
- Mishmash et al. (2011) Ryan V. Mishmash, Matthew S. Block, Ribhu K. Kaul, D. N. Sheng, Olexei I. Motrunich, and Matthew P. A. Fisher, “Bose metals and insulators on multileg ladders with ring exchange,” Phys. Rev. B 84, 245127 (2011).
- Stoudenmire and White (2012) E.M. Stoudenmire and Steven R. White, ‘‘Studying two-dimensional systems with the density matrix renormalization group,” Annual Review of Condensed Matter Physics 3, 111–128 (2012), http://dx.doi.org/10.1146/annurev-conmatphys-020911-125018 .
- Depenbrock et al. (2012) Stefan Depenbrock, Ian P. McCulloch, and Ulrich Schollwöck, “Nature of the spin-liquid ground state of the heisenberg model on the kagome lattice,” Phys. Rev. Lett. 109, 067201 (2012).
- McCulloch (2007) Ian P McCulloch, “From density-matrix renormalization group to matrix product states,” Journal of Statistical Mechanics: Theory and Experiment 2007, P10014 (2007).
- Verstraete et al. (2004) F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, “Matrix product density operators: Simulation of finite-temperature and dissipative systems,” Phys. Rev. Lett. 93, 207204 (2004).
- Crosswhite and Bacon (2008) Gregory M. Crosswhite and Dave Bacon, “Finite automata for caching in matrix product algorithms,” Phys. Rev. A 78, 012356 (2008).
- Fröwis et al. (2010) F. Fröwis, V. Nebendahl, and W. Dür, “Tensor operators: Constructions and applications for long-range interaction systems,” Phys. Rev. A 81, 062337 (2010).
- Pirvu et al. (2010) Bogdan Pirvu, Valentin Murg, J Ignacio Cirac, and Frank Verstraete, “Matrix product operator representations,” New Journal of Physics 12, 025012 (2010).
- Huyghebaert and Raedt (1990) J Huyghebaert and H De Raedt, “Product formula methods for time-dependent Schrodinger problems,” Journal of Physics A: Mathematical and General 23, 5777 (1990).
- Omelyan et al. (2002) IP Omelyan, IM Mryglod, and Reinhard Folk, “Optimized forest–ruth-and suzuki-like algorithms for integration of motion in many-body systems,” Computer Physics Communications 146, 188–202 (2002).
- Vidal (2004) Guifré Vidal, “Efficient simulation of one-dimensional quantum many-body systems,” Phys. Rev. Lett. 93, 040502 (2004).
- Daley et al. (2004) Andrew John Daley, Corinna Kollath, Ulrich Schollwöck, and Guifré Vidal, “Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces,” Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
- Sornborger and Stewart (1999) A. T. Sornborger and E. D. Stewart, “Higher-order methods for simulations on quantum computers,” Phys. Rev. A 60, 1956–1965 (1999).
- Manmana et al. (2005) Salvatore R. Manmana, Alejandro Muramatsu, and Reinhard M. Noack, “Time evolution of oneâdimensional quantum many body systems,” AIP Conference Proceedings 789 (2005).
- García-Ripoll (2006) Juan José García-Ripoll, “Time evolution of matrix product states,” New Journal of Physics 8, 305 (2006).
- Haegeman et al. (2014) Jutho Haegeman, Christian Lubich, Ivan Oseledets, Bart Vandereycken, and Frank Verstraete, “Unifying time evolution and optimization with matrix product states,” arXiv preprint arXiv:1408.5056 (2014).
- Danshita and Naidon (2009) Ippei Danshita and Pascal Naidon, “Bose-hubbard ground state: Extended bogoliubov and variational methods compared with time-evolving block decimation,” Phys. Rev. A 79, 043601 (2009).
- Stoudenmire and White (2010) E M Stoudenmire and Steven R White, “Minimally entangled typical thermal state algorithms,” New Journal of Physics 12, 055026 (2010).
- Gopalakrishnan et al. (2009) Sarang Gopalakrishnan, Benjamin L Lev, and Paul M Goldbart, “Emergent crystallinity and frustration with bose–einstein condensates in multimode cavities,” Nature Physics 5, 845–850 (2009).
- Graß et al. (2015) Tobias Graß, David Raventós, Bruno Juliá-Díaz, Christian Gogolin, and Maciej Lewenstein, “Controlled complexity in trapped ions: from quantum mattis glasses to number partitioning,” arXiv preprint arXiv:1507.07863 (2015).
- Leibfried et al. (2003) Dietrich Leibfried, Brian DeMarco, Volker Meyer, David Lucas, Murray Barrett, Joe Britton, WM Itano, B Jelenković, Chris Langer, Till Rosenband, et al., “Experimental demonstration of a robust, high-fidelity geometric two ion-qubit phase gate,” Nature 422, 412–415 (2003).
- Milburn et al. (2000) GJ Milburn, S Schneider, and DFV James, “Ion trap quantum computing with warm ions,” Fortschritte der Physik 48, 801–810 (2000).
- Sørensen and Mølmer (2000) Anders Sørensen and Klaus Mølmer, “Entanglement and quantum computation with ions in thermal motion,” Phys. Rev. A 62, 022311 (2000).
- Islam et al. (2013) R Islam, C Senko, WC Campbell, S Korenblit, J Smith, A Lee, EE Edwards, C-CJ Wang, JK Freericks, and C Monroe, “Emergence and frustration of magnetism with variable-range interactions in a quantum simulator,” Science 340, 583–587 (2013).
- Magnus (1954) Wilhelm Magnus, “On the exponential solution of differential equations for a linear operator,” Communications on Pure and Applied Mathematics 7, 649–673 (1954).
- Blanes et al. (2009) S. Blanes, F. Casas, J.A. Oteo, and J. Ros, “The magnus expansion and some of its applications,” Physics Reports 470, 151 – 238 (2009).
- James (1998) Daniel FV James, ‘‘Quantum dynamics of cold trapped ions with application to quantum computation,” Applied Physics B: Lasers and Optics 66, 181–190 (1998).
- Wang et al. (2013) C.-C. Joseph Wang, Adam C. Keith, and J. K. Freericks, “Phonon-mediated quantum spin simulator employing a planar ionic crystal in a penning trap,” Phys. Rev. A 87, 013422 (2013).
- Khan et al. (2015) A. Khan, B. Yoshimura, and J. K. Freericks, “Theoretical basis for quantum simulation with a planar ionic crystal in a penning trap using a triangular rotating wall,” Phys. Rev. A 92, 043405 (2015).
- Feiguin and White (2005) Adrian E. Feiguin and Steven R. White, “Finite-temperature density matrix renormalization using an enlarged Hilbert space,” Phys. Rev. B 72, 220401 (2005).
- White (2009) Steven R. White, “Minimally entangled typical quantum states at finite temperature,” Phys. Rev. Lett. 102, 190601 (2009).
- Singh et al. (2010) Sukhwinder Singh, Robert N. C. Pfeifer, and Guifré Vidal, “Tensor network decompositions in the presence of a global symmetry,” Phys. Rev. A 82, 050301 (2010).
- Singh et al. (2011) Sukhwinder Singh, Robert N. C. Pfeifer, and Guifre Vidal, “Tensor network states and algorithms in the presence of a global u(1) symmetry,” Phys. Rev. B 83, 115125 (2011).
- Koller et al. (2016) Andrew P Koller, Michael L Wall, Josh Mundinger, and Ana Maria Rey, “Dynamics of interacting fermions in spin-dependent potentials,” arXiv preprint arXiv:1601.01004 (2016).
- Wall (2016) Michael L Wall, “In preparation,” (2016).
- Kitagawa and Ueda (1993) Masahiro Kitagawa and Masahito Ueda, “Squeezed spin states,” Phys. Rev. A 47, 5138–5143 (1993).
- Sørensen and Mølmer (2001) Anders S. Sørensen and Klaus Mølmer, “Entanglement and extreme spin squeezing,” Phys. Rev. Lett. 86, 4431–4434 (2001).
- Wineland et al. (1992) D. J. Wineland, J. J. Bollinger, W. M. Itano, F. L. Moore, and D. J. Heinzen, “Spin squeezing and reduced quantum noise in spectroscopy,” Phys. Rev. A 46, R6797–R6800 (1992).
- Pezzé and Smerzi (2009) Luca Pezzé and Augusto Smerzi, “Entanglement, nonlinear dynamics, and the heisenberg limit,” Phys. Rev. Lett. 102, 100401 (2009).
- Krauter et al. (2011) Hanna Krauter, Christine A. Muschik, Kasper Jensen, Wojciech Wasilewski, Jonas M. Petersen, J. Ignacio Cirac, and Eugene S. Polzik, “Entanglement generated by dissipation and steady state entanglement of two macroscopic objects,” Phys. Rev. Lett. 107, 080503 (2011).
- Foss-Feig et al. (2012) Michael Foss-Feig, Andrew J. Daley, James K. Thompson, and Ana Maria Rey, “Steady-state many-body entanglement of hot reactive fermions,” Phys. Rev. Lett. 109, 230501 (2012).
- Lin et al. (2013) Y Lin, JP Gaebler, F Reiter, TR Tan, R Bowler, AS Sørensen, D Leibfried, and DJ Wineland, “Dissipative production of a maximally entangled steady state of two quantum bits,” Nature 504, 415–418 (2013).
- Uys et al. (2010) H. Uys, M. J. Biercuk, A. P. VanDevender, C. Ospelkaus, D. Meiser, R. Ozeri, and J. J. Bollinger, “Decoherence due to elastic rayleigh scattering,” Phys. Rev. Lett. 105, 200401 (2010).
- Zwolak and Vidal (2004) Michael Zwolak and Guifré Vidal, “Mixed-state dynamics in one-dimensional quantum lattice systems: A time-dependent superoperator renormalization algorithm,” Phys. Rev. Lett. 93, 207205 (2004).
- Plenio and Knight (1998) M. B. Plenio and P. L. Knight, “The quantum-jump approach to dissipative dynamics in quantum optics,” Rev. Mod. Phys. 70, 101–144 (1998).
- Daley (2014) Andrew J. Daley, “Quantum trajectories and open many-body quantum systems,” Advances in Physics 63, 77–149 (2014), http://dx.doi.org/10.1080/00018732.2014.933502 .
- Foss-Feig et al. (2013a) Michael Foss-Feig, Kaden R. A. Hazzard, John J. Bollinger, and Ana Maria Rey, “Nonequilibrium dynamics of arbitrary-range ising models with decoherence: An exact analytic solution,” Phys. Rev. A 87, 042101 (2013a).
- Foss-Feig et al. (2013b) M Foss-Feig, K R A Hazzard, J J Bollinger, A M Rey, and C W Clark, “Dynamical quantum correlations of ising models on an arbitrary lattice and their resilience to decoherence,” New Journal of Physics 15, 113008 (2013b).
- Sawyer et al. (2015) Brian C. Sawyer, Justin G. Bohnet, Joseph W. Britton, and John J. Bollinger, “Reversing hydride-ion formation in quantum-information experiments with ,” Phys. Rev. A 91, 011401 (2015).
- McAneny et al. (2013) M. McAneny, B. Yoshimura, and J. K. Freericks, “Effect of defects on phonons and the effective spin-spin interactions of an ultracold penning-trap quantum simulator,” Phys. Rev. A 88, 043434 (2013).
- Kevorkian and Cole (May 15, 1996) J.K. Kevorkian and J.D. Cole, Multiple Scale and Singular Perturbation Methods, edited by 1 edition (Springer, May 15, 1996).
- Wall et al. (2016) Michael L Wall, Arghavan Safavi-Naini, and Ana Maria Rey, “In preparation,” (2016).
- Cohen and Retzker (2014) I. Cohen and A. Retzker, “Proposal for verification of the haldane phase using trapped ions,” Phys. Rev. Lett. 112, 040503 (2014).
- Senko et al. (2015) C. Senko, P. Richerme, J. Smith, A. Lee, I. Cohen, A. Retzker, and C. Monroe, “Realization of a quantum integer-spin chain with controllable interactions,” Phys. Rev. X 5, 021026 (2015).
- Gong et al. (2016) Z.-X. Gong, M. F. Maghrebi, A. Hu, M. L. Wall, M. Foss-Feig, and A. V. Gorshkov, “Topological phases with long-range interactions,” Phys. Rev. B 93, 041102 (2016).
- Prior et al. (2010) Javier Prior, Alex W. Chin, Susana F. Huelga, and Martin B. Plenio, “Efficient simulation of strong system-environment interactions,” Phys. Rev. Lett. 105, 050404 (2010).
- Chin et al. (2010) Alex W. Chin, Ángel Rivas, Susana F. Huelga, and Martin B. Plenio, “Exact mapping between system-reservoir quantum models and semi-infinite discrete chains using orthogonal polynomials,” Journal of Mathematical Physics 51, 092109 (2010), http://dx.doi.org/10.1063/1.3490188.
- Woods et al. (2015) M. P. Woods, M. Cramer, and M. B. Plenio, “Simulating bosonic baths with error bars,” Phys. Rev. Lett. 115, 130401 (2015).
- White (2005) Steven R. White, “Density matrix renormalization group algorithms with a single center site,” Phys. Rev. B 72, 180403 (2005).