Anomalous periodicity in superpositions of localized periodic patterns

Interference between overlapping periodic patterns gives rise to important phenomena, such as Moiré fringes, appearing when the patterns have different periods or orientations. Here we present a novel phenomenon, applicable to both the classical and quantum regimes, where two one-dimensional localized periodic patterns with the same period interfere to create fringes with anomalous periodicity. We analyze the effect theoretically and demonstrate it with atomic matter waves. When a central parameter of the system is scanned continuously, we observe a discontinuous but piecewise-rigid periodicity of the resulting fringes. We show that this is a universal phenomenon that emerges from a superposition of two spatially shifted localized periodic patterns of any source or nature when they interfere with a global phase difference. The rigidity of the spectrum becomes even more robust for a coherent superposition of non-overlapping wavepackets, although the conventional interferometric visibility drops to zero. The effect is expected to appear in space and time, as well as in the momentum distribution of quantum particles.

Interference between overlapping periodic structures or waves gives rise to a variety of classical and quantum phenomena. One example are Moiré patterns, which are a result of alternating constructive and destructive interference between periodic patterns with a slightly shifted periodicity, or relative rotation. Such patterns appear in many interesting contexts (e.g., [1]). While interference in general is the basis of classical and quantum wave mechanics, it still remains a rich ground for new experimental findings and theoretical interpretations. This has been made evident, for example, with the recent emerging field of overlapping periodic thin layers and their effect on the electronic and magnetic properties of van der Waals hetero-structures or correlated oxides [2][3][4][5]. Novel points of view regarding interference are also at the heart of studies on the double-slit experiment [6], Bohmian mechanics [7], and weak values [8]. Most relevant to this work, are intricate effects due to the overlap of localized wavepackets, for example in the context of anomalous arrival times, where interference is an underlying mechanism [9][10][11][12][13].
Here we report on anomalous features arising from interference between localized periodic phenomena, which have applications to a wide range of classical and quantum systems. In contrast to Moiré patterns, the anomalous features of our system arise from interference between two patterns with the same period thus having a constant phase difference throughout the patterns. Nonetheless, the phase of the combined pattern shows a phase gradient due to the relative amplitudes of the two patterns which vary along the patterns. This phase gradient gives rise to a new period of the combined pattern which exhibits surprising features.
Our experimental model-system utilizes two parallel interferometers, each creating its own periodic fringe pattern. The interferometers use matter waves of ultracold atoms precisely controlled by an atom chip [14], as described in figure 1.
The general system we have in mind is a superposition of two periodic one-dimensional patterns (hereafter named constituents) displaced from each other by Δz and both having the same wavenumber κ Figure 1. Experimental sequence (schematic representation, not to scale): a BEC is prepared in state |2 ≡ |F = 2, m F = 2 , represented by the red line, and released from the magnetic trap (not shown). After a time-of-flight (TOF) of ∼1 ms an RF π/2 pulse (light blue) transfers each atom into an equal superposition of |1 + |2 , where |1 ≡ |F = 2, m F = 1 is represented by the blue line. The magnetic field gradient pulses are in pink. See text for more details. After a TOF of ∼10 ms, the expanded wavepackets overlap and interfere with one another, forming two interferometric fringe patterns, one of the |1 state and one of the |2 . Our absorption imaging method is invariant to the internal state and so a sum pattern is formed on the charge coupled device (CCD) image. Red and blue represent the two constituent periodic patterns, while purple represents the sum pattern. and finite-size σ of their envelopes, which we describe by Gaussians. The superposition is described by f (z) = e −z 2 /2σ 2 e iκz+iθ 1 + e −(z−Δz) 2 /2σ 2 e iκ(z−Δz)+iθ 2 . (1) The complex function f(z) may represent a wave function of a quantum particle in a superposition state. Alternatively, the real part of f(z) may represent a classical field or density variation, where the phases θ 1 , θ 2 determine the positions of the maxima of the periodic functions (cosines) with respect to their corresponding envelope centers. We are especially interested in the case θ 1 = θ 2 , which appears, for example, in phase-coherent sources or as a result of coherent splitting and non-dispersive propagation, and is also a fundamental feature in some coherent sources, such as our cold-atom model-system. The phase difference between the two constituents Δφ = κΔz, is a global phase independent of position z and hence the interference between them is similarly constructive or destructive everywhere. The localized character of the system is best defined by the number of periods, a dimensionless constant proportional to κσ (typically 3-5 in our model-system).
Our main observable is the periodicity wavenumber K S of the sum pattern in equation (1). In quantum sources K S represents the peak of the momentum distribution (most probable momentum) of the quantum particles, and hence the following features reported for the spatial periodicity of the combined interference pattern apply similarly to the momentum of quantum particles in a superposition state as in equation (1). By varying a system parameter that controls both Δφ and κ, we observe the surprising effect of rigid periodicity, i.e. no variation of K S within certain ranges of the scan in which Δφ does not cross an odd multiple of π, as demonstrated in figure 2. The observed rigidity is a result of an interplay between the relative phase of the constituents and their period, as discussed below. This interplay is shown to be related to a fundamental conservation law in the quantum system that generates the signals in our system. Between these rigidity ranges, sharp transitions (discontinuities, jumps) appear in K S . Below we explain the origin of the effects and show that at in the limit of high Δφ, where the two envelopes are completely separated, a robust and uniform rigidity emerges in Fourier space, whereby all the jumps are of the same height.
The experimental procedure is depicted in figure 1. Our experiment begins by releasing a Bose-Einstein condensate (BEC) of about 10 4 87 Rb atoms from a magnetic trap below an atom chip. We initially prepare the BEC in the state |F, m F = |2, 2 , and then create a superposition of the two spin states |2, 2 ≡ |2 and |2, 1 ≡ |1 by applying a π/2 radio-frequency (RF) pulse. These two states constitute an effective two-level system, as all other states in the F = 2 manifold are pushed out of resonance by the non-linear Zeeman shift generated using an external bias field (see appendix A). A Stern-Gerlach interferometer (SGI) is then Figure 2. Anomalous pattern periodicity-rigidity and jumps: (a) the measured wavenumber K S (purple, as in figure 1) vs the deceleration pulse duration, T 2 (figure 1). The absolute value of the FT of the sum patterns (CCD images shown in the insets) is calculated from the data, and the value of the maxima, K S , is presented. For values of T 2 where a secondary peak is detected with the relative intensity of at least 20%, two points are plotted with the dot size representing the relative intensity of each peak. The error bars are calculated as the standard error of the mean (SEM) over several iterations. As can be seen, although σκ is constant, more fringes appear in the CCD insets as T 2 becomes smaller. This is due to the growing Δφ which is tantamount to a decreasing overlap between the two constituent patterns. (b) Visibility of the interference pattern, V, vs T 2 . The minima in the visibility, emphasized by the vertical dashed gray lines, correspond to the periodicity jump locations in (a). In both (a) and (b) the green line represents the results of a complete numerical analysis of wavepacket propagation based on the exact experimental conditions. For low values of T 2 , the simulation overshoots the observed visibility, as under these conditions 2π/κ is smaller, and as the clouds are moving (free-fall), the limited optical resolution gives rise to smearing and consequently a smaller visibility.
implemented by using a series of two magnetic gradient pulses (gradients along the axis of gravity, z), which are generated by running currents on the atom chip (more details on the setup can be found in reference [15]). The first gradient pulse, of fixed duration T 1 = 4 μs, splits the superposition into two momentum components, which then freely propagate during a delay time T d . During this delay time the spin state is manipulated by a π/2 RF pulse, after which the two wavepackets have equal amplitudes to be in spin state |1 . Then they are decelerated relative to each other by a gradient pulse of duration T 2 . While the first gradient pulse separates the two wavepackets by a spin-dependent force, the second gradient pulse decelerates the relative motion by applying an inhomogeneous force along z on the two wavepackets that have the same spin state but are centered at different positions. The same gradient pulse applies a twice stronger force on the two spin-|2 wavepackets and drives them away from the experimental region, so they are ignored in this experiment. After the deceleration pulse, we apply a third π/2 RF pulse that duplicates the pair of wavepackets of spin-state |1 into two wavepacket superpositions of |1 and |2 . We then apply a third gradient pulse of a fixed duration T 3 = 30 μs. The spin-dependent force of the third gradient pulse gives different momentum kicks to the two spin states such that after free-space expansion and overlap, two interferometric fringe patterns translated from each other in space are formed. As our imaging is insensitive to the spin state, we obtain on the CCD the sum image of the two interference patterns, whose positive spatial-frequency part has the form of equation (1) with θ 1 = θ 2 . Note that the inhomogeneity of the third gradient pulse shifts the periodicities of the interferometric fringe patterns of the two spin states such that they have slightly different wavenumbers κ 1 = κ 2 . However, the wavelength 2π/|κ 1 − κ 2 | corresponding to this difference is much larger than the size σ of the interferometric fringe patterns, such that Moiré effects due to spatially varying phase difference between the two patterns are negligible and we may assume κ 1 ≈ κ 2 = κ (see appendix B). Figure 2 presents the measured wavenumber K S of the sum pattern, which is extracted from the Fourier transform (FT) of the data, as a function of T 2 . The value of K S exhibits a clear rigidity, and singularity points at which this value abruptly changes. The rather good agreement with the data of our numerical simulation, in which care was taken to take into account all the experimental conditions (see appendix C), indicates a good understanding of the experimental apparatus. In what follows we gain insight into the . Emergence of anomalous periodicity: (a) two localized one-dimensional periodic patterns (blue and red), with wavenumber κ and phase difference Δφ, join to create a sum pattern (purple) with wavenumber K S . (b) The phase of the sum pattern (color code) projected on the outline of the two constituent periods. The phase gradient is not constant and has a minimum value at the center, giving rise to a narrow FT peak at K S , which is different from κ (see text). At one edge, the phase of the blue pattern dominates, as the amplitude of the blue pattern dominates, see (a), and at the other edge, the phase of the red pattern dominates. (c) The sum pattern wavenumber K S (z-axis) vs κ (x-axis) and Δφ (y-axis), as calculated from equation (2). The two parameters κ and Δφ may interplay to form the rigidity plateaus of K S around points where Δφ is an integer multiple of 2π. This interplay is represented by the red curve, which shows the actual experimental trajectory giving rise to the data in figures 2 and 4. This line is always parallel to the equi-K S lines except for the discontinuity (jump) points, manifestation of the rigid periodicity. In the model presented here κσ is a constant. Specifically, we take the number of periods in the system (edge-to-edge size 4σ) to have the independently measured value of N p = (2/π)κσ = 5.6. universal origin of the observed effects and in particular how they emerge in our system from the dependence of the model parameters κ and Δφ on the central experimental parameter T 2 .
The principles that stand at the basis of the effect and the general behavior of the model are depicted in figure 3. As demonstrated in figure 3(a), the periodicity K S of the sum pattern is shifted from κ due to a phase gradient which develops along the system. Because of the spatially varying relative magnitude of the two underlying constituent patterns, at one edge of the sum pattern the phase of the first periodic pattern dominates, while at the other edge, the phase of the second periodic pattern dominates. This gives rise to a phase gradient which changes K S as a function of Δφ, for any κ, as demonstrated in figure 3(b) (see a quantitative analysis in appendix D). The dependence of K S on T 2 in figure 2, is understood as an interplay between the fundamental parameters κ and Δφ of equation (1). The specific conditions formed in our experiment and explaining figure 2, are depicted by the red trajectory in the κ-Δφ plane in figure 3(c). In order to understand this interplay we need to elaborate on the crucial role of T 2 in our experiment. On the one hand it determines the periodicity κ of the interferometric fringes by controlling the spatial and momentum difference between the same-spin wavepackets before they expand to form interferometric fringe patterns. On the other hand, as it determines the absolute distance of the atoms from the chip during the translation pulse, it influences the magnetic field gradient (as the gradient is not homogeneous) and hence the differential momentum applied to the two spins, and this in turn, determines the final spatial translation Δz and therefore the relative phase Δφ.
Let us also note that in our experiment κσ is a constant due to a general conservation law concerning the unitary evolution of a pair of Gaussian wavepackets of the same spin in free space or in smooth potentials (see appendix E). This law applies to our interferometric sequence following the splitting pulse whose duration T 1 is kept constant in the experiment. The quantity (κσ) 2 + (Δz/2σ) 2 , where Δz is the distance between the Gaussian centers, is a constant of the evolution, which is approximately equal to κσ at the time of observation. This conservation is most vividly visualized by the evolution of the Wigner function in phase space [16]. The unitary evolution during the interferometric sequence is nothing but a phase space rotation with appropriately scaled phase space coordinates [15], hence the shape of the Wigner function, including the number of periods N p = 2κσ/π, is constant (see appendix F).
Quantitative insight is obtained by comparing the experimental data to the model of equation (1) with the experimentally measured parameters, as presented in figure 4. Good agreement with the data is obtained. The model parameters κ(T 2 ), Δφ(T 2 ) and N p = (2/π)κσ = 5.6, are obtained by independent measurements (appendix B). The Δφ corresponding to each jump are presented in the upper x-axis. The , Δφ(T 2 ) and N p = (2/π)κσ = 5.6, are independently measured (see appendix B). The maxima of both the model and the data were normalized to 1. In the inset, the jump height as a function of the number of periods in our finite system. Each line represents a different jump, classified by their Δφ. The jump corresponding to Δφ = π is expected to appear at T 2 ≈ 300 μs, and is however below our noise limit for detection. inset extends the model to an arbitrary number of periods N p and shows that the magnitude of the jumps in K S for each specific Δφ is determined by N p . For a large N p the jumps disappear, at which point a smooth function of K S ≈ κ as a function of T 2 is expected.
In figures 2-4 we extract K S numerically by finding the maximum of the absolute value of the positive spatial frequency part (K > 0) of the FT (AFT) of the sum pattern. By Fourier transforming the model function f(z) of equation (1) we obtain Δφ .
The peak of the AFT is at K S = κ when Δφ = 2πn is an integer multiple of 2π, where the cosine function is peaked at the same K = κ as the Gaussian peak. In general, when Δφ = 2πn + α, where −π < α < π, the cosine peak shifts to K = κ/(1 + α/2πn), and therefore the peak of the AFT shifts to K S < κ or K S > κ, depending on whether α is positive or negative, respectively (see appendix G for more details). Note that for most of the data points of our experiment a similar wavenumber K S is obtained by fitting the real-space sum pattern (of which equation (1) is only the positive-K part) of the CCD images to the form , wherez,σ andφ are the center position, effective width and phase of the sum pattern, respectively, v is the visibility and A is a normalizable constant. This fitting is valid as long as the two patterns forming the sum pattern significantly overlap. Finally, the effects described above take quite a different and surprising form when applied to quantum particles. Equation (1) may be used to describe the wave function of a quantum particle in a superposition of two wavepackets. As demonstrated in figure 5, such a wave function may be readily obtained at the output of a Mach-Zehnder interferometer. In this case the square of the AFT in equation (2) represents the momentum distribution (or probability) of the particle. The result is especially interesting when we consider a system where the source contains particles with a large range of input momenta but their Figure 5. Discrete peak momentum: (a) a quantum particle with peak momentum κ, represented by a wavepacket with a Gaussian envelope, enters a Mach-Zehnder interferometer and splits at the beam splitter BS1 into two paths with a length difference Δz (for example, due to a different mirror configuration M1 vs M2 and M3). After recombination at BS2 the wave function at each of the output ports is a superposition of two wavepackets whose centers are separated by Δz with a phase difference Δφ = κΔz as in equation (1) with f(z) → ψ(z). This requires θ 1 = θ 2 , which is fulfilled, for example, in propagation in a non-dispersive medium (e.g., photons in vacuum). (b) Origin of peak momentum rigidity: the momentum distribution P(K) ∝ | dz e −iKz ψ(z)| 2 (solid curves for two values of κ) is the square of the AFT in equation (2): a product of a Gaussian envelope of width σ (dashed curves) and a cosine function cos 2 (KΔz/2) with Δz = 4σ. In both the top panel (κΔz = 9.8π) and the bottom panel (κΔz = 10.8π) the dominant peak of the momentum distribution is at K S ≈ 10π /Δz (dotted line)-almost independent of the input peak momentum κ. (c) Momentum distribution (heat map) and most probable momentum (solid curve), as a function of input momentum κ where κσ = 8 is kept constant while κ is scanned. A discrete spectrum of the peak momentum K S emerges when Δz σ (no overlap).
splitting and propagation involves a constant delay of one wavepacket in the superposition with respect to the other, as it occurs in non-dispersive propagation. In this case the momentum distribution at the output of the interferometer is expected to have discrete peak values K S = 2πn /Δz, as shown in figure 5(c) if the delay generated by the interferometer is longer than the wavepacket width such that Δz σ. This effect is the universal and robust limit of the rigidity of periodicity observed in our experiment.
Our model may be used to describe the interplay between any two modulated pulses, with a similar modulation frequency and number of periods. These could be, for example, sound or electromagnetic waves, as well as more exotic phenomena. The anomalous features described above may be observed for any phase-coherent source emitting pairs or trains of pulses, as long as there exist detectors with a bandwidth wide enough to follow the oscillations within the pulses, or detectors that can measure directly the Fourier spectrum of the pulses, or the momentum of the particles.

Acknowledgments
We are grateful to Zina Binstock for the electronics, and the BGU nano-fabrication facility for providing the high-quality chip. This work is funded in part by the Israeli Science Foundation Grant Nos. 856/18, 1314/19, 3515/20, 3470/21, and in part by the program for postdoctoral researchers of the Israeli Council for Higher Education.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).

Appendix A. Detailed experimental scheme
The experiment is based on 87 Rb atoms cooled using a standard reflection-MOT apparatus. The atoms are then loaded into an Ioffe-Pritchard magnetic trap, and forced evaporation cooling, using a RF knife, is applied until degeneracy is reached. The BEC, containing ∼10 4 atoms, is released from the trap in the |2 ≡ |F = 2, m F = 2 state, and the experiment is conducted during free-fall. The initial size and expansion of the BEC after trap release are governed by atom-atom interactions, but as the cloud is relatively dilute these interactions have a minor effect on the interferometric process and may be completely neglected when the atoms expand to form large-scale interference patterns. The whole experiment is performed under a constant bias magnetic field of ∼35 G which isolates an effective two-level system of |2 and |1 ≡ |F = 2, m F = 1 . The wavepacket is transferred into an equal superposition of the two internal states, (|1 + |2 )/ √ 2, using an on-resonance RF π/2 pulse. The transition frequency of this two-level syste is ∼25 MHz. By applying a magnetic field gradient pulse of duration T 1 = 4 μs, the wavepacket is split into two distinct trajectories with different momenta. We note that the values reported here for the duration of the pulses are the values given to the experimental control of the apparatus, while the actual values are systematically shorter by a few hundred nanoseconds. The magnetic field gradient originates from a current of I = 1.1 A flowing along three parallel wires in alternating directions. This configuration of the wires helps reduce the phase noise originating from the chip currents [15].
After the initial splitting pulse, the two wavepackets freely propagate for a time T d1 = 230 μs. During this time we apply a second π/2 pulse which transfers each of the trajectories, previously in a pure state of either |1 or |2 into an equal superposition of (|1 − |2 )/ √ 2 and (|1 + |2 )/ √ 2, respectively. The wavepackets are then decelerated relative to one another by a second magnetic gradient field pulse of varying duration T 2 , which causes the |1 components from each of the wavepackets to have roughly the same momentum. The |2 components are ejected from the interferometer region-of-interest and are ignored for the rest of this experiment. This decelerating is possible due to the non-linearity of the magnetic field; that is, each wavepacket in internal state |1 feels different acceleration since the magnetic field gradient is not constant in space. The same non-linearity also causes the wavepackets to go through a focal point and then expand at an increased rate.
After the deceleration pulse, we apply a third and final RF π/2 pulse, creating two superpositions of |1 and |2 . In total, our quantum system now consists of two spatially separated wavepackets, each in a superposition of |1 and |2 . The atoms freely propagate until we apply a third magnetic field gradient pulse of duration T 3 = 30 μs. The third pulse has the opposite polarity compared to the first two pulses, such that the acceleration due to the magnetic field gradient is directed upward. By reversing the polarity of the third magnetic field, we can achieve slightly longer measurements since the BEC remains in the field of view of the CCD for a longer time. This third pulse imparts a differential phase Δφ and a differential position Δz, on the two spin interference patterns. The timing of the third magnetic field gradient pulse is T d2 = 410 μs after the start of T 2 .
The bias magnetic field is shut down 660 μs after the last magnetic field gradient pulse, after which the wavepackets fall under gravity and expand for another 14 ms. At the end of the experimental cycle, we image the interference pattern using a standard absorption imaging technique. Since the bias magnetic field is turned off before the imaging pulse, our imaging beam is invariant to the two energy levels and can not distinguish between |1 and |2 . This measurement is repeated for different values of T 2 , and its results are presented in the main text.
As auxiliary measurements, we repeat the same experimental sequence described above, where we change the third RF pulse, of duration T R3 and Rabi frequency Ω, from ΩT R3 = π/2 to ΩT R3 = 0 or ΩT R3 = π. Effectively, these two extra measurements, presented in figure 6, measure the periodicity of the single-state interference pattern, κ i .  figure 1 is repeated twice, once without the last RF pulse, ΩT R3 = 0 (blue) and a second time with the last RF pulse twice as long ΩT R 3 = π (red). From these measurements we get the single-state interference pattern. Error bars are calculated from the SEM. The solid lines are phenomenological fits.

Appendix B. Data analysis
An example of the raw data can be seen in the insets of figure 2 and in figure 7. The analysis starts by fitting a Gaussian function of width σ 0 to the interference pattern. While the accuracy of this fit is limited, it gives a rough estimate of the envelope function. Using these results, we subtract the Gaussian envelope from the data and calculate the absolute value of the FT along the elongated axis of the interference pattern (z, gravity). We cut the spatial frequencies lower than 1/0.9σ 0 . This removes the central peak of the FT. While there is some redundancy in filtering the central peak (low-frequency cutoff and subtraction of the envelope from the data), we found this method to give more consistent results. We now locate the two highest maxima of the FT spectrum and calculate their relative intensity. We take the highest value as the main peak, K S , and the second peak is recorded only if its relative intensity is at least 20%. These values are plotted in figure 2(a).
To calculate the visibility of the sum pattern we fit the original data in real space to a sine function multiplied by a Gaussian envelope where A,z,σ, v,φ, and c are the fitting parameters. In this fit, we set the value of K S to be the main peak of the FT. Two images from the CCD and the corresponding fits are presented in figure 7. From these results, we take v to be the visibility of the interference pattern, plotted in figure 2(b). For our model we fit the visibility, figure 2 where c, v 0 , a, and φ 0 are the fitting parameters. From the results of this fit we obtain the differential phase, given by where a = (163 ± 2) × 10 3 μs 2 and φ 0 = 1.3 ± 0.2 rad are the result of the fit to v(T 2 ). We also fit each of the single-state wavenumbers, plotted in figure 6, to the function where a i and b i are the fitting parameters. The results of these fits are a 1 = 0.175 ± 0.004 μm −1 μs 1/2 and b 1 = −56 ± 1 μs for the case of ΩT R3 = 0 and a 2 = 0.183 ± 0.004 μm −1 μs 1/2 and b 1 = −56 ± 1 μs for the case of ΩT R3 = π. For our model, presented in the main text, we use κ(T 2 ) = 1 2 [κ 1 (T 2 ) + κ 2 (T 2 )]. For the independent measurement of N p , presented in figure 8, we use the aforementioned values for κ 1 (T 2 ) and κ 2 (T 2 ) together with the values ofσ(T 2 ) extracted from the fit to equation (B1). We then calculate N p (T 2 ) = (2/π)κ(T 2 )σ(T 2 ), and average the results to get N p = 5.61 ± 0.03.

Appendix C. Numerical simulation
The numerical simulation presented in figure 2, was performed by using the wavepacket evolution method [15,18], for a BEC under the influence of time-dependent potentials. After the final gradient pulse of duration T 3 , there are four wavepackets such that each pair corresponding to the same spin state is summed coherently to yield an interference pattern, and the probabilities of the two interference patterns are summed incoherently to yield the sum pattern. The initial state is a BEC of 10 4 87 Rb atoms in a trap of frequencies ω x = 2π × 38 Hz and ω y = ω z = 2π × 113 Hz. Taking into account the number of atoms and the s-wave scattering length of 87 Rb (taken to be 100× Bohr radius) we obtain an atomic density with sizes (in terms of standard deviation) σ x = 3.45 μm and σ y = σ z = 1.23 μm. Note that these calculated numbers are for a pure BEC model with the given number of atoms. The wavepacket sizes in the trap are not measurable directly. The final clouds' measured sizes after the formation of the interference fringes turn out to be larger than the prediction of the numerical model but almost 20%. This discrepancy may be attributed to effects that are not accounted for in the simulation, such as an incomplete condensation in the trap or unknown magnetic field gradients during the trap release.
The simulation used the experimental parameters. To achieve fine tuning of these parameters within their uncertainty range given by their direct measurement we use values that are consistent with the periodicity measurements, particularly the periodicity of the interference patterns of each of the spin states (κ 1 and κ 2 , given in figure 6), which are measured independently of the combined patterns that are the main result of this work. Fine tuning of the parameters was required as the periodicity of the fringe patterns is sensitive to variations of some of the parameters within the range of uncertainty given by their direct measurement. We take the initial trapping distance from the chip to be z 0 = 89.5 μm, the duration of the splitting gradient pulse to be T 1 = 3.75 μs and the chip current during the gradient pulses to be I = 1.122 A. These values are well within the experimental uncertainty, and they reproduce quite accurately the values of κ that are measured independently. In addition, to reproduce the positions of the spatial frequency jumps in figure 2, we take the bias gradient during the whole evolution (except the last 13 ms when the bias is turned off) to be 90 G m −1 . This value is not completely supported by experimental evidence, but it turns out to be the right value to consider for a few possible sources of the magnetic field's inhomogeneity. Figure 9 shows some properties of the wavepackets during their evolution, which are not directly measured in the experiment.

Appendix D. Phase gradient analysis
In this appendix we provide a more intuitive insight into the effects presented in this work by looking at the variation of the phase ϕ of the sum pattern in real space. Consider two infinite sinusoidal periodic patterns with the same wavenumber and amplitude but different phases φ 1 and φ 2 . Their sum is a similar periodic pattern with the same periodicity: cos(kx + φ 1 ) + cos(kx + φ 2 ) = A cos[ϕ(z)], where A = 2 cos 1 2 (φ 1 − φ 2 ) and ϕ(z) = kz + 1 2 (φ 1 + φ 2 ). The periodicity can be defined as the phase gradient ∂ϕ/∂z = k. If the amplitudes of the two patterns are not equal, then the sum pattern phase is closer to the phase of the pattern with the larger amplitude. Furthermore, if these amplitudes vary with z, then the sum pattern's local phase is not linear in z, and one can define a local periodicity ∂ϕ/∂z = k + δk(z). Figure 10. The phase gradient of a superposition of two displaced fringe patterns as the origin of the periodicity jumps. The phase of the left pattern ϕ 1 (z) = κz is linear with a slope ∂ϕ 1 (z) ∂z = κ and represented by a blue dotted curve. The phase of the right pattern ϕ 2 (z) = κz + Δφ (defined up to a integer multiple of 2π) is also linear with the same slope and represented by two equivalent red dashed curves that differ by 2π. The phase ϕ(z) of the sum pattern (solid black curve) is dominated by ϕ 1 (z) at z < −σ and by ϕ 2 (z) at z > σ, while in the central region it gradually shifts from ϕ 1 (z) to the nearest branch of ϕ 2 (z). (a) For Δφ = π + π/9 the phase ϕ(z) decreases along z from ϕ 1 to the nearest (lower) branch of ϕ 2 and hence the phase gradient of the sum pattern is smaller than the periodicity of both constituent patterns, and we get ∂ϕ(z) ∂z < κ (this holds when π < Δφ + 2πn < 2π). (b) For Δφ = π − π/9 the phase increases from ϕ 1 to the nearest (upper) branch of ϕ 2 and hence ∂ϕ(z) ∂z > κ (this holds when 0 < Δφ + 2πn < π). The deviation of the phase gradient ∂ϕ(z) ∂z from κ is maximal when Δφ is close to an odd multiple of π. When continuously scanning Δφ through such a value the phase gradient and hence the periodicity of the sum pattern jumps discontinuously. For these plots we used κσ = 8.
The sum pattern in equation (1) of the main text can be written as a complex function with a pre-factor e iκz and a sum of two translated Gaussians with two z-independent phase factors The local phase of this function in the complex plane is We define Δφ =φ 1 −φ 2 and obtain ϕ(z) = κz +φ 1 − atan sin Δφ e −(z−z)Δφ/κσ 2 + cos Δφ , ( D 3 ) where Δφ = κΔz if θ 1 = θ 2 = φ in equation (1). The local gradient of the phase is then The phase variation is demonstrated in figure 10. In the transition region between the dominance ranges of the two constituent patterns the phase gradient, which represents an additional effective contribution to the wave number, reduces when Δφ is larger than the closest value of 2πn and increases when Δφ is smaller than the closest 2πn. This gives rise to a jump of K S at Δφ = π(2n + 1). For a qualitative explanation of the effect see also figure 11.

Appendix E. Model and conservation laws
In this appendix we derive the basic properties of our model following the experimental scenario. We show why the interference patterns appearing in the experiment should have a form given by the real part of equation (1) with θ 1 = θ 2 , explain the phase relation Δφ = κΔz and prove the conservation law for κσ. This derivation is based on a general argument that allows for interactions in the initial state and also allows for an arbitrary form of the magnetic field B y (z). Consider an initial time t 0 after the second π/2 pulse and before the deceleration pulse of duration T 2 (see figure 1). We consider only states |1 ignoring states |2 (whose separation from the states |1 is Figure 11. Real-space explanation for the shift and jumps of the sum pattern periodicity K S . The local phase in a sum of two Gaussian envelopes centered at z 1 and z 2 (z 1 < z 2 ) with global phases φ 1 and φ 2 , respectively. The phase of the sum pattern is closer to φ 1 when z < 1 2 (z 1 + z 2 ) and closer to φ 2 when z > 1 2 (z 1 + z 2 ). It follows that if Δφ = φ 1 − φ 2 is in the range 0 < Δφ < π (upper line of circles) the phase decreases along z (from about φ 1 to about φ 2 < φ 1 ), corresponding to a negative contribution δk = ∂φ/∂z < 0 to the wavenumber K S . If π < Δφ < 2π, or equivalently −π < Δφ < 0 the phase increases along z (lower line of circles) so that the contribution to K S is positive δk > 0. When the phase difference Δφ is scanned through π the periodicity wave-vector jumps from below k f to above k f . achieved later with the deceleration pulse), thus we have two states |1 with different momenta and we follow their unitary evolution in time. We define the overlap of these states ψ a (z, t 0 ), ψ b (z, t 0 ) as (using with Γ, χ real and the initial states normalized as z |ψ a (z, t)| 2 = z |ψ b (z, t)| 2 = 1. In our experiment the initial time t 0 is after the splitting pulse and quite shortly after the BEC is released from the trap, so that the wave functions include effects of interactions and may deviate considerably from a Gaussian form. We next consider the full evolution with arbitrary magnetic fields during T 2 and T 3 , however we follow two evolution scenarios: one for states that stayed in the spin state |1 during the third π/2 pulse and one for the states that flipped into the spin state |2 during this RF pulse. The superposition of the two scenarios corresponds to the wave function keeping the normalization of the initial states (ψ a (z, t) + ψ b (z, t))|1 . The pair ψ 1a , ψ 1b forms the interference pattern of spin 1 while the pair ψ 2a , ψ 2b forms that of spin 2, the summation of both interferences gives the observed sum pattern. For either scenario the time evolution represented by the evolution operator U i (t, t 0 , z), is identical for ψ ia (z), ψ ib (z) (common Hamiltonian for a given spin state i = 1, 2), hence the overlap at the final time t is From this conservation law of the overlap integral we derive more explicit conservation laws that involve parameters of the wavepackets. We next assume that the wave-packets after the gradient pulses are Gaussians with equal time-dependent width σ i (t) for the pair of wavepackets with the same spin, but different center positions and momenta. In our experiment the final wavepackets are fairly close to a Gaussian and the two spin states have similar widths, see figures 6 and 8. This is indeed expected as long as the two wavepackets, which originate from the same initial wave-packet before splitting, propagate in free space and in a potential that has a constant curvature ∂ 2 V/∂z 2 over the range occupied by the two wavepackets. This condition is satisfied in the experiment since the potential due to the magnetic field during the gradient pulses is generated by current wires whose distance from the atoms (about 100 μm) is much larger than the distance between the wavepackets or their width (a few μm).
If this condition is satisfied then the interference terms have the form where A i = (2πσ 2 i ) −1/4 is the normalization constant and z ia , z ib are the corresponding Gaussian centers. The explicit forms of κ i and φ i are given below, though for now they are not needed. We is the average position of each pair of states. The interference pattern from (E4) is then where By comparing this form to equation (E1) the conservation law for the overlap integral yields our central result Note that all the parameters δz i , σ i , κ i , z i , φ i are time dependent so that the right-hand sides of equations (E7) and (E8) are conserved in time for either spin state and are in fact spin independent, since Γ, χ are common to both spins. Equation (E5) can now be identified with the real part of each of the two patterns in equation (1), with φ i = κ i z i − 2θ i . By comparing this phase term with equation (E8) we identify θ i = −χ/2. As χ is a conserved quantity that is common to both spins we therefore conclude that θ 1 = θ 2 . This is an exact result exhibiting the quantum signature of our implementation of equation (1) in the main text as it follows the quantum evolution from a common source.
Note that σ i and δz i are nearly the same for i = 1, 2, determined mainly by the size, separation and relative momenta of the common parent state entering the third π/2 pulse. While the curvature of the third gradient pulse of duration T 3 affects the sizes and relative distances and momenta within pairs of states differently for each spin, this effect is quite minor due to the fact that T 3 is much shorter than T 2 . As δz 1 ≈ δz 2 and σ 1 ≈ σ 2 , equation (E7) implies that κ 1 ≈ κ 2 , as observed in the experiment (see appendix B).
In general, the phase difference between the two patterns Δφ = (κ 1 z − φ 1 ) − (κ 2 z − φ 2 ) is, in general z-dependent. However, when we use the approximation κ 1 = κ 2 = κ this phase difference becomes To conclude that κσ is independent of T 2 we need to show that δz 2σ at the time of observation is relatively small. In fact, in order to observe a well defined interference one needs the separation δz to be smaller than the combined width of the wavepackets, thus reliable experimental data must have δz 2σ. In figure 8 we present the product N p = 2 π κσ, found by fitting our real space data. The data is therefore consistent with the conservation law and determines κσ = π 2 5.61 = 8.81. The first term of equation (E7) is then ( 1 8.81 ) 2 ≈ 1% of the second term and can indeed be neglected. Let us now derive the form of the interference equation (E4) and identify explicit forms for the interference wavenumber κ(t) and for the width σ(t) during the evolution in free space. The general shape of two Gaussian wave functions with the same widths at any time is (the following applies to either spin state and we ignore the index i = 1, 2, for simplicity) ψ a (z, t) = A e −(z−z a ) 2 /4σ 2 + 1 2 iα(z−z a ) 2 +ik a (z−z a )+iφ a , where α(t) is a coefficient of the quadratic phase that evolves when the Gaussian expands or shrinks and is the same for the two wavepackets if the width σ(t) is common. Here also the centers z a , z b , the phases φ a , φ b and A = (2πσ 2 (t)) −1/4 are time dependent; k a , k b are (time independent) momenta of the wavepackets.
For a particle with mass m the free space evolution of the Gaussian width is where σ m is a minimum size occurring at time t = t m in the past or future history of the evolution if it occurs in free space. In our case this minimal wavepacket size occurs at some time (focusing time) after the end of the deceleration pulse (usually after the translation pulse T 3 ). The quadratic phase coefficient is given by so that α = 0 at the focusing time t = t m and α(t) → m/ (t − t m ) after a long time where σ(t) σ m . The phase of the interference term of ψ a and ψ b of equation (E10) can be written as wherez = 1 2 (z a + z b ),k = 1 2 (k a + k b ) and δk = k b − k a are the average momentum and momentum difference between the two wavepackets, while δφ = φ b − φ a . We therefore identify Finally we derive an interesting relation between the properties of the state at the observation time t f of our experiment and at the focusing time t m , where t f occurs at a time T f (time of flight) after t m , i.e. t f = t m + T f . The width from equation (E11) and wavelength from equations (E12) and (E14) are given by where d = δz(t m ) is the distance (positive or negative) between the wavepacket centers of the same spin at the focusing time and ξ = 2mσ 2 m / T f . In our experiment the wavepacket that has a larger momentum at t 0 is always further away from the chip so we can choose the indices a and b such that always z b > z a and d > 0 and hence κ(t f ) > 0 in the limit ξ 1. When the Gaussian wavepackets are much larger than their minimal size σ m the factor ξ becomes negligible and the equations for κ and σ simplify considerably. This situation indeed applies to the conditions of our experiment, as predicted by our numerical simulation (see appendix C), but let us show this here by using only direct experimental evidence. Without the focusing effect due to the positive curvature of the magnetic field potential applied during the pulses the expansion of the BEC was measured to obey the free expansion rule σ(t) = σ(0) √ 1 + ω 2 t 2 for any time t after trap release, where ω ≈ 2π × 120 Hz is the trap frequency along z in our experiment. For σ(0) = 1.2 μm [17], free expansion in our experiment would lead after a time t ≈ 16 ms (from trap release to observation) to a final size σ ∼ σ(0)ωt ≈ 15 μm. However, for the observed range of T 2 the measured width in our experiment is in the range σ(t f ) = 30-80 μm, which is at least twice as large as the expected width in free expansion. This enlarged wavepacket size indicates that the magnetic fields focus the wavepackets to a minimal size σ m that is much smaller than 15 μm and this focusing causes the enhanced expansion. As a first step let us only assume σ m < 15 μm. This puts an upper bound to ξ, as (using equation (E15)): This bound allows us to set an upper bound for σ m = T f ξ/2m < 1.7 μm (using T f < 14 ms and the mass of rubidium m = 1.44 × 10 −25 kg). Once this bound on σ m is set, we can estimate that ξ ≈ σ m /σ(t f ) < 1.7/30 1 and therefore we can take the limit ξ → 0 in equations (E15) and (E16), giving rise to the well known relation κ = md/ T f . We then obtain by multiplying the two equations d This result can also be obtained by applying the conservation law equation (E7) at t m , where As σ(t f )) we concluded above, hence equation (E17) is proved.
Having established the relations above, we can obtain estimates of the parameters of the state at the focusing time by assuming that the focusing time t m relative to the end of the translation pulse of duration T 3 is much smaller than the remaining time of flight T f , so that T f ∼ 14 ms. It then follows from κ being in the range 0.3-0.1 μm −1 that d ≈ κT f /m is in the range 3 μm to 1 μm, while σ m is in the range 0.17-0.06 μm. It then follows that ξ is in the range 0.006-0.0007. These values are reproduced in our numerical simulation (see appendix C and figure 9).
In appendix F we present an additional demonstration of the conservation laws, which does not rely on the overlap integral between the wavepackets and uses, instead, arguments based on the evolution of the phase space distribution (Wigner function) of a superposition of two wavepackets.
We have shown in this appendix that an internally coherent system, which we may define as a system having constituents from a single coherent parent that is split into a superposition of two internal states, satisfies Δθ = 0, as each of the constituents preserve the phase χ of the parent state. While this proof is based on properties of quantum superpositions, the same property may also appear in classical electromagnetic pulses (e.g. RF or microwave). When a single pulse with a measurable real field representing an electric or magnetic field, is split into a pair of pulses propagating in different trajectories and then recombined with a time delay ΔT, the resulting field is a superposition ψ(x, t) + ψ(x, t − ΔT). This field has a form similar to equation (1) of the main text with Δθ = 0. Whether such an electromagnetic field is generated by a splitting and recombination process or by an engineered electronic pulse generator, it will conserve the property Δθ = 0 when it propagates through any homogeneous medium. However, if the delay between the two pulses is generated in a dispersive medium where the group velocity v g is not equal to the phase velocity v φ , then the time delay is different for the Gaussian envelope of the pulse and for the phase of the oscillations within this envelope and therefore θ is not conserved and becomes different for the two pulses. Results similar to those obtained in this work can also be obtained in the case where Δθ = constant = 0, as discussed in appendix G.

Appendix F. Demonstration of the conservation laws in phase space
In this appendix we complement the proof of the conservation laws given in appendix E by a demonstration that indicates that the conservation of the number of interference fringes and their phase do not essentially require the overlap integral to be involved. While the proof given in appendix E is general in the sense that it shows that the quantities Γ and χ are conserved in any unitary operation, the specific form of Γ is in terms of parameters of a Gaussian wavepacket and it relies on the overlap integral between the two wavepackets, which may crucially depend on details of the shape of the wave functions at their tails. Here we present a complementary vision that does not rely on an overlap integral and may apply to wavepackets that are completely non-overlapping. On the other hand, it is based on the assumption that the potential acting during the evolution can be fairly well described by a quadratic form over the region occupied by the distribution.
As we have shown in a previous work that analyzed the interferometric sequence in our SGI [15], the evolution of the pair of wavepackets can be viewed as a scaled phase space rotation (see figure 9 in reference [15] and an additional rigorous proof here below). The form of the phase space distribution of a superposition of two wavepackets consists of two peaks at the phase space coordinates where these wavepackets are centered and a fringe pattern that appears in the Wigner distribution in between these centers, whose wave vector points perpendicular to the line that connects between the phase space center coordinates, as demonstrated in figure 12. Under rotation of phase space coordinates this structure conserves the form of this fringe pattern and in particular it conserves the number of fringes along this pattern and the phase, namely the position of the fringes relative to the center of the pattern. When the two coordinate centers are separated by momentum, while the distributions around these centers overlap in space, the fringe pattern appears as a real interference pattern in space. On the other hand, when the two centers are separated in space, such as during the focusing of the wavepackets before they expand, then the interference fringes do not appear in the real space distribution but rather in the momentum distribution.
Let us now explicitly derive the starting point result that evolution in a quadratic potential can be represented by a phase space rotation. The Wigner function for a pure state represented by a wave function ψ(x, t) in one dimension is defined as We have then proved that evolution in a quadratic Hamiltonian is equivalent to phase space rotation of the Wigner function of an arbitrary pure state. This proof can be easily extended to impure states that are given by a density matrix which is a weighted sum over density matrices of pure states.

Appendix G. Fourier analysis of the model
In this appendix we identify the periodicity wavenumber K S of the sum pattern modeled by a sum of two translated localized periodic patterns with a Gaussian profile (equation (1)). We consider the FT of this pattern with Fourier variable K up to a prefactor σ √ 2π, and using Δφ = κΔz yields equation (2) of the main text. The position at the maximum of |F V (+) (K)| is defined as K S and is given by Only when the phase difference Δφ is an integer multiple of 2π, Δφ = 2πn, equation (G2) has the trivial solution K S = κ. Otherwise the period of the sum pattern is different than the period of its constituents: the fundamental result of this work. In figure 13 we present the solution for K S /κ at the peak of the FT in equation (G2) as a function of Δφ and N p where is the number of periods of each of the constituent patterns over the range where their envelopes are larger than 1/e 2 of their maxima. The value of N p , as for κσ, is conserved, i.e. independent of T 2 in our experimental setup. The deviation of K S from κ is larger when the number of periods within the system size is small, while it diminishes with N p so that in the limit of an infinite periodic system, N p → ∞, the period K S of the sum pattern is just the period of the constituent patterns κ.
G.1. Jumps K S has in particular two degenerate solution at Δφ = π(2n + 1) at K S = κ ± 1 2 ΔK S where the jump ΔK S satisfies ΔK S = π(2n + 1) with solutions shown in the inset to figure 4. Figure 14. Mutual cancellation of the variation of system size (2σ) and relative translation Δz leading to the rigidity of sum pattern periodicity in the experiment: numerical results (see appendix C). (a) Average wave-packet size σ increases with T 2 due to increased focusing strength and a resulting increased expansion speed during TOF, while the translation distance decreases with T 2 due to a weaker final gradient pulse given further away from the chip. (b) The root sum of squares √ (2σ) 2 + Δz 2 varies relatively slightly over the range shown, giving rise to a minor deviation from zero slope of K S as a function of T 2 at the center of each of the plateaus.

G.2. Rigidity
We can gain some insight into the formation of plateaus (rigidity of K S between jumps) by looking analytically at the derivative of equation (G2) with respect to Δφ in between the jumps, i.e. at Δφ = 2πn. A plateau requires that this derivative vanishes, This requirement can be satisfied for each n if where κ 0 is a constant. By using the relations Δφ = κΔz and N p = 2κσ/π equation (G5) is satisfied if In figure 14 we show the dependence of Δz and 2σ on the deceleration time T 2 over the range shown in figures 2 and 4. The reason that Δz depends on T 2 is that with a longer T 2 the atoms spend more time before reaching T 3 and also acquire a larger momentum, hence at T 3 the atoms are further away from the chip and experience a smaller gradient. On the other hand the mean wavepacket width σ increases with T 2 because the focusing strength and hence the speed of expansion is larger for longer T 2 . The value of √ (2σ) 2 + Δz 2 in this range has therefore a relatively low variation (standard deviation of 2.4% of the mean). This explains why K S shows rigidity over the range of T 2 variation.
Let us summarize. In our experiment the rigidity, namely insensitivity of K S to the change of the periodicity κ of the constituent patterns, is achieved due to the following properties of the model: • θ 1 = θ 2 .
• κσ = 1 2 πN p constant. • (2σ) 2 + (Δz) 2 constant. While the first two conditions are exact constraints that follow from fundamental properties of the system, as shown in appendix E, the third condition is an approximate result related to the specific choice of the experimental parameters. It is a necessary condition for obtaining zero slopes in the middle of each plateau but not a sufficient condition for obtaining strict flatness over the whole range of each plateau.
In the case of non-overlapping constituent patterns, i.e. Δz 2σ, the third condition for rigidity simplifies to Δz = constant. In this case we can also relax the first two conditions (in an arbitrary system that does not inherently satisfy them): it is sufficient that Δθ is a constant (not necessarily zero) and κσ is not required to be conserved. Let us consider this scenario where Δz is kept constant. In our experimental simulation, a constant Δz (independent of T 2 ) may be realized by adjusting delay times between the pulses. In the context of real electromagnetic pulses, it may be achieved by picking a sub-sample satisfying Δz = constant (or Δt = constant) out of a random flux or it can be generated artificially (see discussion of a mode-locked laser below). In such a case, equation (G2) is written as The solution for K S is plotted in figure 15, clearly there are plateaus at large Δφ. These plateaus become visible at Δφ N p , equivalent to Δz σ. This implies that the two patterns are separated and do not interfere in real space. However, Δφ and the sum pattern wavevector are well defined. We note in particular that the plateaus are at universal values K S Δz = 2πn where n is defined at the center of the plateau as Δφ = 2πn.
We also plot in figure 15 a case with a finite Δθ = θ 2 − θ 1 that results in shifting the argument of the tan function in equation (G8) by −Δθ. The orange-dashed line corresponds to Δθ = π/4, showing a similar structure except that the jump positions and plateau values are shifted.
To understand the rigidity in this case, consider the AFT of the pair of patterns |F V (+) (K)| = e − 1 2 σ 2 (K−κ) 2 | cos( 1 2 KΔz − Δθ)| is a product of a Gaussian depending on κ and a cosine that does not vary when κ is scanned. This AFT is an absolute cosine with Fourier space periodicity 2π/Δz with a Gaussian envelope whose width κ and extends over 1/σ extends over a few periods of the cosine and its center varies with κ. When κ is changed the peaks of the cosine stay static and K S , defined as the largest peak of the AFT, changes only when the peak of the Gaussian function is closer to the next cosine peak.
This model of non-overlapping pulses with Δθ = constant may be implemented with a mode-locked laser, whose output is a train of equally spaced light pulses with a constant phase shift Δθ (so-called carrier-envelope offset) between each pulse and the next one, which can be locked to a value that is independent of other laser parameters. By analyzing the FT of pairs of consequent pulses of this laser while scanning the carrier frequency we can mimic at least some of the features of our model, namely the main peak of the spectrum of each pair will show jumps when the carrier frequency is scanned. In order to obtain exactly the same behavior of K S as in figure 15, where its value is constant within a range of frequencies (i.e., rigidity), giving rise to a quantized spectrum of K S , one has to manipulate the properties of the output such that the time Δt between consequent pulses is also independent of carrier frequency. This is not a straight-forward goal and requires careful engineering of the dispersive properties of the laser medium. Obviously, such tests may only be done with frequencies for which available detectors have enough bandwidth to follow the oscillations.