Sideband thermometry of ion crystals

,


I. INTRODUCTION
Trapped ions are one of the leading platforms for quantum computation [1][2][3], simulation [4][5][6][7], sensing [8,9], and metrology [10][11][12].The exquisite degree of quantum control over electronic and motional degrees of freedom and their unrivalled coherence times make trapped ions excellent qubits [1,5] and enable quantum gates with record fidelity [13][14][15].The unique combination of isolation and controllability also guarantees low systematic effects in precision measurements of the electronic structure of trapped ions, making them perfect systems for frequency and time references in next-generation optical clocks and for searches of physics beyond the Standard Model [16].
Scaling up ion Coulomb crystals (ICCs) is desirable for all these applications, but poses increasing challenges in entropy management as the number of degrees of freedom of motion grows.This concerns first the efficient cooling of ICCs, since thermal excitations are a major limiting factor, e.g. in quantum gates mediated by normal modes of motion and, due to the relativistic Doppler effect, also in the systematics of ion clocks [17,18].Laser cooling of trapped ions has become a rich methodology [19], which in recent years has made it possible to cool even large ICCs consisting of dozens [20][21][22] and even hundreds [23] of ions near their ground states of motion.Beyond cooling, an equally significant challenge is to accurately mea- sure the final state of motion of ICCs achieved by a particular cooling method and to determine its precise effective temperature.Accurate and efficient thermometry is important for evaluating spurious effects in quantum technology, such as quantum gate errors or systematic clock shifts, as well as for evaluating the efficiency of cooling schemes. is the mean excitation probability averaged over the ion crystal under red and blue sideband drive, respectively.The absorption probability is much smaller on the red vibrational sidebands than on the blue sidebands.Yet, the motional modes are still far from being cooled close to their ground states with temperatures ranging at about 5-10 phonons on average.(b) Simulated dynamics of the mean electronic red/blue sideband excitation P r,b of a 19ion ICC (center-of-mass mode, assuming n = 5).(c) Taking P r /(P b − P r ) as an estimate of the mean phonon number, as suggested by the sideband ratio method (Eq.( 3)), yields completely erroneous results and a significant underestimate of the motion temperature.
There exist well-established methods for thermometry of trapped ions, yielding reliable results for ion Coulomb crystals (ICCs) in high-temperature thermal states and for single or few ions near the ground state, as summarized in Fig. 1.Far from the motional ground state, timeof-flight, Doppler lineshape or image analyses achieve satisfactory accuracy [24][25][26][27][28] and suit practically any ion numbers.Close to the motional ground state, the excellent control over the coupling of motional degrees of freedom to internal levels, combined with the extreme precision that can be achieved in measuring the latter, can be used to perform an indirect temperature measurement.For single trapped ions, techniques such as singular value decomposition and numerical fits are employed to analyze motional sideband transitions [29,30].The primary tool for single ion thermometry is the resolved sideband ratio technique [31], which exploits the pronounced asymmetry of motional sidebands near the ground state.
However, these techniques cannot be directly applied to large ICCs near the ground state, posing an open challenge for achieving reliable thermometry in this regime, cf.Fig. 1.The problem is that for globally addressed crystals, driving a first-order (red or blue) sideband transition induces strong and complex spin-spin correlations due to their joint coupling with the driven normal mode, much like in a trapped ion quantum gate.These correlations significantly affect the asymmetry between red and blue sideband transitions, as is illustrated in Fig. 2 at the example of the sideband spectrum of a Dopplercooled 19-ion crystal.The marked asymmetry in the mean excitation probabilities on the red and blue sidebands would naïvely suggest a mean phonon number of approximately n ≲ 1 when applying the sideband absorption ratio technique [31] (cf.Eq. ( 3)) to this scenario.However, this estimate significantly underestimates the actual mean phonon number, expected to be around 5 to 10 after Doppler cooling.This clearly shows that spinspin correlations must be accounted for to accurately determine the temperature of motional states based on a measurement of internal state populations.For small ion crystals these correlations can be described exactly by solving the dynamics numerically or analytically [32][33][34].However for large ICCs consisting of tens or hundreds of ions, an exact consideration of ion-ion correlations is numerically intractable, as in any quantum many-body problem.Despite this complication, the resolved sideband technique (as well as other techniques based on spin-motional coupling) has been used in recent experiments [21,22,35,36] where many-body correlations have been either partially or completely neglected.The occurrence of many-particle dynamics in sideband thermometry can be trivially suppressed if only a single ion in a crystal can be interrogated.In this case, the known single-ion techniques are generally applicable, but become inefficient for larger ICCs due to poor statistics and long interrogation times required.When many or all ions of an ICC are interrogated, many-body dynamics can be easily accommodated in the exceptional case where the symmetric centre-of-mass (COM) mode is interrogated [23,37,38].In this case, the dynamics preserves permutation symmetry and shows a growth of the effective dimension of the Hilbert space that is only polynomial, instead of exponential, in the number of ions.Accurate thermometry of the out-of-phase modes, i.e. all modes except the COM mode, remains generally out of reach.
In this work, we present a new method for the thermometry of states of motion for arbitrarily large, globally addressed ICCs cooled close to the ground state.The method is based on measuring the ratio of the first-order excitation probabilities of the red and blue motional sidebands.Our technique takes into account the entanglement between the individual spins and remains accurate regardless of the number of ions.Using this global sideband thermometry method, we can extract the temperature of each motional mode that is assumed to be in a thermal motional state at the end of the cooling cycle.
Our approach to the many-body aspect of the problem may also serve as a useful reference and possibly complement the other existing thermometry methods.We demonstrate our method on a linear four-ion ICC and verify the estimated result with a full numerical simulation.We also test the global sideband thermometry on a two-dimensional 19-ion ICC; in this case the comparison with the numerical simulation is only possible for the COM mode.
The article is organised as follows: We start with a theoretical description of the resolved sideband thermometry of a single ion in section II A and place it in the framework of parameter estimation.In section II B we present our new global sideband therometry method for ICCs and discuss its advantages and limitations.In section II C we describe an alternative bichromatic technique that could be used if single ion addressing is possible and compare it with our new method.In section III we demonstrate the new global sideband method on four modes of motion of a linear 4-ion crystal and verify the results.In section IV we apply the new technique to a 2D 19-ion crystal and verify it.We summarise the results and give an outlook on open questions and possible further developments in section V.

A. Sideband Thermometry of a Single Ion
We consider a trapped ion with two relevant internal states |↓⟩ and |↑⟩, harmonically bound in a quadrupole trap with quantized motion along the three main trap axes.We restrict ourselves to a single axis and assume the ion is prepared in a thermal state of motion ρ(n) = ∞ n=0 p n (n) |n⟩⟨n| with Fock state occupation probability The mean occupation number is n = Tr ρ(n)a † a where a and a † are bosonic creation and annihilation operators for the motional degree of freedom.The mean occupation number can be inferred by coupling the motional state to the internal states and measuring the excitation probability.This can be realized by initializing the ion in ρ(n) ⊗ |↓⟩⟨↓| and illuminating it for a time t by a laser driving resonant transitions on either the red or blue sideband.The dynamics in these cases is governed by the Hamiltonians respectively.Here, the effective Rabi frequency g = Ωη/2 is obtained by rescaling the carrier Rabi frequency Ω with the Lamb-Dicke factor η ≪ 1 and the spin operators are defined as σ + = |↑⟩⟨↓|.The probability to find the ion in the excited state |↑⟩ after an interrogation time t is given by P α (n, t) = Tr U α (t)ρ(n) ⊗ |↓⟩⟨↓| U † α (t) |↑⟩⟨↑| , where α = r, b for red or blue sideband dynamics, respectively, and the time evolution operators are given by U α (t) = exp(−iH α t).For convenience, the excitation probabilities are given explicitly in Eq. (A2).Several examples of P α (t, n) are plotted at specific values of n in Fig. 3(a).
A measurement of the excitation probability on both the red and blue sideband transitions as a function of interrogation time gives access to the mean occupation number via, e.g., a numerical fit of the data to P α (n, t), [34].Apart from this, one can also infer the temperature from measurements at a single interrogation time by using the convenient formula This identity is at the heart of the sideband ratio thermometry method [31,39,40].It can be easily verfied using the expressions for P α (n, t) in Eq. (A2) and the fact that pn+1(n) pn(n) = n n+1 for thermal states.Eq. ( 3) is formally a correct analytical statement relating the measured excitation probabilities to the motional temperature.It is, however, important to note that this formula implicitly suggests that the values of P α are the precise statistical averages and hence only holds true if the data comes from an infinitely large measurement sample.In reality, the sample size is finite.The real measurement outcomes are the relative statistical frequencies f α , which will inevitably deviate from the true excitation probabilities.As the sample size N increases, this deviation goes down following a Gaussian law: where ξ α is a normally distributed random variable, ξ α ∼ N (0, 1).Since f α slightly differ from the real excitation probabilities, plugging these values into Eq.(3) results in an expression which is not exactly the motional temperature n, but its statistical 'estimator': The values of the estimator (referred to as estimates for short) are distributed with a certain statistical uncertainty and might carry a bias, which has to be accounted for when applying Eq. (3) to experimental data.Using Eq. ( 4) we calculate the statistical bias and the variance of the estimates: Both the variance and the bias converge to zero with the number of measurements.In Fig. 3(b, c) the intrinsic bias and the relative statistical uncertainty of the estimator in Eq. ( 3) are plotted as a function of the interrogation time and rescaled to N .As one would expect, both bias and the statistical uncertainty are large at small interrogation times, where the obtained statistics is poor due to the weak distinguishability of the sidebands.The bias and uncertainty also strongly scale with mean occupation number n.For example, achieving an uncertainty of 3% at n = 0.5 requires N = 10 4 measurements at an optimal interrogation time.An optimal interrogation time is found near the time of the maximal sideband excitation.This is the operating point, where the estimation using Eq.(3) converges the fastest to the true temperature value and has minimal bias.It is also interesting to compare the statistical uncertainty of the estimator with the fundamental bounds given by the (quantum) Cramér-Rao bounds.Since these observations are not central to the thermometry problem we consider below, we defer them to Appendix A.

B. Sideband Thermometry of an Ion Crystal
In an ICC, the motion of ions is collective and occurs in normal modes [39].In the following, we operate under the commonly used premise that after laser cooling, each mode is in a thermal state characterised by a certain mean phonon number n.This is usually a good approximation, although there may well be deviations, especially for short cooling times [30,41,42].Perform-ing ion crystal thermometry thus amounts to estimating the mean phonon number of the collective modes of motion by performing spin measurements on the crystal.
To date, many of the used temperature estimation techniques operating in the nearly ground-state cooled regime (n ≲ 1) rely on knowing the exact system dynamics governed by the first-order red and blue sideband Hamiltonians, arising in the Lamb-Dicke regime of the atom-light interaction for a globally-addressed crystal.These types of coupling entangle the crystal's particular mode of motion with the spins and are given by with (pseudo-)collective spin operators where N is the number of ions in the crystal.Here, g denotes an average Rabi frequency of sideband transitions.The average is chosen such as to have i η 2 i = 1, where the factors η i account for all inhomogeneities in the couplings of the ions to the mode under consideration.Thus, η i comprises the net effect of Lamb-Dicke factors, varying Rabi-frequency etc., which we assume to be calibrated for a given setup.
To obtain the desired exact solution, one needs to timepropagate the Schrödinger equation with Hamiltonians (7).The complexity of this calculation grows exponentially with the number of ions N , making it impractical to extract the exact solution for large ion crystals in reasonable time.The only exceptional case is the symmetric center-of-mass mode, where all the individual coupling constants are equal, η i ≡ 1/ √ N .This allows one to solve the problem in a symmetric Hilbert subspace, thus lifting the exponential scaling, cf.Appendix D.
In order to resolve the large-N bottleneck for out-ofphase modes, a new temperature estimation method is needed, which does not require solving the Schrödinger equation for the dynamics of the coupled N -ion system.We intend to stay as close as possible to the concept of thermometry of a single ion, which we have described in Sec.II A. To extend this approach, we first need to define the excitation probabilities for an ion crystal.A naive way to do this would be to use the mean electronic excitation of the ions, However, adopting this definition of the excitation probability for sideband thermometry leads to an incorrect result, and the temperature is greatly underestimated, as already shown in Fig. 2.
We will show that it is much more efficient to define the global crystal excitation probability as the complement of the probability to remain in the ground state of all ions, i.e.
Here, |0⟩ = |↓ . . .↓⟩ is the collective spin ground state, and propagated full density matrix of the system for red/blue sideband excitation.
Defining the excitation probability in this way has two benefits: firstly, the measurement implied by Eq. ( 8) can be easily performed and does not require single ion resolution.Secondly, for sufficiently small n and t, the functional dependence between the sideband ratio and n can be computed efficiently, taking into account the entanglement between the ions in the crystal.The generalization of Eq. ( 3) for ICCs is where ) .
(9b) The P k (n) are certain known polynomials in n of order k with coefficients depending solely on the mode coupling coefficients η i .Their explicit form and details of their derivation can be found in Appendix B. Further below we will comment on the idea behind this calculation.
In Fig. 4 we plot R t (n)/n as a function of time for temperature values in the regime of interest.The curve's deviation from the values of 1 thus shows the relative temperature estimation error when naïvely applying the single-ion formula (3) to an ICC.As case studies, we plot the curves for a 1D 4-ion crystal as considered in Section III and for a 2D 19-ion crystal considered in Section IV.The chosen modes are the COM-mode and one representative out-of-phase mode.As is evident from Fig. 4, applying formula (3) to an ICC will result in a relative error of roughly 20%, depending on the chosen mode and the crystal interrogated.Comparing this result with Fig. 2(b,c) shows that defining the global excitation probability as per (8) significantly improves the temperature estimation using the single-ion formula (3) compared to using the mean electronic excitation.The new estimator in Eq. ( 9) accounts for the residual discrepancy and thus mitigates the systematic error.
In the spirit of the thermometry of a single ion discussed in Sec.II A, the temperature estimator appropriate for the normal mode of an ICC can be constructed by replacing the excitation probabilities on the left hand side of Eq. ( 9) by measured relative statistical frequencies, and solving the resulting equation for n.Thus, with the (properly chosen) root of the quartic polynomial R Eq. ( 10) is the sought for generalization of Eq. ( 5) to an ICC.The systematic bias and the estimation error for finite sample size can be computed as in the case of a single ion, and are given in Eqs.(B3) and (B4) of the Appendix B. Only sideband excitation probabilities and R t (n) with its derivatives enter the formula Eq. (B3).Hence, no new data needs to be gathered to perform the bias correction.Before discussing the properties and limitations of this estimator, we sketch how Eqs.(9) are derived.Firstly, we exploit that both Hamiltonians in Eqs.(7) have a conserved quantity, namely H r , a † a + J 0 = H b , a † a − J 0 = 0, where J 0 = N i=1 σ i + σ i − measures the number of spin excitations.In both cases, α = r, b, this entails for the probability in Eq. ( 8) of remaining in the spin ground state that where p n (n) is the thermal occupation probability according to Eq. ( 1).Thus, only diagonal components of the time evolution operator enter the excitation probability P α (n, t).Secondly, we use that it is sufficient to describe the time dependence of the excitation probabilities up to their first fringe, as is evident from the discussion in Sec.II.This observation holds true for the ion crystal case as well.To exploit this, the diagonal matrix elements of the evolution operator are expanded in a power series in time.The series is then truncated at eighth order, since it is found to be enough to cover the significant fraction of the first sideband excitation fringe: Note that only even powers of the sideband Hamiltonians make a non-vanishing contribution.The relevant matrix elements ⟨0, n| H 2k α |0, n⟩ are polynomials in n of degree k, and comprise the many-body dynamics generated by the sideband drive.The matrix elements can be evaluated analytically, and the lowest two orders are given in Appendix B. The expressions for the cases k = 3, 4 are rather lengthy, and are evaluated by means of computer algebra, which is available at [43].With the approximation in Eq. ( 12), the matrix elements of the evolution operator become polynomials of order four in n and eight in gt.The average with respect to n in Eq. ( 11) can then be taken exactly, and yields the excitation probabilities in Eq. ( 8) as polynomials in n and gt, still of order four and eight, respectively.Since P α (n, 0) = 0, the final result for the ratio R t (n) in Eq. (9a) is correct within sixth order in gt.
Since the new estimator (10) relies on a timetruncation of the dynamics, cf.Eq. ( 12), the temperature estimation will be reliable only up to a certain time, which should cover a significant fraction of the first fringe in the excitation probabilities.In order to investigate this more quantitatively, we define a 'cutoff time' t * at which the estimated mean phonon number deviates from the true value for more than a predefined error threshold ϵ, which we chose to be ϵ = 5 • 10 −3 .For small ion crystals, t * can be calculated numerically and its dependence on the size of the ICC can be investigated.In Fig. 5 the cutoff time is shown for all motional modes for a case study of ion crystals containing N = 4 to 12 ions.The results show that there is no significant dependence of t * on the motional mode index or on the number of ions, and that the cutoff time is sufficient to measure an excitation signal with good signal-to-noise ratio on both motional sidebands.We also observed no tendency for the cutoff time to significantly decrease when increasing the temperature within the regime n ≲ 1.Since no assumptions on the crystal size were made along the derivation, one may infer that the proposed global sideband temperature estimator of Eq. ( 10) remains valid for large ion crystals.
In summary, the estimator in Eq. ( 10) provides a suitable extension of the well-established sideband thermometry to large, cold ICCs.It allows collective addressing and readout of the ions, providing fast dynamics and a strong signal, and adequately reflects the many-particle correlations involved.In the following, we will indicate yet another approach to ICC thermometry which also allows for collective addressing, but exploits single ion readout in order to avoid the complications connected to many-body dynamics.However, global sideband thermometry as discussed in Sec.II B gives better statistics at low temperatures, as we will show.

C. Thermometry from collective bichromatic drive and single ion readout
When the red and blue sidebands are driven simultaneously, the dynamics of the ICC follows the Hamiltonian x .This can be exploited to avoid the difficulty of dealing with complex many-body interactions, when the readout can be done via a particular single ion of the crystal.For a crystal initially prepared in the state ρ 0 = |0⟩⟨0|⊗ρ(n), the probability to find atom i in the excited state after a time t of bichromatic driving is This shows an exponential loss of contrast at a rate determined by the sought-after mean phonon number n, independent of the exact dynamics of the other ions in the crystal.It can therefore form the basis for an estimator of the motion temperature without having to consider particle correlations.However, the necessary interrogation time will depend strongly on the chosen ion via the mode coefficient η i , and might get large as the crystal size is increased.Moreover, interrogation of a single ion will require larger statistics.
A quantitative comparison between thermometry based on bichromatic drive and the approach to global sideband thermometry can be obtained by considering their statistical uncertainties and their Cramér-Rao (CR) lower bounds.The latter follow directly from the probabilities in Eqs. ( 8) and (13), respectively.In each case, the interrogation time can be optimized to achieve minimal estimation error at a given motional mode temperature n.For the bichromatic drive, the CR bound is independent of the number of ions and the specific mode.For the global sideband thermometry we use as an example the COM mode of a 10-ion crystal, and compare the CR bound also to the error for the specific estimator in Eq. (10).The results are shown in Fig. 6 and demonstrate that the two methods show advantages in complementary regimes: for low temperatures, sideband thermometry will yield better statistics, while for higher temperatures corresponding to n ≳ 1 the bichromatic approach is more efficient.The sideband ratio estimator does not saturate its CRB, yet it gets close to the CRB curve in the limit of small n, while diverging from it for larger n.This is mostly due to the fact that outside of the regime n ≲ 1 the cutoff time starts to decrease significantly with growing temperature, shifting the statistical uncertainty minima to higher values.As the cutoff time does not cover the optimal interrogation time required for the CR bound, there is space for getting closer to the bound with a higher cutoff time, which could be achieved by increasing the truncation order in (12).

III. THERMOMETRY OF A LINEAR 4-ION CRYSTAL
To test the new global sideband themometry method described in Sec.II B we measure the motional temperature of a linear ICC of four 172 Yb + ions.The size of the crystal allows us to evaluate the sideband dynamics exactly and thus benchmark the new method by comparing it to a direct numerical fit.

A. State preparation and cooling
The crystal is stored in a segmented linear radiofrequency (rf) Paul trap [44,45].The radial confinement, i.e. in the xy plane, is set by an rf electric field driven at Ω rf = 2π × 24.4 MHz, which is supplied to the trap electrodes by a resonant circuit.The axial confinement, i.e. along z axis, is set by a combination of dc voltages supplied to the trapping segment and both neighbouring segments.The corresponding secular frequencies are ω x,y,z = 2π × (586, 666, 111) kHz.The 12 motional modes along the x, y and z directions are cooled to the Doppler temperature of about 0.5 mK on the dipole allowed 2 S 1/2 → 2 P 1/2 transition near 370 nm, assisted by a repumper laser near 935 nm.The ions are detected individually by collecting the fluorescence from the decay of the short-lived 2 P 1/2 state.With a high numericalaperture lens of N/A = 0.27 the light from individual ions is imaged onto an electron multiplying charge-coupled device (EMCCD).For more details on the experimental apparatus, see [17,44,46].
The four modes along the strong radial axis of the trap (y) are further cooled to near the ground state using quench-assisted resolved-sideband cooling on the 2 S 1/2 → 2 D 5/2 electric quadrupole transition near 411 nm and the dipole allowed 2 D 5/2 → 2 P 3/2 near 1650 nm.The 411 nm beam is derived from a 822 nm laser that is locked to a cavity with a fractional instability of 5 × 10 −16 at 10 s of averaging time.It propagates parallel (within 3 • ) to the strong radial axis of the trap and efficiently addresses only the corresponding radial y-modes.It is focused down to a waist of 50 µm at the position of the ions and is aligned to the center of the crystal by measuring the carrier Rabi frequency of all four ions individually.After optimization, the Rabi frequencies are measured to be Ω Rabi [ion1, ion2, ion3, ion4] = 2π × [10.66(6), 10.61(6), 10.58(6), 9.88(3)] kHz, which varies by at most 10% over the crystal.
The exact mode frequencies of the four radial y-modes are measured to be ω y [mode1, mode2, mode3, mode4] = 2π × [666.0(1),656.9(1), 643.1(1), 623.6(1)] kHz with resolved sideband spectroscopy on the 411 nm transition.These mode frequencies are used to calculate the Lamb-Dicke factors for the motional modes, where the COM mode is at the highest mode frequency.In order to cool all modes simultaneously, the frequency of the 411 nm laser is set to be 640 kHz red-detuned from the carrier transition, such that it is roughly at the center of the four measured mode frequencies.The 1650 nm laser propagates along the axis of the trap and the power is tuned to reach an effective linewidth of the three level system of 67(2) kHz, see also [42].

B. Sideband thermometry measurement
After ground-state cooling, a thermometry measurement is performed on each motional mode along the y direction.The corresponding red and blue sidebands on the 2 S 1/2 → 2 D 5/2 transition are addressed to measure the excitation probabilities P r and P b , as defined in Sec.II B. For simplicity, the electronic states of individual ions are denoted as 2 S 1/2 i = |↓⟩ i and 2 D 5/2 i = |↑⟩ i .The in- ternal state of each individual ion is measured spatially resolved after a sideband pulse using the electron shelving technique, i.e. fluorescence was only detected on the 2 S 1/2 → 2 P 1/2 transition if the ion was in the |↓⟩ state.
where α = r, b.For a specific interrogation time, the excitation probabilities P r and P b are obtained by averaging over N /2 = 200 measurements and the interrogation time is scanned from 10 µs to 800 µs.As an example, the data obtained for mode 3 is shown in Fig. 7. Since for a 4-ion crystal the dynamics can be obtained numerically, we can obtain a temperature estimation by fitting the experimental data with the simulated curves of the sideband flops and searching for the optimal temperature values (shown in the right panel of Fig. 8).Further description of the fitting method, together with the data of the other motional modes can be found in Appendix D.
For each pulse time below the cutoff, the global sideband ratio is calculated from P r and P b according to equation (9a) and individual estimates for ni are obtained using equation (9b), see left panel of Fig. 8.
The cutoff time, as defined in Sec.II B, for all motional modes in this 4-ion crystal is around 400 µs (gt ≈ 0.22).In order to avoid possible implicit biases in the single measurements of ni at a given interrogation time, we use data from all available points up to the cutoff time to estimate the temperature.This set of m = 6 individual ni estimations, which are bias-corrected for N /2 = 200 according to Eq. (B3) and carry individual error bars σ i , are averaged to obtain the final estimation according to the following weighted sum: The estimations for the extremely small interrogation time of 10µs (the seen very left point for modes 3 and 4) and their uncertainties fall outside of the plot range.The results are compared to the values obtained from a least-square numerical fit from the 'Rabi flops' of the red and blue sidebands.Agreement is found between the two methods within 1σ.
The final estimation of n for each motional mode is n = {0.22 ± 0.05, 0.27 ± 0.02, 0.32 ± 0.03, 0.35 ± 0.04} and is shown in the right panel of Fig. 8.The results are compared with the estimations obtained from the numerical fit of the temperature.
Since a good agreement is found between the theoretical and experimental curves, the results from the fit are expected to give an accurate estimation for n.For all the motional modes, the values of n obtained from the global sideband method agree with those extracted from the fit within a 1σ uncertainty.The global side- band ratio method reaches the same accuracy level of δn ∼ 10 − 20% as the numerical fit, but requires less data points to be taken.We want to explicitly emphasize that the used least-square fit relies on the exact calculation of the sideband dynamics, which scales exponentially with the number of ions and thus cannot be applied for larger ion crystals in general.

IV. THERMOMETRY OF A 2D CRYSTAL A. Setup, state preparation and cooling
To demonstrate the new scheme also on a larger ion crystal, which already imposes challenges in numerically simulating the sideband dynamics, we perform thermometry measurements on the out-of-plane motional modes of a two-dimensional 40 Ca + ion crystal.A planar 19-ion Coulomb crystal is stored in the anisotropic potential of a novel microfabricated monolithic linear Paul trap de-signed for trapping 2D ion crystals [47].The trap is operated at oscillation frequencies of ω s = 2π × 2.189 MHz, ω w1 = 2π × 645 kHz and ω w2 = 2π × 340 kHz where ω s is the secular frequency in the strongly confining direction and ω w1 and ω w2 are the secular frequencies in the two weakly confining directions.The direction of the weakest confinement is aligned with the rf-zero line.Our geometry allows micromotion-free optical access in the plane spanned by the directions of ω s and ω w2 .
Ions are Doppler-cooled on the 4S 1/2 ↔ 4P 1/2 dipole transition at 397 nm.Electromagnetically-induced transparency (EIT) cooling [48,49] is applied for 300 µs after Doppler cooling to simultaneously cool all N out-of-plane secular modes of motion in an N -ion crystal close to the ground state.The direction of the magnetic field, defining the quantization axis, is oriented at an angle of 45 • with respect to the crystal plane and allows for an optimal geometry for EIT cooling.A more detailed description of the beam geometry is given in Ref. [47].
The laser used for EIT cooling is blue-detuned by 110 MHz from the 4S 1/2 ↔ 4P 1/2 transition.The chosen detuning enables efficient cooling over a frequency range large enough to accommodate all transverse modes of motion of a 19-ion crystal, spanning a range of ∼ 300 kHz.The power of the σ − -polarized beam is calibrated such that the induced Stark shift overlaps with the center of the frequency range to be cooled.Further details on the calibration procedure can be found in Ref. [50].
The motional modes are probed via sideband spectroscopy on the 4S 1/2 ↔ 3D 5/2 quadrupole transition with light from a frequency-stable laser (∼ 1 Hz linewidth) at 729 nm.A global beam along the out-ofplane direction excites the individual ions with a maximum variation in the single-ion carrier Rabi frequencies of about 6 % across the 19-ion crystal.For a single ion, we find n = 0.06 for the transverse mode of motion (ω s = 2π × 2.188 MHz) corresponding to the out-of-plane direction of a planar multi-ion crystal.A frequency scan of the red sideband spectrum of a planar 19-ion crystal, once after only Doppler cooling and once after an additional EIT cooling pulse of 300 µs, is shown in Fig. 9.

B. Sideband thermometry of a planar 19-ion crystal
Sideband thermometry based on Eq. ( 10) is applied to the COM mode and the lowest-frequency mode in the out-of-plane direction of a two-dimensional 19-ion crystal.We calculate the normal-mode frequencies and the Lamb-Dicke factors of the individual ions for all out-ofplane modes using simulations within the pseudopotential approximation.The knowledge of both is required for the employed temperature-estimation method.However, the pseudopotential approximation can sometimes fail to reproduce the observed normal mode frequency spectrum, in particular for planar crystals [51,52].We measure the out-of-plane motional-mode spectrum and observe a good match with the simulated mode frequencies (see Fig. 9), providing us with confidence that the pseudopotential approximation yields accurate results for the 19-ion crystal.
Sideband excitation dynamics of the COM and the lowest-frequency mode with a carrier Rabi frequency of 2π × 38.8 kHz are shown in Fig. 10.Single-ion-resolved measurements reveal the structure of the investigated modes, shown as P i r,b in Fig. 10(a,b), and yield the global excitation probabilities, P b and P r , shown in Fig. 10(c,d).
Thermometry measurements are carried out after EIT cooling where we measure the temperature of the planar ion crystal as a function of the carrier Rabi frequency and the interrogation time.The excitation probabilities P b and P r are probed in 4000 individual experiments each.The mean phonon number is then calculated according to Eq. (10).
Using high laser intensities can lead to crosstalk to neighboring modes, which results in inaccuracies in the phonon-number estimates.For the lowest-frequency mode, we observe some crosstalk to the neighboring mode separated by less than 30 kHz.We can circumvent this problem simply by probing the modes at lower Rabi frequencies.Figure 10(e,f) shows the estimated mean phonon numbers for different values of gt.Consistent results are obtained across the probed parameter space.The bias-corrected (Eq.B3) weighted mean phonon number of the COM mode and lowest-frequency mode are determined to be n = 0.149(3) and n = 0.069(3), respectively.
To cross-check the measured mean phonon number of the COM mode, we simulate the sideband dynamics of the COM mode in the symmetric Hilbert subspace and, as for the 4-ion data, perform a weighted least-square fit of the measured data on the red sideband of the COM mode.More details on the numerical calculation for the symmetric COM mode as well as the fit estimator are given in Appendix D. The fit yields a temperature estimate of n = 0.147±0.02,which is in good agreement with the value obtained using the new sideband thermometry technique.The theoretical curves for the red and the blue sideband are shown as solid lines in Fig. 10(c).Up to ∼ 350 µs we find a good agreement between theory and experiment.For longer probe times, we observe deviations from the simulated curve, which we attribute to motional heating as well as instabilities in the trap oscillation frequencies due to power fluctuations.We thus use the first 20 data points (< 350 µs) to fit the dynamics of the COM mode.A numerical simulation of the sideband dynamics of the LF mode, however, is computationally demanding.Therefore, a reference value for the lowestfrequency mode is not given.
In order to test the method with phononically higher excited states, a heating rate measurement of the COM mode is performed in which the probe pulse is delayed by a predefined wait time between 0 and 20 ms after groundstate cooling.The heating-rate curve in Fig. 10(g) shows the bias-corrected estimated values of n.The data is fitted with a linear function by least-squares, weighted with the inverse variance obtained from Eq. (B4).The fit reveals a heating rate of 15.3(1.7)quanta/s per ion consistent with previous measurements with a single ion as well as an 8-ion crystal.In contrast to the COM mode, measurements on the lowest-frequency mode do not indicate significant heating within tens of milliseconds, as expected.

V. DISCUSSION & CONCLUSIONS
We have presented here a method for the thermometry of cold ICCs that generalizes the well-known sideband thermometry of single ions.The effects of manybody quantum dynamics that arise in this process can be taken into account with sufficient accuracy by exploiting conservation quantities of the sideband dynamics and by a suitable truncation of its time series expansion.It turns out that this limitation is not critical, since a sufficient signal can be extracted within a time span in which the truncation still yields reliable results.As we show, the tolerable interrogation time does not change with the number of ions.In principle, if required, a higher truncation order can also be achieved using the methods we have presented here, whose implementation in Python and Mathematica can be accessed at [43].Applications of this methodology to a linear as well as a planar ion crystal give good results, also in comparison to other methods, in cases where such a comparison is possible.A reliable tool for temperature measurement in ultracold ion crystals is an important requirement in experimental quantum metrology and information science.We believe that the thermometry method presented here meets the current needs and can be of practical use for research with cold trapped ions.Moreover, the approach presented in this work could serve as a useful reference for the treatment of many-body effects in similar systems.
As an outlook, we would like to indicate a number of questions and possible further developments that go beyond the results presented here: a central premise of sideband thermometry is the presence of a canonical thermal state.This is a useful and mostly very good approximation, but it will not always be fulfilled for all cooling methods and especially not for short cooling durations, which will occur in quantum technology applications due to time limitations [30,41].Our general approach would also allow to consider more general, non-canonical parameterisations of the occupation probability and to estimate the corresponding parameters systematically.For this, one has to consider a many-parameter estimation problem and in a similar fashion derive the corresponding estimators from the measured observables.Another possible extension would be to consider correlated spin states in order to achieve a quantum metrological improvement of the accuracy of the thermometric measurement.Finally, a way could also be sought for thermometry based on bichromatic driving to exploit measurements of more than one ion and account for the many-body correlations that occur there.
The data underlying the reported measurements are available via Zenodo [53].
An appropriate root of Eq. (B1) gives an estimator for the temperature of the motional mode.Although this equation generally has four roots, in practice it is easy to identify the one corresponding to the temperature estimation.The other roots are typically complex, negative or have values far outside of the n ≲ 1 region.If somehow the ambiguity is still present, dropping all time-dependent terms in the r.h.s. of (B1) and measuring the l.h.s. at the smallest possible time provides the simplest rough estimation helping to choose the correct root.
The asymptotic bias and the variance arising from the finite sampling of N /2 for both P r and P b are calculated as: ) where the derivative of R t (n) is understood as the derivative with respect to n.

Appendix C: Mode-dependent coefficients
The coefficients A, B i , C i and D i given in Eqs.(B2) depend only on the structure of the interrogated motional mode and are defined as follows: Essentially, the coefficients are all the non-vanishing expectation values of strings of operators J ± of a fixed length.Using the definition J ± = N i=1 η i σ ± i each of the coefficients is decomposed into a sum of many strings of particular single-atom Pauli operators.Each string is weighted with the corresponding prefactor, consisting of the mode vector components and sandwiched with the spin ground state.With some combinatorics one can classify and pick out the small number of non-zero terms to ease the calculation.There may be several ways to do so, one of them would be to distinguish the terms based on the number of unique atomic indices appearing in an individual term.For each mode-dependent coefficient we count the number of terms of each class (C j i , D j i ) and then multiply it with the corresponding expectation values.This brings us to the general expression for C-and D-coefficients given by   In order to use the new temperature estimation method, for a given motional mode one needs to com-pute in total 22 mode-dependent coefficients.Using formulas (C1, C2), the calculation boils down to programming a cascade of several FOR-loops, which results in a computationally-friendly polynomial scaling with respect to N .Hence, one may calculate the coefficients entering the estimator easily and fast to apply the new thermometry method for arbitrarily large ion crystals.We provide a program for calculating these coefficients for a given motional mode in the code supplement [43].
Appendix D: Numerical fit estimator A common way for estimating cold ion temperatures is fitting the experimental data with theory curves while using the temperature as a free parameter.Naturally, this technique requires the exact numerical solutions of the Schrödinger equation or its reasonable approximation to be available.This is, however, problematic for large ion crystals due to the exponential scaling of the Hilbert space, and causes the main bottleneck for the numerical fit.For smaller crystals, one can employ this method for temperature estimation using a weighted least square optimization of the model curves with respect to unknown parameter n.The choice of model functions may vary, though in our case the best performance was achieved when using the red sideband collective excitation probability P r (t, n).This choice is further backed by the previous observation from the single-ion case, that the Fisher information on parameter n contained in the red sideband is typically significantly higher than the one contained in the blue sideband or in the combinations of both.This observation also holds in the multi-ion case for relatively short interrogation times.Using the red sideband curves as model functions makes the numerical fit estimator in our case have the following form n = arg min n S(n) = arg min where the sum is taken over m experimental points x i , each carrying a normally-distributed error σ i .The variance of this non-linear weighted least-square estimator is given by [55] ∆n data.The results of the numerical fit together with the experimental data for all motional modes are shown in Fig. 12.
The symmetric center-of-mass mode deserves a special mention.Since all the individual coupling strengths are equal for this mode (η i = 1/ √ N ∀i), the spin dynamics evolves within the symmetric Hilbert subspace and could effectively be described using symmetric spin Dicke basis of states This lifts the exponential scaling of the Hilbert space with respect to the number of ions N and thus makes the numerical fit applicable for the center of mass mode of larger ion crystals.

FIG. 1 .
FIG. 1. Schematic overview of thermometry methods for cold ions.Well-established approaches exist for ion crystals at large temperatures.Close to the ground state, the temperature can be measured of single ions or small crystals.Development of methods applicable to large and cold crystals is an open challenge and subject of the present work.

FIG. 2 .
FIG. 2. (a)Red and blue sideband spectrum of the out-ofplane modes of a 19-ion planar crystal after Doppler cooling (for further experimental details, see section IV).P r,b =1 N N i=1 P i r,bis the mean excitation probability averaged over the ion crystal under red and blue sideband drive, respectively.The absorption probability is much smaller on the red vibrational sidebands than on the blue sidebands.Yet, the motional modes are still far from being cooled close to their ground states with temperatures ranging at about 5-10 phonons on average.(b) Simulated dynamics of the mean electronic red/blue sideband excitation P r,b of a 19ion ICC (center-of-mass mode, assuming n = 5).(c) Taking P r /(P b − P r ) as an estimate of the mean phonon number, as suggested by the sideband ratio method (Eq.(3)), yields completely erroneous results and a significant underestimate of the motion temperature.

FIG. 3 .
FIG. 3. (a)A simulated Rabi flop on blue (blue curves) and red (red curves) sideband for a single ion for several values of the mean phonon occupation number.Bias (b) and the statistical uncertainty (c) of the sideband temperature estimator (3), rescaled to the total number of measurements N .

FIG. 5 .
FIG. 5. Cutoff time gt * for N motional modes of a crystal containing N = 4 to 12 ions (for a linear chain with generic trap parameters).Each point corresponds to a particular motional mode.Due to the symmetries in the mode vectors, certain modes have coinciding cutoff time values.The upper points on the plot correspond to the COM mode.The data refers to n = 0.1.

FIG. 6 .
FIG.6.Cramér-Rao bound for temperature estimation in case of the sideband measurement on the whole crystal (bright green) and in case of the measurement based on bichromatic drive and interrogation of a single ion (dark green).For comparison, the statistical uncertainty of the global sideband temperature estimator(10) is plotted (black curve).Here N = 10, the center-of-mass mode is considered.The curves for other modes overlap almost completely with the COM-curves.

FIG. 7 .
FIG.7.Red (red points) and blue (blue points) global sideband flops of the third motional mode in a 4-ion crystal (see the text for details).The lines show the exact numerical solution generated with the least-square fitted temperature value n.

)FIG. 8 .
FIG. 8. Temperature estimation for motional modes 1-4 shown in subpanels (a)-(d), respectively, of a 4-ion crystal.The individual temperature estimations at interrogation times up to the cutoff time of ∼ 400 µs (left panel) are averaged together to produce the final value of n = {0.22 ± 0.05, 0.27 ± 0.02, 0.32 ± 0.03, 0.35 ± 0.04} shown in the right panel using the presented global sideband method.The estimations for the extremely small interrogation time of 10µs (the seen very left point for modes 3 and 4) and their uncertainties fall outside of the plot range.The results are compared to the values obtained from a least-square numerical fit from the 'Rabi flops' of the red and blue sidebands.Agreement is found between the two methods within 1σ.

FIG. 9 .
FIG.9.Red sideband spectrum of the out-of-plane modes of a 19-ion planar crystal after Doppler cooling (top) and an additional EIT-cooling laser pulse (bottom).The frequency is given as detuning from the 4S 1/2 , m = −1/2 ↔ 3D 5/2 , m = −1/2 carrier transition.The mode frequencies obtained from simulations using pseudopotential theory are displayed as dashed lines.A false-color image of the ion crystal is shown at the top center.The mode structure of the COM mode (top left) as well as the lowest frequency mode (top right) is indicated by red arrows.Each arrow's length is proportional to the magnitude of the respective ion's Lamb-Dicke factor ηi, the direction indicating the sign of ηi.The temperature of both modes is probed in the experiments described further below.

FIG. 10 .
FIG. 10.Sideband dynamics and thermometry of a planar 19-ion crystal.(a,b) Single-ion excitation probabilities P i r,b for the COM mode (a) and the lowest-frequency mode (b).As in Fig. 9, the insets indicate the mode vectors of the investigated modes.(c,d) Global excitation probabilities P r,b according to Eq. (8) for the COM mode (c) and the lowest-frequency mode (d) representing the quantities of interest in the sideband thermometry measurements.The blue and the red points correspond to the measured excitation on the blue and red sidebands, respectively.The solid lines in (c) are obtained from simulations in the symmetric Hilbert subspace using the least-square fitted value of n. (e,f) Sideband thermometry of the COM mode (e) and the lowest-frequency mode (f) for varying carrier Rabi frequency (indicated by color) as a function of the interrogation time.The mean values are shown as dashed lines.For the lowest-frequency mode (f) the mean value is calculated after discarding all data points lying outside a range of 1σ from the mean value obtained from all data points.The discarded data points (5 points with the highest values of n) correspond to measurements with higher Rabi frequency and are biased due to crosstalk to neighbouring modes (see main text for discussion).The error bars of individual thermometry measurements are obtained from Eq. (B4).(g) Measurement of n for the COM mode as a function of the probe pulse delay.The solid line shows a weighted linear fit used to determine the heating rate.
tables I and II are the prefactors completing the expressions (C1, C2), which are needed to evaluate each of the C-and D-coefficients.

1 F ( 1 ,FIG. 12 .
FIG.12.Red (lower points) and blue (upper points) global sideband flops of all 1-4 motional modes shown in (a)-(d), respectively, for the 4-ion crystal.The solid lines show the exact numerical solution, where n was fitted as free parameter with the weighted least-square method.