Time-Resolved Observation of Thermalization in an Isolated Quantum System

Time-Resolved Observation of Thermalization in an Isolated Quantum System

Govinda Clos govinda.clos@physik.uni-freiburg.de Physikalisches Institut, Albert-Ludwigs-Universität, Hermann-Herder-Straße 3, 79104 Freiburg, Germany.    Diego Porras Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, United Kingdom.    Ulrich Warring Physikalisches Institut, Albert-Ludwigs-Universität, Hermann-Herder-Straße 3, 79104 Freiburg, Germany.    Tobias Schaetz Physikalisches Institut, Albert-Ludwigs-Universität, Hermann-Herder-Straße 3, 79104 Freiburg, Germany.
July 31, 2019

We use trapped atomic ions forming a hybrid Coulomb crystal, and exploit its phonons to study an isolated quantum system composed of a single spin coupled to an engineered bosonic environment. We increase the complexity of the system by adding ions and controlling coherent couplings and, thereby, we observe the emergence of thermalization: Time averages of spin observables approach microcanonical averages while related fluctuations decay. Our platform features precise control of system size, coupling strength, and isolation from the external world to explore the dynamics of equilibration and thermalization.

37.10.Ty, 37.10.Jk, 03.65.-w, 05.30.-d

How does statistical mechanics emerge from the microscopic laws of nature? Consider, for example, a finite, isolated quantum system: It features a discrete spectrum and a quantized phase space, its dynamics are governed by the linear Schrödinger equation and, thus, remain reversible at all times. Can such a system equilibrate or even thermalize? Progress in the theory of nonequilibrium dynamics and statistical mechanics sheds light on these fundamental questions. It has been shown that individual quantum states can exhibit properties of thermodynamics depending on entanglement within the system Popescu et al. (2006); Goldstein et al. (2006); Reimann (2008); Linden et al. (2009); Cramer and Eisert (2010); Short and Farrelly (2012); Riera et al. (2012). While the entire system may very well be described by a pure state, any small subsystem and related local observables may be found in a mixed state due to disregarded entanglement with the rest of the isolated system, i.e., the large environment. Further, it is predicted that even any individual many-body eigenstate of a nonintegrable Hamiltonian yields expectation values for few-body observables that are indistinguishable from microcanonical averages Deutsch (1991); Srednicki (1994, 1999); Rigol et al. (2008); Polkovnikov et al. (2011); D’Alessio et al. (2016). This conjecture has been extensively studied by numerical simulations of specific quantum many-body systems of moderate size, exploiting available computational power Rigol (2009); Santos and Rigol (2010); Beugeling et al. (2014); Kim et al. (2014). Recently, there have been first experiments in the context of thermalization in closed quantum systems with ultracold atoms Gring et al. (2012); Trotzky et al. (2012); Langen et al. (2015). However, fundamental questions on the underlying microscopic dynamics of thermalization and its time scales remain unsettled Cazalilla and Rigol (2010); Polkovnikov et al. (2011); Eisert et al. (2015).

Trapped-ion systems are well suited to study quantum dynamics at a fundamental level, featuring unique control in preparation, manipulation, and detection of electronic and motional degrees of freedom Wineland et al. (1998); Leibfried et al. (2003); Blatt and Wineland (2008); Myerson et al. (2008); Kienzler et al. (2015); Tan et al. (2015); Ballance et al. (2015). Their Coulomb interaction of long range permits tuning from weak to strong coupling Porras and Cirac (2004). Additionally, systems can be scaled bottom up to the mesoscopic size of interest to investigate many-body physics Schneider et al. (2012); Monroe and Kim (2013); Ramm et al. (2014); Mielenz et al. (2016).

In this Letter, we study linear chains of up to five trapped ions using two different isotopes of magnesium to realize a single spin with tunable coupling to a resizable bosonic environment. Time averages of spin observables become indistinguishable from microcanonical ensemble averages and amplitudes of time fluctuations decay as the effective system size is increased. We observe the emergence of statistical mechanics in a near-perfectly-isolated quantum system, despite its seemingly small size.

The dynamics of our system are governed by the Hamiltonian Porras et al. (2008); sup ()


The spin is described by Pauli operators and denotes the reduced Planck constant. The first term can be interpreted as interaction of the spin with an effective magnetic field , lifting degeneracy of the eigenstates of , labeled and , while the second drives oscillations between these states with spin coupling rate . The sum represents the environment composed of harmonic oscillators with incommensurate frequencies , and the operators annihilate (create) excitations, i.e., phonons, of mode . The last term describes spin-phonon coupling via spin flips, , accompanied by motional (de)excitation which is incorporated in


at a strength tunable by and the spin-phonon coupling parameters . Expanding in series permits restricting to linear terms for values (weak coupling). For , as in our experiment, higher order terms become significant (strong coupling), allowing the system to explore the full many-body set of highly entangled spin-phonon states. This regime is well described by full exact diagonalization (ED) only, since the discrete nature of the bosonic environment of finite size hinders standard approximations applicable to the spin-boson model considering a continuous spectral density Leggett et al. (1987); sup ().

To study nonequilibrium dynamics of expectation values , consider an initial product state , where the spin is in a pure excited state, and the bosonic modes are cooled close to their motional ground states (average occupation ). With this, we ensure that energies of spin and phonons remain comparable to enable the observation of the coherent quantum nature of the dynamics which creates entanglement of spin and phonon degrees of freedom. Because of the coupling, the spin subsystem is in a mixed state for , even though the entire system is evolving unitarily. When thermalization occurs, any small subsystem of a large isolated system equilibrates towards a thermal state and remains close to it for most times Srednicki (1999); Linden et al. (2009).

The so-called eigenstate thermalization hypothesis provides a potential explanation for the emergence of thermalization in an isolated quantum system. It can be phrased as a statement about matrix elements of few-body observables in the eigenstate basis of a many-body Hamiltonian Deutsch (1991); Srednicki (1994, 1999); Rigol et al. (2008); D’Alessio et al. (2016); Polkovnikov et al. (2011). Within this conjecture, infinite-time averages of expectation values of these observables agree with microcanonical averages. A mathematical definition of this hypothesis and further information are given in the Supplemental Material. Based on Refs. Deutsch (1991); Srednicki (1994), to interpret experimental results, we assume that a coupling distributes any of the energy eigenstates of an uncoupled system over a large subset of the energy eigenstates of the coupled system , i.e., sup (). Further, we consider that these participating states lie within a narrow energy shell around the energy of Rigol et al. (2008); D’Alessio et al. (2016). As introduced in Refs. Srednicki (1999); Popescu et al. (2006); Linden et al. (2009); Short and Farrelly (2012), an effective dimension of the subset, , provides an estimation for the ergodicity of a system. It has been shown that mean amplitudes of time fluctuations of expectation values are bounded by  Linden et al. (2009); Short and Farrelly (2012).

Figure 1: Complexity of the Hamiltonian studied numerically. Parameters are MHz, . (a) Dimension of truncated Hilbert space for corresponding fractions of initial-state population lying within (solid lines). For , for example, 85% lie within (circle). We derive up to (dashed line). (b) Choosing and varying (dashed line), we can tune the spin-phonon coupling into resonance with different modes (sketched at the bottom) and boost the system size. Note, that the actual number of participating states is much larger than the normalized quantity ; see Eq. (5). (c) For fixed [cf. dashed line in (b)], increases significantly with . This enables the systematic investigation of equilibration and thermalization depending on the system size. Error bars show systematic numerical uncertainties sup ().

Correspondingly, for our system, we exploit these predictions for infinite-time averages, both of spin expectation values,


and of their time fluctuations,


To this end, we need to quantify the complexity of the dynamics induced by the coupling. Hence, we extend existing definitions of to a weighted effective dimension sup ()


for . Here, in contrast to , the statistical average over is performed after calculating the inverse participation ratio for each pure state in the mixture sup (). Thereby, also incorporates the number of participating states, but is normalized to for the uncoupled system.

Figure 2: Measured unitary time evolution . Experimental results (black dots, error bars: 1 s.d.) for compared to full ED (solid lines). We exclude numerical results for due to their large systematic uncertainties. Oscillations (time fluctuations) of high amplitude during the transient duration , are driven by the evolution of towards the ground state of . (a) For and increasing , excitation is coherently exchanged with a growing number of modes resulting in revivals at (shaded areas). (b) For , expectation values fluctuate around a negative offset. Revivals and this nontrivial bias emphasize the coherence of the dynamics. (c) Histograms of experimental measurements sample the probability distribution which underlies . Here, we show these for , , and to exemplify the quantities and .

Throughout our Letter, we estimate numerically. depends on , , , , and . We approximate the latter by truncating the Hilbert space to , choosing a phonon number cutoff , such that sup (). For a given computational power and increasing , the description of the initial-state population by becomes less representative, leading to increasing systematic uncertainties, illustrated in Fig. 1(a). Here, the exponentially growing complexity becomes evident: is required to achieve for . Figure 1(b) highlights the experimental controllability of . At , the strong coupling to numerous modes leads to a maximum in . For large , the spin can get close to resonance with few modes only, the latter composing a comparatively small environment. Further, the range of accessible values of grows with increasing ; see Fig. 1(c).

We experimentally implement the single spin by two electronic hyperfine ground states of Mg and add up to four Mg to engineer the size of the bosonic environment spanned by (number of ions) longitudinal (axial) motional modes. For details on our experimental setup, see Refs. Schaetz et al. (2007); Clos et al. (2014). First, we prepare the spin state, , by optical pumping and initialize the phonon state, , by resolved sideband cooling Leibfried et al. (2003) close to the ground state. In calibration measurements we determine that the modes are in thermal states with , which effectively enhances . Next, we apply the spin-phonon interaction by continuously driving Raman transitions with spin coupling rate for variable duration , where is the controllable detuning from resonance Porras et al. (2008). Finally, we detect the spin by state-dependent fluorescence. We choose to record , while we numerically check that feature similar behavior. To study dynamics for a large range of , we choose 95 parameter settings: We set MHz which corresponds to an effective spin-phonon coupling parameter for . For each , we use a fixed MHz and vary from up to  sup ().

In Fig. 2, we present two sets of for . Each data point is obtained by averaging over repetitions yielding an expectation value with statistical uncertainty . We compare with numerical full ED of Eq. (1) with . As increases, the accuracy of numerical results decreases significantly. For , we have . For , since , we exclude numerical results in Figs. 2 and 3; here, even state-of-the-art full ED methods Beugeling et al. (2014) could consider only sup (). For and , we confirm oscillations of high and persisting amplitude due to the coupling to the only mode at . For increasing , the spin couples to modes including higher order processes, such that spin excitation gets distributed (entangled) into the growing bosonic environment. Hence, coherent oscillations at incommensurate frequencies lose their common contrast and appear damped. After the transient duration , with , the spin observable remains close to its time average. Still, the conservation of coherence of the evolution is evident in our measurements: Revivals of spin excitation due to the finite size of the system appear at , where is the mean energy difference between modes. And, for , negative time averages of indicate equilibration of the system to the ground state of , biased by . Both observations present strong independent evidence that our total spin-phonon system is near-perfectly isolated from external baths. Independent measurements yield a decoherence rate of sup ().

Figure 3: Time averages and mean amplitudes of time fluctuations of . These are calculated from experimental traces (black dots, error bars: 1 s.d., derived from the underlying probability distribution of sup ()) for varying and and comparison to full ED (solid lines). (a) Increasing shifts the ground state of , adjusts spin-mode couplings, and varies . Even for small systems, we find agreement of time averages with microcanonical averages, (dashed lines, shaded areas indicate systematic uncertainties). As rapidly increases with , time averages follow microcanonical averages for a larger range of . (b) gradually decreases with and correlated resonances in and fade away, indicating that we tune our system from microscopic to mesoscopic size.

This complements the agreement of experimental with numerical results, where we set .

We analyze all recorded time evolutions, each containing data points in the interval , by deriving time averages


and mean amplitudes of time fluctuations


The quantities are illustrated in two examples in Fig. 2(c). We plot these in Fig. 3 for and as a function of , together with full ED results for (solid lines). Tuning across the maximum of , cf. Fig. 1(b), and comparing to numerically calculated microcanonical averages (dashed lines) sup (), we find agreement for a larger range of when increasing . This indicates an extended regime permitting thermalization. For large , we observe its breakdown as the spin couples to an environment of decreasing complexity. Finite-size effects are prominent in resonances of and for , while their amplitudes gradually fade away for higher .

Figure 4: Scaling of mean amplitudes of time fluctuations with and . (a) We plot (error bars: 1 s.d.) as a function of . The spread for highlights finite-size effects, and we show an average value for each (large symbols, error bars: 1 s.d.). For , we observe a decay that ceases for . (b) as a function of calculated (error bars: systematic uncertainties), which captures the dependence of the effective system size on all experimental parameters. We compare to a scaling (solid line), motivated for our system, and our measurements agree for . Further increasing , the system needs longer durations to resolve decreasing energy differences in the environment, unveiling the importance of time scales.

For further analysis, we postselect data points well described by microcanonical averages, i.e., with a deviation of less than sup (). For those, we show the dependence of on in Fig. 4(a). Although sets the size of the environment, the complexity of the spin-phonon coupling is tuned by , , , , and , cf. Figs. 1(b) and 1(c). Consequently, we study the correlation between and by combining our experimental results with numerical calculations of in Fig. 4(b). In general, mean amplitudes of time fluctuations are predicted to be upper bounded by . For our system, we even find a proportionality, : We motivate this scaling, illustrated by the solid line in Fig. 4(b), by a heuristic derivation considering pure initial states and infinite times, which relies on the eigenstate thermalization hypothesis sup (). Our measurements feature such a scaling for , despite our nonidealized initial states and finite observation duration. We observe that, for , measured mean amplitudes of time fluctuations do not further decrease. We attribute this to the fact that a system of increasing complexity features decreasing energy differences in its spectrum, corresponding to smaller relevant frequencies in the dynamics. Explicitly, the system requires longer durations to approach theoretically predicted values. Here, theory considers averages for infinite time, and does not make any prediction about relevant time scales in the dynamics.

In summary, we scale our trapped-ion system including its engineered environment up to relevant Hilbert space dimensions challenging state-of-the-art full ED. We present time averages and fluctuations of a spin observable and exploit an effective dimension to study their dependence on the size of the system. We observe the emergence of quantum statistical mechanics within our isolated system despite its moderate size. Simultaneously, we monitor the coherent dynamics of thermalization, revealing the importance of initial and transient time scales by direct observation of the evolution towards thermal equilibrium. Thereby, we contribute to open questions in the field of thermalization Popescu et al. (2006); Linden et al. (2009); Eisert et al. (2015). Our approach admits generating a multitude of initial conditions, choosing different system and environment states, and preparing initial correlations Leibfried et al. (2003); Blatt and Wineland (2008); Kienzler et al. (2015). In addition, it allows us to measure a variety of observables Leibfried et al. (1996, 2003); Kienzler et al. (2015). Applying those techniques, we can study, e.g., non-Markovianity of the dynamics, which is evidenced by revivals in the evolution, in detail Breuer et al. (2009); Clos and Breuer (2012). Further, increasing the strength of the spin-phonon coupling, we can effectively expand the observable time span. Possible extensions include incorporating more and larger spins, tuning long-range interactions, adding external baths Porras and Cirac (2004); Porras et al. (2008); Myatt et al. (2000), and propelling experimental quantum simulations. Beyond numerical tractability, our experimental setup can be used as a test bed to assess the validity of approximated theoretical methods that address strong couplings to vibrational baths in a variety of fields, such as molecular and chemical physics.

Recently, we became aware of related studies with trapped ions, superconducting qubits, and ultracold atoms Smith et al. (2016); Neill et al. (2016); Kaufman et al. (2016).

We thank H.-P. Breuer for discussions, M. Enderlein, J. Pacer, and J. Harlos for assistance during the setup of the experiment, and M. Wittemer for comments on the manuscript. This work was supported by the Deutsche Forschungsgemeinschaft [SCHA 973; 91b (INST 39/828-1 and 39/901-1 FUGG)], the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007-2013, REA Grant Agreement No: PCIG14-GA-2013-630955) (D.P.), and the Freiburg Institute for Advanced Studies (FRIAS) (T.S.).


Appendix A Supplemental Material

Appendix B Trapped-ion Hamiltonian

Our system is described by the Hamiltonian presented in Eq. (1), which we recast Porras et al. (2008) from:


As any two-level system, two electronic states can be interpreted as a pseudospin-1/2, represented by the Pauli operators . Controlling the effective energy difference of the states, , via the first term of the Hamiltonian is equivalent to applying an effective magnetic field to the spin yielding a Zeeman shift that lifts the degeneracy of the system in a controlled way.

The second term describes the environment composed of harmonic oscillators, i.e., bosonic modes, with incommensurate frequencies , and their excitations (phonons) constitute the environment of the spin in the otherwise closed quantum system. In Table 1 we list descriptions of all parameters in Eq. (8).

parameter name realization physical interpretation variation in experiment
number of axial motional modes number of ions number of axial harmonic oscillators spanning the bosonic environment controlled by deterministic loading of ions,
axial mode frequency mutual Coulomb repulsion within axial trapping confinement resonance frequencies of the harmonic oscillators (phonon energy ) controlled by axial confinement, here fixed,
spin coupling rate (Rabi frequency) Two-Photon Stimulated Raman transition (TPSR) resonant coupling of the spin states and , driving coherent oscillations between them controlled by laser intensity, chosen for dedicated to maximize
effective magnetic field TPSR detuning from resonant (i.e., carrier) transition lifts degeneracy of coupled spin states and , i.e., introduces a bias and alters the total ground state controlled by variable detuning
spin-phonon coupling (Lamb-Dicke) parameters momentum kick by the two photons of the TPSR momentum conservation when creating and annihilating phonons: (photon recoil energy)/(phonon energy); for : weak coupling (only carrier transition, i.e., spin flip) controlled via (confinement), here fixed,
Table 1: Detailed description of all parameters in the Hamiltonian.

We couple spin and phonons through the photon recoil term . Here, is the position of the ion, carrying the spin, and represents the effective wave vector of the optical field, driving a two-photon stimulated Raman transition with spin coupling rate , the so-called Rabi frequency (see below). The related recoil, , is required for creating (annihilating) phonons within the modes while flipping the spin, that is, assures energy and momentum conservation of the spin-phonon coupling. We want to emphasize that, in principle, the optical drive induces coherent transitions only, and, consequently, does not represent a channel to an external bath. The last term in Eq. (8) incorporates both, the term , which does not affect the motional state, and the spin-phonon coupling , where the spin-flip operators can be related via .

The position is expressed as


with corresponding spin-phonon coupling (Lamb-Dicke) parameters


Here, denotes the mass of the ion chain and is the wavefunction amplitude of mode . Further, it is convenient to use , with describing the longitudinal center-of-mass mode ().

For values , we find linear/weak spin-phonon coupling that is commonly explored in trapped-ion experiments and


We, however, tune to exploit higher-order phonon emission/absorption terms and in our experiments the full expression,


is relevant. Driving at has two main consequences: First, the coupling of the spin states and the spin-phonon coupling get further increased and second, the individual transitions between single and multiple modes cannot get resolved anymore. Note, that the coupling strength can be effectively enhanced, as is proportional to and, thus, scales with the Fock state :


In our work, we explicitly make use of the nonlinear/strong coupling and bring our system into a nonintegrable regime. The system Hamiltonian is exactly solvable in the following cases: (i) and , in which it corresponds to the quantum Rabi model Braak (2011), and (ii) , in which spins and phonons are not coupled at all. In any other parameter regime, it can be considered a nonintegrable Hamiltonian, by the applicable definitions of nonintegrability Caux and Mossel (2011): (i) it does not have an exact analytical solution, and (ii) there is not a number of independently conserved observables that equals the number of degrees of freedom in the system.

The full quantum dynamics of the thermalization process may be treated with quasi-exact methods, since relevant aspects of thermalization, such as time scales and size of fluctuations, critically depend on the set of frequencies involved in the dynamics. In turn, standard approaches for the spin-boson model, such as path-integral methods or generalized master Eq. Weiss (1999), rely on the assumption that the phonon bath is characterized by a continuous spectrum of frequencies, and thus, they cannot be applied to our mesoscopic system. We access the dynamics by full ED, which provides all eigenstates and eigenenergies of the total (truncated) Hamiltonian, and has been done in recent theoretical studies that investigate similar physics in other systems Rigol et al. (2008); Rigol (2009); Santos and Rigol (2010); Beugeling et al. (2014); Kim et al. (2014); Cazalilla and Rigol (2010).

Appendix C Inverse Participation Ratio

The Inverse Participation Ratio (IPR) is used as an ergodicity measure to quantify the ability of a system to thermalize Srednicki (1999); Linden et al. (2009). Here, thermalization describes a process in a closed quantum system, where values of observables (of any subsystem) equal microcanonical averages Deutsch (1991); Srednicki (1994), while the total system remains in an entangled state during the unitary evolution. In general, it requires for (any) few-body observable: (i) the equivalence between time averages and microcanonical averages, and (ii) decreasing mean amplitudes of time fluctuations around a mean value as a function of the system size.

The IPR is defined for a set of energy eigenstates and an initial pure state as


Recent numerical work shows that systems with increasing values of IPR fulfill condition (i) more accurately, see for example Ref. Beugeling et al. (2014), where it is shown that deviations of diagonal ensemble predictions from microcanonical ensemble averages decrease with IPR. Intuitively, we may expect that it also leads to condition (ii): We argue that a large IPR involves a large set of eigenfrequencies in the dynamics of few-body observables and destructive interferences can lead to a suppression of the size of fluctuations on average. However, large time fluctuations can still appear, but during ever shorter durations. Indeed, an upper bound on these fluctuations can be obtained rigorously based on Ref. Linden et al. (2009). For our system, one can find a heuristic derivation that mean amplitudes of time fluctuations monotonically decrease as a function of IPR (see below). The derivation relies on the Eigenstate Thermalization Hypothesis (ETH), plausible assumptions on the many-body eigenstates, and a large number of participating states (parameter regimes with large IPR). We mathematically formulate the ETH and explain it in detail in the last section of this Supplemental Material. To examine the validity of our derivation, we perform numerical calculations and compare them to the derived scaling:


which is valid for initial pure states and infinite times. In Fig. 5 we show numerical results for for a variety of initial conditions and system parameters as a function of IPR, and the estimate based on Eq. (15) is shown as a solid line. We find reasonable agreement even for small values of IPR.

Figure 5: Numerical study of the scaling of mean amplitudes of time fluctuations. Filled and empty symbols correspond to initial states with one and two phonons per mode, respectively. Time fluctuations averaged over infinite time are calculated using Eq. (20). Parameters: , 20, 0.7 MHz. , 10, 1.0 MHz. , 6, 1.3 MHz. , 5, 1.4 MHz. , 4, 1.6 MHz. For all values of , we set MHz and take , in steps of . This numerical study substantiates the validity of our heuristic derivation of Eq. (15) (solid line), and our experimental observation of the scaling with effective dimension.

In our study, the initial state is a product state of a pure spin state and a motional state cooled close to its ground state, i.e., a mixed phonon state, yet with low average occupation. Thus, we need to find a method to average the IPR and get an appropriate measure of the effective dimension of an initial mixed state. In Ref. Short and Farrelly (2012) it is shown that, considering a mixed initial state , an upper bound to mean amplitudes of time fluctuations around time-averaged states can be found from the definition , with . However, this approach has the disadvantage that the IPR takes large values for highly mixed initial states, even in the case of uncoupled systems, since then and . We choose to define an effective dimension by averaging over the IPRs of each state in the mixture, which directly leads to our definition in Eq. (5). This definition yields for uncoupled Hamiltonians. It takes large values only if interactions lead to a large number of eigenstates participating in the dynamics, allowing us to identify regions where the system should be able to thermalize.

The effective dimension is a weighted measure of the number of coupled basis states that are required to express the uncoupled basis states. For example, consider a pure state in the uncoupled basis, represented in the coupled basis , i.e., . Typically, most coefficients , such that a significant value of is obtained only after considering a large number of states , each of them contributing with a small fraction to the final effective dimension. In addition to this, part of our our initial state is a thermal state that incorporates a large number of uncoupled states, . Our definition is chosen such that we do not artificially increase by simply adding these states, but we renormalize according to their contribution to the mixed state.

Appendix D Numerical diagonalizations and their systematic uncertainties

We perform full ED of Eq. (1) to calculate and corresponding microcanonical averages up to . For the calculations, we use Matlab on a workstation with a quad-core 3.7 GHz processor and 64 GB RAM. Besides limitations in computation time, a fundamental issue when calculating the full eigensystem of a matrix lies in the size of the available main memory that is required to store and process this matrix. Our computer is unable to fulfill this task (already) for systems of ions with a cutoff at phonons per mode that corresponds to a matrix of dimension . For our experimental parameters, this forces us to neglect more than half of the initial-state population, , while would be required for and for . As a benchmark of our numerical effort, we compare this to recent theoretical studies that diagonalize maximum sizes of:  Luitz et al. (2015),  Mondaini et al. (2016), and  Beugeling et al. (2014), while the authors state that those are the maximum sizes tractable by their methods. In Table 2 we list parameters used for all presented numerical results.

For calculations of , we use an approximation, i.e., a block-diagonalization procedure that exploits the band structure of our Hamiltonian matrix. We express Eq. (1) in the basis of eigenstates of the uncoupled system , with , and the Fock state with phonons in the th motional mode. The Hamiltonian takes the band matrix form,


with the matrix of nondiagonal elements. We order the states according to their energy such that the diagonal entries grow with index number . The matrix couples only those states that fulfill . Further, we express the initial state , calculate of each of the initial states by diagonalizing , and keep only quantum states and close to such that (, where denotes the number of states kept to calculate the IPR of each of the states . Thus, is projected onto a subspace of , that optimally contributes to of . Finally, we diagonalize , calculate the contribution of , and increase by steps of , until the resulting values of vary by less than for . For example, for , we use , corresponding to , which converges using . However, in the case of and , we are (again) at the limit of our computational power, because we need for converged values of .

Figure 6: Numerical study of as a function of truncated initial-state population. We plot for (details see text). The truncation corresponds to lowest , while is increased in steps of one up to its maximum value (limited by our computer memory; see Tab. 2). The parameters are the same as in Fig. 1(c): MHz, , , MHz, . For and for , we find that agree within their attributed systematic uncertainties (error bars) with numerically most accurate values which have . This suggests that our estimation of yields reasonable results, yet with large systematic uncertainties, even if about of the initial-state population are neglected in the numerical calculations.
Table 2: Parameters used for numerical calculations. Cutoff in Fock space and truncated initial-state population used in the numerical calculation for the different Figures.

In the following, we describe our procedure to yield final values of including a measure of systematic uncertainties accounting for population which has been neglected by the truncation, . We evaluate as a function of truncation for all parameter settings presented in our manuscript (for a list of see Tab. 2). For each parameter setting, we calculate the truncated initial-state population , take as a first bound , linearly extrapolate from it a second bound with the slope extracted from , and plot as a function of . In Fig. 6, we show an example of this. The error bars represent our attributed systematic uncertainties, given by . For and all parameter settings, we find that, for , all extracted agree within their uncertainties with corresponding , which, in turn, we consider most accurate, because for , it is ensured that . This suggests that our procedure yields reasonable values for , yet with large uncertainties, even when up to of the initial-state population is neglected in the numerical calculation, as it is the case for . We use for final values , throughout the manuscript.

To calculate the microcanonical average, , we calculate mean and standard deviation of the energy with respect to the initial state, , and


We define a microcanonical ensemble by using a Gaussian distribution of the form:




The ETH implies that is a smooth function of the energy , and, thus, the microcanonical prediction should not depend on the details of the energy distribution. Further, we apply a procedure similar to the one described above for , but based on ED (no block-diagonalization), to provide final values of including attributed systematic uncertainties. We include numerical calculations of microcanonical averages in Fig. 3 and Fig. 8, even for , where . This estimation is only used to detect regions in parameter space where the measured mean value is close to the microcanonical average (cp. Fig. 8), and is not used in our further analysis. Additionally, numerical calculations show a faster convergence of with respect to the cutoff , since that average is less sensitive to the total number of quantum many-body states than IPR and .

Appendix E Experimental setup

Figure 7: Test for unitarity of the experimentally measured time evolutions. (left) Measured values with statistical uncertainties (black dots) and full ED of the unitary time evolution, i.e., with no decoherence (solid line). The shaded area emphasizes the time span that is considered to derive and (see Figs. 2 and 3), . For , we find our experimental data in good agreement with the numerical calculations. (right) The remaining discrepancy can be captured by implying a weak decoherence indicated by the dashed line , which is in agreement with independent calibration measurements. Parameters for the experiments: (a) , MHz, (b) MHz, (c) MHz, . This comparison further justifies to describe our system by means of unitary time evolution, and, consequently, treat it as fully isolated within our experimental observation.

We employ a linear radio-frequency (rf) Paul trap with a drive frequency MHz to trap Mg isotopes with secular frequencies of MHz (radial direction) and MHz (axial direction) Schaetz et al. (2007). Using a photoionization laser (wavelength nm), we isotope-selectively load Mg and Mg, and prepare identical spatial configurations in all measurements Hume (2010). Two electronic ground states of the hyperfine manifold of Mg (nuclear spin ) constitute the pseudospin, and , where and denote the total angular momentum quantum numbers of the valence electron. Note, that Mg has no nuclear spin () and consequently, in our realization, it constitutes (only) to the motional mode structure.

To address transitions from to either or in Mg at wavelengths around nm, we use dedicated all solid state laser systems Clos et al. (2014); Friedenauer et al. (2006): A laser beam with  nm (labeled BD) is tuned below the to transition with natural line width  MHz Ansbacher et al. (1989) and aligned with a magnetic quantization field mT. The -polarized light with intensity (saturation intensity W/m) is used for Doppler cooling to and optical pumping to . This state provides efficient, state sensitive detection via the closed cycling transition to   to discriminate between the ground state manifolds. When BD is shifted near resonance (detuning  MHz), we record photon-count histograms with averages corresponding to count rates of ms for all states and ms for the states. We can detect experiments in which the ion-chain order is disturbed by a systematic change in histogram distributions, and postselect experiments where the desired ion-order is achieved, which is true in of all our experiments for . Two additional laser beams with nm (labeled RP1 and RP2) are superimposed to BD and used for repumping electronic state population from (RP1) and (RP2) to via (with MHz  Ansbacher et al. (1989)). Further, two perpendicular laser beams at  nm (labeled RR and BR) with a -vector difference, along the axial direction of the trap (ion-chain alignment) enable ground-state cooling of axial motional modes. This is achieved via sideband cooling with two-photon stimulated Raman transitions (TPSR) detuned by GHz from the to transition. In the regime of resolved sidebands, e.g., for MHz, , and MHz, we routinely achieve occupation numbers by an iterative pulsed sideband cooling procedure. For MHz, we use MHz and achieve for by addressing several sideband transitions simultaneously. In Table 3, we summarize relevant experimental parameters that are determined with dedicated calibration measurements, e.g., mode-temperature measurements Wineland et al. (1998).

Table 3: Experimental parameters. Lowest (highest) frequency axial vibrational mode , effective magnetic field varied within interval, spin coupling rate , detuning of the Raman beams from the to transition, measured mode occupation number . The detection duration is s for and s for .

All terms incorporating the spin in the Hamiltonian are implemented via TPSR as well. In order to achieve strong coupling, i.e., Rabi frequencies up to MHz, we use intensities and . We employ two acousto-optical modulators to switch on and off beams (turn coupling on/off), to tune the relative frequency difference between RR and BR, i.e., the detuning of the TPSR (varying ), and to attenuate the intensity of BR (fine tune ).

At intensities and , the contribution of spontaneous Raman scattering Ozeri et al. (2007) of Mg to residual decoherence in our experiments, is calculated to kHz, assuming a relevant beam waist ( radius of intensity) m for RR and BR. This corresponds to less than scattering events in , resulting in a residual heating effect of less than quanta (in ) distributed among motional modes. This is in agreement with dedicated calibration measurements of the decoherence rate induced by RR and BR. In Fig. 7, we compare typical experimental results for to corresponding calculated unitary time evolutions with and without considering a decoherence rate . From this we infer that such that decoherence contributes to a relative uncertainty of (in most cases much less than) to our derived quantities during . We neglect this effect throughout our manuscript and consider our system completely isolated from external baths.

Appendix F Applicability of microcanonical averages

Figure 8: Comparison of experimental results to microcanonical averages. Deviation from (with statistical uncertainties) as a function of (with systematic uncertainties) for . We observe , while ranges of experimentally accessible values of for fixed overlap, and, additionally, grow with increasing . At fixed , approaches with increasing and the points with a deviation (dashed line) are used for evaluation in Fig. 4.

In Figure 3(a), we display time averages as a function of effective magnetic field for and compare them with the microcanonical ensemble averages . We correlate the deviation