Universal atom interferometer simulation of elastic scattering processes

In this article, we introduce a universal simulation framework covering all regimes of matter-wave light-pulse elastic scattering. Applied to atom interferometry as a study case, this simulator solves the atom-light diffraction problem in the elastic case, i.e., when the internal state of the atoms remains unchanged. Taking this perspective, the light-pulse beam splitting is interpreted as a space and time-dependent external potential. In a shift from the usual approach based on a system of momentum-space ordinary differential equations, our position-space treatment is flexible and scales favourably for realistic cases where the light fields have an arbitrary complex spatial behaviour rather than being mere plane waves. Moreover, the solver architecture we developed is effortlessly extended to the problem class of trapped and interacting geometries, which has no simple formulation in the usual framework of momentum-space ordinary differential equations. We check the validity of our model by revisiting several case studies relevant to the precision atom interferometry community. We retrieve analytical solutions when they exist and extend the analysis to more complex parameter ranges in a cross-regime fashion. The flexibility of the approach, the insight it gives, its numerical scalability and accuracy make it an exquisite tool to design, understand and quantitatively analyse metrology-oriented matter-wave interferometry experiments.

The commonly used approach for treating light-pulse beam-splitter and mirror dynamics in matter-wave systems consists in solving a system of ordinary differential equations (ODE) with explicit couplings between the relevant momentum states.
This formulation starts by identifying the relevant diffraction processes and extracting their corresponding coupling terms in the ODE 1,2 . In the elastic scattering case, each pair of light plane waves can drive a set of twophoton transitions from one momentum class j to the next neighboring orders j ± 2 . The presence of multiple couplings allows for higher order transitions and the system is simplified by choosing a cutoff omitting small transition strengths. This ODE approach works well for simple cases leading to analytical solutions in the deep Bragg and Raman-Nath regimes 1,2 . Using a perturbative treatement, it was generalised to the intermediate, so-called quasi-Bragg regime 3 . A numerical solution in this regime has been extended in the case of a finite momentum width 4 . In a different approach, Siemß et al. 5 developed an analytic theory for Bragg atom interferometry based on the adiabatic theorem for quasi-Bragg pulses. Realistically distorted light beams or mean-field interactions, however, sharply increase the number of plane wave states and their couplings required for an accurate description. The formulation of the ODE becomes increasingly large and inflexible, with a set of coupling terms for each relevant pair of light plane waves.
Here, we take an alternative approach and solve the system in its partial differential equation (PDE) formulation following the Schrödinger equation. This time-dependent perspective 6 has several advantages in terms of ease of formulation and implementation, flexibility and numerical efficiency for a broad range of cases. Indeed, this treatment is valid for different types of beam splitters (Bloch, Raman-Nath, deep Bragg and any regime in between) and pulse arrangements. Combining successive light-pulse beam-splitters naturally promoted our solver to a cross-regime or universal atom interferometry simulator that could cope with a wide range of non-ideal effects such as light spatial distortions or atomic interactions, yet being free of commonly-made approximations incompatible with a metrological use.
The position-space representation seems underutilised in the treatment of atom interferometry problems in favor of the momentum-space description although several early attempts of using it were reported for specific Scientific Reports | (2020) 10:22120 | https://doi.org/10.1038/s41598-020-78859-1 www.nature.com/scientificreports/ cases [7][8][9][10][11] . In this paper, we show the unique insights this approach can deliver and, contrary to widespread belief, its great numerical precision and scalability. In addition we illustrate our study with relevant examples from the precision atom interferometry field.

Theoretical model
Light-pulse beam splitting as an external potential. We start with a semi-classical model of Bragg diffraction, where a two-level atom is interacting with a classical light field 1,2 . This light field consists of a pair of two counter-propagating laser beams realised by a retro-reflection mirror setup for example. Assuming that the detuning of the laser light is much larger than the natural line width of the atom, one may perform the adiabatic elimination of the excited state. This yields an effective Schrödinger equation for the lower-energy atomic state ψ(x, t) with an external potential proportional to the intensity of the electric field with the two-photon Rabi frequency and wave vector k = 2π/ in a simplified 1D geometry along the x-direction. For the present study, we consider a 87 Rb atom that is addressed at the D2 transition with = 780 nm resulting in a recoil frequency and velocity 12 of ω r = k 2 /2m = 2π · 3.8 kHz and v r = k/m = 5.9 mm/s, respectively.
In the context of realistic precision atom interferometric setups, it is necessary to include Rabi frequencies �(x, t) and wave vectors k(x, t) which are space and time-dependent. This allows one to account for important experimental ingredients such as the Doppler detuning or the beam shapes including wavefront curvatures [13][14][15] and Gouy phases [16][17][18][19] . Moreover, this generalisation allows one to effortlessly include the superposition of more than two laser fields interacting with the atoms as in the promising case of double Bragg diffraction [20][21][22] , and to model complex atom-light interaction processes where spurious light reflections or other experimental imperfections are present 23 .
Atom interferometer geometries. The light-pulse representation presented in the previous section is the elementary component necessary to generate arbitrary geometries of matter-wave interferometers operating in the elastic diffraction limit. Indeed, since the atom-light interaction in this regime conserves the internal state of the atomic system, a scalar Schrödinger equation is sufficient to describe the physics of the problem in contrast to the model adopted in Ref. 10 .
For example, a Mach-Zehnder-like interferometer geometry can be generated by a succession of π 2 − π − π 2 Bragg pulses (beam-splitter, mirror, beam-splitter pulses) of order n separated by a free drift time of T between each pair of pulses. In the case of Gaussian temporal pulses, this leads to a time-dependent Rabi frequency where bs , τ bs and m , τ m are the peak Rabi frequencies and their respective durations associated to the beamsplitter and mirror pulses, respectively. We numerically solve the corresponding time-dependent Schödinger equation using the split-operator method 24 to propagate the atomic wave packets along the two arms. The populations in the two output ports |+� = |0 k� and |−� = |2n k� are evaluated after the last recombination pulse waiting for a time of flight τ ToF long enough that the atomic wave packets spatially separate. They are obtained by the integration where the integration domains extend over a space interval with non-vanishing probability density of the states |±� . These probabilities are further normalised to account for the loss of atoms to other parasitic momentum classes Using Feynman's path integral approach, the resulting phase shift between the two arms can be decomposed as 25,26 The propagation phase is calculated by evaluating the classical action along the trajectories of the wave packet's centers. The laser phase corresponds to the accumulated phase imprinted by the light pulses at the atom-light interaction position and time. Finally, the separation phase is different from zero if the final wave packets are not overlapping at the time of the final beam splitter, t = 2T.
To extract the relative phase �φ between the two conjugate ports and the contrast C, one can scan a laser phase φ 0 ∈ [0, 2π] at the last beam splitter and evaluate the populations 1 varying as (5) �φ = �φ propagation + �φ laser + �φ separation .
(1 ± C cos(�φ + nφ 0 )). www.nature.com/scientificreports/ The resulting fringe pattern is then fitted with �φ and C ≤ 1 as fit parameters. This method, analogous to experimental procedures, allows one to determine the relative phase modulo 2π.

Results
Raman-Nath beam splitter. The Raman-Nath regime, characterised by a spatially symmetric beam splitting, is the limit of elastic diffraction for very short interaction times of τ ≪ 1 √ 2�ω r . The dynamics of the system can, in this case, be analytically captured following Refs. 1,2 where g n (t) describes the amplitude of the momentum state |2n k� and J n the Bessel functions of the first kind. Such experiments are at the heart of investigations as the one reported in Ref. 27 where a Raman-Nath beam splitter was used to initialise a three-path contrast interferometer offering the possibility of measuring the recoil frequency ω r .
To demonstrate the validity of our position-space approach, we contrast our results to the analytical ones obtained adopting the parameters of Ref. 27 . Figure 1 shows the outcome of a symmetric Raman-Nath beam splitter targeting the preparation of three momentum states: 50% into |0 k� and 25% in each of the | ± 2 k� momentum classes. As a feature of our solver, we directly observe the losses to higher momentum states ( p = ±4 k and p = ±6 k ) due to the finite pulse fidelity. An excellent agreement is found with the analytical predictions (green filled circles) of the populations of the momentum states.

Bragg-diffraction Mach-Zehnder interferometers.
To simulate a Mach-Zehnder atom interferometer based on Bragg diffraction, we consider a pair of two counter-propagating laser beams with a relative frequency detuning �ω = ω 1 − ω 2 = 2nkv r and a phase jump φ 0 ∈ [0, 2π] . This gives rise to the following running optical lattice For sufficiently long atom-light interaction times, i.e. in the quasi-and deep-Bragg regimes 2,3,28,29 , the driven Bragg order n with momentum transfer p = 2n k is determined by the relative frequency detuning �ω of the two laser beams. The relative velocity between the initially prepared atom and the optical lattice is v = nv r . In the rest frame of the optical lattice, the atom has a momentum p = −n k . The difference of kinetic energy between the initial ( p = −n k ) and target state ( p = +n k ) is vanishing and therefore this transition is energetically allowed and leads to a �p = n k − (−n k) = 2n k momentum transfer. www.nature.com/scientificreports/ We now realise beam splitters and mirrors by finding the right combination of peak Rabi frequency and interaction time (�, τ ) , either by numerical population optimisation or analytically, when we work in the deep Bragg regime. Recent advances by Siemß et al. 5 generalise this to the quasi-Bragg regime in an analytical description of Bragg pulses based on the adiabatic theorem. For the pulses used in this paper, the two approaches give the same result for the optimised Rabi frequencies and pulse durations.
In Fig. 2, we simulate a Mach-Zehnder geometry and illustrate the diffraction outcome by showing a spacetime diagram of the density distribution |ψ(x, t)| 2 . For the parameters chosen here, a clear feature of the dynamics is the appearance of additional atomic channels after the mirror pulse, which can be attributed to the velocity selectivity arising from a pulse with a finite duration characterised by τ . The finite velocity acceptance can, indeed, be estimated over the Fourier width σ f of the applied pulse as with σ f = 1/(2πτ ) and f being the frequency variable. This yields the velocity acceptance 30 With an initial velocity width of the atomic probability distribution of σ atom v = 0.1 v r , it is clear that velocity components with |v| = σ pulse v will have a much smaller excitation probability than the components at the centre of the cloud, which leads to the characteristic double well densities of the parasitic trajectories.
With momenta p upper = 0 k and p lower = 2 k , both parasitic trajectories still fulfill the resonance condition with the final Bragg beam splitter, which leads to the emergence of ten trajectories after the exit beam splitter. For a measurement in position space, it is now important that a sufficiently long time of flight τ ToF is applied that the ports of the Mach-Zehnder interferometer do not overlap with the parasitic ports and bias the relative phase measurement. For large densities, the parasitic trajectories at the Mach-Zehnder ports should not overlap since this may already lead to density interaction phase shifts H int ∝ |ψ(x, t)| 2 . To circumvent these problems it , for strongly suppressed parasitic trajectories due to velocity selectivity.
Implementing high-order Bragg diffraction is a natural avenue to increase the momentum separation of an atom interferometer, and therefore its sensitivity. In Fig. 3, we run our solver to observe the population distribution across the different ports of a Mach-Zehnder configuration with Bragg orders up to n = 3 . This is done in a straightforward way by scanning the laser phase φ 0 . We fit the data points corresponding to the population in the fast port |2 k� for the different Bragg orders according to Eq. (6) and observe a clear sinusoidal signal of the simulated fringes, as expected. The resulting contrasts and phase shifts are directly found by our theory model and numerical solver which include the ideal phase shifts commonly found 25,26 , and go beyond to comprise several non-ideal effects as (i) finite momentum widths, (ii) finite pulse timings and (iii) multi-port Bragg diffraction 37,38 and the resulting diffraction phase. The natural occurrence of these effects and the possibility to quantify them are a native feature of our simulator.

Symmetric double Bragg geometry. Scalable and symmetric atom interferometers based on double
Bragg diffraction were theoretically studied 22 and experimentally demonstrated 21 . This dual-lattice geometry has particular advantages, including an increased sensitivity due to the doubled scale factor compared to single-Bragg diffraction, as well as an intrinsic suppression of noise and certain systematic uncertainties due to the symmetric configuration 21 . Combining this technique with subsequent Bloch oscillations applied to the two interferometer arms led to reaching momentum separations of thousands of photon recoils as was recently shown in Ref. 23 .
In double Bragg diffraction schemes, two counter-propagating optical lattices are implemented in such a way that the recoil is simultaneously transferred in opposite directions, leading to a beam splitter momentum separation of p = 4n k 21,22 . To extend our simulator to this important class of interferometers, we merely have to add a term to the external potential The procedures of realising a desired 4n k momentum transfer, as well as mirror or splitter pulses, are identical to the case of single-Bragg diffraction. A simple scan of the Rabi frequency and pulse timings was enough to obtain a full double Bragg interferometer as shown in Fig. 4. The different resulting paths are illustrated in this space-time diagram of the density distribution |ψ(x, t)| 2 . Similarly to the single-Bragg Mach-Zehnder interferometer, we observe additional parasitic interferometers due to the finite velocity filter of the Bragg pulses after the mirror pulse of the interferometer. Due to a finite fidelity of the initial beam splitter, some atoms remain in the |0 k� port and recombine at the last beam splitter with the trajectories of the interferometer. In a metrological The initial momentum width of the atomic sample is σ p = 0.01 k. The corresponding Rabi frequencies for the higher order Bragg transitions were found by optimising for an ideal 50 : 50 population splitting of the π 2 pulse. This leads to 4 k = 3.7 ω r and 6 k = 8.4 ω r . The Rabi frequency for the 2 k transition is 2 k = 1.0573 ω r . The solid lines are the respective fringe scan fits from which the phase shifts and contrasts are directly extracted.

Gravity gradient cancellation for a combined Bragg and Bloch geometry. Precision atom inter-
ferometry-based inertial sensors are sensitive to higher order terms of the gravitational potential, including gravity gradients. In particular, for atom interferometric tests of Einstein's equivalence principle (EP), gravity gradients pose a challenge by coupling to the initial conditions, i.e. position and velocity of the two test isotopes 39 . A finite initial differential position or velocity of the two species can, if unaccounted for, mimic a violation of the EP. By considering a gravitational potential of the form where Ŵ = Ŵ xx is the gravity gradient in the direction normal to the Earth's surface, the relative phase of a freely falling interferometer can be calculated as 40 with k eff = 2nk. In Ref. 40 , it was shown that introducing a variation of the effective wave vector �k eff = Ŵk eff T 2 /2 at the π pulse can cancel the additional phase shift due to the gravity gradient. This was experimentally demonstrated in Refs. 41,42 . The corresponding probability density |ψ(x, t)| 2 is plotted for an initial momentum width of σ p = 0.1 k. The timings of the Gaussian splitter and mirror pulses are set to τ bs = 25 µ s and τ m = 50 µ s, respectively. The corresponding Rabi frequencies are found by optimising the desired population transfer. The first π 2 pulse corresponds to a 2 k transfer in two directions, realised by two counter-propagating optical lattices which results in a 4 k separation between the two interferometer arms. The mirror pulse is a 4 k Bragg transition with a Rabi frequency of = 1.9 ω r such that both arms make a transition from | ± 2 k� → | ∓ 2 k� . The last recombination pulse now realises a 50 : 50 split of the upper trajectory to | − 2 k� and |0 k� and the lower trajectory to | + 2 k� and |0 k� . This leads to a final population of 25% in the | ± 2 k� ports and 50% in the |0 k� port. The separation time between the pulses is T = 10 ms with a final time of flight after the exit beam splitter of τ ToF = 20 ms. Due to velocity selectivity of the Bragg pulses and a non-ideal fidelity of the initial beam splitter pulse, several parasitic interferometers can be observed.

Scientific Reports
| (2020) 10:22120 | https://doi.org/10.1038/s41598-020-78859-1 www.nature.com/scientificreports/ The same principle applies to the gradiometer configuration of left panel in Fig. 5 where the effect of a gravity gradient is compensated by the application of a wave vector correction. This is reminiscent of another experimental cancellation of the gravity gradient phase shifts 41 . In our example, we first consider a set of two Mach-Zehnder interferometers vertically separated by h = 2 m, realised with 4 k Bragg transitions where the atoms start with the same initial velocities v 0 . Choosing a Doppler detuning according to a Bragg = g , the gradiometric phase reads By scanning the momentum of the applied π pulse, one can compensate the gradiometric phase. This is observed in our simulations at the analytically predicted value of �k eff = Ŵk eff T 2 /2 (red dashed curve crossing the zero horizontal line).
It is particularly interesting to use our simulator to find this correction phase in the context of more challenging situations, such as a combined scalable Bragg and Bloch Mach-Zehnder interferometer or a symmetric Bloch beam splitter 43 where analytic solutions are not easily found.
Bloch oscillations can be used to quickly impart a momentum of p = 2n Bloch k on the atoms 44,45 . This adiabatic process can be realised by loading the atoms into a co-moving optical lattice, then accelerating the optical lattice by applying a frequency chirp and finally by unloading the atom from the optical lattice. In our model, this corresponds to the following external potential www.nature.com/scientificreports/ where τ load , τ chirp and τ unload are the durations of the lattice loading, acceleration and unloading, respectively. By ramping up the co-moving optical lattice, the atoms are loaded into the first Bloch band with a quasimomentum q = 0 . An acceleration of the optical lattice acts as a constant force on the atoms which linearly increases the quasimomentum over time. When the criterion for an adiabatic acceleration of the optical lattice is met, the atoms stay in the first Bloch band and undergo a Bloch oscillation, which can be repeated n Bloch times leading to a final momentum transfer of p = 2n Bloch k.
The π pulse correction �k eff = Ŵk eff T 2 /2 is proportional to the space-time area A Bragg = k eff T 2 /m of the underlying 2n k Mach-Zehnder geometry and does not compensate the gravity gradient effects in the Bloch case. Analysing the space-time area A Bragg+Bloch immediately shows a non-trivial correction compared to A Bragg . The suitable momentum compensation factor is, however, found using our solver at the crossing of the dashed blue line and the vertical zero limit ( k Bragg+Bloch eff = 0.932 k Bragg eff ). This straightforward implementation of our toolbox in a rather complex arrangement is promising for an extensive use of this framework to design, interpret or propose advanced experimental schemes.
Trapped interferometry of an interacting BEC. Employing Bose-Einstein condensate (BEC) sources 46,47 for atom interferometry 34,48,49 has numerous advantages such as the possibility to start with very narrow momentum widths σ p 31-36 , which enables high fidelities of the interferometry pulses 4 . For interacting atomic ensembles, it is necessary to take into account the scattering properties of the particles. The Schrödinger equation is not anymore sufficient to describe the system dynamics and the ODE approach becomes rather complex to use as shown in the section on scalability and numerics. We rather generalise our position-space approach and consider a trapped BEC atom interferometer including two-body scattering interactions described in a meanfield framework. The corresponding Gross-Pitaevskii equation reads 50 where the quantum gas of N atoms is trapped is a quasi-1D guide aligned with the interferometry direction and characterised by a transverse trapping at an angular frequency ω ⊥ much stronger than the longitudinal one. These interactions can effectively be reduced in 1D to a magnitude of g 1D = 2 a s ω ⊥ . For our calculation, we set the s-wave scattering length of 87 Rb to one Bohr radius, i.e. a s = a 0 = 5.3 × 10 −11 m. Experimentally, such a value can be realised using a Feshbach resonance technique 51 . This model is well valid in the weakly interacting limit, i.e. when a s N|ψ| 2 ≪ 1 52,53 .
All atom interferometric considerations mentioned earlier, like the Bragg resonance conditions, construction of interferometer geometries, the implementation of Doppler detunings, phase calculations and population measurements are also valid in this case without any extra theoretical effort. The non-linear Gross-Pitaevskii equation is solved following the split-operator method as in the Schrödinger case 24 .
If the atom interferometer is perfectly symmetric in the two directions of the matterwave guide, no phase shift should occur. In realistic situations, however, the finite fidelity of the beam splitters creates an imbalance δN of the particle numbers between the two interferometer arms. The phase shift in this case can be related to the differential chemical potential by We illustrate the capability of our approach to quantitatively predict this effect by contrasting it to the well-known treatment of this dephasing. Following Ref. 48 , we introduce δN = 0 and analyse the dephasing by using the 1D Thomas-Fermi chemical potential of the harmonic oscillator potential The + and − signs refer here to the arms 1 and 2, respectively. We assume the Thomas-Fermi radii before and after the atom-light interaction to be approximately the same. To this end, one needs to introduce the correction factor of 1/ √ 2 which is a direct consequence of Using Eq. (21), one finds a phase shift of  where L denotes the half-width of the BEC, simply by replacing L = R TF , as assumed by Ref. 48 . In Fig. 6a, the mean-field shift is plotted as a function of the atom number imbalance in the two cases of the numerical solution of the Gross-Pitaevskii equation and with the analytical model using the Thomas-Fermi approximation. It is worth noting that the dephasing is accompanied by a loss of contrast consistent with previous theoretical studies 54 . We performed a numerical optimisation to find the maximal particle number N up to which we find a contrast of C > 99 % , which is N ≤ 6 · 10 4 in this case. In Fig. 6b, the absolute value of the difference between the numerical and the analytical solutions of �φ is plotted. For an imbalance of the order of 10 % , we observe an agreement at the mrad level. We could point to different possible sources for the relative phase difference. First, the assumptions of the Thomas-Fermi approximation at the heart of the analytical method are not necessarily satisfied here with N ≤ 6 · 10 4 . Moreover, the analytical treatment neglects all time-dependent effects occurring during the light-atom interactions at the mirror and beam-splitter pulses. These effects, combined with a nonvanishing mean-field would lead to additional phase shifts and shape deformations of the wave functions that are absent from a simple Thomas-Fermi assumption.

Scalability and numerics
Numerical accuracy and precision. To gain a better understanding of the numerical accuracy of the simulations, we plot in Fig. 7 the dependency of the phase shift |�φ| on the momentum width of the atomic sample σ p for a 2 k Bragg Mach-Zehnder interferometer. We study two realisations which differ in the peak Rabi frequency with corresponding pulse lengths to perform beam splitter and mirror pulses. For both cases we observe a similar characteristic qualitative behaviour of |�φ| scaling with σ p . Going to smaller initial momentum widths systematically decreases the phase shift until it reaches a plateau of 1 × 10 −7 rad for � = 1.06ω r and 2.5 × 10 −14 for � = 0.53ω r . www.nature.com/scientificreports/ This qualitative behaviour can be explained by considering the effect of parasitic trajectories. In Fig. 2 it is clearly visible that after the time of flight of τ ToF = 2T , there is no clear separation between the parasitic trajectories and the main ports of the Mach-Zehnder interferometer, which leads to interference between them. We choose the integration borders by setting up a symmetric interval around the peak value of each of the ports (see Eq. (3)), ensuring a minimal influence of the parasitic atoms on the interferometric ports. Nevertheless, the interference between the interferometric ports and the parasitic trajectories modifies the measured particle number and therefore also the inferred relative phase. This effect decreases with smaller initial momentum width since less atoms populate the parasitic trajectories overlapping with the main ports, which explains the decrease of relative phase |�φ| between σ p = 0.1 k to σ p = 0.05 k ( = 1.06 ω r ) and σ p = 0.03 k ( = 0.53 ω r ). Another important contribution to the relative phase |�φ| which is not captured by Feynman's path integral approach 25,26 is the diffraction phase, which is fundamentally linked to the excitation of non-resonant momentum states 37,38 . Using smaller Rabi frequencies leads to a reduced population of non-resonant momentum states (after a beam splitter pulse we find P(−2 k, � = 1.06 ω r ) + P(4 k, � = 1.06 ω r ) = 1.3 × 10 −7 and P(−2 k, � = 0.53 ω r ) + P(4 k, � = 0.53 ω r ) = 1.9 × 10 −18 ) and therefore to a reduced diffraction phase which explains that operating a Mach-Zehnder interferometer at = 0.53 ω r leads to a much smaller residual phase shift than at = 1.06 ω r .
These results indicate that our simulator reaches at least a relative phase accuracy at the level of 2.5 × 10 −14 rad. It is worth mentioning, that the numerical parameters chosen to reach this performance are very accessible on modestly powerful desktop computers. The computation took τ CPUtime = 12.7 s on an Intel Xeon X5670 processor using four cores (2.93 GHz, 12 MB last level cache). Modeling precision atom interferometry problems with this method is therefore a practical, flexible and highly accurate approach. Using improved resolutions in position and time or higher order operator splitting schemes 55 leads to even better numerical precision and accuracy.
Numerical convergence. To analyse the numerical convergence as well as the connected numerical precision and accuracy of the split-operator method applied to the previously presented systems, we simulate three different interferometer settings for different space and time grids. We consider a 2 k Mach-Zehnder interferometer, a 2 k Mach-Zehnder interferometer in a waveguide and as a last example a (2 + 2) k Bragg+Bloch Mach-Zehnder interferometer in order to quantify the necessary resolution and grid sizes that can be derived from these results. In Figs. 8 and 9 we extract the relative phase over successively decreasing spatial and temporal steps and we compare them to simulations with sufficiently fine spatial and temporal resolutions by plotting the  Figure 7. Phase shift of a 2 k Mach-Zehnder interferometer as a function of the initial momentum width of an atomic sample. We evaluate the phase shift for pulse lengths of τ bs = 25 µ s and τ m = 50 µ s (blue dots) and for τ bs = 50 µ s and τ m = 100 µ s (red dots), using peak Rabi frequencies of 25µs = 1.06 ω r and 50µs = 0.53 ω r . The dashed lines are a guide to the eye. We find a systematic decreasing behaviour of the relative phase offset |�φ| starting from an initial momentum width of σ p = 0.1 k (far right) to σ p = 0.05 k (red dots) and σ p = 0.03 k (blue dots). Reaching those critical initial momentum widths both curves show fixed relative phase offsets |�φ| , which in the case of the interferometer with smaller Rabi frequency of � = 0.53ω r (red dots) reaches a value of 2.5 × 10 −14 rad (see text). The numerical simulations were performed with 65,536 grid points, an interaction time step of dt int = 1 µ s and a free evolution time step of dt free = 10 µ s, leading to a computational time of τ CPUtime = 12.69 s on four cores of an Intel Xeon X5670 processor with 2.93 GHz frequency and 12 MB of cache.
For the presented cases we compare to dt = 4 ns and dx = 0.01 . The choice of the steps fulfilling the necessary resolution is motivated in the following, relating them to the physical quantities of the problem (optical lattice and atomic wavepacket). The fast Fourier transform (FFT) efficiently switches between momentum and position representations to apply kinetic and potential propagators. The corresponding position and momentum grids are defined by the number of grid points N grid and the total size of the position grid x as where dp and dx are the steps in momentum and position, respectively, and p the total size of the momentum grid.
To resolve a finite momentum width of the atomic cloud we are restricted to ( = 780 nm) which sets a bound to the size of the position grid. Finally, x has to be chosen according to the maximal separation of the atomic clouds x sep . With this we find To include all momentum orders necessary to simulate the considered atom interferometric sequences, we are naturally bound by Hence, we find that , dp = 2π �x and �p = 2π dx , x x sep ≫ . Bragg MZ, 0hk port Bragg MZ, 2hk port Bragg MZ, waveguide, 2hk port Bragg+Bloch MZ, 2hk port www.nature.com/scientificreports/ which is the natural condition imposed by the necessity of resolving the atomic dynamics in the optical lattice nodes and anti-nodes of the Bragg and Bloch beams. The absence of data points of the (2 + 2) k Bragg+Bloch Mach-Zehnder interferometer in Fig. 8 shows the limits given by Eq. (31). Choosing position steps at dx = 0.07 leads to a maximal computed momentum of ± 7.1 k , which results in the impossibility to find probabilities at 8 k . For this specific interferometer, however, it is critical to resolve those momenta, since they are residually populated during the atom-light interaction processes. Imposing that the position step is roughly one order of magnitude smaller than the wavelength ( dx 0.06 ) results in a reasonable momentum truncation and resolution of the light potential and therefore in the convergence of the numerical routine. Additionally, we can observe that reaching spatial resolutions of dx = 0.01 , one finds, for all three studied cases, satisfying numerical accuracy and precision which in the worst case is approximately 1 × 10 −9 rad for the (2 + 2) k Bragg+Bloch Mach-Zehnder interferometer.
The typical time scales we need to consider are set on the one hand by the velocities of the optical lattice beams and the atomic cloud, and on the other hand by the duration of the atom-light interaction τ . The beams, as well as the atomic cloud, move with velocities which are proportional to the recoil velocity v r . Given that we want to drive Bragg processes of the order of n, we find the following bound on the time step dt The typical duration of a pulse in the quasi-Bragg regime is typically adapted to the momentum width due to the spectral properties of the finite pulse. Here, we assume a lower bound of τ = 10 µ s, which leads to dt < τ . It is worth noting that this time step is only necessary during the atom-light interaction. One can simulate the free evolution between the pulses with a much larger time step (without external and interaction potentials a single step suffices) or using scaling techniques 34,[56][57][58][59] . Figure 9 shows that depending on the specific form of the simulated light potential or the consideration of two-particle interactions, we observe a characteristic convergence behaviour which we directly connect to the propagation error of the split-operator routine 55,60 . For the 2 k Mach-Zehnder interferometer one can observe the already found level of convergence around 1 × 10 −13 rad (see Fig. 7). Interestingly, the 2 k port reaches a level of 1 × 10 −13 rad at a time step of dt ≈ 0.1 µ s, whereas the 0 k port already convergences to that level at a time step of dt = 1 µ s. Note, that the diffraction phase and therefore the relative phase in the slow port vanishes but reaches a finite value of approximately 1 × 10 −7 rad in the fast port 37,38 (see Fig. 7). The next analysed case is the 2 k Mach-Zehnder interferometer in a waveguide whereby introducing the non-linear interaction term (see (31)   Bragg MZ, 0hk port Bragg MZ, 2hk port Bragg MZ, waveguide, 2hk port Bragg+Bloch MZ, 2hk port Figure 9. Numerical convergence analysis of three different interferometer realisations given by a 2 k Mach-Zehnder interferometer (blue and orange dots corresponding to the 0 k and 2 k ports), a 2 k Mach-Zehnder interferometer in a waveguide (green dots) and a (2 + 2) k Bragg+Bloch Mach-Zehnder interferometer (red dots). We analyse the numerical convergence behaviour when changing the temporal step dt of the numerical simulation using the third-order split-operator method with a spatial step of dx = 4 × 10 −2 . The same parameters as Fig. 8 (20)) one can observe a more demanding convergence behaviour, leading to an initial precision of 1 µrad at dt = 1 µ s, which converges to 1 × 10 −10 rad at a time step of dt ≈ 10 −8 s. Introducing additional Bloch oscillations shifts the convergence curve again by two orders of magnitude at a minimal time step of dt = 1 µ s and reaches a level of approximately 1 µrad at dt = 10 −8 s. Note that the precision and accuracy of the split-operator algorithm strongly depends on the potential that is simulated 60 and that in the case of a Bloch oscillation the optical potential linearly changes its velocity, where in the case of a Bragg transition the optical lattice moves with a constant velocity during the atom-light interaction process. Additionally, the atom-light interaction time of a Bloch oscillation is typically one order of magnitude larger compared with a Bragg transition which explains the need for finer temporal grids in order to achieve reasonable precision and accuracy.
Time complexity analysis. In this section, we compare the time complexity behaviour of the commonlyused method of treating the beam splitter and mirror dynamics given by the ODE approach with the PDE formulation presented in this paper, based on a position-space approach to the Schrödinger equation. To assess the time complexity of the ODE treatment, we re-derive it from the Schrödinger equation We decompose the wave function in a momentum state basis as done in Refs. [2][3][4] where j denotes the momentum orders considered and δ the discrete representation of momenta in the interval [k j − k/2, k j + k/2] which captures the finite momentum width of the atoms around each momentum class k j . Making the two exponential terms appear in cos 2 (kx) , one obtains which is a set of N eq coupled ordinary differential equations. This number N eq of equations to solve is equal to N j N δ , set by the truncation condition restricting the solution space to N j momentum classes, each discretised in N δ sub-components. Using standard solvers for such systems as Runge-Kutta, multistep or the Bulirsch-Stoer methods 61 , we generally need to evaluate the right hand side of the system of equations over several iterations. With N eq differential equations, where each one has only two coupling terms, one finds a time complexity of O(N eq ).
In Fig. 10 we present a visualisation of different possible momentum couplings starting from 0 k to other momentum components. One starts with only three momentum states and two coupling elements in Fig. 10a corresponding to a vanishing momentum width ( δ = 0 in Eq. (35)). If the momentum width is introduced (Fig. 10b), the number of coupling elements increases since every δ sub-momentum class of j is connected to the same sub-momentum class of j − 2 and j + 2 as suggested by Eq. (35). In order to reduce visual complexity we are only showing couplings that start from the 0 k wavepacket, while dropping coupling elements starting from ±2 k . We also fixed the number of momentum states per integer momentum class k j to three, which in a realistic example is at least an order of magnitude larger.
In a next step, the coupling terms are calculated for more general potentials with time and space-dependent Rabi frequencies �(x, t) and wave vectors k(x, t). For this purpose, the momentum-space representation of the Schrödinger equation is more appropriate and can be written for the Fourier transform of the atomic wave function g(p, t) where Expressing the wave function in momentum space gives Discretising p → (j + δ) k and p ′ → (l + γ ) k , one finds where l and γ span the same indices ensembles as j and δ . The new equations to solve read (39) V (p) * g(p, t) ≈ 1 2π l,γ F((j + δ) k, (l + γ ) k, t)g l+γ (t), Scientific Reports | (2020) 10:22120 | https://doi.org/10.1038/s41598-020-78859-1 www.nature.com/scientificreports/ which yields the necessary momentum couplings for an arbitrary potential V(x, t). In the worst case, the sum in Eq. (40) runs over N eq nonzero entries ( N l N γ = N j N δ = N eq ) which leads to a time complexity of O(N 2 eq ) . This, however, is an extreme example that contrasts with commonly operated precision interferometric experiments since it would correspond to white light with speckle noise. Realistic scenarios rather involve time-dependent potentials with a smaller number of momentum couplings, i.e. N eq ≫ #coupling terms 2 as would be the case in Fig. 10c. To evaluate the momentum couplings, it is necessary to calculate the integral F(p, p ′ , t) at each time step using the FFT, which leads to a final time complexity class for solving the ODE of O(N eq log N eq ).
The next important generalisation aims to include the effect of the two-body collisions analysed in the meanfield approximation, i.e. H int = g 1D |ψ(x, t)| 2 . In this case, the equation describing the dynamics of the system and the couplings can be written as where ν and o are running indices over the same values as l and γ . One ends up with N eq differential equations where each has more than N 2 eq coupling terms, and finds a time complexity class of O(N 3 eq ) . This shows the growth in numerical operations of the ODE treatment as reflected by the number of couplings in Fig. 10d.
We analyse now the time complexity class for the PDE approach, using the split-operator method 24 . Based on the application of the FFT, it is known that the complexity class of this method is scaling as O(N grid log N grid ) , (40) i ġ j+δ (t) ≈ ((j + δ) k) 2 2m g j+δ (t) + 1 2π l,γ F((j + δ) k, (l + γ ) k, t)g l+γ (t), (41) i ġ j+δ (t) = ((j + δ) 2 ω r + �)g j+δ (t) + � 2 (g j+δ+2 (t) + g j+δ−2 (t)) www.nature.com/scientificreports/ where N grid is the number of grid points in the position or momentum representations. Since the discretisation of the problem for the ODE and PDE (Schrödinger equation) approaches is roughly the same ( N eq ≈ N grid ), a direct comparison between the two treatments is possible. The time complexity analysis is summarised in Table 1. It shows that the standard ODE approach is only better suited in the case of ideal light plane waves. In every realistic case where the light field is allowed to be spatially inhomogeneous, the amount of couplings increases and it is preferable to employ the PDE approach with a scaling of O(N grid log N grid ) , independently of any further complexity to be modelled.

Conclusion
In this paper, we have shown that the position-space representation of light-pulse beam splitters is quite powerful for tackling realistic beam profiles in interaction with cold atom ensembles. It was successfully applied across several relevant regimes, geometries and applications. We showed its particular fitness in treating metrologicallyrelevant investigations based on atomic sensors. Its high numerical precision and scalability makes it a flexible tool of choice to design or interpret atom interferometric measurements without having to change the theoretical framework for every beam geometry, dimensionality, pulse length or atomic ensemble property. We anticipate the possibility of accurately implementing this approach to analyse important systematic effects in the field of precision light-pulse matter-wave interferometry such as the ones related to wavefront aberrations, large momentum transfer and inhomogeneity and fluctuations of the Rabi pulses. Finally, we would like to highlight the possibility to generalise this method to Raman or 1-photon transitions if we account for the internal state degree of freedom change during the diffraction.
Received: 13 July 2020; Accepted: 30 November 2020 Table 1. Comparison of the different time complexity classes of the commonly-used ODE treatment with the position-space approach developed in this work (PDE-based). Including more and more realistic features of the atom-light system leads to an ODE time complexity unfavourably scaling. The PDE formulation, however, routinely scales with O(N grid log N grid ).