Karim Douch
email: douch@ife.uni-hannover.de

Hu Wu
email: wuhu@ife.uni-hannover.de

Christian Schubert
email: schubert@iqo.uni-hannover.de

Jürgen Müller
email: mueller@ife.uni-hannover.de

Franck Pereira
email: franck.pereira@obspm.fr

Dos Santos

Simulation-based evaluation of a cold atom interferometry gradiometer concept for gravity field recovery

Keywords: Cold atom interferometry, space gravity gradiometry, gravity field recovery, closed-loop simulation

Introduction

The GOCE [START_REF] Rummel | GOCE gravitational gradiometry[END_REF] and GRACE [START_REF] Tapley | The gravity recovery and climate experiment: Mission overview and early results[END_REF]) missions have provided unprecedented insights in the static and time-variable gravitational field of the Earth and proved to be very useful in numerous fields, from unification of height systems to the determination of global mass distribution and mass transport in the Earth system [START_REF] Kusche | Mass distribution and mass transport in the Earth system[END_REF].

The two satellites of the gravimetry mission GRACE (Gravity Recovery and Climate Experiment) are orbiting the Earth since 2002. Its major measurement method relies on Satellite-to-Satellite ranging which measures the varying distance between the two satellites with µm precision using a microwave tracking system. GRACE provided monthly solutions of the gravity field during most of its lifetime and therefore enables analysing its temporal variations. The latter are caused, among others, by ice mass loss in Greenland or Antarctica, glacial isostatic rebound or changes in the global hydrological cycle see, e.g. [START_REF] Eicker | Does GRACE see the terrestrial water cycle intensifying?[END_REF]. The GOCE (Gravity field and steadystate Ocean Circulation Explorer) mission delivered data between 2009 and 2013, which allowed to determine the geoid [START_REF] Torge | Geodesy[END_REF] with an accuracy of about 1 to 2 cm for a spatial resolution of 100 km on the Earth surface. This provided the basis to establish global physical height systems (i.e. those related to the gravity field) and -combined with altimetric measurements of the mean sea surface -to estimate global ocean currents. GOCE is also the first mission which employed gradiometry. The core of the GOCE gradiometer consists of 3 pairs of electrostatic accelerometers measuring the differential gravitational accelerations on its sensitive axis with a noise of about 10 mE/ √ Hz (1 E = 10 -9 s -2 ) in a measurement bandwidth (MBW) between 5 and 100 Hz [START_REF] Rummel | GOCE gravitational gradiometry[END_REF]. Hereafter, the noise level will be expressed (if not mentioned otherwise) in terms of amplitude spectral density (ASD) which is the square root of the power spectral density.

Novel science results were obtained thanks to GRACE and GOCE, giving insight in a multitude of processes in the hydrosphere, cryosphere, oceans, atmosphere and solid Earth that could not be assessed before. As stressed in various reports as for instance in [START_REF] Pail | Science and user needs for observing global mass transport to understand global change and to benefit society[END_REF] or [START_REF] Panet | Earth system mass transport mission (e.motion): A concept for future Earth gravity field measurements from space[END_REF], there is a strong demand for sustaining a global, homogeneous and continuous observation of the Earth's gravitational field. Furthermore, the spatial and temporal resolution along with the precision should be improved in order to meet the challenging requirements of the different end users. Therefore, new concepts for follow-on missions have been developed with the aim to achieve a better resolution (ideally below 100 km spatially and below 1 month temporally), a higher accuracy and a longer time series [START_REF] Gruber | e2. motion: Earth system mass transport mission (square)concept for a next generation gravity field missionfinal report of project satellite gravimetry of the next generation (nggm-d)[END_REF][START_REF] Elsaka | Comparing seven candidate mission configurations for temporal gravity field retrieval through full-scale numerical simulation[END_REF].

In this paper we investigate to what extent a new concept of gradiometer based on cold-atom interferometry (CAI) sensors [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF] can reach such objectives. Gradiometry offers several advantages: a gradiometer measures directly a functional of the gravitational field, moreover, it performs like a high-pass filter increasing the relative sensitivity to the high-frequencies of the gravitational signal corresponding to the small wavelengths of the field. The constraints for the measurement bandwidth of a CAI gradiometer might be relaxed and potentially a higher accuracy might be achieved. Last but not least, CAI technology provides a very good stability of the gradiometer bias and scale factor, a high common mode rejection of 140 dB was demonstrated by [START_REF] Mcguirk | Sensitive absolutegravity gradiometry using atom interferometry[END_REF] and is likely to have white noise over a wide frequency range while the GOCE gradiometer noise was approximately white only in its MBW.

After a brief description of the physical principle of the gradiometer concept and the conditions to operate it in space, the method of evaluation of its performance in terms of quality of the recovered gravity field model is described. We have chosen to carry out the study through a closed-loop simulation of the measurement process in the time domain. Finally the results and the technical feasibility of the proposed concept are discussed.

Principle of the CAI gradiometer concept

The gravity gradients V ij are the second-order partial derivatives of the gravitational potential V :

V ij = ∂ 2 V ∂i∂j . (1)
where x i and x j , are Cartesian coordinates of the 3 dimensional space. The resulting 9 components give the GGT (gravitational gradient tensor). In the following, we note V the GGT and V ij i, j ∈ {x, y, z} any component of the tensor. Since the gradiometer is an inertial sensor, it is not possible in a non-inertial frame to measure directly the gravitational gradients, but rather the gradients of the acceleration field. The latter are however connected to the gravitational gradients as follows (for the diagonal elements):

V xx = Γ xx -ω 2 z -ω 2 y V yy = Γ yy -ω 2 z -ω 2 x V zz = Γ zz -ω 2 x -ω 2 y (2)
where Γ ii i ∈ {x, y, z} is the diagonal acceleration gradient along the i-axis and (ω x , ω y , ω z ) t = Ω is the vector of the angular velocity of the GRF (gradiometer reference frame) with respect to the IRF (inertial reference frame). We will refer hereafter to the terms in ω in equations (2) as the centrifugal terms. Finally, two other observables must be given in order to define correctly the gravitational gradient: the position of the measurement point and the orientation of the measurement frame, i.e. the attitude of the GRF with respect to the IRF.

In the case of the CAI gradiometer concept presented here, the acceleration gradient is approximated by the difference of acceleration measured simultaneously at adjacent points. For the x axis, one has for instance

Γ xx = ∆a x l (3)
where ∆a x is the difference of the x-component of the acceleration between two points aligned along the x-direction and separated by a baseline l. Similar expressions hold for the two other diagonal acceleration gradients.

Here we introduce the light-pulse atom interferometry utilized in the acceleration measurement and the CAI gradiometer considered in this study (1) The laser beam features two frequencies ω 1 , ω 2 and is retro reflected at the mirror. This leads to four different components in the light, distinguished by frequency and / or propagation direction (k). Depending on the duration of the interaction of the light field with the atoms, the atoms are coherently excited into a superposition state of initial momentum p 0 and p 0 + (k 1 -k 2 ) ( reduced Planck's constant) by absorption and stimulated emission of a photon from the incoming (wave vector k 1 ) and the retro reflected beam (wave vector k 2 ). During this process, the phase k • d is imprinted onto the atoms, depending on the distance of the atoms to the mirror d and the effective wave vector k = k 1 -k 2 . (2) In the special case of negligible velocity of the atom parallel to the effective wave vector k, a four photon process leads to a symmetric momentum transfer with the superposition of the three momentum states p 0 and p 0 ± (k 1 -k 2 ), effectively doubling the effective wave number.

Figure 2: Atom interferometry scheme for gradiometric measurements as proposed in [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF]. During the preparation time t prep an atomic ensemble is generated, split into two leading to a separation l between the ensembles, and injected into the two atom interferometers with forward drift velocity v. Here, the atomic ensembles are in free fall. The three-pulse sequence with drift times T forms a diamond shaped geometry. The first pulse imprints a reference phase (k • d) onto the atoms to which the phases imprinted at the two subsequent pulses is compared. In a simplified picture, the total phase accumulates to kd(0) -2kd(T ) + kd(2T ), where d(t) denotes the centre of wave packet position (grey dotted line) and other spurious contributions are omitted. With

d(t) = d 0 + v(t) + 1
2 aT 2 the phase shift becomes k • aT 2 . This result is recovered for the upper interferometer with the replacement d → d + l. Differential displacements (d 0 = d 1 = d 2 = d 3 ) or vibration of the mirrors (1 -3) are suppressed in the differential signal. Depending on the total phase shift, the transition probability oscillates between the inner an the outer output port. During detection, this transition probability is recorded by flashing light onto the atoms and either collecting the fluorescence or measuring the absorption. This configuration enables a measurement in 1D and can be reproduced for the other two axis for a 3D measurement. A preparation time t prep shorter than the total interferometer time 2T enables an interleaved mode of operation [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF].

which was proposed in [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF], based on recent developments in microgravity experiments (Barrett et al., 2016;[START_REF] Rudolph | A high-flux BEC source for mobile atom interferometers[END_REF][START_REF] Aguilera | Ste-questtest of the universality of free fall using cold atom interferometry[END_REF][START_REF] Müntinga | Interferometry with Bose-Einstein Condensates in Microgravity[END_REF]. An atom interferometer exploits interference of matter waves which are coherently manipulated by light fields. At each atom-light interaction, a phase is imprinted onto the atoms which depends on the relative position between the light field and the atoms (see 1). The light field, consisting of two frequencies, is retro reflected at a mirror and drives two-photon transitions during each interaction. This couples different momentum states which are required to form the interferometer geometry sensitive to accelerations. Effectively, the mirror serves as a reference plane for the position of the atoms during the interactions. A typical implementation for an atom interferometer sensitive to accelerations is a three pulse sequence as used in various experiments [START_REF] Peters | Measurement of gravitational acceleration by dropping atoms[END_REF][START_REF] Mcguirk | Sensitive absolutegravity gradiometry using atom interferometry[END_REF][START_REF] Louchet-Chauvet | The influence of transverse motion within an atomic gravimeter[END_REF][START_REF] Hu | Demonstration of an ultrahigh-sensitivity atom-interferometry absolute gravimeter[END_REF][START_REF] Rosi | Precision measurement of the newtonian gravitational constant using cold atoms[END_REF]Kovachy et al., 2015a;[START_REF] Berg | Composite-light-pulse technique for high-precision atom interferometry[END_REF][START_REF] Freier | Mobile quantum gravity sensor with unprecedented stability[END_REF][START_REF] Abend | Atom-chip fountain gravimeter[END_REF] and also proposed for the CAI gradiometer [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF] as shown in fig. 2. Initially, a beam splitter pulse creates a coherent superposition between two momentum states. Due to the double diffraction [START_REF] Lévèque | Enhancing the Area of a Raman Atom Interferometer Using a Versatile Double-Diffraction Technique[END_REF] employed here, the trajectories drift apart with four single photon momenta. After a free evolution time T , a second pulse coherently inverts the momenta. This allows closing the atom interferometer with a third beam splitting pulse after another free evolution time T . The detection reads out the relative population of the output ports. More elaborated schemes, spatially resolving structures in the atomic ensembles are available [START_REF] Dickerson | Multiaxis Inertial Sensing with Long-Time Point Source Atom Interferometry[END_REF][START_REF] Müntinga | Interferometry with Bose-Einstein Condensates in Microgravity[END_REF], but in the most simple case, the populations N 1 , N 2 of the two ports are just counted. The normalized output signal then reads

P (φ a ) = N 1 /(N 1 + N 2 ) = P 0 + A • cos(φ a ) (4)
with offset P 0 , amplitude A, contrast C = A/P 0 , and accumulated phase φ a . Neglecting other effects, the low expansion rate enables high beam splitting efficiency, and consequently a contrast close to unity C ∼ 1. Neglecting finite pulse durations [START_REF] Antoine | Rotating matter-wave beam splitters and consequences for atom gyrometers[END_REF], and omitting other / higher order terms, the phase shift due to accelerations a is calculated [START_REF] Bordé | Quantum Theory of Atom-Wave Beam Splitters and Application to Multidimensional Atomic Gravito-Inertial Sensors[END_REF]Arimondo et al., 2009) to

φ a = k • aT 2 (5)
where k denotes the effective wave vector corresponding to the four photon momenta. Here, the three-pulse sequence is preferred to the four-pulse sequence also discussed in [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF], because it offers a higher contrast and a larger scale factor, scaling as ∼ T 2 .

Additional elements important for a light-pulse atom interferometer are the source generating the atomic ensemble, and the detection system, reading out the output ports of the atom interferometer. For each cycle, the source prepares an atomic ensemble with a number of atoms N , and a finite velocity spread σ v which corresponds to an (effective) temperature. In CAI [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF], the generation of 87 Rb Bose-Einstein condensates is proposed for the source to enable very low residual expansion rates at a level of 100 µm/s (Kovachy et al., 2015b;[START_REF] Müntinga | Interferometry with Bose-Einstein Condensates in Microgravity[END_REF].

The fundamental limit in such an interferometer is the quantum projection noise (QPN), implying a minimum phase noise per cycle of

σ φ QP N = 1/(C √ N ) (6)
Consequently, the QPN limited acceleration sensitivity is

σ a = 1 C √ N kT 2 tc τ
for a cycle time t c , an integration time τ , k = |k|, and acceleration a parallel to k. Manipulating two of these atom interferometers, spatially separated by a baseline l, simultaneously with the same beam splitter light fields forms the gradiometer geometry [START_REF] Mcguirk | Sensitive absolutegravity gradiometry using atom interferometry[END_REF]. The gradient information is recovered from the differential signal divided by the baseline. In this case, the QPN limited sensitivity to an acceleration gradient is

σ Γ = √ 2 C √ N kT 2 l t c τ (7)
for an acceleration gradient Γ in direction of the beam splitters, k • l = kl. In [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF], an interleaved measurement mode is proposed. Several atom interferometers are operated concurrently with a cycle time t c = 1 s.

With N = 10 6 , k = 8π/(780 nm), T = 5 s, and l = 0.5, a sensitivity of 3.5 mE/ √ Hz is obtained by [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF], which corresponds to a white noise with a one-sided ASD of 3.5 × √ 2 ≈ 5 mE/ √ Hz. In principle, the trace of the gradient tensor can be measured by reproducing this set-up in all three axes. In a nadir orientation, an additional constraint has to be considered. A typical GOCE-like altitude would imply an orbital frequency of Ω ∼ 1 mrad/s which affects the interferometers oriented along track and in the nadir direction. Due to the rotation, the trajectories do not close anymore at the third beam splitting pulse, implying a significant contrast reduction. The two trajectories will miss each other by ∼ 1.4 mm (in direction of forward drift velocity v) which is significantly larger than the thermal de Broglie wavelength of the atoms of ∼ 19 µm. Without sufficient overlap, the contrast is lost, rendering the gradiometer insensitive. Due to the linear dependency, a reduction of the rotation rate by a factor of ∼ 1000 would be required to recover the contrast. Additionally, a rotation induced phase shift

φ rot = 2(k × v) • ΩT 2 appears if a drift velocity v is present.
In the CAI proposal, the atoms are launched into the interferometer with a forward drift velocity v ≈ 2.4 cm/s, v = |v|, |k × v| = kv. This shift is cancelled in the differential signal between the accelerometers which contains the information of the gradient, but only if the forward drift velocities are the same for both. Due to the finite velocity spread and the finite number of atoms per ensemble, a random velocity uncertainty of σ v / √ N which is uncorrelated between the interferometers will remain. This limits the maximum rotation rate compatible with shot noise limited performance to ∼ 1 µrad/s, including some margin. This value is similar to the requirement for closing the interferometer. Moreover, the finite velocity distribution of the ensemble implies a velocity dependent dephasing. When averaging over an output port by counting the population, the signal vanishes (C → 0), since a fringe pattern inside the output port appears for a rotation rate significantly exceeding ∼ 1 µrad/s leading to 2(k [START_REF] Tackmann | Self-alignment of a compact large-area atomic Sagnac interferometer[END_REF]). Without further modifications, a maximum rotation rate of ∼ 1 µrad/s can only be maintained in the cross-track direction if a nadir pointing mode is chosen. Counter rotating the mirror to cancel the rotation rate [START_REF] Freier | Mobile quantum gravity sensor with unprecedented stability[END_REF][START_REF] Dickerson | Multiaxis Inertial Sensing with Long-Time Point Source Atom Interferometry[END_REF][START_REF] Lan | Influence of the coriolis force in atom interferometry[END_REF] is not an option due to the interleaved mode (several concurrent atom interferometers) and the huge dynamic range requirement of several tens of milli-radians. Note, that a gravity gradient parallel to k also leads to a non-closed atom interferometer, where the displacement at the last beam splitter pulse is parallel to k. This effect can be mitigated as described in [START_REF] Roura | Overcoming loss of contrast in atom interferometry due to gravity gradients[END_REF][START_REF] Roura | Circumventing heisenberg's uncertainty principle in atom interferometry tests of the equivalence principle[END_REF]). An alternative is the operation in a quasi-inertial pointing mode with maximum spurious rotation rates of 1 µrad/s in all three axes. These two pointing modes, nadir and quasi-inertial, motivate the performance discussion in the following sections.

•σ v )•ΩT 2 0.1π ((

Framework and method of the simulation

Objectives of a potential CAI gradiometer mission

Before starting any simulation, it is important to clearly specify the scientific objectives of a potential CAI gradiometer mission, keeping in mind the expected noise level of the sensor. These objectives will be a key factor in the choice of the satellite orbit which is an essential parameter of the simulation study. As previously stated in the introduction, the various proposals for future space gravity missions emphasize the need to sustain a continuous monitoring of the time-variable gravity field. However, and as already noticed in [START_REF] Carraz | A spaceborne gravity gradiometer concept based on cold atom interferometers for measuring Earth's gravity field[END_REF], a CAI gradiometer in the configuration described in the previous section is not sensitive enough to detect the corresponding signals with a sufficient signal-to-noise ratio (SNR) on a monthly basis.

To get an idea of the magnitude of the time-variable gravity signal as compared against the CAI gradiometer sensitivity, the amplitude spectral densities of the time-variable V xx , V yy and V zz in the orbital frame of a circular orbit at 259 km altitude are plotted in Figure 4 together with the expected gradiometer sensitivity. The time-variable gravitational signal is computed from one month (April 2006) of the updated Earth System Model from ESA [START_REF] Dobslaw | The updated ESA Earth System Model for future gravity mission simulation studies[END_REF], which contains the contribution of the atmosphere, oceans, cryosphere, land hydrology and solid Earth. The x, y and z-direction defining the gradiometer reference frame are outlined in Figure 3. One can see that a sensitivity below 0.1 mE/ √ Hz is required if these gradients are to be detected with an SNR larger than 1 in the frequencies between 10 -3 and 10 -2 Hz.

We have therefore decided to focus our study on the ability of a CAI gradiometer mission to recover the static gravity field. As such, the GOCE mission represents a perfect benchmark to evaluate the performance of the CAI gradiometer concept and to quantify precisely how the latter can improve our knowledge of the static field. It has thus been decided to choose for the simulations a GOCE-like orbit so that, in addition to the true error computation in the closed-loop simulations, a comparison of the GOCE-derived gravity models and the ones derived from the simulated data can be carried out.

Figure 3: Orientation of the axes of the gradiometer reference frame for a nadir pointing mode. The x-axis is parallel to the along-track direction, the y-axis is perpendicular to the orbital plane while the z-axis points in the nadir direction.

Description of the closed-loop simulator

The main results of this study are based on the output of a closed-loop simulator designed in MATLAB, which enables to synthesise realistically degraded observables, namely the estimated orbit, the satellite angular velocity and attitude and the gravity gradients from which the gravitational model is finally derived. The same approach has been used for both the nadir and quasi-inertial mode, with different processing in what concerns the attitude. Starting from a global gravity field model and an orbit sampled at 1 Hz, the GGT is computed at each measurement point in the LNOF (Local North Oriented Frame). Then the GGT is rotated into the gradiometer Reference Frame which axes define the measurement axes of the gradiometer. This requires to know for each time the true attitude of the GRF with respect to a reference frame, say the IRF. In the third step, the GGT components (V yy in the nadir case, the three diagonal elements in the quasi-inertial case) are degraded by adding the gradiometer noise and error terms due to the imperfect correction of the centrifugal part. This step also requires to have a realistic model of the estimation error of Ω. Based on the noisy Ω, the attitude quaternions of the GRF with respect to the IRF are computed and transformed into rotation matrices. Combined with the noisy gradients and the satellite estimated position, these data are used in a Least-Squares (LS) adjustment (see section 3.3) to recover gravity field models which are eventually compared to the initial one.

In this study, we have considered that the different observables were only corrupted by additive noise, that is we assumed a white noise floor of the sensors, as it would be the case for a solely shot-noise limited atom interferometer, ideal suppression of spurious accelerations between the two atom interferometers in gradiometer configuration, and neglected higher order terms.

Gravity field recovery

As the comparison is done at the level of gravity field models, we briefly describe the method used to recover them from the simulated gravity gradients. The global gravity field models are usually expressed in terms of Spherical Harmonic (SH) series [START_REF] Hofmann-Wellenhof | Physical geodesy[END_REF]:

V (r, θ, λ) = GM R N n=0 R r n+1 n m=0
C nm cos(mλ) + S nm sin(mλ) P nm (cos θ), (8) where V is the gravitational potential, GM is the gravitational constant of the Earth, R is the equatorial radius of the reference ellipsoid, (r, θ, λ) are the spherical coordinates (r radius, θ co-latitude, λ longitude), n, m denote SH degree and order, N denotes the maximum degree of the model expansion, P nm (cos θ) are the fully normalized associated Legendre functions, and C nm , S nm are the normalized SH coefficients, also known as Stokes coefficients, which are the unknowns of the gravity field solution.

The analytical expression of the GGT V LN OF in the local north-oriented frame can then be directly derived from equation ( 8). However the gravity gradients are measured in the GRF. It is therefore necessary to bring the measured gravitational gradients and the Earth's gravity field model in the same reference frame. Therefore we rotate the analytical expression of the gravitational gradients from the LNOF to the GRF. The rotation is given by:

V GRF = RV LN OF R t , (9)
where R, the rotation matrix between the two reference frames is computed from the attitude quaternions. To resolve the observational equation, the Least-Squares (LS) adjustment method is applied. Concerning the stochastic model, a unity matrix is used mostly because white noise is assumed in the simulation. But in some cases, the spectral density of the post-fit residuals exhibits a coloured behaviour, see more details in Section 5. In such cases, a full variance/covariance matrix (VCM) of observations is considered in the LS adjustment. The VCM is estimated from the post-fit residuals and further treated as the empirical stochastic model in an iterative step. Due to the large dimension, it is impossible to compute and store a complete VCM for the whole data set. Alternatively, the VCM is computed arc-wisely, with an arc length of 15 minutes. Because of the orbit inclination of 96.5 • , there is no observation in both polar regions. This is known as the polar gap problem [START_REF] Sneeuw | The polar gap[END_REF]. The polar gap will severely affect the quality of the zonal and nearzonal SH coefficients (orders around 0). These coefficients are therefore omitted in the computation of the error degree variances in the rest of the study. As the magnitude of the high-degree harmonics of the gravitational gradients decreases rapidly with the altitude, it is crucial to have an altitude as low as possible. For the sake of rigorous comparisons to GOCE, we have chosen 71 days of the GOCE orbit starting from the 2 nd of March 2013. During this period the circular orbit had an average altitude of 239 km and a repeat cycle of 65 days. It must be pointed out that this altitude is lower than the altitude that GOCE had most of its lifetime which was typically of 255 km. The time series of the orbit coordinates have been interpolated at 1 Hz from the GOCE L2 dynamical orbit data sampled at 0.1 Hz using spline functions. The static GGT expressed in the LNOF is computed using the EIGEN-6c4 model [START_REF] Förste | EIGEN-6C4-the latest combined global gravity field model including GOCE data up to degree and order 1949 of GFZ Potsdam and GRGS Toulouse[END_REF] up to degree/order 360.

Satellite angular velocity

Nadir mode. The generation of the angular velocity Ω is split into two steps: first, by definition of the nadir mode, the instantaneous angular velocity ω lorf of the local orbital frame with respect to the IRF has been computed from the orbital position and velocity and then assigned to the component ω y so that the main rotation of the GRF at the orbital frequency is indeed about the y-axis. Second, a residual zero-mean angular velocity δΩ = (δω x , δω y , δω z ) t is added so that we finally have:

Ω =   δω x ω lorf + δω y δω z   (10)
This residual angular velocity consists of five successive harmonics of the orbital frequency with random phases and a 1/f decreasing amplitude, meeting the constraint of a total angular velocity perpendicular to the GRF y-axis smaller than 10 -6 rad/s.

Quasi-inertial mode. Following the definition of a truly inertial mode and keeping in mind the non-perfect attitude control and the notation used in the previous paragraph, the angular velocity vector of the GRF with respect to the IRF only consists of tiny residuals: Ω = (δω x , δω y , δω z ) t . However, we will see in section 5.2 that a better gravitational model can be recovered if the GRF rotates together with the orbital plane. In the typical case of a sun-synchronous orbit, this additional angular velocity is below 2 × 10 -7 rad/s, still within the margin of the requirement when added to a spurious rotation rate of 10 -6 rad/s. If we note Ω op = (ω op x , ω op y , ω op z ) t the angular velocity of the orbital plane with respect to the IRF, expressed in the GRF, we obtain:

Ω =   ω op x + δω x ω op y + δω y ω op z + δω z   (11)
This is the reason why it is more appropriate to call this mode the quasiinertial mode than the inertial mode.

Attitude quaternions

Rotations and attitude are computed using the quaternion formalism. The attitude of the GRF with respect to the inertial frame is calculated directly from the angular velocity described above and an initial condition by integrating the following differential equation:

q = 1 2 Ωq ( 12
)
where q is a unit quaternion describing the rotation and Ω is the angular velocity matrix of the frame. This equation has been numerically integrated with the method given in [START_REF] Zupan | Integrating rotation from angular velocity[END_REF]. For the nadir mode, the attitude is initialized with the quaternions describing the coordinates of the local orbital frame at the initialization time. This way we make sure that the satellite attitude will be a classical Earth-pointing attitude. To avoid any possible divergence of the error of the numerical integration, the attitude is reinitialized every day. This induces a tiny jump in the attitude at midnight but it has no impact on the results of the simulations.

For the quasi-inertial mode, a question arises: is there a best quasi-inertial attitude to recover the gravitational field? For the simulations we have decided to keep an attitude for which one axis, namely the y-axis, is always perpendicular to the orbital plane so that the V yy gradient will be the same as in the nadir mode. Therefore the initial quaternions are chosen so that the y-axis is perpendicular to the orbital frame at the initialization time.

Generation of noisy observables 4.2.1. Error of the orbit determination

For the orbit determination, we assume a zero-mean, normally distributed error with a standard deviation of 2.6 cm added to the orbit inertial coordinates. This precision was typically achieved with GOCE [START_REF] Bock | GPS-derived orbits for the goce satellite[END_REF].

Noise of the angular velocity and attitude quaternions

For the simulation we assume that the satellite payload includes three star-trackers and a gyrometer for the determination of the angular velocity. Both kinds of sensor give estimates of the three components of the angular velocity. Star-trackers and gyrometer outputs are optimally combined at the level of angular velocity and in the frequency domain using a Wiener filter like in the processing of the angular velocity of GOCE (see [START_REF] Stummer | Alternative method for angular rate determination within the GOCE gradiometer processing[END_REF]). This method requires to have a model of the power spectral density (PSD) of the noise of both sensors. Thereafter, we assume we have an exact knowledge of their spectral characteristics. The gyrometer considered here is the space fiber-optic gyrometer ASTRIX R 200 manufactured by Airbus Defence & Space. We considered a white noise of 6.8 × 10 -6 rad/s/ √ Hz, as guaranteed by the constructor. For the star-trackers, we assumed that the noise PSD for the three axis follows the model derived in [START_REF] Stummer | Alternative method for angular rate determination within the GOCE gradiometer processing[END_REF] for the y-component of the angular velocity for GOCE. Both noise spectra are plotted in Figure 5 along with the example of the spectrum of a combined solution. Below 4 × 10 -5 Hz the angular velocity information is essentially given by the star-trackers while above it is only the gyrometer. The estimated attitude quaternions are deduced from the noisy angular velocity by integrating equation (12).

Noise on the gradients in the GRF

The diagonal gravitational gradients are deduced from equations (2). The final error on the gravitational gradients is therefore the sum of the gradiometer noise and the error due to the correction of the centrifugal terms. The expression of the error δV xx is given in equation ( 13) for the term V xx , where the measured quantities are topped with a ∼: where n ωy (resp. n ωz ) is the additive noise on ω y (resp. ω z ) discussed in the previous section and n Γ is the gradiometer white noise. Identical expressions can be derived for the two other gradients, for the nadir and the quasi-inertial mode.

Ṽxx = Γxx -ω2 z -ω2 y = (V xx + ω 2 z + ω 2 y + n Γ ) -(ω z + n ωz ) 2 -(ω y + n ωy ) 2 = V xx -2ω z n ωz -2ω y n ωy -n 2 ωz -n 2 ωy + n Γ δVxx (13)
One must also pay attention to the impact of the error of the attitude determination on the estimated gravitational gradients. We can show for analytical purposes that the additional error at δV on the GGT caused by the attitude determination error is given at first order by the following expression:

at δV = dΘV + VdΘ t ( 14
)
where V is the GGT expressed in the GRF and dΘ is the skew-symmetric matrix of the attitude error.

In particular, we have for the diagonal gradients:

at δV xx = -2dθ z V xy + 2dθ y V xz (15) at δV yy = 2dθ z V xy -2dθ x V yz (16) at δV zz = -2dθ y V xz + 2dθ x V yz ( 17
)
where dθ x dθ y and dθ z are the errors of the Euler angles describing the rotation between LNOF and GRF.

For the sake of simplicity, we approximate here the GGT V lnof expressed in the LNOF along the orbit by a diagonal matrix with V lnof .5 E and V lnof yz ±1 E. In the nadir case, the GGT expressed in the GRF will remain essentially diagonal. As a consequence, we only expect a minor impact of the attitude error on the estimated V yy . Assuming a perfectly circular orbit and under the previous hypothesis, the GGT in the GRF for the quasi-inertial mode is given by

V =   686 -2058 cos(2ω lorf t) 0 -2058 sin(2ω lorf t) 0 -1372 0 -2058 sin(2ω lorf t) 0 686 + 2058 cos(2ω lorf t)   . ( 18
)
In this mode, at δV xx and at δV zz will be dominated by the partial projection of V xz . If we note σ dθy the standard deviation of the angle error dθ y , then the standard deviation of the induced error on the gradient becomes 2|V xz |σ dθy = 4116σ dθy |sin(2ω lorf t)| E. With σ dθy ≈ 1.7 × 10 -6 rad in our case, the error will have an amplitude of 7 mE modulated at twice the orbital frequency. We therefore expect the attitude error to be a critical issue in the quasi-inertial mode and a major contributor to the noise contaminating both V xx and V zz .

Results

In this section, global gravity field models recovered from gravity gradients in both nadir and quasi-inertial mode are resolved up to degree/order 240. Both the true errors (the differences between the recovered gravity field coefficients and the input background model, i.e. EIGEN-6C4) and the formal errors (the accompanied standard deviations of the parameters, obtained in the LS adjustment) are analysed.

Nadir mode

Only the gravity gradient in the cross-track direction V yy is measurable in this case. We first investigate the impact of the gradiometer noise by running simulations for 4 different levels of white noise (cases 1 to 4). The different cases are summarized in Table 1 and the error degree variances are shown in terms of geoid height in Figure 6(a).

As expected, the gravity field model derived from observations with a noise of 7.0 mE/ √ Hz performs worst while the model derived for a noise of 2.5 mE/ √ Hz performs best. In between, the performances of the gravity field models are gradually improved when the noise is decreasing. The model Hz noise is about twice better than that from the larger noise of 5 mE/ √ Hz, especially in the spectral band above degree 50. This is coherent with the levels of input noise. Furthermore, the patterns of the true and formal errors coincide with each other. This indicates that the error propagation of the simulation model is correct. For an identical orbit and timespan, Figure 6(b) shows that the solution derived from the GOCE V yy is slightly better than the CAI gradiometer with a nominal noise (case 2) for degrees larger than 100. The GOCE solution derived from the combination of V xx , V yy and V zz (not plotted here) is even better over degrees larger than 20. We therefore conclude that under the same conditions (orbit and duration), the CAI gradiometer in nadir mode with a noise of 5 mE/ √ Hz does not outperform GOCE. The pyramid representation of the model errors are also shown in Figure 7. We find that the sectorial and near-sectorial coefficients are determined with high precision. This is due to the fact that V yy is approximately oriented in cross-track direction so that the gravity field features in this direction are better mapped. However, some faint stripes appear in the pyramid representation of the true errors at orders of m = 16 × k, k ∈ N. From the ASD of the post-fit residuals of case 2 and 4, displayed in Figure 8, the peaks at frequencies of multiple (with even numbers) cycles per revolution (cpr) are supposed to be responsible for these stripes since the coefficients at orders m = 16 × k are particularly sensitive to those frequencies. A similar conclusion is also drawn in the processing of GOCE gravity gradients, see [START_REF] Stummer | Alternative method for angular rate determination within the GOCE gradiometer processing[END_REF]. To further explore the reason for these stripes, the effect of the quaternion noise is investigated. We generate two new data sets, case 5 and case 6, with the same noisy V yy as in respectively case 2 and case 4 but with noise-free quaternions. In terms of degree variance (see Figure 6(b)), there is no evident difference between the models (case 2 vs. case 5, and case 4 vs. case 6). But in the pyramid representation of the true errors in Figure 7, most stripes disappear in the solutions where the noise-free quaternions are used. Only one tiny stripe is still visible at the order of about 120. It is supposed to be caused by the sub-cycle of the repeat ground track orbit which is about eight days. The ASDs of the post-fit residuals are also compared in Figure 8. With the noise-free quaternions, the ASDs of the residuals at the very low frequencies (below about 2 cpr) are significantly reduced and the amplitude of the peaks at multiples of 2 cpr are slightly attenuated. This agrees with the improvement of the coefficients at orders m = 16 × k.

Quasi-inertial mode

In the quasi-inertial mode the three diagonal gravitational gradients are measured. As a consequence, we would expect better results than in the nadir mode. Here, the most critical point to investigate is the impact of the error of the estimated GRF attitude. Contrary to the nadir mode, a significant part of the energy of the GGT is located in the off-diagonal elements of the tensor, as outlined in section 4.2.3. In particular, any error in the knowledge of the GRF orientation will introduce a partial projection of V xz on the diagonal gradients V xx and V zz , which is much larger than the gradiometer noise.

We carry out two simulations, case 1 and case 2, for which the gradiometer noise is the same and set to its nominal value, i.e., 5 mE/ √ Hz. The difference between both cases is on the noise of the gyrometer angular velocity. In case 1, the angular velocity noise is at its nominal value while in case 2 the noise is divided by 2. The consequence is twofold: first the error on the attitude quaternions is lower in case 2 and second, the error due to the correction of the centrifugal terms in case 2 is also smaller than in case 1. For both cases, we derived solutions from each gravitational gradient taken separately (component-wise solution) and a solution combining the three gradients (combined solution).

The error degree variance of the true errors of the estimated spherical harmonic coefficients in terms of geoid height error for both cases are plotted in Figure 9. As expected, the component-wise solutions from V xx and V zz have a larger error compared to the one derived from V yy while the latter compares quite well with the nadir solution (nadir case 2). The solutions derived from V xx and V zz show a similar error which is not surprising since they are contaminated by the same dominant error related to the attitude error and the projection of V xz . The regular oscillations on both error are due to the modulation of this error at twice the orbital frequency. The comparison between the solutions of case 1 and case 2 confirms however the slight improvement of case 2, as can be seen from their respective combined solutions.

The quality difference of the models recovered from V yy on the one hand, and V xx and V zz on the other hand, is even better illustrated in the pyramid representation of the model errors in Figure 10. In this figure, one can see that stripes at orders m = 16×k, k ∈ N degrade significantly the gravitational models derived from V xx and V zz , and to a much lower extent from V yy . These stripes correspond to the regular oscillations in the degree variance curves in Figure 9. Like in the nadir case, we also observe the presence of peaks at 2 cpr and multiples in the ASDs of the post-fit residuals of V xx and V zz plotted in Figure 11. They correspond to the error contaminating V xx and V zz which is essentially caused by the attitude error and the partial projection of the V xz component also modulated at 2 cpr. This specific error on the gradients is directly mapped onto the Stokes coefficients of order m = 16 × k, k ∈ N.

We can thus conclude that these stripes are directly related to the attitude error.

It is noteworthy that V yy is less affected by these stripes. The reason is that we have forced the y-axis of the gradiometer to be always perpendicular to the orbital plane. In this way, we make sure that the projection of the off-diagonal gradients on V yy due to the attitude error always have a weak amplitude, just like in the nadir mode. This is the reason why the GRF must follow the rotation of the orbital plane: if we keep a truly inertial attitude then the y-axis would slowly rotate with respect to the orbital plane and eventually be parallel to it, like V xx or V zz initially. Accordingly, the noise on V yy due the attitude error would become similar to the one on the 2 other gradients, with the same notable stripes.

A first idea in order to mitigate these stripes would be to have a more accurate and precise attitude determination for the quasi-inertial case than for the nadir case. We will see however that a simple processing can largely reduce these stripes. We introduce empirical parameters to absorb the lowfrequency contributions that are contained in the post-fit residuals. More specifically, a drift and a bias term are added to the functional model and estimated in an arc-wise manner. The length of the arc has been fixed by trial and error to 15 minutes. These parameters are estimated together with the Stockes coefficients. Following this method, the gravity field solutions of case 1 are re-computed and shown in Figure 12. This figure also includes the plots of the ASD of the corresponding post-fit residuals. The amplitude of the peaks at 2 cpr and multiples are dramatically reduced for all gradients, as compared to the previous situation. As expected, the error level and oscillations of the degree variances for the component-wise solutions of V xx and V zz and the combined solution are also attenuated. This is also reflected in the pyramid distribution of the true errors in Figure 10 where the stripes almost totally disappeared.

It can be pointed out that the post-fit residual ASDs drop below 1 mE/ √ Hz for frequencies smaller than 2 × 10 -4 Hz, which -according to Figure 4-is approximately the magnitude of the time-variable gravitational signals in this frequency domain. Yet, this does not allow to detect the time-variable signal. First because the arc-wisely estimated empirical parameters tend to absorb the signal and noise at low frequencies (at least lower than 1 15×60 ≈ 1.1 × 10 -3 /Hz) and this effect increases as the frequency decreases. There- fore this signal would not be captured by the gravitational model. Second, the gravitational signal at frequencies smaller than 1 cpr ≈ 1.8 × 10 -4 /Hz is related to spatial variations of the field whose wavelengths, when sensed along the orbit, is greater than 1 revolution. It is therefore a complex mix of spatially and / or temporally aliased gravitational signal, also modulated by variations of the orbit altitude. In our case, where only a static gravity field is considered, we can see from the example of V xx in Figure 12 that the absorption of these low frequencies does not degrade the recovery of the field. Although white noise is assumed for the gradiometer, we see that the spectral shape of the post-fit residuals exhibits a slight colour behaviour, especially in the low-frequency part. To further improve the estimation of the gravity field solutions, we run a second analysis where the variance/covariance matrix is taken into account to model the correlation of the observations. We notice a tiny improvement for degrees smaller than 100, however the error above degree 100 remains at the same level as before.

Discussion

We analyse now to what extent the CAI gradiometer concept can outperform GOCE when assuming a nominal gradiometer noise of 5 mE/ √ Hz. To this end, we extrapolate the error on these models from 71 days to 8 months assuming the error is purely stochastic and reduces as √ t i where t i is the time of integration. To gain an idea of the GOCE solution error, we first compute a gravity model based on the three diagonal gravitational gradients of the whole GOCE mission period (November 2009 -October 2013, about 47 months). The least-squares adjustment is performed for GOCE data in the same way as for our simulations. Then, we compare it to 2 different reference gravity models: on the one hand, the EIGEN-6C4 model used for the synthesis of the gravitational gradients in the simulation and, on the other hand, the GOCO05S (Mayer-Guerr, 2015). Since the comparison concerns only the contribution of the gradiometer, we have not taken into account the high-low satellite-to-satellite tracking data used to recover the low degrees of the gravity model. The results are shown in Figure 13 in terms of geoid height degree variance. We can note two bumps on the GOCE solution when it is compared to EIGEN-6C4. One appears in the range between d/o 100 and 130; the other one is visible between d/o 170 and 200. These bumps are due to the background model errors. This appears clearly when we compare our derived gravity field solutions to a satellite combined solution, i.e., GOCO05s. The bumps disappear in the GOCE solution but appear in the quasi-inertial solution based on EIGEN-6C4. Besides these two bumps, we can make the following comments:

• As already mentioned in section 5.2, for identical conditions, the solution in the quasi-inertial mode is significantly better than the nadir mode one, especially above d/o 50. Additionally, the nadir mode with For GOCE, the whole period combined solution is based on the gravitational gradients of the whole mission period while the low orbit period solution correspond the 8 last months during which the satellite was at an altitude equal or lower than 239 km.

a nominal noise does not allow to outperform GOCE under identical conditions (orbit and measurement timespan), see section 5.1. As a consequence, the nadir mode should be discarded in favour of the quasi-inertial mode.

• The 8-month quasi-inertial solution is better than the GOCE solution for the whole mission period. We can thus conclude that an 8-month mission at an altitude of 239 km and using a 3-axis CAI gradiometer with a nominal white noise of 5 mE/ √ Hz in the quasi-inertial mode would outperform the full GOCE mission and eventually yields a more precision gravity field model. Nonetheless, it should be borne in mind in this comparison that the GOCE satellite was most of its lifetime at a higher altitude than 239 km, about 3 months at 239 km and 5 months lower than 239 km. This fact shows however that it is technically possible to fly a satellite at an altitude as low as 239 km for a duration of 8 months.

• For a fairer comparison, we have also plotted the GOCE solution for the 8 months during which the satellite was at an altitude equal or lower than 239 km (see Figure 13). Again, the 8-month quasi-inertial solution at a constant altitude yields a slightly better solution despite the fact that the 8-month GOCE solution is partly based on gravity gradients measured at a lower altitude.

Among the assumptions made in this study, the attitude control requirements and the quasi-inertial mode seem particularly challenging, especially at such a low altitude where the atmospheric drag may cause variations on the satellite angular velocity greater than 1 µrad/s. Inertial pointing satellites are rather uncommon although control feedback strategies to maintain such a pointing have been studied [START_REF] Lovera | Global magnetic attitude control of inertially pointing spacecraft[END_REF] and robust solutions using only magnetorquers have been put forward, see e.g. [START_REF] Celani | Robust three-axis attitude stabilization for inertial pointing spacecraft using magnetorquers[END_REF] and [START_REF] Celani | Spacecraft attitude stabilization using magnetorquers with separation between measurement and actuation[END_REF]. The MICROSCOPE space mission [START_REF] Touboul | The MICROSCOPE experiment, ready for the in-orbit test of the equivalence principle[END_REF] launched in April 2016 is an example of a satellite partly operated in inertial pointing with a specification on the angular velocity stability smaller than 10 -9 rad/s [START_REF] Touboul | The MICROSCOPE experiment, ready for the in-orbit test of the equivalence principle[END_REF]. However such an inertial pointing stability is achieved after a long in-orbit calibration and at an altitude of 720 km, where the atmospheric drag is considerably smaller than at the altitude we assumed in our simulation.

Summary, conclusions and recommendations

Recent developments in light-pulse atom interferometry have opened the way to the design of new high-accuracy and high precision inertial sensors, in particular when operated in a microgravity environment where the operating conditions are optimal. Used in an accelerometer configuration and in pairwise assemblies, atom interferometers constitute the backbone of the space gravitational gradiometer concept introduced in this article. A nominal white noise of 5 mE/ √ Hz is expected for the CAI gradiometer. Its performance and ability to recover the Earth's gravitational field have been investigated via closed-loop simulations assuming realistically degraded ancillary data (position, angular velocity, gradiometer noise and attitude quaternions).

Because the atom interferometer architecture requires to limit the angular velocity perpendicular to the measurement axis to 1 µrad/s, only two operating modes are possible for the considered CAI gradiometer concept: the nadir mode and the quasi-inertial mode. In the nadir mode only the V yy component, along the across-track direction, is measured. With one observed gradient, our simulations show that the nadir mode does not allow to clearly outperform GOCE. The results are more promising in the case of the quasi-inertial mode where the gradiometer can measure the three diagonal components of the GGT. To outperform the whole GOCE mission in term of geoid accuracy, we conclude that a CAI gradiometer operated in quasiinertial mode for at least 8 months, on a circular orbit and at an altitude of 240 km is required. Of course, this simple scenario is just an indication and further studies are required to find the best trade-off between the orbit altitude, the mission life-time and the technical limitations.

Although the quasi-inertial mode proved to be the best operating mode for the CAI gradiometer considered in this study, keeping such a pointing seems technically challenging. Therefore we recommend to investigate whether it is possible to modify the design of the CAI gradiometer so that the 3 diagonal elements of the GGT can be measured in the nadir mode. In other words, an improved interferometer architecture that allows the atoms to interfere in an almost uniform rotation reference frame is highly desirable. A solution for the compensation of the Earth rotation in ground based set-ups has already been implemented [START_REF] Lan | Influence of the coriolis force in atom interferometry[END_REF]. It consists in repeatedly counter-rotating the reference mirror using tip-tilt actuators. However, implementing such a solution would require in our case large tilt chirps of the mirror, if a single mirror is to be used. This is not readily implementable in the case of our interleaved interrogation geometry. As an alternative, we are currently studying an architecture based on a set of mirrors with pre-adjusted tilts, which would allow to cancel the effect of rotation at the orbital frequency.

With the 3 diagonal elements of the GGT measured in the nadir mode, better performances are expected, particularly thanks to the V zz component which would be estimated directly along the radial direction. Furthermore, the off-diagonal gradients would keep a small amplitude, leading to a much smaller impact of the attitude error determination on the error of the estimated gradients. like to thank Sven Abend for comments and acknowledges financial support from "Niedersächsisches Vorab" through the "Quantum and Nano-Metrology (QUANOMET)" initiative within the project QT3.

Figure 1 :

 1 Figure 1: Atom-light interaction for beam splitting.(1) The laser beam features two frequencies ω 1 , ω 2 and is retro reflected at the mirror. This leads to four different components in the light, distinguished by frequency and / or propagation direction (k). Depending on the duration of the interaction of the light field with the atoms, the atoms are coherently excited into a superposition state of initial momentum p 0 and p 0 + (k 1 -k 2 ) ( reduced Planck's constant) by absorption and stimulated emission of a photon from the incoming (wave vector k 1 ) and the retro reflected beam (wave vector k 2 ). During this process, the phase k • d is imprinted onto the atoms, depending on the distance of the atoms to the mirror d and the effective wave vector k = k 1 -k 2 . (2) In the special case of negligible velocity of the atom parallel to the effective wave vector k, a four photon process leads to a symmetric momentum transfer with the superposition of the three momentum states p 0 and p 0 ± (k 1 -k 2 ), effectively doubling the effective wave number.

Figure 4 :

 4 Figure 4: Comparison of the Amplitude spectral density of the expected CAI gradiometer noise to the time-variable signal in V xx , V yy and V zz in the orbital frame of a circular orbit at 259 km altitude.

4.

  Description of the simulation data 4.1. Generation of noise-free data 4.1.1. Orbit and gradients in the LNOF

Figure 5 :

 5 Figure 5: Amplitude spectral density of the noise degrading the estimation of one component of the angular velocity of the GRF with respect to the inertial frame.

Figure 6 :

 6 Figure 6: Error degree variance of the geoid height for the gravitational models derived from 71-day measurements in nadir mode. Solid curves represent the true errors and dashed curves are formal errors. In Figure 6(b) the solution without attitude error (cases 5 and 6) and the model error derived from the GOCE V yy for 71 days of data at the same altitude are plotted for comparison.

Figure 7 :

 7 Figure 7: Pyramid representation of relative errors (logarithm scale) of the Stokes coefficients. The left column shows simulations with a nominal gradiometer noise of 5 mE/ √ Hz while in the right column the results for a gradiometer noise of 2.5 mE/ √ Hz are given. For the first and second row, the computation is made with noisy attitude quaternions while noise-free quaternions are used for the third row.

Figure 8 :

 8 Figure 8: Amplitude spectral density of post-fit residuals in the nadir mode.

Figure 9 :

 9 Figure 9: Degree variances of true errors for the component-wise (Fig.9(a) and 9(b)) and combined (Fig.9(c)) gravity field solutions derived from 71-day measurements in quasiinertial mode. For comparison, the solution in the nadir mode with nominal gradiometer and gyrometer noises (Nadir case 2 V yy ) is plotted.

Figure 10 :

 10 Figure 10: Pyramid representation of true errors (logarithm scale) of the Stokes coefficients. The left column shows component-wise and combined solutions of case 1 while the right column shows results of case 2.

Figure 11 :

 11 Figure11: Amplitude spectral density of the post-fit residuals in the quasi-inertial mode for case 1 and 2.

Figure 13 :

 13 Figure13: Extrapolation of the error of the gravity field solutions to 8 months for the CAI gradiometer and comparison to the error of the GOCE gradiometer-only solution. For GOCE, the whole period combined solution is based on the gravitational gradients of the whole mission period while the low orbit period solution correspond the 8 last months during which the satellite was at an altitude equal or lower than 239 km.

Table 1 :

 1 List and attributes of the different simulations run for the nadir case. nom. stands for nominal error level.

	case 1 case 2 case 3 case 4 case 5 case 6 Hz) 7 √ grad. noise (mE/ 5 4 2.5 5 2.5
	attitude err.	nom. nom. nom. nom. none	none
	derived from observations with 2.5 mE/	√

  Figure 12: Degree variances of the true errors of the gravity field model recovered in case 1 (left) and amplitude spectral density of the post-fit residuals (right) when empirical parameters are included in the LS adjustment. For comparison, the same quantities are plotted when no empirical parameters are considered for V xx .

		Obs. Vxx
		Obs. Vyy
		Obs. Vzz
		Res. Vxx w/i emp
		Res. Vyy w/i emp
	1/2 )]	Res. Vzz w/i emp Res. Vxx w/o emp

Acknowledgements

This study was supported by the European Space Agency through Contract No. 4000112677/14/NL/MP. The authors would like to thank all the collaborators of the project, in particular Baptiste Battelier, Andrea Bertoldi, Sven Herrmannn, Merle Cornelius, Naceur Gaaloul, Dennis Becker, Ernst Rasel and Norman Guerlebeck. Comments and feedbacks from Christian Siemes were greatly appreciated. K. Douch, J. Müller, C. Schubert and H. Wu also acknowledge the Collaborative Research Center (SFB) 1128 geo-Q, funded by the Deutsche Forschungsgemeinschaft (DFG). C. Schubert would


A Coustenis
email: athena.coustenis@obspm.fr

D E Jennings

R K Achterberg

G Bampasidis

C A Nixon

P Lavvas

V Cottini

F M Flasar

Seasonal Evolution of Titan's Stratosphere Near the Poles

Keywords: infrared, planetary systemsplanets and satellites, atmospheresplanets and satellites, compositionplanets and satellites, individual (Titan)radiation mechanisms, thermalradiative transfer

In this Letter, we report the monitoring of seasonal evolution near Titan's poles. We find Titan's south pole to exhibit since 2010 a strong temperature decrease and a dramatic enhancement of several trace species such as complex hydrocarbons and nitriles (HC 3 N and C 6 H 6 in particular) previously only observed at high northern latitudes. This results from the seasonal change on Titan going from winter (2002) to summer (2017) in the north and, at the same time, the onset of winter in the south pole. During this transition period atmospheric components with longer chemical lifetimes linger in the north, undergoing slow photochemical destruction, while those with shorter lifetimes decrease and reappear in the south. An opposite effect was expected in the north, but not observed with certainty until now. We present here an analysis of high-resolution nadir spectra acquired by Cassini/Cassini Composite Infrared Spectrometer in the past years and describe the temperature and composition variations near Titan's poles. From 2013 until 2016, the northern polar region has shown a temperature increase of 10 K, while the south has shown a more significant decrease (up to 25 K) in a similar period of time. While the south polar region has been continuously enhanced since about 2012, the chemical content in the north is finally showing a clear depletion for most molecules only since 2015. This is indicative of a non-symmetrical response to the seasons in Titan's stratosphere that can set constraints on photochemical and GCM models.

Introduction

The Cassini mission was in the Saturnian system observing, as a privileged target, Titan for about half of a Saturnian year. The last targeted Titan flyby was on 2017 April 22, followed by the start of the Cassini Grand Finale polar/proximal orbits the next day. Titan's northern hemisphere was in early winter when Cassini arrived in 2004 July and by the end of the extended mission, after 13 years in the system, it was just past northern summer solstice (2017). Within that time, we have had the opportunity to monitor Titan's atmosphere and surface response through two full seasons, including the same time as the Voyager encounter in 1980 November (Coustenis et al. 2013). Since the equinox in mid-2009, where both hemispheres received equal insolation, we have witnessed rapid and strong seasonal effects both in the atmosphere and on the surface of the satellite.

In particular for Titan's stratosphere, roughly between 50 and 400 km in altitude, variations in the thermal and chemical structure were reported from investigators analyzing spectra acquired by the Cassini Composite Infrared Spectrometer (CIRS; Flasar et al. 2004;Jennings et al. 2017). These investigations have shown that a large-scale reversal occurred in the single pole-to-pole circulation affecting the distribution of the atmospheric structure after the northern equinox in mid-2009, with the gases in the summer hemisphere showing upwelling while in the winter pole a corresponding subsidence was found (Achterberg et al. 2008(Achterberg et al. , 2011;;Bampasidis et al. 2012;Teanby et al. 2012). Whereas until 2011 a significant build-up of trace gases, haze, and condensates was still observed over the north pole (West et al. 2011(West et al. , 2016;;Jennings et al. 2012aJennings et al. , 2012b;;Coustenis et al. 2013;Vinatier et al. 2015), toward the end of that period, the vortex and hotspot reported in Jennings et al. (2015) and Achterberg et al. (2008) had practically disappeared. In its place, similar features started to appear near the south pole, which have evolved to what we see today.In the absence of sunlight, gases that subsidize have accumulated since 2012 and become enhanced (e.g., Coustenis et al. 2016 and references within;Teanby et al. 2017), and a subset of those have been suggested to be causing the haze decrease found in the north and its subsequent increase in the south (Jennings et al. 2012b(Jennings et al. , 2015)).

We have already published an analysis of spectra acquired by Cassini/CIRS at high resolution from 2010 onward, both in nadir and limb mode and at high and middle latitudes (Coustenis et al. 2013(Coustenis et al. , 2016)). In our 2016 work, we showed that the hydrocarbons found at high concentrations and longer lifetimes did not show significant latitudinal variability in the past years and until mid-2013. But in more recent dates, these molecules show a strong increase in the south. Notably, the 70°S and mid-latitude abundances show different trends indicating different dynamical processes at play inside and out of the polar vortex region extending down to 50°in latitude. While most trace stratospheric gases in the 70°N data generally show only a small decreasing trend until 2014, the 70°S results on the contrary exhibit a strong enhancement after 2012. This argues for a fast and strong build-up of the gases in the southern pole as it goes deeper into the shadow during the 2013-2014 onset of winter, as predicted by models. The limb spectra also show variations in time with altitude, with, in 1 particular, an area covering the whole northern hemisphere at altitudes between 400 and 500 km strongly depleted in gas and haze in early 2016, indicative of a complex dynamical effect yet to be identified (Vinatier et al. 2015;Teanby et al. 2017).

Here, we report on Cassini/CIRS observations of Titan's stratosphere at both high northern and southern latitudes (between 60°and 90°) from 2012 to 2016 (where we have calibrated high-resolution data) to infer the behavior of the atmospheric temperature and chemical structure as a function of time and the results of such seasonal effects. In particular, this Letter follows the Coustenis et al. (2016) work on seasonal effects in the atmosphere, looking for the expected decrease in abundance of the northern pole constituents, not previously observed.

Observations and Analysis

For the purposes of this Letter, we exploited CIRS high spectral resolution (0.5 cm -1 ) data taken from 2012 to 2016 in the surface-intercepting nadir mode. CIRS is a Fourier transform spectrometer spanning 10-1500 cm -1 in three spectral channels or focal planes (FPs). We use here spectra from the far-infrared FP3 and FP4 that cover the 600-1500 cm -1 range (Jennings et al. 2017).

Figure 1 shows the high-resolution spectra of Titan's stratosphere acquired at high latitudes (from 60°to 90°) from 2014 to 2017. The purpose of this figure is to expose unfortunate gaps in the observations for such high latitudes caused by the nature of Cassini's orbital geometry. Indeed, the north pole is not represented in the data between the fall of 2014 and 2016 May, and for the south pole we do not have enough high-resolution (0.5 cm -1 ) spectra available from fall 2014 onward. We then have to work with the existing data as far as temporal variations are concerned and to take into account these observational gaps.

Table 1 shows the spectra used in this work for inferring the time dependence of Titan's polar thermal and compositional structure in the stratosphere. With respect to previous studies, the data analyzed and interpreted here cover up to Titan Flyby 122 and essentially deal with recent data taken near the northern pole, where observations were missing until recently.

The selections that were already analyzed in August), so 7 out of a total of 12 data sets, are new. The new calibration provides a better sensitivity for line detection and the enhanced radiance in some of the parts of the FP3 spectrum allows for a more precise inference of the abundances and their variation with time for some of the weaker species. The difference in the results still remains within the 10%-15% range with respect to the abundances published in Coustenis et al. (2016).

The FOV of our CIRS observations was restricted to be entirely on Titan's disk and emission angles were limited to 0°-65°to avoid the need for large atmospheric corrections. To secure good sensitivity we only use FP3 averages with more than 500 spectra in a bin (with one exception in 2013 May at high southern latitudes) and with as many as about 1100 spectra in some cases. Available spectra within each time period were zonally averaged in 10°latitude bins for the south and 10°or 20°latitude bins for the north that covered latitudes within the 65°S-85°S range. The averages depend on the data available for a given date/flyby within a latitude range and aim at improving the signal-to-noise ratio (which is at worst 10 and at best around 50 on average; depending on the selection, the lower values of the S/N apply mostly to the southern data). Thus, in some cases, we had to extend the latitude range toward lower or higher latitudes, but we checked in our contribution functions that the atmospheric range probed remains in the same pressure levels. For purposes of clarity, all of these spectral averages are marked "SP" in the plots and discussion hereafter. We apply the same process for the data near the north pole and we mark their averages by "NP".

We use a monochromatic radiative transfer code adapted to Titan's atmosphere (ARTT; Coustenis et al. 2010;Bampasidis et al. 2012), and we extract the temperature profile from the best fit of the 7.7 μm methane band in FP4 (Figure 2 shows some examples of fitting the methane v 4 band in the 1250-1350 cm -1 region). Our model is based on a CH 4 abundance of 1.5% through the stratosphere and increasing below the cold trap to about 5% at the surface, in agreement with CIRS and with Huygens Gas Chromatograph Mass Spectrometer (GCMS) measurements (Flasar et al. 2004;Niemann et al. 2010). The temperature profile is originally set to the one measured at altitudes 0-147 km by the Huygens Atmospheric Structure Instrument (HASI; Fulchignoni et al. 2005). Adjustments were applied to the HASI profile to extend in altitude and to account for latitude dependences and seasonal variations in temperature as found by CIRS emission values in the 7.7 μm CH 4 band and from lower-resolution observations (Coustenis et al. 2007(Coustenis et al. , 2010(Coustenis et al. , 2016;;Achterberg et al. 2008Achterberg et al. , 2011;;Bampasidis et al. 2012). For more details on the method see the literature.

We then implement the temperature profile in our radiative transfer code and solve for the opacity in the rest of the spectrum retrieving the abundances of trace gases and their Note. Titan high-resolution data taken in nadir mode near the south and north pole in focal planes 3 and 4 of the CIRS instrument are analyzed here. Ls is the solar longitude at the time of the observation. The data sets marked with ( * ) were already included in our 2016 work but were re-analyzed here with improved calibration. isotopes in the FP3 region of the CIRS spectra. Figure 3 shows some of our fits in the 600-900 cm -1 range where we observe the emission of several gaseous bands (e.g., HCN, C 2 H 2 , C 3 H 8 , C 2 H 6 ). Our spectra are taken in the nadir mode so we do not have information on the altitude dependence of the mixing ratios, We apply constant-with-height vertical profiles from the higher levels of the atmosphere (around 400 km in altitude) down to the condensation level and then following the saturation curve. This assumption is essentially valid for all of the molecules considered here, due to their weak emission bands, except for C 2 H 2 whose abundance variations with height should be taken into account. Our method for inferring the acetylene mixing ratio values by making use of weighted averages over the P, Q, and R branches of the 729 cm -1 band and by testing vertical profiles is described in detail in Section 2 of Bampasidis et al. (2012).

Results and Discussion

In our recent work (Coustenis et al. 2016), we inferred the temperature profiles and the chemical composition at different dates from 2010 and until 2014 for high northern and southern latitudes (at and beyond 50°N and 50°S).

We had already shown then that as the southern hemisphere moved into winter after 2010, large temperature variations were observed near the south pole (70°S) in the stratosphere (from 0.1 down to 1 mbar pressure levels). Indeed, while a moderate warming is observed in the summer-entering north for the middle and high northern latitudes, a decrease of about 10-15 K in temperature was measured already for the 50°S latitudes. A more spectacular drop in temperature by as much as 25 K at 70°S was measured from 2012 to 2014. These temperature variations were accompanied by a strong enhancement of chemical compounds in the south polar region, while the north failed to show the opposite effect, which indicated a non-symmetrical reaction to the seasonal influence for each pole. We explore here more recent dates seeking to determine longer-term seasonal effects near Titan's poles.

Thermal Structure Variations in the Stratosphere Near the

Poles from 2012 to 2016

In Figure 5, we show the stratospheric temperature on Titan at 70°S with solid lines for different years and with different colors from 2012 to 2014 and at around 70°N with dashed lines for different dates from 2013 to 2016. Our contribution functions indicate that the emission observed in the 7.7 μm methane band originates from levels 0.01-20 mbar (70-400 km in altitude), with a peak around 0.1-1 mbar.

We note that we do not have nadir high-resolution CIRS data for the south pole after mid-2014 (Figure 1). Comparing the earlier dates in the northern polar region (light and dark blue and red dashed lines) with the other dates, we find relatively moderate changes: the new spectral averages from 2016 observations show that the temperature increased by 10 K from 2013 to 2016. The north pole then took three years to gain 10 K, while the south pole lost 2.5 times as much in only two years. There does not therefore appear to be a direct symmetry in the thermal response from each pole on Titan in the stratosphere in terms of temperature. In particular, a 6 K warming is found in the north during the period from 2014 February to September, whereas during approximately the same period of time, the southern stratospheric temperature dropped by an additional 12-15 K near 0.5 mbar, totaling up to a 25 K loss in a little more than two years. From 2014 to 2016, the north pole has gained another 2 K. The real increase in temperature in the north is actually significant only after 2014 September (Figure 4). The decrease in the stratospheric temperature in the south began as early as late 2011 (Coustenis et al. 2016). Our temperature profile for 2014 September is quite compatible with the results of limb observations in 2015 March by Teanby et al. (2017).

We then report here that the north polar region's temperature is beginning to show a reaction to the increasing sunlight as the north moves into summer, but with a three-year delay with respect to the high southern latitudes. This should be accompanied by a decrease in radiative emitters, as already observed for latitudes around for 50°N, a latitude located at the border of the northern vortex, which was shown to be weakening in the past years (Teanby et al. 2012;Jennings et al. 2015).

Chemical Composition Changes in the Stratosphere Near the Poles from 2012 to 2016

We have demonstrated that as the southern hemisphere moved into winter after 2010, strong temperature and composition variations are found essentially in the south pole. The drop in temperature over the south pole is associated with a strong enhancement in chemical content near the south pole. We thus witness in Titan's south pole the appearance and positive evolution of several trace species, such as HC 3 N and C 6 H 6 , previously observed only at high northern latitudes. This is due to the evolution of Titan's seasons, going from northern winter in 2002 to advanced northern spring in 2012 with the opposite effect in the south pole as it moves into winter. During that time, Titan's main atmospheric circulation cell reverses and models predict that it channels newly produced photochemical species from mid-latitudes toward the south pole.

Figure 5 shows our results for the abundances of some of the minor gaseous components near both poles, with a maximum of the contribution functions in the 0.1-10 mbar range (90-250 km). While subtle changes were seen for middle southern latitudes near 50°S (Coustenis et al. 2016), the red lines representing the south polar latitudes around 70°S in Figure 5 (indicated as SP) show an increase in the abundances by one or several orders of magnitude for many gases. For all of the trace molecules shown in this figure, the higher southern latitudes have become as or even more laden in organics than the north by late 2014. We thus note that the SP content was considerably increased from late 2010 until 2014 September, and in particular since 2012 by a factor of 2.5 for C 3 H 4 and almost an order of magnitude for C 4 H 2 . As was already shown in Coustenis et al. (2016), the increase is even more dramatic at SP (by at least three orders of magnitude) for C 6 H 6 , HCN, and HC 3 N. This enhancement has been noticed since late 2010 for HCN and HC 3 N and since late 2011 for C 6 H 6 (it should be noted that we only have upper limits prior to that date for these molecules; see Coustenis et al. 2016), and these species were significantly more enhanced at SP than at NP by late 2014 (by more than one order of magnitude for HCN and C 6 H 6 and by two orders of magnitude for HC 3 N).

When compared to the findings of Teanby et al. (2017) for the southern polar regions from nadir data, we find remarkable agreement for all the molecules in the 2014-2015 period, with the exception of the nitriles (HCN and HC 3 N) where our results indicate higher abundances already by about an order of magnitude. That could be due to the difference in spectral resolution (Teanby et al. use lower-resolution data), to the altitude ranges probed, and to the exact time periods. But all in all, our results are consistent.

The atmospheric constituents that were transported to the northern hemisphere during the previous season remain concentrated there around the north pole (indicated as NP in Figure 5) and are subjected to a slow photochemical destruction. Following the northern winter, those among the species with longer chemical lifetimes are expected to remain in the north for a longer time, while those with shorter lifetimes should be removed, reappearing in the south (Coustenis et al. 2016). The purpose of our investigation here is to determine how quick or slow this process is.

We do not have any high-resolution nadir data at 70°N before 2014, so we cannot infer the evolution of the abundances of these species prior to the past year, but they do not seem to have changed much since their values reported in Vinatier et al. (2015) from limb data in 2011 (their Figure 8). With that caveat in mind, our radiative transfer results show that the high northern latitudes (at about 20°from the pole) only became significantly depleted in trace gases after 2015 (Figure 5). No significant difference is found for CO 2 (not shown in Figure 5), as in previous studies. A limited drop in mixing ratios within the error bars is measured for the other most abundant molecules resisting the photolysis (such as C 3 H 8 , C 2 H 6 , C 2 H 2 , or C 2 H 4 ), given the uncertainties. However, the trend for a drop in content for these molecules is attested by the consistently lower values found for these species in all 2016 spectral averages with respect to 2014. More remarkable decreases in abundance at NP are found for C 3 H 4 and C 4 H 2 (by 50%-70%); 25%-35% decreases are found in the mixing ratios of C 6 H 6 , HCN, and HC 3 N.

The southern and northern polar values seem to cross (become similar) at around late 2014 for the hydrocarbons and earlier (at the beginning of 2014) for the nitriles.

Conclusions

We have evaluated the temperature and chemical composition variations in the stratosphere near Titan's poles from 2012 to 2017. The retrieved thermal profiles show that there is a three-year lag in the northern response with respect to the south.

Following the advent of winter in Titan's north, species with longer chemical lifetimes linger in the north while those with shorter lifetimes decrease and reappear in the south. Where we have data, we searched for signatures of weak stratospheric gases in large spectral averages around the poles and determined their abundances. The results show that the southern pole of Titan became hugely enhanced in trace gases after 2012. The northern composition decrease set in only about three years after the south enhancement had begun (roughly 2015 versus 2012).

The cause of these variations in abundance and temperature could be spatial changes (due to Titan's inclination) in the energy input to Titan's atmosphere, which acts as a driver for changes in the advection patterns, circulation, meteorology phenomena, and other physical processes. Those in turn lead to a stronger variability in the latitudinal abundances of the photochemical species after some time.

Among the most recent GCM works for Titan's atmosphere (Newman et al. 2011;Lebonnois et al. 2012;Lora et al. 2015), only the Institut Pierre-Simon-Laplace (IPSL) Titan GCM is coupled to atmospheric photochemistry, but no complete study of the seasonal variations of atmospheric composition has been published yet. According to a preliminary study (Lebonnois & Vinatier 2013), the decrease of modeled trace gas abundances in the winter polar region at pressures around 1-10 mbar starts as soon as spring equinox, and the enrichment of the opposite polar region begins as soon as the autumn equinox. In this GCM, both variations are synchronized together and with the equinox, which is not what has been observed with the CIRS data set, probably indicating that more work is needed to refine the model. This work then helps to set strong constraints on photochemistry and circulation models that need to account for the asymmetry found in the poles' behavior.

The lack of observations of Titan's atmosphere near the poles over long periods of time begs for a return to Titan with a dedicated mission, after the end of the highly successful Cassini-Huygens mission in 2017 September. Benefiting from the whole pool of CIRS spectra until the end of the mission should allow us to better constrain the degree of complexity attained by the chemistry in Titan's stratosphere and to search for new molecules.

Figure 1 .

 1 Figure 1. Cassini/CIRS high-resolution nadir data available from 2014 January until 2017 January for latitudes 60°-90°in the northern and southern hemispheres of Titan. This figure can be complemented with previous years of data from Figure 1 of Coustenis et al. (2016). It shows, among other effects, the lack of enough highresolution CIRS observations of the north pole between the fall of 2014 and 2016 May and of the south pole from 2014 September onward.

Figure 2 .

 2 Figure 2. Sample of fits in the 1200-1400 cm -1 part of the FP4 Titan spectrum containing the CH 4 emission band that gives access to the temperature profile. The blue curves are the observations by CIRS at different dates, and the red curves are the model simulations that best fit the data: (top) 2013 December near 70°S; (middle) 2013 December near 70°N; (bottom) 2016 August near 70°N.

Figure 3 .

 3 Figure3. Same as Figure2, but for a sample of fits in the 600-900 cm -1 part of the FP3 Titan spectrum: (top) 2013 December near 70°S; (middle) 2013 December near 70°N; (bottom) the same spectral region in 2016 August near 70°N. The feature close to 795 cm -1 not reproduced in the simulation corresponds to an interference pixel in the data that was not completely suppressed in some cases.

Figure 4 .

 4 Figure 4. Temperature evolution in Titan's stratosphere near the poles from 2012 to 2016. The southern polar profiles are shown with solid lines, while the north pole profiles are depicted with dashed lines and the different dates are indicated with different colors. 3σ uncertainties are shown.

Figure 5 .

 5 Figure 5. Abundances of trace hydrocarbon and nitrile gases in the stratosphere of Titan from 2012 to 2016 included for the south polar region (red lines-SP) and for the north polar region (black lines-NP). 3σ error bars are indicated.

  Coustenis et al. (2016) are re-investigated here with spectra from an updated calibration. This concerns the 2012 June and 2013 May southern averages and the 2013 December-2014 September northern averages. All of the other selections in Table 1 are new (three new south pole averages: 2012 September, 2013 December, and 2014 September, and four new north pole averages: 2013 September, 2016 May, 2016 June, and 2016

The authors acknowledge support from NASA's Cassini mission and Cassini data Analysis program. Our thanks to Sebastien Lebonnois for interesting discussions and for insights on the current status of GCM models.

ORCID iDs

A. Coustenis https://orcid.org/0000-0003-3414-3491 D. E. Jennings https:/ /orcid.org/0000-0003-1340-754X R. K. Achterberg https://orcid.org/0000-0002-7643-7626 G. Bampasidis https:/ /orcid.org/0000-0002-0109-7823 C. A. Nixon https:/ /orcid.org/0000-0001-9540-9121 P. Lavvas https:/ /orcid.org/0000-0002-5360-3660 V. Cottini https://orcid.org/0000-0003-0839-5855 F. M. Flasar https:/ /orcid.org/0000-0001-7186-5181