MICROSCOPE mission: statistics and impact of glitches on the test of the weak equivalence principle^*

During its science operations, the MICROSCOPE satellite orbited the Earth on a 7090 km semi-major axis, nearly circular, sun-synchronous (nearly polar) orbit.

The instrument consists of two differential accelerometers (or sensor units). Each differential accelerometer includes two coaxial sensors each consisting in a cylindrical test masses and associated electrodes and electronics [2]. We call SUEP the sensor unit for the equivalence principle test, with test masses made of different composition, and SUREF the reference sensor unit with test masses of the same composition (which should not detect a WEP violation). For each sensor, accelerometry data is sampled at a 4 Hz rate. In this paper, we use the measured acceleration from each sensor, as well as the differential acceleration defined as the difference of acceleration between the two test masses of a given sensor unit (see reference [20] for more details about data levels). Raw data are used in section 3 and calibrated data in section 4. Note that in the MICROSCOPE realm, the acceleration of a test mass is defined as the electrostatic force per unit mass required to keep the test mass motionless [21]. All science sessions (as defined in reference [20]) are used, amounting to 2340 orbits (about 160 days worth of data). Two satellite's spin rates were used: f_spin = 35/2f_orb for some sessions and f_spin = 9/2f_orb for others.

The measurement relies on the intertwined control loops of both sensors and of the drag-free and attitude control system (DFACS). A block diagram of the system is shown in figure 1, where the sensors' transfer functions are labelled A₁ and A₂, while G is that of the DFACS. The drag-free system is controled by only one sensor. The accelerations M₁ and M₂ measured by individual sensors are functions of the external accelerations E (such as atmospheric drag, solar radiation pressure and glitches—we show in section 3.3 that glitches are events external from individual sensors), local accelerations (δ_i represent a composition-dependence of the ratio between the inertial and gravitational masses—${\delta }_{\mathrm{i}}=[{({m}_{\mathrm{g}}/{m}_{\mathrm{i}})}_{\mathrm{i}}-1]g$, where g is the Earth gravity acceleration—and p_i other, generic perturbations), instrumental noise n_i and of the closed-loop transfer functions.

Zoom In Zoom Out Reset image size

Assuming that the transfer functions of the two sensors have been matched up to an accuracy dA, such that A₁(f) = A(f) and A₂(f) = A(f) + dA(f), the acceleration measured by the sensor which controls the drag-free is then

$$\begin{equation}{M}_{1}=\frac{A}{1+GA}(E+{\delta }_{1}+{p}_{1})+\frac{{n}_{1}}{1+GA},\end{equation} \tag{ 1 }$$

where all variables depend on the frequency, e.g. G = G(f), and that measured at the other sensor is

$$\begin{equation}{M}_{2}=\frac{A+\mathrm{d}A}{1+GA}E+(A+\mathrm{d}A)({\delta }_{2}+{p}_{2})-\frac{GA(A+\mathrm{d}A)}{1+GA}({\delta }_{1}+{p}_{1})-\frac{G(A+\mathrm{d}A)}{1+GA}{n}_{1}+{n}_{2}.\end{equation} \tag{ 2 }$$

The differential acceleration, from which we measure the WEP violation, is

$$\begin{align}\hfill {M}_{d}& \equiv {M}_{1}-{M}_{2}=-\frac{\mathrm{d}A}{1+GA}E+A({\delta }_{1}-{\delta }_{2})+A({p}_{1}-{p}_{2})\hfill \\ \hfill & \quad +\mathrm{d}A\left(\frac{GA}{1+GA}{\delta }_{1}-{\delta }_{2}\right)+\mathrm{d}A\left(\frac{GA}{1+GA}{p}_{1}-{p}_{2}\right)+{n}_{1}-{n}_{2}+\frac{G\,\mathrm{d}A}{1+GA}{n}_{1}.\hfill \end{align} \tag{ 3 }$$

If the transfer functions of both sensors are equal and perfectly matched (dA = 0), external forces do not impact the differential acceleration, thereby glitches should not have an effect on the estimation of the WEP violation.

The WEP is estimated at frequency f_EP, where local accelerations p_i ≪ n_i (by design [21]), and where, also by design A(f_EP) ≈ 1, G(f_EP) ≈ 10⁵ [1], and we assume that dA(f_EP) ≪ A(f_EP). In this case, equation (3) becomes

$$\begin{align}\hfill {M}_{\mathrm{d}}({f}_{\text{EP}})& \approx -\frac{\mathrm{d}A({f}_{\text{EP}})}{G({f}_{\text{EP}})A({f}_{\text{EP}})}E({f}_{\text{EP}})+A({f}_{\text{EP}})\left[{\delta }_{1}({f}_{\text{EP}})-{\delta }_{2}({f}_{\text{EP}})\right]\hfill \\ \hfill & \quad +{n}_{1}({f}_{\text{EP}})-{n}_{2}({f}_{\text{EP}}).\hfill \end{align} \tag{ 4 }$$

The impact of glitches Δ, is therefore simply bound by

$$\begin{equation}\vert {\Delta}\vert =\left\vert \frac{\mathrm{d}A({f}_{\text{EP}})}{G({f}_{\text{EP}})A({f}_{\text{EP}})}\pi ({f}_{\text{EP}})\right\vert \leqslant \left\vert \frac{\mathrm{d}A({f}_{\text{EP}})}{G({f}_{\text{EP}})A({f}_{\text{EP}})}E({f}_{\text{EP}})\right\vert ,\end{equation} \tag{ 5 }$$

where π(f) is the frequency signature of glitches (E = π + other external accelerations) and we assumed no cancelation between different external accelerations to provide the upper limit on the impact of glitches.

Before investigating it in the following section, we can note that up to the noise, the acceleration measured at the drag-free point (equation (1)) provides E(f_EP)/G(f_EP), when assuming |δ₁(f_EP)| ≪ |E(f_EP)| and |p₁(f_EP)| ≪ |E(f_EP)|. Then, an estimate of dA(f_EP)/A(f_EP) is enough to provide an estimate of the impact of the total external accelerations.

Equation (5) shows that the impact of glitches on the test of the WEP is related to the power of glitches at the test frequency f_EP. Before specialising to the MICROSCOPE case, we provide a rule of thumb analysis of what the frequency signature of glitches may look like.

We assume that the signal coming from glitches is some distribution of features external to the instrument (section 3.3) and sharing a similar shape localized in time (the mth glitch starting at time t_m ):

$$\begin{equation}{s}_{\mathrm{g}}(t)=\sum\limits _{m}{g}_{m}(t).\end{equation} \tag{ 6 }$$

For simplicity, we assume that all glitches have the same shape χ but distinct amplitudes a_m (this is indeed the case, see section 3.5), such that the mth glitch shape is observed in the measured time series as

$$\begin{equation}{g}_{m}(t)=({a}_{m}{\chi }^{\ast }h)(t-{t}_{m}),\end{equation} \tag{ 7 }$$

where ∗ denotes the convolution product and h(t) is the transfer function of the measuring apparatus (including the sensors and the DFACS).

The shape χ(t) is a priori unknown, as is the exact transfer function, preventing us from recovering χ. However, we can rewrite equation (7) as

$$\begin{equation}{g}_{m}(t)=({a}_{m}\delta \ast k)(t-{t}_{m}),\end{equation} \tag{ 8 }$$

where δ is the Dirac function and k is some kernel that provides the observed shape of glitches (albeit we do not have access to their real shape). We must emphasize that k(t) is not the transfer function, but is a mathematical description that allows us to investigate the effect of glitches on the measurement and data analysis.

Recalling that the power spectral density (PSD) of a stationary function u is, according to Wiener–Khintchine theorem,

$$\begin{equation}{S}_{u}(f)\equiv \mathcal{F}\left\{E[u(t)u(t+{\Delta}t)]\right\}(f)\end{equation} \tag{ 9 }$$

$$\begin{equation}=E[\vert \mathcal{F}\left\{u\right\}(f){\vert }^{2}],\end{equation} \tag{ 10 }$$

where $\mathcal{F}\left\{u\right\}$ is u's Fourier transform, and E[u(t)u(t + Δt)] is the signal's two-point correlation function (E denoting the expectation value). Glitches are assumed to follow a stationary random process of order N, where N is the number of points corresponding to some periodicity. Over one period, their averaged PSD is

$$\begin{equation}{S}_{\mathrm{g}}(f)=\mathcal{F}\left\{E[{\Xi}(t){\Xi}(t+{\Delta}t)]\right\}(f)\times \vert \tilde{k}(f){\vert }^{2},\end{equation} \tag{ 11 }$$

where Ξ(t) = ∑_m a_m δ(t − t_m ), $\tilde{k}(f)\equiv \mathcal{F}\left\{k\right\}(f)$ and we used both forms of the PSD to better highlight the two competing effects of glitches in the frequency domain: their time distribution and their shape. Their time distribution (modulated by their amplitude) enters in the PSD as the Fourier transform of their two-point correlation function while their shape is encoded in the Fourier transform of the kernel k, those two contributions being multiplied in the frequency domain. The kernel can therefore strongly impact the observed effect of glitches.

As an illustration, we can assume that glitches come as a Dirac comb ${\mathcal{I}}_{T}(t)$ of period T and that their observed shape is an exponentially damped sine k(t) = Θ(t)sin(2πf_k t)exp(−t/τ), where Θ is Heaviside step function; as we show below, these assumptions are not far from reality. In this case, the PSD of glitches is a Dirac comb of period 1/T multiplied by a Lorentzian (figure 2),

Note that the central frequency of the Lorentzian, its maximum value and its full width at half maximum depend on the kernel's parameters f_k and τ.

Zoom In Zoom Out Reset image size

In this case, glitches can therefore impact the WEP analysis if

• The power of the Dirac comb is significant at the test frequency f_EP (e.g. if it has a period 1/f_EP).
• And the glitches kernel does not bring the Dirac comb below the noise level at f_EP.

In the example of figure 2, if the dashed line (lower panels) shows the noise level, glitches would affect only the search for a signal of frequency f = k × 0.02 Hz ($k\in \mathbb{N}$) and lower than 0.22 Hz.

In more realistic cases, glitches can come randomly in time, or as sums of approximate Dirac comb (i.e. affected by some jitter), each one with its own amplitude. Some combinations of different periods may create an impact at f_EP, even if no Dirac with period 1/f_EP exists. Nevertheless, we expect that the kernel of glitches originating from similar processes is constant; it may be that different kinds of glitches exist, each with its own kernel. Such realistic cases should be treated numerically and taylored to the actually measured data, as done in the next section.

We shall summarise this section by emphasising that the effect of glitches on the measured PSD comes down to the multiplication of the Fourier transform of their time distribution by the Fourier transform of their observed shape. It is then important to be able to quantify their periodicities and model their shape.

Glitches can be seen in the acceleration data of all sessions and sensors, albeit with a signal-to-noise ratio (SNR) that obviously depends on the local noise level. They occur at a rate ranging from 0.02 s⁻¹ to 0.06 s⁻¹: thus, between 14 000 and 42 000 glitches can be detected in a 120-orbit session. We show in section 3.4 that the geographical and time distributions (the latter being modulated by the orbital period and the satellite's spin rate) of glitches hint at an Earth albedo origin with heating driving the satellite's coating to crackle. However, we have data only in a limited period of the year because of operational constraints; this makes it difficult to reach a robust conclusion, see also appendix A.

Glitches are undersampled: depending on their SNR, they are visible from half a second (two data points) to 5 s (20 data points). As a consequence, only a statistical analysis can shed light on their shape. A principle component analysis allowed us to show that up to their sign, all glitches have the same shape (section 3.5.2), this shape being the instrument's transfer function: glitches are thus shorter than the instrument's response time. However, it is not possible to invert the transfer function to know what the actual excitation looks like at the entry of the system.

Finally, as the test of the WEP is performed along the most sensitive axis of the sensors (x-axis), all accelerations shown in this article are along this axis. Nevertheless, glitches are visible along all axes simultaneously, though their amplitude is usually higher along one or the other. Note that we discuss masking glitches in section 4.1: depending on how conservative masks are (i.e. lasting from 1 s to 10 s), 4% to 50% of the data may be masked.

Figure 3 shows the acceleration of the inner sensor of the SUEP instrument for two sessions, in the time and frequency domains. In those sessions, the satellite's spin rate f_spin = 35/2f_orb which entails a rotation of the satellite on itself of period T_spin = 340 s. Zooms are provided in insets, which last 25 000 s (approximately five orbits). A glitch is also shown in the upper left panel. From this figure, it is clear that even under identical instrumental conditions (albeit the external experimental conditions may vary, e.g. due to the solar weather), the measured acceleration, as well as number and amplitude of glitches, can significantly vary (note that the acceleration range of the time domain is the same for both sessions). Note however that, despite the impression given by the graphes in time domain (upper panels) the difference of noise in the two sessions is not drastic: as shown by the spectra (lower panels) there is about a factor two, and only for high frequencies. Reference [22] discusses the evolution of the noise during the mission.

Zoom In Zoom Out Reset image size

Missing and invalid data (unrelated to glitches) represent only a few data points per session, and their effect is circumvented with the techniques described in [15–19]. In the remainder of this paper, unless stated otherwise, we fill gaps and replace invalid data by a local average of the data. They are so rare that this technique does not bias the results.

Glitches are detected in two steps. First, we perform a recursive σ-clipping (e.g. reference [23]) to extract outliers from the measured acceleration. Second, we segment the acceleration to gather outliers that belong to a single glitch (for instance, in the glitch shown in the inset of figure 3's upper panel, the outliers at t = 4.5 s and t = 5.5 s clearly belong to the same outlier).

Glitches are events external to the instrument, since they appear simultaneously and with similar amplitudes on all four test-masses, as can be seen in figure 4. Otherwise, there would be no reason to systematically see glitches on all sensors at the same time; indeed, the electronics are independent, and the gravitational interaction between sensors is too weak to explain how a sudden motion of a given sensor would affect the others [24].

Zoom In Zoom Out Reset image size

We now look for time and geographic periodicities in the occurrence of glitches. Evidently, such periodicities are linked since MICROSCOPE's orbit is circular, such that the satellite moves with a constant speed around the Earth; nevertheless, it is instructive to investigate them in parallel.

3.4.1. Geographical distribution

We start by investigating how glitches are distributed about the Earth. Note that we do not aim to investigate a potential systematic geographic distribution here. Some hints for a seasonal variation of the geographic distribution are given in appendix A, where 'local' effects (since each side of the satellite has its own propension to crack) are erased by averaging over sessions. Instead, we aim to visualise all effects (local and geographical). Thus, we bin the surface of the Earth in latitude and longitude and compute the mean number of glitches per session in each bin with a grid of resolution 30 × 40 pixels for individual sessions (to avoid erasing local effects when averaging over different sessions). The resolution is enough to see the effects of the satellite spinning as it orbits the Earth, even for the highest spin rate.

Figure 5 shows the mean number of glitches for two sessions with different spin rates: f_spin = 35/2f_orb (T_spin = 340 s) for the left panel, and f_spin = 9/2f_orb (T_spin = 1322 s) for the right panel. A geographic dipole appears in the former (with glitches more likely to occur in the southern hemisphere), modulated by faint stripes corresponding to the distance travelled by the satellite during a rotation; the dipole is more difficult to see in the right panel, where stripes corresponding to the satellite's rotation are instead clearly visible. These geographical distributions hint to two periodicities: one linked to the orbital period and a second one related to the satellite spin rate. We checked that the same patterns appear in all accelerometers, hinting to an external cause for glitches.

Zoom In Zoom Out Reset image size

3.4.2. Time distribution

Glitches follow an approximate Poisson process, albeit with a periodic clustering. Although this process is found to be stationary on the timescale of a given session, its mean rate evolves from session to session. Computing the time distribution of glitches is similar to computing their clustering, which we quantify by means of the two-point correlation function w(Δt): it gives the relative probability of finding a pair of glitches separated by a certain time interval Δt with respect to that of a homogeneous Poisson distribution. We use the estimator

$$\begin{equation}w({\Delta}t)=\frac{\mathrm{D}\mathrm{D}({\Delta}t)}{\mathrm{R}\mathrm{R}({\Delta}t)}-1,\end{equation} \tag{ 13 }$$

which is just a normalised version of the estimator used in section 2. In this expression, DD and RR are the number of data–data and random–random pairs in a time bin Δt, where random glitches are drawn from a homogeneous Poisson law within the same time-span as the data. Were glitches distributed uniformally, w(Δt) = 0; any excursion from w(Δt) = 0 highlights a preferred scale.

We estimate the uncertainty on the two-point correlation function w(Δt) for a given session with the method presented in reference [25]: we divide the session into N subsamples of equal length (long enough to include a significant number of glitches), on which we compute the correlation function w_i (Δt). From those N separate estimations, we calculate the mean $\bar{w}({\Delta}t)$ and error on the mean ${\Delta}\bar{w}({\Delta}t)$:

$$\begin{equation}{\left[{\Delta}\bar{w}({\Delta}t)\right]}^{2}=\frac{1}{N}\sum\limits _{i=1}^{N}{\left[\bar{w}({\Delta}t)-{w}_{i}({\Delta}t)\right]}^{2}.\end{equation} \tag{ 14 }$$

Figure 6 shows the two-point correlation function of glitches appearing in the SUEP internal sensor's acceleration for the sessions whose geographical distribution of glitches are shown in figure 5. For the sake of clarity, we restrict the range of the plot to a few orbit lengths; we checked that the same patterns is repeated all along the sessions durations. In this figure, the space between vertical red solid lines is equal to the orbital period T_orb, while that between vertical orange dashed lines shows the satellite's spin period (T_spin = 340 s in the left panel, T_spin = 1322 s in the right panel). Note that error bars (light grey) exclude 0 only close to time intervals corresponding to multiples of those periods: the correlation between other time intervals is negligible. The correlation function is computed in bins small enough to accommodate ten points between consecutive vertical dotted line.

Zoom In Zoom Out Reset image size

Two modulations are visible in both sessions: one at (approximately) the orbital period, and a second one at the spin period. They are intertwined, in the sense that the former is rather at twice the orbital period. This comes from the fact that the spin rate is a half-integer of the orbital frequency, meaning that the satellite–Earth system comes back in the same configuration (attitude vs latitude) every other revolution.

These observations hint towards correlations between the occurence of glitches with both the satellite's position on its orbit and its orientation with respect to the Earth. They concur with the conclusions drawn from the geographical distribution, and support the hypothesis that glitches originate from crackles of the multi-layer insulator (MLI) coating of the satellite, with one side more sensitive than others. The exact source of excitation of the MLI is not well understood, but may be related to the Earth albedo and/or to the Sun illumination. The hypothesis that crackles of the MLI cause glitches is also supported by our experience with other missions: similar transients affected GRACE but GOCE data was exempt from them; GRACE had an MLI coating but GOCE did not. Most importantly, crackles were identified in several types of MLI during pre-launch test, and the MLI with least crackles was selected [1].

We saw in section 2 (equation (11)) that both the time distribution and the shape of glitches affect their frequency content. The former as the Fourier transform of their two-point correlation function, and the latter as the signature of the kernel transforming an impulse into an observed glitch (recall that this kernel is related to the instrument's transfer function, though does not represent it). We discuss those two contributions here.

3.5.1. Frequency distribution

We first compute the Fourier transform of equation (11)'s Ξ function. For a given measurement session, we extract glitches from the acceleration measured by one sensor, retrieve their amplitude, and define Ξ(t) as the time-series consisting of Dirac impulses at the positions of glitches, with their corresponding amplitude. Figure 7 shows the Fourier transform of the two-point correlation of Ξ. The inset provides a zoom about the spinning and WEP test frequencies f_spin and f_EP.

Zoom In Zoom Out Reset image size

Consistently with our time-domain analysis above, we notice that glitches provide power at the orbital and spinning frequencies. The combination of those two periods brings power to other frequencies, distributed around f_spin and its harmonics, creating a forest of groups of spectral lines, mainly visible between a few mHz and 0.1 Hz. Noticeably, the WEP frequency is affected, with significant power brought up by the time distribution of glitches.

3.5.2. Kernel shape

The kernel k introduced in equation (11) represents glitches. Thus, they allow us to estimate it, though on a limited frequency range limited by the length of glitches (i.e. the observed shape gives no information on the behaviour of the kernel at frequencies smaller than 5 × 10⁻² Hz).

We performed a principal component analysis (e.g. reference [26]) on all glitches of SNR > 3, for several measurement sessions. Only two components were found to be significant, each having the sign opposite to the other's. Therefore, all glitches have the same shape, multiplied by their relative amplitude. The left panel of figure 8 shows the average shape of glitches in the time domain. After a sudden increase, it oscillates but quickly decreases, recalling an exponentially damped sinusoid. The error bars correspond to the error on the mean.

Zoom In Zoom Out Reset image size

The right panel of figure 8 shows the Fourier transform of this kernel. Although averaging glitches removes most of the noise, residual noise affects the Fourier transform, as does the limited number of data points. To bypass those limitations, we fit the average glitch with exponential shapelets [27]: this procedure allows us to denoise the average glitch and resample it at will, and to obtain the smooth Fourier transform shown in figure 8. We can note that this kernel, though peaking about a central frequency, has a more complex shape than a Lorentzian. Glitches are thus more complicated than pure damped sines (equation (12)).

3.5.3. Full frequency content

As shown by equation (11), the frequency signature of the glitches is given by the product of the signature of their distribution with the kernel investigated above. We should emphasise that this kernel contains the contribution of the electrostatic sensor and of the drag-free system and covers the entire frequency range from [10⁻⁵–2] Hz. Thus, the shape shown above by no means represents its integrality. However, no in-flight experiment can give us access to the kernel's behaviour at frequencies lower than 5 × 10⁻² Hz, preventing us from directly quantifying the effect of glitches over the entire frequency range. Thence, we must resort to a numerical model of the transfer function of the instrument.

From equations (7) and (8), the (normalised and denoised) observed shape of glitches is given by χ_obs = χ_true ∗ h = δ*k, where we added the subscript 'true' to emphasise that this function is the shape of the glitch at the entry of the system (though sampled at 4 Hz). Noting that $\mathcal{F}\left\{\delta \right\}=1$, we can estimate the 'true' shape of glitches as ${\hat{\chi }}_{\text{true}}(t)={\mathcal{F}}^{-1}\left\{\tilde{k}(f)/\tilde{h}(f)\right\}(t)$, where $\tilde{k}$ and $\tilde{h}$ are the Fourier transforms of the kernel k and of the instrument's transfer function h, restricted to frequencies available from the glitches shape (0.05 Hz ⩽ f ⩽ 2 Hz). This frequency restriction is prejudicial, and prevents us from recovering the correct 'true' shape of glitches.

Although most glitches presumably originate from crackles of the satellite's MLI, some of them also come from crackles (called clanks in reference [1]) of the gas tanks as their pressure decreases while the gas is consumed. Pre-launch tests showed that those clanks are extremely short events, lasting about a few dozens milliseconds [1]. Since they cannot be singled out in MICROSCOPE data, we conclude that all glitches are extremely short events, much shorter than the response time of MICROSCOPE's electronics. Thus, they can be considered as Dirac impulses when seen in MICROSCOPE's 4 Hz data. In the remainder of this paper, we assume that their true shape is a Dirac impulse, so that we can compute the effect of glitches at all frequencies, from equation (11), as

$$\begin{equation}{S}_{\mathrm{g}}(f)=\mathcal{F}\left\{E[{\Xi}(t){\Xi}(t+{\Delta}t)]\right\}(f)\vert \tilde{h}(f){\vert }^{2},\end{equation} \tag{ 15 }$$

where we also make use of the transfer function model (valid at all frequencies) and do not rely on the empirical kernel anymore.

As shown in section 2, in the linear model considered thus far, the transfer function consists of the contribution of the electrostatic sensor and of the DFACS, such that the transfer function of the mass controlling the drag-free is

$$\begin{equation}\tilde{h}(f)=\frac{A(f)}{1+G(f)A(f)},\end{equation} \tag{ 16 }$$

where we take the drag-free and sensor transfer functions from references [1, 22]. The transfer function is shown in figure 9, and shows a good agreement with the kernel $\tilde{k}$ restricted to the frequencies available from the limited size of glitches.

Zoom In Zoom Out Reset image size

A standard analysis of some measurement sessions provides a statistically unlikely (given the estimates of the Eötvös parameter on other MICROSCOPE measurement sessions [28]) level of a WEP violation. One in particular (session 380—see reference [20] for a list of sessions) yields a 7σ detection on the reference instrument (the test masses of which are made of the same material [2]). This detection hints towards a systematic effect unaccounted for and made apparent by the lower-than-average noise level of that particular session [28]. Additionally, a forest of strong lines can be seen in the differential acceleration measured for this session (black line of figure 10).

Zoom In Zoom Out Reset image size

These spectral lines, as they look alike the glitches' spectral distribution (figure 7), make it tempting to attribute this signal to glitches. To investigate this possibility, we mask glitches and reconstruct the masked data with the modified-expectation-conditional-maximization (M-ECM) algorithm [18]. M-ECM maximises the model likelihood with respect to the available data by estimating the missing data conditional expectation, based on the circulant approximation of the complete data covariance. For this exercise, we look for glitches not only on the x-axis, but also on the y- and z-axes, to make sure that glitches preferentially projected on one of those latter axes are detected, though buried in the noise of the x-axes. Detecting glitches with the σ-clipping technique described above (using a 4.5σ threshold), and masking ten seconds after each outlier, we remove 46% of the data. This shows that glitches appear in considerable amount and may significantly contribute to the measurement noise. Figure 10 compares the spectra of the raw differential acceleration and of the differential acceleration after discarding glitches and reconstructing the data with M-ECM. The line forest between 0.01 Hz and 0.1 Hz is significantly reduced, meaning that glitches contribute to it. Some lines are nonetheless still present, hinting at the presence either of lower SNR glitches or of another unidentified contributor.

No significant WEP violation can be detected in the reconstructed data. The aforementioned 7σ WEP detection when not masking glitches is therefore directly related to their presence. We checked that the M-ECM-reconstructed data would allow us to accurately estimate a WEP violation in this masked data by adding a mock WEP signal to the original (before masking) data, and correctly estimating it.

Equations (1) and (5) readily provide an upper bound on the effect of glitches from the acceleration measured at f_EP by the sensor at the drag-free point (right panel of figure 12), when assuming a linear model for the instrument. Indeed, at f_EP, the acceleration of the drag-free sensor is M_df = E/GA ≈ 10⁻¹³ m s⁻². Assuming moreover that dA is of order the measured scale factor difference |dA| ≈ 10⁻² [8], the impact of external forces (including glitches) on the differential acceleration is ≈10⁻¹⁵ m s⁻². In terms of the Eötvös parameter, this amounts to about 1.2 × 10⁻¹⁶, two orders of magnitude below the 1σ uncertainty, well below what is actually observed.

This result can also be seen graphically. Figure 11 shows the measured differential acceleration (grey) and the differential signal created by glitches, still assuming |dA| ≈ 10⁻². The contribution of glitches at f_EP, in terms of the Eötvös ratio, can be seen to be at the level of 10⁻¹⁶, two orders of magnitude below the noise.

Zoom In Zoom Out Reset image size

Nevertheless, our model recovers the forests of peaks, though with an amplitude lower than observed, except around the f_EP frequency. This can be seen in figure 12, which shows the frequency content of glitches (equation (15)) expected from our linear model for both sensors (red), compared with the measured acceleration (grey). The left panel corresponds to the out-of-drag-free sensor (thence, the frequency-dependent low-frequency noise) and the right panel to the drag-free sensor (thence, the flat low-frequency noise); at low frequency, the higher noise level from glitches is due to a conservative assumption about a stochastic background of low SNR glitches.

Zoom In Zoom Out Reset image size

In summary, the model is satisfactory at the level of individual sensors' acceleration (especially about f_EP), but seems unable to reproduce the signal created by glitches in the differential acceleration. As discussed below, the explanation of this failure can be twofold:

• (a)  
  A linear model is too simple and not representative of the complexity of a real system, especially with a measurement system so sensitive that any perturbation is easily spotted
• (b)  
  The distribution of glitches input to the model is incorrect.

4.2.1. Non-linear model and internal saturations

4.2.2. Glitches distribution

4.2.3. Conclusion