Observation of supersymmetry and its spontaneous breaking in a trapped ion quantum simulator.

Supersymmetry (SUSY) helps solve the hierarchy problem in high-energy physics and provides a natural groundwork for unifying gravity with other fundamental interactions. While being one of the most promising frameworks for theories beyond the Standard Model, its direct experimental evidence in nature still remains to be discovered. Here we report experimental realization of a supersymmetric quantum mechanics (SUSY QM) model, a reduction of the SUSY quantum field theory for studying its fundamental properties, using a trapped ion quantum simulator. We demonstrate the energy degeneracy caused by SUSY in this model and the spontaneous SUSY breaking. By a partial quantum state tomography of the spin-phonon coupled system, we explicitly measure the supercharge of the degenerate ground states, which are superpositions of the bosonic and the fermionic states. Our work demonstrates the trapped-ion quantum simulator as an economic yet powerful platform to study versatile physics in a single well-controlled system.

Introduction

Supersymmetry (SUSY) is a quantum mechanical symmetry that unifies the space-time and the internal degree of freedom of elementary particles[1,2]. In an SUSY theory, bosons and fermions have one-to-one correspondence with the same mass, which leads to a cancellation in the Lagrangian of the Higgs particle and thus helps solve the hierarchy problem or the ultraviolet divergence in Standard Model[2]. However, in reality, this cancellation in Higgs mass is not exact and no SUSY partners of known particles have been discovered at the current energy scales, which requires the spontaneous breaking of SUSY[2,3]. Since it is mathematically daunting to decide if SUSY is spontaneously broken in a quantum field theory, supersymmetric quantum mechanics (SUSY QM)[4,5] has been proposed as a toy model for understanding this crucial property[3]. In an SUSY QM model, Hamiltonians have nonnegative eigenvalues E ≥ 0 with all the E > 0 levels being degenerate, which is a consequence of the boson-fermion correspondence[6]. On the other hand, there can also be E = 0 bosonic or fermionic ground states, which are vacuum states annihilated by supercharges, the generators of SUSY, and do not necessarily have the boson-fermion correspondence. This can be characterized by the Witten index, the number of E = 0 bosonic states minus that of the fermionic states[6]. If such E = 0 levels do not exist, the SUSY is said to be spontaneously broken as the degenerate ground states with positive energies transform between bosons and fermions by the supercharge rather than being invariant[6].

Although the SUSY theory has inspired wide applications in fields like optics[7], condensed matter physics[8], quantum chaos[8] or even outside physics[9], whether it correctly describes the physical world has not yet been determined by the state-of-the-art high-energy physics experiments[10,11] or other indirect experimental evidences[12,13]. However, many theoretical proposals already exist[14-18] to examine its effects by quantum simulation[19,20]. As one of the leading platforms for quantum information processing with long coherence time, convenient initialization and readout, as well as accurate laser or microwave control[21-25], ion trap has demonstrated the quantum simulation of various phenomena such as quantum phase transitions[26,27], many-body dynamics[26], relativistic effects[28] and quantum field theories[29]. In this work, we report experimental realization of a prototypical SUSY QM model[14] in a trapped ion quantum simulator and demonstrate the spontaneous SUSY breaking in this model. The SUSY QM model is realized by tuning suitable parameters in a quantum Rabi model (QRM) Hamiltonian[27,30]. Through a combination of the state-of-the-art manipulation techniques for spin and phonon states in ion trap, and the joint spin–phonon state tomography scheme we develop that has not been achieved before, we explicitly demonstrate the characteristic signatures of the SUSY theory by measuring the energy degeneracy between the bosonic and the fermionic states, and the nonvanishing supercharges of the degenerate ground states in the spontaneous SUSY breaking case.

Results

Our experimental scheme is sketched in Fig. 1. We create a QRM Hamiltonian with bichromatic Raman laser beams on a single trapped 171Yb+ ion[27,30] (see Methods)where σ, σ denote the components of the Pauli spin operator, a (a†) the bosonic annihilation (creation) operator for the phonon mode, and the last term a renormalization constant[14]. This system possesses two supersymmetric points. One is at g = 0 and ω = ω where the system is as simple as a spin and an oscillator without coupling. This becomes evident if we write σ = 2f† f − 1 and σ = f + f† where f (f†) denotes the annihilation (creation) operator of a fermionic mode. The excited states |↑〉|n〉 and |↓〉|n + 1〉 (n ≥ 0) are degenerate fermionic and bosonic states (|↑〉 and |↓〉 have fermionic occupation number of 1 and 0, respectively) with positive energy E = (n + 1)ω, while the unique ground state |↓〉|0〉 has zero vacuum energy E0 = 0 owing to the cancellation between the spin (fermionic) energy ω and the phonon (bosonic) energy ω.

Fig. 1

Quantum Simulation of SUSY and its spontaneous breaking.

We simulate the QRM Hamiltonian by coupling the internal states of an ion with its spatial oscillation using bichromatic Raman laser beams. Under suitable parameters, the system possesses SUSY and there is a unique E = 0 bosonic ground state with all the other excited states being double-degenerate for bosonic and fermionic states. This corresponds to a Witten index of one. When moving to a different set of parameters, the Hamiltonian is still supersymmetric, but this symmetry spontaneously breaks and the ground states are no longer unique but obtains nonvanishing supercharges and positive energies. The Witten index in this case is zero.

The other more interesting SUSY point appears if we turn on the spin–phonon coupling g and set ω = 0 (see Supplementary Note 1). In this special situation, while the Hamiltonian is still supersymmetric, its ground state is not, which indicates a spontaneous SUSY breaking. We connect these two supersymmetric points by the pathwhich corresponds to ω = (1 − r)ω and g = rg in the QRM. The energy spectrum for this system under ω = g = 2π × 5.73 kHz is plotted in Fig. 2a, b with r ranging from zero to 0.8. At r = 0, the first and the second excited states are degenerate, with a finite energy gap of ω from the unique ground state. As r increases, the degeneracy between the two excited states is lifted while the gap between the ground state and the first excited state shrinks. To probe this energy spectrum, we first prepare the ground state at the desired r through an adiabatic passage; then we apply a weak probe pulse H = Ω σ cos ω t with the QRM Hamiltonian turned on and measure the change in the spin population[30] (see Methods). By scanning ω, we obtain the energy gap between the ground and the excited states as shown in Fig. 2c, d as two examples.

Fig. 2

QRM energy spectrum and SUSY.

a Theoretical energy gaps between the ground state and the first two excited states as a function of r in Eq. (2). b Experimentally measured energy gaps. At r = 0, the system is supersymmetric with a unique bosonic ground state and degenerate first and second excited states for the bosonic and the fermionic modes. As r → 1, the Hamiltonian approaches the other supersymmetric point, but the unique ground state disappears and the lowest bosonic and fermionic states become degenerate, which indicates a spontaneous SUSY breaking. c, d The resonant signals probed by a weak pulse for energy gaps at r = 0 and r = 0.3, respectively. The solid blue curves are the theoretical results and the red dots are the measured data, which are fitted by multiple Gaussian functions as the dashed orange curves to extract the peak locations.

In the above process, the energy gap between the ground state and the first excited state closes as r → 1, which seems to violate the adiabatic condition when preparing the ground state. Nevertheless, these two states locate in different symmetry branches of the QRM under parity transform , which commutes with the QRM Hamiltonian. Therefore, the two states will not evolve into each other and the adiabatic condition can still hold. However, there is another problem that as r increases, the coupling between these two lowest levels by H becomes weaker, which leads to a reduction in the resonant signal even for the ideal state as illustrated in Fig. 2a at large r. Consequently, it is still difficult to accurately determine the energy degeneracy near r = 1, the nontrivial spontaneous SUSY breaking point. For this purpose, we explicitly measure the energy expectation values of the two ground states at r = 1.

To prepare these two states, we fix ω = 0 in Eq. (1) and gradually turn up g from zero to g (see Methods). The two initial ground states |↑〉|0〉 and |↓〉|0〉 will then adiabatically evolve into the ground states at the spontaneous SUSY breaking point , which are ideally two Schrödinger’s cat states.

Now to evaluate the expectation value of the Hamiltonian , we only need to measure the average phonon number 〈a†a〉 and a spin–phonon coupling term 〈σ (a + a†)〉. The former can be derived by the standard procedure of first resetting the spin state to |↓〉 and then driving the phonon blue-sideband to fit the phonon number population from the spin dynamics[21], as shown in Fig. 3a, b. As for the second term, we apply a spin-dependent force HSDF = (−Ω/2)σ (a + a†) and measure the evolution of σ (t) by observing that has a linear term in t as Ω tσ (a + a†)[28]. Therefore, after preparing |ψ±〉, we adjust the parameters of the bichromatic Raman laser beams to create this spin-dependent force and fit 〈σ (t)〉 to extract the slope at t = 0 (see Methods), from which we obtain 〈σ (a + a†)〉, as shown in Fig. 3c. Combining these two results, we obtain E+ = 2π × (7.3 ± 1.9) kHz and E− = 2π × (5.2 ± 2.1) kHz for |ψ+〉 and |ψ−〉 respectively at the parameters ω = 2π × 10 kHz and g = 2π × 5.73 kHz. The two levels are degenerate within about one standard deviation, and also agree well with the theoretical value of E0 = ω/2 = 2π × 5 kHz. Note that the error bar here, mainly caused by fitting the phonon distribution, is much smaller than the gap ω to higher levels. Such a double-degeneracy of the ground state with positive energy clearly indicates the spontaneous SUSY breaking.

Fig. 3

Measuring ground state energies at spontaneous SUSY breaking point.

a, b Phonon number distribution of the ground states |ψ+〉 and |ψ−〉, respectively, by fitting the time evolution under a blue-sideband driving as shown in the inset. The estimated average phonon numbers are and , respectively. c Measuring by fitting the spin-up state population under a spin-dependent force. The dynamics can be fitted by an analytic formula (solid curves, see Methods) and then the slope at t = 0 can be extracted as (−25.5 ± 1.9) kHz and (−32.3 ± 1.8) kHz for |ψ+〉 and |ψ−〉, respectively. Each data point is averaged over N = 200 experimental shots. Error bars represent one standard deviation. Combining the two results, we measure the average energy of the two ground states as E+ = 2π × (7.3 ± 1.9) kHz and E− = 2π × (5.2 ± 2.1) kHz with the ideal value E0 = ω/2 = 2π × 5 kHz. The two levels are degenerate within about one standard deviation.

In quantum field theory, E0 > 0 has the nontrivial meaning of a positive vacuum energy. However, in quantum mechanics, a constant in the Hamiltonian only contributes to a global phase and can always be discarded. Therefore, we would prefer more direct evidences of the spontaneous SUSY breaking by measuring the supercharge, which is the generator of the SUSY that transforms between bosons and fermions. An SUSY-invariant vacuum will be annihilated by the supercharge, while in the case of spontaneous SUSY breaking, the degenerate ground states can take nonzero supercharges. For the Hamiltonian of Eq. (2) with r = 1, the supercharge can be defined as[14] (see Supplementary Note 1 for details)where V± ≡ exp[±(g/ω)(a† − a)] are displacement operators, and we define two phonon operators A and B to simplify the expression. We have taken out a coefficient to make the supercharge dimensionless. Ideally, |ψ±〉 are eigenstates of Q with eigenvalues of .

To evaluate 〈σ ⊗ A〉, in principle we can decompose the spin–phonon density matrix as where ρ () are operators on the phonon states. Then we have . Naively, to obtain ρ↑↑ and ρ↓↓ separately, we can first measure σ and then reset the spin state as an ancilla for quantum state tomography on the phonon state[21], conditioned on the measurement outcome σ = ±1. However, in ion trap experiments, this measurement of the spin state typically takes hundreds of microseconds to accumulate sufficient photon counts, during which the decoherence of the phonon state may not be ignored. This suggests that we shall incorporate the spin states into the tomography of the phonon states. Rather than resetting the spin state to |↓〉, we remove its off-diagonal term by applying random σ operations but keep the population in the |↑〉 and |↓〉 basis (see Methods). Then we follow the standard procedure of phonon quantum state tomography by applying displacements in different directions and blue/red-sideband driving[21,31]. Through maximum-likelihood estimation (MLE)[32], we can reconstruct the joint spin–phonon density matrix in the diagonal spin basis, that is, ρ↑↑ and ρ↓↓. Some typical fitting results for various displacement directions and blue/red sidebands are presented in Fig. 4a–d. Similarly, by first rotating the σ basis to the σ basis, we can reconstruct the partial information required to evaluate 〈σ ⊗ B〉, as shown in Fig. 4e–h. Combining the two results, we get Q− = 0.501 ± 0.016 and Q+ = −0.512 ± 0.020 after correcting a phase in the phonon state due to the slow drift of the trap frequency (see Methods and Supplementary Note 2). The measured supercharges are significantly away from zero, which suggests that the ground states are superposition of bosonic and fermionic states rather than an SUSY-invariant vacuum. Note that here the error bars only account for the statistical errors of MLE in finding the best-fitted density matrices for the prepared states (see Methods), while imperfect state preparation and measurement can still result in a systematic deviation from the ideal values of ±0.707. For example, a slow fluctuation in the trap frequency of 2π × 0.5 kHz during the joint tomography, which takes tens of minutes, can already reduce the supercharges to about ±0.5 (see Methods). Also note that, while here we are only interested in the σ and σ parts of the spin for measuring supercharges, our method can be extended to include the σ part as well. Combining these results altogether one will be able to perform full quantum state tomography of the joint spin–phonon system.

Fig. 4

Measuring ground state supercharges.

We reconstruct the phonon density matrix with the spin projected to σ = ±1 or σ = ±1 by maximum-likelihood estimation. a–d The measured spin dynamics under blue-sideband (blue dots) and red-sideband (red dots) driving and those predicted by the best-fitted density matrix (blue and red curves) for |ψ−〉 after four typical phonon displacements with σ projected to ±1. Each data point is the average of 400 experimental shots, half of which have a σ gate before the measurement sequence to remove the off-diagonal spin terms. e–h Similar plots for |ψ−〉 after four typical phonon displacements with σ projected to ±1. Here, we choose ω = 2π × 10 kHz and g = 2π × 5.43 kHz. The displacement operators D(β) are characterized by where β = 0.687, N = 12 and j = 0, 1, …, N − 1. The complete data for all the displacements and both |ψ−〉 and |ψ+〉 are presented in Supplementary Note 2. From the fitted density matrices, we get 〈σ ⊗ A〉 = 0.201 ± 0.013 and 〈σ ⊗ B〉 = 0.300 ± 0.009 for |ψ−〉 and 〈σ ⊗ A〉 = −0.163 ± 0.017 and 〈σ ⊗ B〉 = −0.349 ± 0.011 for |ψ+〉 where error bars are estimated by Monte Carlo sampling (see Methods).

To sum up, we have demonstrated the quantum simulation of an SUSY QM model using a single trapped ion. By measuring the energy spectrum of the QRM along a carefully chosen path, we illustrate the characteristic properties in the degeneracy of the energy levels, whether the SUSY spontaneously breaks or not. We adiabatically prepare the degenerate ground states at the spontaneous SUSY breaking point and explicitly measure the expectation value of the energy and the supercharge, from which the spontaneous SUSY breaking can be verified experimentally. Recently, quantum simulation using up to 16 ions and 16 phonon modes has been demonstrated[33], which should allow the demonstration of some lattice quantum field theory models using spin-boson-coupled systems, similar to ref. [29] for spin models. To follow this direction to further explore the SUSY theory, suitable SUSY lattice models need to be designed to fit into the ion trap setup, and suitable observables should be chosen to characterize its properties: while the joint quantum state tomography scheme we develop here can directly be generalized to the multi-ion case to reconstruct the joint state of any particular collective phonon mode with any particular spin, and provide some useful information, in general the full tomography of all the spins and all the phonon modes will be exponentially more challenging. Nevertheless, as we show in this experiment, partial information about the system may suffice to demonstrate the spontaneous SUSY breaking or other properties of concern. Our work showcases the suitability of ion trap as an economic but powerful platform for simulating diverse physics through its unparalleled controllability.

Methods

Experimental setup

Our experimental setup is a single 171Yb+ ion in a linear Paul trap. The spin state is encoded in the and the levels of the ion with atomic transition frequency around ω = 2π × 12.643 GHz, and we exploit a radial oscillation mode with secular frequency around ω = 2π × 2.35 MHz.

We use counter-propagating 355 nm pulsed laser beams with a specially designed repetition rate around 2π × 118.415 MHz and a bandwidth of about 200 GHz to manipulate the ion through Raman transition. Two acousto-optic modulators (AOMs) are used to fine-tune the frequency, the phase and the amplitude of the laser beams.

Before each experiment, we initialize the phonon state to |0〉 through sideband cooling using 355 nm laser. Then we initialize the spin state by optical pumping of 369 nm laser with a 935 nm laser to repump the population leaked to the 2D3/2 levels. More details about our setup can be found in our previous work[27].

Simulating quantum Rabi model

We follow the scheme of refs. [27,30] to simulate the QRM for an ion with qubit frequency ω in a harmonic trap with trap frequency ω. By applying bichromatic Raman laser beams along the x direction, we get red (blue) sideband Hamiltonian with small detuning δr = ω − ω − ω (δ = ω + ω− ω) to the red (blue) sideband, respectively. By setting δ = ω − ω, δ = ω + ω and Ω = Ω = 2g/η where is the Lamb-Dicke parameter, we get the QRM Hamiltonian of Eq. (1) in the interaction picture of H0 = (ω − ω)σ/2 + (ω − ω)a† a.

Measuring excitation spectrum

After preparing the ground state of the system (see below), we apply a weak probe with Ω = 2π × 0.17 kHz (much smaller than the typical energy scale ω and g of the system) and a duration τ = 1000 μs to measure the excitation spectrum of the system[30], whose peaks give the energy gaps to the excited states. Note that the QRM Hamiltonian is realized in an interaction picture as mentioned above, in which the effective probe frequency is ω − ω + ω. This is what we plot as the horizontal axes in Fig. 2c, d.

Adiabatic evolution

In principle, we can use arbitrary adiabatic passages to prepare the ground states so long as the adiabatic condition is satisfied. However, the change in the AC Stark shift must be taken into account when we tune the laser frequency or intensity, which may accumulate into a non-negligible phase during the slow quench.

To prepare the ground state of Eq. (2), we first set ω = (1 − r)ω and g = rg at the desired values, and choose the phonon frequency to be 10ω, which is significantly larger than the coupling g. We initialize the system to be |↓〉|0〉, which is close to the ground state in this case. Here, if we apply a weak probe to measure the excitation spectrum as mentioned above, since the spin–phonon coupling can be neglected, the resonant signal is expected to locate at ω = ω, that is, the carrier frequency. This allows us to calibrate the AC Stark shift under the coupling strength we use. Now we can adiabatically turn down the phonon frequency following an exponential function[30] (9e− + 1)ω with T = 40 μs and an overall quench time of 400 μs. From numerical simulation, the non-adiabatic excitation can be below 1% for r ranging from 0 to 0.9.

As for the degenerate ground states used for measuring the average energy and the supercharge, we observe that in the case of r = 1, the driving on the blue and the red sidebands are balanced. Our 355 nm pulsed laser has a specially designed repetition rate such that in this situation the AC Stark shift from these two sidebands largely cancel each other[33]. For coupling as large as 2π × 10 kHz, the measured AC Stark shift is still below 100 Hz. Therefore, here we simply perform a linear quench in the coupling g as g(t) = g t/τ with τ = 200 μs. From numerical simulation we see that the error due to the violation of the adiabatic condition is far below 0.1% and is negligible in our experiment.

Measuring expectation value of Hamiltonian

The Hamiltonian at the spontaneous SUSY breaking point r = 1 contains a phonon term ωa† a and a spin–phonon coupling term . The phonon term can be measured by fitting the phonon number distribution through the spin dynamics under a blue-sideband driving[21]. Details can be found in our previous work[27].

To measure the spin–phonon coupling term, we observe that has a linear term in t as , where is the evolution under a spin-dependent force. Therefore, we measure under the spin-dependent force and fit its slope at t = 0 to get . This can in principle be obtained by measuring only the initial part of the evolution, which, however, would require high measurement accuracy. Here, we observe that for the ideal ground states (), we have using the fact that displacement operator and thus . Considering experimental imperfections, we fit the measured spin-up state population as with two fitting parameters A and B, as shown in Fig. 3c. After fitting these two parameters, we can compute the slope at t = 0 and evaluate the spin–phonon coupling term using Ω = 2π × 8.1 kHz.

Partial spin–phonon state tomography

Rather than first projecting the spin state to σ = ±1 and then perform the quantum state tomography for the phonon state conditionally, here we reconstruct the the desired ρ↑↑ and ρ↓↓ altogether. Different from the standard procedure for phonon state tomography where the spin state needs to be discarded[21], here we only remove the off-diagonal terms of the spin state by inserting a σ gate on half of the measurements. Then we follow the steps of the phonon state tomography by applying phonon displacements in various directions and measure the evolution of spin population under a driving on the blue or red sidebands[21]. Theoretically, after removing the spin coherence between |↑〉 and |↓〉, the spin-up state population for any quantum state under the blue or red-sideband driving can be given byand

In principle, by fitting these curves, one can get complete information about and , the diagonal elements in the spin–phonon joint density matrix. Then by applying different phonon displacements, the desired density matrices ρ↑↑ and ρ↓↓ can be obtained[21].

In practice, this procedure of first extracting the coefficients and and then solving the original density matrix can be sensitive to noise in the measurement, and the obtained matrix may not satisfy the physical constraints such as being positive semidefinite with a trace of one. Therefore, here we use maximum-likelihood estimation[32] to find the most probable density matrix that can generate all the measured blue-/red-sideband evolutions. We use the Hamiltonian and to drive the blue or red sidebands along with a phonon displacement β at the same time[31]. We use N = 12 displacements D(β) with , β = 0.687 and , and we set a phonon cutoff of ncut = 7 for the reconstructed density matrices with only 10−9 population outside the truncated Hilbert space for the ideal states.

The projection to σ = ±1 can be obtained in a similar way by first rotating the σ basis to the σ basis. The complete data and fitting results are presented as Supplementary Figs. 1–4 in Supplementary Note 2.

Calibrate phonon state rotation error due to trap frequency drift

In Supplementary Fig. 5 of Supplementary Note 2, we present the fidelity between the reconstructed density matrices and the ideal ones (with the off-diagonal spin terms discarded as described in the main text), versus a rotation angle θ applied on the phonon state . As we can see, while |ψ−〉 has the highest fidelity near θ = 0, the best result for |ψ+〉 is achieved at a finite rotation angle of about θ = 1.74π (or θ = −0.26π). This may come from the slow drift in the trap frequency of our setup between the measurements of |ψ−〉 and |ψ+〉, which in turn would result in a drift in the rotating speed of the interaction picture and thus cause an effective rotation in the phonon state. As we present in Supplementary Fig. 6 of Supplementary Note 2, the supercharges will also be changing under the phonon rotation. The maximizer of the supercharge (strictly speaking, its magnitude) is close to that of the fidelity, but a small discrepancy exists. Here we use the θ for the highest average state fidelity to calibrate the phonon rotation in |ψ+〉, with no such corrections in |ψ−〉 since it is measured earlier and the drift in parameters can be smaller. After this correction, the measured supercharges in the main text can be obtained. Also note that once the optimal θ for |ψ+〉 is determined, it is fixed in the data processing. That is, we keep using the same θ when estimating the error using Monte Carlo sampling as shown in the next section, rather than optimizing a θ for each simulated result.

Error estimation for maximum-likelihood method

We use Monte Carlo sampling to estimate the error in the reconstructed density matrix obtained from maximum-likelihood estimation and that in the computed supercharge. These measurement outcomes are computed from the raw experimental data whether the spin is in |↑〉 or |↓〉 for each displacement and the red/blue-sideband driving at each time point, which is repeated for 400 times (200 times with σ operations and 200 times without to remove the spin coherence). Note that for each data point, the number of spin-up events follows a binomial distribution whose parameter p can be estimated using the measured frequency P↑. Therefore, we can generate new sets of measurement outcomes following these binomial distributions, and then use them to find the most likely density matrix and the corresponding supercharges. We repeat this procedure for 100 times and use the standard deviation of the simulated results to estimate the error of the measurement outcomes.

Estimation of preparation and measurement errors for ground state supercharges

The error bars for the measured supercharges using Monte Carlo sampling only consider the uncertainty from finding the best fit for the experimentally prepared ground states. If the state preparation and measurement by themselves are imperfect, clearly there will be a deviation from the ideal values of . For the adiabatic state preparation, we can numerically integrate the master equation[34] for the linear quench, using a Lindblad term describing the motional dephasing with an empirical coherence time of about τ = 1.2 ms. The supercharges will reduce to ±0.66 due to this error.

A more severe error source is the slow fluctuation of the trap frequency during the joint spin–phonon state tomography. If such a fluctuation Δω occurs when measuring the blue/red-sideband dynamics for the N = 12 displacement directions, effectively the phonon state will experience a random rotation of Δω τ where τ = 200 μs is the adiabatic evolution time. Such a random rotation among different displacement directions for the same ground state |ψ−〉 or |ψ+〉 cannot be corrected by a constant rotation angle θ as done in the previous sections for a fixed drift between the measurements of |ψ−〉 and |ψ+〉. To estimate its effect, we purposely apply a rotation Δω τ to the ideal ground state and calculate the supercharges. It turns out that a fluctuation of 2π × 0.5 kHz already reduces the supercharges to about ±0.52. In principle, such a slow fluctuation can be suppressed if we calibrate the trap frequency before measuring each data point, which however is too time-consuming for this experiment.