Author Correction: Quantum non-demolition measurement of a many-body Hamiltonian

A Correction to this paper has been published: https://doi.org/10.1038/s41467-021-22105-3

R ecent experimental advances provide intriguing opportunities in the preparation, manipulation, and measurement of the quantum state of engineered complex many-body systems. This includes the ability to address individual sites of lattice systems enabling single-shot read-out of single-particle observables, as demonstrated by the quantum gas microscope for atoms in optical lattices 1,2 , single-spin or qubit read-out of trapped ions [3][4][5][6][7] and Rydberg tweezers arrays [8][9][10][11][12] , and superconducting qubits 13,14 . In contrast, we are interested below in developing single-shot measurements of many-body observables such as the HamiltonianĤ of an interacting many-body system. For an isolated quantum system,Ĥ represents a quantum nondemolition (QND) observable, and our goal is to implement a QND measurement of energy of a quantum many-body system in an analog simulator setting. We note that quantum optics provides with several examples of QND measurements; however these have so far been confined to observables representing few quantum degrees of freedom [14][15][16][17][18][19][20] .
Developing QND measurement of a many-body Hamiltonian H provides us with the unique opportunity to distill-in a single run of the experiment-an energy eigenstate ' j i from an initial, possibly mixed, or finite temperature state, by observing in particular run the energy eigenvalue E ℓ . In case of measurement with finite resolution, this will prepare states in a narrow energy window, reminiscent of a microcanonical ensemble. We emphasize that state preparation by measurement is intrinsically probabilistic, that is, will vary from shot to shot, reflecting the population distribution. Furthermore, this provides us with a tool to determine populations and population distributions of (excited) energy eigenstates, as required in, for example, many-body spectroscopy 21 . The ability to prepare and measure (single) energy eigenstates provides us with a unique tool to address experimentally fundamental problems in quantum statistical physics, such as the eigenstate thermalization hypothesis (ETH) [22][23][24] , which asserts that single energy eigenstates of an (isolated) ergodic system encode thermodynamic equilibrium properties. Developing the capability to turn QND measurements on and off allows one to alternate between time periods of free evolution of the unobserved many-body quantum system and energy measurement. This allows quantum feedback in a many-body system conditional to measurement outcomes, and in particular provides a framework to monitor non-equilibrium dynamics and processes in quantum thermodynamics 25 , including measurement of work functions and quantum fluctuation relations (QFRs) [26][27][28] . These relations express fundamental constraints on, for example, the work performed on a quantum system in an arbitrary non-equilibrium process, imposed by the universal canonical form of thermal states and the principle of microreversibility.
Our aim below is to develop QND measurement ofĤ in physical settings of analog quantum simulation, in particular exploring the regime of mesoscopic system sizes. We will demonstrate this in detail as an example of an analog trapped-ion quantum simulator, realizing a long-range transverse Ising Hamiltonian and the associated QND measurement. Our implementation in an analog quantum device should be contrasted to QND measurement ofĤ via a phase estimation algorithm 29 , which however requires a universal (digital) quantum computer.
Results QND measurement ofĤ. On a more formal level, we define QND measurement of a many-body HamiltonianĤ as an indirect measurement by coupling the system of interest S, illustrated in Fig. 1a, to an ancillary system M as meter. In a first step, the system is entangled with the meter according to the time evolution UðtÞ ¼ expðÀiĤ QND tÞ generated by the QND Hamiltonian with coupling strength ϑ (ℏ = 1). To be specific and in light of examples below, we consider here as meter a continuous variable system with a pair of conjugated quadraturesX andP obeying the canonical commutation relation ½X;P ¼ i. Consider now an initial state of the joint system prepared as ; and x 0 j i is an (improper) eigenstate ofX (or squeezed state). We obtain for the time-evolved state Reading the meter as x ℓ ≡ x 0 + ϑE ℓ t, and thus measuring the eigenvalue E ℓ , will prepare the system in ' j i (or in the relevant subspace in case of degeneracies). The probability for obtaining the particular measurement outcome E ℓ is P ℓ = |c ℓ | 2 . Repeating the QND measurement will reproduce the particular E ℓ with certainty, with the system remaining in ' j i. The above discussion is readily extended to mixed initial system states, and to initial meter states, for example, as coherent states.
In an analog quantum simulator setting, QND measurement of the many-body HamiltonianĤ is incorporated by engineering the extended system-meter HamiltonianĤ SM ¼Ĥ Iþ ϑĤ P. In an interaction picture with respect toĤ I, the joint system then evolves according to the HamiltonianĤ int SM H QND realizing the QND measurement discussed above and illustrated in Fig. 1b. On the other hand, by allowing the systemmeter coupling ϑ(t) to be switched on and off in time, we can alternate between the conventional free-evolution simulation and QND measurement mode of the system. In an actual implementation, as discussed below for trapped ions, we will achieve building the extended system-meter Hamiltonian whereĤ 0 andĤ may differ (slightly). We note that the QND measurement ofĤ is obtained by fine-tuningĤ 0 ¼Ĥ. A mismatchĤ 0 ≠Ĥ will be visible as quantum jumps between energy eigenstates in repeated measurements.
Remarkably, in our implementation, the HamiltonianĤ 0 will differ fromĤ just by the transverse field taking on the value h 0 . We will be able to tune h ¼ h 0 , thus achieving the QND condition. As the last step in our formal development, we wish to formulate QND measurement ofĤ as measurement continuous in time [34][35][36][37] . Physically, this amounts to making frequent and, in a continuum limit, continuous readouts X(t) of the meter variableX, with the quantum many-body system evolving according to (2). Following a well-established formalism of quantum optics 38,39 , we write for the system under continuous observation a stochastic master equation (SME) for a conditional density matrixρ c ðtÞ of the many-body system. In our context this SME reads with dW(t) a Wiener increment, to be interpreted as an Itô stochastic differential equation. In a quantum optical setting, as in the ion trap example below, I(t) is identified with photocurrent in homodyne detection of scattered light 38 . Monitoring the photocurrent IðtÞ $ hĤi c thus provides continuous read-out of the many-body HamiltonianĤ with h i c Tr½ ρ c ðtÞ up to shot noise. Thus,ρ c ðtÞ describes the many-body quantum state conditional to observing a particular photocurrent trajectory I (t), as can be observed in a single run of an experiment. In (4) and (5) γ is an effective measurement rate, and ϵ is a measurement efficiency. Furthermore, we have defined a Lindblad superoperator D½ŝρ c ŝρ cŝ y À ðŝ yŝρ c þ H:c:Þ=2 describing decoherence due to the quantum measurement backaction, and the nonlinear superoperator H½ŝρ c ðŝ À hŝi c Þρ c þ H:c:, which updates the density matrix conditioned on the observation of the homodyne photocurrent. Finally, not reading the meter, that is, averaging overall measurement outcomes I(t), the SME (4) reduces to a master equation with Lindblad term $D½Ĥρ, that is, realizing a reservoir coupling with jump operatorĤ, which erases all off-diagonal terms of the averaged density matrixρ in the energy eigenbasis.
Equations (4) and (5) allow us to simulate single measurement runs corresponding to a stochastic trajectory I(t). Figure 1c illustrates ideal QND measurement,Ĥ 0 ¼Ĥ, of the Hamiltonian (3) by plotting a sample trajectory of a filtered photocurrent, obtained by averaging I(t) over a time window τ, I τ ðtÞ ¼ ð2N ffiffiffiffi ffi γϵ p τÞ À1 R 1 0 dt 0 Iðt À t 0 Þe Àt 0 =τ . As initial condition we take all spins pointing against the transverse field. As seen in Fig. 1c the trajectory I τ ðtÞ (red curve) stabilizes on a time scale~γ −1 on a particular energy eigenvalue E ℓ of (3) (up to fluctuations from shot noise). In this figure we consider and show only the eigenstates and eigenenergies (thin horizontal lines) within the symmetry sector containing the ground state of the Ising model with J, h > 0, see Methods. The collapse, and thus preparation of the many-body wavefunction in the corresponding energy eigenstate, is indicated by plotting the populations P ' ðtÞ ' h jρ c ðtÞ ' j i (blue shadings in Fig. 1c). Figure 1d shows quantum jumps between energy eigenstates induced byĤ 0 ≠Ĥ. For weak perturbation (j½Ĥ 0 ;Ĥj ( jĤj 2 ) there are rare jumps between the energy eigenstates, indicated as t 1 and t 2 for the trajectory in Fig. 1d. Finally, Fig. 1e plots the integrated current I τ ¼ ð2N ffiffiffiffi ffi γϵ p τÞ À1 R τ 0 IðtÞdt and its fluctuations as a function of total integration time τ. For N = 8 spins starting in a thermal state, the integrated current I τ (red curve) exhibits a collapse at a rate~γ to a particular energy eigenstate. The insets shows the probabilities P ℓ for various times, and the narrowing of the energy resolution as ΔE=J $ 1= ffiffiffiffiffiffiffi γϵτ p with growing τ (see Methods); first to small energy window containing a few eigenstates as in a microcanonical ensemble, and eventually to a single energy eigenstate.
Implementation with trapped ions. We now provide a trappedion implementation of the system-meter HamiltonianĤ SM (2). As shown in Fig. 2a, we consider a string of N ions in a linear Paul These two-level ions can be driven by laser light # j i ! " j i, where the recoil associated with absorption and emission of photons provides a coupling to vibrational eigenmodes of the ion chain. This includes in particular the center-of-mass motion (COM) withX andP position and momentum operators, respectively, which play the role of meter variables.
To generate inĤ SM both the Ising interaction À P i<j J ijσ x iσ x j I, as well as the Ising term coupled to COM, Àϑ P i<j J ijσ x iσ x j P, we choose a laser configuration consisting of two pairs of counterpropagating laser beams (cf. Fig. 2a). In generalization of refs. 40 1 QND measurement of a many-body HamiltonianĤ in a quantum simulator setting. The many-body spin system S, shown in a, is entangled with an ancillary system M (meter) by the unitaryÛ SM ðtÞ ¼ expfÀi R t 0 dt 0 ϑðt 0 Þ prepares the many-body system in an energy eigenstate ' j i with the eigenvalue E ℓ . c Single trajectory simulation (4) and (5) of an ideal QND measurement for the Ising Hamiltonian (3) for N = 5 spins, α = 1.5, h∕J = 1.5. The window-filtered homodyne current I τ ðtÞ (red curve) fluctuates around a value corresponding to the eigenenergy prepared by the measurement ofĤ. The thin horizontal lines show the system eigenenergies E ℓ and the blue color indicates the conditional populations P ℓ (t) of the corresponding eigenstates. d Observation of quantum jumps due to the mismatch of the transverse fieldŝ H 0 ¼Ĥ þ δh P jσ z j with δh=J ¼ À0:2. The filtered photocurrent (red) clearly shows sudden jumps between eigenstates at times t 1 and t 2 . e Preparation of energy eigenstates or microcanonical ensembles by the ideal QND measurement for N = 8 spins, α = 1.5, h/J = 0.8. The estimate of the system energy given by the cumulative time average of the homodyne current I τ (red line) gradually converges to a single eigenenergy (gray lines) as averaging time τ increases. The corresponding uncertainty (red area) due to shot noise decreases as $ 1= ffiffiffiffiffiffiffi γϵτ p  is detuned by ±Δ from ionic resonance, while the second pair (blue) is detuned by ±Δ 0 . Furthermore, we choose Δ 0 À Δ ¼ ω 0 with ω 0 the COM frequency. These four laser beams give rise to laser-induced two-photon processes involving pairs of ions, which are depicted in Fig. 2b, c). First, as shown in Fig. 2b, absorption of a photon from the one of the amber MS laser beam followed by an absorption from the counterpropagating amber beam gives rise to a two-photon excitation ## j i ! "" j i, which is resonant with twice the (bare) ionic transition frequency of the two-level ion. We emphasize that this process leaves the motional state of the ion chain unchanged, as illustrated by n j i ! n j i for the COM mode, with n the phonon occupation number. This process will thus contribute a term $σ þ iσ þ j to the effective spin-spin interaction. The second pair of MS beams (blue) will again contribute a resonant two-photon excitation, which adds coherently to the first contribution. By considering all possible processes, we obtain the effective Ising interaction À P i<j J ijσ x iσ x j I inĤ SM . An explicit expression for J ij is given in Methods in second-order perturbation theory in the Lamb-Dicke parameter η ¼ k= ffiffiffiffiffiffiffiffiffiffiffi 2mω 0 p ( 1, where m is the ion mass and k is the magnitude of the laser wavevector. Second, with the choice Δ 0 À Δ ¼ ω 0 two-photon processes involving absorption from an amber MS beam and a blue MS beam will be detuned by the COM frequency from two-photon resonance, that is, be resonant with the motional sidebands ±ω 0 (cf. Fig. 2c). These processes will change the phonon number by one, and by considering all possible processes contribute a term Àϑ , and J ij is identical to the couplings obtained above. We note that this term is of order η 3 (for details see Methods).
By considering a (small) imbalance of Rabi frequencies in MS laser configurations, we can create inĤ SM a transverse-field term Àh P jσ z j I, and in addition a term þϑh P jσ z j P (see Methods and Supplementary Note I). Thus, our laser configuration generatesĤ andĤ 0 with the same Ising term, but opposite transverse field ±h. To rectify the transverse-field mismatch, we can offset the detuning of the four lasers by a small amount ±Δ ð0Þ ! ±Δ ð0Þ À 2B. We obtainĤ as in Eq. (3) and The choice B = 2h thus allows us to tune to the QND sweetspotĤ 0 ¼Ĥ as in Fig. 1c, e, while away from this point we obtain H 0 ≠Ĥ as considered in Fig. 1d.
Finally, the homodyne current (5) corresponding to a continuous measurement of the COM quadratureX, and thus of the HamiltonianĤ, can be measured via homodyne detection of the scattered light from an ancillary ion driven by a laser on the red motional COM sideband (cf. Fig. 2a and Methods).
QND measurement protocols. Implementation ofĤ SM with time-dependent system-meter coupling ϑ(t) allows protocols where we switch between time windows of unobserved quantum simulation, and measurement of energy, and thus preparation of energy eigenstates, which is verified by observing convergence of the filtered photocurrent. In addition, the Hamiltonians (3) and (6) can be made time dependent, for example, with a timedependent magnetic field. This allows us to perform work on the system, and measure work distribution functions via measurement of energy 25 . Our QND toolbox thus opens up the door to address experimentally fundamental problems of (non-equilibrium) statistical mechanics in analog quantum simulation. We apply the QND toolbox below first to ETH 42,43 and then we consider testing QFRs 25 in interacting many-body systems. We emphasize that our setting explores naturally the interesting regime of mesoscopic particle numbers from a few to tens of spins.
Thermal properties of energy eigenstates. Single energy eigenstates ' j i can encode thermal properties, which we typically associate with a microcanonical or canonical ensemble describing systems in thermodynamic equilibrium. This eigenstate thermalization concerns, on the one hand, expectation values of few-body observables, leading to the remarkable prediction of the ETH that diagonal matrix elements 'jÔj ' have to agree with the microcanonical average at energy E ℓ , 'jÔj' E ' is the microcanonical density operator as a mixture of energy eigenstates within a narrow range centered around E ℓ . On the other hand, ETH imposes constraints on dynamical properties for diagonal and off-diagonal matrix elements ' 0 jÔj' ; for example, two-time correlation functions and dynamical susceptibilities have to be related by the fluctuation-dissipation theorem 42 . To be more specific, ETH suggests a structure 44 ' 0 jÔj' ¼ SðEÞ is the thermodynamic entropy at the mean energy E, and R ' 0 ' is a random number with zero mean and unit variance. An experimental test of ETH, therefore, requires the ability to measure both diagonal and offdiagonal elements, something that is provided by our ion toolbox. The transverse Ising model (3), as realized with ions, provides a rich testbed for ETH 45 . For 1 < α ≤ 2, this model features a ferromagnetic transition at finite temperature or energy density, in a canonical or microcanonical description, respectively. As illustrated above in Fig. 1, our trapped-ion QND toolbox enables the preparation of microcanonical ensembles of variable width ΔE. According to ETH, the ferromagnetic transition persists even in the limit of vanishing ΔE, which corresponds to the preparation of a single energy eigenstate.
For reference, the microcanonical phase diagram at finite ΔE is shown in Fig. 3a for an experimentally accessible system size of N = 14 spins and α = 1.5 (see Supplementary Note I for experimental parameters, we use QuSpin for the exact diagonalization 46 ). The ferromagnetic transition is clearly manifest in the distribution P(m x ) of the magnetizationm x ¼ N À1 P jσ x j , which is bi-modal in the ferromagnetic phase, see the inset in Fig. 3a. Consequently, fluctuationsm 2 x are finite (vanish) in the ferromagnetic (paramagnetic) phase, and indicate order even in the absence of symmetry-breaking fields. A trapped-ion quantum simulator provides the ability to perform single-site resolved read-out of spins, thus giving direct access to the distribution P(m x ) and, consequently, the fluctuationsm 2 x . Due to the quasidiagonal structure of ' 0 jÔj' for ETH-satisfying observablesÔ, the hypothesis is expected to hold for any power of such observables and, in particular, also for the full probability distribution function P(m x ) 44 . Indeed, as we show in Fig. 3b-d, respectively, we find clear signatures of the transition in P(m x ) for individual energy eigenstates both for N = 14 and even much smaller system of only N = 5 spins, in which single eigenenergies can be resolved with current experimental technology (see Supplementary Note I).
The observation of the Ising transition in single eigenstates gives a qualitative indication of eigenstate thermalization in diagonal matrix elements. A stringent quantitative assessment requires to show that fluctuations of single-eigenstate expectation values 'jÔj' around the microcanonical average trðÔρ mc E ' Þ are suppressed with increasing system size 42 . We discuss experimental requirements for such a test, along with a protocol to measure off-diagonal matrix elements, in Supplementary Note III.
Work distribution function and QFRs. Projective measurements in the energy eigenbasis are the key ingredient for the long-sought experimental verification of QFRs 25 . The challenging requirement to measure changes in the energy of the system on the level of single energy eigenstates has been achieved only recently in single-particle systems 47,48 . Our QND measurement scheme opens up the possibility to probe QFRs in a true many-body setting.
As an illustration we consider the celebrated Jarzynski equality, which describes the mean value of the exponentiated work performed on a system in an arbitrary non-equilibrium process defined by a time-dependent HamiltonianĤðtÞ 25 . The equality relates the work to the difference between the free energies ΔF of equilibrium systems described by the Hamiltonian at the initial t 0 and the final t 1 times: Here β is the inverse temperature specifying an initial canonical thermal state of the systemρ β ¼ e ÀβĤðt 0 Þ =Z t 0 with Z t 0 ¼ The average on the left-hand side of Eq. (7) is performed with respect to the distribution of work PðWÞ (see Methods). The work itself is determined as the difference between the outcomes of two energy measurements before and after the time-dependent protocol . Remarkably, while the work distribution does depend on details of the time evolution given byĤðtÞ, the average is defined only by the initial and final Hamiltonians. Therefore, the QFR enables experimental measurements of the equilibrium property ΔF via measurement of work in a non-equilibrium process.
Our scheme provides the required ingredients for probing the QFR in the interesting regime of intermediate system sizes, which is dominated by quantum fluctuations: as presented above, the scheme allows for independent temporal control of parameters of the spin Hamiltonian as well as the system-meter coupling in Eq.
(2). Further, single energy levels are well resolved for a system of five interacting spins as we consider in the following (see Supplementary Note I for experimental parameters). This enables the protocol shown in Fig. 4a-c, in which a first measurement of the energy is carried out while the magnetic field h in Eq. (3) is kept constant at a value h t 0 . Then, the system-meter coupling is switched off, and the magnetic field is linearly ramped to a final value h t 1 , followed by another measurement. The statistics of corresponding measurement outcomes, E t 0 ' and E t 1 ' 0 , determines the distribution of work PðWÞ performed on the system during the magnetic field ramp. Initialization of the system at arbitrary temperatures can be emulated by weighting different runs of the protocol according to the Gibbs distribution e ÀβE t 0 ' =Z t 0 with the initial energy E t 0 ' . The resulting work distributions PðWÞ for various quench durations t Q are shown in Fig. 4d-f. While the probability distribution for fast (blue dots) and slow (green dots) quenches differs significantly the estimated free energy difference ΔF est ¼ Àlnhe ÀβW i=β approximately (due to the finite number of simulated experimental runs) matches the true value of ΔF for all three quench speeds; thus, verifying the equality (7).

Discussion
We have developed a QND toolbox in analog quantum simulation realizing single-shot measurement of the energy of an isolated quantum many-body system, as a key element towards experimental studies in non-equilibrium quantum statistical mechanics. This comprises ETH and quantum thermodynamics, including quantum work distribution, and Jarzynski and Crooks fluctuations relations 25 in mesoscopic quantum many-body systems. The present work outlines an ion-trap implementation with COM phonons as meter. However, the concepts and techniques carry over to other platforms including CQED with atoms 49 and superconducting qubits 50 , where the role of the meter can be represented by cavity photons read with homodyne detection, and Rydberg tweezer arrays 8-12 by coupling to a small atomic ensemble encoding the continuous meter variables 51 , respectively. Finally, while the present work considers QND measurement of the total HamiltonianĤ of an isolated system, our approach generalizes to measuring HamiltoniansĤ A of subsystems, as is of interested in quantum transport of energy, or energy exchange in coupling the many-body system of interest to a bath.

Methods
System-meter coupling Hamiltonian. We choose for the four lasers in our double MS configuration the detuning and the Rabi frequency as (Δ, Ω), (− Δ, Ω + δΩ), Ω is a small imbalance we use to generate the transverse-field term in the spin model. We are interested in the regime of sufficiently large detunings compared to the Rabi frequency Ω, such that single lasers only virtually excite the ions and the phonon modes, Ω ( Δ ð0Þ , η q Ω ( jΔ ð0Þ À ω q j, where ω q is the oscillation frequency of the qth phonon mode and η q η ffiffiffiffiffiffiffiffiffiffiffiffi ffi ω 0 =ω q q . On large timescales t ) 1=Δ ð0Þ ; 1=ω 0 , we obtain an effective HamiltonianĤ SM describing the coupled dynamics of the system and the meter, that is, the spins and the COM phonon mode, by performing the Magnus expansion 52 to the time evolution operator in the interaction picture (see Supplementary Note I). We further expandĤ SM in terms of η. In second order in η we recover the transverse-field Ising Hamiltonian [30][31][32] , x iσ x j þ h P jσ z j Þ I Ĥ 0 I, where the spin-spin couplings include contributions from the two MS configurations independently with M iq denoting the distribution matrix element of the qth phonon mode. The transverse- Crucially, under the condition Δ 0 À Δ ¼ ω 0 , the crosstalk between the two MS configurations leads to an extra resonant processes as exemplified by Fig. 2c. These are described by expanding the effective HamiltonianĤ SM to third order in η, Our double MS configuration can be implemented with both axial and transverse phonon modes. The implementation with transverse modes gives rise to power-law spin interactions J ij = J ∕ |i − j| α with 0 ≤ α ≤ 3, which are considered in the rest of this paper. Experimental considerations and scalability are discussed in the Supplementary Note I.
Continuous read-out ofX. We assume that in Fig. 2 the ancillary ion does not see the four MS lasers (amber and blue) and, similarly, the system ions do not couple to the read-out laser (red), that is, we assume single-ion addressability 53 or with mixed species 5 . The read-out laser is tuned in resonance with the red sideband of the COM mode, Δ e = ω 0 , under the resolved-sideband condition ω 0 ≫ Γ e , Ω 0 ,  Fig. 4 Verification of Jarzynski equality (7) for the transverse-field Ising model with five spins. A single realization of the proposed protocol consists of a the projection of the initial thermal state of spins to an eigenstate E t 0 ' of the initial Hamiltonian with transverse-field value h t 0 , b the free evolution under a linear change of the transverse field h during time t Q = t 1 − t 0 (probabilities P ℓ (t) to populate the instantaneous energy eigenstates are shown in shades of blue), c The final energy read-out E t 1 ' of the system at transverse-field h t 1 (photocurrent realizations for simulated experimental runs are shown in red), resulting in the work W (see text) performed during the non-equilibrium process. d-f The resulting work probability distribution PðWÞ for h t 0 ¼ 2J, h t 1 ¼ 0:5J, α = 2, β = 0.5∕J, and various quench times t Q . The gray bars show theoretical probability as a function of work W. The color dots with the corresponding vertical error bars (one standard deviation) show the estimated probabilities for 1000 simulated experimental runs. Independent of the quench duration the Jarzynski relation (7) yields similar estimations of the free energy difference ΔF est (fluctuating due to the finite number of runs) with the true value given by ΔF ≈ 4.35.
where Γ e is the spontaneous emission rate of the cooling transition e j i ! g j i, while Ω 0 and Δ e are the Rabi frequency and the detuning of the cooling laser, respectively. In this regime, the emitted electric field is proportional to hâ 0 i withâ 0 the annihilation operator of the COM mode. Homodyne detection then directly reveals the quadrature of the COM phonon (the meter). The homodyne current can be written as (see Supplementary Note I) where ϵ is the photon detection efficiency, γ s ' k 2 0 Ω 2 0 =ð2Γ e Nm 0 ω 0 Þ is the measurement rate with k 0 the cooling laser wavevector and m 0 the ancillary ion mass, and we have chosen the frequency and phase of the homodyne local oscillator to maximize the homodyne current. Correspondingly, the evolution of the conditional state ρ SM c ðtÞ of spin system plus the meter is described by a SME dρ SM Eliminating the meter under the condition γ s ≫ |ϑJ|, we realize continuous QND read-out of the spin Hamiltonian as described by Eqs. (4) and (5) with γ = 2(ϑJ) 2 ∕γ s . We further emphasize that the read-out laser, which is tuned to the red sideband, also acts as cooling of the COM mode. Furthermore, the read-out signal can be enhanced with several ancilla ions.
Energy measurement resolution. Here we estimate the signal-to-noise ratio (SNR), which allows us to distinguish two adjacent energy levels separated by ΔE. The difference of photocurrents (5)  Considering the shot noises W 1,2 (t) of two measurements as uncorrelated and using the Wiener increment property dW 2 1;2 ðtÞ ¼ dt we obtain SNR = 2γϵ(ΔE/J) 2 τ. For a given averaging time τ, the condition SNR ≫1 provides us with the minimal energy difference we can distinguish ΔE=J ) 1= ffiffiffiffiffiffiffiffiffi 2γϵτ p .
Symmetries of the long-range transverse-field Ising model. The transversefield Ising model (3) is invariant under the reflection and spin inversion symmetry transformations. We now provide an operational definition of these symmetries and the corresponding symmetry sectors. Consider a product state vector in the σ x basis ϕ j i ¼ s x 1 ; ; s x N . The reflection operator can be defined by its action on the ϕ j i state as R s x 1 ; ; s x N s x N ; ; s x 1 . Analogously, the spin inversion operator can be defined as P s x 1 ; ; s x N Às x 1 ; ; Às x N . Both operators have two eigenvalues ± and commute with each other and the Hamiltonian Eq. (3), thus representing QND observables, which can also be measured in the non-destructive way as presented in the paper.
The Hamiltonian can be independently diagonalized in each of the subspaces corresponding to eigenvalues of the R and P operators. The ground state of the Ising model with J, h > 0 belongs to the þ1; þ1 f gsymmetry sector. For the test of ETH in Fig. 3 we consider the symmetry sector with eigenvalues of R and P given by þ1; À1 f g , respectively. The subspace can be reached from the þ1; þ1 f gsector by flipping odd number of spins (along the σ z direction) in the limit of strong transverse field.
Interaction renormalization. In numerical simulations in Fig. 3, we renormalize the interaction strength coefficient J such that the average interaction strength matches its value in thermodynamic limit. More precisely, for the N-spin Ising model (3) we rescale J → J N ≡ J ⋅ S N ∕ S ∞ , with S N 1 N P N i;j¼1 1= i À j j j α . The results are then expressed in units of J 14 .
Work distribution function. The work distribution of a process defined by a timedependent HamiltonianĤðtÞ (with the corresponding instantaneous energy eigenvalues and eigenstates E t ' and eigenstates ψ t ' ) is defined as follows 25 : where P 0 is the occupation probabilities of the initial state and P t ' 0 ' ¼ ψ t ' 0 Uðt; 0Þ j j ψ 0 ' 2 is the transition probabilities between initial ℓ and final ' 0 states with Uðt; 0Þ T e Ài R t 0Ĥ ðt 0 Þdt 0 the evolution operator. The average of the exponentiated work in Eq. (7) is readily defined as an integral with the work distribution function Eq. (11): he ÀβW i R dWe ÀβW PðWÞ.

Data availability
The data sets generated and analyzed in the current study are available from the corresponding author upon reasonable request.

Code availability
The codes used to generate the numerical data in the current study are available from the corresponding author upon reasonable request.