Entanglement growth in quench dynamics with variable range interactions

Studying entanglement growth in quantum dynamics provides both insight into the underlying microscopic processes and information about the complexity of the quantum states, which is related to the efficiency of simulations on classical computers. Recently, experiments with trapped ions, polar molecules, and Rydberg excitations have provided new opportunities to observe dynamics with long-range interactions. We explore nonequilibrium coherent dynamics after a quantum quench in such systems, identifying qualitatively different behavior as the exponent of algebraically decaying spin-spin interactions in a transverse Ising chain is varied. Computing the build-up of bipartite entanglement as well as mutual information between distant spins, we identify linear growth of entanglement entropy corresponding to propagation of quasiparticles for shorter range interactions, with the maximum rate of growth occurring when the Hamiltonian parameters match those for the quantum phase transition. Counter-intuitively, the growth of bipartite entanglement for long-range interactions is only logarithmic for most regimes, i.e., substantially slower than for shorter range interactions. Experiments with trapped ions allow for the realization of this system with a tunable interaction range, and we show that the different phenomena are robust for finite system sizes and in the presence of noise. These results can act as a direct guide for the generation of large-scale entanglement in such experiments, towards a regime where the entanglement growth can render existing classical simulations inefficient.


I. INTRODUCTION
Advances with atomic molecular and optical (AMO) systems including cold atoms, entangled photons, and trapped ions has rapidly opened possibilities to explore many-body physics in a highly controllable way [1][2][3]. A key example of this is the new possibility to explore coherent non-equilibrium dynamics in a closed many-body system, e.g., the dynamics induced by quantum quenches [4][5][6][7][8][9][10][11]. There have been several recent quench experiments with cold atoms in optical lattices, which not only probe the microscopic behavior of the system, e.g., the propagation of quasi-particles [12], but also indicate the possibility to probe dynamics beyond the regimes that are currently accessible to classical simulations [13][14][15]. In this context, the growth of entanglement in the system underlies the complexity of simulating the dynamics classically. Demonstration of large entanglement growth after a quantum quench would be a crucial step in demonstrating the possibility to use these systems as controllable quantum simulators, effectively using experimental systems to compute dynamics in a way that exceeds the capabilities of classical computations [16].
Systems of trapped ions are a very promising candidate for realizing a quantum simulator, due to the control already demonstrated in the development of gate-based quantum computation and simulation [17][18][19] with these systems, and the ability to make measurements by state tomography [20]. Recently, analogue quantum simulation of interacting spin systems [21,22] was also realized in ion traps [23][24][25] with a key novel element being the possibility to realize variable-range interactions [26][27][28][29][30][31], as shown in Fig. 1, in contrast to the short-range interactions of neutral atoms, or the dipole-dipole interactions possible with polar molecules. So far, these variable-range interactions were discussed primarily in the case of ground-state calculations and near-adiabatic dynamics. Here we explore nonequilibrium coherent dynamics after a quantum quench in these systems, identifying qualitatively different behavior as the exponent α of algebraically decaying spinspin interactions is varied. Beginning with all spins aligned with a transverse field, we use a combination of analytical and numerical methods to compute the dynamics after the Ising interactions are quenched on, incorporating matrix product operator techniques [32][33][34][35][36][37][38][39] to treat variable long-range interactions with up to 50 spins.
In particular, we investigate the build-up of bipartite entanglement in the chain as well as mutual information between distant spins. For interactions with α 1 we show that the behavior is qualitatively similar to nearestneighbor interactions, with correlation build-up well described by the propagation of quasi-particles at a rate equal to or slower than the Lieb-Robinson bound [40][41][42]. This leads to a linear increase in bipartite entanglement in time, so that the dynamics cannot be efficiently computed in existing classical simulations beyond short times [14,15]. Interestingly, in this limit we find that the maximum growth rate of bipartite entanglement even in small systems occurs when we quench the interaction strength to the value corresponding to the quantum phase transition point, shifting accordingly for varying α. For interactions with α 1, we observe qualitatively different behavior. Counter-intuitively, quenches above the critical point for these long-range interactions lead only to logarithmic increase of bipartite entanglement in time, so that in this regime, long-range interactions produce a slower growth of entanglement than short-range interactions. For the extreme case of infinite range interactions we show that the bipartite entanglement is bounded by a constant value.
Finally we discuss specific experimental parameters for the realization of different regimes in ion traps with finite chain lengths, and experimental measurement protocols for these effects, opening possibilities for the regimes considered here to be observed in the laboratory. We take typical experimental noise sources into account and show that the observable features are robust against these. The result that long-range interactions do not always give rise to strong entanglement in quench dynamics has implications for the realization of large scale-entanglement in quantum simulations in general systems with longrange interactions. This paper is organized as follows. In Sec. II we introduce the setup and model, as well as the entanglement measures we compute. In Sec. III we show how the entanglement growth depends on the model parameters and how the entanglement distribution mechanisms can be understood. In Sec. IV we show entanglement growth for typical experimental parameters and how the effects can be measured in noisy experiments. The final Sec. V provides a conclusion and an outlook.

II. MODEL FOR A QUENCH WITH LONG-RANGE INTERACTIONS
In this paper we study the non-equilibrium dynamics of spatial entanglement in systems with long-range interactions, especially as they are realizable with variable range in ion traps. In this section we introduce the long-range transverse Ising model governing the time-evolution, and the measures of entanglement we compute.

A. Transverse Ising model
We consider the transverse Ising model with long-range interactions, described by the Hamiltonian Here, theσ α i denote the local Pauli matrices (α = x, z), J i,j is a general interaction matrix with potentially longrange interactions, and B is the transverse field. This Hamiltonian can be realized experimentally, e.g. with a string of trapped ions that are harmonically confined in a linear trap as depicted in Fig. 1. Using two stable (or metastable) electronic states of these ions as local spin representation at site i, |↑ i and |↓ i , it has been shown [21] that one can use collective couplings of these local states to motional degrees of freedom of the whole chain to produce the effective spin-model (1) (an example of J i,j for the ion trap experiment"case B" of Sec. IV is shown in Fig. 1b/c). Note that we will deal throughout with open boundary conditions, which are typical in ion chain experiments.
We define the local eigenstates ofσ z i as |0 i ≡ |↓ i and |1 i ≡ |↑ i with eigenvalues −1 and 1 respectively. We consider a quench experiment (see Fig. 1a), where the system starts in the fully polarized state |ψ 0 = M i |0 i , which is the ground-state for B(t = 0) → ∞. We are interested in the non-equilibrium dynamics of the many body-quantum state under a coherent evolution, i.e. ( = 1): We will concentrate on the case J i,j > 0 for all i, j. However, we note that nevertheless we obtain solutions for both the ferromagnetic and the anti-ferromagnetic case. Since we start in a state with a real probability amplitude in the spin-basis, the evolution for any ob-servableÔ z (withÔ † z =Ô z ), such as a density matrix of any sub-system is completely symmetric under the timereversal transformation t ↔ −t. This can be seen by the fact that where the cross terms have to vanish due to the real coefficients of the initial state and the fact that the expectation value must be real. Thus, the evolution of any observable is identical under both HamiltoniansĤ and −Ĥ. Therefore, the results we obtain for B > 0 with the anti-ferromagnetic (J i,j > 0) model are identical to the ferromagnetic model (J i,j < 0) with either negative field −B or a rotated initial state. In this paper we will show how the entanglement growth behavior changes with the strength of the magnetic field and the range of the interactions. Initially, we will idealize the interaction matrix, taking the form whereJ denotes the nearest-neighbor interaction strength and α is the decay exponent. This gives a good representation of the basic behavior of the interactions, but in real experiments there are typically small deviations from purely algebraic behavior of the interactions. In Sec. IV we will consider a full interaction matrix J i,j for real experimental parameters, as well as considering the effects of noise in the experiment.

B. Spatial entanglement
In characterizing the growth of spatial entanglement in the spin chain, we will make use of two complementary measures: The von Neumann entropy for a bipartite splitting of the chain in the center of the system, and the quantum mutual information between two distant spins. The former gives a measure of the overall entanglement buildup, and gives an idea of the complexity of the state being generated. The latter measure will give more detailed information as to how correlations propagate spatially, and will help us to characterize what part of the entanglement build-up is due to propagation of quasiparticles produced in the quench, and which part is due to direct interactions through long-range interactions. Both of these measures are accessible in experiments, though the mutual information is substantially less costly to measure (see Sec. IV B for more information).

Half-chain von Neumann entropy
Consider a chain of M spins as depicted in Fig. 1a,b. We can split this system into to halves L and R in the center of the system. In the case that the (pure) state of the composite system |ψ cannot be written as a product state of two states on the sub-systems L and R, i.e. |ψ = |ψ A |ψ B , we call the state entangled. The reduced density matrix of the sub-system L is defined via ρ L ≡ tr R (|ψ ψ|), where tr R denotes the partial trace over the system R. This density matrix will only be pure for a product state, and in the case of an entangled state the amount of bipartite entanglement is quantified by the von Neumann half-chain entropy of this matrix, which is defined as The time-dependent growth of the half-chain entropy summarizes the build-up of quantum correlations between two halves of the system. In a sense, it also underlies the complexity of numerical simulations of the quench when using matrix product state (MPS) representations. As we show in appendix A, the size of the MPS, represented by the bond dimension D has to grow exponentially as a function of time in the case that S vN grows linearly. Regimes of linear entanglement growth in time are thus important in demonstrating the power of a quantum simulator, in that realizing such regimes is a necessary requirement in order to observe dynamics that cannot be captured by state of the art numerical techniques over long timescales [13,14].

Quantum mutual information
An alternative measure, which gives more information on the distance of correlations, is the quantum mutual information between two distant spins i and j. In an experiment, this is also more straight-forward to measure than the von Neumann entropy for a bipartite splitting into two large blocks (see Sec. IV B), and clearly allows to distinguish different regimes of entanglement growth. The quantum mutual information is defined as Here, ρ i = tr k =i (|ψ ψ|) and ρ j = tr k =j (|ψ ψ|) denote the reduced density matrices of the single spins (obtained by tracing over all other spins k), and ρ ij = tr k =i,j (|ψ ψ|) the reduced density matrix of the composite system of the two spins.
Note that one has to be careful when interpreting the half-chain entropy and the quantum mutual information in an experiment in which the quantum state of the whole chain is in general mixed due to coupling to the environment and classical noise. In general the von Neumann entropy for each reduced density matrix is expected to increase compared to the zero temperature case [43]. We will consider these imperfections in Sec. IV.

III. ENTANGLEMENT GROWTH DYNAMICS
In this section we study the evolution of the entanglement after the quench. We identify three very different regimes: i) For relatively short-range interactions α 1 (depending on the system size) we find a linear growth of the half-chain entropy as a function of time, which we can understand in terms of free quasi-particle propagation within an effective Lieb-Robinson light-cone. ii) For long-range interactions α ∼ 0.8, 0.9, 1 we find a regime, where the half-chain entropy grows logarithmically. iii) For nearly infinite range interactions with α 0.2 we find rapid oscillations of the half-chain entropy around small values, which we can understand in an effective Dickestate model [44]. We treat the case i) in Sec. III A, ii) in Sec. III B 1, and case iii) in Sec. III B 2.
A. Entanglement dynamics for short-range interactions

Nearest-neighbor interactions
To understand the entanglement entropy growth behavior in this regime it is instructive to first revise the the case of nearest neighbor interactions, i.e. an interaction matrix (4) with a decay-exponent α → ∞, and discuss the dynamics of the quantities we study here. In this limit, the model Hamiltonian (1) becomes a standard transverse Ising model of the form which has been well studied in the literature. Note that since the spectrum of the Hamiltonian is symmetric under the exchange B ↔ −B, the dynamics will not only be identical under a change of the sign of the total Hamiltonian, but also under a change of the sign of B, and we therefore focus on B > 0 here. The model (7) can be diagonalized analytically [45] (see appendix B for more details). After performing a Jordan-Wigner transformation and diagonalizing the quadratic Hamiltonian in quasi-momentum space, the resulting diagonal model is a model of free fermions γ q , The γ q (γ † q ) are the annihilation (creation) operators for a fermionic quasi-particle with quasi-momentum q, which obey the anti-commutation relations {γ q , γ † p } = δ q,p . In the thermodynamic limit, i.e. for a chain of infinite length, the quasi-momenta become continuous −π < qa < π (a is the spatial separation between the spins) and the dispersion relation of the free particles is two-fold degenerate for q = −q = 0 and given by q = 2 (J − B) 2 + 4JB sin 2 (qa/2). The group velocity of quasi-particle excitations in this system is given by v g (q)/a = d q /d(qa). The maximum velocity of the quasi-particles gives rise to the Lieb-Robinson bound, which defines an effective light-cone for spatial correlations, outside of which the correlations are exponentially suppressed [40]. This sets an upper linear bound on the block entropy growth as we will see below. It is straight forward to calculate that the fastest particles move at a Lieb-Robinson velocity v R = max |v g | = 2aJ for B ≥ J and v R = 2aB for B < J.
Following [41], we can understand the entanglement distribution mechanism in model (8) as follows: In a coherent time-evolution, the initially excited state acts as a source for quasi-particle excitations. Pairs of the free fermions with quasi-momenta p and −p, which have been created at a certain point in space are entangled pairs. These pairs move freely through the system with a corresponding group-velocity v g and −v g , respectively. Parts of pairs which have been produced in block L and arrive in block R entangle the two blocks. An illustration of this mechanism is given in Fig. 2a. Thus, the arrival rate in block R for quasi-particles belonging to a pair created in block L is constant. Therefore, the increase of halfchain entropy is linear and we expect S vN = ηv g t, with some constant η. Since the group velocity is limited by the Lieb-Robinson bound, S vN ≤ ηv R t.
We can test this mechanism explicitly by making use of the boundary effects with open boundary conditions. Consider a quasi-particle pair, which has been created at the left edge of block L and moves at the Lieb-Robinson velocity. As soon as the right-moving quasi-particle arrives in block R, the linear entanglement increase has to break down since there are no more entangled pairs available to the left, which could further entangle blocks L and R. We can estimate the time at which this happens as t * = (M/2)/v R , which corresponds to t * = (M/4J) for B ≥ J and to t * = (M/4B) for B < J. In Fig. 2b we plot a comparison of this critical time with a numerical exact diagonalization simulation (ED, see appendix A for details) of the half-chain entropy evolution for B = J for increasingly large system sizes of 10 ≤ M ≤ 20 spins. We see that as expected exactly at the critical time, for each system size, the entropy starts to level off and remarkably reduces again after this maximum peak.
The quench experiment for B = J is special in the sense that this is a quench to the critical point of the quantum phase transition of Hamiltonian (7). We will now ask how the entanglement growth depends on B. In the limit of B → ∞ the initial state becomes an eigenstate of the system and no evolution will take place, i.e. S vN (t) = 0. On the other hand, in the limit of B → 0, the Hamiltonian has a spectrum with M degenerate levels, which are separated by an energy of ∼ 2J (spin flips). Thus in the latter case, we expect dynamics which is dominated by oscillations between those levels at a frequency scale given by J. Indeed for B = 0 it is straightforward to calculate analytically that S vN (t) = − cos 2 (Jt) log 2 (cos 2 (Jt)) − sin 2 (Jt) log 2 (sin 2 (Jt)). In Fig. 2c we find the expected behavior for B < J in an exact diagonalization simulation of a system with 20 spins. For increasing B, the oscillatory behavior of S vN (which fits the analytical result) breaks down and changes into a linear increase before boundary effects become important. Since the Lieb-Robinson velocity decreases with B for B < J, correspondingly the boundary effects shift to later times for smaller B (critical times t * are indicated as vertical arrows in Fig. 2c). It is interesting to note that in this case the maximum value of S vN can be actually larger than for B = J. In Fig. 2d we analyze the opposite case of B > J. Remarkably, we find that also away from the critical point the half-chain entropy growth becomes slower with increasing B. Below we will find that this also holds for finite-range interactions with α 2 and that the fastest entanglement growth follows precisely the point of the phase transition.

Finite-range interactions (α 1)
We now investigate the same behavior, but with algebraically decaying interactions. In Fig. 3 we demonstrate that the picture of entanglement distribution via entangled quasi-particle pair propagation also holds for a large range of finite α 1 (depending on the system size). We find that despite the existence of direct spin-spin interactions over all distances, it is the propagation of quasiparticles that dominates the dynamics of entanglement growth over this range of α values.
One marked signature of the linear half-chain entropy growth is that before boundary effects become important, the rate of the growth is essentially independent of the size of the system as shown in Fig. 3a. There we show the time-evolution of S vN after the quench for various decay exponents in the range ∞ ≥ α ≥ 1.5. The solid lines show an exact diagonalization (ED) simulation for a system of 20 spins, and the dashed lines are for M = 50 spins. We obtain the results for large systems using t-DMRG methods. Specifically, we use a MPO [32][33][34] of the Hamiltonian to time-evolve a MPS via a Runge-Kutta type method [46] (see appendix A for more details). In all cases with α > 2 we find that the two lines coincide and that the increase is linear. For α = 1.5 we find a slight change in the behavior in the sense that the results for M = 20 and M = 50 start to differ. In Fig 3b, we show an overview over the regime of linear half-chain entropy growth and its finite size scaling. Therefore as a function of system size and α we plot the error of a linear fit. In large systems we find that the linear increase (small error) breaks down at α ∼ 1. For smaller systems, boundary effects and finite-size effects become significant and the linear regime breaks down at larger α. Note that the change in behavior for large systems at α ∼ 1 can in a sense be understood, since this marks the point at which in the thermodynamic limit the sum in the interaction term in the Hamiltonian begins to diverge with increasing system size. In Sec. III B 1 we will show that the growth of S vN becomes logarithmic for α 1.
We can identify the regime of linear growth of S vN also by looking at the mutual information of distant spins. In Fig. 3c we plot the mutual information I 1,j , between site 1 and j, while increasing j = 4 . . . 10 for nearest neighbor interactions and for α = 2. We find that in both cases, for two particular spins that are separated by some distance, the mutual information remains nearly zero for a longtime, until it suddenly peaks at a time corresponding to the arrival of a quasiparticle pair originally produced on a site between the two spins, which then entangles the two spins. After the quasiparticles pass, the mutual information remains at a value slightly greater than zero (barely visible). This is consistent with the nearest-neighbor model, where entanglement is also distributed by freely propagating quasi-particles. Furthermore, we find that the time of the arrival of the wave at site j is consistent with quasi-particles moving at the Lieb-Robinson velocity for nearest neighbor interactions. Since the two sites become entangled once a quasi-particle that has been created in the middle of the two sites first arrives at both spins, this time is given by t j = (j −2)/2v R = (j −2)/2aJ and is shown as black bars in Fig. 3c). In the lower panel we find that even for α = 2, despite the rather long-range interactions, the characteristic behavior is the same. For longer-range interactions we find that the wave moves more slowly. Furthermore we note that for α → ∞ one finds a much more "diffusive" behavior in the sense that the peaks of mutual information broaden and become smaller with distance.
The peaks in Fig. 3c allow us to extract an effective velocity of the "wave" of entangling quasi-particles, v I . We do this by fitting a line to the position of the peaks (in time) as function of the distance j to obtain 1/v I . As we show in Fig. 3d, the rate of the half-block entropy increase (linear fit 1 ≤ tJ ≤ 3) is directly related to v I . We find that for α 5, S vN (t) = ηv I t/a with η ≈ 0.253. For α 3, the proportionality constant starts to depend on α and η < 0.253. This means that for increasing range of interaction both the half-chain entropy growth rate and v I reduces, and that v I decreases less strongly. We note that the in general η also depends on B.
We emphasize that the fact that the entanglement growth mechanism is directly reflected in the timedependence of the mutual information between two distant spins is very important for experimental observations. Instead of having to reconstruct 2 (M/2) × 2 (M/2) density matrix elements of a large block via quantum state tomography, the growth behavior of the half-chain entropy can be directly verified by measuring only 4 × 4 density matrices for a system of two composite spins. In Sec. IV B we will show how the measurement further simplifies for the particular quench we consider here.

Ground-state phase transition point and entanglement growth
In this section we study how the linear growth of S vN depends on the value of B, and the decay exponent α. We focus on the caseJ = 1, B > 0. For ground states, as for the nearest neighbor transverse Ising model, the long-range model undergoes a quantum phase transition from an antiferromagnetic to a paramagnetic phase at a critical field B c for all decay exponents α [47,48]. For example, for α = 3 it has been calculated that B c ≈ 0.83 [49]. As shown in Ref. [49], we can estimate the point of the phase transition by locating the discontinuity of the transverse magnetization, i.e. the jump in d 2 m z /dB 2 , where m z = m σ z m . In order to do this for moderately large systems, we use a MPO Lanczos diagonalization for 100 spins. In Fig. 4, we compare this ground-state phase transition point for a moderate system size to the linear growth rates of the von Neumann entropy in a small systems of only 20 spins as function of α and B. It is remarkable that the point of the ground-state phase transition is not only reflected in the scaling behavior of the block entropy in the ground-state [48], but we also find that in the evolution following the quench from B(t = 0) = ∞ → B, the growth rate of the von Neumann entropy as a function of time is largest at the critical points B c . Since the entangling quasi-particles are bounded by the Lieb-Robinson bound, this effect is independent of the system size up to times t * when boundary effects limit the quasiparticle propagation.

B. Entanglement dynamics for long-range interactions
In this section we study the entanglement growth for very long-range interactions with α ≤ 1. In this regime the picture entangling quasi-particles that move freely within a light-cone breaks down, and instead distant parts of the system can become almost instantaneously entangled based on direct interactions. Here we observe that for α ∼ 0.8, 0.9, 1, the half-chain entropy still can increase steadily as a function of time for our quench, but that the increase becomes logarithmic instead of linear. When further increasing the range of interactions for α 0.2, we find a regime where S vN oscillates rapidly around small values. We understand this behavior via an effective model in a basis of Dicke states [44] for infinite range interactions α = 0.

Logarithmic entropy growth
When increasing the range of interactions, eventually the linear growth of S vN breaks down, and the growth becomes slower than linear as depicted in Fig. 5a. In Fig. 5b, we find that this growth becomes logarithmic.
For very long-range interactions, the timescale of the dynamics is dominated by the interaction energy term  Fig. 5b we plot the evolution of the von Neumann entropy in these units for system sizes of M = 30, 40, 50, for α = 0.8, 0.9, 1 and for B = 0.7J, 1J, 1.3J. If ideally, the entropy increased logarithmically without oscillation, S vN = log 2 (Ct + 1), with some constant C. On an exponential scale, i.e. by plotting 2 SvN(t) , we would see a straight line. In Fig. 5b we see that indeed after some initial behavior, the von Neumann entropy shows oscillations, around a straight line. It is remarkable that for a fixed value of B, independently of the system size and for all α = 0.8, 0.9, 1 we find roughly the same constant C in units of ∆τ −1 . With decreasing B, i.e. for quenches, which put an increasing amount of energy into the system, the constant increases. For interactions with α 0.7, the oscillations become more dominant, so that the logarithmic increase is hard to verify.
Also the time evolution of the mutual information between distant spins shows a completely different qualitative behavior for α 1 than for α 1. In Fig. 5c we show the evolution of I 1,8 after the quench for a system of 20 spins (B =J) for 2 ≥ α ≥ 0.2. In the upper panel of Fig. 5c we find that the incoming wave picture breaks down when decreasing α from α = 2 to α = 1. In the latter case we find a mixed behavior, where a wave peak is still roughly visible around tJ ∼ 3, however also peaks for very short times appear. These indicate that due to the long-range part of interaction, distant parts of the system become rapidly entangled with an immediate increase in correlations after the quench. When further increasing the interaction range (lower panel), these contributions to I 1,8 become dominant and the quasi-particle peak disappears completely. While for α ∼ 1, I 1,8 still shows some slow overall increase as function of time, in the case of nearly infinite range interactions (α 0.2) we only find rapid oscillations around a constant value of I 1,8 ∼ 0.1. In Fig. 5d we show the crossover from the short-range physics regime to the long-range physics regime, by plotting the maximum of the short-time mutual information I ≡ max[I 1,8 (tJ < 1)], as a function of α. We show results for different system sizes. Remarkably, the curves are identical. What we find is that the frequency of the oscillations of the mutual information for α 1 does depend on the system size, but the height of the short time peak depends only on the range of the interactions.

Entanglement dynamics for infinite-range interactions
To understand the oscillatory behavior for decay exponents of α 0.2, it is instructive to consider the case of infinite range interactions, i.e. α = 0. In this limit each spin interacts with equal strength with all others and the ("mean-field") Hamiltonian reads: As in the nearest-neighbor case, this limit is analytically exactly tractable. We can introduce effective spin-M/2 operators S x,y,z ≡ M iσ x,y,z i . With these, we can we rewrite Hamiltonian (9) aŝ  The eigenstates of Hamiltonian (10) are given by Dicke spin-M/2 states, which are defined as Here n 0 , n 1 denotes the number of spins down/up, respectively, and S is the symmetrization operator. Best illustrated by an example: For a system of 4 spins, we would have an effective spin-2 model, and the Dicke state |S = 2, m S = −1 = (|0001 +|0010 +|0100 +|1000 )/2. In this picture the quench experiment is equivalent with a free evolution under Hamiltonian (10) with the initial state |ψ 0 = |S = M/2, m S = −M/2 . It is straightforward to calculate the half-chain von Neumann entropy for an arbitrary Dicke state [50] (see appendix B for more details). One finds that where p l are combinatorical factors depending only on the number of single up spins of the corresponding Dicke state and 0 ≤ l ≤ M/2. Example results are plotted in Fig. 6a, and it is important to note that simply due to the fact that the sum in (12) contains a maximum of M/2+1 terms, the entropy is bounded by S vN ≤ log 2 (M/2 + 1). In our quench, the time-dependent state will assume a superposition of Dicke states, |ψ(t) = m c m (t)|M/2, m with c m (t = 0) = δ m,−M/2 . In the small Dicke Hilbert-space the time-evolution can be easily numerically simulated, and we can at any instance in time construct the reduced density matrix for the left system. From this density matrix, the von Neumann entropy can be readily calculated and examples are shown in Fig. 6b. Again, we find simply by noting that the dimension of the M/4 Dicke subspace is M/2+1, that for an arbitrary Dicke state superposition, the von Neumann entropy is limited by S vN ≤ log 2 (M/2 + 1). However, for our specific experiment a tighter bound can be found. Since the Hamiltonian (10) only couples spin-states with m S = ±2, the time-dependent coefficients c m (t) can only be non-zero for states with m = −M/2, −M/2+2, . . . , M/2. Assuming an even number of spins, the entropy is then limited by S vN ≤ log 2 (M/4+1), which is in agreement with the exact von Neumann entropy evolution as shown in Fig. 6b.

A. Entanglement growth for realistic experimental parameters
In this section we ask to what extent the effects shown in the previous sections are experimentally observable in ion traps. Therefore, we consider two experimentally realistic full interaction matrices J i,j , which show the characteristic behavior of linear entanglement growth and logarithmic one. In a case A, over short distances, the averaged interactions decay as α < 1 (logarithmic growth regime), and in a case B, they decay as α ∼ 2 (linear growth regime), as depicted in Fig. 7a. We define the energy unit by the largest element of J i,j , which we denotē J in this section.
We consider a linear chain of 20 ions which interact via the mechanism described in [29,30]. In summary, a force is applied which couples the electronic 'spin' state of the ions to the spectrum of closely spaced (nearly frequency degenerate) vibrational modes transverse to the ion string. By setting the driving force to be far offresonant the phonon states can be adiabatically eliminated, allowing an analogue simulation. While the simulations presented are for 40 Ca + ions driven with bichromatic laser fields, similar interaction matrices can be derived in a number of other systems. The exact experimental parameters considered are given in appendix C.
We show the evolution of the half-chain entropy and the mutual information I 1,5 in Fig. 7b and 7c, respectively. As expected from the averaged decay of the interactions in both cases A and B and from the previous discussions, we find that for the case B, S vN increases linearly in time, while for A, the growth behavior is logarithmic. Accordingly, we find the behavior of the mutual information as we have found it for the case of homogeneously decaying interactions: In the case B, the initial mutual information is zero and a peak appears at a time ∼ 2.8/J. In case A, in the logarithmic regime, we find an instantaneous increase of I 1 due to the long-range part of the interactions, which entangles distant parts instantaneously.
In a realistic experiment setup, the string of ions will be subjected to noise. Here, we consider the two most significant imperfections: i) Fluctuations in the energy splitting of the electronic states used to encode a spin (e.g. due to ambient magnetic field fluctuations) and ii) fluctuations in the coupling strength between the spindependent force mechanism and the ions (e.g. due to laser intensity fluctuations). Noise of the form i) will lead to a correlated rotation of the qubits around the zaxis, while noise of type ii) causes a stochastic fluctuation of the overall interaction strengthJ. We idealize both cases as white noise fluctuations ξ(t), ξ(0)ξ(t) = s α δ(t) with strength s α (the bar denotes the time-average). The evolution of the state can then be described with a Stratonovich stochastic Schrödinger equation [51,52] where dW (t) is the Wiener increment, andL α is the noise (jump) operator. For case i) Equivalently we can derive a master equation for the evolution of the full density matrix In the long-time limit the master equation drives the system into a state that commutes with the jump operators. For noise of type i) this is a state diagonal in the z basis, and for type ii) a state diagonal in the x basis. For the time-scales of the experiment, the dynamics consists of a complicated interplay between the coherent evolution and the dissipative part. In general, we expect that noise of type i) leads to a global dephasing in the sense that it will reduce the purity of the full state and thus result in a slightly higher measured entropy, whereas noise of type ii) can be more complicated because an overall fluctuation ofJ acts with different strength between different spins according to J i,j . We can simulate the evolution of the master equation numerically by evolving the stochastic Schrödinger equation in time using a first-order semiimplicit method with strong order 1.0 convergence [51] and statistically averaging over a large amount of trajectories.
In Fig. 7b,c we find that as expected the noise adds an additional entropy growth as function of time for both, the half-chain entropy and the mutual information. However, the underlying entropy features, which arise from the entanglement-build up in the coherent evolution remain clearly visible. For example, instead of the shorttime initial mutual information being zero in the regime of the linear S vN growth regime, for noise of type i) we find an overall increase of I 1,5 as function of time. Nevertheless, the quasi-particle peak clearly remains observable even in the presence of the noise. In general we find that the overall entropy growth, which is induced by the fluctuations on B is larger than the one induced by fluctuations on J i,j . Especially we find that the mutual information between distant spins is very robust against noise on the coupling strength due the decay of J i,j with distance. In the case A we find that for long times, the characteristic logarithmic growth eventually breaks down due to the entropy increase from the noise, however for  tJ 3 it remains observable. In general these results suggest that the mutual information is a very robust experimental measure for the entanglement growth behavior of the systems, which furthermore can be easily extracted from experimental data as we show in the next section.

B. Measurements of block entropies and mutual information
We will now briefly review, how the entanglement measures we used in this paper can be measured experimentally. To calculate the von Neumann entropy of a subsystem A of l spins (which do not have to be next to each other), one can simply measure the reduced density matrix of that block. While the process of measuring this matrix is known as quantum state tomography, we show that for our particular experimental setup this tomography simplifies significantly.
The reduced density matrix of A, after tracing over the remaining system is where bold greek symbols denote the set of indices for the sub-system of spins, i.e. all 2 M A binary representations of M A spins: α = (α 1 , α 2 , . . . , α M A ) with α k ∈ {0, 1}. The diagonal elements ofρ A can be easily measured. They are the probabilities for finding the spin-combination α, p α . The off-diagonal elements are more challenging but reduce to the measurement of spin-correlations. For our experimental situation we can make use of the fact that for for any time-evolved state |ψ(t) = α c α (t)|α where denotes the sum modulo two. This is easily verified by the fact that the matrix elements of any power of Hamiltonian (1), α|Ĥ n |β = 0 for M k α k = M k β k . Thus, since we start in the state with c α (t = 0) = δ α,0 , the time-evolution operator exp(−itĤ) = n=0 (−itĤ) n /n! can only produce states with non-zero coefficients c α (t) for which M k α k = 0. Thus, half the elements of any reduced density matrix calculated from the time-evolved state |ψ(t) will always remain zero. The remaining spin-spin correlations consist only of σ x and σ y terms.
We illustrate this here for the example of a sub-system of a single (l = 1) and two (l = 2) spins, and show how to reconstruct the corresponding density matrix. In the case of a single-spin, the density matrix will be diagonal: Here the off-diagonal part completely vanishes since trivially for allρ β α with α = β, α ⊕ β = 1. For a block of two spins, the density matrix becomes The six density matrix values are all that have to be measured. To experimentally obtain the off-diagonal elements, one has to measure the real and imaginary parts of ρ 11 00 and ρ 10 01 . This can be done by expanding those elements into spin-spin correlations via Re(ρ 11 00 ) = (σ x ⊗σ x −σ y ⊗σ y ) /4 Im(ρ 11 00 ) = (σ x ⊗σ y +σ y ⊗σ x ) /4 Re(ρ 10 01 ) = (σ x ⊗σ x +σ y ⊗σ y ) /4 Im(ρ 11 00 ) = (σ y ⊗σ x −σ x ⊗σ y ) /4, which means that 4 spin correlations have to be measured. Thus, for the whole density matrixρ 2 , a total of 8 expectation values have to be experimentally determined. Note that this is a simplified version of full quantum state tomography for two qubits [20], however a full state tomography, i.e. a measurement of all density matrix elements might be useful to quantify the deviations of the other elements from zero and therefore get a measure for experimental deviations.
From the density matrices one can directly extract the corresponding von Neumann entropy. Note however, that one can also calculate Renyi entropies of arbitrary order n. These entropies are defined as S [n] (ρ) = log 2 [tr(ρ)]/(1 − n) and can serve as lower bound to the von Neumann entropy. In an analogue fashion, one can use these entropies to define a mutual information A generalization of the density matrix measurement to larger blocks is straightforward. In particular for a block length of l sites, one has to measure a total amount of [2 2l−2 + 2 (l−1) ] density matrix values. For example, to clearly observe a linear entanglement increase for α = 2, one would have to simulate the system until a time of ∼ 2/J (cf. Fig. 3a). Before finite size effects of the block size play a role (cf. Fig. 3c) one therefore has to consider a block with l = 8. However, it is important to note that the mutual information can always be extracted by only using blocks of l = 1 and l = 2.
Alternatively, entanglement entropies could be directly measured by incorporating copies of the same state and using "beam-splitter" operations between those copies [53]. These multiple copies could for example be realized by performing the quench identically for co-trapped strings in micro-traps. Alternatively one could use a single large string and hiding pulses to effectively realize two copies.

V. CONCLUSION & OUTLOOK
We have studied the dynamical evolution of entanglement in quench experiments in ion traps. We found that a regime of linear half-chain von Neumann entropy growth is present even for relatively long-range interactions with interaction decay exponents of α 1. The growth rate is sensitive to the underlying low-energy spectrum and notably largest at the point of the quantum phase transition, which varies with changing α. For longer-range interactions, we find a regime of logarithmic entropy growth (α 1), and for infinite range interactions (α = 0) the entropy remains bound by a small constant and oscillates rapidly. We showed that mutual information between distant parts of the system can be experimentally measured and used to distinguish the different regimes.
The entanglement entropy growth behavior has important implications for t-DMRG/MPS/MPO algorithms on classical computers, since it underlies the complexity of these algorithms. For a specific experimental example we have shown that this regime is in reach for an experimental quantum simulator. This means that we can find a regime in these systems that fulfill a necessary condition realization of a quantum simulator regime where state of the art numerical simulations on classical computers can become inefficient.
We emphasize that ion trap experiments are not the only experimental realization in which the entanglement dynamics studied here could be observed. For example, exactly the same spin model [54] and also more complicated models with long-range interactions can be realized in system with polar molecules [55,56] or Rydberg atoms [57][58][59][60] in optical lattices. These systems would have the disadvantage that the decay exponent is not directly tunable, however they might have the advantage that halfchain entropies might be easier to measure directly, by employing schemes which rely on the preparation of multiple copies in an optical lattice [53].
In the final stages of preparing our manuscript we became aware of some related work [61], in which a local quench in an Ising model with algebraically decaying interactions is studied, together with the resulting spread of quasiparticles. Though our quenches are qualitatively different, and large-scale entanglement growth cannot be observed in a local quench, Ref. [61] has interesting parallels with what we observe here, and identifies similar parameter regimes for dynamics with long-range and shortrange interactions. ble for large system sizes. The von Neumann entropy is defined as where in the second equation we defined the eigenvalues of ρ L as λ α and introduced the number of non-zero eigenvalues, χ L ≤ dim(ρ L ). χ L is called the Schmidtrank and can be considered as an entanglement measure itself. Note that also the maximum possible Schmidt rank grows exponentially as a function of M , and correspondingly the maximum possible von Neumann entropy grows linearly as S [max] vN = M/2. However the Schmidt rank for states as they occur in typical experiments can be much smaller than the dimension of ρ L . The approximation, which is made in numerical DMRG/MPS algorithms (see below) consists therefore of truncating the Schmidt rank (for all possible bipartite splittings) at a maximum value, which is called the bond dimension D, thus effectively limiting the von Neumann entropy to S [MPS] vN ≤ log 2 (D). The error made is called the truncated weight, D = χ L i=D+1 λ α . The truncation made in typical DMRG simulations is very significant. For example, considering a string of 50 spins, the dimension of ρ L is given by 2 25 . If one represents a state of this system by a MPS with a bond dimension of D M/2 = 1024, the size of the effective Hilbert space for one half of the system is still only ∼ 0.003% of the full Hilbert space. However, it is quite remarkable that this approximation is in many cases quasi exact. The reason for that is that most low-energy states of physical systems in nature turn out to be in fact very slightly entangled. If for example in Hamiltonian (1) we restrict the interaction to nearest neighbors (standard transverse Ising model), then it can be proven that for the critical model, i.e. for J = B where for an infinite system the energy gap from the ground state to excited states disappears, the entanglement entropy scales with the block length as S [GS] L ∼ log 2 (L). If the system is gapped, the entanglement entropy scales as S L ∼ const., i.e. it obeys an "area law" [62,63]. Thus, in the worst case (critical model), the bond dimension only has to grow linear with the system size, and therefore ground-states can be easily calculated up to a quasi exact precision with DMRG/MPS for systems of hundreds of spins.
It is obvious that therefore in a time-evolution simulation, whether the simulation of the system over long times with t-DMRG methods is in practice possible or not depends on how fast the entanglement grows. If in our quench experiment with ions the von Neumann entropy grows linearly as a function of time, in order to keep D small, D has to grow exponentially as a function of time. The computational resources to store a MPS grows thus exponentially with the time the system is to be simulated. This becomes prohibitively expensive for large system sizes.

Exact diagonalization
Despite the exponential growth of the Hilbert space, quantum systems of moderate size can still be diagonalized exactly, simply by exploiting the sparseness of typical Hamiltonians. For example, for 20 spins the full Hamiltonian is a 2 20 × 2 20 matrix, which is too large to even store in the memory of current computer hardware. However the amount of non-zero elements of Hamiltonian (1) is only O(2 20 ), even for a full interaction matrix. Therefore, one can use Krylov subspace projection techniques [64] to evaluate the matrix exponential of the Hamiltonian matrix, as well as semi-implicit first order methods to propagate a state vector for systems of 20 spins in time.
A MPS is defined as decomposition of the complex amplitudes of the full quantum state of a lattices system, |ψ = {i k } c i1,i2,...,i M |i 1 |i 2 . . . |i M (M sites with local basis states {|i k }) into a matrix product. Specifically, we define a MPS in its canonical form as Here, the A  [38]). In general the size of the matrix dimensions that are required to represent a certain state exactly are given by the Schmidt rank χ i for the bipartite splitting between sites i and i+1. Limiting all matrix dimension by the bond-dimension D limits the maximum allowed von Neumann entropy to S = 0, and it is readily seen that the MPS reduces to a non-entangled product state.
We can study systems with long-range interactions numerically by making use of matrix product operators (MPO), which are as MPS a decomposition of the now real 4 Mdimensional operator tensor of the full Hamiltonian H = {i k },{j k } o j1,j2,...,j M i1,i2,...,i M |i 1 |i 2 . . . |i M j 1 | j 2 | . . . j M | into a matrix product. The long-range interaction Hamiltonian (1) with a decaying interaction J ij =J/|i − j| α can be up to a very good approximation implemented as MPO with relatively small bond-dimension, which can be achieved by expanding the power-law decay function into a sum of exponentials [32][33][34]. With the Hamiltonian in MPO form it is then possible to implement time-evolution algorithms, using e.g. a Runge-Kutta type evolution scheme [46]. Alternatively one can also use the original adaptive t-DMRG methods [35][36][37] and introduce swap-gates, which interchange indices in a MPS. This has the advantage that arbitrary interaction matrices J ij can be implemented.
To calculate ground-states, one can either evolve a MPS in negative imaginary time, or one can construct the local representations of the Hamiltonian expectation value by leaving the indices on a particular site open and contracting the remaining tensor network (see e.g. [38]). We then use a local iterative Lanczos solver to find the local MPS matrix, which minimizes the corresponding energy. By sweeping through the system from site to site, this is a very efficient method to find the overall MPS ground-state for large systems.
We checked the validity of all our results by comparing different methods and confirm convergence in the bond dimension by running multiple calculations with increasingly large D.
The p l can be found using combinatorical arguments, . (B7) The von Neumann entropy of half the chain can be trivially extracted from the p l , since they are simply the eigenvalues of the reduced density matrix for half the chain. Therefore, S vN = − l p l log 2 (p l ). In our quench experiment, initially the ion chain is in the Dicke state |M/2, −M/2 , and the subsequent timeevolution will rotate the state vector in the Dicke manifold, and therefore prepare a time-dependent superposition Numerically determining the c m (t) is easily achieved for large systems, since the dimension of the Hilbert space is only given by M + 1. However, when inserting the formula (B6) into (B8), the resulting decomposition is not a proper Schmidt-decomposition anymore, since the"Schmidt values" can be complex now. Therefore one has to construct the full reduced density matrix by tracing over one half of the system, This matrix can then be easily diagonalized numerically, which allows us to time-dependently calculate the halfchain entropy. For the same reason as discussed in Sec. IV B, we find that half the coefficients c m (t) will always remain zero (The Hamiltonian only couples terms with m ↔ m ± 2. Thus the amount of non-zero eigenvalues of the matrix (B10) is limited by M/4 + 1 (for even M ). This gives the upper bound of the entropy, S vN ≤ log 2 (M/4 + 1).
Appendix C: Realistic experiment with 20 ions In section IV the entanglement growth is analyzed in experimentally realistic situations. Here we give details on the exact experimental parameters considered. In each case we consider a string of 20 ions in a 3D harmonic trap with highest transverse motional frequency ω t =2π×4.9 MHz, and variable lowest axial frequency 0.1<ω z /2π<0.5 MHz. A state-dependent driving force is simulated at a fixed detuning δ=+2π×80 kHz from ω t . These three parameters are sufficient to completely determine the form of the interaction matrix [29,30].
In order to vary the interaction strength we choose to vary ω z , which has the effect of changing how closely spaced the transverse mode are.
The absolute value of the interaction strengthJ is determined by the specific choice of driving force mechanism, driving strength, ionic species, and electronic transition used to encode the spin. We consider an optical transition, at 729 nm, between the S 1/2 , m=1/2, ground and metastable D 5/2 , m=3/2, states in 40 Ca + . Copropagating bichromatic laser fields at 729 nm drive the interaction, with a Rabi frequency of Ω=2π × 0.5 MHz (if put on resonance with the electronic spin). For the highest transverse mode cooled to the ground state, the coupling strength on the upper vibrational sideband of the transition is therefore η t Ω=2π × 22 kHz, where η t =0.044 is the Lamb-Dicke parameter. In the case δ η t Ω, the assumption that the phonon states can be adiabatically eliminated holds. In our case δ ≈ 4η t Ω. This approximation could be improved at the expense of a slower overall interaction strengthJ.
For the cases shown in figure 7a, b, and c, we find overall spin-spin coupling ratesJ/2π=2π × 0.5, 2π × 0.4 and 2π×0.3 kHz. These rates correspond to quantum dynamics observable on time scales of a few ms, which compares favorably with typical decoherence times of several 10's of ms that are typically achieved in these experiments.