# MICROSCOPE mission: Data analysis principle

Joel Bergé1, Quentin Baghi1,2111Current address: NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA, Emilie Hardy1, Gilles Métris2, Alain
Robert3, Manuel Rodrigues1, Pierre Touboul1, Ratana Chhun1, Pierre-Yves
Guidotti3222Current address: Airbus Defence and Space, 31 rue des Cosmonautes,
31402 Toulouse, France, Sandrine Pires4, Serge Reynaud5, Laura Serron2, Jean-
Michel Travert6 1 DPHY, ONERA, Université Paris Saclay, F-92322 Châtillon,
France  2 Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, IRD,
Géoazur, 250 avenue Albert Einstein, F-06560 Valbonne, France  3 CNES
Toulouse, 18 avenue Edouard Belin - 31401 Toulouse Cedex 9, France  4
Laboratoire AIM, CEA/DSM–CNRS, Université Paris Diderot, IRFU/SAp, CEA Saclay,
Orme des Merisiers, 91191 Gif-sur-Yvette, France  5 Laboratoire Kastler
Brossel, UPMC-Sorbonne Université, CNRS, ENS-PSL University, Collège de
France, 75252 Paris, France  6 Altran technologies, 17 Avenue Didier Daurat,
31700 Blagnac, France<EMAIL_ADDRESS>
###### Abstract

After performing highly sensitive acceleration measurements during two years
of drag-free flight around the Earth, MICROSCOPE provided the best constraint
on the Weak Equivalence Principle (WEP) to date. Beside being a technological
challenge, this experiment required a specialised data analysis pipeline to
look for a potential small signal burried in the noise, possibly plagued by
instrumental defects, missing data and glitches. This paper describes the
frequency-domain iterative least-square technique that we developed for
MICROSCOPE. In particular, using numerical simulations, we prove that our
estimator is unbiased and provides correct error bars. This paper therefore
justifies the robustness of the WEP measurements given by MICROSCOPE.

††: Class. Quantum Grav.

  * December 2020

Keywords: Experimental Gravitation, Data Analysis

##  1 Introduction

MICROSCOPE’s test of the Weak Equivalence Principle (WEP) is based on the
comparison of the acceleration of two concentric cylindrical test masses of
different composition as they orbit the Earth [1, 2, 3]. Thus, the gravity
that pulls the test masses being sourced by the Earth, any Equivalence
Principle violation (EPV) signal will be proportional to the Earth gravity
acceleration. More precisely, as it is measured along the main axis of the
cylindrical masses, the signal we look for can be expected to be proportional
to the Earth’s gravity acceleration modulated by the motion and the attitude
of the satellite around the Earth, resulting in a periodic signal with a well-
known frequency (noted $f_{\rm EP}$ in this paper).

The main MICROSCOPE data consists in time series of accelerations measured by
two concentric accelerometers [4]. To look for an EPV, it is then enough to
look for a non-zero signal at the $f_{\rm EP}$ frequency in the difference of
those two time series. Obviously, this process is impacted by instrumental
noise and systematics (either instrumental or environmental). Thus, it
consists in correcting the measured time series from –calibrated and/or
modelled– instrumental and environmental systematics, before seeking a
possible periodic signal amounting to a violation of the WEP in coloured-
noise-dominated data.

As described in Ref. [5], the MICROSCOPE mission is divided in different
measurement sessions. Sessions represent a time span during which the
satellite and the instrument keep the same configuration (spin, drag-free
control law…). Some of these sessions are directly devoted to the WEP test
(called “EP sessions” in this paper) while others (“calibration sessions”) are
used to calibrate or characterise the instrument. EP sessions are the longest,
most of them lasting 120 orbital periods (about 8 days), while calibration
sessions typically last a few orbits.

In this paper, we present MICROSCOPE’s data analysis pipeline. After recalling
the measurement equation in Sect. 2, we present our frequency-domain iterative
ordinary least squares (OLS) algorithm and its mathematical background in
Sect. 3. We then use “worst-case” simulations from a hybrid software-hardware
MICROSCOPE simulator to discuss the optimal way to correct for gaps in
MICROSCOPE data, in Sect. 4. Finally, in Sect. 5, we show with well-controlled
numerical simulations that our pipeline provides unbiased estimates and
reliable error bars in the presence of instrumental systematics. Combined with
previous works where we showed how we can successfully deal with missing data,
those results prove that our data analysis pipeline allows us to reliably
measure an EPV, if any. We conclude in Sect. 6. A detailed appendix discusses
uncertainty propagation. Note that this paper is only about the MICROSCOPE
data analysis methodology, and does not present real data. For the real data
analysis (using the methods described in this paper), see Refs [3, 6, 7, 8].

The terminology and notations used in this paper, as well as all observables,
are defined in Ref. [4]. In particular, we define the common-mode (resp.
differential-mode) of a given observable or parameter as the half-sum (resp.
half-difference) of this observable/parameter for both test masses,
$o^{(c,d)}=(o^{(1)}\pm{}o^{(2)})/2$, and we use the convention $\mathbf{o}$ to
denote a vector and $[\mathbf{o}]$ to denote a second-order tensor,
respectively. We note the time derivative of observable $o$ as
$\accentset{\mbox{\large.}}{o}$.

##  2 Measurement principle

###  2.1 Measurement equation

MICROSCOPE looks for an EPV by monitoring the difference in accelerations
potentially undergone by the two test masses of a differential accelerometer
[1, 2, 3, 4]. In an ideal case, the measurement equation is then
straightforward to establish, since the difference of acceleration of the two
test masses is expected to be proportional to Earth gravity field
${\mathbf{g}}$,
${\mathbf{\Gamma}^{(d)}}\equiv{\mathbf{\Gamma}^{(1)}}-{\mathbf{\Gamma}^{(2)}}=\delta{\mathbf{g}}$,
where $\delta$ is the (approximate) Eötvös parameter that we wish to estimate.
However, the instrument is not perfect: for instance, scale factors are not
exactly unit and test masses are imperfectly centered and aligned with respect
to each other and to the satellite’s frame of reference, such that the common-
mode and differential-mode sensitivity matrices $[\mathbf{a_{c}}]$ and
$[\mathbf{a_{d}}]$ are not the identity and null matrices, respectively.
Furthermore, since the satellite rotates it imparts a Coriolis acceleration on
the test masses.

The measurement equation is therefore much more complicated than in the ideal
case. Ref. [4] establishes it in full generality, introducing and taking into
account all instrumental defects and their notation. We measure the
differential acceleration along the test masses’ sensitive axis
($x$-coordinate in the instrument’s frame, see Fig. 1 of Ref. [9]), such that

$\Gamma_{x}^{(d)}=2\tilde{b}_{x}^{(d)}+a_{c11}{\delta}{g_{x}}+a_{c12}{\delta}{g_{y}}+a_{c13}{\delta}{g_{z}}+{\Delta^{\prime}_{x}}{{S_{xx}}}+{\Delta^{\prime}_{y}}{{S_{xy}}}+{\Delta^{\prime}_{z}}{{S_{xz}}}\\\
+{\left(a_{c13}\Delta^{\prime}_{y}+a_{c12}\Delta^{\prime}_{z}\right)}{{S_{yz}}}+{a_{c12}\Delta^{\prime}_{y}}{{S_{yy}}}+{a_{c13}\Delta^{\prime}_{z}}{{S_{zz}}}\\\
+\left({-a_{c13}\Delta^{\prime}_{y}+a_{c12}\Delta^{\prime}_{z}}+2c_{d11}\right){\accentset{\mbox{\large.}}{\Omega}_{x}}-\left({\Delta^{\prime}_{z}-2a_{c13}\Delta^{\prime}_{x}}-2c_{d12}\right){\accentset{\mbox{\large.}}{\Omega}_{y}}\\\
+\left({\Delta^{\prime}_{y}-2a_{c12}\Delta^{\prime}_{x}}+2c_{d13}\right){\accentset{\mbox{\large.}}{\Omega}_{z}}+2\left({a_{d11}}{\tilde{\Gamma}^{(c)}_{x}}+{a_{d12}}{\tilde{\Gamma}^{(c)}_{y}}+{a_{d13}}{\tilde{\Gamma}^{(c)}_{z}}\right)\\\
+2\accentset{\mbox{\large.}}{\Delta}^{\prime}_{x}\Omega_{x}-2\accentset{\mbox{\large.}}{\Delta}^{\prime}_{z}\Omega_{y}+2\accentset{\mbox{\large.}}{\Delta}^{\prime}_{y}\Omega_{z}-a_{c11}\accentset{\mbox{\large..}}{\Delta}_{x}-a_{c12}\accentset{\mbox{\large..}}{\Delta}_{y}-a_{c13}\accentset{\mbox{\large..}}{\Delta}_{z}\\\
+{\tilde{K}^{(1)}_{2xx}}\left({\tilde{\Gamma}^{(1)}_{x}}\right)^{2}-{\tilde{K}^{(2)}_{2xx}}\left({\tilde{\Gamma}^{(2)}_{x}}\right)^{2}+2n_{x}^{(d)},$
(1)

where $\left[{\mathbf{T}}\right]$ the Earth gravity gradient tensor (GGT) in
the instrument’s frame,
$\left[{\mathbf{In}}\right]=\left[\accentset{\mbox{\large.}}{\bm{\Omega}}\right]+\left[\bm{\Omega}\right]\left[\bm{\Omega}\right]$
the gradient of inertia tensor, $\left[{\mathbf{S}}\right]$ the symmetric part
of the $\left[{\mathbf{T}}\right]-\left[{\mathbf{In}}\right]$ matrix,
$\left[\bm{\Omega}\right]$ the angular velocity tensor of the satellite, and
$\mathbf{\Delta}$ the vector between the center of the two test masses (called
“offcentering vector” hereafter).

In Eq. (1),
$\tilde{\mathbf{\Gamma}}^{(c)}=\mathbf{\Gamma}^{(c)}-\mathbf{n}^{(c)}$ is the
noise-free common-mode measured acceleration,
$\left[a_{c11},a_{c12},a_{c13}\right]$ is the first row of the common-mode
sensitivity matrix, $\left[a_{d11},a_{d12},a_{d13}\right]$ is the first row of
the differential-mode sensitivity matrix,
$\left[c_{d11},c_{d12},c_{d13}\right]$ is the first row of the differential-
mode sensitivity matrix to the angular acceleration,
$\tilde{b}_{x}^{(d)}=b_{0x}^{d}+a_{c11}b_{1x}^{(d)}+a_{c12}b_{1y}^{(d)}+a_{c13}b_{1z}^{(d)}$
(with $\mathbf{b_{0}}^{(d)}$ the differential electrostatic bias and
$\mathbf{b_{1}}^{(d)}$ the difference of mechanical perturbations acting on
the two test masses). For convenience, in the remainder of this paper, we
denote $\delta_{x}\equiv a_{c11}\delta$, $\delta_{y}\equiv a_{c12}\delta$,
$\delta_{z}\equiv a_{c13}\delta$, but we warn the reader that $\delta_{i}$
should not be confused with the component of a tensor.

We should stress that the measurement is not directly sensitive to the actual
offcenterings $\mathbf{\Delta}$, but to the following combinations of
instrumental parameters:

$\displaystyle\Delta^{\prime}_{x}$ $\displaystyle\approx$ $\displaystyle
a_{c11}\Delta_{x}-a_{c12}\Delta_{y}-a_{c13}\Delta_{z}$ (2)
$\displaystyle\Delta^{\prime}_{y}$ $\displaystyle\approx$ $\displaystyle
a_{c11}\Delta_{y}+2a_{c12}\Delta_{x}-a_{c23}\Delta_{z}$ (3)
$\displaystyle\Delta^{\prime}_{z}$ $\displaystyle\approx$ $\displaystyle
a_{c11}\Delta_{z}+2a_{c13}\Delta_{x}+a_{c23}\Delta_{y}.$ (4)

Only those “derived” offcenterings $\Delta_{x,y,z}$ can be estimated. Finally,
note that the quadratic factors $K_{2xx}^{(j)}$ were neglected in Ref. [4]
since they are estimated to be negligible in the real MICROSCOPE data. In this
paper, we aim to be as general as possible, so that we take them into account,
and set them to non-negligible values in the simulations below. Similarly,
Ref. [4] ignored the motion of the test masses (thus, setting the time
derivatives of $\mathbf{\Delta}$ to 0); here, we consider it since it should
be taken into account in some calibration sessions.

###  2.2 Modelled, estimated and ignored (combinations of) parameters

Each term of Eq. (1) is the product of a time-varying signal by parameters
related to the experiment (called “instrumental parameters” in the remainder
of this paper). Their effects are clear from Eq. (1): the differential bias
shifts the measured acceleration from the true acceleration; the common-mode
sensitivity matrix mixes components from all axes; the differential-mode
sensitivity matrix projects common-mode acceleration into the measured
differential acceleration; the coupling with the angular acceleration makes
the measurement depend on the satellite’s motion; and ${n}^{(d)}$ adds
stochastic instrumental noise to the measurement.

The varying signals are either directly measured or are derived by analyzing
tracking data, as described below and shown in Refs. [10, 6]: several
satellite manoeuvres (summarised in A) allow us to discriminate between
instrumental defects and estimate them separately.

Those terms come in four distinct categories and can be either corrected for,
or ignored in the measurement equation (1):

  1. 1.

Terms computed through models and/or on-ground restitutions:

     * •

The position of the satellite as determined in the celestial reference frame
J2000 [11] is expressed in the Earth’s reference frame using the appropriate
rotations [12]. From this position, the components of the gravity acceleration
$\mathbf{g}$ and of the GGT $\left[{\mathbf{T}}\right]$ are computed in the
Earth’s reference frame [13]; afterwards they are transformed to the J2000
celestial frame and then to the instrument frame using the attitude of the
satellite [11].

     * •

Like the attitude, the angular velocity ${\mathbf{\Omega}}$ and acceleration
$\accentset{\mbox{\large.}}{\mathbf{\Omega}}$ result from an on-ground
restitution [11].

     * •

$\tilde{\Gamma}^{(c)}_{x}$, $\tilde{\Gamma}^{(c)}_{y}$ and
$\tilde{\Gamma}^{(c)}_{z}$ are directly approximated by the measured
quantities $\Gamma^{(c)}_{x}$, $\Gamma^{(c)}_{y}$ and $\Gamma^{(c)}_{z}$.

  2. 2.

Terms corrected using results from calibration sessions
($\Delta^{\prime}_{y}$, $a_{d1j}$):

     * •

${\Delta^{\prime}_{y}}\left({{S_{xy}}}+\dot{\Omega}_{z}\right)$

     * •

$2\left({a_{d11}}{\tilde{\Gamma}^{(c)}_{x}}+{a_{d12}}{\tilde{\Gamma}^{(c)}_{y}}+{a_{d13}}{\tilde{\Gamma}^{(c)}_{z}}\right)$

  3. 3.

Terms dependent on parameters estimated with the current session:

     * •

$2\tilde{b}_{x}^{(d)}$: in principle $\tilde{b}_{x}^{(d)}$ is a readily
measurable constant, although it can drift in time (section 3.3).

     * •

$\delta_{x}g_{x}+\delta_{z}g_{z}$: we take advantage from the fact that
$g_{x}$ and $g_{z}$ vary at the same frequency $f_{\rm EP}$ but in quadrature
to estimate both $\delta_{x}$ and $\delta_{z}$ almost without correlation.
Since $a_{c11}\simeq 1$ and $a_{c13}\ll{}1$, $\delta_{x}\approx\delta$
(estimating $\delta$ being MICROSCOPE’s main objective). Given the upper bound
on $\delta$ from previous experiments [14], $\delta_{z}=a_{c13}\delta$ is in
principle not observable but a statistic over its value estimated in different
sessions could give interesting indications on the quality of the experiment.

     * •

$\Delta^{\prime}_{x}S_{xx}+\Delta^{\prime}_{z}S_{xz}$: the components
$\Delta^{\prime}_{x}$ and $\Delta^{\prime}_{z}$ of the offcentering can be
estimated very accurately thanks to the strong variations of $S_{xx}$ and
$S_{xz}$ at $2f_{\rm EP}$

  4. 4.

Terms neglected due to their very small magnitude [6]:

     * •

$\delta_{y}g_{y}$

     * •

$\left(a_{c13}\Delta^{\prime}_{y}+a_{c12}\Delta^{\prime}_{z}\right)S_{yz}+a_{c12}\Delta^{\prime}_{y}S_{yy}+a_{c13}\Delta^{\prime}_{z}S_{zz}$

     * •

$\left(-a_{c13}\Delta^{\prime}_{y}+a_{c12}\Delta^{\prime}_{z}+2c_{d11}\right)\dot{\Omega}_{x}+\left(2a_{c13}\Delta^{\prime}_{x}+2c_{d12}\right)\dot{\Omega}_{y}+\left(-2a_{c12}\Delta^{\prime}_{x}+2c_{d13}\right)\dot{\Omega}_{z}$

###  2.3 Sessions dedicated to the test of the EP

The measurement equation (1) is valid for all scientific sessions (calibration
and EP sessions). However, it can be simplified for EP sessions where test
masses are kept motionless, so that the velocities
$\accentset{\mbox{\large.}}{\Delta}_{x}$,
$\accentset{\mbox{\large.}}{\Delta}_{y}$,
$\accentset{\mbox{\large.}}{\Delta}_{z}$ and the accelerations
$\accentset{\mbox{\large..}}{\Delta}_{x}$,
$\accentset{\mbox{\large..}}{\Delta}_{y}$,
$\accentset{\mbox{\large..}}{\Delta}_{z}$ vanish. Further applying the
corrections described in Sect. 2.2, Eq. (1) then simplifies to

$\Gamma^{(d)}_{x,{\rm
corr}}=2\tilde{b}_{x}^{{}^{\prime}(d)}+\delta_{x}g_{x}+\delta_{z}g_{z}+\Delta^{\prime}_{x}S_{xx}+\Delta^{\prime}_{z}S_{xz}+2n_{x}^{(d)},$
(5)

which is the core model fitted to the data after applying the calibration
parameters: in addition to the Eötvös parameter $\delta_{x}$ we also estimate
$\delta_{z}$ which quantifies the amplitude of a signal proportional to
$g_{z}$ (varying also at the $f_{\rm EP}$ frequency but in quadrature with
$g_{x}$) and the components $\Delta^{\prime}_{x}$ and $\Delta^{\prime}_{z}$ of
the apparent offcentering.

The next section provides more conceptual and mathematical background on the
(instrumental and Eötvös) parameters used in the MICROSCOPE data analysis.

##  3 Parameters estimation

###  3.1 Iterative weighted least square fit

Each in-flight calibration session is dedicated to estimating one (or two)
parameters and designed so that the signals sourced by those parameters have a
favourable signal-to-noise ratio. Although it is theoretically possible to
cumulate all calibration sessions and estimate all parameters simultaneously
from Eq. (1), we use the fact that they are almost independent from each other
to simplify and better control the estimation process via an alternative
method. We thus devised a technique to estimate each parameter iteratively,
refining and updating the estimation of a given parameter using the estimation
of the other parameters until some convergence criterium is reached.

The measurement equation Eq. (1) is of the form
$\Gamma_{x}^{(d)}=f(p_{k},t)+n_{x}^{(d)}$, where $p_{k}$ are parameters and
the time dependence is linked to measured or modeled signals $s_{i}(t)$. For
each session, the data provides us with $\Gamma_{x}^{(d)}$ and all $s_{i}(t)$.
It is then possible to perform a least-square (or similar) fit to estimate the
parameters $p_{k}$ from a given model.

Moreover, for a given calibration session, we have a priori values $p_{k,0}$
for the parameters $p_{k}$, as some of them have been measured on ground, or
as others may have been estimated during an earlier in-flight calibration
session (Sect. 2.2). It is then possible to correct the measurement for the
corresponding signals, and use an updated version of the measurement equation,

$\Gamma^{(d)}_{x,{\rm corr}}(t)=\Gamma_{x}^{(d)}(t)-f(p_{k,0},t).$ (6)

This equation can finally be used to refine the estimation of some parameters
$p_{ke}$, with a least-square method using the model

$\Gamma^{(d)}_{x,{\rm corr}}(t)=\Sigma_{ke}\frac{\partial f(p_{k},t)}{\partial
p_{ke}}\left(p_{ke}-p_{ke,0}\right).$ (7)

The actual estimation of a parameter depends on the technique used to deal
with missing data (see Sect. 4 for the introduction of the techniques we use
to deal with them –KARMA [15], M-ECM [16] and inpainting [17, 18]). M-ECM
estimates the parameters and the noise and deals with missing data all at
once; however, inpainting only fills in missing data, and must be augmented by
a least square estimate.

We use iteratively the Adam (Accelerometric Data Analysis for MICROSCOPE)
Fortran code to estimate parameters in the frequency domain. Once this
iterative process has converged (typically in two to three iterations), we use
Eq. (5) to measure the Eötvös parameter $\delta_{x}$ on calibrated data. The
procedure is summarised in Algorithm 1. The data analysis processes underlying
Adam are described below.

Algorithm 1 Iterative least-squares estimation of $N$ instrumental parameters
from $M$ independent calibration sessions.

Initial prior:
$\Pi=\\{\pi_{0}(p_{1}),\dots,\pi_{0}(p_{i}),\dots\pi_{0}(p_{N})\\}$

while not converged do

for $i=1$ to $M$ do

Correct measurement from $i$th session with priors $\Pi_{i-1}$ (Eq. 6)

Estimate $i$th parameter (Eq. 7): estimator $\hat{p}_{i}$

end for

Update prior: $\Pi=\\{\hat{p}_{1},\dots,\hat{p}_{i},\dots\hat{p}_{N}\\}$

end while

Estimate the Eötvös parameter $\delta_{x}$ on calibrated data (Eq. 5)

###  3.2 Frequency domain least square analysis: Adam

To allow for a least-squares analysis [19], the corrected measurement equation
(7) can be formally written

${\mathbf{Y}}={\left[{\mathbf{A}}\right]}{\bm{\theta}}+{\mathbf{n}},$ (8)

where ${\mathbf{Y}}$ is the vector of $N$ measurements, ${\bm{\theta}}$ is a
vector of $q$ unknown parameters to estimate (e.g. Eötvös parameter or test
masses’ offcentering), $\left[{\mathbf{A}}\right]$ is the design matrix and
${\mathbf{n}}$ is the noise vector. Since the models (given by Eq. 1 in the
most general case and by Eq. 5 for EP sessions) are linear with respect to the
estimated parameters, the columns of $\left[{\mathbf{A}}\right]$ simply
correspond to the signal associated to each parameter, sampled at the epochs
of the measured acceleration. The $N$ measurements are assumed to be regularly
sampled at a frequency $f_{\rm e}$ over a duration $T$. In case of missing
data the analysis takes place after reconstruction of these data using
inpainting or M-ECM algorithms (Sect. 4).

####  3.2.1 Transformation of the problem from the time domain to the
frequency domain

In order to solve the problem in the Fourier domain, we take the Fourier
transform of Eq. (8). To this aim, we make use of the Discrete Fourier
Transform Operator $\left[{\mathbf{F}}\right]$. The DFT operator being
unitary, the signal energy content is preserved by the transformation.

The new system can be simply written

$\hat{\mathbf{Y}}={\left[\hat{\mathbf{A}}\right]}{\bm{\theta}}+\hat{\mathbf{n}}.$
(9)

The original quantities being real, the new system can be reduced to $N$
useful real equations. These new equations can be grouped by pair (related to
real and imaginary parts of the DFT), corresponding to frequencies
$f_{k}=\frac{k}{T},k=1\cdots\lfloor\frac{N-1}{2}\rfloor$.

Moreover, the discrete Fourier transform drastically decreases the
correlations between the measurements projected in the frequency space (Fig.
1): the covariance matrix associated to
$\left[{\mathbf{F}}\right]{\mathbf{n}}$ is diagonal dominant. This is
beneficial as a diagonal weight matrix leads to a quasi-optimal solution.

Figure 1: Covariance matrix of the measured differential acceleration of an EP
session. Working in the frequency domain allows us to deal with a diagonal
covariance matrix.

####  3.2.2 Weighting in the frequency domain

Since each measurement projected in the Fourier domain can be associated to a
discrete frequency, the corresponding weight is

$w(f_{k})=\frac{1}{\sqrt{\gamma(f_{k})}},$ (10)

where $\gamma(f_{k})$ is the Power Spectral Density (PSD) of the noise at the
frequency $f_{k}$. In practice, the PSD is estimated by smoothing the residual
noise in the frequency domain resulting from a first estimation and removal of
the signal. The estimation–correction–smoothing–weighting process can be
iterated until convergence; in practice, two iterations are sufficient.

####  3.2.3 Restriction to the frequency bands containing the dominant signals

As shown in Ref. [4], the MICROSCOPE mission was designed to concentrate
useful signal on specific frequencies (i.e., a potential EPV signal peaks at
$f_{\rm EP}$, the GGT signal at $2f_{\rm EP}$ and calibration signals at
$f_{\rm cal}$). This is so true in the real data that a very simple analysis
such as synchronous detection could lead to reasonable results. However, we
use a more flexible method: we limit our least square inversion to the bands
of frequency containing the relevant signals. In practice, this is equivalent
to extracting a subsystem of Eq. (9) by selecting the relevant equations to
get the truncated system

${\left[\hat{\mathbf{A}}_{r}\right]}{\bm{\theta}}+\hat{\mathbf{n}}=\hat{\mathbf{Y}}_{r}.$
(11)

This trade-off between synchronous detection and the inversion of the full
system brings several advantages:

  * •

it is more robust than synchronous detection in case of small fluctuation of
the frequencies of the signals; it is not even necessary to know precisely the
value of these frequencies;

  * •

the choice of large enough bands containing also a substantial sample of noise
allows us to compute consistent values of the goodness of fit;

  * •

contrary to the inversion of the full system (9) the solution of the truncated
system (11) is immune to possible unmodelled perturbations in frequency bands
containing no useful information (especially high frequency bands). Taking
these bands into account would not change the parameter estimations (because
equations at different Fourier frequencies are uncorrelated) but could
decrease the global goodness of the fit;

  * •

the number of observation equations is neatly decreased and the analyses are
faster.

A final clarification is in order: the very low frequency part of the signal,
which contains in particular the zero frequency component, is not used in the
analysis. Indeed, not only does the zero frequency contain well-known signal
such as some components of the gravity gradient but it also contains unknown
signals such as the bias of the accelerometers. All these contributions are
not separable, so that the zero frequency is difficult to exploit.

###  3.3 Bias drift

Strictly speaking, the bias $\vv{b_{0}}^{(d)}$ included in Eq. (1) should be
constant. But in practice a very low frequency evolution can appear, in
particular because of thermal variations in the instrument (left panel of Fig.
2). This evolution can be efficiently corrected by fitting and subtracting a
polynomial (right panel of Fig. 2). As illustrated by Fig. 3, the polynomial
affects all frequencies in a wide low frequency band; consequently, the
correlation with the EPV signal (which is concentrated at the $f_{\rm EP}$
frequency) is very low. But to be sure that fitting a polynomial does not
affect the signal around $f_{\rm EP}$, we take advantage of our estimation in
the Fourier space to remove the contribution of these frequencies in the
estimation of the polynomial. More precisely:

  1. 1.

we first estimate simultaneously (i) the Eötvös parameter using a narrow
frequency band around $f_{\rm EP}$, (ii) the components $\Delta^{\prime}_{x}$
and $\Delta^{\prime}_{z}$ of the offcentring using a narrow frequency band
around $2f_{\rm EP}$, and (iii) the coefficients of the polynomial using a
wide band from 0 to 1 Hz but excluding the $f_{\rm EP}$ and $2f_{\rm EP}$
bands.

  2. 2.

we then use the first estimations of all these parameters to correct the
measurements according to Eq. (6); since this correction is performed in the
time domain, we implicitly correct all the frequencies and in particular the
polynomial is also applied to the $f_{\rm EP}$ and $2f_{\rm EP}$ bands.

  3. 3.

finally, we use the corrected acceleration to re-estimate the same parameters
in the same way as in step (i).

Another option is to correct low frequency variations with a model of
temperature sensitivity. Ref. [8] shows that the substraction of a model
linear with respect to the temperatures produces results similar to those
obtained with a polynomial fit.

Figure 2: Time evolution of the differential acceleration. Left: evidence for
a polynomial-like tendency. Right: after correcting for a polynomial of degree
3. Figure 3: Same as Fig. 2 but in the frequency domain. As expected, the
subtraction of the polynomial (blue) clearly reduces the low frequency
contribution.

###  3.4 Combining sessions

In order to decrease the stochastic error in the estimation of (instrumental
and Eötvös) parameters, it is interesting to combine measurements from several
sessions with the same configuration. This combination can be done either in
the time or in the frequency domain, as shown below.

####  3.4.1 Time domain

Suppose that session $k$ leads to the following linear system of $N_{k}$
equations (of the form of Eq. 8):

${\left[{\mathbf{A}}_{k}\right]}{\bm{\theta}}+{\mathbf{n}}_{k}={\mathbf{Y}}_{k}$
(12)

To carry on with a classical inversion in the time domain, one can just gather
the above matrices and vectors:

$\left[{\mathbf{A}}\right]=\begin{bmatrix}{\left[{\mathbf{A}}_{1}\right]}\\\
{\left[{\mathbf{A}}_{2}\right]}\\\ \vdots\\\
{\left[{\mathbf{A}}_{m}\right]}\end{bmatrix},\quad{\mathbf{n}}=\begin{bmatrix}{\mathbf{n}}_{1}\\\
{\mathbf{n}}_{2}\\\ \vdots\\\
{\mathbf{n}}_{m}\end{bmatrix},\quad{\mathbf{Y}}=\begin{bmatrix}{\mathbf{Y}}_{1}\\\
{\mathbf{Y}}_{2}\\\ \vdots\\\ {\mathbf{Y}}_{m}\end{bmatrix}$ (13)

where $m$ is the number of sessions considered. Then an appropriately weighted
least-squares technique could be used to solve for the concatenated system.
Given the advantages with doing an analysis in the frequency domain as
outlined above, we have not used this time domain analysis.

####  3.4.2 Fourier domain

First, DFT are applied separately to the systems corresponding to each session
$k$, as described in Sect. 3.2.1, to get equivalent systems in the frequency
domain:

${\left[\hat{\mathbf{A}}_{k}\right]}{\bm{\theta}}+\hat{\mathbf{n}}_{k}=\hat{\mathbf{Y}}_{k}.$
(14)

We use a weighting based on the PSD, as described in Sect. 3.2.2. Finally, all
considered weighted systems are gathered as shown above for the time domain
case.

Like in the case of a unique session, we can use a Generalised Least-Squares
regression to solve either the whole system or, more efficiently, a selected
subset of equations corresponding to the desired bands of frequency.
Nevertheless, the estimation process is flexible enough to allow us to take
into account the drift of some parameters (e.g. offcenterings, which vary with
respect to temperature).

###  3.5 Error propagation

Uncertainties in estimated parameters (or in the a priori knowledge of
parameters that cannot be finely estimated) propagate when correcting for them
in the measured differential acceleration, which involves three effects when
estimating the EPV signal: (1) residual systematic errors coming from the bias
of each parameter estimate; (2) statistical errors coming from the variance of
each parameter estimate and (3) statistical errors coming from the fact that
we use noisy signals to calibrate the measurement, whose noise is coupled with
the parameters. The first type of errors is assessed by upper bounds given by
a detailed performance analysis [20, 6]. The second and third type of errors
are discussed in this section: they can be assessed through Monte Carlo
simulations or directly included in the least-square covariance matrix.

Using an iterative estimation technique involves propagating the uncertainties
of the $i$th iteration estimates to the $(i+1)$th iteration. Consequently, as
can be seen from Eq. (6), the variance of the calibrated differential
acceleration is not only that of the raw measured differential acceleration,
but also involves a contribution from the (imperfect) knowledge of the
corrected instrumental parameters.

Error propagation is derived rigorously in B. Eq. (20) gives the expected
value of the calibrated differential acceleration used to estimate the Eötvös
parameter when assuming that the estimates of all instrumental parameters are
unbiased. Its variance is given by Eq. (23). Those equations can easily be
generalised to the case where a given instrumental parameter is estimated
given some others.

Two effects coming from uncertainties in the iterative correction of
instrumental parameters (through errors in the estimated amplitudes of the
signals that it contains) can be noted: (i) statistical errors in the
calibrated differential acceleration, because of the imprecise estimation of
corrected parameters (during previous iterations), and (ii) uncertainty on the
calibrated differential acceleration residuals (due to the fact that we
subtract noisy quantities multiplied by estimated parameters). Therefore, we
can expect that the least-square estimator of a parameter measured from the
calibrated differential acceleration will be affected by statistical errors
from instrumental parameters, as well as by statistical errors coming from the
measurement noise, and errors on statistical errors coming from the
uncertainties on the noise.

In practice, in order to rigorously propagate uncertainties in our iterative
least-square estimates, Eq. (23) can be used as the input variance, for each
parameter, at each iteration. Another possibility is to run Monte Carlo
simulations that use the variance as measured from the noise of several
calibrated differential accelerations created from drawings of instrumental
parameters (that can be assumed normally distributed, with mean and variance
defined by the previous estimations). Although more CPU- and time-consuming,
this technique is more sound than using Eq. (23), since it does not involve
“educated guesses” on some parameters (note that some terms appearing in Eq.
(23) are the “true” values and not their estimates).

We show below that given the smallness of the uncertainty on the estimated
instrumental parameters, we can actually safely ignore propagating errors on
the estimation of instrumental parameters.

##  4 Glitches and missing data

###  4.1 (Masked) glitches and missing data (in practice)

Data “gaps” come in two flavours: (i) glitches, where data is available but
contaminated by physical (e.g. impact with a micrometeorite or a satellite
crack) or measurement (e.g. an internal saturation in the instrument’s servo-
loop command) processes uncorrelated with testing the WEP, and (ii) telemetry
losses, where data is missing during some time intervals. Glitches may bias
the measurement or introduce spurious signals at specific frequencies, or
within some frequency range [7]. Similar events are present e.g. in LIGO/Virgo
data, where they hamper the detection and characterisation of candidate
gravitational waves signals [21], and several techniques have been developed
to correct for them (see e.g. Refs [22, 23, 24, 25, 26]).

To counteract the direct effect of glitches on MICROSCOPE’s WEP measurement,
we can model and subtract them as discussed in Ref. [27]; here, we choose to
remove glitches, hence giving them the same status as telemetry losses
(“missing” data). Missing data break the even sampling of data in time,
thereby preventing us from using standard OLS in the frequency domain and,
most critically creating an important spectral leakage due to noise colour,
potentially burying the EPV signal in MICROSCOPE data [15, 17].

We have developed and adapted methods to cope with such missing data, that
allow us to reach the required accuracy on the EPV measurement:

  * •

KARMA [15]: the Kalman Auto-Regressive Model Analysis is a generalized least-
square technique based on an autoregressive model of the (unknown)
instrumental noise, whitened with a Kalman filter.

  * •

M-ECM [16]: the Modified-Expectation-Conditional-Maximization technique allows
us to maximize the likelihood of available data through the estimation of
missing data by their conditional expectation, based on the circulant
approximation of the complete data covariance.

  * •

inpainting [17, 18]: originally developed for 2D observational cosmology [28]
and 1D asteroseismology [29], the inpainting technique uses a sparsity prior
to estimate the most probable value of missing data, therefore allowing us to
use an ordinary least-square fit (Adam –see below).

Those techniques have already been tested on numerical simulations with the
simplifying assumptions that instrumental defects were perfectly known and
corrected for [15, 16, 17, 18]. Note that Eq. (5) is equivalent to the form
used by Refs. [15, 16]; when further correcting for offcenterings, it gives
the measurement equation used by Refs. [17, 18]. Ref. [1] used Adam with no
correction of transients nor missing data (taking advantage of quiet and clean
measurement sessions). Ref. [2] uses Adam on data corrected for missing data
with inpainting.

Figure 4: Example of a transient and its correction. The black curve shows the
raw measurement. The gray zone is the region discarded by the masking
procedure. The red curve is the acceleration corrected by inpainting.

To search for glitches, we use a standard recursive $\sigma$-clipping
technique, in which we flag as outliers every point that deviates from the
moving average of the data by more than a specific number of times (4 when
defining masks on the differential-mode acceleration, or 5 when masks are
defined on the common-mode acceleration) the standard deviation of the data
[7]. We then mask a specific time interval after each detected outlier to make
sure that the transient regime is always removed. The gray area in Fig. 4
shows such a mask after an outlier to remove the corresponding glitch, visible
in black, in SUSON-simulated data (see below). Additional data points are
flagged by the instrument’s electronics, if an internal saturation has been
detected by the accelerometer digital electronics (and eventually smoothed
out, thereby invisible in the acceleration data [9]): we decide to reject
those points. The red curve in Fig. 4 shows the inpainting correction to the
masked glitch (black).

We discuss below (Sect. 4.3) how missing data couple with MICROSCOPE
instrumental parameters. This allows us to define an optimal way to deal with
missing data within the iterative process described in Sect. 3, by considering
at which level of the data processing the correction for missing data should
appear. This discussion is based on numerical simulations created by a hybrid
software-hardware simulator developed by CNES (Centre National d’Etudes
Spatiales) to validate MICROSCOPE’s ground segment operations before flight.
Those simulations correspond to worst-case scenarios, with exaggerated number
of glitches and maximum allowed values for instrumental parameters. Although
those simulations consider instrumental defects higher than the in-flight
MICROSCOPE data, they are valuable to amplify the effect of instrumental
parameters and of missing data, and allow for a pedagogical discussion
hereafter (all figures in this section were produced using those simulations).

###  4.2 Hybrid simulations: SUSON

The mock data used in this section were produced by the BVSS (Banc de
Validation Système et SCAO –Drag Free and Attitude Control System Validation
Bench), a CNES real time test bench which serves different purposes in project
development plans. By simulating very accurately the spacecraft behavior, it
aims at validating all the onboard functional chains, in particular the
attitude and acceleration control. It also aims at performing all the
technical and operational system tests while connected to the ground entities
(operational and scientific centers).

The MICROSCOPE test bench relies on the other CNES microsatellite simulation
tools, and shares many hardware and software commonalities, which gives good
confidence in the product. Three benches were developed to enable the various
and intensive tests performed by the project. Each bench is constructed as
follows. A Sun Microsystems station, running at 64Hz, includes various models
of the spacecraft environment (Earth gravitational potential, Earth atmosphere
and albedo, Sun illumination) and of its dynamics. It also includes physical
models of all the embarked equipment (actuators and sensors), with major
evolutions for the MICROSCOPE peculiarities (see below). Those equipment
models are coupled with a rack of interface electronic boards, emulating the
real hardware communication protocol for each equipment. A harness then
connects the interface boards to a real onboard computer, on which the onboard
software is running. Finally, electrical ground support equipment simulates
the ground stations and the radiofrequency link between stations and
spacecraft, so that the users receive the telemetry data output by the onboard
computer as in the real life.

A peculiarity of the MICROSCOPE bench is the T-SAGE333T-SAGE (Twin Space
Accelerometers for Gravitation Experiments) is MICROSCOPE’s instrument core
[9] simulation. A very representative model was developed by CNES, called
SUSON. It is composed of a real Interface Control Unit (payload computer),
running at 1027Hz and interfaced on one side with the main onboard computer
and on the other side with a specific board simulating both the Front End
Electronic Unit and the Sensing Unit [9]. The latter is fed with 64Hz data
coming from the dynamics models of the spacecraft, so as to inject outside
acceleration and gravity as sensed by the test masses. Within SUSON, the
1027Hz control loop of two tests masses is simulated with the real T-SAGE
computer and software in the loop, down to the electrode level. Many
parameters enable to tune the characteristic of the measurements: mass
misalignments, mass centering, scale factor, quadratic scale factor,
stiffness, and biases. A theoretical spectrum of T-SAGE noise is also
injected.

Another interesting aspect of this test bench is the simulation of micro-
perturbations. A complementary model was developed to inject random impacts
within the dynamics model, with the possibility to tune the time, direction,
and momentum random distributions. The lack of flight data at the time of the
first simulations campaign imposed a conservative parametrization444When
compared to flight data (MICROSCOPE was not launch at the time the simulations
were created, but at the time of writing, we now have in-orbit estimates of
the amount of crackles), it shows that the real level of micro-perturbation is
far below hypothesis..

Finally, we should emphasize that it was the first time a BVSS test bench was
pushed so far in the simulation process. Among its notable characteristics is
the wide range of timecales used: 1027Hz for the payload frequency, 64Hz for
dynamics simulation, 4Hz for the onboard activation of the SCAA (Attitude and
Orbit Control System) tasks, the whole process running during 10 days to
obtain the final simulation outputs. Such a long run time is required to
simulate 120-orbit sessions during which the hardware elements of the
simulator ensure a very realistic simulation of the spacecraft and the
possibility to validate in depth the functional chains.

Figure 5: Spectra of the simulated differential acceleration for two different
measurement sessions, simulated with SUSON. Left: inertial WEP-measurement
session (shown are spectra with glitches present (gray) and masked then
replaced by inpainted signal (black)). Right: $\Delta_{y}$ calibration
session.

Fig. 5 shows the spectra of the simulated differential acceleration for two
different measurement sessions, simulated with SUSON. The left panel shows the
noise and spectral leakage induced by glitches and missing data (gray) and
their correction with inpainting (black) in the frequency domain. The right
panel shows the differential acceleration during a simulated $\Delta_{y}$
calibration session. The bump between 0.01 Hz and 0.1 Hz is not caused by
transients but is due to worst-case inaccuracies simulated in the attitude
control coming from the star sensor [30], and is linked to the angular
velocity of the satellite. This bump is not prejudicial since it occurs at
frequencies different than those of the excitations and estimated signals.
Furthermore, it can be corrected for when estimating parameters in the
iterative way presented in Sect. 3.

###  4.3 Effect of the data level on glitches masking and missing data
correction

As presented in Ref. [5], we use different data levels during data processing
and analysis. Those relevant for our current discussion are:

  * •

N1a data: raw science data sorted by inertial sensor

  * •

N2a data: raw science differential- and common-mode accelerations

  * •

N2b data: derived from N2a data after correcting for instrumental parameters.
They are used to estimate the EPV signal.

We can wonder at which stage we should define the mask and correct for
glitches and missing data. We show here that the answer depends on what task
we aim to accomplish.

####  4.3.1 Instrumental parameters estimation

We deal with calibration at the N2a level; thus, two routes can be taken to
define masks: either we mask N1a accelerations, then logically add those masks
when creating N2a data, or we ignore masking N1a data and directly mask at the
N2a level. One could naively think that when creating the differential
acceleration, glitches would cancel out (both sensors see them in the same
way, since they come from external perturbations); then, the mask needed at
the N2a level should be less conservative than the logical addition of the N1a
masks. However, sensors are not exactly identical, therefore leaving
significant (though attenuated) imprints of glitches and other invalid data at
N2a level. We find that defining the mask either at N1a level or N2a level
lead to approximately the same final N2a mask, although we find that it is
best to apply inpainting on N2a data level rather than on the N1a level. When
applying inpainting to N1a data, reconstruction errors add up in the N2a level
data, thereby providing a less optimal result. We also find that KARMA and
M-ECM are more stable than inpainting with respect to the mask definition.

####  4.3.2 WEP measurement: non-commutativity of inpainting and systematics
subtraction

We look for an EPV at the N2b level. Once instrumental parameters are
estimated, subtracting their contribution from Eq. (1) can a priori be done on
data still plagued by missing or invalid data, to which simple subtractions
should be immune. We then have once again two routes to choose between in
order to go from N2a to N2b data: either correct for instrumental parameters
on already inpainted N2a data, or correct for them on raw N2a data, then apply
inpainting, KARMA or M-ECM on those raw N2b data level. KARMA and M-ECM do not
rely on filling missing data, but only on the mask definition. Therefore,
their precision does not depend on whether we define masks at the N2b level or
we create N2b level masks from N2a masks. Hence, the following discussion
pertains only to methods that fill data gaps and we specialise on the
inpainting case.

The instrumental parameters are coupled with the measured common-mode
acceleration (through the $a_{d}$ parameters) and the square of the
acceleration of individual sensors (through the quadratic factors $K_{2}$).
Their subtraction will then evidently depend on whether we subtract them (i)
before applying inpainting or (ii) after correcting them for invalid and
missing data. In the former case, we create raw N2b data from raw N2a data,
that we must eventually mask and inpaint; in the latter case, we create
inpainted N2b data from inpainted N2a data. The guiding principle to choose
between those two possibilities is to maximise the consistency between N2a
data (inpainted or raw) and systematic/instrumental effects.

We start by computing the systematic/instrumental parameters that are to be
subtracted from N2a data. Fig. 6 shows the different contributors to
systematic effects when the common-mode and individual accelerations are not
inpainted (left panel) and are inpainted (right panel). In those figures, the
blue curve shows the main effect of the Earth gravity gradient combined with
the offcentering of the test masses; the green curve shows the effect of the
common-mode misalignments $a_{cij}$; the purple curve shows the effect of the
test masses’ motion and satellite’s attitude; the orange curve shows the
effect of the projection of the common-mode accelerations; the cyan curve
shows the effect of the angular to linear couplings; the red curve shows the
effect of the quadratic factors; the transparent black curve is the total
systematics.

The main deterministic systematics come from the coupling between the test-
masses offcentering and the gravity gradient. The 2$f_{\rm EP}$ line can be
clearly seen, but the $f_{\rm EP}$ peak is strong enough to mimic an EPV
violation of about $10^{-14}$ m/s2 if not corrected for. Furthermore, the
spectral leakage from the mean of the quadratic factor dominates at almost all
frequencies. It is responsible for the flat spectrum at low frequency, at a
level of $10^{-13}$ m/s2 before inpainting and $10^{-15}$ m/s2 after
inpainting. This plateau can be seen clearly on the N2a differential
acceleration with the gray curve of the left panel of Fig. 5.

Figure 6: Systematics and instrumental parameters contributions on a SUSON
simulation. Left panel: with no correction of missing data. Right panel: after
correction of missing data. See main text for the description of each curve.

The two options to create N2b data corrected for invalid and missing data can
be seen graphically as (i) subtracting the black curve (at the data level) of
Fig. 6’s left panel to Fig. 5’s gray curve, then mask the resulting data and
apply inpainting to it or (ii) just subtracting the black curve (at the data
level) of Fig. 6’s right panel to Fig. 5’s black curve. The resulting N2b
spectra are shown in Fig. 7.

The differences are striking. On the one hand, subtracting raw systematics to
raw N2a differential accelerations, then mask and inpaint them (left panel of
Fig. 7), provides a very clean N2b differential acceleration, with no
significant deterministic signals. On the other hand, subtracting inpainted
systematics to inpainted N2a differential acceleration provides a poor N2b
acceleration (right panel of Fig. 7), with a higher low-frequency noise and
still plagued by deterministic signals, with peaks at $f_{\rm EP}$ and
2$f_{\rm EP}$ and the “attitude” bump about $10^{-2}$Hz (right panel of Fig.
5): although this bump is not problematic to estimate an EPV signal, the peak
at $f_{\rm EP}$ would lead us to incorrectly conclude for the existence of an
EPV.

Figure 7: Left panel: inpainted N2b differential acceleration created from raw
N2a differential acceleration and systematics. Right panel: N2b differential
acceleration created from inpainted N2a differential acceleration and
systematics.

This can be explained in the following way. Systematics and instrumental
parameters are coupled with the common-mode acceleration and with the
individual test mass squared accelerations; those are to be subtracted from
the (N2a) differential-mode acceleration. Although there seems to be four
different accelerations involved in this process, there are actually only two
(whose sum and difference make up the common-mode and differential-mode
accelerations), and hence all accelerations are significantly correlated (up
to the differences of both sensor’s transfer function). In particular,
transients are seen for both test masses (albeit not in exactly the same way
due to their different transfer functions). Therefore, when subtracting raw
systematics to raw N2a differential acceleration, we consistently remove most
of the invalid data, and we almost get clean N2b data (see the gray spectrum
in Fig. 8); masking the remaining outliers and inpainting them then improves
the spectrum (to obtain the clean N2b data of the black curve of Fig. 8).

In contrast, when masking and inpainting N2a data and N1a data (needed to
subtract the quadratic factors, that couple to the individual accelerations
–see Eq. 1), inpainting introduces small errors (although not visible on
individual spectra) that break the consistency between the four accelerations
involved. In particular, if some bias (that can be different for the
differential-mode and common-mode acceleration) is introduced by inpainting in
the inpainted deterministic signals, their subtraction will be imperfect,
resulting in an incorrect N2b differential acceleration. Such biases are the
likely cause of the remaining peaks in the spectrum of the lower panel of Fig.
7. Furthermore, if transients are periodic, those inconsistent errors will
create a periodic pattern at the N2b differential acceleration, that can
plague the inpainting reconstruction, thereby introducing an artefact at the
frequency of this pattern. Such an artefact can increase the appearance of a
sub-optimally subtracted deterministic systematics.

To summarise this discussion: the invalid/missing data correction (using
inpainting) and systematics substraction do not commute. Therefore, we first
create raw N2b-level differential accelerations from raw N2a accelerations and
raw systematics, before masking them and filling their gaps. This allows us to
apply inpainting only once (thus using it in the regime explored in [17, 18]
and minimising the bias observed in [18]), and to take advantage of the
correlations between all the accelerations involved in the process.

Figure 8: Spectrum of N2b level differential acceleration for the inertial
session used to estimate the EPV in SUSON simulations, before (grey) and after
(black) glitches and missing data correction. This spectrum can be compared
with that shown (in black) in Fig. 5 before correcting for instrumental
parameters.

##  5 Validating iterative least squares through numerical simulations

Although the SUSON simulations used above are sufficient to investigate
qualitatively the impact of instrumental parameters and the effect of missing
data, they are not precise enough to test our full pipeline down to the
precision required for MICROSCOPE. This is because they were not designed for
such a high-precision task, but rather to make sure that all ground-segment
systems worked as expected. Consequently, we developed Monte Carlo numerical
simulations dedicated to test the data analysis technique presented in Sect.
3, to show the internal consistency of our iterative least-square estimations
in the absence of data gaps. We already showed that KARMA, M-ECM and
inpainting allow us to correctly deal with missing data when estimating (in a
non-iterative way) parameters [15, 16, 17, 18], thereby allowing us to ignore
missing data in the remainder of the paper.

We assess the accuracy and precision of (i) ADAM’s least-squares estimator and
of (ii) the iterative calibration of instrumental parameters. We investigate
the first point on focusing on how well Adam recovers the
$\Delta^{\prime}_{x}$ offcentering input in the simulations, when all other
parameters are either ignored or perfectly corrected for. The second point is
investigated by simulating 200 sets of calibration sessions and running our
iterative least-square on them. In particular, such Monte Carlo simulations
allow us to easily look into error propagation. All the simulations use the
same instrumental parameters and noise characteristics, each of them with its
own noise realisation (simulated from the power-law noise power spectral
density discussed in Ref. [31]), so that the dispersion of Adam’s outputs are
solely statistical and representative of the data analysis process.

###  5.1 Simulations

We use Simula, a numerical fortran simulation code developed specifically for
MICROSCOPE, to simulate the acceleration measured by both test masses. The
simulator takes into account instrumental parameters (Eq. 1) the satellite’s
orbit and attitude (either measured from real data or simulated by any
external orbit simulator –e.g. GINS [32] for the orbit) and computes the Earth
GGT in the instrument’s frame. Finally, the motion of each test mass is
simulated individually along its orbit (see Ref. [10, 4] for the acceleration
of individual test masses). Measurement errors and further systematic effects
(such as, but not limited to, attitude errors, datation errors, drifts,
periodic perturbations, missing data, scale factor stability, stability of
coupling between axes) can be taken into account, but the electronic servo-
loop used to control the test masses cannot easily be simulated, and is
therefore ignored in this paper. Although we simulate the effect of
instrumental parameters (misalignments, scale factors, couplings…), we ignore
systematic effects not linked to the instrument itself (e.g. drag-free
residual or attitude control imperfections) as well as transients and missing
data. However, we take care to define all calibration sessions as they are
performed in reality (i.e. with the same duration and manoeuvres
characteristics, see Ref. [6]). Using those simplified simulations allows us
to check that the least-squares estimator is unbiased.

The values of the parameters that can be estimated are the same as in the
SUSON simulations of Sect. 4, and are listed in the second column of Table 1.
We use exaggerated values, so that the instrumental defects’ effect is most
easily visible. Besides instrumental parameters, we simulate an EPV
$\delta=10^{-14}$. Our goal is thus to recover all instrumental parameters,
and estimate the EPV with a $7\times{}10^{-15}$ precision (given the length
and noise level of the corresponding simulated session) after correcting them
from the simulated differential acceleration.

###  5.2 Results

####  5.2.1 Adam estimates’ bias and variance

We first check that the Adam least-square estimates are not intrinsically
biased. To this aim, we measure the $\Delta^{\prime}_{x}$ offcentering, first
with no prior on other instrumental parameters, i.e. without correcting for
them (which is equivalent to the first iteration of the iterative
calibration). The left panel of Fig. 9 shows the distribution of the
$\hat{\Delta}^{\prime}_{x}$ estimator, based on 2000 simulations; the black
dot and associated error bar show the mean and rms of the distribution; the
input offcentering (as should be recovered by Adam) is shown by the red dotted
vertical line on the right hand side of the figure. This estimation is clearly
biased; although the bias is small ($\approx$ 1%), it is highly significant.

The observed bias is due to other uncorrected instrumental parameters
affecting the measurement, which justifies the use of an iterative technique.

Figure 9: Distribution of $\hat{\Delta}^{\prime}_{x}$ estimated with adam,
from 2000 simulations. The red dotted vertical line is the value to be
recovered. The black dot and error bar show the mean and rms of the
distribution. Left: without correcting for other parameters. Right: with
perfect prior knowledge of other parameters

To make sure of this assertion, we now consider the estimation of
$\Delta^{\prime}_{x}$ when other instrumental parameters are perfectly known,
and can therefore be corrected for in the measured differential acceleration
before estimating $\Delta^{\prime}_{x}$. The right panel of Fig. 9 shows the
distribution of the $\hat{\Delta}^{\prime}_{x}$ estimators in this case:
adam’s estimate is now unbiased, proving that the bias observed above was only
due to ignoring the effect of other parameters. We can therefore claim that it
is essential to iterate across instrumental parameters when measuring them and
that the iterative Adam’s estimation is unbiased. Nevertheless, we will show
below that the Eötvös parameter estimation is insensitive to small biases on
instrumental parameters, and therefore does not require an iterative process.

We show the distribution of the estimators’ variance in Fig. 10. We can first
notice that adam’s estimated variance is distributed as a
$\chi^{2}$-distribution, as expected for a least-square estimator’s variance.
In this figure, the red dotted vertical line is the “true” variance, as
computed from the estimator’s distribution (Fig. 9). The black dot and
associated error bars give the mean and rms of the variance of Adam’s
individual estimator. Since the “true” variance is within those error bars, we
can conclude that the variance provided by Adam for each estimation is
reliable. This is particularly important in the MICROSCOPE landscape, and we
can safely conclude that the upper bounds statistical errors on the Eötvös
parameter provided in Refs. [1, 2] are correct.

We can then safely conclude that adam is intrinsically unbiased and provides
correct error bars. We checked that our conclusions about
$\Delta^{\prime}_{x}$ still hold for different values of
$\Delta^{\prime}_{x}$, as well as for other instrumental parameters.

Figure 10: Distribution of Var($\hat{\Delta}^{\prime}_{x}$) estimated with
adam. The red dotted vertical line is the “true” variance as computed from the
distribution of the estimator (Fig. 9). The black dot and error bar show the
mean and rms of the estimated variance distribution.

####  5.2.2 Full iterative calibration

Table 1: Results of the Monte Carlo simulations: measured value of
instrumental parameter at each iteration. For each parameter, the first line
provides the results when propagating errors from one iteration to the next;
the results listed in the second line ignore the error propagation. The second
column gives the input values with the precision required to meet MICROSCOPE’s
objective in brackets.

| Expected | Iter 1 | Iter 2 | Iter 3 | Iter 4 | Iter 5
---|---|---|---|---|---|---
$\Delta^{\prime}_{x}$ [$\mu$m] | 33.275 | $33.172\pm{}0.025$ | $33.303\pm{}0.023$ | $33.266\pm{}0.023$ | $33.264\pm{}0.024$ | $33.269\pm{}0.028$
| (0.1) | $33.171\pm{}0.025$ | $33.301\pm{}0.026$ | $33.269\pm{}0.026$ | $33.267\pm{}0.028$ | $33.271\pm{}0.029$
$\Delta^{\prime}_{y}$ [$\mu$m] | -16.900 | $-17.552\pm{}0.010$ | $-17.268\pm{}0.012$ | $-17.368\pm{}0.011$ | $-17.368\pm{}0.013$ | $-17.367\pm{}0.013$
| (2) | $-17.550\pm{}0.011$ | $-17.265\pm{}0.011$ | $-17.368\pm{}0.012$ | $-17.367\pm{}0.013$ | $-17.367\pm{}0.012$
$\Delta^{\prime}_{z}$ [$\mu$m] | -26.559 | $-26.561\pm{}0.026$ | $-26.559\pm{}0.027$ | $-26.559\pm{}0.027$ | $-26.557\pm{}0.023$ | $-26.554\pm{}0.023$
| (0.1) | $-26.557\pm{}0.021$ | $-26.557\pm{}0.026$ | $-26.556\pm{}0.025$ | $-26.556\pm{}0.024$ | $-26.557\pm{}0.027$
$a^{\prime}_{d11}$ | 9.877 | $9.8866\pm{}0.0004$ | $9.8865\pm{}0.0004$ | $9.8865\pm{}0.0003$ | $9.8865\pm{}0.0004$ | $9.8866\pm{}0.0004$
$[\times{}10^{-3}]$ | (0.15) | $9.8866\pm{}0.0004$ | $9.8865\pm{}0.0004$ | $9.8865\pm{}0.0003$ | $9.8865\pm{}0.0003$ | $9.8865\pm{}0.0004$
$a_{d12}$ | 1.5983 | $1.5941\pm{}0.0004$ | $1.5940\pm{}0.0004$ | $1.5941\pm{}0.0004$ | $1.5940\pm{}0.0003$ | $1.5941\pm{}0.0003$
$[\times{}10^{-3}]$ | (0.05) | $1.5940\pm{}0.0004$ | $1.5940\pm{}0.0004$ | $1.5941\pm{}0.0004$ | $1.5940\pm{}0.0004$ | $1.5940\pm{}0.0003$
$a_{d13}$ | -1.4026 | $-1.4067\pm{}0.0004$ | $-1.4068\pm{}0.0003$ | $-1.4068\pm{}0.0003$ | $-1.4067\pm{}0.0003$ | $-1.4068\pm{}0.0004$
$[\times{}10^{-3}]$ | (0.05) | $-1.4068\pm{}0.0003$ | $-1.4068\pm{}0.0004$ | $-1.4068\pm{}0.0004$ | $-1.4068\pm{}0.0004$ | $-1.4068\pm{}0.0004$
$a_{c12}$ | 2.99 | $4.39\pm{}0.05$ | $2.84\pm{}0.05$ | $2.83\pm{}0.05$ | $2.83\pm{}0.04$ | $2.83\pm{}0.05$
$[\times{}10^{-3}]$ | (0.9) | $4.39\pm{}0.04$ | $2.83\pm{}0.04$ | $2.82\pm{}0.05$ | $2.83\pm{}0.04$ | $2.84\pm{}0.04$
$a_{c13}$ | -3.75 | $-4.95\pm{}0.04$ | $-3.55\pm{}0.05$ | $-3.55\pm{}0.05$ | $-3.56\pm{}0.05$ | $-3.56\pm{}0.04$
$[\times{}10^{-3}]$ | (0.9) | $-4.94\pm{}0.05$ | $-3.55\pm{}0.05$ | $-3.54\pm{}0.05$ | $-3.56\pm{}0.04$ | $-3.56\pm{}0.05$
$K_{2d}$ [s2/m] | 5612.6 | $5746.6\pm{}1.2$ | $5746.7\pm{}1.2$ | $5746.7\pm{}1.2$ | $5746.6\pm{}1.1$ | $5746.7\pm{}1.2$
| (250) | $5746.7\pm{}1.2$ | $5746.7\pm{}1.1$ | $5746.7\pm{}1.2$ | $5746.9\pm{}1.1$ | $5746.8\pm{}1.2$
$\delta$ [$\times{}10^{-15}$] | 10 | $7.95\pm{}7.39$ | $10.31\pm{}8.07$ | $11.63\pm{}7.97$ | $11.29\pm{}6.58$ | $9.52\pm{}6.68$
| (7) | $8.45\pm{}7.48$ | $9.58\pm{}7.00$ | $10.81\pm{}7.80$ | $9.50\pm{}7.34$ | $8.94\pm{}7.30$

After focusing on the $\Delta^{\prime}_{x}$ estimation above, we now consider
all parameters to mimic a real in-flight iterative calibration. We estimate
the parameters listed in the first column of Table 1, with input values listed
in the second column; the required combined accuracy and precision of the
estimation of each parameter is given in Table 2. The remaining columns list
the values measured for each parameter at each iteration of the process. For
each parameter, the first line gives the value estimated when considering
error propagation from one iteration to the other; in the second line, we
ignore error propagation and use for the $i$th iteration the best fit of each
parameter obtained at the $(i-1)$th iteration. The best fit and uncertainty
for each parameter are computed as the mean and rms of the distribution of 100
estimates from 100 simulations for each parameter and iteration. For the
estimation of a given parameter, all simulations share the same deterministic
signals but have a different realisation of the instrumental noise. Note
however that because of computational limitations, we could not use the same
set of simulations in both cases. This led to different best fits when
propagating uncertainties or not (nevertheless, the differences are well
within the error bars).

It is clear that after the first iteration where the estimation of some
parameters is somewhat biased (most notably $\Delta^{\prime}_{x}$, $a_{c12}$
and $a_{c13}$), the iterative correction of instrumental parameters allows us
to recover the input values, well within the required accuracy. In practice,
based on these results, we iterate on instrumental parameters while estimates
from one iteration to the next do not vary significantly, with a maximum
number of iterations set to five [6]. Note also that the $a_{c12}$ and
$a_{c13}$ parameters, though correctly estimated here, are not measured in
real data [6, 8] since they are close to their required value from
construction, and the noise (higher than in these simulations) does not allow
for a precise measurement.

The discussions in Sect. 3.5 and B show that error propagation is not trivial
in our iterative calibration. The comparison of each instrumental parameter’s
lines in Table 1 (the first one corresponding to the case where we propagate
errors, the second one to the case where we ignore error propagation) shows
that we can ignore error propagation when estimating instrumental parameters.

The last line of Table 1 gives the measurement of the Eötvös parameter after
correction of the estimated instrumental parameters. It is clear that we are
able to accurately recover the input EPV signal with the expected precision.
Just like in the case of instrumental parameters, the uncertainty on the
Eötvös parameter is not impacted by the uncertainty propagation from one
iteration to the next. We can therefore reliably ignore the error propagation
in our iterative estimations.

##  6 Conclusion

In this paper, we summarised MICROSCOPE’s data analysis process. Building on
the measurement equation, we showed how we can either estimate, model or
ignore instrumental parameters. The estimation of those that cannot be ignored
nor modelled is performed through an iterative weighted least square fit in
the frequency domain. We provided an extensive characterisation of inpainting,
a gap-filling technique that we adapted to MICROSCOPE. This exercice, based on
worst-case scenario numerical simulations, shows the behavior of the algorithm
in a complex data-processing pipeline, and can be useful to other experimental
data analyses.

Using well-controlled numerical simulations, we then showed that our iterative
least-square method is robust to estimate instrumental parameters and correct
for them in order to reliably measure the Eötvös ratio. In particular, we
showed that our main least-square estimator is intrinsically unbiased.
Combined with our previous works [15, 16, 17, 18], those results prove that
our data analysis pipeline allows us to reliably measure an EPV even in the
presence of missing data and instrumental imperfections.

Finally, we discussed the non-trivial, but important, problem of uncertainties
propagation in our pipeline. Like in any precision experiment, the problem of
uncertainty estimation is central to the MICROSCOPE data analysis. We showed
that not only does our iterative least-square technique provide correct error
bars, but we can also safely ignore the problem of uncertainty propagation
throughout iterations. This paper therefore justifies the results provided in
Refs. [1, 2, 3, 8].

This work makes use of technical data from the CNES-ESA-ONERA-CNRS-OCA
MICROSCOPE mission, and has received financial support from ONERA and CNES. We
thank Bruno Christophe, Bernard Foulon and Isabelle Petitbon, as well as the
members of the MICROSCOPE Science Working Group for useful discussions.
Special thanks go to Pierre Fayet and Pieter Visser for detailed comments on
the manuscript. JB and SP acknowledge the financial support of the UnivEarthS
Labex program at Sorbonne Paris Cité (ANR-10-LABX-0023 and
ANR-11-IDEX-0005-02).

##  Appendix A In-flight instrumental parameters estimation

Table 2 lists the parameters that we can calibrate, and summarizes the
techniques used to perform their estimation. The third column of the table
lists the precision and accuracy required on the estimation of each parameter
to satisfy the overall MICROSCOPE goal to reach the $10^{-15}$ level for the
Eötvös parameter; the fourth column shows the maximum value allowed by design
of the instrument; those numbers were obtained with a performance analysis,
whose description is beyond the scope of this paper. See Ref. [6] for a
complete description of which parameters can be estimated in flight.

Table 2: Parameters that can be calibrated in orbit, method to perform their
calibration, and precision required to satisfy MICROSCOPE’s objective.
Oscillations are given along or about axes in the instrument frame. $f_{\rm
EP}$ is the frequency of the test, where we may expect to detect a WEP
violation. The fourth column gives the maximal allowed value by design. See
Ref. [10, 6] for details.

Parameter | Calibration technique | Aimed precision | Max. allowed value
---|---|---|---
$\Delta^{\prime}_{x}$ | No satellite maneuver; take advantage of Earth’s GGT signal at $2f_{\rm EP}$ | $0.1\mu m$ | 20 $\mu$m
$\Delta^{\prime}_{y}$ | Angular oscillation of the satellite about the $z$-axis | $2\mu m$ | 20 $\mu$m
$\Delta^{\prime}_{z}$ | No satellite maneuver; take advantage of Earth’s GGT signal at $2f_{\rm EP}$ | $0.1\mu m$ | 20 $\mu$m
$a_{c12}$ | Angular oscillation of the satellite about the $x$-axis and linear oscillation |  $9\times{}10^{-4}$rad |  $2.6\times{}10^{-3}$rad
| of the test mass along the $z$-axis |  |
$a_{c13}$ | Angular oscillation of the satellite about the $x$-axis and linear oscillation |  $9\times{}10^{-4}$rad |  $2.6\times{}10^{-3}$rad
| of the test mass along the $y$-axis |  |
$a^{\prime}_{d11}$ | Linear oscillation of the satellite along the $x$-axis | $1.5\times{}10^{-4}$ | 0.01
$a_{d12}$ | Linear oscillation of the satellite along the $y$-axis |  $5\times{}10^{-5}$rad |  $1.5\times{}10^{-3}$rad
$a_{d13}$ | Linear oscillation of the satellite along the $z$-axis |  $5\times{}10^{-5}$rad |  $1.5\times{}10^{-3}$rad
$K_{2xx}^{(d)}/\left(K_{1x}^{(c)}\right)^{2}$ | Linear oscillation of the satellite along the $x$-axis | 250 s2/m | 14000 s2/m

##  Appendix B Error propagation through least-square estimation

This appendix presents a pedagogical derivation of error propagation when
correcting the measured differential acceleration from estimated instrumental
parameters. We show two different approaches: the first one is a traditional
variational approach, while the second one computes the full variance of the
calibrated differential acceleration. We show that the latter provides a more
complete view of the effects of propagating uncertainties when correcting for
(estimated) instrumental parameter.

###  B.1 Uncorrected differential acceleration and estimated parameters

The uncorrected differential acceleration (Eq. 1) can be re-written as:

$\Gamma^{d}(t)=2\Gamma_{\delta}(t)+\sum_{j,k=1}^{3}(T_{jk}(t)-In_{jk}(t))a_{c1j}\Delta_{k}(t)\\\
+2\sum_{k=1}^{3}(-1)^{k}\left[a_{c1i}\dot{\Delta}_{j}(t)-a_{c1j}\dot{\Delta}_{i}(t)\right]_{i<j,i,j\neq
k}\Omega_{k}\\\
-\sum_{k=1}^{3}a_{c1k}\ddot{\Delta}_{k}(t)+2\sum_{k=1}^{3}a_{d1k}\tilde{\Gamma_{k}^{c}}(t)\\\
+K_{21}\left(\tilde{\Gamma}^{(1)}(t)\right)^{2}-K_{22}\left(\tilde{\Gamma}^{(2)}(t)\right)^{2}+2n_{d}(t),$
(15)

where $\Gamma_{\delta}(t)$ gathers terms depending on the Eötvös ratio
$\delta$, and where the indices $i,j,k$ are either in (1,2) or ($x$,$y$,$z$)
depending on the variable they index. In Eq. (15), we ignore angular-to-linear
couplings. In this equation, the common-mode and individual accelerations,
marked by a tilde, are noise-corrected (e.g.
$\tilde{\Gamma_{k}^{c}}(t)=\Gamma_{k}^{c}(t)-n_{c}(t)$).

The common-mode and differential-mode noises are assumed normally distributed
with variance $\sigma_{d}^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}$ and
$\sigma_{c}^{2}=(\sigma_{1}^{2}+\sigma_{2}^{2})/4$, where $\sigma_{k}^{2}$ is
the variance for the $k$th ($k=1,2$) sensor’s noise.

The parameters estimated before correction are the following: $a_{c1k}$,
$a_{d1k}$, $\Delta_{k}$, $K_{2k}$. In the following, we assume that their
estimators $\hat{a}_{c1k}$, $\hat{a}_{d1k}$, $\hat{\Delta}_{k}$,
$\hat{K}_{2k}$ are unbiased and of variance $\sigma_{ac1k}$, etc.

Note that although $a_{c11}$ cannot be estimated, it is possible to consider
it as an unknown parameter, and assume it is a random variable of mean 1 and
variance given by its specifications (or an educated guess). For the sake of
clarity, we shall not pursue this possibility here.

###  B.2 Synthetic uncorrected and calibrated accelerations

In order to simplify the computation of the calibrated differential
acceleration, of its expectation value and of its variance after propagating
the errors on the estimation of instrumental parameters, we define a synthetic
(uncorrected) differential acceleration as:

$\Gamma^{d}(t)=\Gamma_{\delta}(t)+\kappa_{a}aY(t)+\kappa_{bc}bcZ(t)+\kappa_{k}k\tilde{\Gamma^{c}}(t)+l_{1}\left(\tilde{\Gamma}^{(1)}\right)^{2}-l_{2}\left(\tilde{\Gamma}^{(2)}\right)^{2}+n_{d}(t),$
(16)

where $a$, $b$, $c$, $k$, $l_{1}$ and $l_{2}$ are estimated, $\kappa_{i}$ are
numerical (constant) factors, $Y(t)$ and $Z(t)$ are deterministic signals that
can be measured or modeled, and we kept the common-mode and individual
accelerations unchanged. Although $b$ and $c$ may be degenerate, we will
assume that they are not correlated. This equation encompasses all types of
terms found in Eq. (15).

The calibrated synthetic acceleration can then be shown to be:

$\Gamma_{\rm
cal}^{d}(t)=\Gamma_{\delta}(t)+\kappa_{a}(a-\hat{a})Y(t)+\kappa_{bc}(bc-\hat{b}\hat{c})Z(t)+\kappa_{k}(k-\hat{k})\tilde{\Gamma}^{c}(t)\\\
+(l_{1}-\hat{l_{1}})\left[\tilde{\Gamma}^{(1)}(t)\right]^{2}-2\hat{l_{1}}\tilde{\Gamma}^{(1)}(t)n_{1}(t)-(l_{2}-\hat{l_{2}})\left[\tilde{\Gamma}^{(2)}(t)\right]^{2}+2\hat{l_{2}}\tilde{\Gamma}^{(2)}(t)n_{2}(t)\\\
+\frac{1}{2}\left[(n_{1}(t)-n_{2}(t))-\kappa_{k}\hat{k}(n_{1}(t)+n_{2}(t))\right]-\hat{l_{1}}n_{1}^{2}(t)+\hat{l_{2}}n_{2}^{2}(t)$
(17)

In what follows, we assume that we wish to estimate the Eötvös ratio $\delta$
(hence, the assumption that the estimates of the instrumental parameters are
unbiased). The discussion can readily be generalized to the estimation of any
instrumental parameter given priors on the others.

###  B.3 Error propagation. Method 1: variational approach

We first use a variational approach to compute the variance of the calibrated
differential acceleration (see e.g. [33]).

###  B.4 Expectation value of calibrated differential acceleration

Under the assumption that all estimates are unbiased, the expectation value of
Eq. (17) is

$\mathrm{E}[\Gamma_{\rm
cal}^{d}(t)]=\Gamma_{\delta}(t)-\hat{l_{1}}\sigma_{1}^{2}+\hat{l_{2}}\sigma_{2}^{2}.$
(18)

It should be noted that the calibration entails a constant non-zero bias.
However, the same bias is applied at all times $t$, and therefore does not
affect the estimation of the amplitude of the WEP violation signal.

####  B.4.1 Error propagation

Taking the sum of partial derivatives of Eq. (17) and assuming that
$|\Gamma^{(1)}|\ll{}1$ and $|\Gamma^{(2)}|\ll{}1$, we get:

$\sigma_{\Gamma_{\rm
cal}}^{2}=\kappa_{a}^{2}Y^{2}(t)\sigma_{a}^{2}+\kappa_{bc}Z^{2}(t)\hat{b}^{2}\hat{c}^{2}\left(\frac{\sigma_{b}^{2}}{\hat{b}}+\frac{\sigma_{c}^{2}}{\hat{c}^{2}}\right)+\kappa_{k}^{2}\left[\Gamma^{c}(t)\right]^{2}\sigma_{k}^{2}\\\
+\left[\frac{1-\kappa_{k}\hat{k}}{2}-2\hat{l_{1}}\Gamma^{(1)}(t)\right]^{2}\sigma_{1}^{2}+\left[2\hat{l_{2}}\Gamma^{(2)(t)}-\frac{1+\kappa_{k}\hat{k}}{2}\right]^{2}\sigma_{2}^{2}$
(19)

Note that in Eq. (19), the common-mode and individual accelerations that
appear are the effectively measured (hence, noisy) ones, contrary to those
which appeared in the uncorrected acceleration.

####  B.4.2 Actual calibrated acceleration: Expectation value and variance

Returning to the actual MICROSCOPE accelerations, Eqs. (18) and (19) give:

$\mathrm{E}[\Gamma_{\rm
cal}^{d}(t)]=\Gamma_{\delta}(t)-\frac{1}{2}(\hat{K}_{21}\sigma_{1}^{2}+\hat{K}_{22}\sigma_{2}^{2})$
(20)

and

$\sigma_{\Gamma_{\rm
cal}(t)}^{2}=\frac{1}{4}\sum_{j,k=1}^{3}\left[T_{jk}(t)-In_{jk}(t)\right]^{2}\hat{a}_{c1j}^{2}\hat{\Delta}_{k}^{2}\left(\frac{\sigma_{\Delta_{k}}^{2}}{\hat{\Delta}_{k}^{2}}+\frac{\sigma_{ac1j}^{2}}{\hat{a_{c1j}}^{2}}\right)\\\
+\sum_{k=1}^{3}\left[\dot{\Delta}_{j}(t)^{2}\sigma_{ac1i}^{2}+\dot{\Delta}_{i}(t)^{2}\sigma_{ac1j}^{2}\right]_{i,j\neq
k,i<j}\Omega_{k}(t)^{2}\\\
+\frac{1}{4}\sum_{k=1}^{3}\ddot{\Delta}_{k}(t)^{2}\sigma_{ac1k}^{2}+\sum_{k=1}^{3}\left[\Gamma^{c}(t)\right]^{2}\sigma_{ad1k}^{2}\\\
+\left[\frac{1-\sum_{k=1}^{3}\hat{a}_{d1k}}{2}-2\hat{K}_{21}\Gamma^{(1)}(t)\right]^{2}\sigma_{1}^{2}\\\
+\left[2\hat{K}_{22}\Gamma^{(2)}(t)-\frac{1+\sum_{k=1}^{3}\hat{a}_{d1k}}{2}\right]^{2}\sigma_{2}^{2}.$
(21)

Note that the uncertainties on the quadratic factors $\sigma_{K2i}^{2}$ do not
appear in this equation. This is because we assumed that
$|\Gamma_{(i)}|\ll{}1$.

Eq. (21) shows that the instrumental parameters themselves and the
uncertainties on their estimation bring up two different effects on the
variance of the calibrated differential acceleration:

  * •

extra-contributors to total variance: uncertainties on the estimated
instrumental parameters add up to the noise, therefore increasing the
calibrated differential acceleration’s variance, and increasing the
uncertainty on the Eötvös parameter estimated with a least-square fit of the
calibrated differential acceleration. Those extra contributions can be taken
into account simply by adding them in the data covariance passed to the least-
square fit. Another method (more computationally expensive) is to use Monte
Carlo simulations where we vary the correction terms within their allowed
bounds.

  * •

modification of the noise: the calibrated differential acceleration noise
variance is not simply $\sigma_{d}^{2}=\sigma_{1}^{2}+\sigma_{2}^{2}$ as that
of the uncorrected differential acceleration noise. Instead, it is modified by
the presence of non-zero differential parameters ($a_{d1k}$) that couple to
the common-mode acceleration and of non-zero quadratic factors that couple to
the individual accelerations. The knowledge of the best estimates for those
instrumental parameters allows us to quantify the modification of the noise
(which, given the values of the involved parameters, remains negligible).
However, Eq. (21) is valid only for a given set of instrumental parameters’
estimators, and does not tell us anything about the distribution of the
corrected differential acceleration noise, and hence on the uncertainty on the
calibrated differential acceleration noise. A more general approach is
necessary (see below).

###  B.5 Error propagation. Method 2: general variance analysis of the
calibrated differential acceleration

For this analysis, we first go back to our synthetic model (17), whose
variance we directly compute.

####  B.5.1 Synthetic model

Under the same assumptions as for Eq. (19), we find (with a straightforward
but tedious algebra):

$\mathrm{Var}[\Gamma_{\rm
cal}(t)]=\kappa_{a}^{2}Y^{2}(t)\sigma_{a}^{2}+\kappa_{bc}^{2}Z^{2}(t)\hat{b}^{2}\hat{c}^{2}\left[\frac{\sigma_{b}^{2}\sigma_{c}^{2}}{\hat{b}^{2}\hat{c}^{2}}+\frac{\sigma_{b}^{2}}{\hat{b}^{2}}+\frac{\sigma_{c}^{2}}{\hat{c}^{2}}\right]+\kappa_{k}^{2}\left[\tilde{\Gamma^{c}}(t)\right]^{2}\sigma_{k}^{2}+l_{1}^{2}\sigma_{1}^{4}+l_{2}\sigma_{2}^{4}\\\
+\frac{1}{4}\left[1+\kappa_{k}(\kappa_{k}\sigma_{k}^{2}-2k+k^{2})+16(\sigma_{l1}^{2}+l_{1}^{2})\left[\tilde{\Gamma}^{(1)}(t)\right]^{2}-2(1-\kappa_{k}k)l_{1}\tilde{\Gamma}^{(1)}(t)\right]\sigma_{1}^{2}\\\
+\frac{1}{4}\left[1+\kappa_{k}(\kappa_{k}\sigma_{k}^{2}+2k+k^{2})+16(\sigma_{l_{2}}^{2}+l_{2}^{2})\left[\tilde{\Gamma}^{(2)}(t)\right]^{2}-2(1-\kappa_{k}k)l_{2}\tilde{\Gamma}^{(2)}(t)\right]\sigma_{2}^{2}.$
(22)

Note that although we assume $\sigma_{1,2}\ll{}1$, we do not ignore the
$\sigma_{1,2}^{4}$ terms since the $l_{1,2}$ terms may be significant (they
encode the quadratic factors).

We can observe that the contributors from the error propagation (those which
do not contribute to the noise) are the same as in Eq. (19), except for the
$\sigma_{b}^{2}\sigma_{c}^{2}/\hat{b}^{2}\hat{c}^{2}$ term, that pops-up here
since we did not assume that $b$ and $c$ were independent. However, for
MICROSCOPE, we always have $\sigma_{b,c}^{2}\ll{}(b,c)^{2}$, so we will ignore
this term in the following, so that the contribution to the uncertainty on the
Eötvös parameter from error propagation is the same in our two analyses (Eqs.
(19) and (22)).

However, the noise term is much more complex in Eq. (22) than in Eq. (19).
This is because it encompasses the uncertainty on the estimated instrumental
parameters. In that sense, whereas Eq. (19) tells us how the noise is modified
when we correct the differential acceleration for a given set of instrumental
parameters, Eq. (22) tells us how the modified noise is distributed. It
provides a more conservative uncertainty propagation, and should be favored
when analyzing data.

####  B.5.2 Actual calibrated differential acceleration

Ignoring the $\sigma_{b}^{2}\sigma_{c}^{2}/\hat{b}^{2}\hat{c}^{2}$ term from
Eq. (22), the total variance for the actual MICROSCOPE acceleration is

$\mathrm{Var}[\Gamma_{\rm
cal}(t)]=\frac{1}{4}\sum_{j,k=1}^{3}\left[T_{jk}(t)-In_{jk}(t)\right]^{2}\hat{a}_{c1j}^{2}\hat{\Delta}_{k}^{2}\left(\frac{\sigma_{\Delta_{k}}^{2}}{\hat{\Delta}_{k}^{2}}+\frac{\sigma_{ac1j}^{2}}{\hat{a}_{c1j}^{2}}\right)\\\
+\sum_{k=1}^{3}\left[\dot{\Delta}_{j}(t)^{2}\sigma_{ac1i}^{2}+\dot{\Delta}_{i}(t)^{2}\sigma_{ac1j}^{2}\right]_{i,j\neq
k,i<j}\Omega_{k}(t)^{2}\\\
+\frac{1}{4}\sum_{k=1}^{3}\ddot{\Delta}_{k}(t)^{2}\sigma_{ac1k}^{2}+\sum_{k=1}^{3}\left[\Gamma^{c}(t)\right]^{2}\sigma_{ad1k}^{2}\\\
+\frac{1}{4}\left\\{1+\sum_{k=1}^{3}\sigma_{ad1k}^{2}-2\sum_{k=1}^{3}a_{d1k}+\sum_{k=1}^{3}a_{d1k}^{2}+16\left(\sigma_{K21}^{2}+K_{21}^{2}\right)\left[\tilde{\Gamma}^{(1)}\right]^{2}\right.\\\
\left.-2\left(1-\sum_{k=1}^{3}a_{d1k}\right)K_{21}\tilde{\Gamma}^{(1)}\right\\}\sigma_{1}^{2}\\\
+\frac{1}{4}\left\\{1+\sum_{k=1}^{3}\sigma_{ad1k}^{2}+2\sum_{k=1}^{3}a_{d1k}+\sum_{k=1}^{3}a_{d1k}^{2}+16\left(\sigma_{K22}^{2}+K_{22}^{2}\right)\left[\tilde{\Gamma}^{(2)}\right]^{2}\right.\\\
\left.-2\left(1+\sum_{k=1}^{3}a_{d1k}\right)K_{22}\tilde{\Gamma}^{(2)}\right\\}\sigma_{2}^{2}\\\
+K_{21}^{2}\sigma_{1}^{4}+K_{22}^{2}\sigma_{2}^{4}$ (23)

The noise term in Eq. (23) makes it clear that the $a_{d1k}$ parameters and
quadratic factors, as well as the uncertainty on their estimation, modify the
measured noise.

Although in practice the uncertainty on those parameters can be seen as
bringing an uncertainty on the noise (the noise is corrected by some uncertain
estimates of the parameters), Eq. (23) combines those uncertainties to provide
an upper bound of the measured noise. That is, it provides the most
conservative variance of the calibrated differential acceleration. Therefore,
it should be used to completely take into account the errors in the estimate
of instrumental parameters, and should be included in the least-square fit.

Note that in this equation, the values of $a_{d1k}$ and $K_{2i}$ are not their
estimates, but their “real” value. In practice, assuming that our estimates
are unbiased, or at least give a correct order of magnitude of the real value,
we can replace them by their estimated values.

## References

## References

  * [1] Touboul P, Métris G, Rodrigues M, André Y, Baghi Q, Bergé J, Boulanger D, Bremer S, Carle P, Chhun R, Christophe B, Cipolla V, Damour T, Danto P, Dittus H, Fayet P, Foulon B, Gageant C, Guidotti P Y, Hagedorn D, Hardy E, Huynh P A, Inchauspe H, Kayser P, Lala S, Lämmerzahl C, Lebat V, Leseur P, Liorzou F, List M, Löffler F, Panet I, Pouilloux B, Prieur P, Rebray A, Reynaud S, Rievers B, Robert A, Selig H, Serron L, Sumner T, Tanguy N and Visser P 2017 Physical Review Letters 119 231101 (Preprint 1712.01176)
  * [2] Touboul P, Métris G, Rodrigues M, André Y, Baghi Q, Bergé J, Boulanger D, Bremer S, Chhun R, Christophe B, Cipolla V, Damour T, Danto P, Dittus H, Fayet P, Foulon B, Guidotti P Y, Hardy E, Huynh P A, Lämmerzahl C, Lebat V, Liorzou F, List M, Panet I, Pires S, Pouilloux B, Prieur P, Reynaud S, Rievers B, Robert A, Selig H, Serron L, Sumner T and Visser P 2019 Classical and Quantum Gravity 36 225006
  * [3] Touboul P and Microscope team in prep Class. Quant. Grav.
  * [4] Touboul P, Rodrigues M, Métris G, Chhun R, Robert A, Baghi Q, Hardy E, Bergé J, Boulanger D, Christophe B, Cipolla V, Foulon B, Guidotti P Y, Huynh P A, Lebat V, Liorzou F, Pouilloux B, Prieur P and Reynaud S 2020 arXiv e-prints arXiv:2012.06472 (Preprint 2012.06472)
  * [5] Touboul P, Métris G, Rodrigues M, Bedouet J, Bergé J, Boulanger D, Chhun R, Christophe B, Cipolla V, Damour T, Dittus H, Fayet P, Foulon B, Guidotti P Y, Hardy E, Huynh P A, Lämmerzahl C, Lebat V, Liorzou F, List M, Pouilloux B, Prieur P, Reynaud S, Rievers B, Robert A, Serron L, Sumner T and Visser P in prep Class. Quant. Grav.
  * [6] Hardy E and Microscope team in prep Class. Quant. Grav.
  * [7] Bergé J, Baghi Q, Robert A, Rodrigues M, Foulon B, Hardy E, Métris G, Pires S and Touboul P 2020 arXiv e-prints arXiv:2012.06485 (Preprint 2012.06485)
  * [8] Métris G and Microscope team in prep Class. Quant. Grav.
  * [9] Liorzou F, Touboul P, Rodrigues M, Métris G, André Y, Bergé J, Boulanger D, Bremer S, Chhun R, Christophe B, Danto P, Foulon B, Hagedorn D, Hardy E, Huynh P A, Lämmerzahl C, Lebat V, List M, Löffler F, Rievers B, Robert A and Selig H 2020 arXiv e-prints arXiv:2012.11232 (Preprint 2012.11232)
  * [10] Hardy E, Levy A, Métris G, Rodrigues M and Touboul P 2013 Sapce Science Reviews 180 177–191 (Preprint 1707.08024)
  * [11] Robert A, Cipolla V, Prieur P, Touboul P, Métris G, Rodrigues M, André Y, Bergé J, Boulanger D, Chhun R, Christophe B, Guidotti P Y, Hardy E, Lebat V, Lienart T, Liorzou F and Pouilloux B 2020 arXiv e-prints arXiv:2012.06479 (Preprint 2012.06479)
  * [12] Petit G and Luzum B (eds) 2010 IERS Conventions (2010) International Earth Rotation and Reference Systems Service (IERS) (Verlag des Bundesamts für Kartographie und Geodäsie)
  * [13] Métris G, Xu J and Wytrzyszczak I 1998 Celestial Mechanics and Dynamical Astronomy 71 137–151
  * [14] Will C M 2014 Living Reviews in Relativity 17 4 (Preprint 1403.7377)
  * [15] Baghi Q, Métris G, Bergé J, Christophe B, Touboul P and Rodrigues M 2015 Phys. Rev. D 91 062003 (Preprint 1503.01470)
  * [16] Baghi Q, Métris G, Bergé J, Christophe B, Touboul P and Rodrigues M 2016 Phys. Rev. D 93 122007 (Preprint 1608.08530)
  * [17] Bergé J, Pires S, Baghi Q, Touboul P and Métris G 2015 Phys. Rev. D 92 112006 (Preprint 1512.00492)
  * [18] Pires S, Bergé J, Baghi Q, Touboul P and Métris G 2016 Phys. Rev. D 94 123015 (Preprint 1612.05452)
  * [19] Lupton R 1993 Statistics in theory and practice (Princeton University Press) ISBN 9780691074290
  * [20] Hardy E, Levy A, Rodrigues M, Touboul P and Métris G 2013 Advances in Space Research 52 1634–1646 (Preprint 1707.07630)
  * [21] Abbott B P et al. (LIGO Scientific, Virgo) 2020 Class. Quant. Grav. 37 055002 (Preprint 1908.11170)
  * [22] Cornish N J and Littenberg T B 2015 Classical and Quantum Gravity 32 135012 (Preprint 1410.3835)
  * [23] Zackay B, Venumadhav T, Roulet J, Dai L and Zaldarriaga M 2019 (Preprint 1908.05644)
  * [24] Venumadhav T, Zackay B, Roulet J, Dai L and Zaldarriaga M 2019 Phys. Rev. D 100 023011 (Preprint 1902.10341)
  * [25] Wei W and Huerta E 2020 Phys. Lett. B 800 135081 (Preprint 1901.00869)
  * [26] Torres-Forné A, Cuoco E, Font J A and Marquina A 2020 Phys. Rev. D 102 023011 (Preprint 2002.11668)
  * [27] Bergé J, Massey R, Baghi Q and Touboul P 2019 Monthly Notices of the Royal Astronomical Society 486 544–559 (Preprint 1903.05837)
  * [28] Pires S, Starck J L, Amara A, Teyssier R, Réfrégier A and Fadili J 2009 Monthly Notices of the Royal Astronomical Society 395 1265–1279 (Preprint 0804.4068)
  * [29] Pires S, Mathur S, García R A, Ballot J, Stello D and Sato K 2015 Astronomy and Astrophysics 574 A18 (Preprint 1410.6088)
  * [30] Prieur P, Lienart T, Rodrigues M, Touboul P, Denver T, Jorgensen J L, Bang A M and Metris G 2017 Microscope mission: on-orbit assessment of the drag-free and attitude control system ISTS-2017-d-038/ISSFD-2017-038 1-10, 26th Int Symp Space Flight Dynamics
  * [31] Touboul P 2009 Space Science Reviews 148 455–474
  * [32] Bourda G, Charlot P and Biancale R 2008 ArXiv e-prints (Preprint 0803.1432)
  * [33] Ku H H 1966 Journal of research of the National Bureau of Standards 70 4