Lattice-depth measurement using multi-pulse atom diffraction in and beyond the weakly diffracting limit

Precise knowledge of optical lattice depths is important for a number of areas of atomic physics, most notably in quantum simulation, atom interferometry and for the accurate determination of transition matrix elements. In such experiments, lattice depths are often measured by exposing an ultracold atomic gas to a series of off-resonant laser-standing-wave pulses, and fitting theoretical predictions for the fraction of atoms found in each of the allowed momentum states by time of flight measurement, after some number of pulses. We present a full analytic model for the time evolution of the atomic populations of the lowest momentum-states, which is sufficient for a"weak"lattice, as well as numerical simulations incorporating higher momentum states for both relatively strong and weak lattices. Finally, we consider the situation where the initial gas is explicitly assumed to be at a finite temperature.


I. INTRODUCTION
Precision measurement of optical lattice [1] depths is important for a broad range of fields in atomic and molecular physics [2,3], most notably in atom interferometry [4,5], many body quantum physics [6,7], accurate determination of transition matrix elements [8][9][10][11][12], and, by extension, ultraprecise atomic clocks [13,14]. Lattice depth measurement schemes include methods based on parametric heating [15], Rabi oscillations [16], and sudden lattice phase shifts [17]. The most commonly used scheme is Kapitza-Dirac scattering [18], where an ultracold atomic gas is exposed to a pulsed laser standing wave and theoretical predictions for the fraction of atoms found in each of the allowed momentum states are fitted to time of flight measurements [7,[19][20][21]. However, when determining the matrix elements of weak atomic transitions, the lattice depths involved are correspondingly small (V 0.01E R for any atom, here V is the lattice depth and E R is the laser recoil energy), such that signal-to-noise considerations become an issue [22].
Recently, the work of Herold et al. [23] and Kao et al. [24] has suggested that this complication can be mitigated by using multiple laser standing wave pulses, alternating each with a free evolution, such that each alternating stage has a duration equal to half the Talbot time [25][26][27]. With each pulse, population in the first diffraction order is coherently increased, improving contrast relative to the zeroth order. 1 The modeling approach taken in [23,24] is valid for a weak lattice which is pulsed a small number of times, corresponding to the "weakly-diffracting limit". Following description of our model system and its general time evolution in section II, in section III we present a full analytic model for the time evolution of the atomic populations of the zeroth and first diffraction orders; this is sufficient for a "weak" lattice.
In section IV we present numerical simulations incorporating higher momentum states at both large and small lattice depths V ("small" is taken to mean when V is less than a tenth of the recoil energy E R ), which we compare for typical experimental values. We also explore the role of finite-temperature effects in such experiments (section V), and present our conclusions in section VI.

A. Alternating Hamiltonian evolutions
We consider an atomic Bose-Einstein condensate (BEC) with interatomic interactions neglected. 2 This can be achieved experimentally by exploiting an appropriate Feshbach resonance [28][29][30][31], or by allowing the cloud to expand adiabatically [32]. Working in this regime means that we need only consider the single-particle dynamics of each atom. The optical lattice laser is far off resonance such that we consider the atomic center of mass motion only [33], and we consider the atoms to be periodically perturbed by a 1d optical lattice, alternated with a free evolution [34]. The atomic center of mass dynamics are then alternatingly governed by the following Hamiltonians:Ĥ wherep is the 1d momentum operator in the x direction (see Fig 1),x is the associated position operator, M is the atomic mass, and V the lattice depth 3

(dimensions of energy) of a
Mirror Laser

Atomic gas
Laser standing wave shows the modulation of the lattice depth in time, where V is the lattice depth (dimensions of energy) when the standing wave pulse is on, and T 1/2 is the Talbot time as defined in Eq. (2). For simplicity, the laser standing wave has been oriented orthogonally to the gravitational direction, however we note that this is equivalent to a vertically oriented system in which a phase-shifter element is used to introduce a time dependent phase on the standing wave, which is tuned to cancel out gravity [34,37].
lattice with wavenumber K (K = 2K L , where K L is the laser wavenumber) [35,36]. As stated in the introduction, Herold et al. [23] and Kao et al. [24] proposed that when measuring very small lattice depths (V∼0.01E R , here E R = 2 K 2 /8M), the signal can be optimized by both the lattice pulse and free evolution having a duration equal to the half Talbot time [38], This is half the full Talbot time, which is the elapsed time for which the free evolution operator [generated by Eq. (1b)] collapses to the identity when applied to a momentum state that is an integer multiple of K. 4

B. Time evolution
The time-periodicity of the system admits a Floquet treatment [44]; the time evolution of an initial state |ψ(t = 0) for N successive lattice-pulse sequences is given by repeated applications of the system Floquet operatorF to the initial state, i.e., |ψ(t = N) =F N |ψ(t = 0) .
We determine the relevantF, governing a lattice pulse of duration T 1/2 [Eq. (2)], followed by a free evolution of the same duration, straightforwardly from the time evolution operators generated by Eqs. (1a) and (1b). The spatial periodicity of the laser standing wave also enables us to invoke Bloch theory [47]. Recasting the momentum operatorp such that: with k ∈ Z and β ∈ [−1/2, 1/2) [48], we elucidate that the total dimensionless momentum ( K) −1 p associated with a single plane wave is the sum of k, the discrete part, and β as the continuous part or quasimomentum, which is a conserved quantity. Hence, only momentum states separated by integer multiples of K are coupled [34,49]. Within a single quasimomentum subspace, the system Floquet operator can therefore be written: where V eff = V M/ 2 K 2 is the dimensionless lattice depth,θ = Kx and the rescaled half Talbot time is equal to 2π. 5 Using Eq. (4) to calculate |ψ(t = N) = j c j (N)|k = j , the population in each discrete momentum state |k = j after N pulses is given by the absolute square of the individual coefficients P j (N) = |c j (N)| 2 . In this paper we employ both the well-known splitstep Fourier approach [34,50], and matrix diagonalization in a truncated basis [23,51] to determine |ψ(t = N) beyond the weakly-diffracting limit, as well as an analytic approach in the weakly-diffracting case.

III. ANALYTIC RESULTS IN A TWO-STATE BASIS
For an initially zero-temperature gas (β = 0) subjected to a small number of pulses from a shallow lattice, a useful approximation is to assume that no population is diffracted into momentum states with |p| > K, the so-called "weaklydiffracting limit". Mathematically, this regime corresponds to the time evolution of an initial state |ψ(t = 0) = |k = 0 in a space spanned only by the |k = −1 , |k = 0 and |k = 1 states of the β = 0 quasimomentum subspace.
The symmetry of the lattice and free evolution Hamiltonians about |k = 0 guarantees that, for our chosen initial state, the population diffracted into the |k = 1 state is identical to that diffracted into the |k = −1 state. We therefore express the system Hamiltonians (1a), (1b) as matrices in the truncated momentum basis: yielding the following 3×3 matrix representation of the lattice Hamiltonian: There is no coupling between the |0 state and the antisymmetric |− state. Hence, for an initially zero-temperature gas, there is no population transfer into the |− state for all time. The relevant basis is therefore two-dimensional, with basis states |0 2 ≡ 0 1 and |+ 2 ≡ 1 0 . We use these to represent Eq. (1a) as the 2 × 2 matrix: We recognize Eq. (7) as a Rabi matrix, the eigenvalues and normalized eigenvectors of which are well known [52]. We use these to calculate the populations after N pulses of the |0 and |+ states [P 0 (N, V eff ) and P + (N, V eff ), respectively]: A = 8V 2 eff sin 2 π 1 + 8V 2 eff /2 as explicitly derived in Appendix A. From Eqs. (8a) and (8b), we see that in the weaklydiffracting limit P 0 and P + oscillate sinusoidally with the number of pulses N, and are entirely characterized by an amplitude A and a "frequency" φ, both of which depend solely on the dimensionless lattice depth V eff . We note the similarity to the result reported in [21] for single pulse diffraction. We display the variation of A and of φ versus V eff in φ initially increases approximately linearly with V eff , meaning that over a sufficiently small range of lattice depths, we should expect to see an approximate universality in the population dynamics when the time axis is scaled by V eff (we explore this scaling in Section IV). In the limit where V eff → 0, it follows that φ = 4 √ 2V eff (see Appendix B), depicted by the solid straight line plotted in Fig. 2(a). Substituting this result into Eq. (8b) and expanding the corresponding Taylor series to leading order, we recover the familiar quadratic dependence of Herold et al. [23,24] (see Appendix B 3): The validity of this result is subject to Nφ(V eff )/2 1. Increasing V eff beyond this regime, A, which decreases steadily in the range of linearity of φ, first reaches a node at Physically, these values of V eff correspond to there being no pulse-to-pulse population transfer out of the |k = 0 state, at least in the weakly-diffracting limit. As shown in Appendix B, φ = π at those values of V eff where A has a node, visualised by the intersection of the vertical and horizontal dashed lines in Fig 2(a).
with an overall oscillatory behavior of ever-decreasing amplitude around this value, while A takes on the form of a sinusoidal oscillation: A = sin 2 ( √ 2πV eff ).
"Atan2" numerical routine in, e.g., Python. This ensures that the sign of the argument is taken into account, which avoids singularities in the frequency.

A. Numerical simulations for a large momentum basis
Having obtained analytic results for the time-evolved populations in the weakly-diffracting limit, we test their domain of validity by using standard numerical techniques to compute the full momentum distribution of the system, and sampling the population in the |k = 0 state, P 0 . We follow the same approach as [44,50] and work within the momentum basis. The action of the Floquet operator (4) on the total state of the system, |ψ , is calculated by a split-step Fourier method, on a basis of 2048 momentum states, which is exhaustive for any practical purpose.
In Fig. 3 we compare the analytic results of Eqs. (8a,8b,8c,8d) to this exact numerical calculation for fixed values of the effective lattice depth V eff . From Fig. 3(a) we see that the sinusoidal character of the analytic result for P 0 is revealed for higher values of V eff , as well as a similar oscillatory behavior in the exact numerics. Naively, we may say that increasing V eff gives rise to a greater deviation of the exact numerics from the analytics. This is true when comparing over a fixed number of pulses, however we can use our argument that there is an approximate universality in V eff and the number of pulses (see section III) to clarify this statement by means of the universal curve displayed in Fig. 3(b). This clearly shows that the universality holds approximately for the exact numerics also, and that the analytics cease to agree with the exact numerics at approximately the same point on the universal curve, regardless of the value of V eff in the chosen range. Hence, more completely, the analytics are sufficient to understand the system provided the product of the number of pulses and effective lattice depth is sufficiently small. We note specifically that there is a frequency drift which increases along the curve, and a marked reduction in amplitude of the exact numerics as compared to the analytics at its first revival. Both features appear due to leakage of population into momentum states with |p| > K, and inform our discussion of the range of validity of the weakly-diffracting limit taken in previous work. Indeed, the quadratic result of Herold et al. [Eq. (9), shown as dashed lines in Fig. 3] deviates from the exact numerics at a significantly smaller value of NV eff than our exact analytic result for two diffraction orders.
In [23,24], the regime in which the weakly-diffracting limit is satisfied (recast in our system of variables) is given by NV eff 1/4. Though this inequality places an upper bound on the allowed value of NV eff , it is reasonable to ask at what point is NV eff "much smaller" than 1/4? By inspection of Fig. 3(b), we can see that at NV eff = 1/4, there is still excellent agreement between our analytics and exact numerics. We calculate the RMS difference between our analytics and full numerics [53] at this point over the range of chosen lattice depths (defined as RMS=  9)]. Each set of markers corresponds to a fixed value of the effective lattice depth ranging from the slowest-oscillating curve at V eff = 0.01 to the fastest oscillating one at V eff = 0.11 in steps of 0.02. (b): Reproduction of (a), with the number of pulses axis scaled by the dimensionless lattice depth V eff to reveal an approximate universal curve both in the analytics and the numerical simulations. The data have been extended to span the full range of the horizontal axis. The universal curve reveals a drop in the amplitude of P 0 as calculated by the full numerics at the first revival, which is not reproduced by the analytics, but is reproduced in the truncated momentum basis. In (b), the oscillation frequency of the numerical curve increases compared to that of the analytic result as the number of pulses or the lattice depth is increased. After three half-oscillations on the universal curve, the truncated basis result begins to deviate appreciably from the full numerics. corresponding quadratic result deviates at the 42% level. 7 The point at which leakage into higher momentum-states first becomes appreciable is NV eff ∼ 1/2, with an RMS of 0.0043 (deviation at the 0.4% level). Though this is clearly suffi- Each false-color plot shows the time evolved population in the first 13 momentum states (|k| < 6), to be read on the colorbar to the right. A cutoff population value of P cutoff = 10 −11 has been applied to each population distribution to accommodate the log scale. This illustrates that for this choice of parameters, the amount of population diffracted into momentum states with |p| > 3 K is negligible. Row 2 [(b), (d), (f)] shows firstly, slices through the momentum distribution corresponding to the population in the k = 0 state, P 0 , (red circles) and the |p| = K states, P ±1 , (blue squares), to which our two-state analytic model is compared (red and blue solid lines respectively). To clarify the drop in amplitude in the first revival of P 0 , the green triangles have been added, which correspond to 1 − P ±2 and almost intersect the red circles corresponding to P 0 , indicating that the overwhelming majority of the population which has left P 0 at this point, has in fact been diffracted into the |k| = 2 states. At the second revival, the two sets of points are further apart. Population leakage into the |p| = 3 K states, corresponding to the magenta diamonds, which represent 1 − P ±3 , explains this effect. Solid lines have been added as a guide to the eye. Each column corresponds to a fixed value of ciently small to still be considered within the range of validity of the weakly-diffracting limit, beyond NV eff ∼ 1/2, where the RMS becomes larger, we must incorporate higher momentumstates. This motivates the question of how many momentum states are necessary to include for such a model to be useful for a reasonable choice of experimental parameters. Figures 4 (a,c,e) show a selection of momentum distributions for a range of values of V eff as calculated by the full numerics, showing momentum states up to |p| ≤ 6 K, with Figs. 4 (b,d,f) showing corresponding slices through the momentum distributions. The log scale makes clear that there is very little population leakage into momentum states with |p| > 3 K for the chosen values. Instead we see that there are pronounced oscillations in population between the |p| = 0 and |p| = K states, which are modulated by population leakage into the |p| = 2 K states, and to a lesser extent the |p| = 3 K states. By inspection of the lattice Hamiltonian in the momen-tum basis, this can be explained by the decrease in magnitude of the off-diagonal coupling terms with state number. In fact, the decrease in amplitude at the first revival in Fig. 3(b) is almost entirely due to population leakage into the |p| = 2 K states, suggesting that a model incorporating only n = 5 momentum states should be sufficient to capture the dynamics, up to at least V eff = 1.1.

B. Small momentum bases of dimension > 2
To incorporate higher momentum-states we numerically diagonalize Eqs. (1a) and (1b), in a truncated basis of n momentum states, and propagate the time-evolution using the procedure described in Appendix C. Our analysis in the previous section suggests that simulations using a basis of n = 5 momentum states ought to be sufficient for practical purposes.
Corresponding results are shown by the hollow markers in Fig. 3(b). The five state model is an order of magnitude more accurate than the analytics at NV eff = 1/4 and NV eff = 1/2, with RMS differences with respect to the full numerics of 0.00018, and 0.00011 respectively. As expected, the decrease in amplitude at the second revival on the universal curve is reproduced by this approach, but is clearly also valid over a larger range, up to the fourth turning point (NV eff ∼ 1.6, RMS deviation 0.0022), beyond which the model begins to overestimate and then underestimate the exact numerical result.
This difference appears as a result of the basis truncation, as population leakage into states with |p| ≥ 5 K is explicitly not possible in this model, though it should be noted that this effect would only be relevant to experiments performed using a very large effective lattice depth. An attractive feature of the five state model is that it can in principle be solved analytically for the time-evolution of the populations, which could be fit to experimental data to extract more accurate lattice depths.

V. FINITE-TEMPERATURE RESPONSE
The results presented in the previous sections are valid for a gas which is assumed to be initially at zero temperature; in practice this regime is never fully achieved, even for a BEC. To find the response of P 0 versus the number of pulses for a finite-temperature gas, we calculate the time evolution of P 0 for an ensemble of initial momentum states |ψ(t = 0) = |( K) −1 p = k + β according to Eq. (4), where the initial momentum is defined in a Bloch framework with k and β as free parameters. For a sufficiently cold gas [temperature T w ( 2 K 2 /64k B )K] 8 we need only consider initial states with k = 0 in order to capture the essential features. In this regime we choose a fixed value of the lattice depth and scan across the full range of the quasimomentum β as the only free parameter, to find the momentum dependence in the first Brillouin zone [47] displayed in Fig. 5. Figure 5 clearly shows the central resonance at β = 0, where our zero-temperature analysis is applicable. Increasing the quasimomentum to |β| = 0.0625, we see that the oscillation in P 0 has an amplitude of less than 50% of that at β = 0, and a substantially different frequency. Hence, the width of the central resonance is relatively narrow compared to the full width of the Brillouin zone. For an initial momentum distribution of appreciable width we must consider the surrounding structure when calculating the population dynamics, as the zerotemperature behavior will be washed out over time, or even be unresolvable altogether if the temperature is sufficiently high.
Note that for broader initial momentum distributions the dynamics will include the secondary resonances at |β| = 0.5, which have a periodicity of the form P 0 (N) = cos 2 (πV eff N), such that P 0 varies between 0 and 1 for all V eff .
Having characterized the first Brillouin zone, we calculate the full finite-temperature response of P 0 by performing Gaus- 8 This rule of thumb is chosen such that the initial width of the momentum distribution is at most one quarter that of the first Brillouin zone. sian weighting in momentum space according to a rescaled Maxwell-Boltzmann distribution: where the dimensionful temperature is given by T w = 2 K 2 w 2 /Mk B [44], and k B is Boltzmann's constant. Figure 6 shows the variation of P 0 with the number of pulses, including both the strong and weak lattice regimes, and three different values of the initial momentum distribution width w. In overview: in regimes where we have a weak lattice and low temperature the analytic formula is adhered to almost perfectly; in regimes where we have a weak lattice and a higher temperature we begin to see noticeable deviations, which occur for a smaller number of pulses as the temperature is increased; in regimes where we have a strong lattice and low temperature, although the analytic formula is not strongly adhered to as the lattice depth increases, the oscillation frequency appears to be reasonably robust as V eff increases and the amplitude of oscillation consequently decreases; finally in the regime of strong lattice and higher temperature, the analytic formula is again only adhered to for relatively short times, with that time being dependent on the temperature.

VI. CONCLUSIONS
We have a zero-temperature analytic formula which yields significant insight assuming that we are working in the weakly-diffracting limit. We have shown that at zero temperature, very small basis sizes are sufficient to capture the essential features of the population dynamics outside the weaklydiffracting limit. We have explored the effects of finite temperature initial distributions, and elucidated regimes from which the lattice depths can be determined from the observed dynamics in the lowest diffraction order.

Floquet operator in two-state basis
We may calculate the time evolution of the |0 and |+ state populations by first diagonalizing Eq. (7) (reproduced here for convenience) using the well known eigenvalues and normalized eigenvectors of a Rabi matrix, E ± = (1 ± 1 + 8V 2 eff )/4, and respectively, where α = arctan(2 √ 2V eff ). H 2×2 Latt can then be written: such that R is the matrix of normalized eigenvectors. This leads directly to the part of the Floquet operator governing the lattice evolution: Expressing F Free in the truncated momentum basis, |0 2 ≡ 0 1 ; |+ 2 ≡ 1 0 , we can represent the total Floquet operator in matrix form thus: 2. Floquet evolution for a general two-level system Any time-evolution operator associated with a two-level system can be expressed as a 2 × 2 unitary matrix, and all unitary matrices are diagonalizable, hence we may represent such a time-evolution operator thus: Here S is a matrix composed of the normalized eigenvectors of U: and λ ± are the corresponding eigenvalues of U, which have unit magnitude and so can be expressed as: where θ + and θ − are phase angles to be determined. The matrix which produces N successive evolutions can therefore be written: Suppose that the initial state of the system can be represented by |0 2 ≡ 0 1 , and the excited state by |+ 2 ≡ 1 0 , the probability of the system occupying the |0 state after N evolutions can be written: which is the absolute square of the top-left matrix element of Eq. (A9). The corresponding probability of the system being in the |+ state is simply P + (N) = 1 − P 0 (N). Since S is a unitary matrix, v + 0 and v − 0 must satisfy |v + 0 | 2 + |v − 0 | 2 = 1, using this identity and inserting Eq. (A8), P 0 (N) and P + (N) can be written: By finding v ± 0 and θ ± for our specific Floquet operator (A5), we explicitly determine Eq. (A11a) and (A11b), in terms of the number of pulses N and the effective potential depth V eff , this is the origin of Eq. (8a) and (8b).

Quadratic approximant to Equation (8b)
Equation (8b) can be rewritten in terms of the first few orders of a Taylor expansion: with x ≡ Nφ/2, in a regime where x 1. Further, assuming that V eff is near zero, we may replace φ and A with our leading order approximations of Eqs. (B7,B9), with x ≈ 2 √ 2NV eff .
where α, γ ∈ Z. Equation (C3) can then be expressed in matrix form, and numerically diagonalized in order to find the time evolution of an initial momentum eigenstate.
By expressing Eq. C3 in matrix form thus: we can construct the matrix P n×n diagonalizing H n×n latt , such that H n×n latt,diag = (P † ) n×n H n×n latt P n×n . We are led to the expression: |ψ(t = N) n×1 = [H n×n free P n×n H n×n latt,diag (P n×n ) † ] N |K = α n×1 , (C5) for |ψ(t = N) n×1 , the time evolution due to N pulse sequences of an initial eigenstate |K = α n×1 , where α ∈ [−(n − 1)/2, (n − 1)/2]. The n × 1 superscript denotes that the ket should be understood as an n-dimensional column vector.