Characterization of Errors in Interferometry with Entangled Atoms

Recent progress in generating entangled spin states of neutral atoms provides opportunities to advance quantum sensing technology. In particular, entanglement can enhance the performance of accelerometers and gravimeters based on light-pulse atom interferometry. We study the effects of error sources that may limit the sensitivity of such devices, including errors in the preparation of the initial entangled state, imperfections in the laser pulses, momentum spread of the initial atomic wave packet, measurement errors, and atom loss. We determine that, for each of these errors, the expectation value of the parity operator $\Pi$ has the general form, $\overline{\langle \Pi \rangle} = \Pi_0 \cos( N \phi )$, where $\phi$ is the interferometer phase and $N$ is the number of atoms prepared in the maximally entangled Greenberger--Horne--Zeilinger state. Correspondingly, the minimum phase uncertainty has the general form, $\Delta\phi = (\Pi_0 N)^{-1}$. Each error manifests itself through a reduction of the amplitude of the parity oscillations, $\Pi_0$, below the ideal value of $\Pi_0 = 1$. For each of the errors, we derive an analytic result that expresses the dependence of $\Pi_0$ on error parameter(s) and $N$, and also obtain a simplified approximate expression valid when the error is small. Based on the performed analysis, entanglement-enhanced atom interferometry appears to be feasible with existing experimental capabilities.


I. INTRODUCTION
Neutral atoms are used in some of the most precise, stateof-the-art quantum sensors, including atomic clocks [1][2][3][4], optical atomic magnetometers [5][6][7], and atom interferometers (AIs) [8][9][10][11]. Looking forward, an ambitious goal is to harness the power of quantum entanglement to decrease the phase uncertainty in these atomic sensors [12], first, beyond the standard quantum limit (SQL), N −1/2 , which arises in measurements with independent atoms, and, ultimately, as close as possible to the Heisenberg limit (HL), N −1 , where N is the number of atoms used in the measurement.
There are two main approaches to generating metrologically useful entanglement in atomic systems. One approach utilizes spin squeezing [13][14][15] generated in an ensemble of ultracold neutral atoms [16][17][18][19]. It has been shown that spin squeezing makes it possible to surpass the SQL in atomic clocks [18,20] and atomic magnetometers [21,22]. In these experiments, spin-squeezed states were generated in large ensembles of atoms (10 5 -10 6 ), but only a small fraction of the atoms (∼ 0.1%) have been actually entangled [18]. Also, strong spin squeezing is generated via an optical-cavity-based measurement, and releasing the atoms from the cavity into free space (which is typically required in order to use them in an AI) results in a fast degradation of squeezing [23].
An alternative approach is to produce maximally entangled atomic spin states such as the Greenberger-Horne-Zeilinger (GHZ) state [24], which are known to achieve the HL under ideal conditions [25]. The most promising method for generating high-fidelity entangled states of atomic spins is by using Rydberg-mediated interactions in arrays of ultracold, optically trapped neutral atoms [26]. In particular, this method was used in a number of experiments [27][28][29][30][31] to generate * cnbrif@sandia.gov † bruzic@sandia.gov ‡ biedermann@ou.edu maximally entangled two-spin states (i.e., Bell states [32]), with fidelity as high as 0.97 [29]. This method, combined with optimal control [33], also produced GHZ spin states in arrays of 4 to 20 atoms with fidelities of 0.85 to 0.54, respectively [34]. Additional proposed advances such as the use of rapid adiabatic Rydberg dressing [35] and in situ adaptive optimal control [36], along with various technical improvements [34], are likely to further increase the number of entangled atoms and enhance fidelities in these systems. Although maintaining a high fidelity as the number of entangled atoms increases is an outstanding challenge, the possibility of producing highly entangled spin states of ultracold neutral atoms is a promising avenue for quantum sensing. One of the most prominent sensing techniques for neutral atoms is light-pulse atom interferometry [37][38][39][40]. Inertially sensitive AIs have many important practical applications [8,9,11], including accelerometers and gyroscopes [41][42][43][44][45], and gravimeters [46][47][48][49]. Promising advances also demonstrate gravity measurements using spatially separated atomic wave packets suspended in a lattice [50][51][52] and most recently with ultralong interrogation times [53].
Feeding atoms prepared in a highly entangled spin state into an AI [54][55][56] opens up an exciting path towards inertially sensitive measurements with a phase uncertainty beyond the SQL or even close to the HL. The unparalleled control achieved in ultracold neutral atom experiments [26,34] makes them an ideal platform to both reach the fundamental limits of AI-based sensors and understand key errors that affect their performance.
In this paper, we provide a detailed analysis of a light-pulse AI utilizing ultracold atoms that have been prepared in the GHZ spin state. We devise a protocol for the operation of the entanglement-enhanced AI and discuss the relevant physical parameters that affect its performance. In the ideal case of an AI without any errors, the maximally-entangled GHZ state achieves a phase uncertainty at the HL, similarly to the analogous result in spectroscopy [25]. Of course, in any physical implementation of the AI, there will be many sources of noise that reduce its sensitivity from the ideal case. Our main goal is to characterize the non-negligible sources of noise and develop mitigation strategies for them if necessary. First, we consider the effect of imperfections in the initial many-atom GHZ state. Next, we investigate how various types of noise affect the operation of the AI after the GHZ state has been prepared with a given fidelity. We characterize the relevant noise sources in terms of error parameters and analyze how the phase uncertainty scales as a function of the number of entangled atoms N in the initial GHZ state. We use the error parameters and the scaling laws we derive to learn about the fundamental limitations of this type of sensor.
The non-negligible errors we consider come from noise in the initial state preparation, laser intensity fluctuations, laser phase noise, the initial momentum spread of the atoms, measurement errors, and atom loss. The initial momentum spread leads to a detuning error in the light pulses of the AI, and we find that this error is dominant for AI operation. This error grows linearly with the atom's vibrational energy in the trap and decreases with the square of the effective Rabi frequency for the two-photon Raman transition. Consequently, the effect of the initial momentum spread can be reduced by cooling the atom to the ground state of the trap [57], by decreasing the trap frequency through adiabatic lowering of the trap depth [58], and by increasing the Rabi frequency through the use of a high-intensity laser with tight focusing. We forecast that a realistic implementation of an entanglement-enhanced AI can achieve inertial sensing near the HL for ∼ 100 atoms.
We also develop a detection scheme that allows for a measurement of the interferometric phase φ, with a number of measurements that scales linearly with N , in contrast to the exponential scaling of the Hilbert space dimension for a system of N entangled atoms. This is possible because a measurement of the parity of the final N -atom spin state, which provides a sufficient amount of information to determine φ, can be performed via a state-selective detection [59][60][61] with only N + 1 possible outcomes, instead of needing to distinguish between the 2 N possible final states. We estimate the effect of measurement error for this detection scheme.
Finally, we quantify the effect of atom loss. By using the capability of detecting whether an atom was lost, we can postselect only lossless outcomes, which will result in a reduced data-acquisition rate. If the total number of experiments (AI cycles) is fixed, a decrease in the number of lossless experiments can be interpreted as an effective deterioration of phase uncertainty per one experiment.

II. BACKGROUND
Precision AI experiments [38][39][40] typically employ the clock transition between the ground-state hyperfine levels with m F = 0 (this transition is magnetic field insensitive). In particular, for 133 Cs atoms, the electronic ground state is 6S 1/2 , and the clock transition is between the hyperfine levels |g = |F = 3, m F = 0 and |e = |F = 4, m F = 0 . The standard experimental approach [38,39] is based on driving the two-photon stimulated Raman transition between these levels, via an intermediate level |i , as schematically shown in where m is the atom's mass, ω g , ω e , and ω i are the energies of the three atomic levels, d is the atom's dipole operator, and E is the electric field of two Raman beams: Here, x is the atom's position, and each Raman field is characterized by its amplitude E j , wave vector k j , frequency ω j and phase φ j (j = 1, 2). A solution to the Schrödinger equation governed by the Hamiltonian of Eq. (1) can be generally expressed as a superposition: where the basis state |α, p = |α ⊗ |p corresponds to an atom in its internal state |α (α = {g, e, i}) and in its momentum eigenstate |p , whose position-space representation is ψ p (x) = e ip·x/ . The time dependence in Eq. (4) captures the field-free evolution. Due to momentum conservation, for a given p, the Hamiltonian of Eq. (1) only couples the states |g, p , |i, p + k 1 , and |e, p + K , where K = k 1 − k 2 . The momentum kick K is maximized by using counter-propagating Raman beams, for which K ≈ 2k 1 . Adiabatic elimination of the intermediate level |i results in a model for an effective two-level system, whose Hamiltonian couples the states |g, p and |e, p + K . In the fieldfree evolving basis {|g, p (t), |e, p + K (t)}, this effective Hamiltonian is [40] where , Here, ω 12 = ω 1 − ω 2 is the frequency difference between the two Raman fields, ω eg = ω e − ω g is the frequency difference between the hyperfine levels, and ∆ = ω 1 − (ω i − ω g ) is the detuning from the optical resonance, as shown in Fig. 1. The diagonal elements of the Hamiltonian in Eq. (5) are the ac Stark shifts of levels |g and |e , and the relative ac Stark shift of the two levels is δ AC ≡ Ω AC e − Ω AC g . The effective detuning from the Raman resonance, δ 12 in Eq. (8), includes the two-photon detuning, δ = ω 12 − ω eg , the detuning due to the Doppler shift, −p · K/m, and the detuning due to the atom's recoil energy, − |K| 2 /2m. The phase φ 12 is chosen to make the effective Rabi frequency Ω eff a positive real number.
The Schrödinger equation governed by the Hamiltonian of Eq. (5) can be solved analytically to obtain the evolution operator for a single atom [40]. If the Raman fields are turned on at time t and act for duration τ , the resulting evolution operator, denoted as U t (τ ), is represented in the basis {|g, p (t), |e, p + K (t)} as where we omitted a global phase factor e −i(Ω AC g +Ω AC e )τ /2 and defined If the atom moves at a constant velocity, then δ 12 is constant, and φ t = δ 12 t. However, if the atom accelerates, δ 12 is time dependent, and Eq. (11) should be used. Note that Eq. (10) neglects the time dependence of δ 12 during the pulse, assuming that the pulse duration τ is very short compared to periods of free evolution. Since the time dependence for the fieldfree evolution is included in Eq. (4), the evolution operator for a period of field-free evolution is represented in the basis {|g, p (t), |e, p + K (t)} by the two-by-two identity matrix. We will use the general form (10) in Sec. VII where we explicitly take into account the effect of the initial momentum spread of the atoms, which leads to a detuning uncertainty via the Doppler shift term in Eq. (8). However, in Secs. III-VI we assume that the Raman detuning is much smaller than the effective Rabi frequency, |δ 12 | Ω eff . This assumption is satisfied if three conditions are met: (1) short, intense driving fields are used to make Ω eff large, (2) the momentum distribution of the atoms is sufficiently narrow to make |δ 12 | small for all relevant momentum components, and (3) laser frequency chirping is used to compensate for the evolving Doppler shift due to the acceleration of the atom (see [39,40] for more details). Under this assumption, we neglect terms on the order of |δ 12 − δ AC |/Ω eff and |δ 12 |τ in Eq. (10) to obtain: where A ≡ Ω eff τ is the pulse area and, since τ enters only via A, we changed the notation from U t (τ ) to U t (A). We will use this approximation for the pulse evolution operator throughout this paper, except for Sec. VII and Appendix B.
Since atoms interact independently with light fields, the evolution operator for N atoms is where U (k) t (A) is the evolution operator for the kth atom. A typical AI operation includes three pulses: π/2-π-π/2, with two periods of field-free evolution, each of duration T , between the pulses. The evolution operator for one entire cycle of the AI operation is where t 1 is the starting time of the first π/2 pulse and, correspondingly, t 2 = t 1 + τ /2 + T and t 3 = t 1 + 3τ /2 + 2T , assuming that the π pulse has a duration τ and each of the π/2 pulses has a duration τ /2. In what follows, without loss of generality, we set t 1 = 0. We also assume that the duration of each pulse (typically, τ ∼ 1 µs) is negligible compared to the duration of the free evolution (typically, T ∼ 1 ms). Under this assumption, we neglect terms on the order of |δ 12 |τ [which is consistent with the approximation we made in deriving Eq. (12)] and obtain: Immediately after the final pulse, a state-dependent detection of the atoms is used to measure the interference. Specifically, we assume that one measures the expectation value of the parity operator, where σ (k) z = |g k k g| − |e k k e| is the Pauli z matrix for the kth atom's spin. Since the parity operator (or any operator corresponding to a measurement of the populations of the atomic levels and hence diagonal in the basis {|g k , |e k }) is invariant under the field-free evolution, its expectation value at the final time t = 2T is given by where 1 1 p = N k=1 d 3 p k |p k p k | is the identity operator for the motional degrees of freedom and |Ψ(0) is the initial state of the system of N atoms.

III. THE IDEAL CASE
As a reference point, we first consider the case of an AI without errors. In this ideal case, we ignore all physical errors in the Raman pulses, which means that both A and φ 12 in the evolution operator of Eq. (12) are precisely known.
The evolution operator U t (A) of Eq. (12) represents the rotation by an angle A about the axisx cos(φ t + φ 12 ) − y sin(φ t + φ 12 ). The parity expectation value is independent of the specific choice of φ 12 , so, for the sake of simplicity, we use φ 12 = −π/2 throughout the rest of this paper. With this choice, the evolution operator of Eq. (12) becomes Specifically, the evolution operators for the three AI pulses are 2T 1 .
Note that in Eq. (19) we used the fact that φ t=0 = 0.

A. Ideal interferometer with independent atoms
For an AI with N independent atoms, each atom interferes with itself [44]. Let us first consider one atom which is initially in the state  2. A scheme of an AI with the π/2-π-π/2 sequence of pulses acting on the initial state |g,p . Trajectories of the wavefunction components that differ by the momentum K are indicated.
denotes a motional state of the atom with the average momentump = d 3 p |ψp(p)| 2 p. The AI operation with the pulse sequence π/2-π-π/2 acting on the initial state |g,p is shown schematically in Fig. 2. Since this initial state is a linear superposition of the basis states |g, p , we first calculate the evolution for a basis state with a given p and integrate over the momentum distribution only at the final step when computing the parity expectation value via Eq. (17). Using Eqs. (19)-(21), we find that the evolution operator U tot of Eq. (15) has the following matrix form: and where By noting that we see that φ is the difference between the phases accumulated via the Raman detuning in the second and first halves of the atom's trajectory. The phase φ is zero if the atom moves at a constant velocity. If the atom moves with a constant acceleration a, then δ 12 depends linearly on time. In the laboratory frame, the atom's momentum changes with time as p(t) = p(0) + mat. In the frame that accelerates with the atom, the momentum p is constant, but the Raman frequency experiences a time-dependent Doppler shift: Regardless of which frame is used, the time-dependent part of the Raman detuning is −K · at, which yields Since φ is independent of p, we do not need to know the specific momentum distribution |ψp(p)| 2 of the initial state |ψ(0) to compute the expectation value of the parity operator via Eq. (17). Using Eqs. (17) and (25), we obtain: This result holds for any momentum distribution of the initial state as long as it is sufficiently narrow for the approximation used for Eq. (12) to be valid. For the AI with N independent atoms, the initial state is Note that, for the sake of generality, we assume that the wavefunctionsψp k can be different for different k, although, as we will see below, this assumption does not affect the result for the parity expectation value. The evolution operator for one atom is given by Eq. (24) or, explicitly for the kth atom, by Eq. (A7), where For the kth atom, the phases φ (k) 2T depend on the momentum p k through the Doppler term −p k · K/m in the Raman detuning δ (k) 12 . However, if the kth atom moves with a constant acceleration a k , then, analogously to Eq. (28), φ k = −K · a k T 2 is independent of p k . Therefore, we once again do not need to know the specific momentum distributions |ψp k (p k )| 2 to compute the expectation value of the parity operator via Eq. (17). It also does not matter whetherψp k is the same for all k or not. The derivation of the parity expectation value is described in detail in Appendix A. Specifically, for the initial state of Eq. (30), we obtain: If all atoms experience the same constant acceleration a k = a, then φ k = φ = −K · aT 2 , and we find: The uncertainty of the measured φ value is where we obtain: A scheme of an entanglement-enhanced AI, where the generation of the GHZ state replaces the first π/2 pulse. Trajectories of the wavefunction components that differ by the momentum K are indicated.
It is easy to verify that the phase uncertainty of Eq. (35) is minimized for φ = nπ, where n = 0, ±1, ±2, . . .. Then we obtain: As expected for the AI with independent atoms, this phase uncertainty scales according to the SQL.

B. Ideal interferometer with entangled atoms
In order to surpass the SQL, we propose an AI using an entangled state of atoms. Specifically, we consider a collection of N atoms in the GHZ state [24]: where for α = {g, e} with K g = 0 and K e = K. Normally, the term "GHZ state" is used only for N ≥ 3, however, the state of Eq. (37) is defined for any integer N ≥ 1. For N = 1, this is the single-atom state U 0 (π/2)|g,p = 1 √ 2 (|g,p + |e,p + K ). For N = 2, this is the two-atom Bell state |Φ + [32].
As shown in Fig. 3, the proposed entanglement-enhanced AI scheme differs from the standard operation described in Sec. III A above by the replacement of the first π/2 pulse at time t = 0 by a process that generates the GHZ input state of Eq. (37). The rest of the AI operation, including the π pulse at time t = T and the π/2 pulse at time t = 2T , remains the same. Correspondingly, U T (π), and, using Eqs. (20) and (21), we obtain: The derivation of the parity expectation value is described in detail in Appendix A. Specifically, for the initial GHZ state of Eq. (37) and the evolution operator of Eq. (39), we obtain: where φ k is given by Eq. (31). Since φ k is independent ofp k , the result of Eq. (40) holds for any momentum distribution of the initial state as long as it is sufficiently narrow for the approximation of Eq. (12) to be valid. If all atoms experience the same constant acceleration a, then φ k = φ = K · aT 2 , and we obtain: The phase uncertainty is obtained by substituting Eq. (41) into Eq. (34), which yields, for any value of φ, This is the HL, the ultimate scaling of the phase uncertainty with respect to the number of particles allowed by the laws of quantum mechanics.

IV. EFFECT OF IMPERFECT PREPARATION OF THE INITIAL ENTANGLED STATE
In practice, inevitable noise and imperfections can preclude the achievement of the HL. Our goal is to quantify the deviation from the HL due to these practically relevant errors.
First, we investigate the effect of imperfections in the GHZ state at the input to the AI. We consider two types of state imperfections: a random relative phase between the two components of the GHZ state and the admixture of a noise state. This input state can be represented by the density matrix: Here, |Ψ(β) is the GHZ state with a random relative phase β between the two components, and |ζ is a noise state admixed with probability q ζ (0 ≤ q ζ ≤ 1). With the initial state of Eq. (43) and the same AI operation as described in Sec. III B, i.e., the evolution operator U (k) tot of Eq. (39), the expectation value of the parity operator is (see Appendix A for details of the derivation) where φ k is given by Eq. (31). If we assume, as we have done previously, that all atoms experience the same constant acceleration a k = a, then φ k = φ = K·aT 2 , and N k=1 φ k = N φ in Eq. (46). We also assume that {ϑ k , ϕ k } and β are independent random phase variables. Therefore, we are interested in the expectation value of the parity operator averaged over these random variables, Π . If |ζ represents a contribution from completely random noise, each ϑ k has a uniform distribution with the probability density function P (ϑ k ) = 1/π on [0, π], and each ϕ k has a uniform distribution with P (ϕ k ) = 1/(2π) on [0, 2π]. Since 2π 0 cos(φ − ϕ k )dϕ k = 0 for any value of φ, the contribution from the second term in Eq. (46) is zero. With these assumptions, we obtain: and the amplitude of parity oscillations, Π 0 , is given by where P (β) is the probability density function of a distribution for the random relative phase β. We assume that β has the wrapped normal distribution [62] with zero mean and variance σ 2 β , whose probability density function is given by Consequently, we obtain: If the variance of the relative phase is small, σ 2 β 1, then Eq. (50) can be approximated as Π 0 ≈ (1 − q ζ )(1 − σ 2 β /2), and if both errors are small, i.e., σ 2 β 1 and q ζ 1, then A set of numerical simulations we performed for N = 2 confirms the accuracy of the approximate result in Eq. (51). The fidelity of the initial state is For the initial state of Eq. (43) and using the same averaging over the random variables as in the derivation of Π above, we obtain the average fidelity: Accordingly, the amplitude of parity oscillations can be expressed in terms of the initial state fidelity: The phase uncertainty is obtained by substituting Eq. (47) into Eq. (34), which yields This uncertainty is minimized for the measurement at a dark fringe, where cos(N φ) = 0. Then we obtain: In terms of the initial state preparation error, the phase uncertainty can be expressed as where the approximate result is valid when the errors q ζ and σ 2 β are small. Recent experiments reported the fidelity of the generated two-atom Bell state ranging from 0.89 to 0.97 [28][29][30][31]. This corresponds to Π 0 values from 0.78 to 0.94, respectively. As reported in [34], the fidelity of the generated GHZ state decreases as N increases (specifically, inferred fidelity values are 0.852 for 4 atoms, 0.745 for 8 atoms, 0.643 for 12 atoms, 0.582 for 16 atoms, and 0.542 for 20 atoms). It is expected that the GHZ state fidelity can be increased and generation of even larger GHZ states should be feasible with additional proposed advances [35,36] and technical improvements [34]. However, at this point, it is too early to speculate about how the fidelity will scale with N for very large (100 atoms or more) GHZ states. In any case, the fundamental value of the result of Eq. (56) is the demonstration that, as long as the state preparation errors are small, the deviation from the HL is also small.

V. EFFECT OF LASER INTENSITY FLUCTUATIONS
Random fluctuations in the laser intensity will produce an error in the pulse area: where A 0 is the nominal (errorless) pulse area and δA is the pulse-area error. We assume that δA is a random variable that has the wrapped normal distribution with zero mean, whose probability density function is P (δA) = f WN (δA; 0, σ 2 A ). We denote the errors in the π pulse and the π/2 pulse as w and v, respectively. Using Eq. (18), we obtain the matrices that represent the evolution operators for these pulses: where we used notation C ≡ cos(w/2), S ≡ sin(w/2), and Q ± ≡ cos(v/2) ± sin(v/2). Using Eqs. (58) and (59), we obtain the matrix elements of the evolution operator U With the initial GHZ state of Eq. (37) and the evolution operator matrix elements of Eqs. (60), the expectation value of the parity operator is (see Appendix A for details of the derivation) Under the usual assumption of a constant uniform acceleration, we set φ k = φ = K · aT 2 and N k=1 φ k = N φ in Eq. (61). It then follows from Eq. (61) that the expectation value of the parity operator averaged over the pulse-area error variables, Π , has the general form of Eq. (47), i.e., Π = Π 0 cos(N φ), where the amplitude of parity oscillations, Π 0 , is given by The averages in Eq. (62) are over the wrapped normal distributions with the probability density functions The sums in Eqs. (63) are difficult to compute numerically for large N . However, we can approximate them analytically by using an asymptotic formula for the binomial coefficients [63]: Phase Uncertainty (rad.) which is valid when n is large and m is linear in n. By substituting the asymptotic formula (64) into Eqs. (63) and approximating the sums by respective integrals, we obtain an approximate expression for Π 0 , which is valid for N 1: For the important case of small pulse-area errors such that N σ 2 v 1 and N σ 2 w 1, we expand the exponentials in Eqs. (63) to derive: which yileds another approximation for Π 0 : Analogous to the case of a one-dimensional random walk, the variance σ 2 A is proportional to the pulse duration τ and, correspondingly, to the nominal pulse area A 0 . Specifically, σ 2 v = ξ 2 π/2 and σ 2 w = ξ 2 π, where the parameter ξ is on the order of 10 −3 or better in state-of-the-art experiments. Therefore, Eq. (65) is transformed into which is valid for N ξ 2 π/2 1 (this condition is satisfied for N 10 5 when ξ = 10 −3 ).
The phase uncertainty has the general form of Eq. (54) and, for the measurement at a dark fringe, of Eq. (55). The explicit form of the phase uncertainty is, for N 1, and, for small pulse-area errors (N ξ 2 π/2 1), (71) Figure 4 shows the phase uncertainty as a function of the number of entangled atoms, N , for various values of ξ, with solid lines for ∆φ of Eq. (70) and dashed lines for ∆φ of Eq. (71). The term responsible for the deviation from the HL, N ξ 2 π/2, increases linearly with N . Fortunately, due to the parameter ξ being so small in state-of-the-art experimental conditions, the deviation from the HL is insignificant for at least N 10 5 .

VI. EFFECT OF LASER PHASE NOISE
As shown in Sec. II, the phases of the two Raman fields affect the atomic dynamics only through their difference φ 12 = φ 1 − φ 2 . If the two Raman fields are produced by splitting the light from a single laser, the laser phase noise will mostly cancel out in φ 12 , and the only remaining contribution will be from spurious noise due to technical imperfections, e.g., vibrations and noise in active optical elements such as modulators. While this phase noise could be very small (especially when the two optical paths are well balanced), its effect on the AI performance warrants consideration.
Our analysis is based on using the form (12) for the evolution operator U t (A) and assuming φ 12 = −π/2 + ϑ, where ϑ is a random phase variable that represents the phase noise. Correspondingly, the evolution operators for the two AI pulses have the following matrix forms: where ϑ π and ϑ π/2 denote the random phases in the π and π/2 pulses, respectively. Consequently, the evolution operator U 2T + ϑ π/2 . Therefore, with the initial GHZ state of Eq. (37), the expectation value of the parity operator has the form of Eq. (40) with the replacement φ k → φ k −θ, i.e., whereθ ≡ 2ϑ π − ϑ π/2 . Assuming, as usual, a constant uniform acceleration, resulting in φ k = φ = K · aT 2 , we obtain: Under a general assumption that ϑ π and ϑ π/2 are symmetrically distributed around zero, the expectation value of the parity operator in Eq. (75) averaged over the random phase variables has the general form Π = Π 0 cos(N φ), where Next, we assume that each of the random phase variables ϑ π and ϑ π/2 has the same distribution, which is the wrapped normal distribution with zero mean and variance σ 2 ϑ , whose probability density function is P (ϑ µ ) = f WN (ϑ µ ; 0, σ 2 ϑ ), where µ = {π, π/2}. Thenθ is also a random phase variable that has the wrapped normal distribution with P (θ) = f WN (θ; 0, r corr σ 2 ϑ ), where the value of the numerical factor r corr depends on the degree of correlation between ϑ π and ϑ π/2 . This correlation is determined by the ratio between the characteristic time scale of the spurious noise and the time separation T between the π and π/2 pulses. If ϑ π and ϑ π/2 are completely correlated then r corr = 1, if ϑ π and ϑ π/2 are completely independent then r corr = 5, and in general 1 ≤ r corr ≤ 5. With this assumption, we obtain: where the approximate expression is valid for small phase errors such that N 2 σ 2 ϑ 1.
where the approximation in Eq. (79) is valid for small phase errors (N 2 σ 2 ϑ 1). Figure 5 shows the phase uncertainty as a function of the number of entangled atoms, N , for various values of σ ϑ and r corr = 5, with solid lines for ∆φ of Eq. (78) and dashed lines for ∆φ of Eq. (79). The term responsible for the deviation from the HL, N 2 σ 2 ϑ r corr /2, increases quadratically with N . While we do not know the exact value of σ ϑ , it should be very small based on considerations described above, which ensures that the deviation from the HL is insignificant even for large values of N .

VII. EFFECT OF INITIAL MOMENTUM DISTRIBUTION
In this section we investigate the effect that the momentum uncertainty of the atoms has on the phase uncertainty of the AI. This momentum uncertainty arises because each atom is initially trapped in individual optical tweezers. In the anticipated AI protocol, the tweezers are extinguished, and short pulses of the Raman and Rydberg lasers generate the GHZ spin sate, while also imparting a state-dependent momentum kick to the atoms. The GHZ state preparation time can be made as short as ∼ 1 µs using optimally shaped pulses [34]. Hence we can neglect atomic motion during this step. After the GHZ state has been prepared, the first free evolution step of the AI operation begins, and we assume that the atoms are in the same motional state in which they existed in the traps, plus the state-dependent momentum kick from the entangling procedure.

A. Description of the trapped atoms
The trapped atoms are confined in three dimensions and have three components of vibrational motion, but the only motion relevant to the interferometer is that along the direction of the momentum kick K from the two-photon stimulated Raman transition. This is clear because the atomic momentum p always appears in the dot product p · K in the effective Hamiltonian for the Raman transition, as shown in Eq. (8). Therefore, we consider motion of the atoms in only one dimension, along the coordinate parallel to K. We denote the position and momentum of the kth atom along this coordinate as x k and p k , respectively.
We assume that each atom is initially in its own harmonic trap with energy eigenstates |n k , E n k = (n k + 1 2 ) ω trap , n k = {0, 1, 2, . . .}, where ω trap = 2πν trap is the trap frequency which is assumed to be identical for all traps. The position-space and momentum-space representations of the state |n k are, respectively: ψ n k (p k ) = i n k 2 n k n k ! πσ 2 p −1/2 where H n k is the Hermite polynomial of degree n k . The parameter σ x = /mω trap is the length scale of the trap, and σ p = /σ x = mω trap , where m is the atom's mass. In this paper, we consider 133 Cs atoms for which m ≈ 132.90545 u. The position and momentum uncertainties of the ground vibrational state (for one atom) are (∆x) 0 = σ x / √ 2 and (∆p) 0 = σ p / √ 2. The vibrational motion of the trapped atoms is decoupled from their spin state. As a result, the total density matrix of the atoms has the form: We assume that each atom is in a thermal vibrational state, with temperature T , identical for all atoms. The thermal vibrational state of N atoms has the form: where n = [exp( ω trap /k B T ) − 1] −1 is the average vibrational excitation number, which is identical for all atoms, and z ≡ n /(1 + n ) = exp(− ω trap /k B T ).
As the spins of the atoms can be entangled, we express the spin density matrix in terms of the general N -atom spinor X α : where α stands for the set of indices {α 1 , α 2 , . . . , α N } and α k = {g, e}. Immediately after the traps are switched off, the atoms undergo an entangling process that prepares them in the GHZ spin state and, due to the use of counter-propagating Raman beams, provides a momentum kick of K, where K = |K| ≈ 4π/λ D2 ≈ 1.4743 × 10 7 m −1 , to the excited spin state only. In the GHZ state, the components of X α acquire the following values: X α = 1/ √ 2 if α k = g for all k, X α = 1/ √ 2 if α k = e for all k, and X α = 0 otherwise. The spin-dependent momentum kick for the kth atom is denoted as K α k , where K e = K and K g = 0.
To represent the effect of this entangling procedure, we modify the density matrix ρ vib in Eq. (82) by expanding the energy eigenstates |n k in the basis of momentum eigenstates and shifting their momenta by K α k . Hence we can express the total density matrix at t = 0 as where P (p k , p k ) ≡ p k | ρ (k) vib |p k is given by The sum in Eq. (85) can be evaluated analytically. Using Eq. (80b), we rewrite P (p k , p k ) as is the generating function for the Hermite polynomials [64]. Substituting Eq. (87) into Eq. (86), we derive the following expression for P (p k , p k ):

B. Parity expectation value
Analogously to Eq. (17), the expectation value of the parity operator is where ρ(0) is the initial density matrix, given by Eq. (84), and U tot = N k=1 U (k) tot is the evolution operator for the AI operation. With the initial GHZ state, the AI performs the π-π/2 pulse sequence, and therefore U (k) tot is given by where the evolution operator for each of the pulses, U (k) t (τ ), is given by Eq. (10), where t = T , τ = τ π for the π pulse and t = 2T , τ = τ π /2 for the π/2 pulse.
By substituting the initial density matrix ρ(0) of Eq. (84) into Eq. (89), we obtain: where P (p k ) = P (p k , p k ) = p k | ρ (k) vib |p k is the momentum distribution for the thermal vibrational state of one atom. Using Eq. (88), we obtain: which is the probability density function of the normal distribution for variable p k , with zero mean and variance σ 2 th , where σ th ≡ (∆p) th = p 2 = σ p n + 1/2 is the momentum uncertainty of the thermal state for one atom.
In contrast to the analysis in Secs. III-VI, in this section we do not make the approximation used to obtain Eq. (12), i.e., we do not neglect terms on the order of |δ 12 − δ AC |/Ω eff and |δ 12 |τ in Eq. (10). As mentioned above, in practice, laser frequency chirping is used to compensate the evolving Doppler shift due to the acceleration of the atom, and thereby keep |δ 12 | small enough for this approximation to be valid [39,40]. However, if the momentum spread of the initial state ρ(0) is significant, then it might be impossible to make the aforementioned terms negligible for all momentum components of the atomic wave packet simultaneously. As we mentioned above, in the frame that accelerates with the atom, the Raman frequency experiences a time-dependent Doppler shift: ω 12 (t) = ω 12 (0) − Ka k t, where a k = a k · K/K is the component of a k parallel to K. If one also implements a linear chirp bt of the Raman frequency ω 12 (t), the full time dependence is ω 12 (t) = ω 12 (0)+(b−Ka k )t [39]. Correspondingly, an expression for the Raman detuning that explicitly takes into account the time dependence is It is customary to set the Raman frequency at t = 0 to be ω 12 (0) = ω eg + K 2 /2m (i.e., at the resonance for p k = 0), which yields δ (k) In this section, we assume that the value of the frequency chirp rate b is chosen such that the term (b − Ka k )t is much smaller than σ th K/m at the times when the pulses are applied (i.e., for t = T and t = 2T ), and therefore we set δ We also make the choice φ 12 = −π/2 and, for the sake of simplicity, set δ AC = 0.
By using the evolution operator of Eq. (90) in Eq. (91), the expectation value of the parity operator is (see Appendix A for details of the derivation) where Under the usual assumption of a constant uniform acceleration, i.e., φ k = φ = (b − K · a)T 2 , the parity expectation value of Eq. (94) has the general form Π = Π 0 cos(N φ), where Π 0 = (1 − η) N . In the regime where N η 1, we can use the approximation Π 0 ≈ 1 − N η.
The value of η can be easily computed via numerical integration in Eq. (95). However, it is instructive to obtain an approximate expression for η, which is valid when the uncertainty of the Doppler shift term, (∆δ 12 ) th = Kσ th /m, is small compared to Ω eff . We expand the integrand in Eq. (95) in the powers of r 2 and neglect all terms on the order of r 4 or smaller, to obtain: where κ = 1 2 (π 2 + 3 − 2π) ≈ 3.2932 is a numerical factor. Using the explicit form σ 2 th = mω trap ( n +1/2), we rewrite Eq. (97) as where is the average vibrational energy of the atom in the trap.
In the limit of high temperature, ω trap k B T , we have E vib ≈ k B T , and η is independent of the trap frequency. In the limit of low temperature, ω trap k B T , we have E vib ≈ ω trap /2 (the ground state energy), and η scales linearly with the trap frequency.
The phase uncertainty has the general form of Eq. (54) and, for the measurement at a dark fringe, of Eq. (55). The explicit form of the phase uncertainty is From Eq. (98), the strategy to minimize the error associated with the momentum uncertainty is very straightforward: minimize the vibrational energy by cooling the atoms as close as possible to the ground state and lowering the trap frequency, and increase the Rabi frequency for the Raman transition by using a high-intensity laser with tight focusing. Ideally, we would prefer the regime in which atoms are cooled to submicrokelvin temperatures, the trap frequency is lowered to ν trap < 10 kHz, and the Rabi frequency is increased to Ω eff > 2π × 500 kHz. A recent experiment [31] reported a Rabi frequency of Ω eff ≈ 2π × 250 kHz for a two-photon Raman transition driven by a laser field tuned near the D 1 (5 2 S 1/2 −→ 5 2 P 1/2 ) transition of 87 Rb. Since the dipole moment for the D 2 (6 2 S 1/2 −→ 6 2 P 3/2 ) transition of 133 Cs is about 1.5 times larger than that for the D 1 transition of 87 Rb, one can expect that a laser system with the same intensity and focusing as the one used in Ref. [31] would produce a Rabi frequency of Ω eff ≈ 2π × 560 kHz for a two-photon Raman transition driven by a laser field tuned near the D 2 transition of 133 Cs.

C. Numerical results
To differentiate between exact (numerical) and approximate (analytical) results, we denote the values of η obtained via numerical integration in Eq. (95) as η num and the values calculated using the approximate analytical formula in Eq. (98) as η approx . The relative error of the approximate value of η is |η approx − η num |/η num . Similarly, we denote the values of the phase uncertainty obtained by substituting η num into Eq. (100a) as (∆φ) num and the values calculated using the approximate analytical formula in Eq. (100b) as (∆φ) approx . The relative error of the approximate value of the phase uncertainty is |(∆φ) approx − (∆φ) num |/(∆φ) num . Figure 6 shows η num and |η approx − η num |/η num as functions of the Rabi frequency Ω eff /2π for various values of the atom temperature T and a fixed value of the trap frequency (ν trap = 10 kHz). As expected from Eq. (98), the scaling is η ∝ 1/Ω 2 eff . Figure 7 shows η num and |η approx − η num |/η num as functions of the trap frequency ν trap for various values of the atom temperature T and a fixed value of the Rabi frequency (Ω eff /2π = 300 kHz). As expected from Eq. (98), η is independent of the trap frequency in the high-temperature regime ( ω trap k B T ) and scales linearly with ν trap in the low-temperature regime ( ω trap k B T ). We observe that the value of the relative error is very close to the value of η num and scales in exactly the same way. This is not surprising, since the leading term in |η approx − η num | is on the order of r 4 and hence |η approx − η num |/η num ∼ r 4 / r 2 ∼ r 2 ∼ η num . Figure 8 shows the phase uncertainty (∆φ) num and the relative error |(∆φ) approx − (∆φ) num |/(∆φ) num as functions of the number of entangled atoms, N , for three different parameter combinations that are described in Table I. These three parameter combinations correspond to η ≈ 5.0 × 10 −3 (yellow curve), η ≈ 2.0 × 10 −3 (red curve), and η ≈ 5.0 × 10 −4 (blue curve). Figure 8(b) demonstrates that the difference between the numeric and approximate values of ∆φ is comfortably small for all considered values of N .  To quantify how close a phase uncertainty curve is to the HL, we use the number N * , which is defined as the number of atoms, for which the deviation from the HL is less than 10%, i.e., N * is the largest N for which (1 − η) N ≥ 0.9. The values of N * for the curves in Fig. 8 are reported in Table I. We see that, with optimistic parameter values, the deviation from the HL is small for N ∼ 100.

A. Parity measurement protocol
The analysis of the AI performance in this paper is based on the assumption that the expectation value of the parity operator Π of Eq. (16) is measured for a system of N atoms. Therefore, an important question is how the number of measurements required to evaluate Π scales with N . We present here a protocol for parity measurement that scales linearly with N despite the fact that the Hilbert space dimension grows exponentially.
The Hilbert space of a system of N two-level atoms is spanned by the set S of 2 N basis states, which are defined in Table II.
Here, M j denotes the number of atoms in the |e level for the state |j . From Table II, it is easy to see that the entire set S consists of N + 1 subsets S M such that the expectation value of the parity operator is Using the partition of the basis states |j into the subsets {S M |M = 0, 1, . . . , N }, we obtain: where is the probability that the system is in the subset S M of states.
Consider now a state-selective measurement, for example, via state-dependent fluorescence imaging (SFI) [59,60] or via coherent spatial splitting (CSS) in a state-dependent optical lattice [61], which detects whether an atom is in the |e level or in the |g level. For example, in SFI, by driving a resonant cycling transition, fluorescence is produced if the atom is in the |e level and no fluorescence is produced if the atom is in the |g level. Similarly, in CSS, by applying a sequence of pulses and lattice transformations, atoms in the |e level will be shifted to the left and atoms in the |g level will be shifted to the right. Hence, for a system of N atoms, any state in the subset S M will produce M fluorescence images (in SFI) or M atoms shifted to the left (in CSS). Correspondingly, for an arbitrary state |ψ , the probability to detect M fluorescence images (in SFI) or M atoms shifted to the left (in CSS) is P M . Therefore, the measurement of the N + 1 probabilities {P M } directly yields the parity expectation value via Eq. (105).

B. Measurement error
For this measurement protocol, a conservative estimate of the effect of measurement error is based on assuming that a random wrong value of Π is obtained if at least one atom is detected in a wrong level (|g instead of |e or |e instead of |g ). By averaging over all these error outcomes, we obtain: where q det is the probability of erroneous state detection for one atom. Correspondingly, the phase uncertainty, minimized for the measurement at a dark fringe, where cos(N φ) = 0, is Equations (107) and (108) comply, respectively, with the general forms Π = Π 0 cos(N φ) and ∆φ = (Π 0 N ) −1 , with Π 0 ≈ (1 − q det ) N . A further approximation Π 0 ≈ 1 − N q det holds when N q det 1. In an SFI experiment with 10 atoms, an average state detection fidelity of 0.987 was reported [60], which corresponds to (1 − q det ) N ≈ 0.877 for N = 10. While fidelity scaling to a much larger N is not yet known in SFI experiments, the same fidelity for N = 100 would result in (1 − q det ) N ≈ 0.27. In a CSS experiment with 160 atoms, an average state detection fidelity of 0.9994 was reported [61], which corresponds to (1 − q det ) N ≈ 0.94 for N = 100. Furthermore, in CSS, fidelity is essentially independent of the number of atoms measured, and the measurement causes negligible atom loss. However, in order to perform a CSS experiment, atoms need to be loaded into a three-dimensional optical lattice. Figure 9 shows the phase uncertainty ∆φ of Eq. (108) as a function of the number of entangled atoms in the GHZ state, N , for q det = 1.3 × 10 −2 (corresponding to the average state detection fidelity of 0.987 reported in the SFI experiment [60]), q det = 5.0 × 10 −3 (a hypothetical value), and q det = 6.0 × 10 −4 (corresponding to the average state detection fidelity of 0.9994 reported in the CSS experiment [61]). If the state detection error is sufficiently small (one per cent or lower), the phase uncertainty surpasses the SQL for N ∼ 100.

IX. EFFECT OF ATOM LOSS
The loss of an atom during the AI operation can be detected, and the outcomes of the respective experiments can be eliminated from the data. Hence, the possibility of losing an atom leads to a reduced data-acquisition rate. At the end of our AI protocol, each atom is imaged individually by SFI [59,60], and the number of atoms in the bright state is recorded. Atom loss is then detected by transferring the dark-state atoms to the bright state, repeating the SFI measurement, and comparing the total number of atoms detected with the initial number of atoms. Therefore, post-selecting for the measurement outcomes without atom loss eliminates the error introduced by counting a lost atom as a dark-state measurement.
This method of post-selection of lossless outcomes is not always possible in other systems. For example, in optical interferometry the loss of a single photon transforms the maximally entangled NOON state into an incoherent mixture [12,65]. Without the ability to discriminate between lossless and lossy operations, the outcomes of measurements corresponding to incoherent mixture states add to the phase uncertainty of an optical interferometer.
The effect of a reduced data-acquisition rate due to atom loss on the phase uncertainty can be quantified as follows. If an AI experiment with the N -atom GHZ state is repeated M times, then, in the ideal case, the phase uncertainty after M experiments is If the loss probability for one atom in one AI experiment is q loss , then the number of post-selected lossless outcomes is M (1 − q loss ) N . The resulting phase uncertainty after M experiments is From Eq. (110), the scaling of the phase uncertainty per one experiment, ∆φ, with N can be interpreted differently depending on the type of application, for which the AI is used. For example, in a stationary gravity measurement, it might be possible to repeat the experiment as many times as needed to achieve the desired number of lossless outcomes. This can be interpreted as if the atom loss prolonged the effective time that it takes to perform one lossless experiment but did not affect the scaling of ∆φ. However, in many cases (in particular, when the AI is used in an inertial navigation system), the total measurement time and hence the total number of experiments (AI cycles) are fixed. In such a case the result of Eq. (110) can be effectively interpreted as the phase uncertainty for M lossless experiments with the uncertainty per one experiment given by   Figure 10 shows the phase uncertainty ∆φ of Eq. (111) as a function of the number of entangled atoms in the GHZ state, N , for q loss = 2.0 × 10 −2 , q loss = 1.0 × 10 −2 , and q loss = 0.5 × 10 −2 . We see that if atom loss is sufficiently low (two per cent or smaller), the phase uncertainty surpasses the SQL for N ∼ 100.

X. SUMMARY AND CONCLUSIONS
Within the protocol for AI operation that we have described, non-negligible errors come from the imperfect initial state preparation, intensity and phase fluctuations of the Raman beams, the initial momentum spread of the atoms in the optical traps, imperfect measurement, and atom loss. A key finding is that each of these errors preserves the general forms for the parity oscillations, Π = Π 0 cos(N φ), and for the minimum phase uncertainty, ∆φ = (Π 0 N ) −1 , and manifests itself through a reduction of the amplitude Π 0 below the ideal value of Π 0 = 1. We have derived analytical results that express the dependence of Π 0 on error parameters and N for all of these error sources. When an error is small, Π 0 is close to 1 and can be approximated as Π 0 ≈ 1 − ε, where ε 1. The obtained results indicate that an entanglement-enhanced AI with a phase uncertainty close to the HL is feasible for N ∼ 100 with state-of-the-art experimental capabilities.
If we use an index to enumerate different uncorrelated error sources, then the total effect of all errors on the amplitude of parity oscillations is given by When all of these errors are small (ε 1 ∀ ), the total error is simply a sum of all the individual errors: ε total ≈ ε . Table III reports a summary of all the non-negligible error sources considered in this work. This table lists exact and approximate (valid for small error) expressions for Π 0 for each error source, along with descriptions of respective error parameters. Currently, error in the initial state preparation is not prohibitive for N 20. For this error, an explicit scaling with N is not available, and generating larger high-fidelity entangled states is an active field of study. Errors arising from intensity and phase fluctuations of the Raman beams are insignificant for large numbers of atoms (N 10 4 ), assuming stateof-the-art optical technology. The error due to the initial momentum distribution of the atoms scales linearly with the average vibrational energy of the atom in the trap and inversely proportional to the square of the effective Rabi frequency for the two-photon Raman transition. For optimistic (yet still realistic) parameter values (atoms cooled to sub-microkelvin temperatures in shallow traps and driven with very intense Raman fields), this error can be made insignificant for N 100. Optimal shaping of Raman fields is potentially useful for mitigating this error [66,67]. Measurement error and detection loss can also be significant for large numbers of atoms, but attaining a phase uncertainty near the HL still appears to be feasible for N ∼ 100. The CSS method of state-selective measurement, while challenging in an AI setup, promises superior scaling with the number of atoms [61].
Interferometers with highly entangled atomic spin states show great promise as inertial sensors, and an array of singleatom optical tweezers provides a suitable platform for these sensors. When the atoms are cooled to ultracold temperatures and initiated in the GHZ state, only a few error sources stand in the way of reaching a phase uncertainty near the HL. We have shown how these errors scale with the number of entangled atoms and identified parameters that quantify these errors for a realistic system of entangled Cs atoms. III. Summary of error sources that affect interferometry with entangled atoms. Each error source results in a reduction of the parity oscillation amplitude, Π0, below the ideal value of Π0 = 1. N is the number of atoms in the initial GHZ state. fWN(θ; 0, σ 2 ) is the probability density function of the wrapped normal distribution for variable θ, with zero mean and variance σ 2 .

Error source Π0
Approximate Π0 for small error

Initial momentum spread
(1 − η) N 1 − N η η -parameter that quantifies the error due to the initial momentum spread of trapped atoms; the exact expression for η is given by Eq. (95) and an approximate expression (valid for small errors) is given by Eq. tot for an entire sequence of the AI pulses acting on the kth atom, in the basis {|g, p k k , |e, p k + K k }. Note that we omitted the time dependence from the basis states, since, according to Eq. (17), the field-free evolution has no effect on the parity expectation value. In this basis, U (k) tot has the general matrix form: (A1) In general, the matrix elements in Eq. (A1) depend on the phases φ (k) 2T which, in their turn, depend on the momentum p k through the Doppler term −p k · K/m in the Raman detuning δ (k) 12 experienced by the kth atom. Following Eq. (13), the evolution operator for the system of N atoms is given by the tensor product of all one-atom evolution operators: First consider an initial state for N atoms which are not entangled. The general form for an unentangled state is where c   |α,p k + K α k = d 3 p kψp k (p k )|α, p k + K α k , (A3) for α = {g, e} with K g = 0 and K e = K. The output state |Ψ out = U tot |Ψ in is given by By substituting this result into Eq. (17) and using the unitarity of the matrix U (k) tot , we obtain: where we introduced the notation: for averaging over the momentum distribution of the kth atom.
For the three-pulse sequence π/2-π-π/2, the evolution operator for one atom is given by Eq. (24) or, explicitly for the kth atom: Using the matrix elements from Eq. (A7), we obtain: Since these terms depend only on φ k , and φ k is independent of p k , each of the averages in Eq. (A5) is equivalent to multiplication by 1, regardless of the specific forms of the momentum distributions |ψp k (p k )| 2 . Therefore, by substituting Eqs. (A9) into Eq. (A5), we obtain: In particular, if all atoms are initially prepared in the ground state, c For the two-pulse sequence π-π/2, the evolution operator acting on the kth atom is given by Eq. (39). Using the matrix elements from Eq. (A7), we obtain: Once again, these terms are independent of p k , and each of the averages in Eq. (A5) is equivalent to multiplication by 1, which produces In particular, for the state |ζ of Eq. (45), the coefficients are c where c g and c e are complex coefficients (subject to the normalization condition |c g | 2 + |c e | 2 = 1). The output state |Ψ out = U tot |Ψ in is given by × U (k) gg |g, p k k + U (k) eg |e, p k + K k + c e N k=1 d 3 p kψp k (p k ) × U (k) ge |g, p k k + U (k) ee |e, p k + K k . (A16) By substituting this result into Eq. (17) and using the unitarity of the matrix U (k) tot , we obtain: For the two-pulse sequence π-π/2, we use the terms in Eqs. (A12) and substitute them into Eq. (A17) to obtain: In particular, for the GHZ state of Eq. (37), c g = c e = 1/ √ 2, which yields Π = cos N k=1 φ k . Similarly, for the state |Ψ(β) of Eq. (44), c g = 1/ √ 2 and c e = e iβ / √ 2, which yields Π = cos N k=1 φ k − β . Next, we consider the two-pulse sequence π-π/2 in the presence of pulse-area errors due to laser intensity fluctuations. In this case, U Substituting these terms into Eq. (A17), we encounter factors These averages are integrals over very rapidly oscillating functions and therefore they are extremely small. Specifically, for an atom in a harmonic trap, the averages in Eqs. (A20a) and (A20b) scale as e −γ and the one in Eq. (A20c) scales as e −4γ , where γ = (KT ∆p/m) 2 and ∆p is the atom's momentum uncertainty. Even with the minimum uncertainty, (∆p) 0 = mω trap /2, we find γ ∼ 10 3 for typical values of T (∼ 1 ms) and ν trap (∼ 10 kHz). After neglecting all terms that include these extremely small factors, we obtain: Π = 2 cos N v cos 2N (w/2)Re c * g c e e −i N k=1 φ k + |c g | 2 + (−1) N |c e | 2 sin N v cos N w.
Since the pulse-area error v is small, v 1, the second term in Eq. (A21) includes a factor that scales as v N . This factor is extremely small in the regime of large atom numbers, N 1, in which we are interested, and therefore the second term in Eq. (A21) can be safely neglected. Then we find: In particular, for the GHZ state of Eq. (37), c g = c e = 1/ √ 2, which yields Π = cos N v cos 2N (w/2) cos N k=1 φ k .
Finally, we consider the case where the pulse detuning error due to the initial momentum uncertainty of the atoms is taken into account, as described in Sec. VII. In this case, the initial state of the system is described by the density matrix ρ(0) of Eq. (84), and the parity expectation value is given by Eq. (91). With the coefficients X α = c g if α k = g for all k, X α = c e if α k = e for all k, and X α = 0 otherwise, Eq. (91) takes the form of Eq. (A17), but now the average over the momentum distribution of the kth atom is where P (p k ) = p k | ρ (k) vib |p k is the momentum distribution for the thermal vibrational state of one atom, given by Eq. (92). Since this momentum distribution is the same for all atoms, in what follows we will rename the integration variable from p k to p in the averages of the form (A24), to simplify the notation.
In deriving Eqs. (A25), we neglected all terms that include the factors shown in Eqs. (A20) since, as we showed above, these averages are integrals over very rapidly oscillating functions and therefore they are essentially zero. It is easy to see that, for small r, the average in Eq. (A25b) scales as |U (k) gg | 2 − |U (k) eg | 2 ≈ −(1 − π/4) r 2 , where r 2 = (Kσ th /mΩ eff ) 2 , and its contribution to Π [i.e., the second term in Eq. (A17)] scales as (1 − π/4) N r 2 N . Since we are interested in the parameter regime, in which r 2 1 and N 1, the contribution of this term to Π can be safely neglected. Also, the imaginary part of the average in Eq. (A25a) is an integral over an odd function of p and therefore it is equal to zero. Using these facts, we obtain: where η ≡ 1 − sin 3 2λ (1 + r 2 ) 3/2 cos πr + r tan λ sin πr √ 1 + r 2 . (A28) In particular, for c g = c e = 1/ √ 2, Eq. (A27) yields fringes, which needs to be taken into account for an accurate measurement of φ. The ratio is independent of . Typically, θ φ, for example, with Ω eff ∼ 2π × 10 5 Hz and T ∼ 10 −3 s, we find θ/φ ∼ 10 −3 . This ratio can be further reduced by increasing Ω eff (which is also useful for reducing η and η ) and/or T .