Many-body effects in quantum metrology

We study the impact of many-body effects on the fundamental precision limits in quantum metrology. On the one hand such effects may lead to non-linear Hamiltonians, studied in the field of non-linear quantum metrology, while on the other hand they may result in decoherence processes that cannot be described using single-body noise models. We provide a general reasoning that allows to predict the fundamental scaling of precision in such models as a function of the number of atoms present in the system. Moreover, we describe a computationally efficient approach that allows for a simple derivation of quantitative bounds. We illustrate these general considerations by a detailed analysis of fundamental precision bounds in a paradigmatic atomic interferometry experiment with standard linear Hamiltonian but with both single and two-body losses taken into account---a model which is motivated by the most recent Bose-Einstein Condensate (BEC) magnetometry experiments. Using this example we also highlight the impact of the atom number super-selection rule on the possibility of protecting interferometric protocols against decoherence.


I. INTRODUCTION
The research effort in the field of quantum metrology [1][2][3][4] goes along two main directions. The first is the experimental one, with the ultimate goal of designing metrological setups that harness quantum behavior of light and matter while at the same time minimize unwanted effects of decoherence in order to reach the precision regime inaccessible within the classical paradigm [5][6][7][8][9]. The theoretical effort, on the other hand, aims at identifying new promising proposals for quantum metrological protocols as well as fundamental bounds on precision that can in principle be achieved in quantum systems. The bounds are derived using simplified models that are intended to capture the essential features of the relevant physical systems and take into account the dominant sources of decoherence [10][11][12][13][14][15][16][17].
Typically, a significant discrepancy is observed between the fundamental bounds arising from theoretical considerations and the performance of practical setups. This is due to either oversimplified theoretical modeling of important physical effects present in the experiment or experimental inability of preparing the optimal quantum states and implementing the optimal quantum operations required by the theory. A noticeable exception from this rule is the case of quantum enhanced optical interferometry, where the performance of squeezed-light enhanced gravitation wave detector devices [18,19] approaches closely the fundamental limit predicted by the theory [20].
Atomic systems are inherently more complex than purely optical ones, and hence the theoretical analysis of fundamental metrological bounds is much more challenging. While in optical interferometry the dominant source of decoherence is loss, which has single-body character, in atomic systems interactions between atoms naturally lead to various manybody effects. In cold atom experiments, apart from single-particle losses, which in principle might be described using a similar formalism as in the optical case, two-and three-body losses play a crucial role in the system dynamics. This poses a serious challenge to theorists, as the most powerful theoretical methods developed in the field of quantum metrology are naturally tailored to deal with uncorrelated, or in other words, single-body noise [12,13,21].
In this paper, we substantially develop the ideas from [16], where it was proposed to view the evolution of an N -particle permutationally invariant system with at most n-body interactions as equivalent to an application series of n-particle quantum channels to n-element subsets of N particles. In this way, one is able to make use of theoretical methods, developed with single-body noise models in mind, in many-body noise cases, without the need to consider the full N -particle space but rather a much smaller n-particle space. In what follows we will call this the reduced particle number (RPN) approach. In [16] this idea was used to obtain scalings of precision bounds in cases of a k-body Hamiltonian and l-body noise processes. Here we generalize these considerations to a situation where both Hamiltonian and the noise part contains many terms with different level of non-linearity and provide a clear recipe for how to proceed in order to obtain the eventual fundamental precision scaling. Furthermore, we discuss the way to obtain explicit quantitative bounds using this approach, and apply it to find limits on precision of atomic interferometry in presence of both arXiv:1810.03651v2 [quant-ph] 17 Oct 2018 single-and two-body losses.
Interestingly, according to [16,17] in case of a linear Hamiltonian the two-body losses are in principle correctable, and therefore one should only focus on the effects of single-body losses. This is true, but we show that this is true only under an assumption that one is able to prepare states with indefinite particle numbers. Conversely, in case one is limited to definite particle number states, as is the case of atomic interferometry where the particle number super-selection rule holds, this statement is no longer valid. We show this fact analytically, and also derive an explicit bound on precision of linear interferometry in presence of two-body losses. We also show how this super-selection rule should be applied on the level of the RPN approach to obtain explicit quantitative bounds that include the effects of both single and two-body losses that are meaningful for atomic interferometric experiments. Finally, we compare the obtained bounds with simulations of an explicit atom interferometric protocol as well as data from a Bose-Einstein Condensate (BEC) magnetometry experiment.
The paper is organized as follows. In Section II we set up the metrological model that serves as a basis for our considerations and review the fundamental tools of quantum estimation and metrology theory that we use in our analysis. In Section III, we introduce the RPN approach and provide a general recipe to obtain asymptotic scaling of precision with the number of atoms in models where multiple processes with different degrees of non-linearity act simultaneously. We also show how quantitative bounds can be obtained within the RPN approach. In Section IV, we focus on an atomic interferometry model with linear Hamiltonian in presence of both single-and two-body losses. We show similarities and differences between the result obtained within the RPN approach and the approach that is based on the study of the properties of operators acting in the full N -particle bosonic space, and point to the consequences of application of the atom number super-selection rule. We compare the bounds with numerical simulations and experimental data. In Section V we conclude the paper and provide an outlook on possible future research directions.

II. METROLOGICAL MODEL
Let ρ describe the state of an atomic system undergoing an evolution governed by a general quan- The general scheme of an adaptive quantum metrological experiment assuming arbitrary fast control. E ω t represents physical evolution of the N sensing atoms, while Ui are the control operations that may entangle/disentangle the atoms with arbitrary number of ancillary systems. The goal is to estimate ω with lowest possible uncertainty under given total time T of the experiment. tum master equation [22]: where ω, multiplying the generator of the unitary part of the dynamics H (which in what follows we will refer to as the Hamiltonian), is the parameter to be estimated, while L(ρ) represents decoherence processes, which are assumed here to be Markovian and have no inherent time dependence-formally this means that the dynamics can be described using a semi-group. The goal is to estimate ω under a fixed total probing time T which is treated here as a resource. This problem can be viewed as a generalized quantum frequency estimation problem [11,[23][24][25]. Let E ω t be the quantum map representing the above dynamics integrated over time t. Hence, if no other operations are applied to the state we can write the evolved state as ρ ω t = E ω t (ρ 0 ), where ρ 0 is the initial state. Similarly as in [15-17, 26, 27] we look for the fundamental bound and therefore consider the most general quantum adaptive metrological protocol, which is depicted in Fig. 1.
In the scheme in Fig. 1 we have explicitly indicated N lines entering each channel, which represent N atoms being used for sensing. The general scheme allows the atoms to be entangled with an arbitrary number of ancillary systems, which may contribute nontrivially to the sensing process thanks to the unitary control operations U i -we denote the state of the system + ancillas by in contrast to ρ which describes the atomic system alone. The scheme, with N separate lines, suggest that we deal with distinguishable particles. Indeed, our method based on the RPN approach refers to particles as formally distinguishable. Many important atom interferometric experiments, however, are performed with Bose-Einstein condensates, with only the symmetric states involved. Still, since distinguishable particles in principle offer a greater metrological potential, as symmetric states form only a subset of all N particle states, any fundamental bounds derived for distinguishable particles will also be valid for indistinguishable ones. Moreover, these bounds often prove to be saturable by strategies based on the use of symmetric states [2] and hence tightness of the bounds is typically not affected by turning to distinguishable particle paradigm.
Note also, that in the scheme the number of particles N is unchanged from one protocol step to the other. When we deal with losses, this should be understood in the spirit that there are additional internal degrees of freedom of each particle-line that represent the particle being lost. Hence, while the formal number of particles remains unchanged, the loss process can still be described within this framework. This also points to the possibility that control operations U i may effectively "reload" new particles in place of the old lost ones, by a formal change of internal states of each particle-line. We discuss this issue in detail in Sec. IV, where we also show how one can tighten the bounds by assuming impossibility of reloading the lost particles on the course of estimation process.
The bound on precision of estimating ω can be obtained using the Quantum Cramér-Rao bound [28]: where F Q is the Quantum Fisher Information (QFI), dot signifies d dω , while |a , λ a are eigenvectors and eigenvalues of the final state of the protocol ω T . Taking QFI as a figure of merit, the goal of determining the fundamental bound can be viewed as the problem of maximizing QFI over input state and control operations U i . In particular the key qualitative question is whether it is fundamentally possible to protect coherence of the sensing system on long time scales keeping the T 2 scaling of QFI-the so called Heisenberg scaling achievable in noiseless quantum metrology-or is the scaling of the QFI bound to be linear in T for sufficiently long times.
This in itself is an extremely complicated mathematical problem, mainly due to the fact that the presence of decoherence causes the state ω T to be mixed. Fortunately a powerful theoretical methods based on the idea of minimization over different representations of quantum channels originated in [29] and developed in [12,13,21,26] has finally led to a concise and efficiently computable solution of the problem [16,17], where it was shown that provided the following condition is satisfied where H , AH denote the Hermitian and anti-Hermitian part of an operator, QFI scales at most linearly with T . We will refer to this condition as the Hamiltonian in the Lindblad Span (HLS) condition. In this case the explicit upper bound reads [16,17]: where h is a real variable, h a complex vector of length J (J is the number of noise operators L j appearing in Eq. (1)), h is a J × J Hermitian matrix, · is the operator norm, while operators α and β read where the noise operators have been collected in an operator valued vector L = [L 1 , . . . , L J ] T . If, on the other hand, the HLS condition is not satisfied, one can construct a quantum error-correction protocol that effectively protects the system from decoherence while preserving some part of the unitary dynamics, which results in an effective T 2 scaling of QFI [17]. In most realistic cases, such as cold atom metrology with losses, the HLS condition is satisfied and therefore in this paper we will focus mostly on the bound (4).
If the structure of H and L is simple enough the above formulation often allows to find an analytical form of the bound [12,13,16,21]. Otherwise, provided the number of noise operators L is sufficiently small the above problem can be cast as a simple semi-definite program and the bound can be found numerically [16].
Since we will typically be interested in the large particle number limit, both H and the noise operators L will act on a very large Hilbert space and hence in general will not be directly suited to plug into e.g. the semi-definite program mentioned above. However, in cases where we deal with permutationally invariant systems (such as BEC), where both H and L can be described via at most k-body interactions, one can effectively calculate the bound using only properties of a subchannel that describes the dynamics on a k-particle subset of all the particles. Since the dynamics of most cold atom systems is effectively described with terms that never involve k larger than 3-4, this method will be very efficient numerically and may also help to obtain an analytical form of the bound. This approach was introduced in [16], and we will refer to it here the reduced particle number approach. In the next section we will describe this approach systematically and in full generality and show how to obtain the fundamental precision scalings in the general quantum metrology models with many-body effects-the task that was only sketched and studied in some special cases in [16].

III. PRECISION BOUNDS VIA THE REDUCED PARTICLE NUMBER (RPN) APPROACH
Consider a permutationally invariant system of N particles, where the Hamiltonian H and the noisy part of the dynamics L can be split according to different level of non-linearity of interactions they represent: where H (k) N part contains a contribution to the Hamiltonian arising from k-body interactions, and similarly L (l) N represents the l-body contribution to the noisy part. More explicitly, permutational invariance of the problem implies that H . . , i k )} represents all k element combinations of the N element set, and the operator index ν ∈ Υ k N denotes the set of particles that a given operator acts on. Similarly we can write L (l) ν represent the noise part acting on an l-particle subset ν.
Since, a duration of a single step in the most general protocol as depicted in Fig. 1 can without loss of generality be chosen as arbitrarily small, we may equivalently view the action of H and L as a sequence of operations acting on particular sets of particles. According to the Trotter expansion this may at most introduce errors of the order t 2 . Consider for example a situation where we have a threebody H = H   N . We can model the dynamics as effectively consisting of subchannels ε ω t acting on all 3-particle subsets of all N particles, as depicted in Fig. 2.
The tilde symbols above the H and L in Fig. 2 are to remind that due to combinatorial reasons it may be necessary to rescale the operators describing the evolution. Since we consider a three-particle elementary channel, the three-body Hamiltonian is properly accounted for, there is no need for rescaling  term acting on particles ν = (1, 2, 3) as depicted in the picture, and we omitted the ν subscript for brevity. In what follows whenever we write H (k) or L (l) without any subscript we understand this as a term representing action of Hamiltonian or noise on k or l particle subsets. Now, if we consider L (1) term, we see that the number ofL gates is 3 1 N 3 , and since there is N = N 1 uses of the L (1) operator in the proper evolution we need to rescale and set We can now formulate a general recipe. When considering an RPN model with n-particle subchannels, we need to replace the original H (k) operators with and similarly for L (l) ν . Note that, since the noisy part L (l) ν is rescaled by a factor of χ l the corresponding noise operators L (l) ν,j that enter quadratically in L (l) are rescaled as:L where ν is a vector containing all noise operators L (l) ν,j of a given many-body character acting on a given subset of particles ν.
We are now ready to describe a general procedure that allows us to obtain fundamental bounds on precision of estimating ω in such models. First, we make the assumption that neither H (k) nor L (l) has any dependence on N , so in other words the rates of elementary processes do not depend on the total number of particles. We will relax this assumption in Sec. IV, when discussing the effect of single and two-body losses, where we will take into account the possibility of dependence of two-body loss coefficient on the total particle number via effective changes in the size of the atomic cloud in the trap with the increasing number of atoms. (3) is not satisfied, this implies that there is at least one operator H (k) N for some k which is not in space S. Hence, one may in principle remove all the noise present in the system via appropriate error-correction protocols and be left with the N 2k T 2 scaling, known from the models of noise-less non-linear metrology [30][31][32][33][34][35][36][37][38]. Motivated by the physical situation in cold atomic systems, in this paper we assume that this is not the case and all H (k) N satisfy the HLS condition. Let us reconsider now the HLS condition in the context of the RPN approach. When considering an n-particle subchannel, as in Fig. 2, we can formulate an analogous HLS condition for this subchannel, where the H (k) and L (l) terms appear in different variants where they act on different k or l element subsets of n particles. For example, in the case depicted in Fig

First of all, if the HLS condition
We will refer to the HLS condition on the level of n-particle sub-channels as the HLS n condition.
Formally, the effective Hamiltonian in the RPN approach can be written as, while the set of all noise operators combined in a single vector reads: Hence, the HLS n condition is given by the same formula as the standard HLS condition given in Eq. (3), but with H replaced byH n and the space S replaced by S(L n ), for which the set of noise operators is taken from the vectorL n . If the HLS n condition is satisfied for some n then this implies that the HLS condition is satisfied for the whole system. This simply stems from the fact that full H and noise operators L j consist of all possible variants of elementary operators acting on different subsets of particles. Therefore if there is linear dependence between operators within a given sub-channel then there will still be a linear dependence if we add all their variants acting on all different subsets. The opposite, however, may not necessarily be true. It might be the case that for some particular n the HLS n is not satisfied while if we look at the whole system it is. This may be due to the fact that non-trivial new operators appear when we consider products of elementary noise operators L (l) ν,j , L (l ) ν ,j with ν and ν only partially overlapping; since such products enter into the definition of the S(L n ) space, they may help to span the relevant Hamiltonian. Still, there is no point in increasing the space of the sub-channel too much, as no qualitatively new operators will appear in the definition of S(L n ), if we take n > 2l-all possible overlaps of noise operators are then taken into account. Hence when considering the HLS condition using sub-channels one may always limit the considerations to max(k, l) ≤ n ≤ max(k, 2l)-to be precise, here by k, l we mean the maximum values of k and l that are relevant in the considered model.
In order to obtain a quantitative bound within the RPN approach using Eq. (4), note that the number of sub-channels that appear in the decomposition of E ω t as depicted in Fig. 2 is N n , and therefore recalling the reasoning that leads to derivation of the bound (4) [16,17]-where the number of channel uses is entered as the proportionality factor-the corresponding bound reads: This problem can now be cast as a semi-defnite program as described in [16] and solved efficiently provided that n is reasonably small. First we construct the matrix where J is the length of theL n vector. Minimizing the operator norm α n is now equivalent to minimizing λ subject to A n ≥ 0 with the additional constraint coming from the equationβ n = 0. The bound on QFI can therefore be written as which is indeed a semi-definite program. Note that the earlier argument stating that increasing n above max(k, 2l) is not necessary in order to determine the validity of the HLS condition, has to be carefully reanalyzed when thinking of the actual quantitative form of the bound. Indeed, increasing n above max(k, 2l) will not change the scaling character of the bound, but may sometimes allow us to tighten the coefficient appearing in the bound, see the example of two-body losses in Sec. IV C.
A. Scaling in the limit of large particle number Without going into numerics, we can get a deeper understanding of the qualitative behavior of the above bound if we consider the asymptotic limit of large number of particles N → ∞. Recall that according to Eq. (7), the relevant operators scale as H , hence the terms with larger non-linearity are dominating.
First, for simplicity let us assume thatH n consist of only k = k * -body operatorsH . Let us also assume that we can satisfy the HLS n condition using only noise operators with a given non-linearity l = l * . Then according to Eqs. (12), in order to satisfyβ n = 0, the correspondingh coefficients need to scale as N k * −l * /2−n/2 (in order to haveh †L n to scale as N k * −n -the same scaling as theH n ), while according to the same reasoningh needs to scale as N k * −l * andh needs to scale as N k * −n . Consequently, the resultingα n will scale as N 2k * −l * −n and finally the bound (11) implies From the above reasoning it follows that ifH n can be spanned by operators with a given type of nonlinearity l i on their own, then we get the tightest bound, in the sense of scaling with N , by considering the noise operators with the highest nonlinearity. Hence in the above formula we can write l * = max({l i }), where the set {l i } contains nonlinearity levels, where for each of them there are corresponding noise operators that allow to satisfy the HLS n condition using just operators of a given non-linearity type. Note that when we have an interplay of decoherence processes with different non-linearities and we investigate the overall scaling of the bound N 2k * −l * we are left with decoherence processes with relative strength scaling as √ N li−l * . The tightest bound on QFI will come from α n dependent on the strongest decoherence; where by strongest we mean such that scales with l * = max({l i }). The reason why we pick this strongest decoherence as the main scaling factor is that only then α n depends on asymptotically weakening or constant decoherence processes and those processes that get weaker with increasing N are not crucial to span the Hamiltonian (as all L (i) are sufficient to span the Hamiltonian on their own).
Consider now, a more complicated situation where in order to satisfy the HLS n condition we need to use noise operators with two (or more) different types of non-linearities {l i } simultaneously, since none of them on their own is sufficient to setβ n = 0. In this case, the coefficients h (and h) which multiply operators with a given non-linearities l i will scale as N k−li/2−n/2 (N k−li ). When inspecting the asymptotic behavior ofα n we then note that the dominant scaling ofα n with N will be determined by the coefficients h (and h) which have the highest scaling exponent-the terms corresponding to noise operators with the lowest non-linearity level. Therefore, in the scaling formula (15), we need to take l * = min({l i }). The presence of higher order nonlinearities, may influence the bound as well, but only in via modification of the multiplying constant and not the scaling-we will see this behavior in IV when discussing the impact of single and two-body losses on fundamental precision bounds in cold atom metrology.
The argument that bases on the relative strength of decoherence processes does not work in this case. Even though the relative strength grows with N it does not dominate the norm ofα n . The part with lower non-linearity is crucial to span the Hamiltonian and so no matter how much the strength of decoherences with higher non-linearities grow, the norm will depend mostly on presence of this lower degree.
So, in short, given Hamiltonian with a fixed nonlinearity k we can say that l * in bound (15) is equal to the highest level of non-linearity that when combined with higher non-linearity terms allows to satisfy the HLS n condition: Finally, let us allow the Hamiltonian to consist of many terms with different non-linearity levels k. In this case for each k-body Hamiltonian part there will be a corresponding l * (k) providing the tightest bound as given by (16). Since all the Hamiltonian parts are present simultaneously, the resulting scaling will be determined by the part that dominates the formula forα n . This will be the term with nonlinearity k for which 2k − l * (k) is maximal. Hence, in general, the scaling of precision resulting from the RPN approach is given by the formula (15) with This bound resolves the open problem of scaling of metrological bounds analyzed in [39,40].

B. Adapting the bound to experimental realities
While the preceding considerations take into account experimental imperfections (general Markovian noise), there are some implicit assumptions that were made on the way to the formula for the fundamental bound (4) that might be modified in order to take into account additional experimental limitations present in the physical systems considered.
First of all, throughout the derivation we have assumed a fixed number of particles N being present in the system. It does not mean that we exclude losses (they can be formally incorporated via additional internal degrees of freedom of the particles), but it means that in principle, the adaptive strategies are allowed to reload the particles continuously during the course of the experiment. In real cold atom experiments, this is not the case, and only after the experiment with a given cloud of atoms is finished a new sample is prepared. If we want our bound to represent this experimental limitation we may do it by taking into account the time dependence of the number of atoms in the experiment N (t). Given this function, we may easily adapt the reasoning leading to the bound (4) following the general recipe discussed for time dependent models in [16], and obtain the corresponding bound where the dependence ofα n ,β n on t is due to the rescaling procedure that leads from non-tilded quantities to tilded quantities involving the particle number parameter N , which now is time dependent.
The other potential modification in the preceding derivation may be due to the fact that the noise parameters have some additional nontrivial dependence on the particle number on their own. As discusses in detail in Sec. IV this is a typical situation in cold atom systems, where with increasing number of particles the size of the atomic cloud increases and effectively the two-body loss coefficient decreases. In such situations, this effect needs to be taken into account in order to provide the proper scaling exponent. Still provided that noise coefficient go down slower than N −1 the hierarchy of non-linearities remains valid, in the sense that higher order non-linear terms dominate over lower order non-linear terms in the limit of large N -see Sec. IV for more detailed discussion in the context of cold atom interferometry experiments.

IV. QUANTUM INTERFEROMETRY WITH SINGLE-AND TWO-BODY LOSSES
In this section we present the full analysis of the simultaneous effects of single-and two-body losses on quantum interferometric protocols, keeping the cold atom physical context in mind as the main motivation for this study. Consider an atomic interferometry model with a system of N two-level atoms where the dynamics is described by the following master equation: where is the standard linear interferometry Hamiltonian, with a i (i = 1, 2) being the annihilation operators removing an atom from level i, while the noise operators corresponding to the single and two-body losses parts L (1) and L (2) read respectively In [16] it was stated that since two-body loss operators and their products are linearly independent from H (1) , the system may be effectively protected from two-body losses via some kind of quantum error-correction protocol, and the resulting bound on precision will be determined solely by single-body losses in the case of the linear Hamiltonian. This reasoning made an implicit assumption, however, that we put no additional restriction on the allowed class of states. We show below that if we impose a superselection rule that forbids preparation of superpositions of states with different total atom numbers, and if at least two out of the three two-body loss coefficients γ ij are nonzero, it is no longer possible to protect the system against two-body losses. As a result, in the large particle limit the actual bound will be determined by the two-body effects which leads to a much less favorable scaling of precision than the one resulting from the impact of single-body losses only.
A. Impossibility to protect against two-body losses under atom number superselection rule Following the general considerations presented in [16,17], in order to protect the sensing system from two-body losses we need to construct at least a twodimensional code space, {|ψ 1 , |ψ 2 }, where the following conditions need to be satisfied: where the bold index j ∈ {(1, 1), (1, 2), (2, 2), 0} represents the two-body loss operator double indices, where value of the index 0 returns the identity matrix L (2) 0 = 1 1, and µ jj is some Hermitian matrix. Furthermore, in order to have a nontrivial effective action of the Hamiltonian in the code space, and hence be able to sense the parameter of interest, we require that Consider now We now assume that the code states have a definite atom number, and hence (a † 1 a 1 + a † 2 a 2 )|ψ k = N |ψ k -note that the atoms may still be entangled with some additional ancillas, and hence |ψ k may in principle live on a larger Hilbert space, but the assumption simply means that within the system Hilbert space we will only make use of states with a fixed number of particles in both modes. Using this assumption we get Provided γ 12 > 0 the error-correction condition for the L 12 operator yields ψ k |a † 1 a 1 a † 2 a 2 |ψ k = 1 γ 12 δ kk µ 12,12 .
As a result we have: and analogously for ψ k |a † 2 a 2 |ψ k . This implies that and hence the effective Hamiltonian is trivial on the codespace. The above derivation is valid provided all γ ij > 0. However, since (a † 1 a 1 + a † 2 a 2 )|ψ k = N |ψ k even if one of the coefficients is zero two out of three error correction (25,26,27) conditions are sufficient to arrive at the same conclusion.
If on the other hand two coefficients are zero the reasoning fails and this is the only situation when one can construct an error correcting code to protect the sensing process from two-body losses with the atom number superselection rule imposed.
To give a concrete example. Assume γ 11 > 0 while γ 12 = 0, γ 22 = 0. In this case a possible codespace (which may be shown to be optimal, using the criterion given in Eq. (36) in [17]) consists of vectors , and we assumed N > 4. The only nontrivial error correcting condition to check is ψ 1 |a †2 1 a 2 1 |ψ 1 = ψ 2 |a †2 1 a 2 1 |ψ 2 which indeed holds thanks to the choice of the parameter s as above. Furthermore, the Hamiltonian is nontrivial in the codespace, since ψ 1 |H (1) |ψ 1 = − N 2 4(N −1) , while all other matrix elements of H (1) in this subspace are zero. In the end this leads to QFI F Q = T 2 N 4 (N −1) 2 16 , which yields the scaling as in the ideal lossless case but with a quantitative reduction by approximately a factor of 16. Similar constructions can be performed when a different loss coefficient is nonzero.
As mentioned above, if two or more two-body loss coefficients are non-zero it is not possible to construct the correcting protocol with the superselection rule imposed. If one insists on constructing such a code one must violate the particle number superselection rule. For example, a code that works in case all types of two-body loss processes are present, can be written as a simple generalization of the code presented above: |ψ 1 = (s|N + √ 1 − s 2 |0 )|N/2 , |ψ 2 = |N/2 (s|N + √ 1 − s 2 |0 ). It can be easily checked that, provided N > 4, it satisfies errorcorrecting conditions with respect to all types of twobody loss processes and yields a non-trivial Hamiltonian within the codespace.

B. Asymptotic bounds due to two-body losses
Assuming that at least two γ ij coefficients are nonzero, we will now use the formula (4) to derive the bound on performance of linear interferometry due to the presence of two-body losses. Assuming that we are in the large N limit, and γ ij coefficient are constant with growing N (or declining slower than N −1 ) this will also be a valid bound in presence of both single-and two-body losses, since as discussed in detail in Sec. II the two-body effects will dominate over single-body ones in the limit of large N .
The procedure will be similar to that presented in [16] where the bounds due to two-body losses have been derived, but with a notable exception that here we focus on the linear Hamiltonian rather than the quadratic one. where we neglected the terms for which the coefficients are set to zero as they do not help in satisfying the condition β = 0. We keep in mind, that the operator will be acting within the space of states with total fixed atom number N , hence we can use the relation a † 1 a 1 + a † 2 a 2 = N 1 1 provided these operators will act directly on the potential states (in other words are not acted upon by other operators in the formula). As a result we can write: Consequently, we can express all coefficients as a function of a single parameter h: h 11,11 = − (1/2 + ξ) /(N γ 11 ), h 22,22 = (1/2 − ξ) /(N γ 22 ), h 12,12 = −2ξ/(N γ 12 ), where ξ = h N , and we approximated N − 1 ≈ N as we are focused on the large N regime. In order to obtain the bound we now need to minimize the operator norm of over the parameter ξ. All involved operators are diagonal in the atom number basis |n |N − n . Hence the operator norm of α equals the largest absolute value of its diagonal elements in this basis and as a result the formula for the bound reads (again approximating N − 1 ≈ N when necessary): where we will treat x = n/N as a continuous parameter in what follows. This optimization can be easily done numerically, as this function involves at most products of quadratic terms in x and ξ. We can also extract an analytical formula from this expression with some additional assumptions on the relation between γ ij coefficients.
With fixed ξ the bound is a quadratic function of x. Hence the maximum over x is achieved on either the boundaries of the region x = 0, x = 1 or at the extremum point. If for a given ξ the function x is convex in x then the maximum is indeed achieved on the boundaries. This is the case provided 4ξ 2 /γ 12 ≥ (ξ + 1/2) 2 /γ 11 + (1/2 − ξ) 2 /γ 22 . Assuming this condition is satisfied we can write: where the minimum is achieved for ξ = 1 √ γ11 . The above reasoning is valid provided the function is indeed convex in x for the optimal ξ parameter found above. Plugging it into the condition for convexity we arrive at the condition for validity of the above bound γ 12 ≥ 1 2 √ γ 11 − √ γ 22 2 , and the final form the bound reads: (35) In particular, the above abound is always valid in the case of symmetric losses γ 11 = γ 22 . Note that the bound is similar to the one derived for singlebody losses [16,21], but differs by the missing factor N , which is in agreement with our general scaling considerations which imply that for k = 1 Hamiltonian and l = 2 noise-type QFI should approach a constant value for large N .
The other case that is easy to do analytically and is also relevant in experimental context (see Sec. IV C) is the case where one of the two-body loss coefficients is zero. Let us take γ 22 = 0 and γ 12 , γ 11 > 0. In this case the minimization over ξ becomes trivial because γ 22 = 0 implies ξ = −1/2, note the discussion below Eq. (31). We are left with the task of calculating max 0≤x≤1 γ12 . This quadratic function decreases monotonously on 0 ≤ x ≤ 1 for γ 11 ≤ 2γ 12 , hence the maximum is attained for x = 0. In the case when γ 11 > 2γ 12 the function is concave with the extremum point in x = 2γ12−γ11 2(γ12−γ11) . Finally the bound reads C. The reduced particle number (RPN) approach We will now show how the problem of deriving the bound-taking into account single and two-body losses-can be analyzed within the RPN approach, and discuss the advantages compared with the reasoning presented in the previous subsection. If we would like to derive a bound in the regime of N where contribution of single and two-body is comparable we would not be able to follow the simplified reasoning presented in the previous subsection. Instead we would need to take into account all single and two-body loss operators, consider the condition for β = 0 and minimize the operator norm of α. This would be cumbersome, first of all due to a larger number of independent parameters over which we need to optimize and furthermore when calculating the operator norm we need to face the fact that the operators act in a very large Hilbert space. The RPN approach will allow us to reduce the problem to a small dimensional Hilbert space, which allows us to effectively perform numerical optimization required for the calculation of the bound even if the procedure involves optimization over many free parameters.
We follow here the RPN approach using twoparticle n = 2 sub-channels. Apart from the two internal states |1 , |2 representing the particle being in mode 1 or 2 respectively, we also introduce a vacuum state to represent the particle being lost |v , and hence the two-particle space is 9 dimensional. Within the space of two particles we will also restrict ourselves to symmetric states, as this is the class of states that we encounter in cold atom systems, and moreover, we have numerically checked that relaxing this requirement did not have any impact on the derived bounds. The symmetric subspace is 6 dimensional and is spanned by H S 2 := {|0, 0 , |0, 1 , |1, 0 , |0, 2 , |1, 1 , |2, 0 }, where we use the standard occupation number notation with numbers representing the number of particles occupying mode 1 and 2 respectively. The relation between these basis states and the states of distinguishable particles with three internal degrees is straightforward. For example: |0, 0 = |v ⊗ |v , We can now write the appropriately rescaled Hamiltonian and noise operators where the annihilation operators restricted to the considered subspace are given by: Since all operators now act in a 6-dimensional space one can calculate the bounds efficiently using the semi-definite program defined in (14).
Without resorting to numerics, using simple algebra of 6 × 6 matrices one can also analyze the implications ofβ 2 = 0 condition. First, of all one notices that whileH 2 ∈ S(L (1) ) is inside the space spanned by single-body loss operators, it is not spanned by the the two-body loss operatorsH 2 / ∈ S(L (2) ). This is analogous to the situation we have encountered when analyzing the full operators acting on the whole Hilbert space without the atom-number superselection rule imposed. Indeed, this is to be expected, since in the way we approached the problem using the RPN formalism we in fact allow for states which are superposition of different atom number states up to the total atom number n = 2. If we want to derive tighter bounds using the RPN approach, that take into account the atom number superselection rule we must formulate the problem mathematically in a way that forbids the use of superposition of different atom number states.
The simplest approach, is to modify the operators in Eqs. (37) in a way that their input space is restricted to the subspace with exactly two-atoms, so H S 2! = {|0, 2 , |1, 1 , |2, 0 }. As discussed earlier, physically this would correspond to the situation where, at every adaptive step of the protocol lost atoms are being replaced with new ones, so that the total number of atoms is unchanged and we can consider only states with definite number of atoms throughout the protocol. Formally, this assumption can be incorporated in the RPN approach by restricting the input space of noise operatorsL 2 P H S 2! , where P H S 2! is the projection on H S 2! , as well as projecting the Hamiltonian on this subspaceH Note, that we are multiplying noise operators only from one side, as this is how they act on the state of the system, while we project the Hamiltonian from both sides. By doing so, we get thatH (1) 2! is effectively a 3 × 3 matrix acting on H S 2! and so are the products of noise operators and their Hermitian conjugationL † 2!,jL 2!,i (the single noise operators that also appear in the definition of the S(L 2! ) space will not be relevant as they introduce terms that are beyond the space spanned by the Hamiltonian and hence need to be set to zero in order to satisfyβ 2! = 0).
The operators restricted to the H S 2! subspace which are relevant for derivation of the bound read explicitly: It is now clear, that provided two of the two-body loss coefficients are non-zeroH 2! ) and hence we recover the observation that we have arrived at while considering the operators on the whole N -particle Hilbert space: with the super-selection rule imposed the two-body losses will determine the asymptotic scaling in the limit of large N . Here we see it since the two-body noise terms are not renormalized by 1/(N − 1) factor and hence will be dominating in the large N limit.
In the case considered in Eq. (35) (so in articular symmetric losses) the bound achieved in the RPN model is the same no matter what n we take. However, if we consider an experimentally motivated case γ 22 = 0, γ 12 > 0, and γ 11 > 0, we see that we can tighten the bound in some regime of parameters by increasing n, see Fig. 3. The bigger space gives us freedom to consider states with the different ratio of atom occupation numbers in modes 1 and 2, which is useful to tighten the bound in some parameter regions. So in this case, it appears that while, qualitative scaling of the bound is seen in the RPN approach already at the smallest n = 2 subchannel considerations, the tightness of the bound may sometimes be improved by increasing n.

D. Comparison of the bounds with simulations and experimental results
The physical system particularly relevant to theoretical considerations is the Bose-Einstein conden- sate. Here we briefly recall basics of the bimodal Bose-Einstein condensates in the context of metrology. It will give us better understanding of the experimental situation necessary to make a meaningful comparison of the presented bounds with experimental results.
We focus on the case with |1 and |2 corresponding to internal atomic states [35,41], although there are also experiments on different spatial modes of a double well potential [42][43][44] and proposals to use different orbitals of a scalar fragmented Bose-Einstein condensate [45]. It is instructive to consider first two clouds of atoms, one with N 1 atoms in an internal state |1 and the second with N 2 atoms in another internal state |2 . The Bose-Einstein condensation means that the total state of the system is close to a pure separable state, built up from two macroscopically occupied orbitals, ψ The formula relies on the fact, that the atoms in the internal states |1 are distinguishable from the atoms in the state |2 . The symmetrization has to be performed between atoms in |1 , and between atoms in |2 , but not necessarily between atoms in different internal states. In practice, when the evolution starts with a single Bose-Einstein condensate with all atoms initially in |1 , shined with the Rabi pulse to obtain the spin coherent state, the state is totally symmetric. These subtle differences, do not impact the further considerations. It follows from the theory of condensation, that the two orbitals are minimal energy solutions of coupled Gross-Pitaevskii equations [46,47] where coupling constants g refer to the interaction between two atoms, one of which is in the state | , and the other in | and , m is the atomic mass of the atoms. The eigenvalues µ are Lagrange multipliers, which physically are the chemical potentials of condensates. Finally, functions U (r) are the statedependent potentials trapping atoms where ω σ are the trap frequencies, and r ( ) σ are trap centers for the atoms in the state | .
The energy of the two condensates is, up to a constant, given by the formula In a typical experiment the initial state is a spin coherent state which is a superposition with the total number of atoms differently distributed between two condensates. The evolution of such system (in the states for which dispersions of occupations are small, precisely ∆N 1 N and ∆N 2 N ) is approximately generated by the operator whereŜ z := N 1 −N 2 /2 is the imbalance between condensates' occupations. Parameters of the Hamiltonian (47) have to be derived from the solutions of the Gross-Pitaevskii equations and the corresponding energies according to [48]: The derivatives with respect to occupations N 1 and N 2 are evaluated around the average occupations, i.e.N 1 andN 2 , which in the initial state are equal to N/2. The central term of the Hamiltonian (47) is the one-axis twisting model, χŜ 2 z . Hence, one expects that the system can be used to obtain entangled states, in particular states useful in metrology. Indeed, the spin-squeezed state of two Bose-Einstein condensates have been produced already a decade ago [35,41,42] and used to demonstrate a gain in precision of an interferometer.
On the other hand, the potential benefits from condensation may be limited due to particle losses, which are an important source of decoherence for Bose-Einstein condensates. Particle losses are usually divided into three classes, respectively of the number of atoms which are lost in a single event.
One-body losses are caused by collisions between condensed atoms and other particles, from the residual gas unpumped from the vacuum chamber in which the condensate is kept.
The term χŜ 2 z in the Hamiltonian (47), which is responsible for generation of entanglement, comes in fact from elastic collisions between atoms. It happens, however, with little probability, that the collision is not elastic -in consequence of such inelastic collision both atoms are lost from the trap. The rate of two-body losses depends on the rate of the binary collisions, which stems from both Bose-Einstein wave-functions. Atomic constants K are associated with probabilities that the collision will be inelastic. In the limit of large number of atoms the condensates enter the so called Thomas-Fermi regime, in which one finds excellent analytical approximations for the solutions of the Gross-Pitaevskii equations. Using this Thomas-Fermi approximations one may evaluate Eq. (50) to find analytical approximation for two-body loss rates. In case of atoms in the spherically symmetric harmonic trapping potential, assuming equal coupling strengths g 11 = g 22 = g 12 =: g one gets where a = g m 4 2 π is called scattering length, l osc = mω is the oscillatory length and N is the initial number of atoms. Due to the factor N −3/5 the rate of two-body losses is getting lower when the number of atoms increases. This is caused by the van der Waals repulsion between condensed atoms, which results in the spatial broadening of the condensate wave function. Finally, the integrals GPE | 2 decays, so as the twobody rates as given in Eq. (50). This has far going consequences for scaling of the optimal Quantum Fisher Information F Q . One has to remember that the total number of lost atoms per unit of time increases with the total number of atoms, but with slower rate compared to the case with constant γ. In Fig. 4 the approximation (51) is compared to the exact result, obtained directly from Eq. (50), but using the exact numerical solutions of equations (43) for the orbitals ψ  [49,50] in which the potential trapping atoms is practically equal to a box. In this case the condensate wave-functions are practically constants ψ GPE = 1/ √ V , where V is the volume of the trap. Then Eq. (50) leads to rate of two-body losses independent on the number of particles. The influence of the two cases-the Thomas-Fermi approximation and the box trap-on bounds on QFI is depicted in Fig. 5. This figure also depicts the main result of this paper concerning scaling of the fundamental bounds with the number of particles. Constant decoherence parameters result in constant precision bound, as we have already showed in Eqs. (35), (36). If on the other hand γ scales with N the fundamental bound changes significantly: to O(N 3/5 ), which shows that manipulation of density of the atomic cloud offers trade-offs visible also in fundamental bounds.
Another experimental complication is the fluctuation of the total number of atoms in the initial state. Usually N is a stochastic variable, with distribution close to the Poissonian one. Due to the term χNŜ z , the fluctuations of total N leads to a phase noise. This effect is known under the name "collisional shift".
As it is impossible to fully avoid particle losses, hence the entangled state has to be produced and used on faster time-scale than the ones connected with particle losses. This demands large value of the non-linearity parameter χ. From the equation for χ (48) and energy (46) one finds the approximation In the case of the widely used 87 Rb atoms, the interaction strength parameters g are close to each other. In consequence, if atoms in |1 and |2 share a common spatial mode, i.e. ψ GPE , then the nonlinearity practically vanishes χ ≈ 0. One of the way to increase χ is to separate both clouds using the state dependent trapping potentials [41] and varying their centers r ( ) σ , the other -to change one of the interaction strengths using Feshbach resonances [35]. The entangling dynamics is initiated by taking the two atomic clouds apart to maximize non-linearity. In this preparation stage the useful state, as a spin squeezed state can be obtained. Due to asymmetry in the loss rates, the state leaves the equator as shown schematically in the pseudo Bloch sphere. To use the state in linear interferometry, first the two clouds have to be merged into a common spatial mode, and then the quantum state encoded in the internal degrees of freedom is appropriately rotated. Afterwards the state should accumulate a phase.
In both cases, once the target entangled state is prepared at T PREP during the nonlinear dynamics, the further entangling evolution has to be stopped. In the case of two separated condensates it is sufficient to bring both clouds together, ideally to a common spatial mode. The common spatial mode is also necessary to perform any SU(2) rotations, usually required to adjust the state appropriately to the quantum information task. In particular, Ramsey interferometry requires appropriate orientation of the state with respect to the phase-imprinting Hamiltonian. If phase will be accumulated in the dynamics generated by the HamiltonianĤ ∝Ŝ z , as discussed here, and the entangled state is squeezed, then the state on the Bloch sphere has to be elongated along the meridian, see Fig. 6.
We were simulating numerically the scheme presented in Fig. 6. For each chosen Ramsey time T , we were optimizing the following Fisher Information with respect to the preparation time T PREP F s (T ) = max The state of the system at time t was computed from the master equations using quantum trajectory method [51]. The parameters of the master equations were computed numerically, by solving the coupled Gross-Pitaevskii equations (43), then evaluating Eqns. (46), (48) to find non-linearities χ andχ and (50) to find the rate of two-body losses. Three body losses were also included, but they do not played any important role. The parameters were computed separately for the nonlinear evolution, during which the two condensates are taken  (53), "shallow" denotes the simulation with a modulated trap that decreases the impact of two-body losses. The green points are the result of the experiment [7]. In our numerical simulations we assumed the trap frequencies as in the experiment [7], equal to (νx, νy, νz) = (540, 540, 115)Hz with in average N = 1400 atoms initially. The loss rates were computed from Gross-Pitaevskii equation applied to two clouds of N1 = N2 = 700 atoms, specified in the caption of Fig. 4.
apart, and separately for the interferometry , during which the two clouds should stay together (to avoid further nonlinear dynamics). The initial number of atoms was stochastic with the Poissonian distribution to refer to the experimental situation. As long as no special tricks are used, like the compensation method [52,53], then the optimal preparation times are short, resulting in the weakly squeezed states. In this case there is no substantial difference between the full Quantum Fisher Information and its estimate from the error propagation formula based on the measurements of the spin components.
In Fig. 7 we present a comparison of precision of estimation coming from fundamental bounds, numerical simulations, and an actual experiment. The two curves representing fundamental bounds show that in the case of the experimental loss rates of [7] the influence of two-body losses on the best achievable precision is marginal. This small discrepancy between fundamental bounds is mainly caused by the small number of particles. Note that in Fig. 5 for N = 1400 the bound in the Thomas-Fermi approximation scales as QFI for single-body losses. There is a subtle difference, between the experiment and theory: the experimental frequency estimation at time T , was not based on stopping the experiment at an intermediate time t and then repeating the experiment T /t times, as it is assumed in Eq. (53). The experimentally found ∆ω should be rather benchmarked with 1/ F Q (T ), where F Q (T ) comes directly from the definition (2). As in the following expression: The deviations due to the subtleties in the definition of the uncertainty bound are definitely not sufficient to explain the huge difference between the experimental results and the ultimate bound.
We investigate the origins of the theory-toexperiment difference using the numerical simulations described above. The brown curve in Fig. 7 shows the numerical results for the trap geometry, the separation between clouds, and the fixed preparation time mimicking the experimental situation [7]. On the other hand, the dynamical spatial phenomena and the phase noise coming from the thermal effects are omitted in the simulation. Especially the former factor may be important-in the preparation stage of the very experiment [7] the clouds' centers of masses were not stationary, but oscillating with amplitude around 140nm. Although the simulation is not comprehensive, still it stays relatively close to the experiment. Then we repeated the simulation but fully separating both clouds (which is still reasonable, compare with [54]), optimizing over the preparation time and using the definition (53). Due to larger separation the non-linearity parameter χ increases, whereas losses slightly decrease. By this, the initial states for Ramsey interferometry are less deteriorated by decoherence within the preparation stage. Precision increases by order of magnitude. Especially important is the optimization over the preparation time-by adjusting T PREP separately to each Ramsey time T we have found that experimentally strongly squeezed states are in fact a bad choice for long interferometry durations. Finally we considered the next scenario, with the atomic clouds placed for the Ramsay interferometry into a very shallow trap. In this case, there are no two-body losses, neither the residual non-linearities effect, i.e.χ = 0, within interferometry. This last simulation, shows how much one can gain by avoiding the two-body losses and the collisional shift in the real experiment.

V. CONCLUSIONS
We have shown how to use the recently developed tools of theoretical quantum metrology to analyze in great detail interferometric experiments performed with cold atoms. We provide the fundamental and more realistic bounds on precision and complete the ongoing discussion on achievable scaling of QFI with the number of particles. Future work might involve adopting presented methods in multi-parameter metrology. On the other hand, as the field of experimental non-linear metrolgy grows rapidly, it will be necessary to perform similar analysis for setups involving the quadratic Hamiltonian. We analyzed theoretical bounds in the context of interferometry based on Bose-Einstein condensates. We benchmarked our theoretical results with experimental data [7] and numerical simulations. First, we have shown that our simulations are close to the experimental results presented in [7]. Then, we focused on ways to improve experimental precision by tuning "simple" parameters, like separation between atomic clouds in the preparation stage, duration of the entangling dynamics and frequency of the final trap. One could in principle increase the precision by an order of magnitude. For longer Ramsey times, of the order of 100ms, it is actually beneficial to stop the "state preparation" stage early, much before reaching the maximally squeezed states. Numerical analysis should be extended by a better model of the spatial dynamics and the non-adiabatic effects -the factors influencing strongly the quantum engineering using the Bose-Einstein condensates [55]. As the benchmark between simulations for the optimized parameter and fundamental bounds shows a big room for improvement, one should perform full optimization of the tunable parameters to either direct the experiments on the fastest track or to find the most important limiting factors.

VI. ACKNOWLEDGMENTS
We thank Philipp Treutlein for giving us access to cold atom magnetometry experimental data. JC acknowledges support by the Netherlands Organi-