blue span represent respectively the excess $S_{\Delta g,e}^{1/2}$ and the
total modeled noise, as in Fig. 22. Adding deviations of
$1\text{\,}\mathrm{\SIUnitSymbolMicro V}$ and
$2\text{\,}\mathrm{\SIUnitSymbolMicro V}$ RMS, respectively, the total modeled
noise becomes the one represented by the green and yellow spans.

Potentially, an inaccuracy of $2\text{\,}\mathrm{\SIUnitSymbolMicro V}$ RMS,
or even less, could result in non-negligible force noise, explaining a
relevant fraction of the detected excess noise. At the same time, it would go
undetected in in-flight calibration experiments. Preliminary measurements on
LPF-prototype FEE models show that the presence of residual nonlinearities of
the magnitude above can not be excluded, i.e., measurements are compatible
with the presence of a nonlinearity of $1\text{\,}\mathrm{\SIUnitSymbolMicro
V}$ RMS. However, we are currently planning deeper and more systematic tests,
to better characterize nonlinearities and assess their impact.

FIG. 23:  Total modeled noise, including potential inaccuracies in the applied
actuation voltage amplitudes. Red points and blue span: respectively, excess
noise $S_{\Delta g,e}^{1/2}$ and total modeled noise, as in Fig. 22. Green and
yellow spans: modeled noise, summed to voltage amplitude nonlinearities with
$1\text{\,}\mathrm{\SIUnitSymbolMicro V}$ and
$2\text{\,}\mathrm{\SIUnitSymbolMicro V}$, respectively. The additional PSD
due to such inaccuracies is the joint posterior distribution of 10 independent
realizations with the given RMS, as described in the text.

Another mechanism through which voltage disturbances could lead to increased
in-band force noise and, at the same time, go undetected in dedicated
measurement campaigns, is through down-conversion of high-frequency spurious
disturbances. If voltage anomalies should be present at frequencies outside
the measurement band, they would down-convert into the band because of the
quadratic nature of the electrostatic force. In addition, due to the lack of
significant associated torque, the voltage anomaly should have involved one or
more pairs of electrodes facing the same face of the same TM.

To fix the amplitude scale of such a disturbance, the force due to the voltage
$V$, when applied to both electrodes in one such pair, exerts a force
$\lvert\Delta g\rvert=(1/M)\lvert\partial C_{x}/\partial x\rvert V^{2}\simeq
0.15\,\text{pm\,s}^{-2}\,(V/100\,\text{mV})^{2}$. Thus a line with a mean
amplitude of 100 mV, and a noisy relative amplitude fluctuation with ASD
$\simeq 5\times 10^{-3}/\sqrt{\text{Hz}}\,(1\,\text{mHz}/f)$ would produce the
observed noise. We note that the actuation circuitry connected to the
electrodes includes a passive two-stage lowpass filter with a bandpass near 2
kHz, thus this putative high-frequency disturbance would have needed to be
accordingly higher at the amplifier outputs for frequencies larger than 2 kHz.

No such line has been detected either during the ground testing of the FEE
flight hardware, or during laboratory testing of its various prototypes, but
we cannot exclude for instance some damage due to the launch stresses.

###  IX.3  Unmodeled gravitational noise

We have discussed the gravitational noise due to propellant tank depletion and
that due to LTP distortion. We have accounted for the first and estimated that
it is unlikely that the second may have caused more than 10% of the
unaccounted noise.

Besides those effects, any other mass motion, either because of distortion of
solid parts or because of evaporation of volatile fractions, may cause
gravitational force noise and may have contributed to excess noise. We discuss
here a few possibilities, first for the case of distortion and then for that
of evaporation.

The Tungsten balance mass is the most intense source of gravitational field
gradient at the TM location. The LTP distortion, already discussed, moves the
balance mass, relative to the TM, _together with the entire GRS_. However, in
addition to that, any _internal GRS distortion_ may also have moved the
balance mass relative to the TM.

With a gravitational gradient $\simeq\partial g_{x}/\partial x_{s}\simeq
5\times 10^{-7}\,\text{s}^{-2}$, a random jitter $\delta x$ of the balance
mass position with ASD $S_{\delta
x}^{1/2}\simeq$1.5\text{\,}\mathrm{n}\mathrm{m}\mathrm{/}\sqrt{\text{Hz}}$\,\left($1\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$/f\right)$,
would account for the entire unexplained excess. Such a jitter, if thermally
induced, given the construction materials and the geometry of the GRS, would
amount to thermal fluctuations with an ASD of
$S_{T}^{1/2}\simeq$4\text{\,}\mathrm{m}\mathrm{K}\mathrm{/}\sqrt{\text{Hz}}$\,\left($1\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$/f\right)$,
definitely larger than the measured one [8].

As for non-thermal deformation, an obvious example would be long-term, noisy
mechanical secondary creep due to stress release, like that due to the unlock
of the TM on orbit. Crudely approximating this creep as a Poisson sequence of
steps with root mean square (rms) amplitude $\delta x$, and average relaxation
rate $\delta\dot{x}$, one would need $\delta
x\delta\dot{x}\simeq$0.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\mathrm{/}\mathrm{yr}$$.
As realistically
$\delta\dot{x}\ll\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{2}\mathrm{/}\mathrm{yr}$$,
this would require steps of $0.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
happening then at a rate of less than 5 per year, very different from the
observed time series.

Thermal or non-thermal distortion resulting in the motion of massive
components farther away from the TMs may also exert significant forces.
Calculations [47] show for instance that the spacecraft alone, without the
LTP, exerts a static difference of force on the TM of $\Delta g_{sc}\simeq
5\times 10^{-9}\,$\mathrm{m}\,\mathrm{s}^{-2}$$, very similar to that of a
homogeneous square toroid with an inner diameter of one meter, an outer one of
two meter, and a mass of $\simeq$300\text{\,}\mathrm{k}\mathrm{g}$$, a crude
approximation to the spacecraft shape.

One can calculate that a homogeneous relative distortion $\delta$ of any of
the toroid dimensions, causes a variation of differential force
$\delta\times\Delta g_{sc}$. It would thus take $\delta$ to be a random
process with ASD $S_{\delta}^{1/2}\simeq 1.4\times
10^{-7}\,$\mathrm{Hz}^{-1/2}$\,($1\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$/f)$
to justify the unexplained noise.

We have no way of assessing if such mechanical distortions of the spacecraft
(about $0.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ root-mean-square variation
of the corresponding dimension over 1 day) did take place during operations. A
correlation analysis of $\Delta g$ with the solar power hitting the spacecraft
gave no significant results. A simulation of thermal induced distortion and of
the resulting gravitational noise performed during mission development [47],
pointed to one order of magnitude smaller figures. This simulation was done
though on a simplified model and was significant only for
$f\gtrsim$2\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$$.

All that said, detailed thermo-mechanical simulations are standard practice in
aerospace engineering. In addition, for LPF rather complete tools have been
set up to interface the thermo-mechanical model of the system with the
gravitational one. We calculate that a thorough dynamical simulation of the
LISA spacecraft’s gravitational field should have the accuracy and the
sensitivity to keep gravitational noise due to mechanical distortion under
control.

Outgassing of volatile fractions from spacecraft has been observed, for
instance in the Rosetta mission, to continue years after launch [48]. In
particular, after a water desorption-dominated phase lasting 100-200 days, a
longer, possibly diffusion-dominated phase was observed, with a very slowly,
if at all decaying evaporation rate, at least for the following 500 days.

Given the non-symmetric geometry of the LPF spacecraft, even a spatially
homogeneous _noisy_ outgassing would have caused some $\Delta g$ noise.
Definitely more so if the outgassing had also some patchy pattern, that would
have further reduced the spatial symmetry.

To fix the scale of the effect, we have considered the spacecraft structure,
which is made out of sandwich panels of Carbon Fiber Reinforced Plastic (CFRP)
with an aluminum honeycomb core. We have considered the case of N outgassing
patches of size small enough to approximate sums with integrals, having noisy
evaporation rates all with same ASD, incoherently adding up to a total
evaporation rate with ASD $S_{\dot{M}}^{1/2}$.

We calculate that $S_{\Delta g}^{1/2}\simeq
0.5\left(\mathcal{G}/L_{0}^{2}\right)S_{\dot{M}}^{1/2}/(2\pi f)$, with
$L_{0}=37.6\text{ cm}$ the distance between the centers of the two TMs, and
$\mathcal{G}$ the gravitational constant. To match our observation of
$S_{\Delta
g}^{1/2}\simeq$0.7\text{\,}\mathrm{f}\mathrm{m}\,\mathrm{s}^{-2}\mathrm{/}\sqrt{\text{Hz}}$\,(1\,\text{mHz}/f)$,
we need a frequency independent value
$S_{\dot{M}}^{1/2}\simeq$0.02\text{\,}\mathrm{m}\mathrm{g}\,\mathrm{s}^{-1}\mathrm{/}\sqrt{\text{Hz}}$$.

Just for the sake of discussion, we note that such a “white” evaporation rate
behavior would naturally be obtained if the outgassing on LPF consisted of a
Poisson succession, with mean evaporation rate $\langle\dot{M}\rangle$, of
discrete outgassing events from any of the patches, with a rms value $\delta
m$. This process would indeed give $S_{\dot{M}}^{1/2}=(\delta
m\,\langle\dot{M}\,\rangle)^{1/2}$.

Ref. [48] estimates the mass loss of Rosetta in several hundred grams per
year. LPF was a lighter satellite, by almost a factor 3, but with more plastic
and closer to the sun. Thus a direct projection of the outgassing properties
may be rather speculative.

Nevertheless, a crude estimate based on the outgassing properties of LPF
materials, gives a lower limit for $\langle\dot{M}\rangle$ at
$\langle\dot{M}\rangle\simeq$0.1\text{\,}\mathrm{k}\mathrm{g}\,\mathrm{y}^{-1}$=$3\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{g}\,\mathrm{s}^{-1}$$,
not far from what one would extrapolate from Rosetta based just on the mass
ratio of the two satellites. Using this figure for $\langle\dot{M}\rangle$ in
the Poisson noise ASD formula, one gets $\delta
m\simeq$0.2\text{\,}\mathrm{m}\mathrm{g}$$ and
$\lambda=\langle\dot{M}\rangle/\delta m\simeq$0.01\text{\,}\mathrm{s}^{-1}$$.

To our knowledge, there are no studies on the spatial distribution of
outgassing in spacecraft, nor on its dynamics. Thus this crude scenario of
noisy outgassing remains rather speculative. We note however that nothing, in
the figures above, or in the current knowledge, allows us to rule it out. We
conclude that a cautious approach for LISA would be to stay as close as
possible to the LPF design in terms of quantity and distribution of materials
with significant volatile components.

###  IX.4 Pressure fluctuation

As said in Sec. V.6, the complex geometry of the TM environment may create
quasi-static pressure gradients. Any in-band fluctuation of such gradients
would directly translate into an in-band acceleration fluctuation.

To account for the measured excess $S_{\Delta g_{e}}^{1/2}$ of about
$0.7\,$\mathrm{f}\mathrm{m}\,\mathrm{s}^{-2}\mathrm{/}\sqrt{\text{Hz}}$\,($1\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$/f)$,
the ASD of these fluctuations should be $S_{\Delta p}^{1/2}\sim 0.6\times
10^{-12}\,$\mathrm{Pa}\mathrm{/}\sqrt{\text{Hz}}$\,($1\text{\,}\mathrm{m}\mathrm{H}\mathrm{z}$/f)$.

We note that the process at the basis of this hypothetical fluctuating
pressure should necessarily be different from that generating the static
gradient discussed in Sec. V.6. Indeed, (see Sec. IV) the latter steadily
decreased during the mission, following the general decay of the pressure,
whereas the former (Sec. V.2) did not.

Hence, if the process responsible for the $1/f$ noise has been due to the
outgassing of some species for the TM environment, it must have been
substantially stable over the entire mission and independent of the outgassing
rate of the main fraction. This observations parallels that on the possible
outgassing origin of force glitches in LPF [20] that also had properties that
have been stable over the course of the entire mission.

It does not require much outgassing to produce the $1/f$ noise. If for
instance such hypothetical gas phase were hydrogen diffusing out of the
various elements of the TM environment, and if the noisy outgassing took place
as the series of undetectable glitches discussed in Sec. VII.2 and in the App.
G, it would only take a mean outgassing rate of $\simeq 6\text{ pg/d}$, that
is some 3 ng of total emission over the course of the entire mission, to
explain the noise. We have no specific piece of experimental evidence to
support this hypothesis, neither we could trace any relevant study on pressure
fluctuations in vacuum systems. Neither we have, however, any evidence proving
the model false or unlikely, neither from our own experiments nor from the
literature. It is also unlikely that we could devise a laboratory experiment
able to detect such tiny pressure fluctuations.

For this reason, we conclude again that the LISA GRS outgassing environment
should be kept as close as possible to that of LPF, to retire the risk of
unwanted large pressure gradient fluctuations.

###  IX.5  High frequency magnetic field noise

In addition to low-frequency effects, discussed in Sec. VIII, magnetic fields
at high frequency may induce eddy currents within the test masses, and then
exert Lorentz forces on them [49]. The effect is thus quadratic and would
convert the low frequency amplitude fluctuations of a high frequency magnetic
spectral line into a corresponding low frequency force.

To give a scale of the effect, a recent finite-element electromagnetic
calculation by the LISA project [50], has shown that the effect of a dipole of
$1\text{\,}\mathrm{m}\mathrm{A}\,\mathrm{m}^{2}$ located at a distance
$d=$20\text{\,}\mathrm{c}\mathrm{m}$$ from the test mass and oscillating at
the frequency of 100 Hz, would cause a force of $\Delta
g\simeq$4\text{\,}\mathrm{f}\mathrm{m}\mathrm{/}\mathrm{s}^{2}$$. The effect
reaches its peak at $100\text{\,}\mathrm{Hz}$, while at lower frequency the
induced current decreases, and above that, the screening effect of the
metallic electrode housing attenuates the oscillating field. The effect of a
dipole source decreases with $d^{-7}$, so that at the closest distances of
about $0.4\text{\,}\mathrm{m}$ between the test mass and any active device on
the LPF spacecraft the effect might be $\sim 100$ times smaller.

The spacecraft prime contractor, during LPF development, performed a test
campaign on ground against audio frequency magnetic lines [51]. A few, barely
significant lines have been identified, with peak amplitudes
$\ll$1\text{\,}\mathrm{nT}$$ at the position of the test masses. In the point
dipole model at a distance of $d\simeq$0.4\text{\,}\mathrm{m}$$, a
$1\,\text{nT}$ line would be generated by a dipole of
$\simeq$0.3\text{\,}\mathrm{m}\mathrm{A}\,\mathrm{m}^{2}$$ at most, and would
exert a static force $\Delta
g\simeq$4\text{\times}{10}^{-3}\text{\,}\mathrm{f}\mathrm{m}\,\mathrm{s}^{-2}$$.
To reproduce the observed excess noise with a noisy amplitude modulation of
one of these lines $1\text{\,}\mathrm{nT}$, one would need a relative
amplitude modulation with ASD $\simeq 2.5\times
10^{2}\text{\,Hz}^{-1/2}\left(1\,\text{mHz}/f\right)$.

One can calculate that the rms fluctuation of a process $x$ with ASD given by
$S_{x}=S_{0}^{1/2}(f_{0}/f)$, over a data stretch of duration $T$, is $\langle
x_{\text{rms}}\rangle=S_{0}^{1/2}\pi f_{0}\sqrt{T/3}$. With $S_{0}^{1/2}\simeq
2.5\times 10^{2}\text{\,Hz}^{-1/2}$, it would take about 5 s to make the
relative rms fluctuation of the line amplitude 100%. Such a large fluctuation
would have been noticed during testing.

The required amplitude fluctuation for such lines would become smaller in the
presence of more than one line. However consider that 1 nT is a real
conservative upper limit for the observed lines, and that the effect is
quadratic in the amplitude of the lines, so that the combined effect of many
lines would be dominated by the few brightest among them.

Despite these considerations, that would rule out this source of noise, it
must be stressed that unfortunately we had no magnetometer on board sensitive
to the audio band. Thus, as the operating conditions may have been different
from those during testing, we cannot exclude that additional, more intense,
amplitude-modulated lines had been generated once on orbit.

For instance, a single line generated by a dipole of
$\simeq$2\text{\,}\mathrm{m}\mathrm{A}\,\mathrm{m}^{2}$$ located at a distance
$d=$20\text{\,}\mathrm{c}\mathrm{m}$$ from the test mass, with a more
reasonable relative amplitude fluctuation ASD
$\simeq$0.2\text{\,}\mathrm{Hz}^{-1/2}$\left(1\,\text{mHz}/f\right)$, would
explain the excess noise. This is a relatively large magnetic dipole, for
instance 32 mA in a fully uncompensated square loop of 25 cm size, of the
existence of which we are not aware. Nevertheless we certainly recommend that
in LISA thorough testing is performed on ground, and that on-board diagnostic
magnetometers with sensitivity up to, at least, $1\text{\,}\mathrm{kHz}$ are
seriously considered.

##  X Concluding remarks

LISA Pathfinder reached an acceleration noise performance achieving the LISA
requirements with margin, and better than both its requirements and what had
been estimated before launch [49].
This last achievement was mostly allowed by two facts:

  * •

A much better self-gravity cancellation than what had been very cautiously
estimated on ground. This in turn allowed using much less electrostatic force
authority than predicted [2], and suppressing accordingly the actuation noise.

  * •

The achievement, over the course of the mission, of a lower base pressure, and
thus a lower Brownian noise, than had been assumed in pre-launch estimates
[49].

Both these facts revealed the existence of an excess noise above the Brownian
noise level, with a $1/f^{2}$ PSD.
Based on this starting framework, in this paper we have:

  * •

shown that the Brownian noise evolved in agreement with the outgassing of a
single gaseous species diffusing out of the immediate, complex surroundings of
the TM, where such outgassing also maintains a quasi-static pressure gradient
across the TM;

  * •

shown that the temperature stability of the system was good enough that
temperature fluctuation played a significant role in the acceleration noise
only at the lowest analyzed frequency $18\text{\,}\mathrm{\SIUnitSymbolMicro
Hz}$, well below the LISA lower frequency of $0.1\text{\,}\mathrm{mHz}$;

  * •

shown that the intrinsic stability of the $1/f^{2}$ excess PSD (referred to in
the article as $1/f$ excess for ASD) was $\pm 20\%$ in amplitude over more
than 16 months, with the residual fluctuations being independent of any
traceable change in the operational conditions that were needed to run the
mission;

  * •

analyzed all sources of noise for which we had a verified model, either from
correlation analysis or from dedicated experiments, and concluded that these
sources account for a fraction $0.46\pm 0.09$ of the total power of the
excess;

  * •

finally, we have discussed all possible explanations we could trace for the
unaccounted part of the excess, with patch potentials, actuation electronics
nonlinearity, audio frequency magnetic fields, and pressure fluctuations as
the prominent ones, and identified possible measures to keep them under full
control during the implementation of LISA.

## Acknowledgments

This work has been made possible by the LISA Pathfinder mission, which is part
of the space-science program of the European Space Agency.
We thank Paolo Chiggiato and his team at CERN, for very helpful discussions
about the LPF outgassing environment.
The Italian contribution has been supported by Istituto Nazionale di Fisica
Nucleare (INFN) and Agenzia Spaziale Italiana (ASI), Project No.
2017-29-H.1-2020 “Attività per la fase A della missione LISA”. The UK groups
wish to acknowledge support from the United Kingdom Space Agency (UKSA), the
Scottish Universities Physics Alliance (SUPA), the University of Glasgow, the
University of Birmingham, and Imperial College London. The Swiss contribution
acknowledges the support of the Swiss Space Office via the PRODEX Programme of
ESA, the support of the ETH Research Grant No. ETH-05 16-2 and the support of
the Swiss National Science Foundation (Projects No. 162449 and No. 185051).
The Albert Einstein Institute acknowledges the support of the German Space
Agency, DLR. The work is supported by the Federal Ministry for Economic
Affairs and Energy based on a resolution of the German Bundestag (No. FKZ
50OQ0501, No. FKZ 50OQ1601, and No. FKZ 50OQ1801). J.I.T. and J.S. acknowledge
the support of the U.S. National Aeronautics and Space Administration (NASA).
The Spanish contribution has been supported by Contracts No. AYA2010-15709
(Ministerio de Ciencia e Innovación MICINN), No.ESP2013- 47637-P, No.
ESP2015-67234-P, No. ESP2017-90084-P (Ministerio de Asuntos Económicos y
Transformación Digital, MINECO) and No. PID2019-106515GB-I00 and
PID2022-137674NB-I00 (MCIN/AEI/10.13039/501100011033). Support from AGAUR
(Generalitat de Catalunya) Contracts No. 2017-SGR-1469 and 2021-SGR-01529 is
also acknowledged. It has also received partial support from the program
Unidad de Excelencia María de Maeztu CEX2020-001058-M (Spanish Ministry of
Science and Innovation). The French contribution has been supported by the
CNES (Accord Specific de projet No. CNES 1316634/CNRS 103747), the CNRS, the
Observatoire de Paris and the University Paris-Diderot. E.P. and H.I. would
also like to acknowledge the financial support of the UnivEarthS Labex program
at Sorbonne Paris Cité (No. ANR-10-LABX-0023 and No. ANR-11-IDEX-0005-02).
N.K. would like to thank for the support from the CNES Fellowship.

## Appendix

The following appendices contain some rather important information, tables,
and calculations.
In App. A, we give details about the analyzed noise runs.
In App. B, we describe our spectral analysis tools and decorrelation methods.
In App. C, we list the frequencies used in PSD estimation.
In App. D, we provide the ASD of $\Delta g(t)$ for all analyzed noise runs.
In App. E, we present the calculations for instrument distortion analysis and
the drift evaluation procedure.
In App. F, we present the calculations for in-band noise of instrument
distortion and thermal effects.
In App. G, we study whether a Poisson flow of low-amplitude glitches could
produce noise with the observed properties.

##  Appendix A  Experimental configurations for the 13 runs

There were minor differences in the operating conditions of the 13 runs of
Table 1. We describe them in the following and summarize the different
configurations in Table 2.

Run |  Propellant |  Thruster |  Actuation |  TM |  Voltage |  Heaters |  ST7
---|---|---|---|---|---|---|---
|  tank |  branch |  authority |  alignment |  compensation |  configuration |  state
1 |  2 |  A |  URLA |  1 |  1 |  1 |  OFF
2 |  2 |  A |  UURLA |  1 |  1 |  1 |  OFF
3 |  2 |  A |  UURLA |  1 |  2 |  1 |  OFF
4 |  3 |  A |  UURLA |  2 |  3 |  1 |  OFF
5 |  3 |  A |  UURLA |  3 |  3 |  2 |  DIAG
6 |  1 |  A |  UURLA |  3 |  3 |  2 |  DIAG
7 |  3 |  B |  UURLA |  3 |  3 |  2 |  DIAG
8 |  1 |  B |  UURLA |  3 |  3 |  3 |  DIAG
9 |  1 |  B |  UURLA |  3 |  4 |  4 |  OFF
10 |  1 |  B |  UURLA |  3 |  4 |  5 |  OFF
11 |  3 |  A |  UURLA |  3 |  4 |  6 |  OFF
12 |  3 |  A |  UURLA |  3 |  4 |  1 |  OFF
13 |  3 |  A |  UURLA |  3 |  4 |  1 |  OFF
TABLE 2: Experimental configuration for the 12 different runs of Table 1. The
meaning of the numeric labels is explained in the text.

  * •

Thruster propellant was stored in three different tanks. The gravitational
signal from the depletion of these tanks was different, due to the different
positions. In Table 2, we name these three tanks as 1, 2, and 3.

  * •

For redundancy reasons, LPF carried two independent branches of micro-
thrusters [13]. In Table 2, we call them A and B.

  * •

The electrostatic controllers needed a setting for the maximum force/torque
authority they could deliver within a linear regime [19]. These settings
determine the zero-actuation voltages commanded to the various electrodes. The
settings were identical for all the considered runs except for the first one.
Values are reported in Table 3.

|  URLA |  UURLA
---|---|---
$F_{x,1}$ [pN] |  0 |  0
$F_{x,2}$ [pN] |  50 |  50
$N_{x,1}$ [pN m] |  16.37 |  4
$N_{x,2}$ [pN m] |  16.37 |  4
$F_{y,1}$ [pN] |  3670 |  1000
$F_{y,2}$ [pN] |  3670 |  1000
$N_{y,1}$ [pN m] |  13.32 |  4
$N_{y,2}$ [pN m] |  13.32 |  4
$F_{z,1}$ [pN] |  5820 |  500
$F_{z,2}$ [pN] |  5820 |  500
$N_{z,1}$ [pN m] |  1.5 |  1.5
$N_{z,2}$ [pN m] |  1 |  1
TABLE 3: Maximum actuation authorities, $F_{a,i}$ and $N_{a,i}$, respectively
the maximum force and the maximum torque applicable on TM$i$ along axis $a$.

  * •

In the attempt to reduce cross-talk between the differential interferometer
readout and the motion of the remaining TM degrees of freedom [16], we
performed three different adjustments of the zero set points of all control
loops. The most relevant ones were those for the angular orientations of the
TM along $z$ and $y$. The three set points, named 1, 2, and 3, are listed in
Table 4.

Angle |  Set point |  Set point |  Set point
---|---|---|---
$($\mathrm{\SIUnitSymbolMicro rad}$)$ | 1 | 2 | 3
$\phi_{1}$ | -59.25 | -56.32 | -61.2
$\phi_{2}$ | -21.35 | -33.01 | -9.7
$\eta_{1}$ | -3.5 | -2.14 | -4.9
$\eta_{2}$ | 3.5 | 10.3 | -3.3
TABLE 4: The three different set points for the angular orientations of the
TMs along $z$ ($\phi$) and $y$ ($\eta$), used during the noise runs. The index
indicates which TM each angle refers to. Angles are relative to a reference
frame defined by the OMS.

  * •

Random TM charging due to cosmic rays couples to DC voltage differences
between the TM and the surrounding metal surfaces to produce random force and
torque [22]. This noise may be suppressed by purposely applying some DC
voltages to the various electrodes, in order to compensate the parasitic ones
that are found on metal surfaces because of work function differences [22]. In
particular, the force along $x$ and the torque along $\phi$, may be
compensated by a proper combination of voltages on the electrodes facing the
$x$-faces of the TMs. These compensation voltages may be described by just
voltage parameters $\Delta_{x}$ and $\Delta_{\phi}$ as explained in [22]. In
the course of the mission we have made some adjustments of these parameters
for both TMs. The various configurations are listed in Table 5.

Voltage |  Setting |  Setting |  Setting |  Setting
---|---|---|---|---
(mV) | 1 | 2 | 3 | 4
$\Delta_{x,1}$ | 0 | +24 | +12 | +24
$\Delta_{x,2}$ | 0 | 0 | 0 | 0
$\Delta_{\phi,1}$ | 0 | 0 | 0 | +32
$\Delta_{\phi,2}$ | 0 | 0 | 0 | -116
TABLE 5: Compensation voltages for the 4 adopted settings. The index indicates
the TM.

  * •

While the noise runs in Table 1 were all performed by using the satellite cold
gas thrusters for drag-free control [13], LPF also carried a set of
alternative thrusters based on colloidal propellant technology, in the
framework of NASA ST7 mission [52]. These thrusters were used intermittently,
leaving them, for purpose of diagnostics, in some activated state even when
not in use, in the two epochs from June 27, 2016 to December 7, 2016, and from
March 18, 2017 to April 29, 2017. In Table 2 we indicate such state as DIAG,
while the completely off state is indicated with OFF.

##  Appendix B Spectral estimation methods

###  B.1 Periodograms and their spectral properties

We use, for the elementary periodogram $X(k)$ of the $N$-sample series $x[n]$,
sampled with sampling time $T$, the standard definition:

$X(k)=\sqrt{\frac{T}{N}}\sum_{n=0}^{N-1}x[n]w[n]\,e^{-2\pi ikn/N},$ (18)

where $w[n]$ are the coefficients of a Blackman-Harris spectral window, which
gives good side-lobe suppression [53].

Following the Welch method [54], we section our data series in as many 50%
overlapping data stretches of length $N$ as they fit into the length of the
data series, and define the average, one-sided experimental PSD at frequency
$f$:

$\Pi\left(f=\frac{k}{NT}\right)=\frac{2}{M}\sum_{l=1}^{M}X_{(l)}(k)X_{(l)}^{*}(k),$
(19)

where the index $l$ runs over the $M$ data stretches, and $k$ is the frequency
index.

For multiple, synchronously measured data series $x_{i}[n]$, with $1\leq i\leq
p$, Eq. (19) generalizes to the (one-sided) experimental cross-spectral
density (CPSD) matrix $\bm{\Pi}$, with elements:

$\Pi_{ij}\left(f=\frac{k}{NT}\right)=\frac{2}{M}\sum_{l=1}^{M}X_{i,(l)}(k)X_{j,(l)}^{*}(k),$
(20)

a complex and Hermitian matrix. It is also positive-definite if $M\geq p$. Our
founding point is that, if the $x_{i}(t)$ are Gaussian, zero-mean stationary
stochastic processes, with a theoretical cross-spectral density matrix
$\bm{\Sigma}$, then the matrix $\bm{W}$, whose elements are
$W_{ij}=M\Pi_{ij}$, is distributed like a complex Wishart distribution [55,
56]:

$p\left(\bm{W}\big{|}\bm{\Sigma}\right)=\frac{\left|\bm{W}\right|^{M-p}}{\widetilde{\Gamma}_{p}(M)\left|\bm{\Sigma}\right|^{M}}\
\mathrm{etr}\left[-\bm{\Sigma}^{-1}\bm{W}\right]$ (21)

Where $\left|\cdot\right|$ is the determinant, $\mathrm{etr}$ the exponential
trace $\mathrm{etr}(\cdot)=\exp\left(\mathrm{tr}\left(\cdot\right)\right)$ and
$\widetilde{\Gamma}_{p}(M)$ is the multivariate complex Gamma function,
defined by

$\widetilde{\Gamma}_{p}(M)=\pi^{\frac{1}{2}p(p-1)}\prod_{i=1}^{p}\Gamma(M-i+1)$

We denote this distribution as $\mathcal{CW}(\bm{\Sigma},M)$, and indicate
that $\bm{W}$ follows this distribution as
$\bm{W}\sim\mathcal{CW}(\bm{\Sigma},M)$.

Note that in Eq. (21) it is required that $M\geq p$, otherwise the matrix
$\bm{W}$ is singular. This distribution is the basis of our spectral analysis
method.

###  B.2 PSD estimation

For $p=1$,
$\mathcal{CW}(\bm{\Sigma},M)=\mathcal{CW}(S_{x},M)=\Gamma(M,M/S_{x})$, with
$\Gamma(\alpha,\theta)$ the gamma distribution with shape parameter $\alpha$
and scale parameter $\theta$, and with $S_{x}$ the theoretical PSD of the sole
series $x(t)$.

This distribution can be used, along with a proper prior, to build the
Bayesian posterior for $S_{x}$. According to standard, physically sound
practice, we take an uninformative flat prior in the logarithm of $S$, which
is a prior proportional to $1/S$. For $p=1$ this also coincides with the
Jeffreys, re-parametrization independent prior [26].

Such choice gives a posterior for $S$ that is distributed as an inverse-Gamma
distribution, with shape parameter $M$ and scale parameter $M\Pi$:

$p(S|\Pi,M)\sim\text{inv}\Gamma\left(M,M\Pi\right)$ (22)

The expected value of this posterior, $\langle S\rangle=M\Pi/(M-1)$, is
slightly biased and diverges for $M\to 1$. However, the posterior probability
density function remains well-behaved even in the limit case; we have checked,
by a simple numerical simulation, that for all values of $M$ and for Gaussian
data, the theoretical value lies in the $p$-credible interval with probability
$p$.

Based on this posterior distribution, we compute the equally-tailed
$p$-credible interval from the $(1\pm p)/2$ quantiles of the distribution.
Unless otherwise specified, error bars in our plots correspond to $p=68.3\%$,
as for the $\pm 1\sigma$ interval in a normal distribution. The dots in those
same plots represent the median of the distribution, which remains well
behaved also for $M=1$, and has a smaller bias than the mean.

###  B.3  Decorrelation methods

Here we describe the method used to assess the contribution to $S_{\Delta g}$
of the disturbances for which we have a synchronously measured time series. As
said in the main text, we had to face two separate cases: one in which the
susceptibilities $\alpha_{i}$ in Eq. (14) are reasonably known, but the time
series of the disturbance carries some significant readout noise; in the
second case, on the contrary, the susceptibilities are unknown, except for
some rather loose order of magnitude estimate, but the readout noise can be
safely assumed to be negligible. These two cases are separately treated in the
following.

####  B.3.1  Decorrelation of noisy series with known susceptibilities

This is the case when we have a relatively narrow posterior distribution
$f(\alpha_{i})$ for the susceptibility $\alpha_{i}$, but $S_{n_{i}}$ may be
large, actually may even dominate $S_{z_{i},z_{i}}$ in Eq. (15), and refers to
the joint analysis of fluctuations of the average temperature and of LTP
distortion in Sec. VIII.2.1, and to that of the tank depletion gravitational
signal in Sec. VIII.2.2.

Though the two cases differ substantially, they both share the step of
estimating the Bayesian posterior for $\bm{\Sigma}$, the theoretical CPSD.

To that purpose, our starting point is again the distribution in Eq. (21) that
can be directly used to estimate the Bayesian posterior for $\bm{\Sigma}$, or
its inverse $\bm{Q}=\bm{\Sigma}^{-1}$, from the observation of $\bm{W}$, once
a prior distribution for either $\bm{\Sigma}$ or $\bm{Q}$ has been assumed.

Natural, non-informative choices are either a uniform prior on all elements of
$\bm{Q}$, or the Jeffreys prior [26], which, for a $p\times p$ dimensional
$\bm{Q}$ matrix, is $\propto\lvert\bm{Q}\rvert^{-p}$ [57, 58]. For the uniform
prior, the posterior is $\bm{Q}\sim\mathcal{CW}(\bm{W}^{-1},M+p)$, while for
the Jeffreys prior $\bm{Q}\sim\mathcal{CW}(\bm{W}^{-1},M)$

Both choices have limitations. Of those of the uniform prior when $p=1$ we
have discussed in App. B.2. On the other hand, the Jeffreys prior is affected
by bias, and by a significant inconsistency when $p>1$.

The bias consists of the fact that while the mean value of $\bm{Q}$ is
unbiased, the mean value of $\bm{\Sigma}$ is equal to $\bm{\Pi}M/(M-p)$. Thus
the bias depends on $p$ and may become large at low values of $M$.

The inconsistency stays in the fact that, when $p>1$, the marginal posterior
of any of the diagonal elements of $\bm{\Sigma}$, is not the same as the
posterior in Eq. (22), which one would calculate from the Jeffreys prior with
$p=1$.

The bias can be reduced and made independent of $p$, and the inconsistency
solved, if one takes the “intermediate” prior
$\propto\lvert\bm{Q}\rvert^{-1}$, that still coincides with the logarithmic
prior for the case $p=1$.

Note that, also with this choice, $\bm{Q}$ is Wishart distributed,
$\bm{Q}\sim\mathcal{CW}(\bm{W}^{-1},M+p-1)$, which guarantees, without
imposing any further prior constraint, that $\bm{Q}$ and its inverse are
positive definite matrices, a fundamental constraint for CPSD. We use this
posterior to estimate both the entire matrix, in the case of temperature/LTP
distortion, or its $S_{z_{i},\Delta g}(f)$ elements, in the case of the tank
depletion noise.

As we have a closed formula, we generate random samples both from such
posterior for $\bm{Q}$ and from the posterior for the $\alpha_{i}$, and
calculate then the samples of the quantities needed in the two cases, as
discussed in their relative sections.

####  B.3.2 Decorrelation of disturbances with unknown susceptibilities and
negligible readout noise

In the case of unknown susceptibilities and negligible readout noise, we
evaluate $S_{\Delta g_{0}}(f)$, and the susceptibilities assuming the
physically realistic model of real, and frequency-independent
susceptibilities.

The starting point is that one can transform the distribution in Eq. (21) to
separate its dependence on $S_{\Delta g_{0}}$, and the $\alpha_{i}$, from that
on the CPSD of the disturbances $S_{y_{i},y_{j}}$.

To do that consider that Eq. (14) defines a linear transformation $\Delta
g_{0}\to\Delta g,y_{i}\to y_{i}$ whose matrix is

$\bm{U}=\left(\begin{matrix}1&\alpha_{1}&\dots&\alpha_{r}\\\ 0&1&\dots&0\\\
\vdots&&\ddots&\vdots\\\ 0&\dots&\dots&1\end{matrix}\right)$ (23)

Note that $\lvert\bm{U}\rvert=1$, and that the inverse of $\bm{U}$ is obtained
with the transformation $\alpha_{i}\to-\alpha_{i}$.

As $\Delta g_{0}$ is independent of the $y_{i}$, before the transformation,
the processes have a CPSD matrix with block representation

$\bm{\Sigma}^{\prime}=\left(\begin{matrix}S_{\Delta g_{0}}&0\\\
0&\bm{\Sigma}_{yy}\end{matrix}\right)$ (24)

where $\bm{\Sigma}_{yy}$ is the $r\times r$ CPSD of the $y_{i}$. Note that

$\bm{\Sigma}^{\prime-1}=\left(\begin{matrix}1/S_{\Delta g_{0}}&0\\\
0&\bm{\Sigma}_{yy}^{-1}\end{matrix}\right),$ (25)

and that $\bm{\Sigma}=\bm{U}\bm{\Sigma}^{\prime}\bm{U}^{T}$, with $\bm{U}^{T}$
the transpose of $\bm{U}$, so that
$\lvert\bm{\Sigma}\rvert=\lvert\bm{\Sigma}^{\prime}\rvert=S_{\Delta
g_{0}}\lvert\bm{\Sigma}_{yy}\rvert$.

Furthermore, in Eq. (21) one can calculate that
$\mathrm{etr}\left[-\bm{\Sigma}^{-1}\bm{W}\right]=\mathrm{etr}\left[-(\bm{U}^{T})^{-1}\bm{\Sigma}^{\prime-1}\bm{U}^{-1}\bm{W}\right]=\mathrm{etr}\left[-\bm{\Sigma}^{\prime-1}{\bm{W}}^{\prime}\right]$,
having defined ${\bm{W}}^{\prime}=\bm{U}^{-1}\bm{W}(\bm{U}^{T})^{-1}$.

We need one more step to get our separation. $\bm{W}$ and ${\bm{W}}^{\prime}$
have block representation

$\bm{W}=\left(\begin{matrix}W_{1,1}&\bm{W}_{1,y}\\\
\bm{W}_{1,y}^{\dagger}&\bm{W}_{yy}\end{matrix}\right)\qquad{\bm{W}}^{\prime}=\left(\begin{matrix}{W}^{\prime}_{1,1}&{\bm{W}^{\prime}}_{1,y}\\\
{\bm{W}^{\prime}}_{1,y}^{\dagger}&\bm{W}_{yy}\end{matrix}\right)$ (26)

where one can calculate that

$\begin{split}&{W}^{\prime}_{1,1}=W_{1,1}-2\bm{\alpha}\cdot\text{Re}\bm{W}_{1,y}+\bm{\alpha}\cdot\text{Re}\bm{W}_{yy}\cdot\bm{\alpha},\end{split}$
(27)

with $\alpha$ the r-vector the components of which are the susceptibilities.

The consequence of Eqs. 25, 26, and 27 is that

$\begin{split}&\mathrm{etr}\left[-\bm{\Sigma}^{-1}\bm{W}\right]=\mathrm{etr}\left[-\bm{\Sigma}^{\prime-1}\tilde{\bm{W}}\right]=\\\
&=e^{-\frac{{W}^{\prime}_{1,1}}{S_{\Delta
g_{0}}}-\text{tr}\left[\bm{\Sigma}_{yy}^{-1}\bm{W}_{yy}\right]}=\\\
&=e^{-\frac{W_{1,1}-2\bm{\alpha}\cdot\text{Re}\bm{W}_{1,y}+\bm{\alpha}\cdot\text{Re}\bm{W}_{yy}\cdot\bm{\alpha}}{S_{\Delta
g_{0}}}-\text{tr}\left[\bm{\Sigma}_{yy}^{-1}\bm{W}_{yy}\right]}.\end{split}$
(28)

Putting all together, the distribution of $\bm{W}$, conditional on the
susceptibilities, the PSD of the residual noise, and the CPSD of the
disturbances, becomes then:

$\displaystyle p\left(\bm{W}\big{|}\alpha_{i},S_{\Delta
g_{0}},\bm{\Sigma}_{yy}\right)=\frac{\left|\bm{W}\right|^{M-p}}{\widetilde{\Gamma}_{p}(M)}\times$
(29) $\displaystyle\times\frac{1}{S_{\Delta
g_{0}}^{M}}e^{-\frac{W_{1,1}-2\bm{\alpha}\cdot\text{Re}\bm{W}_{1,y}+\bm{\alpha}\cdot\text{Re}\bm{W}_{yy}\cdot\bm{\alpha}}{S_{\Delta
g_{0}}}}\times$
$\displaystyle\times\frac{1}{\left|\bm{\Sigma}_{yy}\right|^{M}}\mathrm{etr}\left[-\bm{\Sigma}_{yy}^{-1}\bm{W}_{yy}\right]$

The distribution in Eq. (29) can be used to build a joint posterior for
$S_{\Delta g_{0}}^{-1}$ and the $\alpha$’s, and an independent posterior for
$\bm{\Sigma}_{yy}$.

Treating the data at the different frequencies as independent, and assuming
that $S_{\Delta g_{0}}$ and $\bm{\Sigma}_{yy}$ depend on the frequency, while
the $\alpha$’s do not, the joint likelihood for the Bayesian inference of
$S_{\Delta g_{0}}$ and $\alpha_{i}$ becomes:

$\begin{split}&p(\alpha_{i},S_{\Delta g_{0}}(f_{1}),S_{\Delta
g_{0}}(f_{2}),\dots|\bm{W})\propto\\\
&\propto\prod_{k}\frac{\widetilde{p}(S_{\Delta g_{0}}(f_{k}))}{S_{\Delta
g_{0}}(f_{k})^{M(f_{k})}}e^{-\sum_{k}\frac{W_{1,1}(f_{k})}{S_{\Delta
g_{0}}(f_{k})}}\times\\\ &\times
e^{2\bm{\alpha}\cdot\sum_{k}\frac{\text{Re}\bm{W}_{1,y}(f_{k})}{S_{\Delta
g_{0}}(f_{k})}-\bm{\alpha}\cdot\left(\sum_{k}\frac{\text{Re}\bm{W}_{yy}(f_{k})}{S_{\Delta
g_{0}}(f_{k})}\right)\cdot\bm{\alpha}},\end{split}$ (30)

having taken a flat prior for the $\alpha$’s, and the prior
$\widetilde{p}(S_{\Delta g_{0}}(f_{k}))$ for $S_{\Delta g_{0}}(f_{k})$. As for
this one, to take into account the presence of the Brownian noise, we split it
as $S_{\Delta g_{0}}(f_{k})\to S_{\Delta g_{0}}(f_{k})+S_{\text{Brown}}$, and
we take a uniform prior for the logarithms of $S_{\Delta g_{0}}(f_{k})$. For
$S_{\text{Brown}}$ we use as prior the posterior obtained for its Bayesian
estimate (see Sec. III). At frequencies low enough that the Brownian noise is
negligible, which is where the excess is best measured, this is fully
equivalent to taking the Jeffreys prior for $S_{\Delta g_{0}}(f_{k})^{-1}$,
which is indeed a uniform prior on its logarithm.

To estimate the posterior of the parameters, we employ a parallel-tempering
Monte Carlo Markov chain algorithm (PT-MCMC [59]). We always find smooth
distributions, no bi-modalities and no strong cross-correlation among
parameters. To implement the prior, the MCMC algorithm explores the parameter
space of the susceptibilities, and of the logarithms of $S_{\Delta
g_{0}}(f_{k})$ and of $S_{\text{Brown}}$.

As for $\bm{\Sigma}_{yy}$, Eq. (29) shows that, using the intermediate prior,
$\bm{\Sigma}_{yy}^{-1}(f_{i})\sim\mathcal{CW}(\bm{W}_{yy}^{-1}(f_{i}),M(f_{i})+r-1)$.
Thus the evaluation of the statistical properties of the posterior can be
obtained numerically from the relative Wishart distribution.

Once the posteriors for $S_{\Delta g_{0}}(f_{k})$, $S_{\text{Brown}}$,
$\alpha_{i}$, and possibly for $\bm{\Sigma}_{yy}$, have been obtained, one can
calculate the posterior for
$S_{c}=\sum_{i,j}\alpha_{i}\alpha_{j}S_{y_{i},y_{j}}=\sum_{i,j}\alpha_{i}\alpha_{j}\left(\Sigma_{yy}\right)_{i,j}$.
However, for small values of $S_{c}/S_{\Delta g_{0}}$, the statistics of such
a posterior is biased toward high values, with the bias increasing with $r$.
This is due to the positive definite quadratic nature of $S_{c}$, which gets
rapidly dominated by susceptibility fluctuations.

To deal with this problem we have resorted to calibrating the bias via
simulations. Our simulation consists of the following steps:

  * •

We extract a sample of $\bm{\Sigma}_{yy}(f)$ from its calculated posterior,
and form a corresponding sample for $\bm{\Sigma}(f)$ assuming $S_{\Delta
g_{0}}=0$.

  * •

We generate a random vector of susceptibilities $\alpha_{i}$ from a normal
distribution with zero mean, and standard deviation adjusted such that, at
some preferred frequency, $S_{c}\simeq\rho S_{\Delta g_{e}}$, with $\rho$ a
desired value of the noise power fraction.

  * •

From the susceptibility vector we form the matrix $\bm{U}$ of Eq. (23), and
with it we form the sample of the (provisional) matrix
$\bm{\Sigma}(f)=\bm{U}\cdot\bm{\Sigma}^{\prime}(f)\cdot\bm{U}^{T}$. This
matrix has $\Sigma_{1,1}=S_{\Delta g}(f)=S_{c}(f)$, a number we store as we
need for our calibration procedure.

  * •

We then form the final value of the matrix $\bm{\Sigma}$ by substituting the
value of $\Sigma_{1,1}$ above, with a random sample extracted from the
posterior distribution of $S_{\Delta g}(f)$ (see Eq. (22)). This our
theoretical CPSD matrix, which has a value of $S_{c}(f)$, coherent with its
elements $\Sigma_{1,j}(f)$, a value for $S_{\Delta g}=S_{\Delta
g_{0}}(f)+S_{c}$ coherent with its posterior, and a value of
$\bm{\Sigma}_{yy}(f)$ also coherent with its posterior.

  * •

From the matrix $\bm{\Sigma}(f)$ above we generate a random sample of the
simulated $\bm{W}$ matrix, using the Wishart distribution in Eq. (21). On each
of these simulated samples of $\bm{W}$ we perform our entire Bayesian
procedure thus getting, among others, a posterior for $S_{c}(f)$

We repeat the procedure 100 times for different values of the power fraction
$\rho$ and we get then, at each value of $\rho$ and at each of the considered
frequencies, a distribution of the ‘true’ values of $S_{c,t}(f)$ and a global
distribution of the estimated values $S_{c,e}(f)$ for the same quantity (see
Fig. 24). We use mean and standard deviations of these distributions, to do a
linear, weighted least square fit $S_{c,t}(f)=AS_{e,t}(f)+B$, that we use to
correct the observed data. The procedure gives a result close to what one
would obtain by just subtracting, from the estimated value of $S_{c}(f)$, the
median from the simulation for $\rho=0$. Indeed in some cases, to avoid
repeating time-consuming simulations at different values of $\rho$, we adopted
this second method.

FIG. 24: $\pm 1\sigma$ intervals for the distributions of the ‘true’ (Tr.)
values of $S_{c}(f)$ and for their global estimated posteriors (Est.) over 100
simulations for each case. The power fraction $\rho$ is an approximation of
the average at the three highest frequencies. In particular $\rho=0$ is in
reality calculated at $\rho\simeq 4\times 10^{-6}$. Data are plotted for
clarity at slightly shifted frequency coordinates, but have all been
calculated at one of the frequencies $f_{i}$ with $4\leq i\leq 8$.

With simulations, we have also checked that the method is unbiased for the
estimate of both $S_{\Delta g_{0}}$ and the susceptibilities.

As a final note, the case of the evaluation of the effect of the temperature
on the first bin involves just one frequency and one disturbance. In such case
the likelihood in Eq. (30) can be integrated analytically, to give a marginal
probability density for $S_{\Delta
g,0}\sim\text{inv}\Gamma\left(M-1,M/(\bm{\Pi}^{-1})_{1,1}\right)$, for the
general case of a complex $\alpha$, and $S_{\Delta
g,0}\sim\text{inv}\Gamma\left(M-1/2,M/(\text{Re}\bm{\Pi}^{-1})_{1,1}\right)$
if alpha is real. We have used these formulas to evaluate the effect of the
temperature.

##  Appendix C Choice of quasi-independent frequencies

As presented in [1, Supp.mat.], we evaluate the CPSD at log-spaced
frequencies, such that the correlation among adjacent frequencies is kept
below 5%. The number of averaged periodograms $M$ varies with frequency, so
that the variance is optimally reduced as more periodograms are available.
However, $M$ is the same for the first and second frequencies: even though
this induces a slightly higher correlation, it allows to analyze the lowest
frequency $18\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$. Frequencies are chosen
so that the lowest bins always have a fair amount of periodograms, and the
fourth one ($0.1\text{\,}\mathrm{mHz}$) is the lower bound of the official
LISA frequency band. The frequencies are listed in Table 6.

TABLE 6: List of frequencies used for spectral estimation.

##  Appendix D ASD of $\Delta g(t)$ for all runs

The following figures, 25 and 26, show the ASD for all runs of Table 1, except
for run 10 which is shown in Fig. 2.

FIG. 25: PSD of $\Delta g(t)$ for runs from 1 to 6. The meaning of all symbols
is the same as in Fig. 2. FIG. 26: PSD of $\Delta g(t)$ for runs from 7 to 13,
except for run 10 which has already been shown in Fig. 2. The meaning of all
symbols is the same as in Fig. 2.

##  Appendix E  Calculation of gravitational long-term variations

###  E.1 Instrument distortion

As already said in Sec. V.6, the capacitive motion sensor gives a measurement
of the relative motion of the GRS relative to its own test mass,
$x_{\text{GRS},i}=x_{i}-X_{i}$. Here the suffix $i$ indicates the test mass,
and $X_{i}$ is the coordinate of the GRS surrounding TM$i$.

As we are only interested in time variations of differences of coordinates,
here the reference frame may be taken that of the spacecraft, and each
coordinate is zero at the nominal position of the body it belongs to, that is,
at the position of the body in the absence of distortion.

The variation of gravitational field due to such distortion-induced motion of
both GRS and the TMs can be written as:

$\begin{split}&\Delta g_{\text{Dist},0}=\\\
&-\omega_{2,2}^{2}\left(x_{2}-X_{2}\right)-\omega_{\text{SC},2}^{2}x_{2}-\omega_{2,1}^{2}\left(x_{2}-X_{1}\right)+\\\
&+\omega_{1,1}^{2}\left(x_{1}-X_{1}\right)+\omega_{\text{SC},1}^{2}x_{1}+\omega_{1,2}^{2}\left(x_{1}-X_{2}\right)\end{split}$
(31)

Here $-\omega_{i,j}^{2}$ is the force gradient (per unit mass) caused by GRS
$j$ on TM $i$ and $-\omega_{\text{SC},i}^{2}$ is that caused by all the
remaining parts of the spacecraft, mostly the OMS, on TM $i$. Notice that
$\omega_{i}^{2}$ in Eq. (1) is given by
$\omega_{i}^{2}=\left(\omega_{i,i}^{2}+\omega_{\text{SC},i}^{2}+\omega_{i,j}^{2}\right)$.

Within this model the $\Delta g$ time series is expected to be affected by a
distortion-induced component:

$\begin{split}&\Delta g_{\text{Dist}}=\Delta
g_{\text{Dist},0}+\omega_{2}^{2}\left(x_{2}-x_{1}\right)+\left(\omega_{2}^{2}-\omega_{1}^{2}\right)x_{1}=\\\
=&\left(\omega_{2,2}^{2}-\omega_{2,1}^{2}\right)X_{2}-\left(\omega_{1,1}^{2}-\omega_{1,2}^{2}\right)X_{1}\simeq\\\
&+\simeq\left(\omega_{2,2}^{2}-\omega_{2,1}^{2}\right)\left(X_{2}-X_{1}\right)\end{split}$
(32)

The last simplification comes from the fact that, due the symmetry of the
instrument,
$\left(\omega_{2,2}^{2}-\omega_{2,1}^{2}\right)\simeq\left(\omega_{1,1}^{2}-\omega_{1,2}^{2}\right)$
within the same accuracy within which $\omega_{1}^{2}\simeq\omega_{2}^{2}$, as
discussed in Sec. II.2.

The relative motion of the two GRS, $\Delta x=X_{2}-X_{1}$, can be obtained as

$\begin{split}&\Delta X=X_{2}-X_{1}=\Delta x_{\text{OMS}}-\Delta
x_{\text{GRS}}=\\\
&=x_{2}-x_{1}-\left(\left(x_{2}-X_{2}\right)-\left(x_{1}-X_{1}\right)\right)\end{split}$
(33)

Thus, in conclusion,

$\begin{split}&\Delta g_{\text{Dist}}=\\\
&=\left(\omega_{2,2}^{2}-\omega_{2,1}^{2}\right)\left(\Delta
x_{\text{OMS}}-\Delta x_{\text{GRS}}\right)\equiv\\\
&\equiv\omega_{d}^{2}\left(\Delta x_{\text{OMS}}-\Delta
x_{\text{GRS}}\right)\end{split}$ (34)

$\omega_{2,2}^{2}$ is dominated by both gravitational and electrical effects
[18], while the origin of $\omega_{2,1}^{2}$ is just gravitational. We can get
an estimate of the electrical terms from [18], while the gravitational
components have been estimated in [60], within the work to suppress the
gravitational field at the TM location [11].

Based on these references, we calculate $\omega_{d}^{2}=(-3.32\pm 0.05)\times
10^{-7}\leavevmode\nobreak\ \text{s}^{-2}$.

###  E.2 Propellant tank depletion

The sign and the magnitude of $\Delta g_{\text{Tanks}}$ depends on the
specific tanks from which the propellant is taken. Tank 1 and 2 of Table 2 in
App. A, produce, respectively, a $\Delta g_{\text{Tanks}}$ of
$\lambda_{t,1}\simeq+39\text{\,pm s}^{-2}\text{ kg}^{-1}$ and
$\lambda_{t,2}\simeq+37\text{\,pm s}^{-2}\text{\,kg}^{-1}$ per unit mass of
contained propellant, while for tank 3 this figure is instead
$\lambda_{t,3}\simeq-43\text{\,pm\,s}^{-2}\text{\,kg}^{-1}$. These figures are
the result of a numerical calculation on the geometry and location of the
tanks, and are affected by a relative error of $\simeq 5\%$.

The mass loss is monitored by set of flow-meters, one for each of the
thrusters, that measure the instantaneous flow of mass through that thruster.
As the thrusting system contains two branches (see App. A), the measurement of
the total flow of mass from the tank in use during any given run, $\dot{m}$,
is the sum of the readings of the six flow-meters of the thrusters belonging
to the branch in use during that specific run. We estimate an absolute
accuracy on $\dot{m}$, for both branches, of the order of 10%.

From the figures above, one can estimate $\Delta g_{\text{Tank}}$, during a
specific run that uses tank $i$ and branch $j$, as

$\begin{split}&\Delta
g_{\text{Tank}}(t)=\kappa_{t,i}\;\lambda_{t,i}\;\kappa_{b,j}\left(m_{j}(0)-\int_{0}^{t}\dot{m}_{j}(t)dt\right)\equiv\\\
&\equiv\kappa_{t,i}\kappa_{b,j}\Delta g_{\text{Tank},0,i,j}(t)+c\end{split}$
(35)

Here we have introduced:

  * •

a calibration factor $\kappa_{t,i}$ for the gravitational gradient of tank
$i$.

  * •

a calibration factor $\kappa_{b,j}$ for the reading $\dot{m}_{\text{j}}(t)$ of
the mass outflow from tank $j$,

  * •

the initial propellant mass $m_{j}(0)$,

  * •

the nominal signal $\Delta
g_{\text{Tank},0,i,j}(t)=-\lambda_{t,i}\int_{0}^{t}\dot{m}_{j}(t)dt$,

  * •

the constant $c=\lambda_{t,i}\kappa_{t,i}\kappa_{b,j}m_{j}(0)$ that has no
relevance for the discussion.

In Sec. V.6 we discuss both the nominal case $\kappa_{t,i}=\kappa_{b,j}=1$ and
some other options.

###  E.3 Drift evaluation procedure

On the long time scale of days or more we are considering here, all data
series have strong auto-correlation. Thus the linear square fitting to data,
that we do to estimate the residual drift, does not allow a consistent and
unbiased estimate of errors. To partially circumvent this limitation we have
resorted to the following procedure.

  * •

We chose the values for $\kappa_{t,i}$, $\kappa_{b,j}$ and $\omega_{d}^{2}$.
One obvious choice is $\kappa_{t,i}=\kappa_{b,j}=1$ for all values of i and j,
and $\omega_{d}^{2}=-3.32\times 10^{-7}\text{\,s}^{-2}$ that we call the
nominal calibration.

  * •

For each run, we partition all relevant data series into one day-long non-
overlapping stretches. Let us call $N_{k}$ the resulting number of stretches
for run $k$

  * •

For each of the stretches in one run, we form the ‘corrected’ data series
$\Delta g_{c}(t)=\Delta g(t)-\kappa_{t,i}\kappa_{b,j}\Delta
g_{\text{Tank},0,i,j}(t)-\omega_{d}^{2}\left(\Delta x_{\text{OMS}}-\Delta
x_{\text{GRS}}\right)$, using the proper values for tank $i$ and branch $j$
used in that same run.

  * •

Again for each of the stretches in one run, we perform a least square fitting
of $\Delta g_{c}(t)$ with the model

$\Delta g_{c}(t)=c_{T}T+c_{t}t+\text{c}$ (36)

  * •

From the $N_{k}$ long sample of fitting coefficient pairs $c_{T,j},\;c_{t,j}$,
with $1\leq j\leq N_{k}$, obtained by fitting the $N_{k}$ stretches of run
$k$, we form the average partial derivatives $\left(\partial\Delta g/\partial
T\right)_{k}\equiv(1/N)\sum_{k=1}^{N}c_{T,k}$ and $\left(\partial\Delta
g/\partial t\right)_{k}\equiv(1/N)\sum_{k=1}^{N}c_{t,k}$, and the
corresponding variances $\sigma_{T,k}^{2}$ and $\sigma_{t,k}^{2}$.

  * •

As $N_{k}$ is quite small for many of the runs, to get a more solid estimate
of the mean fluctuation of a single determination of coefficients, we make a
weighted average over the variance on the entire set of run:

$\begin{split}&\sigma_{T}=\sqrt{\frac{\sum_{k=1}^{13}\left(N_{k}-1\right)\sigma_{T,k}^{2}}{\sum_{k=1}^{13}\left(N_{k}-1\right)}}\\\
&\sigma_{t}=\sqrt{\frac{\sum_{k=1}^{13}\left(N_{k}-1\right)\sigma_{t,k}^{2}}{\sum_{k=1}^{13}\left(N_{k}-1\right)}}\end{split}$

  * •

Finally, we take $\sigma_{T}/\sqrt{N_{k}}$ and $\sigma_{t}/\sqrt{N_{k}}$ as
the uncertainty on $\left(\partial\Delta g/\partial T\right)_{k}$ and
$\left(\partial\Delta g/\partial t\right)_{k}$ respectively

It is worth stressing that such a procedure remains largely empirical, and
that coefficients from stretches within the same run may still be correlated
so that the errors may be still underestimated. However the main results of
the discussion in Sec. V.6 should depend only weakly on the details of the
procedure.

##  Appendix F  In-band noise from instrument distortion and thermal effects

In Sec. VIII.2.1, we discussed the in-band contribution to $\Delta g$ of the
instrument distortion and thermal effects. Here we give details and
calculations.

_Susceptibilities_ — A set of thermal experiments [61] has shown the
possibility of delay effects between $T$ and the other time series, reasonably
represented by a simple pole filter with a cut-off frequency in the
$\mathrm{\SIUnitSymbolMicro Hz}$ range. We summarise all this information, by
taking the following distributions for the susceptibilities:

  * •

$\omega_{d}^{2}$ is frequency-independent and normally distributed with
$\omega_{d}^{2}=(-0.33\pm 0.03)\times 10^{-6}\,$\mathrm{s}^{-2}$$.

  * •

$\alpha_{T}(f)$, the susceptibility of $\Delta g$ to $T$, is given by
$\alpha_{T}(f)=\alpha_{T,0}\left(1+if/f_{0}\right)^{-1}$, with $\alpha_{T,0}$
normally distributed with $\alpha_{T,0}=(0.4\pm
0.2)\,$\mathrm{p}\mathrm{m}\,\mathrm{s}^{-2}\,\mathrm{K}^{-1}$$, and $\log
f_{0}$ uniformly distributed between of
$\log($1\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$)$ and
$\log($100\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$)$.

  * •

$\gamma(f)$, the susceptibility of $\Delta X$ to $T$, is given by
$\gamma(f)=\gamma_{0}\left(1+if/f_{1}\right)^{-1}$, with $\gamma_{0}$ normally
distributed with $\gamma_{0}=(1.04\pm
0.07)\,$\mathrm{\SIUnitSymbolMicro}\mathrm{m}\,\mathrm{K}^{-1}$$, and $f_{1}$
with the same distribution as for $f_{0}$.

_First approach_ — The first obvious, naive way to put an upper limit to the
combined contribution of $\Delta g_{\text{Dist}}$ and $\Delta g_{T}$ to the
PSD of $\Delta g$, is to: assume that these are not affected by any readout
noise; take the posterior for the joint CPSD of $\Delta X$, $\Delta T$ and
$\Delta g$ from the proper distribution, as explained in App. B.3.1; take the
posterior for the susceptibilities as explained above. From all this, the
total contribution to $S_{\Delta g}$ of $\Delta X$ and $\Delta T$ would be

$\omega_{d}^{4}S_{\Delta
X}(f)+\left(2\text{Re}\left[\omega_{d}^{2}\alpha_{T}(f)\gamma^{*}(f)\right]+\lvert\alpha_{T}(f)\rvert^{2}\right)S_{T}(f),$
(37)

_Second approach_ — A more informative bound can be obtained explicitly
including the presence of the readout noise affecting the time series. We
include the following considerations:

  1. 1.

The measured $\Delta X$ and $T$ are affected by unmodeled and uncorrelated
readout noises. Moreover, the noise affecting $\Delta X$ correlates to $\Delta
g_{0}$ in Eq. (15).

  2. 2.

$\Delta\Phi$ is also affected by readout noise, which correlates with $\Delta
X$.

  3. 3.

The angular distortion $\Delta\Phi$ has a negligible contribution to the
linear acceleration noise $\Delta g$.

We define the matrices $\bm{\Sigma}_{\Delta X}$, $\bm{\Sigma}_{T}$, and
$\bm{\Sigma}_{0}$ so that the CPSD matrix of $(\Delta g,\Delta
X,T,\Delta\Phi)$ can be written as

$\bm{\Sigma}=\bm{\Sigma}_{0}+S_{\Delta X}\bm{\Sigma}_{\Delta
X}+S_{T}\bm{\Sigma}_{T}+S_{\Delta\Phi}\bm{\Sigma}_{\Delta\Phi}$ (38)

Here, $\bm{\Sigma}_{0}$ represents the CPSD between $\Delta g_{0}$ and the
readout noise part of $\Delta X$, $T$, and $\Delta\Phi$. According to the
previous assumptions, the matrices read

$\begin{split}&\bm{\Sigma}_{T}=\left(\begin{matrix}|\alpha_{T}|^{2}+2\text{Re}\left(\alpha_{T}\omega_{d}^{2}\gamma^{*}\right)&\alpha_{T}\gamma^{*}&\omega_{d}^{2}\gamma+\alpha_{T}&0\\\
\alpha_{T}^{*}\gamma&0&\gamma&0\\\
\omega_{d}^{2}\gamma^{*}+\alpha_{T}^{*}&\gamma^{*}&1&0\\\
0&0&0&0\end{matrix}\right)\\\ &\bm{\Sigma}_{\Delta
X}=\left(\begin{matrix}\omega_{d}^{4}&\omega_{d}^{2}&0&0\\\
\omega_{d}^{2}&1&0&0\\\ 0&0&0&0\\\ 0&0&0&0\end{matrix}\right)\end{split}$ (39)

Within this model, we estimate the contribution to $S_{\Delta g}$. First, we
take a sample from the posterior of $\bm{\Sigma}$ (see App. B.3.1) and one
from that of the susceptibilities, from which we build $\bm{\Sigma}_{0}$ using
Eqs. (38) and (39). Then, we search the values of $S_{\Delta X}$, $S_{T}$, and
$S_{\Delta\Phi}$ that give the maximum noise contribution $\left(S_{\Delta
X}\bm{\Sigma}_{\Delta X}+S_{T}\bm{\Sigma}_{T}\right)_{1,1}$, subject to the
condition that $\bm{\Sigma}_{0}$ is a positive definite matrix, as a CPSD
should always be. We repeat the calculation over a thousand different
posterior samples, deriving a posterior also for this upper limit. The result
is in Fig. 16 on page 16.

##  Appendix G Noise from Poisson flow of glitches

A very simple model for a stochastic process that shows a $1/f$ ASD, is a flow
of steps $\Delta g(t)=g_{i}\Theta(t-t_{i})$, with $\Theta(t)$ the Heaviside
theta function, $g_{i}$ a random amplitude, and with the arrival times $t_{i}$
following Poisson event statistics. Such a process has indeed an ASD given by
$S_{\Delta g}^{1/2}=\sqrt{2\lambda\langle g_{i}^{2}\rangle/(2\pi f)^{2}}$,
with $\lambda$ the event rate. 555Strictly speaking the calculation must be
done for steps $\Delta g(t)=g_{i}e^{-(t-t_{i})/\tau}\Theta(t-t_{i})$, and the
result is valid for $f\ll 1/(2\pi\tau)$

The process can be generalized. Specifically, one can calculate the ASD of a
flow of Poisson events $\Delta g(t)=g_{i}h_{i}(t-t_{i})\Theta(t-t_{i})$, with
$h_{i}(t)$ the unit peak amplitude version of a randomly selected glitch among
those observed during ordinary LPF noise runs [20]. As we only observed
$N_{g}=98$ glitches, to generate an arbitrary long time series, in this
process glitch shapes need to be repeated.

The PSD of such a process would be

$S_{\Delta g,gl}(f)^{1/2}=\sqrt{2\lambda\langle
g_{i}^{2}\rangle\frac{1}{N_{g}}\sum_{l=1}^{N_{g}}\lvert h_{l}(f)\rvert^{2}},$
(40)

Here $h_{l}(t)$ is the shape of the $l$-th observed glitch, and $h_{l}(f)$ its
Fourier transform. For the sake of the current discussion all 98 glitches had
$h_{l}(t)\propto e^{-t/\tau_{2,l}}-e^{-t/\tau_{1,l}}$ with $\tau_{2,l}$ and
$\tau_{1,l}$ two time constants666Note that here we normalize $h_{l}(t)$ to
have unit peak force-per-unit-mass amplitude, at variance with the
normalization adopted in Ref. [20], where $h_{l}(t)$ had unit impulse per-
unit-mass. [20].

Figure 27 compares the ASD calculated from Eq. (40) taking $2\lambda\langle
g_{i}^{2}\rangle=9.8\times 10^{-5}\text{\,fm}^{2}\text{s}^{-5}$ (the dashed
black line), to data from run 10 (the black dots). To ease the comparison, the
figure reports $\sqrt{S_{\Delta g,gl}+S_{\text{Brown}}}$, with
$S_{\text{Brown}}$ the maximum likelihood value for run 10.

FIG. 27: ASD from a Poisson flow of glitches. Dashed black line: calculated
ASD of a flow of glitches taken at random from the sample of the observed ones
and with $2\lambda\langle g_{i}^{2}\rangle=9.8\times
10^{-5}\text{\,fm}^{2}\text{s}^{-5}$ (see text for definition of symbols).
Blue dot-dashed line: calculated ASD of a flow of glitches with time
parameters taken at random from the extended distribution described in the
text and with $2\lambda\langle g_{i}^{2}\rangle=8.0\times
10^{-5}\,\text{fm}^{2}\text{s}^{-5}$. Red dots: ASD of a simulated time series
with the same length as that of run 10, with glitch parameters extracted from
the extended simulation, $g_{i}$ extracted from a uniform distribution in the
interval $\pm 0.9\,\text{fm\,s}^{-2}$, and with a rate of
$\lambda=105\,\text{d}^{-1}$.

It must be said the indistinguishable results, except for the numerical value
of $2\lambda\langle g_{i}^{2}\rangle$, are obtained by taking just the 54
glitches for which $\tau_{2,l},\tau_{1,l}\geq 10\,\text{s}$.

The figure shows that such a simple model may well be fit to the observed
ASDs, at least for frequencies $f\gtrsim$30\text{\,}\mathrm{\SIUnitSymbolMicro
Hz}$$. However below that limit the predicted ASD saturate, as none of the
observed time constants ever exceeded $\simeq 8\,\text{ks}$ [20]. To achieve a
better result one needs then to slightly extend the family of templates
$h(t)$.

To do that, we first note that more than 90% of the glitches for which
$\tau_{2,l},\tau_{1,l}\geq 10\,\text{s}$, had
$\tau_{2,l}=\tau_{1,l}\equiv\tau_{l}$, a case in which $h_{l}(t)\to
e(t/\tau_{l})e^{-t/\tau_{l}}$. For these glitches the distribution of
$\log_{10}(\tau_{l}/1\,\text{s})$ is quantitatively compatible with a uniform
one in the range 1 to 4.

The simplest model extension is then to a flow of glitches of the kind $\Delta
g(t)=g_{i}e(t/\tau_{i})e^{-t/\tau_{i}}\Theta(t-t_{i})$ with the distribution
of $\log_{10}(\tau_{i}/1\,\text{s})$, extended up to 5, instead of being
limited to 4.

The results for such an extended model, with a fitting constant
$2\lambda\langle g_{i}^{2}\rangle=8.0\times
10^{-5}\,\text{fm}^{2}\text{s}^{-5}$ are reported again in Fig. 27 (blue dot-
dashed line).

We now discuss the possibility that the glitches that form this hypothetical
random flow, the _noise glitches_ , may be the same that occasionally become
large enough to be detected as isolated signals in the data series, the
_isolated glitches_. To do that let us first discuss the possible parameter
distribution of the noise glitches.

The model above would reproduce our observations if the values for
$2\lambda\langle g_{i}^{2}\rangle$ would only refer to glitches with amplitude
$g_{i}$ small enough to be undetectable against the background noise. Indeed
the data series on which we have calculated noise ASD have been purged of any
detectable glitch [2], and thus of their contribution to the ASD. This puts a
constraint on the distribution of $g_{i}$ for the noise glitches, that depends
on our ability to detect a glitch against the background noise.

Our empirical glitch detection method has been described in [20], and its
detection ability found in substantial agreement with the prediction of a
search based on a matched filter. We have repeated such an analysis for the
amplitude normalisation we use here, calculating the joint Fisher matrix for
$g_{i}$, $\tau_{1,i}$, $\tau_{2,i}$ and the glitch arrival time, and confirmed
such an agreement.

Specifically all detected glitches have signal to noise ratio $\text{SNR}\geq
3.6$, with $\text{SNR}=\lvert g_{i}\rvert/\sigma_{g}$, and $\sigma_{g}$ the
uncertainty on $g_{i}$ predicted by the Fisher matrix, and all but 4 have a
signal to noise ratio $\text{SNR}\geq 5$. Indeed visual inspection confirms
that our empirical method would detect almost certainly any glitch with
$\text{SNR}=5$, while would almost certainly not detect a glitch with
$\text{SNR}\leq 3.5$.

For a given amplitude $g_{i}$, SNR depends on $S_{\text{Brown}}$, but also on
both $\tau_{1,i}$ and $\tau_{2,i}$. For run 10, the run with best noise, the
SNR reaches a maximum for
$\tau_{1,i}=\tau_{2,i}\simeq$85\text{\,}\mathrm{s}$$.

Putting all these elements together, and reminding the stability of the
observed noise in the course of the mission, the distribution of $g_{i}$ and
$\tau_{i}$ for the noise glitches should fulfill the following conditions:

  * •

$g_{i}$ and $\tau_{i}$ should be independent;

  * •

their distributions should both be independent of the run;

  * •

$\log_{10}(\tau_{i}/1\,\text{s})$ should be uniformly distributed between 1
and 5

  * •

the distribution of $g_{i}$ should fulfill: $2\lambda\langle
g_{i}^{2}\rangle\simeq 8.0\times 10^{-5}\text{\,fm}^{2}\text{s}^{-5}$

  * •

the distribution of $g_{i}$ should assign a probability $p_{t}$ to
$g_{i}+s_{i}\geq 5\sigma_{g}$, with $s_{i}$ the random amplitude measurement
error, such that $\lambda p_{t}$ is much less than the observed rate, of
$\simeq 1\text{\,d}^{-1}$ of the isolated glitches. Here all quantities refer
to run 10 and to $\tau_{1,i}=\tau_{2,i}\simeq 85\text{s}$ which is the most
favorable case for glitch detection.

Thus the question if the isolated glitches are just a sample of large
amplitude tail of an overall distribution of glitch parameters, translates in
that if the amplitude distribution of the noise glitches, fulfilling the
constraints above, merges at high amplitude into that of the isolated
glitches, without the need of too many pathological assumptions to fill up the
gaps among the two.

In order to start with the minimum such gap, in the condition for
$2\lambda\langle g_{i}^{2}\rangle$ one should pick the lowest possible value
for $\lambda$ corresponding to the largest for $\langle g_{i}^{2}\rangle$.
Barring the non-physical assumption of all equal $g_{i}$’s, the simple
distribution that maximizes $\langle g_{i}^{2}\rangle$ for a given upper bound
$g_{\text{max}}$, is a uniform distribution $-g_{\text{max}}\leq g_{i}\leq
g_{\text{max}}$ for which $\langle g_{i}^{2}\rangle=(1/3)g_{\text{max}}^{2}$.

To decide the value of $g_{\text{max}}$ we take the conditions of run 10, with
$\tau_{1,i}=\tau_{2,i}\simeq$85\text{\,}\mathrm{s}$$, recast
$g_{\text{max}}=n_{\text{max}}\sigma_{g}$, and calculate the probability that
a sample from a uniform distribution $-n_{\text{max}}\sigma_{g}\leq g_{i}\leq
n_{\text{max}}\sigma_{g}$, plus an independent sample from a zero-mean normal
distribution with variance $s_{g}^{2}$, representing the measurement error,
exceeds in absolute value, the threshold $5\sigma_{g}$. We pick
$n_{\text{max}}$ when this probability is low enough to detect at most one
detectable glitch in a run with the same duration of run 10. The procedure
requires a few trial and error loops involving adjusting the value of
$\lambda$ that try to minimize.

Figure 27 shows the ASD (red dots) of simulated data assuming
$\lambda=102\text{\,d}^{-1}$, $n_{\text{max}}=3$,
$\sigma_{g}=0.1\text{\,fm\,s}^{-2}$, and the Brownian noise level and the time
duration of run 10. Inspection of the data shows no glitch we would have
detected with our empirical method.

We have simulated noise glitches to reproduce the noise for all runs in which
we have observed the glitches [20], by generating amplitudes and time
constants from the same distributions, independently of the run. For each of
the simulated glitches, we have also generated the amplitude that would have
been measured by the matched filter, by adding a simulated measurement error,
and tested if the glitch would have been detected within the specific run to
which it belonged. The results of such collective simulation are shown in Fig.
28.

FIG. 28: Measured amplitude $g_{i,m}$ versus time constant $\tau_{i}$ for both
simulated noise glitches and observed isolated ones. Black dots: all observed
isolated glitches with $\tau_{2,i}=\tau_{1,i}\equiv\tau_{i}$ and
$\tau_{i}>10\text{\,s}$. The remaining dots refer to all noise glitches from a
simulation aimed at reproducing the noise ASD of all runs in which isolated
glitches have been observed [20]. Red dots refer to glitches with
$\text{SNR}<5$, while blue dots refer to those with $\text{SNR}\geq 5$ that
would have almost certainly been detected in real data.

The parameter space for these simulations is rather vast, and the criterion
for detectability is not free of ambiguity. Nevertheless, the picture strongly
suggests that:

  1. 1.

there is an abrupt, order of magnitude large jump in the probability density
between the upper edge of the noise glitches distribution and that of the
isolated ones;

  2. 2.

there is a large void in the population of isolated glitches for
$\tau_{i}\gtrsim$200\text{\,}\mathrm{s}$$, between the amplitude of the
observed glitches and that of noise one. This void is real, glitches with
amplitude in the void would be observable with high SNR.

  3. 3.

As already noted, glitches with $\tau_{i}\geq 10^{4}\text{\,s}$, necessary to
reproduce the noise, are missing in the observed sample.

These features are unlikely to disappear as a result of the tuning of any of
the model details.

Thus it may be that the noise is due to a flow of undetectable glitches, but
these appear to belong to a distribution so significantly separated from that
of the observed ones, to make the conclusion that we may be witnessing two
aspects of the same physical phenomenon, highly speculative.

It is also worth noting that, already at the minimal rate of $\lambda\simeq
100-200\text{\,d}^{-1}$, the simulated time series appear Gaussian and
stationary within the statistical uncertainty, making their Poisson nature
undetectable from the data. Thus again the hypothesis that $\Delta g_{e}(t)$
is Poisson in origin, remains a speculation.

## References

  * [1] M. Armano et al.  Beyond the Required LISA Free-Fall Performance: New LISA Pathfinder Results down to $20\text{\,}\mathrm{\SIUnitSymbolMicro Hz}$.  Phys. Rev. Lett., 120:061101, Feb 2018.  DOI: 10.1103/PhysRevLett.120.061101.
  * [2] M. Armano et al.  Sub-Femto-$g$ Free Fall for Space-Based Gravitational Wave Observatories: LISA Pathfinder Results.  Phys. Rev. Lett., 116:231101, 2016.  DOI: 10.1103/PhysRevLett.116.231101.
  * [3] Amaro-Seoane et al.  Laser Interferometer Space Antenna, 2017.  DOI: 10.48550/arXiv.1702.00786.
  * [4] P. McNamara, S. Vitale, K. Danzmann and LISA Pathfinder Sci Working Team.  LISA Pathfinder.  Classical and Quantum Gravity, 25(11), 2008.  DOI: 10.1088/0264-9381/25/11/114034.
  * [5] M. Armano et al.  Sensor Noise in LISA Pathfinder: In-Flight Performance of the Optical Test Mass Readout.  Phys. Rev. Lett., 126:131103, Apr 2021.  DOI: 10.1103/PhysRevLett.126.131103.
  * [6] LISA Pathfinder Collaboration et al.  Capacitive sensing of test mass motion with nanometer precision over millimeter-wide sensing gaps for space-borne gravitational reference sensors.  Physical Review D, 96(6):062004, September 2017.  DOI: 10.1103/PhysRevD.96.062004.
  * [7] M. Armano et al.  Spacecraft and interplanetary contributions to the magnetic environment on-board LISA Pathfinder.  Monthly Notices of the Royal Astronomical Society, 494(2):3014–3027, 04 2020.  DOI: 10.1093/mnras/staa830.
  * [8] M. Armano et al.  Temperature stability in the sub-milliHertz band with LISA Pathfinder.  Monthly Notices of the Royal Astronomical Society, 486(3):3368–3379, 04 2019.  DOI: 10.1093/mnras/stz1017.
  * [9] P. Canizares et al.  The LISA Pathfinder DMU and Radiation Monitor.  Classical and Quantum Gravity, 28(9):094004, 05 2011.  DOI: 10.1088/0264-9381/28/9/094004.
  * [10] M. Armano et al.  Characteristics and Energy Dependence of Recurrent Galactic Cosmic-Ray Flux Depressions and of a Forbush Decrease with LISA Pathfinder.  The Astrophysical Journal, 854(2):113, February 2018.  DOI: 10.3847/1538-4357/aaa774.
  * [11] M. Armano et al.  Constraints on LISA Pathfinder’s self-gravity: design requirements, estimates and testing procedures.  Classical and Quantum Gravity, 33(23):235015, 2016.  DOI: 10.1088/0264-9381/33/23/235015.
  * [12] S. Tibbitts.  High-output paraffin linear motors: utilization in adaptive systems.  page 388, San Diego, CA, 01 1992.  DOI: 10.1117/12.51194.
  * [13] M. Armano et al.  LISA Pathfinder micronewton cold gas thrusters: In-flight characterization.  Phys. Rev. D, 99:122003, Jun 2019.  DOI: 10.1103/PhysRevD.99.122003.
  * [14] L. Ferraioli, M. Hueller and S. Vitale.  Discrete derivative estimation in LISA Pathfinder data reduction.  Classical and Quantum Gravity, 26(9), 2009.  DOI: 10.1088/0264-9381/26/9/094013.
  * [15] LISA Instrument Group.  LISA Performance Model and Error Budget.  Technical note. LISA-LCST-INST-TN-003, 2020.
  * [16] M. Armano et al.  Tilt-to-length coupling in LISA Pathfinder: A data analysis.  Phys. Rev. D, 108:102003, Nov 2023.  DOI: 10.1103/PhysRevD.108.102003.
  * [17] M.S. Hartig and G. Wanner.  Tilt-to-length coupling in LISA Pathfinder: Analytical modeling.  Phys. Rev. D, 108:022008, Jul 2023.  DOI: 10.1103/PhysRevD.108.022008.
  * [18] M. Armano et al.  Calibrating the system dynamics of LISA Pathfinder.  Phys. Rev. D, 97:122002, Jun 2018.  DOI: 10.1103/PhysRevD.97.122002.
  * [19] M. Armano et al.  Analysis of the accuracy of actuation electronics in the laser interferometer space antenna pathfinder.  Review of Scientific Instruments, 91(4):045003, 2020.  DOI: 10.1063/1.5140406.
  * [20] M. Armano et al.  Transient acceleration events in LISA Pathfinder data: Properties and possible physical origin.  Phys. Rev. D, 106:062001, Sep 2022.  DOI: 10.1103/PhysRevD.106.062001.
  * [21] M. Armano et al.  Charging of free-falling test masses in orbit due to cosmic rays: Results from LISA Pathfinder.  Phys. Rev. D, 107:062007, Mar 2023.  DOI: 10.1103/PhysRevD.107.062007.
  * [22] M. Armano et al.  Charge-Induced Force Noise on Free-Falling Test Masses: Results from LISA Pathfinder.  Phys. Rev. Lett., 118:171101, Apr 2017.  DOI: 10.1103/PhysRevLett.118.171101.
  * [23] M. Armano et al.  Measuring the Galactic Cosmic Ray flux with the LISA Pathfinder radiation monitor.  Astroparticle Physics, 98:28–37, 2018.  DOI: 10.1016/j.astropartphys.2018.01.006.
  * [24] M. Armano et al.  Precision charge control for isolated free-falling test masses: LISA pathfinder results.  Phys. Rev. D, 98:062001, Sep 2018.  DOI: 10.1103/PhysRevD.98.062001.
  * [25] A. Cavalleri et al.  Increased Brownian Force Noise from Molecular Impacts in a Constrained Volume.  Phys. Rev. Lett., 103:140601, Sep 2009.  DOI: 10.1103/PhysRevLett.103.140601.
  * [26] H. Jeffreys.  An invariant form for the prior probability in estimation problems.  Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 186(1007):453–461, 1946.  DOI: 10.1098/rspa.1946.0056.
  * [27] K. Jousten, editor.  Handbook of vacuum technology.  Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, second, completely revised and updated edition edition, 2016.  DOI: 10.1002/9783527688265.
  * [28] P. Chiggiato.  Outgassing properties of vacuum materials for particle accelerators, 2020.  DOI: 10.48550/arXiv.2006.07124.
  * [29] LISA Science Study Team.  LISA Science Requirements Document. ESA-L3-EST-SCI-RS-001 1.0, May 2018.
  * [30] H. Akaike.  A new look at the statistical model identification.  IEEE Transactions on Automatic Control, 19(6):716–723, 1974.  DOI: 10.1109/TAC.1974.1100705.
  * [31] A. Gelman et al.  Bayesian Workflow, 2020.  DOI: 10.48550/arXiv.2011.01808.
  * [32] M. Armano et al.  NanoNewton electrostatic force actuators for femtoNewton-sensitive measurements: system performance test in the LISA Pathfinder mission.  2024\.  In print, Phys. Rev. D, DOI: 10.48550/arXiv.2401.00884.
  * [33] P.A. Redhead.  Modeling the pumpdown of a reversibly adsorbed phase. I. Monolayer and submonolayer initial coverage.  Journal of Vacuum Science & Technology A, 13(2):467–475, 1995.  DOI: 10.1116/1.579381.
  * [34] R. Kersevan and M. Ady.  Recent Developments of Monte-Carlo Codes Molflow+ and Synrad+.  In Proc. 10th International Particle Accelerator Conference (IPAC’19), Melbourne, Australia, 19-24 May 2019, number 10 in International Particle Accelerator Conference, pages 1327–1330, Geneva, Switzerland, Jun. 2019. JACoW Publishing.  DOI: doi:10.18429/JACoW-IPAC2019-TUPMP037.
  * [35] M. Armano et al.  Sensor noise in LISA Pathfinder: An extensive in-flight review of the angular and longitudinal interferometric measurement system.  Phys. Rev. D, 106:082001, Oct 2022.  DOI: 10.1103/PhysRevD.106.082001.
  * [36] S. Vitale.  Possible excess GRS sensing noise with TM in science mode. S2-UTN-TN-3151 1.0, 2018.
  * [37] Y. Levin.  Creep events and creep noise in gravitational-wave interferometers: Basic formalism and stationary limit.  Phys. Rev. D, 86:122004, Dec 2012.  DOI: 10.1103/PhysRevD.86.122004.
  * [38] F. Antonucci et al.  Interaction between Stray Electrostatic Fields and a Charged Free-Falling Test Mass.  Physical Review Letters, 108(18):5, 2012.  DOI: 10.1103/PhysRevLett.108.181101.
  * [39] LISA Pathfinder collaboration.  Magnetic-induced force noise in LISA Pathfinder free-falling test masses.  In preparation, 2024.
  * [40] LISA Pathfinder collaboration.  Precision measurements of the magnetic parameters of LISA Pathfinder test masses.  In preparation, 2024.
  * [41] L. Carbone et al.  Thermal gradient-induced forces on geodesic reference masses for LISA.  Phys. Rev. D, 76:102003, Nov 2007.  DOI: 10.1103/PhysRevD.76.102003.
  * [42] B. Kaune.  In-orbit stability analysis of the LISA pathfinder optical metrology: photoreceivers and polarisation.  PhD thesis, Hannover: Gottfried Wilhelm Leibniz Universität, 2021.  DOI: 10.15488/10887.
  * [43] T.W. Darling, F. Rossi, G.I. Opat and G.F. Moorhead.  The fall of charged particles under gravity: A study of experimental problems.  Rev. Mod. Phys., 64:237–257, Jan 1992.  DOI: 10.1103/RevModPhys.64.237.
  * [44] S. Vitale, L. Sala, V. Ferroni and W.J. Weber.  Estimate of force noise from electrostatic patch potentials in lisa pathfinder, 2024.  DOI: https://doi.org/10.48550/arXiv.2404.08340.
  * [45] M. Brownnutt, M. Kumph, P. Rabl and R. Blatt.  Ion-trap measurements of electric-field noise near surfaces.  Rev. Mod. Phys., 87:1419–1482, Dec 2015.  DOI: 10.1103/RevModPhys.87.1419.
  * [46] G. Russano et al.  Measuring fN force variations in the presence of constant nN forces: a torsion pendulum ground test of the LISA Pathfinder free-fall mode.  Classical and Quantum Gravity, 35(3):035017, 02 2018.  DOI: 10.1088/1361-6382/aaa00f.
  * [47] H.D. Wealthy and W. C.  Spacecraft-Level Gravity Model Results. S2.ASU.TN.2270.ISS02, 2008.
  * [48] B. Schläppi et al.  Influence of spacecraft outgassing on the exploration of tenuous atmospheres with in situ mass spectrometry.  Journal of Geophysical Research: Space Physics, 115(A12), 2010.  DOI: 10.1029/2010JA015734.
  * [49] F. Antonucci et al.  From laboratory experiments to LISA Pathfinder: achieving LISA geodesic motion.  Classical and Quantum Gravity, 28(9):094002, 2011.  DOI: 10.1088/0264-9381/28/9/094002.
  * [50] C. Zanoni and V. Ferroni.  Electromagnetic shielding effect of the Electrode Housing. LISA-UTN-INST-TN-022 Iss. 1.1, 2024.
  * [51] H.D. Wealthy and T. C.  Prediction of spacecraft magnetic field based on recent testing. S2.ASU.TN.2523.ISS03, 2011.
  * [52] G. Anderson et al.  Experimental results from the ST7 mission on LISA Pathfinder.  Phys. Rev. D, 98:102005, Nov 2018.  DOI: 10.1103/PhysRevD.98.102005.
  * [53] S. Vitale et al.  Data series subtraction with unknown and unmodeled background noise.  Phys. Rev. D, 90:042003, 2014.  DOI: 10.1103/PhysRevD.90.042003.
  * [54] P.D. Welch.  The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms.  IEEE Trans. Audio and Electroacoustics, 15(2):70–73, 1967.  DOI: 10.1109/TAU.1967.1161901.
  * [55] N.R. Goodman.  Statistical Analysis Based on a Certain Multivariate Complex Gaussian Distribution (An Introduction).  The Annals of Mathematical Statistics, 34(1):152–177, March 1963.  DOI: 10.1214/aoms/1177704250.
  * [56] P. Shaman.  The inverted complex Wishart distribution and its application to spectral estimation.  J. Multivar. Anal., 10(1):51–59, March 1980.  DOI: 10.1016/0047-259X(80)90081-0.
  * [57] L. Svensson and M. Lundberg.  On posterior distributions for signals in Gaussian noise with unknown covariance matrix.  IEEE Transactions on Signal Processing, 53(9):3554–3571, 2005.  DOI: 10.1109/TSP.2005.853102.
  * [58] C. Villegas.  On the A Priori Distribution of the Covariance Matrix.  The Annals of Mathematical Statistics, 40(3):1098 – 1099, 1969.  DOI: 10.1214/aoms/1177697617.
  * [59] J. Ellis and R. van Haasteren.  jellis18/PTMCMCSampler: Official Release, October 2017.  DOI: 10.5281/zenodo.1037579.
  * [60] P. Sarra.  LISA Technology Package Core Assembly Gravitational Analysis Report. S2-CGS-RP-3001. Iss10, 2015.
  * [61] M. Armano et al.  Temperature-induced forces on geodesic reference masses for LISA: results from LISA Pathfinder.  In preparation, 2024.