Incommensurable matter-wave jets in quasi-1D geometry
Introduction
Incommensurability is encountered in aperiodic crystals with incommensurable phases1,2,3, charge and spin structures4,5,6,7,8, twisted moiré bilayer experiments9 with Hofstadter butterfly energy structure10, and surface layers11. In the context of cold atoms, incommensurable optical lattices have been used to study wavefunction localization phenomena in disordered potentials within the Aubry-André model12,13,14,15,16,17,18,19,20.
Matter-wave jets occur when a Bose-Einstein condensate (BEC) is subjected to a modulation of interaction21. The modulation leads to the formation of density waves inside the BEC followed by the emission of jets22. The microscopic process behind jet emission is photon-stimulated two-atom collisions21,23, where a quantum of the magnetic field modulation energy is transferred to the atom pair. First-order jets form from pairs of atoms that each absorb half of the energy of the photon. Further collisions lead to higher-order jets having velocities that can either be rational multiples of first-order jet velocity (commensurate jets) or not (incommensurate jets).
In our previous work we studied the formation of first- and second-order commensurate jets in a quasi-one-dimensional geometry24. We observed the jets experimentally and in numerical simulations, however, no incommensurate jets were seen, because they only occur at higher modulation amplitudes. Previous studies of matter-wave jets in two dimensions (2D) have revealed the formation of density waves inside the condensate prior to the emission of jets and explained the process of formation of higher-order jets22,23. In ref. 23 higher-order jets with velocities (sqrt{2}) and 2 times the velocity of first-order jets were directly observed, commensurate with the first order. Additionally, a hint of incommensurate (sqrt{3}) jets was reported. It is important to note that in 2D geometry the (sqrt{2}) jets actually result from a density modulation with the velocity vector at a 45° angle with respect to the primary density modulation, so that each component of the velocity vector is, in fact, commensurate with the primary modulation. In contrast, (sqrt{2}) jets in 1D geometry result from a genuine incommensurability in the BEC.
Here, we experimentally and numerically demonstrate the formation of incommensurable “golden” jets and other higher-order matter-wave jets in a quasi-1D BEC under a single-frequency interaction modulation. The incommensurable golden jets are precursors of the supercontinuum formation, which remains to be experimentally confirmed.
Results
Formation of matter-wave jets
We start by preparing an almost pure BEC of approximately 4000 cesium atoms in a crossed dipole trap. We release it into a channel with radial frequency ωr = 2π ⋅ 90 Hz and simultaneously change the interaction between the atoms to a value of adc that is close to zero via the wide s-wave Feshbach resonance24 with zero-crossing at 17.1 G. This results in the formation of a BEC soliton. For 4000 atoms the BEC forms a soliton at interaction adc = − 2.7a025, where a0 is the Bohr radius. To induce the emission of matter-wave jets, we modulate the magnetic field and therefore the interaction between the atoms as (a(t)={a}_{{{{rm{dc}}}}}+{a}_{{{{rm{ac}}}}}sin (omega t)) for a time tp, where aac and ω are the amplitude and the frequency of the modulation, respectively21 (see Methods). We let the system evolve for te in the channel and take an absorption image after a 15 ms time-of-flight expansion. The results are shown in Fig. 1a. In addition to previously observed jets of first order (J1) and second order (J2)24, we observe two pairs of jets with dimensionless velocities
where ({v}_{{{{rm{0}}}}}=sqrt{hslash omega /m}) is the velocity of the first-order jets (J1), ℏ is the reduced Planck constant, and m is the mass of an atom. We name these “golden jets” (G2) for their golden ratio dimensionless velocities. The irrational velocity ratio between first order and golden jets implies the coexistence of two incommensurable density waves inside the condensate prior to jet emission. Figure 1a shows three experimental realizations showing golden jets. As can be seen in Fig. 1a the golden jets are not necessarily symmetric around J0 and change in each experimental run despite the supposed symmetry of the involved formation processes. We attribute this to the asymmetry in the initial BEC which is slightly inhomogeneous and fluctuates between experimental runs.

a Absorption images of three different realizations of golden jets. All measured at the same experimental parameters: modulation frequency 4 kHz, amplitude aac = 90a0 (which equals to interaction parameter amplitude kac = 21.9), modulation time tp = 10 ms, and evolution time te = 40 ms. The images show the original BEC (J0), first order jets (J1), second order jets (J2) and as well as two pairs of “golden” jets (G2), but not (sqrt{2}) (D2) jets. b Schematic of possible jets for increasing number of subsequent collisions (order encoded by color).
The origin of golden jets can be understood from the conservation of energy (frac{m{v}_{1}^{2}}{2}+frac{m{v}_{2}^{2}}{2}+hslash omega =frac{m{v}_{1}^{{prime} 2}}{2}+frac{m{v}_{2}^{{prime} 2}}{2}), and momentum: (m{v}_{1}+m{v}_{2}=m{v}_{1}^{{prime} }+m{v}_{2}^{{prime} }), where vi (({v}_{i}^{{prime} })) are the initial (final) velocities of the atoms and the momentum of the photon is neglected. A pair of atoms in J1 is formed from a pair of colliding atoms from the original condensate (J0) by absorbing one quantum of energy from the modulating magnetic field. Similarly, a pair J0-J2 forms from a pair of J1 atoms with the same velocity and a pair G2-G2 from a pair J0-J1. In general the final velocities of the resulting jets are
where u1,2 are the initial dimensionless velocities of the atom pair. With each absorbed photon the number of possible jets increases. By the number of photons that an atom in a jet has absorbed we define the order of the jets: the initial BEC is zeroth order, J1 are first order, J2, G2 and D2 are second order and so on, see Fig. 1(b) and Table 1. The total number of jets increases as a double exponential with order, going as 1, 3, 11, 109, 11605,…, and can be approximately calculated as (1.794{9}^{{2}^{n}}) (from the fit up to n = 4, where n is the order). Experimentally, we observe J1, J2 and G2, but we have not been able to observe jets with dimensionless velocity (pm sqrt{2}) (D2) or any third order jets. J2 are formed from J1 which move with the same velocity so the process is equivalent to the primary process in a moving frame, while the atoms that form G2 move relative to each other, which reduces the interaction time, making golden jets rare. For D2 the velocity difference is two times higher, explaining why we do not observe them. Additionally, the formation of J2 is more probable due to bosonic stimulation, since they form in pair with J0, which is present in the system beforehand23.
Numerical simulation of Gross-Pitaevskii equation
We turn to numerical simulations of 1D Gross-Pitaevskii equation (GPE) to better understand the conditions for the formation of each order of jets. As in our previous work24, we start with a soliton as the initial state and modulate the interaction parameter (k(t)=Na(t)/{a}_{r}={k}_{{{{rm{dc}}}}}+{k}_{{{{rm{ac}}}}}sin (omega t)) for a time tp = 8π with different frequencies ω and amplitudes kac. We leave the wavefunction to evolve for te and examine the resulting density profile n(x) = ∣ψ(x)∣2. The details about the simulation are provided in Methods and in our previous work24. In the simulations, the time t is given in units 1/ωr, ω in units ωr, and distances in units ({a}_{r}=sqrt{hslash /m{omega }_{r}}).
An example of the final density profile can be seen in Fig. 2a where all visible jets are marked. Figure 2b shows the time evolution of the density profile from the start of modulation. In contrast with the experiment, we observe additional jets of third order J3 and G3 which form from J2-J2 and J1-J2 pairs, respectively. Surprisingly, G2 and G3 have comparable amplitude even though G2 jets are second order whereas G3 only arise in third order; this seems to be due to the large amplitude of J2 jets (from which G3 are formed) in the simulation.

a The final density profile for interaction parameter amplitude kac = 27.1. b The time evolution of the density profile for the same modulation amplitude. c Fraction of atoms in J0, J1 and J2 in linear scale as a function of the modulation amplitude. d Fraction of atoms in jets J1, J2, J3, G2 and G3 in logarithmic scale as a function of the modulation amplitude. On both (c, d) the horizontal dashed line at 0.006 represents the threshold for a presence of each jets used to construct the phase diagram. All simulations in this figure have evolution time t = 40, modulation frequency ω = 33, interaction parameter offset kdc = −0.56, and modulation time tp = 8π.
To quantify the size of jets, we calculate the integrals over an interval around each peak in the density profile in momentum space ∣ψ(q)∣2, which gives the fraction of the atoms in each type of jets. For J0 we integrate from −0.56 to 0.56 and for other jets we integrate an area ± 0.16 around each jet. The total integral is normalized, ∫∣ψ(q)∣2dq = 1, and ψ(q) is the Fourier transform of ψ(x). The fraction in the central peak J0 only starts decreasing after the threshold amplitude for J1 formation is reached, as can be seen in Fig. 2c. The J1 first increases steeply with amplitude and then levels out due to J0 depletion and exitations to the higher-order jets. It can be seen that only the first order has a sharp threshold (Fig. 2c), while the higher orders increase smoothly from zero once the first order exists (Fig. 2d). D2 are not observed even in simulations due to the onset of a wide background at high amplitudes needed to create them (discussed below).
Based on calculated atom fractions, we construct a phase diagram of visibility of different sets of jets, shown in Fig. 3a. A new region is declared to exist when the fraction for the next jet exceeds 0.006. We choose this value so that it is much larger than the noise in the simulation, but smaller than the atom fraction of the jets at high modulation amplitudes. As the amplitude kac increases, additional jets appear, represented by new regions in the phase diagram. As discussed above, the transitions are actually continuous and not discrete, except from the “no jets” to the “J1” region. G3 jets occur in parts of the region marked with “J1+J2+J3+G2″ in Fig. 3a.

a Phase diagram of jets from GPE simulations showing which jets emerge for a range of modulation frequencies and amplitudes. Different colors represent regions where different types of jets are above the threshold. b Final density profiles in momentum space ∣ψ(q)∣2 (t = 8π) for modulation frequency ω = 33 [along the dotted line in (a)]. c, e Experimental phase diagram of jets, measured at modulation times tp = 30 ms and 10 ms. The grey lines indicate where the total number of atoms falls by 25% and 50% from the initial number of atoms in the condensate. The dashed black line in (c) is a fit to the boundary between the “no jets” and “J1″ region, proportional to ω1/2. The error bars correspond to the variance of the atom number and the precision of determining the amplitude corresponding to the boundary between the regions. d, f Averages of 50 absorption images as examples of different points in the phase diagram (marked with black crosses in the diagrams above). From bottom to top the images represent regions with: no jets, only J1, J1 and J2, and J1, J2 and G2 (top two images). The images were taken for evolution time te = 42.7 ms. The interaction parameter amplitude is kac = Naac/ar, where N = 4000 ± 200 is the number of atoms in the BEC, before we start modulation. aac is calculated from the amplitude of the pulse by a calibration determined using the Feshbach resonance (see Methods).
Figure 3b shows the final density profiles in momentum space for a fixed frequency ω = 33. The onset of J1, J2 and other jets can be clearly seen as the modulation amplitude increases. For very high amplitudes the discrete jets start disappearing in the increasingly uniform background of a multitude of higher-order jets. The observed spread from the single momentum mode to a wide momentum distribution resembles the generation of supercontinua in nonlinear optical media, where various instabilities cause an initially narrow frequency spectrum to disperse into a wide continuum26. The boundary for the supercontinuum in Fig. 3a is determined with an entropy-like measure (S=-int| psi (q){| }^{2}log left(| psi (q){| }^{2}right)dq), where we integrate over the whole interval. We scale the threshold with the square root of the modulation frequency, to take into account the scaling of the jet velocities with modulation frequency. When (S > 0.38sqrt{omega }) we say that it has transitioned into the supercontinuum region (the transition is actually gradual, the threshold value was chosen to be above numerical noise and visually without separate jets). The threshold boundary for J1 formation is proportional to ω1/2, as previously reported21,24.
Experimental phase diagram of jet formation
We explore the phase diagram of jet formation experimentally as well. Figure 3c, e show the experimental phase diagrams of jet formation for modulation times tp = 10 ms and 30 ms, respectively. We determined the regions of the phase diagram where we observe J1, J1 and J2 (J1+J2), and J1, J2 and G2 (J1+J2+G2).
In contrast to simulation, in the experiments, the number of atoms in jets varies from shot to shot; for example, see Fig. 1a. This was previously observed in 2D27, where the distribution of population in angular modes resembles the distribution of thermal radiation. More information on the statistics of the number of atoms in each type of jet is available in the Supplementary Note 2 and Supplementary Fig. 3. Because of this variance, gathering statistics is necessary to determine the presence of the phase at each point in the phase diagram. We choose the boundary to be at the point where at least one image out of ten shows the corresponding jets. Examples of averages of 50 absorption images of jets from different parts of the phase diagram are shown in Fig. 3d, f and in the Supplementary Fig. 1. Both Fig. 3d, f show experimental images for modulation frequency 3 kHz, but (f) is for the shorter modulation time tp = 10 ms and (d) is for tp = 30 ms. Note that in Fig. 3d the distance between the jets and the initial condensate varies for different modulation amplitudes. This is because jets are emitted sooner if the modulation amplitude is higher (see Supplementary Fig. 6). This effect is not as pronounced in Fig. 3f due to the shorter modulation time.
The range of frequencies in the experimental and the theoretical phase diagram are similar (the frequencies shown in Fig. 3a range from 0 to 7.2 kHz (from 0 to 80 in units of ωr = 2π ⋅ 90 Hz), whereas the experimental phase diagrams in Fig. 3c, e were measured for frequencies between 1 kHz and 7 kHz). On the other hand, the values of kac required for excitation of jets are approximately two times lower in the experiment than in simulation. The difference is most probably because the simulation describes a true 1D system, while in the experiment, the BEC is a 3D object with jets confined into quasi-1D system by a dipole trap beam.
For tp = 30 ms, we are not able to observe golden jets for modulation frequencies higher than 3 kHz, even though they are observable using the shorter tp = 10 ms. In the simulation, the modulation time is tp = 8π, which corresponds to 44 ms. In the experiment, we used shorter modulation times, because longer modulation times cause more atom losses and increase the incoherent background, formed by the excited atoms that are not part of the jets24, making the occurrence of jets rarer and more difficult to observe (see Supplementary Figs. 1 and 2). However, we see from the experimental phase diagrams that the shape of the threshold does depend on tp. For a short tp there is no strong dependence on the modulation frequency, whereas for the longer tp the threshold boundary for the formation of jets is proportional to ω1/2, as expected from the simulation and previous experiments24. Additionally, the longer modulation times sometimes cause the jets to be emitted twice, making the interpretation of the absorption images less reliable.
As we increase the amplitude of modulation in the experiment, the total number of atoms starts decreasing and the incoherent background between the jets becomes stronger. This can be seen in Fig. 3c and e, where we plot the thresholds at which the total number of atoms is 25% and 50% smaller than the number of atoms in the initial BEC. The decrease in the number of atoms is likely due to partial collapses of the BEC, caused by the strong attractive interactions the BEC is subjected to during the modulation, leading to three-body losses. The incoherent background is formed by the excited atoms that are not part of the jets. Their velocities have random directions and that appears as a homogeneous background between J0 and J1, and to a lesser extent between J1 and J2 (best seen in Supplementary Fig. 2). The momentum distribution of the incoherent background remains of the same width for increasing amplitude and does not increase beyond the position of J2 jets as the supercontinuum would. Atom losses and increased fraction of incoherent background, both not accounted for in the simulation, decrease the effective modulation amplitude and prevent the supercontinuum from forming.
To summarize, the experimental phenomena that the 1D GPE simulation does not account for are the atom losses into the incoherent background and three-body losses. For a complete description of the experiment, one would not only need to take into account the three-dimensionality of the experiment, but also the higher order interaction between atoms (e. g. three-particle interaction) and incoherent processes, that are not described by the mean field picture. Despite all of these differences, the discrepancy between the simulation and the experiment only becomes significant at higher modulation amplitudes.
Discussion
In this work, we show the formation of matter-wave jets with golden ratio dimensionless velocity ({u}_{{{{rm{G2}}}}}=(1+sqrt{5})/2), which form when a BEC in a quasi-1D geometry is subjected to a single-frequency interaction modulation. These jet velocities and the corresponding density waves are not commensurate with the previously observed jets with dimensionless velocities uJ1 = 1 and uJ2 = 2. A hint of incommensurate (sqrt{3}) jets has previously been observed in similar experiments with 2D geometry23, but, to the best of our knowledge, not in 1D. It needs to be noted that the occurrence of incommensurate density waves can easily be achieved with two-frequency driving (like in ref. 28, but with the frequencies that produce first-order jets that are not commensurable, see Supplementary Note 3 and Supplementary Fig. 4), but in the case of golden jets incommensurability emerges from single-frequency driving.
Numerically and experimentally, we explore the formation of incommensurable jets for a wide range of modulation frequencies and amplitudes, resulting in a phase diagram of possible jets, while also uncovering an unexpected supercontinuum region. The supercontinuum likely arises due to the double-exponential growth of the number of possible modes, while the wavenumber only grows linearly, resulting in a densely populated momentum space. The supercontinuum is potentially experimentally observable if the incoherent background could be reduced by increasing the radial confinement and suppressing three-body losses.
The demonstrated formation of incommensurable matter-wave jets could, for example, be used in multispecies experiments. The incommensurable density waves inside one species could be used as a disordered potential for the other species (provided there exist appropriate Feshbach resonances). This would eliminate the need for an external disordered potential or a regular lattice that is usually needed to capture impurities, which act as a source of disorder12,13,14,15,16,17,18,19,20.
Methods
Preparation of the BEC
We start the preparation of the BEC with a collimated beam of atoms slowed down using a Zeeman slower. The slowed atoms are captured by a magneto-optical trap in the main experimental chamber. We usually load ~ 60 million atoms with temperature ~ 70 μK into the magneto-optical trap over 10 s. After loading we compress the atoms and use optical molasses to cool the atoms to ~13 μK and transfer them to the internal ground state F = 3. The atoms are then transferred into the Raman lattice. With degenerate Raman sideband cooling we cool about 20 million atoms in the absolute ground state (leftvert F=3,{m}_{F}=3rightrangle) to 1 μK, which we then transfer into the large crossed dipole trap. After 0.5 s of plain evaporation in the dipole trap we start ramping up the dimple trap. The dimple beams are linearly ramped to their maximum powers in 1.5 s and held for another 400 ms to let the atoms thermalize. Then the large dipole trap is turned off and evaporation in the dimple trap begins. The powers of dimple beams are exponentially ramped down in 6 s and after 2 s we start with additional exponential rampdown of the magnetic field gradient. The combination of the two techniques leads to efficient evaporation and at the end of it we produce a BEC of approximately 4000 atoms. More details can be found in refs. 25,29.
Numerical simulation
We use the split-step spectral method to solve the Gross-Pitaevskii equation
in the dimensionless form30, where time t is given in units 1/ωr, ω in units ωr, and distances in units ({a}_{r}=sqrt{hslash /m{omega }_{r}}). k is the interaction parameter k(t) = Na(t)/ar. We solve the equation on an uniform 1D grid of 65536 points in the interval [−1200:1200]. The initial state of the BEC is determined by an imaginary-time evolution, while the calculation of the real-time dynamics is performed with a time-step of τ = 0.001. The code is written in Wolfram Language (Mathematica) and it can solve for the ground state through imaginary-time propagation of an initial approximation, and determine the time-dependence in real-time. Most of the computation time is spent in performing the (inverse) Fourier transform, which is performed highly efficiently in Mathematica, thus the performance of the code is comparable to that of compiled codes while it is significantly more flexible. Further details and examples can be found in ref. 31.
Calibration of modulation amplitude
We measured the amplitude of the modulation of interaction by measuring the splitting of a Feshbach resonance at 19.9 G32. In a typical measurement of Feshbach resonance, we measure the number of atoms remaining in the trap after switching the magnetic field to a chosen value for 100 ms. Around the magnetic field corresponding to the Feshbach resonance we observe enhanced atom losses. We repeat the same measurement with the addition of an oscillating magnetic field at the various frequencies used for jet formation. In this case, we observe that the dip in the number of atoms is widened (see Supplementary Note 4 and Supplementary Fig. 5). The width of the dip corresponds to the peak-to-peak amplitude of the magnetic field oscillations. For jet formation we modulate the magnetic field around 17.1 G, where the wide s-wave Feshbach resonance of cesium has the zero crossing32. The interaction as a function of the magnetic field around the zero crossing is known32, and we use it to determine the amplitude of the interaction modulation. By measuring the effect on the atoms, we directly measure the magnetic field at the position of the atoms. The measurements are done for the same amplitudes and frequencies as for jet formation, meaning that all the non-nonlinearities and other effect of the electrical circuit producing the oscillating magnetic field are taken into account, including the amplifier (Lab Gruppen IPD 1200).
Responses