Bloch Oscillations Along a Synthetic Dimension of Atomic Trap States

Synthetic dimensions provide a powerful approach for simulating condensed matter physics in cold atoms and photonics, whereby a set of discrete degrees of freedom are coupled together and re-interpreted as lattice sites along an artificial spatial dimension. However, atomic experimental realisations have been limited so far by the number of artificial lattice sites that can be feasibly coupled along the synthetic dimension. Here, we experimentally realise for the first time a very long and controllable synthetic dimension of atomic harmonic trap states. To create this, we couple trap states by dynamically modulating the trapping potential of the atomic cloud with patterned light. By controlling the detuning between the frequency of the driving potential and the trapping frequency, we implement a controllable force in the synthetic dimension. This induces Bloch oscillations in which atoms move periodically up and down tens of atomic trap states. We experimentally observe the key characteristics of this behaviour in the real space dynamics of the cloud, and verify our observations with numerical simulations and semiclassical theory. This experiment provides an intuitive approach for the manipulation and control of highly-excited trap states, and sets the stage for the future exploration of topological physics in higher dimensions.

Synthetic dimensions provide a powerful approach for simulating condensed matter physics in cold atoms and photonics, whereby a set of discrete degrees of freedom are coupled together and re-interpreted as lattice sites along an artificial spatial dimension. However, atomic experimental realisations have been limited so far by the number of artificial lattice sites that can be feasibly coupled along the synthetic dimension. Here, we experimentally realise for the first time a very long and controllable synthetic dimension of atomic harmonic trap states. To create this, we couple trap states by dynamically modulating the trapping potential of the atomic cloud with patterned light. By controlling the detuning between the frequency of the driving potential and the trapping frequency, we implement a controllable force in the synthetic dimension. This induces Bloch oscillations in which atoms move periodically up and down tens of atomic trap states. We experimentally observe the key characteristics of this behaviour in the real space dynamics of the cloud, and verify our observations with numerical simulations and semiclassical theory. The Bloch oscillations thus act as a smoking gun signature of the synthetic dimension, and allow us to characterise the effective band structure. Our methods provide an efficient approach for the manipulation and control of highly-excited trap states, and set the stage for the future exploration of topological physics in higher dimensions through the use of a tunable artificial gauge field and finite-range interactions.

I. INTRODUCTION
Synthetic dimensions provide a powerful approach for simulating condensed matter physics in cold atoms 1-14 and photonics [16][17][18][19] , and they are opening up many new avenues for simulating and exploring exotic physics, including quantum Hall ladders 3,4,10 , non-Hermitian topological bands 22 , topological Anderson insulators 23 , and even lattice physics in four dimensions or higher 1, 11,24 .
A key reason that this framework is so powerful is that it is very general, and can be applied to a wide-range of very different physical systems. For example, synthetic dimensions have so far been realised experimentally in cold atoms with hyperfine 1,3,4 , magnetic 7 , Rydberg 8,9 , and clock states 5,6 , as well as with momentum 10,11 , orbital 14 and superradiant states 12 . However, in these experiments, the size of the synthetic dimension has been so far limited by the number of states that can feasibly be coupled. Indeed, the largest momentum state lattice 15 employed so far consists of a 1D lattice of 21 sites.
Notably, it has been recently realised that such limitations can be lifted if external degrees of freedom associated with trapping potentials are used to generate the synthetic dimensions 16,20,21 . Indeed, trapping potentials typically allow for tens or hundreds of trapped states in each direction, and by suitably coupling them one could implement very long synthetic dimensions, unleashing the full potential of this technique for quantum simulation. Additionally, this kind of synthetic dimension is extremely appealing because it provides a framework for the manipulation and control of trap states. A range of applications including quantum simulations in optical lattices [25][26][27] , trapped and guided atom interferometry [28][29][30] and quantum thermodynamics [31][32][33] require the use of highly-excited trapped states, which are generally difficult to realise with a good degree of precision and control.
In this work, we experimentally engineer a very long synthetic dimension of many tens of atomic trap states by dynamically modulating the harmonic trap of an ultracold atomic sample 20,21 . By controlling the driving frequency we generate a force along the synthetic dimension that induces Bloch oscillations, which act as a smoking gun signature that the synthetic dimension behaves as expected. Bloch oscillations were first famously predicted for electrons moving in a crystal under an electric field, and have since been observed in various setups, including optical lattices for cold atoms 34 , as well as synthetic dimensions of photonic frequency modes 35 and of room-temperate molecular angular momentum states 36 . However, Bloch oscillations in our experiment are physically very different from previous realizations, as they correspond to atoms periodically oscillating between lowand high-energy states of the harmonic trap. As such, another benefit of the synthetic dimension Bloch oscillations implemented here is that they allow us to explore highly-excited harmonic states, and thus can lead towards a novel approach for quantum engineering of external atomic states. More generally, this work paves the way for the exploration of higher-dimensional quantum Hall physics with artificial magnetic fields, and opens new opportunities in quantum simulations more widely. For example, a tunable artificial gauge field can be implemented by spatial modulation of the phase of the driving potential, allowing access to 2D quantum Hall physics. Furthermore, mean-field interactions between the atoms in real space should give rise to exotic interactions that are long-range and decay with distance along the synthetic dimension 20 , in contrast to the usual interactions in atomic gases and to the SU (N ) interactions in some other atomic synthetic dimension schemes. We therefore expect interesting ground state physics under the inclusion of interactions which will be of interest for future study.

II. OVERVIEW OF THE SCHEME
To introduce our scheme, let us consider an atomic cloud in a cigar-shaped harmonic trap with trap frequencies ω x = ω z ≫ ω y . In order to realise the synthetic dimension, we couple together the atomic trap states of the strong trapping potential along x with the spatiallyand temporally-varying driving potential given by: where V 0 is the driving amplitude, Θ(x) is the Heaviside step function, ω D is the driving frequency and φ is the initial driving phase. Physically, this corresponds to illuminating the upper half of the atomic cloud with constant power for the first half of the period, T D = 2π/ω D , before illuminating the lower half with the same constant power over the second half of the period (see Fig. 1(a)). This driving protocol is chosen because it is simple to implement, robust and it leads to a simple Floquet Hamiltonian with near-constant nearest-neighbour hoppings, corresponding to a textbook 1D tight binding model, as discussed below. Combining Eq. 1 with the 1D harmonic oscillator Hamiltonian along x gives the time-dependent Hamiltonian: written in the eigenstate basis of the "strong" trap along x, as indexed by λ = 0, 1, 2.... The stroboscopic dynamics of this system is captured by an effective timeindependent Floquet Hamiltonian, which we can approximate over a large number of trap-states by (see Appendix F): where ∆ ≡ ω x − ω D is the (small) drive detuning and J is a uniform hopping amplitude, which itself depends on V 0 and ω x and is calculated using Floquet theory (Appendix F). Note that this description is valid for nearresonant driving in a deep harmonic trap, i.e. such that ω x ≃ ω D ≫ ∆, J/ℏ. As depicted in Fig. 1(a), Eq. 3 describes a particle hopping between nearest-neighbour sites along a 1D tight-binding lattice with unit spacing in a synthetic dimension. The detuning plays the role of a constant force, F ≡ −ℏ∆, which therefore can induce Bloch oscillations. Note that the shaking phase φ appears in the effective Hamiltonian as a hopping phase, which we will exploit below to average over unwanted micromotion effects.
In the absence of a force along the synthetic dimension, i.e. when ∆ = 0, the effective model in Eq. 3 is translationally invariant along the synthetic dimension and has a single energy band in the Brillouin zone. Applying a nonzero force (i.e. ∆ ̸ = 0) accelerates a semiclassical wavepacket formed in the synthetic dimension bulk, such that it undergoes Bloch oscillations across the Brillouin zone, with a center-of-mass (COM) position along the synthetic dimension, λ com , that varies as 20 : from the initial position λ com (t = 0) = λ 0 com and where f B is the Bloch oscillation frequency. As we can set the spacing between the fictitious synthetic lattice sites equal to one, the periodic Brillouin zone covers [−π, π]. The Bloch oscillation frequency is then set by the magnitude of the applied force divided by the length of this Brillouin zone, i.e. f B = |∆|/2π, and so is entirely controlled by the detuning. Conversely, the amplitude of the Bloch oscillations is proportional to the bandwidth divided by the force, i.e. 4J/ℏ∆, and so depends on the detuning but also, through J, on the trap frequency and shaking power (Appendix F). The Bloch oscillations therefore provide a way to transport atoms between different trap states, with independent control over both the timescale and number of trap states explored. Experimentally, we use a thermal cloud of 87 Rb atoms in a harmonic trap with trapping frequencies ω x = ω z ≃ 2π × 160 Hz and ω y ≃ 2π × 10 Hz. The cloud was measured, both in-situ and with standard time-of-flight techniques, to have an initial temperature of T ≃ 20 nK. We used a rapid evaporation ramp to prevent the sample from condensing at this temperature, thus reducing the effect of mean-field interactions 20 , which may complicate the dynamics and will be of interest in future investigations. To realise the driving potential of Eq. 1, we utilise a digital micromirror device (DMD) that allows us to dynamically and spatially control the intensity profile of a laser beam with wavelength 800 nm (see Appendix A). We verify that the trapping and driving potentials are aligned to within ≃ 1µm. Our driving potential was chosen because it was found to be the most effective and robust, in the sense of not being sensitive to misalignments and other imperfections. This is in addition to the favourable theoretical properties discussed above. We then perform absorption imaging of the atomic cloud in position space after a very short time-of-flight (TOF) expansion of t tof = 5 ms, chosen to increase the visibility of the dynamics. This is demonstrated in Fig. 1(b) for a detuning of 9.84 Hz, where we plot the real space cloud density as a function of time for an example Bloch oscillation, showing the (a): Schematic of the DMD pattern (red ), which shakes the harmonic trap (blue/purple) and couples together nearest-neighbour trap states with an approximately uniform hopping amplitude, J, in order to create the synthetic dimension. Bloch oscillations (green) can be driven along the synthetic dimension by applying a force along λ, corresponding to detuning the shaking frequency from the trap frequency. (b): Experimentally, the cloud is imaged after a short time-of-flight expansion, as demonstrated here for a detuning of ∆ = 9.84 × 2π Hz. The colour scale is the column density in arbitrary units. The full dynamics includes micromotion within each driving period, so the experimental data are averaged over several values of φ. This averaging procedure makes the cloud appear to widen, although the true cloud width is approximately constant. (c): The real-space center-of-mass position, xcom, of the atomic cloud is extracted, as shown here for the same data as in (b). The experimental data (orange) are fitted (black) with a function motivated by Bloch oscillations in the synthetic dimension (see Appendices B and D), which captures the real-space dynamics well. Also shown are numerical results (blue) with a suitable TOF correction, discussed in the Main Text. The micromotion within each driving period can be seen in the inset for an initial driving phase of φ = 0. To reduce these unwanted micromotion effects, we average both the experiment and numerics over several values of φ (e.g. φ = 0 and π/2) to obtain the results shown in the main panel. As can be seen numerically, the residual micromotion oscillations have a small amplitude, which can be further reduced by averaging over more initial driving phases. Panels (b) and (c) use experimental parameters of V0 = 4.16 nK, T = 20 nK, ωy = 10 × 2π Hz, and ωx = 166.5 × 2π Hz, where the latter is determined experimentally by shifting the oscillation frequency data to pass through (|∆|, f ) = (0, 0). Experimental errorbars are 1σ statistical errors.
cloud COM being displaced away from the origin. Note that the density is averaged over several values of φ in order to reduce the effect of micromotion. As such, the cloud appears to widen along x in time, although the cloud width is actually approximately constant. The initial temperature of T ≃ 20 nK, corresponds to λ 0 com ≈ 2 for an initial Maxwell-Boltzmann distribution of atomic energies. This means that, initially, the atoms start near one "edge" of the synthetic dimension (at λ = 0), and so atoms must move up the synthetic dimension, irrespec-tive of the sign of the detuning, i.e. the direction of the force. The cloud also does not have a Gaussian distribution with respect to the synthetic dimension as the above semiclassical theory implicitly assumes; nevertheless, as we shall show, the prediction in Eq. 4 works well once appropriate corrections are included (see Appendix G).  Fig. 1(c). Note that the numerical fit parameters also include error bars, but these are smaller than the datapoint size. We use the same parameters as Fig. 1. Experimental errorbars are 1σ statistical errors.

III. EXPERIMENTAL AND THEORETICAL RESULTS
Bloch oscillations in λ-space naturally translate into atomic motion along x in real space, as different harmonic oscillator eigenstates have different real-space profiles. The resulting motion can be seen in Fig. 1, where we report the measured dynamics of the real-space COM position under the action of the shaking potential, as a function of time. It is important to notice that the full dynamics also includes micromotion within each driving period 20 ; for the real-space COM, as shown numerically in the inset of Fig. 1(c), the micromotion corresponds to large and fast oscillations as the atoms slosh backwards and forwards in the trap, while the stroboscopic Bloch oscillations translate into variations in the envelope of the dynamics. We are not able to reliably achieve the high time resolution to observe the micromotion in experiment, so we apply an averaging procedure to remove it. This is because of drifts in experimental parameters on a timescale of hours, which would make the large number of measurements required to reconstruct fast dynamics impractical. To overcome this, we root-mean-square (RMS) average over different experimental runs with suitablychosen different initial driving phases, φ. The raw cloud images in Fig. 1(b) have themselves been averaged in this way. This has the effect of making the cloud appear to widen significantly over the oscillation, although in each single shot the width is approximately constant. Note that this averaging procedure slightly lowers the apparent amplitude of the motion. However, by reducing the micromotion effects, we can then clearly observe the realspace signatures of the synthetic-dimension Bloch oscillations, as reported in Fig. 1(c) for the same parameters as Fig. 1(b) (orange), where the experimental data is fitted by a function (black ) motivated by synthetic-dimension Bloch oscillations (see Appendices B and D): where we fit to find the amplitude A and frequency f . We introduce the additional fit parameters g and ϕ to capture some details of the experimental data (see Appendices B and D). As can be seen, this fit captures the behaviour of the data very well, with agreement between the experiment and numerical simulations (blue curve) of a noninteracting 2D thermal cloud (see Appendix C). Our numerical simulations use the time-dependent Hamiltonian to evolve an ensemble of states, each of which is a superposition over the eigenstates of the 2D trap with random phase factors to destroy the phase coherence. Physical observables such as the cloud density are found by averaging over these phases. To account for the fact that we do not measure the true position of the cloud due to the TOF expansion, we include an approximate TOF correction to the cloud centre-of-mass position in our numerical simulations. In particular, we use the simulated  3) is no longer a good description (see Appendix F), although the numerics and experiment still appear to exhibit Bloch-oscillation dynamics. This is shown further in (b), where we plot the amplitude fit parameter (inset) and the frequency fit parameter, for V0 = 11.96 nK, with ωx = 142.1 × 2π Hz, ∆ < 0, and other parameters as in Fig. 1 and 2. Despite the large shaking power, we still observe similar trends to that at low power [c.f. Fig 2], and still with agreement in panel (b) to the fB = |∆|/2π analytical result from the simple tight-binding model. Experimental errorbars are 1σ statistical errors. Note that all numerical datapoints include error bars, but these are smaller than the datapoints. The green error band on the analytics is calculated from the errors on J and other numerical parameters.
COM momentum p com to find the semiclassical cloud velocity, and then shift the simulated cloud COM at each timestep (Appendix D).
To further characterise our experimental results, we plot in Fig. 2(a) the values of the oscillation frequency obtained by fitting our data for different detunings. As can be seen, for both experiment and numerics, the frequency increases linearly with detuning, as expected from the analytical Bloch-oscillation frequency (green line) given by f B = |∆|/2π in both real and synthetic space [c.f. Eq. 4]. The trapping frequency is determined by shifting a linear fit to the measured oscillation frequencies to pass through (|∆|, f ) = (0, 0). This provides a straightforward way to measure the trapping frequency, but does mean that any systematic uncertainty would not be detected.
In Fig. 2(b) we show how the amplitude of the realspace motion depends on the detuning by plotting the amplitude fitting parameters. As can be seen, the experiment (orange/black) and numerics (blue) both clearly show the expected growth in the real-space amplitude as the detuning decreases and higher-excited atomic trap states are explored. To make a quantitative comparison with semiclassical Bloch oscillations (Eq. 4), we have derived an analytical expression (see Appendix D) that converts the expected Bloch oscillation amplitude from synthetic space to real space under appropriate assumptions, including a correction for the finite fraction of atoms participating in the dynamics, as discussed further below. The expression is based on the formula: which connects the real space cloud COM x com and width σ x to the synthetic space COM λ com under certain assumptions. This result is derived in Appendix E. This analytical prediction is plotted in Fig. 2(b) (green), with errors calculated from our numerical parameters. As can be seen, there is agreement between the experiment, numerics and the analytics, demonstrating that we have achieved good control of the synthetic-dimension Bloch oscillations. We can also independently increase the number of atomic trap states explored (i.e. the Bloch oscillation amplitude) while keeping the oscillation frequency constant, by increasing the shaking power, V 0 , and hence the hopping parameter J [c.f. Eq. 3]. The dependence of the real-space COM amplitude on V 0 is shown in Fig. 3(a) for a fixed detuning ∆ = 8.3 × 2π Hz, while the inset shows the variation of the hopping J, as calculated with Floquet theory (Appendix F). As can be seen, the amplitude in experiment (orange), numerics (blue) and analytics (green) all increase as the hopping increases [c.f. Eq. 4]. Note that the analytical result is only plotted for V 0 ≤ 4 nK, as at higher shaking powers, our simple analytical model (Eq. 3) breaks down (see Appendix F). Despite this, we still observe clear Bloch oscillation dynamics at high power. This is further demonstrated in Fig. 3(b), where we use a very strong potential of V 0 = 11.96 nK, and still observe the same qualitative trends (i.e. the amplitude decreasing with the detuning and the frequency being equal to the detuning) as in the low power regime. Finally, we can visualise the Bloch oscillations along the synthetic dimension more directly, by translating our real-space experimental measurements [from Fig. 1(c)] into λ-space, under suitable approximations (see Appendix G). These experimental results are plotted in Fig. 4(a) (orange), along with numerical simulations (blue curve). Note that the latter is not averaged over different initial driving phases; however, in synthetic space, the micromotion is only a fraction of a "lattice spacing" and becomes negligible as the trapping frequency ω x increases 20 . The extrapolated maximum value of λ reached by the oscillating part of the cloud as a function of ∆, demonstrating that we have a long, controllable synthetic dimension. As discussed in Appendix H, we convert the experimental fit parameters in Fig. 2(b) to synthetic space (orange/black ), correcting for cloud splitting. We also plot our numerical maximum λ values (blue) for φ = 0, again corrected for cloud splitting. Furthermore, we derive appropriate analytics (green) for comparison (see Appendix G). Error bars on experimental points were calculated by converting the real space errors to synthetic space, while the error bars on the analytics were found by propagating errors on J and other numerical parameters. We use experimental parameters of ωx = 166.5 × 2π Hz, V0 = 4.16 nK, T = 20 nK, and ωy = 10 × 2π Hz.
We also compare our experimental results directly with the semiclassical analytical predictions [Eq. 4] with (green curve) and without (red curve) multiplying by a constant numerical re-scaling factor to account for the fraction of atoms contributing to the dynamics (Appendices G and H). This correction was also used in Fig. 2 and 3 and is necessary because the COM is skewed downwards by the thermal cloud splitting into two distinct parts, with approximately half of the atoms remaining in low λ-states during the oscillation, as can be seen in the numerical density distribution heatmap in Fig. 4(b). This effect is likely caused by small oscillations in the Floquet Hamiltonian matrix elements, and we discuss methods to reduce this splitting effect in Appendix I.
Importantly, we can also convert our experimental amplitude fit parameters ( Fig. 2(b)) to synthetic space and divide by the constant numerical cloud splitting factor in order to extrapolate how far up λ the oscillating part of the cloud is actually exploring (see Appendix H). This is plotted in Fig. 5 (orange/black ), where we can see that the fraction of the cloud involved in the Bloch oscillations is moving up several tens of Bloch states, demonstrating that we are creating a very long synthetic dimension. This is supported by our numerical (blue) and analytical results (green) (see Appendices C and D). Note that the number of sites along the synthetic dimension will eventually be limited experimentally by anharmonic-ities in the trapping potential 20 ; however, in this experiment, this does not play an important role. We included the appropriate quartic anharmonic terms in our numerical simulations and found that they had no effect in the range of trap states that are relevant here. We discuss the properties of the excited states we create, and methods to improve their fidelity, in Appendix I. Finally, we note that our methods can be extended towards singlesite resolution to on the order of or below ℏω x /k B . That would start the dynamics with only the λ = 0 state populated. The subsequent dynamics would maintain the single λ-state character. This point is discussed further in Apprendix I.

IV. CONCLUSIONS
Our experimental results show that we have engineered a synthetic dimension of atomic trap states by projecting a time-dependent shaking potential onto an atomic cloud via a digital micro-mirror device. Through control of the shaking potential's detuning, we induced Bloch oscillations along the synthetic dimension, observing the key characteristics of these dynamics in the real-space motion of the cloud. Our experiment demonstrates that a long and controllable synthetic dimension can be created. This opens up the way towards the exploration of topological physics using a synthetic dimension of harmonic trap states 20,21 by introducing a controllable artificial gauge field using a spatially-varying shaking phase. The spatio-temporal control of the shaking potential can also allow for future investigations of phenomena such as magnetic barriers, as well as the controlled population of excited atomic trap states, including direct imaging of the states 37 and single-site resolution of our current methods (Appendix I). Moreover, the addition of mean-field interactions in the cloud should lead to exotic interactions along the synthetic dimension 20 and, in turn, interesting ground state physics.

APPENDIX A: EXPERIMENT
In our experimental sequence we load 87 Rb atoms from a 3d MOT into a crossed optical dipole trap and then perform forced evaporative cooling 38 . The final trapping frequencies are f x = f z ≃ 160 Hz, and f y ≃ 10 Hz, where z is the vertical axis, resulting in a cloud elongated along y, with N ≃ 2 × 10 4 atoms at ≃20 nK. We conclude a posteriori that the degeneracy of the x and z trapping directions does not affect the dynamics because of the good agreement between theory and experimental data. We also verified using horizontal imaging that there are no significant dynamics along the z-direction. The optical setup to realise the dynamical potential and high resolution imaging is described in detail in Ref. 39. In brief, the light produced by a 800 nm laser is reflected by a DMD and then sent onto the atoms along the vertical direction, using an optical setup that produces a demagnification of a factor 100. The DMD is an 2d array of 1920 × 1080 micromirrors, each with size 10.8 µm. Each micromirror can be individually tilted every 100 µs, allowing us to produce dynamical optical potentials. The atoms are imaged on a CCD mounted in the vertical direction using a 20× magnifying system with a resolution of ≃ 2 µm. The numerical aperture of our imaging system is 0.28. We verify the alignment between the driving and trapping potentials to a precision of ≃1 µm by imaging both the cloud and the DMD pattern at once. This was done every 20 runs to avoid slow drifts. We also note that we do not observe significant heating of the cloud due to the driving potential. We determine this by observing the width of the cloud in the weak trapping direction after a short TOF expansion and finding that it does not change, as shown by experimental data in Fig. 16. This is caused by any "heated" atom spilling out of the weak trap, leading to some atom loss but a fixed temperature.

APPENDIX B: REAL SPACE EXPERIMENTAL DATA ANALYSIS
To reduce micromotion effects, the experimental data for the real-space COM position is RMS-averaged over initial driving phases drawn randomly from 2πn/30, with n = 0, 1, ...30. (The effects of RMS-averaging are discussed further in Appendix D) The resulting data are then fitted to the function: where A is the amplitude, f is the frequency, g is a damping factor and ϕ is a phase offset which accounts for random variation in the state preparation. The functional form of this fitting function is motivated by translating the semiclassical prediction for synthetic-dimension Bloch oscillations (Eq. 4) into real space. As shown in Appendix E, under certain approximations, this conversion can be achieved by taking: where x com and σ x are, respectively, the COM and width of the cloud (in harmonic oscillator lengths) with respect to the real position coordinate, x. Note that in choosing the form of the fitting function, we assume that the cloud width σ x is approximately constant as a function of time, as has also been verified numerically (Appendix D). The fitting parameters f and A are then plotted, for example, in Fig. 2(a) and (b) respectively.

APPENDIX C: DETAILS OF NUMERICAL SIMULATIONS
As shown in the main text, we numerically simulate the motion of a thermal cloud in two dimensions under the time-dependent HamiltonianĤ(t) (Eq. 2 in the main text, with the y-terms restored). In so doing, we choose to work in the λ − y basis; this avoids discretising the x-direction and hence reduces the size of the matrix representing the Hamiltonian in the numerics. From our numerical simulation of the wave-function, we then calculate the cloud density ρ(λ, y, t), and convert this into real space to yield ρ(x, y, t). The time evolution of the wave-function is done by numerically time-evolving an appropriate initial state as: where dt is a sufficiently small timestep. In order to simulate the dynamics of a non-interacting thermal cloud (i.e. a non-interacting gas which is distributed over the levels of the trap according to a classical Boltzmann distribution), we choose an initial state of the form 41,42 : where A is a normalisation factor, p i = exp (−βE i )/Z is the Boltzmann weight for the i th eigenstate of the 2D harmonic trap |ϕ i ⟩ with energy E i , β is the inverse temperature, Z is the partition function for the 2D harmonic trap and θ i is a random phase drawn from a uniform distribution between 0 and 2π. Note that a finite number, N , of harmonic trap states is used in this construction; in order to safely neglect higher-energy trap states, we check numerically that the Boltzmann weight has decayed to a sufficiently small value.
In the numerics, we then sequentially generate N ′ such random phase states, each with a different set of random phases. Each state is then time-evolved and the resultant densities are averaged together. By averaging over random phase factors, we destroy the phase coherence of the state; omitting this step would correspond to selecting a particular (and likely unphysical) coherent superposition of trap states. To further illustrate the importance of the random-phase averaging, we show how this can reproduce the expected density for a thermal cloud at t = 0. The latter can be found from the system's density matrix ρ = i p i |ϕ i ⟩ ⟨ϕ i | , as: where ϕ i (r) = ⟨r|ϕ i ⟩ 43 . We can also calculate the density of the random phase state as: which, as can be seen, involves a double sum over the harmonic trap states. However, averaging over the random phases gives: as desired, where we used the identity: The random phase state therefore reproduces the density for a thermal cloud at t = 0 under suitable averaging. We can also demonstrate the effects of phase-averaging numerically as shown in Figure 6, where this technique is applied to a 1D harmonic trap for (a) only five random phase states and (b) for 100 random phase states. In both cases, the blue curve is the density calculated via the above method and the orange curve is the expected thermal cloud density calculated explicitly as: where A is a normalisation constant and σ = k B T /mω 2 x is the cloud width, controlled by the trap frequency and temperature 43 . For large enough numbers of random phases states included in the average (Fig. 6(b)), we see good agreement with the expected thermal cloud density. Note that in so doing, we have assumed that our cloud is non-interacting, because if interactions are present, we can no longer use a Boltzmann-weighted superposition over single-particle trap states. Throughout this work, we use 50 random phase states in our numerics.

APPENDIX D: DETAILS OF NUMERICAL DATA ANALYSIS
Here we describe the data analysis steps carried out on the simulated cloud density (Appendix C) in order to extract the Bloch oscillation frequency and amplitude, which we have then compared to experiment and analytical results in the main text.
Firstly, the cloud center-of-mass (COM) and width are calculated from the real space cloud density ρ(x, y, t), found using the method in Appendix C. An example is shown in Fig. 7, where we see oscillations in both the synthetic and real space COM ((a) and (b) respectively), including high-frequency micromotion, as discussed in the main text. Note that the λ-space COM is calculated as λ com = λ λρ(λ, t), where ρ(λ, t) is the probability density with respect to the synthetic dimension, calculated numerically as in Appendix C, with the y-dependence integrated out. We also see that, although the cloud does visit higher harmonic oscillator states, the cloud width (panel (c)) is approximately constant in time. We use this observation in Appendix G to simplify our analysis. Throughout this section, we use the typical experimental parameters: V 0 = 4.16 nK, ω x = 166.5 × 2π Hz, ∆ = 9.84 × 2π Hz, ω y = 10 × 2π Hz, T = 20 nK and the initial driving phase φ = 0 for any data including micromotion and φ = 0 and π/2 for any data where micromotion has been averaged out, as discussed further below. These are the same parameters as in Fig. 1, 2, 3 and 4 in the main text. As an aside, we mention the effect of different signs of detuning (i.e. different signs of the effective force) on the dynamics in synthetic space. In Fig. 8, we plot the synthetic space dynamics for both φ = 0 (blue and yellow ) and π/2 (red and purple) for opposite detuning signs. We see that, for φ = 0, there is a small FIG. 8. The effect on the synthetic space dynamics of different signs of the detuning for the two initial driving phases used for all the numerics. For φ = 0 (blue and yellow ), we see a small difference in oscillation amplitude of around 1 synthetic lattice site for different signs of ∆, whereas for φ = π/2 (red and purple) the amplitudes are unchanged. We use the same parameters as Fig. 7, but with |∆| = 5 × 2π Hz. amplitude difference of around 1 synthetic lattice site, whereas for φ = π/2 the two amplitudes are the same. This means that we expect our oscillation amplitudes to be similar regardless of the sign of ∆, as we find in our results in the main text. Note that the presence of a hard wall boundary at λ = 0 makes this different to the expected result for a wavepacket prepared in the synthetic dimension bulk. In that case, we expect the wavepacket to move in opposite directions for opposite signs of ∆.
Secondly, we need to take into account the time-offlight expansion (TOF) carried out in the experiment, as this is not included in the numerical method described in Appendix C. If this is not accounted for, then the experimental COM oscillations will be significantly larger than in the simulation. To correct for this in the numerics, we use the simulated momentum distribution to calculate the COM momentum in real space, p com , and hence find the semiclassical cloud COM velocity v com = p com /m at each time-step, such that we can correct the COM position from the numerics as: where t tof is the TOF expansion time, corresponding to 5ms in this experiment. This method applied to our example oscillation (Fig. 7) is shown in Fig. 9(a). We see the same qualitative form as Fig. 7(b), but with an amplitude around five times larger. Thirdly, as discussed above and in the main text, we need to remove the micromotion before we can extract the amplitude and frequency of the Bloch oscillation. Experimentally, this is done by repeating the experiment for multiple different starting phases φ of the driving potential, as discussed in Appendix B. The micromotion is approximately removed by taking the root-mean-square (RMS) average over the M chosen driving phases: wherex (i) com is the COM for the i-th initial driving phase. Physically, the choice of φ controls the phase of the micromotion oscillations. This approach is applied directly to the experimental data, as well as to the simulated data after the TOF correction.
If the driving phases are chosen randomly or if the micromotion is very complicated, then we expect to take the large M limit in Eq. 17 and average over many experimental or numerical runs. However, if the micromotion were described by a perfectly sinusoidal function, such as e.g. f (t) = sin ω D t, then we would in fact only need to average over any two values of φ that are separated by (2k + 1)π/2, with k = 0, 1, 2..., using the property that (f (t)) 2 + (f (t + (2k + 1)π/2)) 2 = 1.
In practice, our numerical simulations ( Fig. 9(a) and 7) suggest that the micromotion is close to being sinusoidal and so we try averaging over only one pair of phases related by (2k + 1)π/2. Indeed, Fig. 9(a) shows an example of this RMS-averaged oscillation (red curve) over two phases (φ = 0, π/2), together with the φ = 0 unaveraged data (blue curve). As can be seen, averaging over only two runs is already sufficient to remove the majority of the micromotion, with a small residual that could be removed by using more pairs of phases. However, the amplitude of the averaged curve is smaller than the un-averaged result, and we will return to this point later when discussing corrections to the analytical results in Appendix G. Finally, to extract the oscillation ampli- tude and frequency, we fit the function: to the real space COM in both simulation and experiment, as introduced in the main text. This functional form is motivated by the analytical results; if we convert the expected oscillation in synthetic space (Eq. 4 in the main text) to real space (using Eq. 6 in the main text which is derived in Appendix E below), we obtain: (1 − cos(∆t)).
FIG. 10. Example of the conversion formula (Eq. 29, red ) applied to the same numerical simulation as Fig. 7 (blue), showing agreement up to a small offset. The formula is applied using the RMS time-averaged data in Fig. 7(b) and (c).
Since, in our numerics, x com (t = 0) = 0, we have λ 0 com − σ 2 x + 1/2 = 0, where we assume that σ x is constant in time as justified previously. This then suggests the functional form of the above fitting function (Eq. 18). By hand, we then add in the exponential damping factor, to capture wavepacket splitting effects (see below) and other sources of damping in experiment (e.g. heating), as well as the phase ϕ to account for random variation in the position of the cloud after state preparation in the experiment. (Note that for the numerics we have defined x(t = 0) = 0, and so we set ϕ = 0 to reduce the number of fitting parameters.) Figure 9(b) shows an example fit (red curve) using Eq. 18 applied to our numerical data (blue curve) from Fig. 9(a). The error bars on the fit parameters are from the least-squares fit. The fitted amplitude A and frequency f are plotted in the main text for varying detuning and shaking power.

APPENDIX E: MAPPING FROM THE SYNTHETIC DIMENSION TO REAL SPACE
Analytical results for our synthetic dimension scheme are expressed in terms of the synthetic dimension λ, but the experiment naturally probes real space x. We therefore now derive a formula linking these, which is given as Eq. 6 in the Main Text, and which is used to justify the fitting function as discussed above, and to process the analytical results in Appendix G. Our aim is to link the COM of a state |ψ⟩ with respect to the synthetic dimension, λ com , to the COM and width of the state in real space, x com and σ x respectively. To start, we will expand our state in the harmonic-trap basis as: where c λ are complex coefficients with λ |c λ | 2 = 1. In terms of the COM variables, it is straightforward to show that where: are the COM in real and momentum space respectively.
Here, x com is measured in units of ℏ/mω x , and p com in units of √ ℏmω x . We can also write which follows from the usual expression for the variance of a random variable. Now writingx = 1 √ 2 (â+â † ), wherê a † andâ are the usual harmonic oscillator raising and lowering operators, and substituting into Eq. 24, yields: where we identified λ com = λ |c λ | 2 λ and where: Finally, re-arranging for λ com gives: To gain some intuition for the S terms, consider the limit of preparing the system in a particular eigenstate of the harmonic trap |λ 0 ⟩, so c λ = δ λ,λ0 . This makes S = 0, and we have: If, instead, the state is a semiclassical Gaussian wavepacket with c λ ∼ exp(−(λ − λ 0 ) 2 /2σ 2 λ ), we expect Eq. 29 to hold approximately when σ λ is sufficiently small. Figure 10 shows an example of this result (red curve) compared against the COM calculated directly from a numerical simulation (blue curve), with agreement up to a small offset. The numerical curve displays small micromotion oscillations as expected, as does the curve calculated from our formula. The formula shows these oscillations because using only a single pair of phases during micromotion averaging does not perfectly remove FIG. 11. The momentum space conversion formula (Eq. 30) applied to the same oscillation as Fig. 10. The result of applying Eq. 30 to the RMS phase-averaged momentum COM and width (blue) is compared to the synthetic space COM calculated numerically (red ) and we see agreement up to a small offset.
the micromotion, as discussed in Appendix D. Note also that the example data and parameters used in this section are the same as the previous one, although the level of agreement is similar in all cases studied. This result demonstrates that the derived conversion formula still holds in the case of the thermal cloud. We can repeat the above calculation but working in terms of momentum rather than position to find: where σ p is the width of the state in momentum space, and the momenta are measured in units of √ ℏmω x . An example of this formula applied to the oscillation in Fig. 10 is shown in Fig. 11, showing similar agreement to the real space result.

APPENDIX F: EFFECTIVE TIME-INDEPENDENT HAMILTONIAN DESCRIPTION
In this section, we show how we arrive at the effective time-independent Hamiltonian shown in the main text: (31) starting from the full time-dependent Hamiltonian (Eq. 2 in the main text). Note that this effective Hamiltonian is different to that of 20 in the sense that we can assume that our nearest-neighbour hoppings are constant in λ, whereas those of 20 scale as √ λ; this is due to the choice of driving potential. As the Hamiltonian is periodic in time, we can use Floquet theory to define the Floquet Hamiltonian,Ĥ F , according to: whereÛ (T D ; 0) is the time-evolution operator over a full period of the driving, T D = 2π/ω D , where ω D is the driving frequency 44 . This Hamiltonian can be calculated numerically by splitting the time-evolution operator (Eq. 32) into sufficiently many small timesteps dt.
For simplicity, we have neglected terms in the Hamiltonian that depend on the y -direction to work with a 1D Hamiltonian in the λ -basis. We first investigate the matrix elements of Eq. 32 for the low-V 0 case (V 0 = 4.16 nK), as considered in Fig.  1, 2, 3 and 4 in the main text. The Floquet Hamiltonian matrix elements are plotted in Fig. 12 for several detunings and for parameters listed in the caption. More precisely, we plot the real (top row) and imaginary (middle row) parts of the first five diagonals ofĤ F , as well as |Ĥ F | as a heat map to show the longer-range structure of the Hamiltonian (bottom row). For all figures in this section, the on-site terms are in blue, the nearest-neighbour (NN) hoppings are in red, and longer-range hoppings in other colours. For this low-power figure, we have applied a constant offset to the time-dependent Hamiltonian to ensure that the on-site matrix elements are zero for λ = 0. This shifts the location of the instability regions (see below) and allows the behaviour for ∆ < 0 to be seen more clearly.
As can be seen, for many values of λ, the Hamiltonian can be approximated by the form used in the main text. Firstly, the detuning induces a tilt on the on-site matrix elements (blue curve) ∼ ∆λ, which allows us to interpret this detuning as a constant force along the synthetic dimension (see e.g. Fig. 12 columns (b), (c) between λ ≈ 0 and λ ≈ 20). Note that we verified that the slope of the on-site terms is equal to the detuning by fitting a straight line to these plots. Secondly, across the same regions, we also find nearly-flat NN hoppings (red curve), such that the NN hopping energy J can be calculated by taking the average of the NN matrix elements up to the onset of the instability region (see below). This hopping energy does not vary significantly as a function of detuning, and we find the value J/ℏ = 106 ± 8Hz for V 0 = 4.16nK. To calculate this error bar on J, we add the standard deviations of the NN hoppings for each of the N ∆ detunings in quadrature to obtain σ, and then use the error bar σ/ √ N ∆ . Thirdly, we also see that we have a small longrange hopping (purple curve) which we neglect. Note that all the matrix elements show small oscillations with respect to λ, which may act as a potential minimum around λ = 0 and be the cause of the cloud splitting effects discussed in the main text and below. This is discussed further in Appendix I. Finally, we see the "odd" diagonals of the Hamiltonian (nearest-neighbour, NNNN, etc.) are purely imaginary and the others are purely real, which is caused by the initial driving phase φ playing the role of a Peierls phase in the effective model.  Fig. 1, 2, 3 and 4 in the main text. The top row shows the real part of the first five diagonals, the middle row shows the imaginary part, and the final row shows a heat map of |ĤF |. The on-site terms are in blue, the NN hoppings in red and longer-range hoppings in other colours. We see that we can approximateĤF by a nearest-neighbour tight-binding model where ∆ plays the role of a force along the synthetic dimension. The onset of instability regions (IR), marked by the black dotted lines, show long-range hoppings for some λ and are a numerical artifact and are not physical.
For non-zero detuning in Fig. 12, there also appear to be regions of λ where we no longer have linear on-site terms and flat NN hoppings, but instead have significant long-range hoppings (e.g. between λ ≈ 20 and λ ≈ 60 in column (b)). However, these are an artifact of our numerics and should not be physical. These regions arise because the Floquet Hamiltonian is not unique, but depends on the branch of the matrix logarithm (Eq. 32). In the numerics, the principal branch is always taken, meaning that the Floquet Hamiltonian is constructed to have eigenvalues that only lie between −π/T D and π/T D . This leads to the apparent"wrapping around" of onsite terms and an associated variation in off-diagonal terms when the on-site shift due to the detuning becomes large. This can be seen by noting that for larger detuning (column (d)), the apparent breakdown happens earlier and more frequently, and between these regions, the matrix elements look regular and well-behaved. We have also checked our interpretation numerically by adding a constant offset to the time-dependent HamiltonianĤ(t) that changes the size and location of the apparent breakdown regions while leaving the on-site slope and NN hoppings otherwise unchanged. Finally, we also do not observe any qualitative change in behaviour of our numerical simula-tions (Appendix C) when the cloud moves in the instability regions, further confirming that these are not a physical effect.
It is important to distinguish between the above numerical artifact and a genuine breakdown of the effective Hamiltonian employed in the low-V 0 case, caused by the failure of the rotating-wave approximation. This is apparent e.g. in Fig. 12 (a) above λ ≈ 55, where we observe that the matrix elements begin to deviate from their previous values.
For larger values of V 0 (such as considered in Fig. 2 in the main text), we also observe numerically that the matrix elements become less uniform and more long-ranged. This is shown in Figure 13, where we plot the Floquet Hamiltonian matrix elements for V 0 = 11.96 nK for two detunings, with parameters as listed in the figure caption. As can be seen, in this case, we find significant long-range hoppings for all λ, and a non-zero detuning does not simply add a slope to the on-site terms. This therefore precludes building a simple analytical model for this behaviour, although we still find that our numerical simulations agree well with experimental results, as in Fig. 3 in the main text.
FIG. 13. Floquet Hamiltonian matrix elements for a higher shaking power than Fig. 12. We have ωx = 2π × 142.1 Hz, V0 = 11.96 nK and φ = 0, and ∆ = 2π × 0 and −2 Hz in columns (a) and (b) respectively. These parameters correspond to the high-power data shown in Fig. 3(a) in the main text. Unlike the low-power case, we find significant long-range hoppings for all λ, and that a detuning does not simply induce a slope on the on-site terms. This is not a numerical artifact, but a result of the failure of the rotating-wave approximation, and precludes building a simple effective Hamiltonian. We use the same layout and colour schemes as Fig. 12.

APPENDIX G: DETAILS OF ANALYTICAL RESULTS
Here we describe in detail how the analytical results for Bloch oscillations in both real and synthetic space are obtained, including corrections to make them comparable to numerical and experimental data, as plotted and discussed in the main text.
As stated in the main text, we expect the cloud COM to oscillate with respect to the synthetic dimension as: (1 − cos(∆t)).
We therefore expect an oscillation frequency of ∆/2π and a maximum λ of λ max = λ 0 com + 4J/ℏ∆. We can then use our result connecting the real and synthetic dimension, Eq. 29, to calculate the maximum displacement of the cloud from x = 0. This gives: as discussed in the main text, and where we measure x max and σ x in units of ℏ/mω x . We can analytically calculate the width for a thermal cloud σ x = k B T /mω 2 x 43 , and approximate the cloud width as constant in time, as justified in Appendix D.
We now take into account cloud splitting, as described in the main text. Note that we discuss methods to reduce this effect in Appendix I. The COM with respect to the synthetic dimension is calculated as: where ρ(λ, t) is the synthetic space density, calculated numerically as in Section I, where we have integrated out the y-dependence. The presence of a split wavepacket component skews this average downward, and the analytical result (Eq. 33) therefore overestimates numerical and experimental amplitudes.
To correct for this, we define a cutoff in λ, λ c , that cleanly separates the two wavepacket components near the peak of the oscillation. In particular, we choose λ c = 2J/ℏ∆, because this is about half of the maximum λ coordinate at the oscillation peak (Eq. 33). Near the oscillation peak, we can then write: where λ < com (λ > com ) is the centre of mass of the lower (upper) wavepacket respectively, and r is the amount of wavepacket above the cutoff. We can calculate r numerically and find that it is approximately constant with respect to the detuning, with an average value of r = 0.52 ± 0.03, where the error bar is the standard deviation of the r values over the detunings. This value is for the low-power data, with V 0 = 4.16 nK, ω x = 166.5 × 2π Hz, ω y = 10 × 2π Hz and T = 20 nK. To make our analytical result for λ max comparable to the numerics and experiment, we hence correct it as: This result is plotted in Fig. 14 (red ) and compared against the maximum λ com in the numerics (blue). We see good agreement, up to a roughly constant offset of around two synthetic lattice sites. We can also use this corrected expression for λ max in our real space result (Eq. 34) to get: FIG. 14. Cloud-splitting corrected analytical synthetic space amplitude (red ) compared to the maximum λcom obtained numerically (blue). We see good agreement, up to a constant offset of around two synthetic lattice sites. Error bars on analytical results are calculated by propagating errors on numerical parameters. Lines are a guide to the eye. Here, we use parameters for the low-power data, as in, for example, Fig. 9.
We now need to apply three further corrections to this real space result in order to compare to the experimental and numerical results including TOF. Firstly, we need to estimate the effect of TOF expansion. To do this, we numerically calculate the ratio of the maximum of the real space oscillation including TOF (Appendix D) to the maximum of the oscillation without. We find that this is approximately independent of detuning, and has a value of α tof = 5.25±0.08, where the error bar is the standard deviation of the α tof values over the detunings. As for r above, this is calculated for the low-power data. We then correct our real space result as: Alternatively, we can consider correcting our result analytically by calculating the COM velocity v com at the oscillation peak. To do this, we can use the momentum space version of the λ formula (Eq. 30) to calculate p com at the oscillation peak, including the wavepacket splitting correction, and then calculating v com = p com /m, to yield: where σ p is calculated from the numerics and assumed constant in time, and we measure σ p in units of √ ℏmω x . We can then apply this time-of-flight correction to x max : The effect of doing this, together with the two corrections described below, is shown in Fig. 15  where the blue curve is the amplitude fit parameter for the numerical simulations, and the orange points are the fit parameter for the experimental data, as discussed in Appendix D and the main text, and where we use the same low-power parameters as Fig. 1, 2, 3 and 4 in the main text. We see that the level of agreement between these analytics and the numerics and experiment is comparable to the method using a numerical TOF correction, although this method clearly overestimates the effect of TOF. Next, we correct for the effect of RMS averaging. As described in Appendix D, the RMS average over initial driving phases to remove micromotion slightly underestimates the true size of the oscillation ( Fig. 9(a)). If we model the dynamics of the COM as a beat: where ω B is the slow Bloch oscillation frequency, we can then RMS-average out the micromotion using two initial driving phases of 0 and π/2: (A sin(ω D t + π/2) sin (ω B t/2)) 2 ). (44) Here, we have shifted the micromotion part of the beat by π/2. Simplifying gives ⟨x com ⟩ 2 = A 2 sin(ω B t/2) 2 /2. We can then calculate a scale factor between the averaged and un-averaged results, α rms ≡ ⟨x com ⟩ max /A = 1/ √ 2. We then apply the correction: Note that we can calculate a scale factor based on an infinite number of driving phases by doing the RMS timeaverage, and we find the same result. As discussed in Appendix D, our numerical results can be well-described by such a sinusoidal beat, so this approach works well for our data. Our final correction is due to the fitting function used (Eq. 18). Neglecting damping, our fitting function will have a maximum value of x max = √ 2A, and we analytically calculate x max . We therefore correct our analytical result so far as: so that we can compare directly to the fit parameters obtained from the numerics and experiment. Our final analytical result for x max , as shown in the main text, is then: where we have restored SI units. Finally, note that these same analytical results are used for the variable-power data in Fig. 3(b) in the main text, where the numerical parameters α tof , J and r were re-calculated for this dataset. To convert our real space experimental oscillation into synthetic space (Fig. 4(a)), we use our conversion formula (Eq. 6), assuming that the real space width σ x (measured in harmonic oscillator lengths) remains constant in time and set by the trap frequency and initial temperature. We then obtain: where x com (t) is the experimental real space COM measured in harmonic oscillator lengths. The conversion for the data in Fig. 5 is similar, but here we convert the experimental amplitude fitting parameters, so we use the conversion: This is then compared against the corresponding numerics, where we extract the maximum value of λ com and rescale by 1/r to account for the cloud splitting. We also compare against an analytical expression, λ com = λ 0 + 4J/ℏ∆, as shown in Fig. 5.

APPENDIX I: CREATION OF HIGHLY-EXCITED STATES
Bloch oscillations along the synthetic dimension offer a way to controllably populate highly-excited atomic trap states. Indeed, by stopping the shaking at a certain given time it is possible to 'freeze' the dynamics. If this is done when the Bloch oscillation is exploring high λ states, the resulting state will have a considerable admixture of highly-excited harmonic trap states. In Fig. 16 we show an example of states obtained after a 5ms TOF expansion with such a technique with ∆/2π ≃ 17 Hz, V 0 ≃ 1.2 nK and stopping the Bloch evolution when the oscillation has reached the peak (40 ms) and half way (20 ms). By averaging over several different phases of the shaking potential, we can reconstruct the whole set of states that can be obtained for a given Bloch evolution time, this is shown as shaded areas in Fig. 16 (a) and the column density profiles in Fig. 16 (b) and (c). We observe that the resulting density profiles acquire the characteristic double-lobe structure typical of highly-excited harmonic trap states. Additionally, as the Bloch dynamics evolve, we observe that the distance between the two lobes increases, as expected when populating increasingly higher states. Notably, by controlling the phase of the shaking potential, it is possible to accurately control the shape and position of the final state. This is shown in Fig. 16  FIG. 17. (a): Time evolution of the integrated synthetic space density for a thermal cloud under our usual driving potential using typical parameters, for a detuning of ∆ = 5 × 2π Hz. We see significant cloud splitting and a relatively wide oscillating component. We use identical parameters to our low-power data elsewhere in this work, with φ = 0. (b): The result in (a) but with increased trapping frequency and decreased temperature, showing that the oscillating part of the cloud becomes much narrower in synthetic space. Note the different colour scale to (a). We use the same parameters as (a) but with T = 5 nK and ωx = 2π × 500 Hz. (c): Time evolution of the integrated synthetic space density of a thermal cloud for a modified driving potential, showing that the entire cloud moves, with no split component. Note that the linear driving potential together with parameters as in (a) produces a wide oscillating cloud, so we also increase the trapping frequency and decrease the temperature here to compensate for this. We use: κx qho = 1 nK, with x qho = ℏ/mωx, ωx = 500 × 2π Hz, ∆ = 5 × 2π Hz, T = 10 nK and ωy = 10 × 2π Hz. The insets of each panel shows the corresponding λ COM as a function of time. Note that the results for (a) and (b) are not rescaled to take into account cloud splitting.
(a), where the two solid lines correspond to averages over several runs with the same phase both for the 20 and 40 ms cases.
We have additionally measured the lifetime of the states created by measuring the number of atoms as a function of time after the Bloch evolution is stopped. We then performed an exponential decay fit, whose decay time sets the lifetime of the state. Concerning the states shown in Fig. 16 in particular, we have measured that the lifetime for the states produced after 20 ms of shaking potential is ≃1 s, while for those produced after 40 ms is ≃600 ms. Therefore the lifetime of these highly-excited states is sufficiently long to allow one to practically use them. As an example, they would provide a good overlap with a double-well potential enabling new possibilities for trapped atom interferometry 30 .
The fidelity of the excited states could be further improved by decreasing the proportion of the cloud that remains in the low-lying λ states, and by decreasing the width with respect to λ of the portion that does oscillate. One approach to achieving this is by a combination of using a stronger trap (i.e. increasing ω x ) and/or a lower temperature, both of which reduce the width of the initial Maxwell-Boltzmann distribution. Fig. 17(a) and (b) show the synthetic space density profile, with ydependence integrated out, for a numerical simulation of a thermal cloud with (a) the typical parameters used elsewhere in our work, and (b) a larger trapping frequency and lower temperature, over a single Bloch oscillation period. As can be seen, in panel (b) the part of the cloud that oscillates is narrower with respect to λ than with our more typical parameters in panel (a). Note that, for a small enough temperature, a significant initial condensate fraction may also affect the dynamics. However, we have verified numerically that an initial state with the whole cloud in the λ = 0 state (i.e. a non-interacting BEC) still undergoes Bloch oscillations under the driving potential.
We also note that the width of the oscillating part of the cloud also depends upon the detuning. In Fig. 18, for the low-power parameters used in the Main Text and a detuning of ∆/2π = 2Hz (panel (a)) and ∆/2π = 5Hz (panel (b)), the density with respect to λ is shown close to the peak of the oscillation. The oscillating part of the cloud can be seen to have a full-width at half-maximum (FWHM) of around 8 states in (a). This width depends upon the detuning, with larger detunings yielding smaller widths. For example, in panel (b), the FWHM is around 5 states. Note that these methods to reduce the width of the cloud with respect to λ could be extended towards single-site resolution.
Another strategy for improving the fidelity is to optimize the driving potential. As an example, Fig. 17(c) shows the effect of a different driving potential, V (x, t) = κx cos(ω D t), upon a lower-temperature thermal cloud (T = 10 nK), as calculated numerically using the technique described in Appendix C. Note that this driving potential is the same as that of 20 (see Eq. 2 therein). As can be seen, the entire cloud undergoes the Bloch oscillation, with no split component. This difference is likely because for the experimental driving potential (Eq. 1 in the Main Text), there are oscillations in the Floquet FIG. 18. Cloud density with respect to λ for the low power parameters used in the Main Text and a detuning of ∆/2π = 2Hz (panel (a)) and ∆/2π = 5Hz (panel (b)). The density is extracted from close to the peak of the oscillation. We clearly see the oscillating part of the cloud, which has a full-width at half-maximum of around 8 states in (a) and around 5 states in (b).
Hamiltonian matrix elements (e.g. Fig. 12) at low λ, which include, for example, an effective potential minimum at λ = 0 which may trap some atoms, leading to the cloud splitting. For the linear driving potential, on the other hand, the Floquet Hamiltonian does not exhibit such oscillations, although there is then a strong λ-dependence in the hopping elements 20 . Going further, quantum control approaches could be used to further numerically optimize the driving potential, in order to improve the fidelity of desired state preparation. Finally, another possible strategy for improving the fidelity of excited states would be to employ a two-step protocol. Firstly, the cloud could be placed into an excited trap state without using the driving potential 26 . Secondly, the driving potential could be activated, causing Bloch oscillations in the synthetic dimension bulk.
Note Added: We note that Bloch oscillations of driven dissipative solitons were recently detected in a photonics synthetic dimension 40 .