Polarization-gradient cooling of 1D and 2D ion Coulomb crystals

We present experiments on polarization gradient cooling of Ca$^+$ multi-ion Coulomb crystals in a linear Paul trap. Polarization gradient cooling of the collective modes of motion whose eigenvectors have overlap with the symmetry axis of the trap is achieved by two counter-propagating laser beams with mutually orthogonal linear polarizations that are blue-detuned from the S$_{1/2}$ to P$_{1/2}$ transition. We demonstrate cooling of linear chains of up to 51 ions and 2D-crystals in zig-zag configuration with 22 ions. The cooling results are compared with numerical simulations and the predictions of a simple model of cooling in a moving polarization gradient.


Introduction
Laser cooling is an important technique to prepare neutral atoms and trapped ions at low temperature [1,2,3]. Doppler cooling [4] enables the formation of ion Coulomb crystals [5,6], which are at the heart of many experiments on quantum computation, simulation, and precision measurements [7,8,9,10,11,12,13]. In order to coherently manipulate ions with high fidelity, their motion in the trap should ideally be completely frozen out; hence, various laser cooling methods have been developed and demonstrated in the past three decades [14]. Resolved sideband cooling [15,16], Raman sideband cooling [17] and cooling by electromagnetically induced transparency (EIT) [18,19] are routinely used to prepare trapped ions in their motional ground state.
While cooling to the ground state might be considered the ultimate goal of laser cooling, a cooling scheme should also be evaluated against other important metrics such as cooling speed, cooling range in frequency space, or the initial energy from which an ion can be cooled down to low temperatures. Conventional ground state cooling techniques, such as Raman and resolved-sideband cooling techniques are known to be well suited for achieving high ground state occupancy, but are in general slow in terms of cooling multiple modes. On the other hand, the EIT cooling scheme, which is the fastest method of cooling to the ground state, maintains that cooling rate only within a limited range of motional frequencies and becomes slow for low-frequency modes, which tend to disturb laser-ion interactions more strongly due to their larger Lamb-Dicke parameters.
The need for scaling up quantum information and quantum simulation experiments leads to new challenges in controlling both electronic and motional states of a large number of trapped ions. In trapped-ion clocks, the goal to increase the signal-to-noise ratio in a limited time interval leads to the development of multi-ion clocks, which require sub-Doppler cooling of ion strings to maintain their precision.
Laser cooling schemes often under-perform when dealing with large-size ion Coulomb crystals. A linear chain of ions can be created by operating a linear Paul trap at low axial confinement, which however gives rise to low-frequency collective motional modes with a strong occupation of high-lying phonon states. Planar crystals in Paul traps, which are of interest in the context of quantum simulation [20,21,22], present a different kind of laser cooling challenge due to the limited number of laser beam directions for which the laser does not couple to the radio-frequency-driven ion motion. The quest for high-fidelity quantum operations in a trapped-ion-based quantum processor motivates us to explore novel approaches for cooling Coulomb crystals containing tens of ions.
Polarization gradient cooling (PGC) [3] is a Sisyphus cooling technique, which has been widely used for cooling neutral atoms to very low temperatures [23]. However, despite early theoretical works investigating Sisyphus cooling of strongly bound atoms [24,25,26], the cooling scheme hasn't yet been extensively applied to cooling trapped ions. After an early demonstration [27], polarization gradient cooling of trapped ions has been reported only recently in a regime where the ions were strongly bound to the trap [28,29]. In experiments with up to four ions [28], the authors reported cooling ions to near the motional ground state from far outside the Lamb-Dicke regime (LDR) over a few hundred microseconds. In the current paper, we demonstrate polarization gradient cooling of long ion strings and of planar ion Coulomb crystals to sub-Doppler temperatures.
Our results demonstrate that quantum computation, quantum simulation, and precision spectroscopy experiments can benefit from polarization gradient cooling. Specifically, for initializing the state of a trapped-ion quantum processor, many collective vibrational modes need to be cooled to low vibrational quantum numbers in order to avoid coupling strength fluctuations in quantum gates based on laser-ion interactions. Here, polarization gradient cooling can improve the gate fidelity without significantly increasing the length of the experimental duty cycle. Similarly, multi-ion clock experiments [30] require ion strings with as little secular energy as possible in order to reduce second-order Doppler shifts [31,32]. In this context, polarization gradient cooling can be seen as a technique bridging the gap between Doppler cooling and ground-state cooling techniques targeting a small number of vibrational modes.
The current article is structured as follows. Section 2 reviews and extends a simple model of polarization gradient cooling that predicts the achievable cooling limit and cooling rate. Our experimental apparatus and measurement techniques are described in section 3. In section 4, we first present measurements carried out with a single trapped calcium ion before discussing results obtained with multiple ions in a 1D chain and a 2D crystal.

Theory of polarization gradient cooling
The theory of polarization gradient cooling of a bound atom has been developed previously [25,26], with a focus on j = 1/2 ↔ j = 3/2 transitions. Here we will briefly present the main elements of the semiclassical theory developed in ref. [25] and derive the cooling rate and cooling limit for the case of a j = 1/2 ↔ j = 1/2 transition.
We consider an ion of mass m confined in a harmonic potential with trap frequency ω z . The ion has a dipole-allowed transition connecting the j g = 1/2 ground state to an excited j e = 1/2 state. The transition is off-resonantly excited by two counterpropagating laser beams with linear-perpendicular-linear (lin-⊥-lin) polarizations. This arrangement of laser polarizations forms a periodically varying polarization gradient along the direction of the laser propagation [3]. For simplicity, we assume that the cooling laser beams propagate in the direction of the quantization axis (ẑ) such that the ion interacts only with σ-polarized laser photons.
For an ion moving along theẑ axis, its Zeeman ground states |± = |J g , m = ±1/2 experience periodically varying ac-Stark shifts that are π out of phase with respect to each other. In conjunction with the axial trapping potential U trap = 1 2 mω 2 z z 2 , these light shifts give rise to a state-dependent total potential energy given by where k is the wave number, ∆ the detuning from the dipole-allowed transition of the laser field and φ determines the position of the polarization gradient with respect to the trap center (for φ = 0, an ion at the trap center is driven by linearly polarized light). The intensity of the cooling beams is expressed in terms of the saturation parameter, s = Ω 2 /2 Γ 2 /4+∆ 2 , with Γ the linewidth and Ω the Rabi frequency that the ion would experience on a |+ ↔ |J e = 3/2, m = 3/2 transition when interacting with only one of the two laser beams whose polarization was changed to being σ + -polarized. In this case, the excited-state population would be given by p = s/2 1+s for the definition of the Rabi frequency that we have chosen (i. e., a π-pulse requiring a duration of τ = π/Ω).
In addition to creating state-dependent potentials, the polarization gradient also gives rise to spatially dependent optical pumping rates between the ground states given by Γs(1 ∓ sin(2kz + 2φ)). (2) Transition probability z Figure 1. Principle of polarization gradient cooling. A trapped ion with an S 1/2 ground state in a polarization gradient coupling the two Zeeman ground states to an excited state. The polarization gradient makes the trapping potential state-dependent, by displacing it into opposite directions for the two Zeeman ground states. It also induces spatially varying pumping rates between the states. These two effects enable cooling a hot ion by favouring energy-reducing transitions (black squiggly arrows) over the reverse processes.
It is the interplay of spatially dependent light shifts and pumping rates that gives rise to an efficient cooling mechanism, sketched in Fig. 1. We assume that the potential of the ion trap is much stiffer than the potential of the polarization gradient so that the latter creates state-dependent potentials U ± (z) for the two Zeeman ground states |± that can be treated as being harmonic. If the potential minima of U ± are spatially displaced, a hot ion can be cooled because the spatially varying pumping rate will favour state changes giving rise to an energy loss (as indicated in the figure) over the reverse processes. Note that for this picture to hold, the pumping rates Γ ± have to be smaller than the ion's oscillation frequency. The ion motion will come to a steady state when the cooling is balanced by heating by spontaneously scattered photons and dipole force fluctuations, and additionally by scattering processes at random times in which the ion is pumped to the other Zeeman ground state. Following ref. [25], we calculate the cooling limit and cooling rate by a simple model, in which an ion having a motional energy E stochastically switches between its ground states with rates Γ ± . By calculating the average motional energy change induced by transitions to the other Zeeman state and accounting for motional heating by absorption and emission processes, one obtains the following differential equation describing the cooling dynamics: where p ± = 1 2 (1 ± sin(2kz + 2φ)) are the ground state populations in steady state, and ... + and ... − indicate the averages over the ion's spatial probability distribution in the potential U + and U − , respectively. The term H sc ω z accounts for the momentum diffusion due to spontaneously scattered photons and fluctuating optical dipole forces [33]. In the LDR, the resulting heating is conveniently described in terms of carrier and sideband processes. In this picture, absorption on the carrier transition followed by spontaneous emission on a first-order vibrational sideband accounts for heating by spontaneously scattered photons, whereas absorption on a first-order sideband followed by emission of a photon on the carrier transition accounts for dipole force fluctuations. At z = 0, these two processes lead to a heating rate H sc = H carr + H sb given by where η = k 2 /(2mω z ) is the Lamb-Dicke parameter and α = 1/3 results from the spatially isotropic spontaneous emission on the S 1/2 ↔ P 1/2 transition that we are considering. Note that the dependence on the phase φ of the polarization gradient differs for these processes as the carrier coupling is maximum in the anti-nodes of the standing waves whereas the sideband coupling is maximum in the nodes [34]. After expanding the pumping rates Γ ± and the potentials U ± to first order in kz, an evaluation of the first two terms of equation (3) gives rise to a term ∼ −W E cooling the ion if the polarization gradient is blue-detuned with respect to the atomic transition. Additionally, there is a heating term proportional to the energy that an ion at the minimum of U ± gains when being excited to the other potential U ∓ (and vice versa). This leads tȯ with cooling and heating rates H(φ) = 2 9 η 2 Γs(8ξ 2 cos 4 2φ + 2 + sin 2 2φ), and ξ = ∆s 3ωz . If the ion experiences perfectly circularly polarized light (φ = ±π/4), the cooling rate vanishes. The most efficient cooling occurs at φ = 0 when the ion is placed at the steepest slope of the optical potential. Then, the mean phonon number in steady state is given by A minimum phonon number of min( n 0 ) = 1/2 is obtained for ξ = 1 2 , when the optical well depth (2∆s/3) becomes equal to the trap frequency.
Experimentally, it is very challenging to achieve the sub-wavelength stability required to position the ion at a precise phase of the polarization gradient. However, there are two ways of overcoming this problem: Cooling of many-ion crystals. When cooling motional modes of ion crystals containing many ions, the ions will roughly uniformly sample the phases of the polarization gradient so that the exact positioning of the polarization gradient no longer matters. One can then average the heating and cooling rates, eqs. (7), (8), over φ in order to calculate the cooling limit which is minimized for ξ min = 5/6 ≈ 0.91, resulting in min n = 15 Cooling in a moving polarization gradient. For a single ion or a few-ion crystal, the cooling scheme can be made robust by means of a moving polarization gradient [35] that can be created by introducing a small frequency detuning δ between the two counterpropagating beams. This detuning should be higher than the cooling rate; otherwise the mean phonon number would adiabatically follow the steady-state value n φ which, on average, would lead to an increased motional energy and large shot-to-shot fluctuations caused by random values of φ at the end of the cooling pulse. On the other hand, δ < ω z is required to make the ion sample the potential well before the polarization gradient changes considerably. We have validated the predictions of the simple cooling model described above by numerically solving the master equation describing polarization gradient cooling. The simulations, presented in Appendix A, demonstrate that eqs. (7)-(10) accurately describe the cooling of an ion in the LDR.

Experimental apparatus and measurement techniques
A linear Paul trap with blade-shaped RF electrodes [36] is used for confining both linear strings and planar crystals of 40 Ca + ions. A detailed description of the experimental apparatus can be found in Ref. [37]. A magnetic field of B ≈ 4.18 G pointing along the trap's rf-zero line defines the ions' quantization axis and leads to a splitting of 11.714 MHz between the two Zeeman sublevels of the S 1/2 manifold. Ions are precooled using an elliptically shaped Doppler cooling laser beam that illuminates the ion chain under an angle of 45 degrees, and has an overlap with all collective modes of motion.
In the current article, we present experimental results of polarization gradient cooling along the axis of the ion chain (ẑ axis). Holes through the trap's tip electrodes, with a diameter of 0.5 mm, allow two counter-propagating laser beams to pass through alongẑ, thus coupling to the axial motion only. Polarization gradient cooling is performed on the S 1/2 ↔ P 1/2 transition at 397 nm. The laser used for Doppler cooling also provides the two counter-propagating beams that are blue-detuned by 210 MHz from the cooling transition. The beams are linearly polarized with orthogonal polarizations, creating a polarization gradient that alternates between right circular and left circular polarizations along the direction of light propagation. The detuning of the beams can be controlled independently by means of two acousto-optical modulators. In most experiments, their frequency difference was set to δ = 2π × 60 kHz to create a polarization gradient moving at a rate greater than the cooling rate.
We calibrate the intensity of the two cooling beams in the trap center separately by optically pumping the ions to one of the Zeeman ground states and subsequently measuring the equilibration time constant of the two ground state populations after one of the cooling beams is switched on. For this, Zeeman ground state populations are measured by transferring the respective population to the D 5/2 state with a narrow linewidth laser beam at 729 nm, which is also used for motional state analysis.

Results and discussion
We investigate polarization gradient cooling using the following experimental sequence: In the first step, calcium ions are prepared in the S 1/2 state and Doppler cooling is performed for 3 ms. Next, a polarization gradient cooling pulse is employed for up to 1 ms. To prevent optical pumping to the D 3/2 state by spontaneous decay of the P 1/2 state, a repumper at 866 nm is switched on during the Doppler and polarization gradient cooling steps. In the third step, ions are prepared in the |S 1/2 , m = 1/2 state by frequency-resolved optical pumping on the |S 1/2 , m = 1/2 ↔ |D 5/2 , m = 3/2 transition, in conjunction with the 854 nm laser beam, preparing the ions in the |S 1/2 , m = 1/2 sublevel with a probability greater than 99.9%. The final step involves mapping the information about the motional state onto the electronic states, via a 729 nm laser pulse on either a carrier or a first-order sideband of the S 1/2 ↔ D 5/2 transition followed by quantum state detection via fluorescence detection at 397 nm. Spatially resolved measurements of the fluorescence of an ion crystal are recorded with a CCD camera in order to assign the quantum state of individual ions.

Cooling of a single ion
We start by investigating polarization gradient cooling of a single ion. The cooling time constant is measured by applying cooling pulses of variable duration τ to a Dopplercooled ion, followed by motional state analysis. For the latter, we drive carrier or first-order sideband Rabi oscillations, which we fit with a thermal phonon distribution. An example of cooling an ion is displayed in Fig. 2 (a), which shows that the ion motion equilibrates to a steady-state value, which is a bit above the cooling limit predicted by eq. (10), and with a slightly lower rate. The resulting data is fitted with a function f (τ ) = n(∞) + ( n(0) − n(∞) ) exp(−W τ ) in order to extract the cooling rate, where parameters n(0) and n(∞) correspond to mean phonon numbers at τ = 0 and τ = ∞, respectively. Here, for ω z = 2π × 1088 kHz and ξ = 1.35 (14), we find a cooling rate W ≈ 66000 s −1 . Taking into account the slightly non-exponential decay that we observe, the data shows that equilibrium is achieved after less than 200 µs.
For a cooling pulse of a duration τ much longer than the cooling time constant 1/W , we measure the mean phonon number n in steady state as a function of the  (14). The simulated cooling dynamics (dashed black curve) for the experimental parameters is based on averaging eqs. (7) and (7) over all phases φ. The dashed red curve is a fit to the experimental data. (b) Mean phonon number n in steady state as a function of the normalized saturation parameter ξ = ∆s/3ω z at 1088 kHz (in red), 501 kHz (in black) and 217 kHz (in blue) axial trap frequencies, together with the cooling limit predicted by eq. 9 (dashed line).
(c) n as function of trap frequency ω z at s ∼ 0.0085 (15) , which is predicted to yield optimum cooling for ω z = 2π × 650 kHz. The dashed black line is predicted by eq. 9. The blue and red dashed lines represent the cooling limits of one-and three-dimensional Doppler cooling, respectively.
trap frequency and the intensity of the polarization gradient. Figure 2 (b) shows the measured n as a function of the dimensionless parameter ξ = ∆s/3ω z for trap frequencies of 1088 kHz, 501 kHz, and 217 kHz, while varying the intensities of both laser beams to achieve the desired values of the saturation parameter s. For the highest trap frequency, optimum cooling is achieved for ξ ∼ 1 in agreement with the prediction of eq. (10), which is indicated by a dashed red line. The measured mean phonon number is slightly above this limit (see also Fig. C1 in the appendix showing the same data on an expanded scale). Master equation simulations confirm that this discrepancy can be explained by the non-zero Lamb-Dicke parameter of η = 0.17. For the data taken at low trap frequency, the differences between measured and predicted mean phonon numbers become very pronounced, in particular at low values of ξ. For ω z = 2π × 217 kHz, the minimum motional energy ( n = 7.7(1.1) phonons) is obtained for an optical potential with a depth that is significantly higher than the optimum value ξ min = 5/6 predicted by eq. (10). This observation could be explained by electric field noise giving rise to motional heating at a rate that is much higher than the heating rate H intrinsic to polarization gradient cooling. Alternatively it could be due to a reduction of the cooling rate outside the LDR. In our experiment, we observed a heating rate of 1350(85) quanta/s at 217 kHz trap frequency. However, we infer from master equation simulations a cooling limit n > 5 phonons and estimate the cooling rate to be higher than 30000 quanta/s. Consequently, we conclude that electric field noise is not to blame and that the raised cooling limit is due to leaving the LDR.
Despite the increased cooling limit at very low trap frequencies, polarization gradient cooling is capable of cooling an ion below the Doppler limit (i.e. n = Γ 4ωz (1+α) for optimum cooling in one dimension) over a wide range of trap frequencies. Figure 2 (c) shows a measurement of n as a function of the trap frequency, which ranged from 127 kHz up to 1200 kHz, with optimum cooling ξ min achieved for ω z = 2π × 650 kHz. Here, sub-Doppler cooling is achieved over the entire frequency range.

Cooling of 1D ion chains
In this section, we present experimental results of PGC of linear ion strings trapped in a potential with transverse trapping frequencies of ω x = 2π ×2.67 MHz and ω y = 2π ×2.64 MHz. The axial trap frequency was set to ω z = 2π × 217 kHz except for experiments with 51 ions, where it was lowered to ω z = 2π × 127 kHz.
Spatially resolved sideband spectroscopy of a 22-ion string demonstrates the advantage of polarization gradient cooling over Doppler cooling as shown in Fig. 3 in a visually compelling way: a laser beam propagating along the symmetry axis of the ion trap is used to excite the |S 1/2 , m = 1/2 ↔ |D 5/2 , m = 3/2 transition including all first-order axial sidebands. The figure shows spectra obtained after Doppler cooling, (a), and with an additional polarization gradient cooling pulse of 1 ms duration, (b). Clearly, the excitation strength on the sidebands is strongly reduced by polarization gradient cooling. Moreover, an asymmetry between red and blue sidebands appears, however this does not indicate ground state cooling for collectively excited ion crystals. The spatially heterogeneous excitation probability on the sidebands reflects the different Lamb-Dicke parameters the ions have [38].
Another way of globally assessing the influence of polarization gradient cooling on the motional state of the ion crystal is to drive carrier oscillations. The carrier coupling strength Ω j of ion j depends on the phonon number n k of mode k as where n = (n 1 , . . . , n N ) is the vector of phonon numbers, L n (η 2 ) a Laguerre polynomial, η jk the Lamb-Dicke parameters, and Ω 0,j the bare Rabi frequencies, which might not be the same on all ions. The higher the motional energy is, the slower the Rabi oscillations will be and the faster they will dephase. This effect is clearly visible in the Rabi oscillations shown in Fig. 4 after Doppler cooling (a) and polarization gradient cooling (b). While a single Rabi flop is barely visible for Doppler-cooled ions, the oscillations persist for many periods after polarization gradient cooling. As η ∼ ω −1/2 , carrier Rabi oscillations mostly probe the occupation of low-frequency modes. The bigger Lamb-Dicke parameter of low-frequency modes also explains that ions at the ends of the string experience a stronger Rabi frequency reduction because the motion of these ions has a stronger overlap with low-frequency modes (see also Fig. 3). It should be noted that the  observed variation of the average Rabi frequency from ion to ion cannot be explained by effects such as the spatial profile of the laser beam or coupling strength variations by axial micromotion. A quantitative characterization of polarization gradient cooling is achieved by measuring the mean phonon numbers of individual collective modes by sideband spectroscopy. With the exception of the center-of-mass mode, it is very hard to extract mean phonon numbers from sideband measurements of collectively excited ions, such as in Fig. 3. We therefore use our single-ion addressing capability to transfer all ions except one to a different Zeeman state so that the axial laser beam at 729 nm can be used to drive sideband transitions coupling to a single ion only. For a given mode k of interest, the ion with the largest Lamb-Dicke parameter η jk is chosen in order to minimize the probe pulse length.  With the exception of the center-of-mass mode, all motional modes are cooled down to a few motional quanta. We attribute the rather high number of quanta in the center-of-mass mode to a competition of laser cooling with motional heating by electric field noise. Given the measured heating rate of 1350(85) quanta/s for a single ion, we expect the heating rate by electric field noise to increase by a factor of N , leading to about 11000 quanta/s for an 8-ion and 30000 quanta/s for a 22-ion string. Moreover, the measurement might slightly overestimate the cooling limit due to motional heating occurring within the time period (typically 250 µs) between the end of the laser cooling pulse and the start of the sideband probe pulse. While polarization gradient cooling cannot prepare the center-of-mass mode in very low phonon number states, it nevertheless substantially reduces the mode occupation after Doppler cooling that is estimated to be 158 (18), 184 (17) for ion strings with 8 and 22 ions, respectively.
Single-ion addressing is a technique that might not be available in every experiment with laser-cooled ion strings. For this reason, a semi-quantitative analysis of the cooling performance, that can be carried out with a laser exciting all ions equally, would be a useful tool. Towards this aim, we analyzed the excited state probabilities p j (t) of the j th ion when driving carrier Rabi oscillations for a thermal state n = ( n 1 , ..., n N ), which, through the dependence of Ω n j on the motional state, provide information about the phonon distribution p n of all motional modes involved. As the damping and frequency shift of Rabi oscillations is predominantly caused by thermally populated lowfrequency modes, one cannot hope to extract the mean phonon numbers of all motional modes from such a data set. Therefore, for the analysis of a polarization-gradient cooled ion string, the mean phonon number n i , (1 ≤ i ≤ N ) of the individual modes oscillating at frequency ω i are parametrized by a model with only three parameters, where n c is the mean phonon number of the center-of-mass mode, n 0 the lowest mean phonon number, and ω 0 the frequency at which n 0 is reached. The choice of this model is motivated by the fact that for a single ion the observed mean phonon numbers conform sufficiently well to the functional dependence on the trap frequency predicted by eq. (9), except for the center-of-mass mode which is affected by motional heating. Fitting to the experimental data is computationally expensive as expression (12) can be neither analytically calculated nor exactly evaluated numerically due to the enormous state space of the N harmonic oscillators involved. For a given thermal state n , we therefore have to sample from the thermal distributions many times in order to arrive at a numerical estimate of p j (t). This presents an obstacle to using standard minimization algorithms in the data fitting routine as the fit function now has become a stochastical variable. To overcome this problem, we used a dividing rectangles algorithm, developed in the context of variational quantum simulation (see [39] and references therein), for finding the parameter set (n c ,n 0 ,ω 0 ) that optimally fits a given data set. Figure 5 shows mean phonon numbers obtained by fitting Rabi oscillations of polarization-gradient cooled ion strings of 8, 22, and 51 ions. When fitting the data, we take into account that axial micromotion caused by an axial q-parameter of q z ≈ 0.0013 reduces the bare Rabi frequencies Ω 0,j (c.f. eq.(11)) of the outermost ions with respect to the ones of the innermost ions by about 2% (11%) for an ion string with 22 (51) ions.
For the 8-ion Rabi oscillations, the extracted mean phonon numbers are in reasonably good agreement with the single-ion sideband spectroscopy measurements. For the 22ion Rabi oscillations, the extracted mean phonon number of the center-of-mass mode and another low-frequency mode also agrees with the sideband measurements, whereas carrier fitting seems to overestimate the populations in the high-frequency modes, whose thermal occupation only weakly affects the shape of the carrier oscillations. At the time the measurement was carried out, we did not have the ability to carry out sideband spectroscopy with 51 ions. The 51-ion carrier oscillations indicate a strongly occupied center-of-mass mode, which is probably caused by a very high heating rate, as creating a 51-ion linear string necessitated dropping the axial trap frequency to 127 kHz. For a Doppler-cooled string, we estimated a mean phonon number of 760(100) for the axial center-of-mass mode.

Cooling of planar ion crystals
Operating our linear trap with a higher axial confinement and an increased frequency splitting of the transverse normal modes enables the creation of planar crystals in a plane that is normal to the direction having the highest frequency. For center-of-mass mode frequencies of ω z = 2π × 438 kHz, ω x = 2π × 2.76 MHz and ω y = 2π × 2.51 MHz, we were able to create the 22-ion zig-zag crystals displayed in Fig. 6. This trap configuration opens up the possibility of cooling the in-plane modes of motion in an axial polarization gradient without suffering from micromotion-induced line broadening that those ions experience which are not located directly on axis. Fig. 6 shows a spectrum of the blue sideband transitions of the planar modes, obtained by polarization-gradient cooling the crystal for 1 ms and mapping motional state excitation to the D 5/2 state by a sideband pulse with a 729 nm laser beam propagating along the axial direction. Here, the 2D-plot shows the spatially resolved excitations of all the ions, which are numbered according to their horizontal position in the crystal. The average excitation p is displayed on top of the plot as a black line. We were able to identify sideband resonances (gray arrows) of 40 out of the 44 in-plane collective modes of motion. The measured sideband frequencies are mostly in good agreement with mode frequency calculations based on both pseudopotential theory and Floquet-Lyapunov theory [40] (for the frequencies of the in-plane modes, which we observe, the predicted values by the two approaches are pretty close). Moreover, we also see a good match between ions with the highest D 5/2 state excitation and the projection of their normal mode vector onto the axial direction. Two example sets of normal mode vectors (white arrows) are graphically shown on top of images of the ion crystal.
A comparison of Rabi oscillations for a Doppler-cooled and a polarization gradientcooled 2D-crystal is shown in Fig. 7. Due to poor Doppler cooling results because of the micromotion-broadened cooling transition, no Rabi oscillations are visible in panel (a). In contrast, persistent Rabi oscillations are observed for polarization-gradient cooled ions. Fitting these oscillations with the model of eqs. (12) and (13), we estimate the center-of-mass mode's mean phonon number to be about 15 and all the other modes to be populated by six or fewer phonons on average.

Discussion and conclusion
We have investigated polarization gradient cooling of long ion strings and planar crystals in a linear Paul trap. Tests with a single ion show that a simple cooling model for a moving polarization gradient accurately predicts the cooling limit and cooling rate for experiments carried out deep in the LDR. At low frequencies, when the Lamb-Dicke parameter becomes big, we observe increased cooling limits and a shift of the optimum cooling conditions towards higher intensities of the lasers realizing the polarization gradient.
With long ion strings, we observe multi-mode cooling over a wide range of frequencies, with mean phonon numbers significantly below the ones that can be achieved by Doppler cooling on the same transition. Towards the low-frequency end, and in particular for the center-of-mass frequency, the cooling limit predicted by the simple model is not reached; in addition, it appears that these modes are optimally cooled for saturation intensities that are higher than predicted. The reason for this behaviour is currently not well understood. As η ∼ N −1/2 , the lowest-frequency mode has a much lower Lamb-Dicke parameter than a single ion would have. The presence of spectator modes might lead to an increased cooling limit and possibly to a reduced cooling rate; as a consequence, electric field noise might impede cooling of the center-of-mass mode.
For the planar crystals that we investigated, nearly all of the 2N in-plane modes of the crystal had ions with sufficiently big Lamb-Dicke parameters along the trap axis to enable polarization-gradient cooling of these modes. This makes polarization gradient cooling a valuable tool for cooling many modes of the crystal below the Doppler limit. Achieving sub-Doppler cooling at high rates is another asset of polarization gradient cooling. For Doppler cooling, the rate at which phonons can be extracted scales as W Doppler ∼ η 2 ω and therefore does not depend on the trap frequency ω (because η ∼ ω −1/2 ). For polarization gradient cooling, we have W P GC ∼ η 2 ω(Γ/∆)ξ 2 . When the depth of the optical potential is chosen to yield the lowest mean phonon number, i.e. at ξ = 5/6, the cooling rate achieved by PGC is expected to be smaller than the Doppler cooling rate if Γ ∆. On the other hand, one can achieve higher cooling rates for PGC than with Doppler cooling when ξ > 5/6 but at the expense of somewhat higher mean phonon numbers [24].
Another interesting application of polarization gradient cooling is sub-Doppler multi-mode cooling over a large frequency range. Polarization gradient cooling cannot compete with EIT cooling when it comes to cooling multiple modes close to the ground state [41,42,43]. EIT cooling, however, only achieves high cooling rates over a comparatively small frequency range due to the narrow dressed state facilitating the cooling. Moreover, since the highest EIT cooling rate scales as W EIT ∼ ω [44], EIT cooling becomes slow for low-frequency modes in contrast to polarization gradient cooling, for which the cooling rate can stay high at low confinement, where the statedependence of the potential seen by the ion becomes more pronounced.
The experimental setup used for the measurements described in this paper is dedicated to carrying out quantum simulations with long ion strings via engineered spin-spin couplings, which are mediated by transverse motional modes cooled close to the ground state. While the axial modes of motion do not need to be cooled to very low quantum numbers, hot axial modes can nevertheless give rise to spurious couplings inducing decoherence in the spin-spin couplings, as well as to coupling strengths variations in addressed single-ion operations with a tightly focused beam. Earlier attempts failed to cool these modes using an EIT cooling setup, which was used for ground-state cooling the transverse modes [42]. In contrast, polarization gradient cooling has been demonstrated to be able to prepare the axial modes at lower phonon numbers than the ones achievable by Doppler cooling. For long strings with a very large ratio between the lowest and the highest collective mode frequency ω, cooling range limitations imposed by n ∼ ω 0 /ω + ω/ω 0 might become a problem. Here, use of concatenated cooling pulses in a polarization gradient with an optical depth decreasing over time might extend the cooling range in a way similar to what has been demonstrated with sideband and EIT cooling [17,43].

Acknowledgements
We thank Haggai Landa for sharing with us his code for normal mode calculations using the Floquet-Lyapunov approach. CR acknowledges useful discussions with Alex Retzker. MKJ would like to acknowledge Pavel Hrmo for discussions and careful reading of the manuscript. The project leading to this application has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No 741541), and from the European Unions Horizon 2020 research and innovation programme under grant agreement No 817482. Furthermore, we acknowledge support by the Austrian Science Fund through the SFB BeyondC (F71) and funding by the Institut fr Quanteninformation GmbH.
H = H trap + H al consists of two parts, with H trap = ω z (a † a + 1 2 ) describing the ion's motion in the trap, where a † and a denote the motional raising and lowering operators, and the atom-light interaction Hamiltonian H al . In the frame rotating at the laser detuning ∆ from the atomic transition, it reads where σ + and σ − are the raising and lowering operators on the transitions driven by the two circularly polarized standing waves, each of which is described by two counterpropagating beams with wave vectors ±ke z . The symbol z denotes the position operator, which we express as kẑ = η(a + a † ). Furthermore, Ω str is the Rabi frequency of the stretched transition |S 1/2 , m j = +1/2 ↔ |P 3/2 , m j = +3/2 and φ is the laser phase at the ion position. We do not model the frequency shifts of the atomic states in a non-zero magnetic field, which we assume to be small in comparison to ∆. The Liouvillian operator L diss accounts for spontaneous decay during the laser-ion interaction. It is constructed using the jump operator describing the decay process from the required transitions, which reads as where J m = q p mq C m Γ(e −ikqz σ − m ). Here m stands for an index indicating transitions associated with 4 decay channels that exist between the Zeeman levels of the excited state and the ground state. We take into account the spatial dipole radiation pattern by subdividing each decay channel into three terms which transfer a momentum of k q with k q ∈ {−k, 0, +k} to the ion along its direction of oscillation [45]. The C m are the Clebsch-Gordan coefficients, which are 1/3 and 2/3 for ∆m j = 0 and ±1, respectively. The values of the pre-factors p mq are given in table A1.  We validate the predictions of the simple cooling model presented in section 2 by comparison with master equation simulations. For a numerical determination of the cooling rate and cooling limit, we time-evolve an ion initially prepared in the motional ground state under eq. (A.1) until the motional state reaches its steady state. Figure B1 shows the results for a Ca + ion oscillating at ω z = 2π × 1088 kHz for different phases φ of the polarization gradient. In order to achieve a comparison deep inside the LDR, we artificially increase the ion's transition wavelength by a factor of 10 (so that η = 0.017). As a consequence, the dynamics is slowed down a by a factor of 100. For determining the heating rate H(φ), we fit a tangent at t = 0 to the mean phonon number n(t) using a motional state space comprising the lowest m = 24 Fock states. For determining the steady-state phonon number and the cooling rate, we fit the data by a function f (t) = n m (1 − exp(−W m t)). For phases φ approaching π/4, the extracted steady-state phonon number and also the cooling rate depend on the state space cut-off m. Therefore, we carried out these simulations for m ∈ {4, 5 . . . 24} and used a finite-size scaling to extrapolate n m to infinity to obtain the true steady-state mean phonon number n ∞ . There are indications that this procedure slightly overestimates n ∞ . In Fig. B1(b), n ∞ (φ) is indicated by filled symbols whereas the open symbols represent n m=24 (φ). The cooling rate W (φ) is extracted by combining the results for H(φ) with n ∞ . In a second investigation shown in Fig. B2, we looked into the cooling limit in a moving polarisation gradient and its dependence on the relative frequency detuning δ between the two counter-propagating beams. The simple model predicts a constant cooling limit provided that δ is considerably bigger than the cooling rate. Again, we set ω z = 2π × 1088 kHz, but replace η → η/10 in order to be deep in the LDR. The data shown in Fig. B2 were obtained for optimum cooling conditions, ξ = 5/6 using a state space comprising the twenty lowest Fock states. After temporally averaging the steady-state n(t) over a period τ = 2π/δ, we find a good agreement between the master equation simulation and the simple model for sufficiently fast values of δ. If δ becomes smaller than the cooling rate, n goes up and the variance of n(t) (indicated by the error bars in the figure) strongly increases.

Appendix C. Comparison of master equation simulations with experimental data
The motional energy of a single ion shown as red data points in Fig. 2 of the main text are slightly higher than the prediction of the simple cooling model, which is shown in Fig. C1 as a pink line. Numerical simulations of the cooling limit (blue curve) agree well with the experimental data and show that this discrepancy can be explained by the fact that the ions are not deep in the LDR (η = 0.17). For small values of ξ, the simulation shows non-thermal phonon number distributions with tails that are more heavily occupied than in a thermal state. For this reason, we carried out the simulation once more for motional state spaces of different size including up to 24 Fock states and used a finite size scaling for extracting the mean phonon number.   Figure B2. Simulated mean phonon numbers (black dots) of a single ion as a function of frequency detuning δ between the two counter-propagating laser beams creating a moving polarisation gradient at the ion position. The error bar indicate the width of the distribution of mean phonon numbers that result from stopping the cooling pulse at a random phase of the moving gradient. The prediction of the simple cooling model is indicated by the red line. Again, for a comparison with the experiment, the horizontal axis should be multiplied by a factor of 100. Numer-Mast-Eq Simple-Model Exp Figure C1. Cooling limit as a function of ξ = (∆s/(3ω z )). Blue curve represents the simulation results with the master equation and magenta represents the results from the semi-classical treatment.