aimed at estimating the TM couplings, the sensing cross-talk and other
relevant parameters needed for system calibration. Such calibration is also
critical for the unbiased estimation of the total equivalent acceleration
noise. The same model is planned to be employed during the identification
experiments of the LPF mission.

###  3.4 Cross-talk from degrees of freedom other than the optical axis

Section 3.1 has introduced a general formalism – the three master equations of
(3.4) – capable in describing the evolution of LPF as a physical system
continuously subjected to a digital control. Subsequently, in Section 3.3 the
formalism has been applied in the derivation of a model of LPF along the two
main optical degrees of freedom. However, there are sources of non-idealities
in this description. The knowledge of the operators $\bm{D}$, $\bm{A}$ and
$\bm{S}$ might lack because of a poor calibration, misalignments, pointing
inaccuracies as in Section 2.4, force gradients and force noise along
different axes. Such errors couple with the main interferometric axis and
manifest themselves as a cross-talk from other degrees of freedom. Despite
describing the system in its full complexity with all degrees of freedom at
once, in the hypothesis of small motion, the absence of strong non-linearities
assures that the contribution from other degrees of freedom is a small
perturbation to the nominal dynamics along the optical axis.

For the sake of clarity, the approach used here in describing such effects is
based on a first-order perturbation theory, where the zeroth-order dynamics,
named the nominal dynamics, is fully known. The readouts $\bm{o}$ can then be
split into a nominal response of the system to the $x$-dynamics, say
$\bm{o}_{0}$, and a small perturbation coming from other degrees of freedom
($y$, $z$ or some torsional angles), say $\delta\bm{o}$, therefore

$\bm{o}\simeq\bm{o}_{0}+\delta\bm{o}\leavevmode\nobreak\ .$ (3.26)

The leading correcting terms can be embedded into the dynamics as
imperfections to the operators introduced in the previous section. For
example, if $\bm{D}_{0}$ is the nominal dynamics operator, $\delta\bm{D}$ is
the relative imperfection.

To first order the dynamics is a generalization of (3.4) [45]

$\displaystyle(\bm{D}_{0}+\delta\bm{D})\,\bm{q}$
$\displaystyle=\bm{g}\leavevmode\nobreak\ ,$ (3.27a) $\displaystyle\bm{g}$
$\displaystyle=\bm{f}_{\text{n}}+(\bm{A}_{0}+\delta\bm{A})\left[\bm{f}_{\text{i}}-\bm{C}_{0}\,(\bm{o}-\bm{o}_{\text{i}})\right]\leavevmode\nobreak\
,$ (3.27b) $\displaystyle\bm{o}$
$\displaystyle=(\bm{S}_{0}+\delta\bm{S})\,\bm{q}+\bm{o}_{\text{n}}\leavevmode\nobreak\
.$ (3.27c)

where the control laws in $\bm{C}_{0}$ are exactly known from the original
design. Figure 3.6 shows a block diagram of the above set of equations, in
where a dashed arrow stands for a cross-talk contribution. Every time a
physical quantity must be passed throughout an operator, then the relative
imperfection mixes it up among many degrees of freedom.

Figure 3.6: An update of the block diagram shown in Figure 3.4 where the
cross-talks (dynamics, sensing and actuation) are introduced as imperfections
(dashed arrows) to the nominal operators. For example,
$\bm{D}\simeq\bm{D}_{0}+\delta\bm{D}$.

The nominal solution $\bm{o}_{0}$ is provided by (3.6) in which all operators
have the proper subscript indicating a zeroth-order term. The differential
operator is obtained from the definition (3.7) by means of simple algebra and
using the inversion lemma
111$\left(\bm{S}_{0}+\delta\bm{S}\right)^{-1}=\bm{S}_{0}^{-1}-\bm{S}_{0}^{-1}\,\delta\bm{S}\left(\bm{1}+\bm{S}_{0}^{-1}\,\delta\bm{S}\right)^{-1}\,\bm{S}_{0}^{-1}\simeq\bm{S}_{0}^{-1}-\bm{S}_{0}^{-1}\,\delta\bm{S}\,\bm{S}_{0}^{-1}$.

$\begin{split}\bm{\Delta}&\simeq\left(\bm{D}_{0}+\delta\bm{D}\right)\left(\bm{S}_{0}+\delta\bm{S}\right)^{-1}+\left(\bm{A}_{0}+\delta\bm{A}\right)\bm{C}_{0}\\\
&\simeq\left(\bm{D}_{0}+\delta\bm{D}\right)\bm{S}_{0}^{-1}\left(\bm{1}-\delta\bm{S}\,\bm{S}_{0}^{-1}\right)+\left(\bm{A}_{0}+\delta\bm{A}\right)\bm{C}_{0}\\\
&\simeq\left(\bm{D}_{0}\,\bm{S}_{0}^{-1}+\bm{A}_{0}\,\bm{C}_{0}\right)+\left(\delta\bm{D}\,\bm{S}_{0}^{-1}-\bm{D}_{0}\,\bm{S}_{0}^{-1}\,\delta\bm{S}\,\bm{S}_{0}^{-1}+\delta\bm{A}\,\bm{C}_{0}\right)\end{split}$
(3.28)

where the nominal operator and its imperfection are defined as

$\displaystyle\bm{\Delta}_{0}$
$\displaystyle=\bm{D}_{0}\,\bm{S}_{0}^{-1}+\bm{A}_{0}\,\bm{C}_{0}\leavevmode\nobreak\
,$ (3.29a) $\displaystyle\delta\bm{\Delta}$
$\displaystyle=\delta\bm{D}\,\bm{S}_{0}^{-1}-\bm{D}_{0}\,\bm{S}_{0}^{-1}\,\delta\bm{S}\,\bm{S}_{0}^{-1}+\delta\bm{A}\,\bm{C}_{0}\leavevmode\nobreak\
.$ (3.29b)

There are two ways to obtain the evolution of the first-order correction
$\delta\bm{o}$. The first method involves a direct computation on (3.27)
similar to the reasoning of Section 3.1. The idea is to combine the three
equations, substitute back the zeroth-order and keep only first-order terms.
The second method is more straightforward and is based on the elaboration of
the equation of motion (3.6). Following this idea, the expansion of (3.6) is

$\begin{split}\left(\bm{\Delta}_{0}+\delta\bm{\Delta}\right)\left(\bm{o}_{0}+\delta\bm{o}\right)&\simeq\bm{f}_{\text{n}}\\\
&\quad+\left(\bm{D}_{0}+\delta\bm{D}\right)\left(\bm{S}_{0}+\delta\bm{S}\right)^{-1}\,\bm{o}_{\text{n}}\\\
&\quad+\left(\bm{A}_{0}+\delta\bm{A}\right)\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)\leavevmode\nobreak\
,\end{split}$ (3.30)

Now, to first order the equation becomes

$\begin{split}\bm{\Delta}_{0}\,\bm{o}_{0}+\bm{\Delta}_{0}\,\delta\bm{o}+\delta\bm{\Delta}\,\bm{o}_{0}&\simeq\bm{f}_{\text{n}}\\\
&\quad+\left(\bm{D}_{0}\,\bm{S}_{0}^{-1}+\delta\bm{D}_{0}\,\bm{S}_{0}^{-1}-\bm{D}_{0}\,\bm{S}_{0}^{-1}\,\delta\bm{S}\,\bm{S}_{0}^{-1}\right)\bm{o}_{\text{n}}\\\
&\quad+\left(\bm{A}_{0}+\delta\bm{A}\right)\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)\leavevmode\nobreak\
,\end{split}$ (3.31)

which is further simplified by assuming the knowledge of the nominal dynamics.
Indeed, $\bm{\Delta}_{0}\,\bm{o}_{0}$ on the left-hand side is canceled out by
$\bm{f}_{\text{n}}+\bm{D}_{0}\,\bm{S}_{0}^{-1}\,\bm{o}_{\text{n}}+\bm{A}_{0}\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)$
on the right-hand side. Therefore, the first-order dynamics is given by

$\begin{split}\bm{\Delta}_{0}\,\delta\bm{o}&\simeq-\delta\bm{\Delta}\,\bm{o}_{0}+\left(\delta\bm{D}\,\bm{S}_{0}^{-1}-\bm{D}_{0}\,\bm{S}_{0}^{-1}\,\delta\bm{S}\,\bm{S}_{0}^{-1}\right)\bm{o}_{\text{n}}\\\
&\quad+\delta\bm{A}\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)\leavevmode\nobreak\
,\end{split}$ (3.32)

where the operator in front of $\bm{o}_{\text{n}}$ is exactly
$\delta\bm{\Delta}-\delta\bm{A}\,\bm{C}_{0}$ (it can be checked by comparing
to (3.29b)). The equation also shows that the first-order dynamics is given in
terms of the zeroth-order $\bm{o}_{0}$. To write down an explicit form of the
above formula, it is necessary to substitute the expression of $\bm{o}_{0}$ in
(3.6) and get

$\begin{split}\bm{\Delta}_{0}\,\delta\bm{o}&\simeq-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\left[\bm{f}_{\text{n}}+\bm{D}_{0}\,\bm{S}_{0}^{-1}\,\bm{o}_{\text{n}}+\bm{A}_{0}\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)\right]\\\
&\quad+\left(\delta\bm{\Delta}-\delta\bm{A}\,\bm{C}\right)\bm{o}_{\text{n}}+\delta\bm{A}\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)\\\
&=-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\,\bm{f}_{\text{n}}+\left(\delta\bm{\Delta}-\delta\bm{A}\,\bm{C}_{0}-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\,\bm{D}_{0}\,\bm{S}_{0}^{-1}\right)\bm{o}_{\text{n}}\\\
&\quad+\left(\delta\bm{A}-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\,\bm{A}_{0}\right)\left(\bm{f}_{\text{i}}+\bm{C}_{0}\,\bm{o}_{\text{i}}\right)\leavevmode\nobreak\
.\end{split}$ (3.33)

The operator in front of $\bm{f}_{\text{n}}$ is the force noise cross-talk,
whose transfer to total equivalent acceleration is denoted with

$\delta\bm{T}_{f_{\text{n}}\rightarrow
f}=-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\leavevmode\nobreak\ ;$ (3.34)

the operator in front of $\bm{o}_{\text{n}}$ is the readout noise cross-talk
and can be rewritten as

$\begin{split}\delta\bm{T}_{o_{\text{n}}\rightarrow
f}&=\delta\bm{\Delta}-\delta\bm{A}\,\bm{C}_{0}-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\,\bm{D}_{0}\,\bm{S}_{0}^{-1}\\\
&=\delta\bm{\Delta}-\delta\bm{A}\,\bm{C}_{0}-\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\left(\bm{\Delta}_{0}-\bm{A}_{0}\,\bm{C}_{0}\right)\\\
&=\delta\bm{\Delta}\,\bm{\Delta}_{0}^{-1}\,\bm{A}_{0}\,\bm{C}_{0}-\delta\bm{A}\,\bm{C}_{0}\\\
&=-\left(\delta\bm{A}+\delta\bm{T}_{f_{\text{n}}\rightarrow
f}\,\bm{A}_{0}\right)\bm{C}_{0}\leavevmode\nobreak\ ;\end{split}$ (3.35)

the operator in front of $\bm{f}_{\text{i}}$ is the force actuation cross-talk

$\delta\bm{T}_{f_{\text{i}}\rightarrow
f}=\delta\bm{A}+\delta\bm{T}_{f_{\text{n}}\rightarrow
f}\,\bm{A}_{0}\leavevmode\nobreak\ ;$ (3.36)

and the operator in front of $\bm{o}_{\text{i}}$ is the bias actuation cross-
talk

$\delta\bm{T}_{o_{\text{i}}\rightarrow
f}=\left(\delta\bm{A}+\delta\bm{T}_{f_{\text{n}}\rightarrow
f}\,\bm{A}_{0}\right)\bm{C}_{0}\leavevmode\nobreak\ ,$ (3.37)

that is exactly the readout noise cross-talk modulo a sign. As reasonable to
expect, the symmetry between the various cross-talk terms comes from the
linearization of the problem.

The equation of motion for the first-order cross-talk is finally given by

$\bm{\Delta}_{0}\,\delta\bm{o}\simeq\delta\bm{T}_{f_{\text{n}}\rightarrow
f}\,\bm{f}_{\text{n}}+\delta\bm{T}_{o_{\text{n}}\rightarrow
f}\,\bm{o}_{\text{n}}+\delta\bm{T}_{f_{\text{i}}\rightarrow
f}\,\bm{f}_{\text{i}}+\delta\bm{T}_{o_{\text{i}}\rightarrow
f}\,\bm{o}_{\text{i}}\leavevmode\nobreak\ ,$ (3.38)

that enlightens the various contributions to the overall cross-talk from other
degrees of freedom. On one side, in a pure noise measurement during the LPF
mission (no application of forces or biases) a noise cross-talk sums up to the
nominal dynamics along the main degrees of freedom. On the other side,
whatever a bias is applied to the system, a perturbation is produced along the
sensitive axis. As it is clear, a non-negligible cross-talk at any level of
the system (dynamics, sensing and actuation) actually breaks the nominal
orthogonality between different degrees of freedom.

###  3.5 Dynamical model for $xy$ cross-talk

The application of the idea of the preceding section is a further development
of the model along the optical axis described in Section 3.3, i.e., the cross-
talk from other nominally orthogonal degrees of freedom to the optical axis.
The example discussed in this section is the cross-talk from the $y$ degree of
freedom to $x$.

Referring to Figure 3.2 for the usual coordinate naming convention, the $xy$
cross-talk can be viewed as a first-order perturbation of the dynamics in the
$xy$ plane to the zeroth-order dynamics along $x$. Obviously, the dynamics in
$xy$ contains the rotation $\phi$ about the $z$ axis.

The control design of the main science mode for this simplified model requires
a minimal number of drag-free, electrostatic suspension and attitude
coordinates as inputs to the same control loops. Table 3.2 presents a list of
those coordinates relevant for the $xy$ cross-talk as taken from Table 3.1.

Table 3.2:  List of the 7 controlled degrees of freedom for the $xy$ cross-talk of the LPF mission in the main science mode. Refer to Table 3.1 for a comprehensive description of all coordinates. Notice that $o_{\phi_{2}}=o_{\phi_{1}}+o_{\phi_{12}}$ is used in the equations for clearness and simplicity. Coordinate | Control | Sensor | Actuator
---|---|---|---
$x_{1}$ | Drag-free | $o_{1}\,=\,\text{IFO}[x_{1}]$ | FEEP
$y_{1}$ | Drag-free | $o_{y_{1}}\,=\,\text{GRS}[y_{1}]$ | FEEP
$\phi_{1}$ | Elect. suspension | $o_{\phi_{1}}\,=\,\text{IFO}[\phi_{1}]$ | GRS
$x_{2}$ | Elect. suspension | $o_{12}\,=\,\text{IFO}[x_{12}]$ | GRS
$y_{2}$ | Drag-free | $o_{y_{2}}\,=\,\text{GRS}[y_{2}]$ | FEEP
$\phi_{2}$ | Elect. suspension | $o_{\phi_{12}}\,=\,\text{IFO}[\phi_{12}]$ | GRS
$\phi_{\text{SC}}$ | Attitude | $o_{\phi_{\text{SC}}}\,=\,\text{ST}[\phi_{\text{SC}}]$ | GRS

The control is such that:

  1. 1.

along $x$: guided by the optical $x_{1}$, the SC is forced to follow the
reference TM through thruster actuation; guided by the optical $x_{12}$, the
second TM is forced to follow the reference TM through capacitive actuation;

  2. 2.

along $y$: guided by the capacitive $(y_{1}+y_{2})/2$, the SC is forced to
follow both TMs through thruster actuation; guided by the star-tracker
$\phi_{\text{SC}}$, the TMs are oriented along $\phi$;

  3. 3.

along $\phi$: guided by the capacitive $(y_{1}-y_{2})/2$, the SC is forced to
follow both TMs through thruster actuation; guided by the optical $\phi_{1}$
and $\phi_{2}$, both TMs are oriented along $\phi$ through capacitive
actuation.

As already pointed out, the cross-talk can be described by a first-order
perturbation of the nominal dynamics along $x$. Three different types of
cross-talks can be identified:

  1. 1.

dynamical cross-talk;

  2. 2.

actuation cross-talk;

  3. 3.

sensing cross-talk.

All equations must be written to within linear terms of an imperfection or
noise contribution.

Concerning the dynamics along the $x$ axis, in the approximation of small
motion, the stiffness constant has been introduced to model residual
oscillator-like couplings between the TMs and the SC. The generalization in
three dimensions is straightforward. Since there are electrodes all around the
TMs and the most important coupling is indeed due to electrostatics, in place
of a single oscillator, 6 coupled harmonic oscillators along the translational
and rotational degrees of freedom must be considered. Therefore, the stiffness
constant becomes a $6\times 6$ quasi-diagonal matrix, the stiffness matrix.
For the $xy$ cross-talk, since it makes sense to inspect the cross-stiffness
from $y$ to $x$ or $\phi$ to $x$, the structure of the matrix for the first TM
is

$\bm{\kappa}_{1}=\begin{pmatrix}m_{1}\,\omega_{1,x}^{2}&\delta_{1,xy}\,m_{1}\,\omega_{1,y}^{2}&\delta_{1,x\phi}\,I_{1,z}\,\omega_{1,\phi}^{2}\\\
0&m_{1}\,\omega_{1,y}^{2}&0\\\
0&0&I_{1,z}\,\omega_{1,\phi}^{2}\end{pmatrix}\leavevmode\nobreak\ ,$ (3.39)

where $m_{1}=1.96\,\mathrm{kg}$ is the TM mass,
$I_{1,z}=\nicefrac{{1}}{{6}}\,m_{1}\,l^{2}\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}7\\!\times\\!10^{-4}\,\mathrm{kg\,m^{2}}$ is the
inertia matrix about $z$ and $l=46\,\mathrm{mm}$ is the TM size. A
$\delta$-coefficient, with the number of the TM and the names of two
coordinates as subscripts, denotes a dynamical cross-talk imperfection
typically $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$
}$$}1\\!\times\\!10^{-3}$, an assumption based on on-ground measurements and
theoretical models of the GRS.

The second type is the actuation cross-talk due to misalignments in the setup
of the thrusters and electrodes. The idea is that every time a force is
actuated along a nominally orthogonal degree of freedom it couples with $x$. A
$\delta$-coefficient, with the names of the actuation and the coordinate,
denotes an actuation cross-talk imperfection.

Finally, the third type is the sensing cross-talk due to miscalibrations and
misalignments in the sensors. For the capacitive sensing, misalignments in the
electrodes produces a mixing in the sensed coordinates. For the optical
sensing, misalignments in the optical elements produce an analogous result. As
in Section 3.3, the 7 physical coordinates

$\bm{q}=\scalebox{0.8}{$\begin{pmatrix}x_{1}\\\ x_{12}\\\ y_{1}\\\ y_{2}\\\
\phi_{1}\\\ \phi_{2}\\\ \phi_{\text{SC}}\end{pmatrix}$}\leavevmode\nobreak\ ,$
(3.40)

are converted into the sensed coordinates

$\bm{o}=\scalebox{0.8}{$\begin{pmatrix}o_{1}\\\ o_{12}\\\ o_{y_{1}}\\\
o_{y_{2}}\\\ o_{\phi_{1}}\\\ o_{\phi_{2}}\\\
o_{\phi_{\text{SC}}}\end{pmatrix}$}\leavevmode\nobreak\ ,$ (3.41)

with an inevitable mixing. Assuming a sensing operator which is nominally
identity, the relative imperfection operator has the following matrix
representation

$\delta\bm{S}=\left(\begin{smallmatrix}0&0&\delta_{S,1y_{1}}&\delta_{S,1y_{2}}&\frac{l}{2}\,\delta_{S,1\phi_{1}}&\frac{l}{2}\,\delta_{S,1\phi_{2}}&\frac{L}{2}\,\delta_{S,1\phi_{\text{SC}}}\\\
S_{21}&0&\delta_{S,2y_{1}}&\delta_{S,2y_{2}}&\frac{l}{2}\,\delta_{S,2\phi_{1}}&\frac{l}{2}\,\delta_{S,2\phi_{2}}&\frac{L}{2}\,\delta_{S,2\phi_{\text{SC}}}\\\
0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0\end{smallmatrix}\right)\leavevmode\nobreak\ ,$ (3.42)

where $l=46\,\mathrm{mm}$ is the TM size and $L=38\,\mathrm{cm}$ is the
nominal separation between the TM centers of mass. $S_{21}$ is the sensing
cross-talk between $o_{1}$ and $o_{12}$, already discussed, and a generic
$\delta$-coefficient denotes a sensing cross-talk imperfection from an
orthogonal degree of freedom to $x$. Since the target of the investigation is
the cross-talk to the two optical degrees of freedom, no cross-talk to other
coordinates is considered as this is second order.

In the derivation of the equations of motion the gravity/torque gradients
between the TMs can be neglected without loss of generality, as well as the
back-reactions, since they are all second order effects. The first set of
equations is along the $x$ degree of freedom. With the same considerations of
Section 3.3, i.e., assuming small motion, small forces and null initial
conditions, the linear equations of motion, per unit mass, in frequency domain
are

$\displaystyle\begin{split}s^{2}\,x_{1}+s^{2}\,x_{\text{SC}}+\omega_{1,x}^{2}\,x_{1}\quad&\\\
\quad+\,\delta_{1,xy}\,\omega_{1,y}^{2}\,y_{1}+\delta_{1,x\phi}\,\frac{l}{2}\,\omega_{1,\phi}^{2}\,\phi_{1}&=f_{1,x}\\\
&\quad+\delta_{\text{sus},\,y_{1}}\,\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad+\delta_{\text{sus},\,\phi_{1}}\,\frac{l}{2}\,C_{\text{sus},\,\phi_{1}}(s)\,o_{\phi_{1}}\leavevmode\nobreak\
,\end{split}$ (3.43a)
$\displaystyle\begin{split}s^{2}\,x_{2}+s^{2}\,x_{\text{SC}}+\omega_{2,x}^{2}\,x_{2}\quad&\\\
\quad+\,\delta_{2,xy}\,\omega_{2,y}^{2}\,y_{2}+\delta_{2,x\phi}\,\frac{l}{2}\,\omega_{2,\phi}^{2}\,\phi_{2}&=f_{2,x}-C_{\text{sus},\,x}(s)\,o_{12}\\\
&\quad-\delta_{\text{sus},\,y_{2}}\,\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad+\delta_{\text{sus},\,\phi_{2}}\,\frac{l}{2}\,C_{\text{sus},\,\phi_{2}}(s)\,o_{\phi_{2}}\leavevmode\nobreak\
,\end{split}$ (3.43b)
$\displaystyle\begin{split}s^{2}\,x_{\text{SC}}&=f_{\text{SC},x}+C_{\text{df},\,x}(s)\,o_{1}\\\
&\quad+\delta_{\text{df},\,y_{\text{SC}}}\,\frac{1}{2}\left[C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}+C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}\right]\\\
&\quad-\delta_{\text{df},\,\phi_{\text{SC}}}\,\frac{1}{2}\left[C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}-C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}\right]\leavevmode\nobreak\
,\end{split}$ (3.43c)

The second row of the left-hand side of the equations for the TMs contain the
dynamical cross-talk due to the stiffness matrix (imperfections
$\delta_{1,xy}$, $\delta_{1,x\phi}$, $\delta_{2,xy}$ and $\delta_{2,x\phi}$).
Instead, the last two rows of the right-hand side of the equations are
actuation cross-talk. Notice that for the TMs there is a cross-talk from the
SC inertial attitude control ($\delta_{\text{sus},\,y}$) and one from the TM
attitude control ($\delta_{\text{sus},\,\phi}$). For the SC there is a cross-
talk from the $y$ and $\phi$ drag-free actuation
($\delta_{\text{df},\,y_{\text{SC}}}$ and
$\delta_{\text{df},\,\phi_{\text{SC}}}$). Other relevant quantities have been
already defined in the previous section and are not further discussed here.

The second set of equations describes the dynamics along the nominally
orthogonal $y$ axis

$\displaystyle
s^{2}\,y_{1}+s^{2}\,y_{\text{SC}}+\omega_{1,y}^{2}\,y_{1}-\frac{L}{2}\,s^{2}\,\phi_{\text{SC}}$
$\displaystyle=f_{1,y}+\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\leavevmode\nobreak\
,$ (3.44a) $\displaystyle
s^{2}\,y_{2}+s^{2}\,y_{\text{SC}}+\omega_{2,y}^{2}\,y_{2}+\frac{L}{2}\,s^{2}\,\phi_{\text{SC}}$
$\displaystyle=f_{2,y}-\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\leavevmode\nobreak\
,$ (3.44b) $\displaystyle s^{2}\,y_{\text{SC}}$
$\displaystyle=f_{\text{SC},y}+\frac{1}{2}\left[C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}+C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}\right]\leavevmode\nobreak\
.$ (3.44c)

The inertial attitude control is implemented as electrostatic suspension
actuation on the TMs along $y$ (notice the opposite signs) through the
$C_{\text{sus},\,\phi_{\text{SC}}}(s)$ control law. $C_{\text{df},\,y_{1}}(s)$
and $C_{\text{df},\,y_{2}}(s)$ are the two drag-free control laws along $y$.
On the left-end side of the equations the SC absolute angular acceleration
$s^{2}\,\phi_{\text{SC}}$ also appears as a strict consequence of the coupling
between the $\phi$ motion with $y$.

The third set of equations describes the dynamics along the nominally
orthogonal $\phi$ angle

$\displaystyle
s^{2}\,\phi_{1}+s^{2}\,\phi_{\text{SC}}+\omega_{1,\phi}^{2}\,\phi_{1}$
$\displaystyle=\tau_{1,z}-C_{\text{sus},\,\phi_{1}}(s)\,o_{\phi_{1}}\leavevmode\nobreak\
,$ (3.45a) $\displaystyle
s^{2}\,\phi_{2}+s^{2}\,\phi_{\text{SC}}+\omega_{2,\phi}^{2}\,\phi_{2}$
$\displaystyle=\tau_{2,z}-C_{\text{sus},\,\phi_{2}}(s)\,o_{\phi_{2}}\leavevmode\nobreak\
,$ (3.45b) $\displaystyle s^{2}\,\phi_{\text{SC}}$
$\displaystyle=\tau_{\text{SC},z}+\frac{1}{L}\left[C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}-C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}\right]\leavevmode\nobreak\
,$ (3.45c)

where $\tau$ denotes a generic component of torque per unit of inertia
$I_{\text{TM},z}\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$
}$$}7\\!\times\\!10^{-4}\,\mathrm{kg\,m^{2}}$ and $I_{\text{SC},z}\raise
0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$
}$$}1\\!\times\\!10^{3}\,\mathrm{kg\,m^{2}}$ [46]. On the left-end side, the
angular stiffness constants are evident. On the right-end side, there are the
$C_{\text{sus},\,\phi}(s)$ electrostatic suspension control law and the
$C_{\text{df},\,y}(s)$ drag-free control law. The SC attitude is controlled by
actuating along the differential TM positions along $y$.

As said, along $y$ the SC follows the average $y$ motion of the TMs sensed
with the capacitive $o_{y_{1}}$ and $o_{y_{2}}$; the TMs are oriented
following the star-tracker $o_{\phi_{\text{SC}}}$. Analogously, along $\phi$
the SC is oriented following the differential $y$ motion of the TMs sensed
with the capacitive $o_{y_{1}}$ and $o_{y_{2}}$; the TMs are oriented
following the optical $o_{\phi_{1}}$ and $o_{\phi_{2}}$.

The SC absolute linear acceleration is not measurable, whereas the SC absolute
attitude is measured by the ST with respect to the celestial inertial frame.
Therefore, the system of 9 equations turns into 7 equations because the SC
acceleration along $x$ and $y$ must be canceled out. By doing so in (3.43) and
defining the differential TM displacement, the dynamics along $x$ becomes

$\displaystyle\begin{split}s^{2}\,x_{1}+\omega_{1,x}^{2}\,x_{1}\quad&\\\
\quad+\,\delta_{1,xy}\,\omega_{1,y}^{2}\,y_{1}+\delta_{1,x\phi}\,\frac{l}{2}\,\omega_{1,\phi}^{2}\,\phi_{1}&=f_{1,x}-f_{\text{SC},x}-C_{\text{df},\,x}(s)\,o_{1}\\\
&\quad+\delta_{\text{sus},\,y_{1}}\,\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad+\delta_{\text{sus},\,\phi_{1}}\,\frac{l}{2}\,C_{\text{sus},\,\phi_{1}}(s)\,o_{\phi_{1}}\\\
&\quad-\delta_{\text{df},\,y_{\text{SC}}}\,\frac{1}{2}\left[C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}+C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}\right]\\\
&\quad+\delta_{\text{df},\,\phi_{\text{SC}}}\,\frac{1}{2}\left[C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}-C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}\right]\leavevmode\nobreak\
,\end{split}$ (3.46a)
$\displaystyle\begin{split}s^{2}\,x_{12}+\omega_{2,x}^{2}\,x_{2}-\omega_{1,x}^{2}\,x_{1}\quad&\\\
\quad+\,\delta_{2,xy}\,\omega_{2,y}^{2}\,y_{2}+\delta_{2,x\phi}\,\frac{l}{2}\,\omega_{2,\phi}^{2}\,\phi_{2}\\\
\quad-\,\delta_{1,xy}\,\omega_{1,y}^{2}\,y_{1}-\delta_{1,x\phi}\,\frac{l}{2}\,\omega_{1,\phi}^{2}\,\phi_{1}&=f_{2,x}-C_{\text{sus},\,x}(s)\,o_{12}\\\
&\quad-\delta_{\text{sus},\,y_{2}}\,\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad+\delta_{\text{sus},\,\phi_{2}}\,\frac{l}{2}\,C_{\text{sus},\,\phi_{2}}(s)\,o_{\phi_{2}}\\\
&\quad-\delta_{\text{sus},\,y_{1}}\,\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad-\delta_{\text{sus},\,\phi_{1}}\,\frac{l}{2}\,C_{\text{sus},\,\phi_{1}}(s)\,o_{\phi_{1}}\leavevmode\nobreak\
.\end{split}$ (3.46b)

Analogously, substituting the SC acceleration along $y$ and $\phi$ into the
dynamics along $y$

$\displaystyle\begin{split}s^{2}\,y_{1}+\omega_{1,y}^{2}\,y_{1}&=f_{1,y}-f_{\text{SC},y}+\frac{L}{2}\tau_{\text{SC},z}+\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad-C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}\leavevmode\nobreak\ ,\end{split}$
(3.47a)
$\displaystyle\begin{split}s^{2}\,y_{2}+\omega_{2,y}^{2}\,y_{2}&=f_{2,y}-f_{\text{SC},y}-\frac{L}{2}\tau_{\text{SC},z}-\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\,o_{\phi_{\text{SC}}}\\\
&\quad-C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}\leavevmode\nobreak\ .\end{split}$
(3.47b)

$y_{1}$ and $y_{2}$ are named drag-free coordinates as they guide a drag-free
actuation; $\phi_{\text{SC}}$ is an electrostatic suspension coordinate as it
guides a capacitive actuation. At the same time the dynamics along $\phi$ is

$\displaystyle\begin{split}s^{2}\,\phi_{1}+\omega_{1,\phi}^{2}\,\phi_{1}&=\tau_{1,z}-\tau_{\text{SC},z}-C_{\text{sus},\,\phi_{1}}(s)\,o_{\phi_{1}}\\\
&\quad-\frac{1}{L}\left[C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}-C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}\right]\leavevmode\nobreak\
,\end{split}$ (3.48a)
$\displaystyle\begin{split}s^{2}\,\phi_{2}+\omega_{2,\phi}^{2}\,\phi_{2}&=\tau_{2,z}-\tau_{\text{SC},z}-C_{\text{sus},\,\phi_{2}}(s)\,o_{\phi_{2}}\\\
&\quad-\frac{1}{L}\left[C_{\text{df},\,y_{2}}(s)\,o_{y_{2}}-C_{\text{df},\,y_{1}}(s)\,o_{y_{1}}\right]\leavevmode\nobreak\
.\end{split}$ (3.48b)

and the one for $\phi_{\text{SC}}$ which remains unchanged. $\phi_{1}$ and
$\phi_{2}$ are electrostatic suspension coordinates as they guide a capacitive
actuation.

The equations of motion presented above can now be mapped to the formalism of
the previous section. Therefore, the $7\times 7$ nominal dynamics operator has
the following matrix representation

$\bm{D}_{0}=\scalebox{0.8}{$\begin{pmatrix}s^{2}+\omega_{1,x}^{2}&0&0&0&0&0&0\\\
\omega_{12,x}^{2}&s^{2}+\omega_{1,x}^{2}+\omega_{12,x}^{2}&0&0&0&0&0\\\
0&0&s^{2}+\omega_{1,y}^{2}&0&0&0&0\\\ 0&0&0&s^{2}+\omega_{2,y}^{2}&0&0&0\\\
0&0&0&0&s^{2}+\omega_{1,\phi}^{2}&0&0\\\
0&0&0&0&0&s^{2}+\omega_{2,\phi}^{2}&0\\\
0&0&0&0&0&0&s^{2}\end{pmatrix}$}\leavevmode\nobreak\ ,$ (3.49)

which is a natural generalization of the same operator (3.18) written for the
model along $x$. The first-order perturbation due to the dynamical cross-talk
is given by

$\delta\bm{D}=\scalebox{0.8}{$\begin{pmatrix}0&0&\delta_{1,xy}\,\omega_{1,y}^{2}&0&\delta_{1,x\phi}\,\frac{l}{2}\,\omega_{1,\phi}^{2}&0&0\\\
0&0&-\delta_{1,xy}\,\omega_{1,y}^{2}&\delta_{2,xy}\,\omega_{2,y}^{2}&-\delta_{1,x\phi}\,\frac{l}{2}\,\omega_{1,\phi}^{2}&\delta_{2,x\phi}\,\frac{l}{2}\,\omega_{2,\phi}^{2}&0\\\
0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0\end{pmatrix}$}\leavevmode\nobreak\ .$ (3.50)

Analogously, the control operator can be identified as

$\bm{C}_{0}=\scalebox{0.8}{$\begin{pmatrix}C_{\text{df},\,x}(s)&0&0&0&0&0&0\\\
0&C_{\text{sus},\,x}(s)&0&0&0&0&0\\\
0&0&C_{\text{df},\,y_{1}}(s)&0&0&0&-\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\\\
0&0&0&C_{\text{df},\,y_{2}}(s)&0&0&\frac{L}{2}\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\\\
0&0&-\frac{1}{L}\,C_{\text{df},\,y_{1}}(s)&\frac{1}{L}\,C_{\text{df},\,y_{2}}(s)&C_{\text{sus},\,\phi_{1}}(s)&0&0\\\
0&0&-\frac{1}{L}\,C_{\text{df},\,y_{1}}(s)&\frac{1}{L}\,C_{\text{df},\,y_{2}}(s)&0&C_{\text{sus},\,\phi_{2}}(s)&0\\\
0&0&-\frac{1}{L}\,C_{\text{df},\,y_{1}}(s)&\frac{1}{L}\,C_{\text{df},\,y_{2}}(s)&0&0&0\end{pmatrix}$}\leavevmode\nobreak\
,$ (3.51)

where possible actuation gains can be intended as multiplicative factor of
each single control law. The first-order perturbation due to the control
actuation cross-talk is given by

$\delta\bm{A}\,\bm{C}_{0}=\frac{1}{2}\scalebox{0.7}{$\begin{pmatrix}0&0&\delta_{\text{df}}^{-}\,C_{\text{df},\,y_{1}}(s)&\delta_{\text{df}}^{+}\,C_{\text{df},\,y_{2}}(s)&-\delta_{\text{sus},\,\phi_{1}}\,l\,C_{\text{sus},\,\phi_{1}}(s)&0&-\delta_{\text{sus},\,y_{1}}\,L\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\\\
0&0&0&0&\delta_{\text{sus},\,\phi_{1}}\,l\,C_{\text{sus},\,\phi_{1}}(s)&-\delta_{\text{sus},\,\phi_{2}}\,l\,C_{\text{sus},\,\phi_{2}}(s)&\delta_{\text{sus}}^{+}\,L\,C_{\text{sus},\,\phi_{\text{SC}}}(s)\\\
0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0\end{pmatrix}$}\leavevmode\nobreak\ ,$ (3.52)

where
$\delta_{\text{df}}^{-}=\left(\delta_{\text{df},\,y_{\text{SC}}}-\delta_{\text{df},\,\phi_{\text{SC}}}\right)$,
$\delta_{\text{df}}^{+}=\left(\delta_{\text{df},\,y_{\text{SC}}}+\delta_{\text{df},\,\phi_{\text{SC}}}\right)$
and
$\delta_{\text{sus}}^{+}=\delta_{\text{sus},\,y_{1}}+\delta_{\text{sus},\,y_{2}}$
are three new definitions of effective cross-talk coefficients.

The imperfection matrices $\delta\bm{D}$, $\delta\bm{A}\,\bm{C}_{0}$ and
$\delta\bm{S}$, together with the nominal matrices $\bm{D}_{0}$, $\bm{C}_{0}$
and $\bm{S}_{0}=\bm{1}$, allows for a simplification of the nominal
differential operator and its imperfection in (3.29)

$\displaystyle\bm{\Delta}_{0}$
$\displaystyle=\bm{D}_{0}+\bm{C}_{0}\leavevmode\nobreak\ ,$ (3.53a)
$\displaystyle\delta\bm{\Delta}$
$\displaystyle=\delta\bm{D}-\bm{D}_{0}\,\delta\bm{S}+\delta\bm{A}\,\bm{C}_{0}\leavevmode\nobreak\
.$ (3.53b)

The second equation explains the fact that the imperfection to the
differential operator converting sensed coordinates into total equivalent
acceleration is given by three terms: the dynamical cross-talk, the sensing
cross-talk and the control actuation cross-talk. The application of the first
matrix to the first-order correction to the sensed coordinates finally gives
the various cross-talk contributions in the equation of motion (3.38).
Inverting $\bm{\Delta}_{0}$, and applying it to the various cross-talk
contributions of (3.38), it allows for a modeling of the response of the
system along the optical axis to noise sources affecting the nominally
orthogonal degrees of freedom.

##  4 System identification

This chapter focuses on a topic that can be considered the core of the whole
LPF mission in view of characterizing the total equivalent acceleration noise
affecting each single LISA arm. In system identification LPF is modeled as a
matrix of parametric transfer functions. Targeted experiments where the system
is stimulated on each degree of freedom can be used to infer the values of the
critical parameters contained in those functions.

The preceding chapter described the closed-loop dynamics underlying LPF, the
methods to handle and subtract the applied control forces, the sensing and the
dynamical couplings between the TMs and the SC, the extent to which system
transients can be suppressed and the estimation of the equivalent out-of-loop
acceleration noise can be made possible.

This chapter shows an application of the ideas in a mission-like fashion with
numerical applications. It assumes a model for LPF along the optical axis,
which gives the description of the dynamics to first approximation. The aim is
to simulate and analyze the data as they will be released during the mission.
To simplify the discussion, only two experiments are considered, allowing for
a complete identification of the system along the optical axis. As the methods
developed in this chapter are general, they can also be applied to the study
of more sophisticated experiments. Examples are the cross-talk experiments
from orthogonal degrees of freedom to the optical axis: the modeled transfer
functions are different, the dimensionality of the system is different, but
the approach is the same. In the end, all experiments analyzed with the
methods described in this chapter will hopefully provide a coherent
understanding of the system, contributing to the final success of the LPF
mission.

In turn, this chapter discusses: the dynamical model assumed for simulations
and analysis; the noise characterization of the system; the simulated
identification experiments; the parameter estimation method, the validation
and the robustness to non-standard scenarios. Finally, it demonstrates the
impact of system identification on the estimation of the residual equivalent
acceleration noise and the suppression of transients in data produced by a
simulator provided by ESA.

###  4.1 Dynamical model

Section 3.3 provided a model of LPF along $x$, the optical axis. In the main
science mode, the reference TM is in free fall along $x$. The other TM and the
SC are, respectively, forced by capacitive and thruster actuation to follow
the reference TM along the same axis. The interferometer keeps track of the
relative motion between the reference TM and the SC ($o_{1}$) and between the
two TMs ($o_{12}$). The two readouts expressed in displacement are fed into
the DFACS to command force actuation, hence minimizing the relative motion.

Figure 4.1 and Figure 4.2 show the frequency dependence of the two control
laws converting the sensed displacements to commanded forces to the thrusters
(drag-free loop) and the electrostatic suspensions (electrostatic suspension
loop), respectively. At low frequency, the drag-free gain is very high due to
the need for suppressing the SC jitter that dominates $o_{1}$. Instead, the
electrostatic suspension gain is designed to suppress the force couplings
between the TMs and the SC that dominate $o_{12}$. The control laws used in
this thesis are provided by ASTRIUM [35] – the main industry contractor of
LPF.

Figure 4.1: Frequency dependence of the drag-free loop controller per unit SC
mass. The very high gain at low frequency is explained by the need for
removing the thruster noise. Following (3.9), $1\,\mathrm{\mu m}$ sensed
displacement of the first TM relative to the SC produces a thruster actuation
of $\raise 0.58554pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}0.02\,\mathrm{\mu N}$
at $1\,\mathrm{mHz}$. Figure 4.2: Frequency dependence of the electrostatic
suspension loop controller per unit TM mass. Here there is no such a huge
variation in the order of magnitude as for the drag-free controller in Figure
4.1. The control law is designed in particular to suppress the force couplings
between the TMs and the SC at low frequency. Following (3.9), $1\,\mathrm{\mu
m}$ sensed displacement of the second TM relative to the first one produces a
capacitive actuation of $\raise 0.58554pt\hbox{$$ \mbox{\scriptsize$\sim$
}$$}0.1\,\mathrm{nN}$ at $1\,\mathrm{mHz}$.

As described in Section 3.3, the system can be modeled by the operators
$\bm{D}$ (dynamics), $\bm{S}$ (sensing) and $\bm{A}$ (actuation) representing
different non-idealities in the practical implementation of the closed-loop
LISA arm. The operators contain all system parameters describing the dynamics
along the optical axis. One last source of indetermination introduced here is
a delay in the application of the guidance signals

$\bm{T}=\begin{pmatrix}e^{-s\,\Delta t_{1}}&0\\\ 0&e^{-s\,\Delta
t_{2}}\end{pmatrix}\leavevmode\nobreak\ ,$ (4.1)

whose possible causes may be either due to the digitalization of the
continuous control laws or to bus delays, a possibility not considered in a
previous model [47]. With the introduction of the delays, the model (3.8b)
becomes now

$\bm{T}_{o_{\text{i}}\rightarrow
o}=\bm{\Delta}^{-1}\bm{A}\,\bm{C}\,\bm{T}\leavevmode\nobreak\ ,$ (4.2)

where the differential operator $\bm{\Delta}$, defined in (3.8a), converts the
sensed motion into total equivalent acceleration.

Figure 4.3 shows the transfer gains of the model
$\bm{T}_{o_{\text{i}}\rightarrow o}$, whereas the dynamical cross-talk from
the differential channel to the first one is definitely negligible with peak
gain of about $4\\!\times\\!10^{-6}$ at $30\,\mathrm{mHz}$. The diagonal
elements have respectively peak gains of almost 3 at $0.1\,\mathrm{Hz}$ and
about 2 at $0.8\,\mathrm{mHz}$. The dynamical cross-talk from the first
channel to the differential one has peak gain of about $5\\!\times\\!10^{-2}$
at $0.5\,\mathrm{mHz}$. The above transfer matrix is used to both model the
outputs of the system subjected to bias injections and perform system
identification.

Figure 4.3: Frequency dependence of the transfer matrix
$\bm{T}_{o_{\text{i}}\rightarrow o}$ used for system identification. The
transfer function $T^{11}_{o_{\text{i}}\rightarrow
o}=T_{o_{\text{i,1}}\rightarrow o_{1}}$ has peak gain of almost 3 at
$0.1\,\mathrm{Hz}$. The transfer function $T^{22}_{o_{\text{i}}\rightarrow
o}=T_{o_{\text{i,12}}\rightarrow o_{12}}$ has peak gain of about 2 at
$0.8\,\mathrm{mHz}$, then it quickly decays. The dynamical cross-talk
$T^{12}_{o_{\text{i}}\rightarrow o}=T_{o_{\text{i,1}}\rightarrow o_{12}}$ has
peak gain of about $5\\!\times\\!10^{-2}$ at $0.5\,\mathrm{mHz}$. The other
dynamical cross-talk is negligible since has peak gain of about
$4\\!\times\\!10^{-6}$ at $30\,\mathrm{mHz}$.

Throughout this chapter bias injections at the level of controller guidance
signals $\bm{o}_{\text{i}}$ 111Following (3.9), a bias in the guidance signals
is equivalent to a commanded force bias directly applied onto the SC through
thruster actuation and the second TM through capacitive actuation. are
considered and the transfer matrix in (4.2) models the response of the system
to those signals. As the modeled system parameters appear in the operators,
$\bm{T}_{o_{\text{i}}\rightarrow o}$ is parameter-dependent. The modeled
system response is then parameter-dependent. The parameters can be arranged in
a vector that will be abstractly referred to $\bm{p}$

$\bm{p}=\scalebox{0.8}{$\begin{pmatrix}\omega_{1}^{2}\\\ \omega_{12}^{2}\\\
S_{21}\\\ A_{\text{df}}\\\ A_{\text{sus}}\\\ \Delta t_{1}\\\ \Delta t_{2}\\\
\end{pmatrix}$}\leavevmode\nobreak\ ,$ (4.3)

where Table 4.1 provides a description of the above system parameter with
initial plausible estimates coming from on-ground measurements and theoretical
modeling.

Table 4.1: List of the modeled system parameters, introduced in Section 3.3,
except for $\Delta t_{1}$ and $\Delta t_{2}$, with descriptions and initial
estimates. The parameters that are fitted to data are $\omega_{1}^{2}$,
$\omega_{12}^{2}$, $S_{21}$, $A_{\text{df}}$, $A_{\text{sus}}$, $\Delta
t_{1}$, $\Delta t_{2}$.

Parameter |  Description |  Note |  Estimate
---|---|---|---
$\omega_{1}^{2}$, $\omega_{12}^{2}$ |  parasitic stiffness constants modeling residual oscillator-like couplings between the SC and the reference TM and between the two TMs |  must be estimated from experiments |  $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}-1\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$
$S_{21}$ |  sensing cross-talk between $o_{1}$ and $o_{12}$ interferometric readouts |  must be estimated from experiments |  $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}1\\!\times\\!10^{-4}$
$A_{\text{df}}$, $A_{\text{sus}}$ |  actuation gains for the application of forces by the thrusters and the electrostatic suspensions |  must be estimated from experiments |  $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}1$
$\Delta t_{1}$, $\Delta t_{2}$ |  delays in the application of biases to the controller computing the actuation |  must be estimated from experiments |  $\lesssim 1\,\mathrm{s}$
$\Gamma_{x}$ |  gravity gradient between the two TMs |  could be estimated from experiments with different actuation stiffness but difficult, considered fixed |  $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}4\\!\times\\!10^{-9}\,\mathrm{s^{-2}}$
$m_{1}$, $m_{2}$, $m_{\text{SC}}$ |  masses of TMs and SC |  considered fixed |  $1.96\,\mathrm{kg}$, $422.7\,\mathrm{kg}$

The aim of system identification, as thoroughly described in this chapter, is
the estimation of these system parameters with targeted experiments.

####  4.1.1 Anelasticity and damping

The parameters defined above are implicitly assumed to be independent from
frequency. For example, the parasitic stiffness constant may show a frequency
dependence due to anelasticity (an “internal” dissipation of the string
constant) or a damping effect.

For the sake of clarity, $\omega_{0}^{2}$ is the (negative) stiffness constant
not to be confused with the Fourier angular frequency $\omega$. An
anelasticity can be modeled as a frequency dependence in the imaginary part
[48] of a complex stiffness constant

$\tilde{\omega}^{2}(\omega)=\omega_{0}^{2}\left[1+i\phi(\omega)\right]\leavevmode\nobreak\
,$ (4.4)

where $\phi$ is named the loss angle modeling the dissipation. Sources of
dissipation are the dielectric losses in the surface of the electrodes facing
the TM that can be modeled by a constant

$\phi_{\epsilon}(\omega)=-\delta_{\epsilon}\leavevmode\nobreak\ ,$ (4.5)

such that it produces a force proportional to displacement and in phase with
velocity 222Thanks to the imaginary unit. The minus sign is due to the fact
that the stiffness constant is usually negative.. The other source of
dissipation is the residual gas damping that can be modeled by a function
proportional to frequency

$\phi_{\text{g}}(\omega)=\frac{\omega}{\omega_{0}^{2}\tau}\leavevmode\nobreak\
,$ (4.6)

where $\tau$ is the damping characteristic time. The above produces a force
proportional to velocity 333In fact, the damped harmonic oscillator in
frequency domain is $(-\omega^{2}+i\gamma\omega+\omega_{0}^{2})x=f$, where
$\gamma=1/\tau$ is the damping coefficient. Then, the complex stiffness
constant is given by
$\tilde{\omega}^{2}=\omega_{0}^{2}+i\gamma\omega=\omega_{0}^{2}+i\omega/\tau$..

The loss angle function may show other interesting features beyond the ones
reported here. To first approximation, the following analysis assumes that all
parameters are independent from frequency, at least within the frequency band
of interest.

###  4.2 Noise characterization

One of the objectives of the LPF mission is to provide a full noise projection
of the total equivalent differential acceleration noise between the TMs. As
this is well beyond the scope of this thesis, the following presents a hint of
the problem. Moreover, a theoretical projection of the observed displacement
noise is needed in advance in order to identify the dominant effects in the
noise and produce the generating filters used for all simulations. The noise
projections shown in this section are given by plausible noise shapes
implemented in the simulator provided by ESA (that will be specifically
introduced in the first paragraph of Section 4.4.5).

Figure 4.4 shows a theoretical noise projection of the equivalent acceleration
noise affecting the $x_{1}$ degree of freedom. Evidently, the thruster
actuation noise dominates the total noise budget in the frequency band of
interest. Other important noise sources are the infrared thermal emission of
the SC external surface and the $o_{1}$ sensing noise.

Figure 4.4: Theoretical noise projection of the residual equivalent
acceleration noise of the relative motion between the SC and the first TM for
the nominal dynamics along $x$. The thruster actuation noise dominates the
total noise budget (dashed line) in the frequency band of interest. Other
important noise sources are the infrared thermal emission of the SC external
surface and the $o_{1}$ sensing noise.

Figure 4.5 shows the second and most important projection of the equivalent
differential acceleration noise affecting the $x_{12}$ degree of freedom. A
turning point around $6\,\mathrm{mHz}$ between two regimes is clearly evident.
At high frequency, the $o_{1}$ sensing noise dominates the total noise budget.
At low frequency, $\nicefrac{{2}}{{3}}$ of the total noise budget (in units of
$\sqrt{\text{PSD}}$) is due to force couplings between the SC and the TMs.
Other important noise sources, intervening at low frequency, are the
capacitive actuation noise on the second TM, forces on the TMs coming from
outside the SC and the $o_{12}$ and $o_{1}$ sensing noises.

Figure 4.5: Theoretical noise projection of the residual equivalent
acceleration noise of the relative motion between the TMs for the nominal
dynamics along $x$. At high frequency, the $o_{12}$ sensing noise dominates
the total noise budget (dashed line). At low frequency, $\nicefrac{{2}}{{3}}$
of the total noise budget is due to force couplings between the SC and the
TMs. Other important noise sources are the capacitive actuation noise on the
second TM, forces on the TMs coming from outside the SC and the $o_{12}$ and
$o_{1}$ sensing noises.

The above acceleration noise projections are the equivalent acceleration
inputs to LPF coming from reasonable noise shapes, producing a characteristic
output in the interferometric readouts. Figure 4.6 and Figure 4.7 contain the
relative projections for the two interferometric readouts, $o_{1}$ and
$o_{12}$ along $x$, produced with a plausible transfer model. Analogously to
the equivalent acceleration noise, for $o_{1}$ the thruster actuation noise
dominates the total noise budget in the frequency band of interest. As
previously pointed out, for $o_{12}$ there is a turning point around
$6\,\mathrm{mHz}$. At high frequency, the $o_{1}$ sensing noise dominates the
total noise budget. At low frequency, $\nicefrac{{2}}{{3}}$ of the total noise
budget is due to force couplings between the SC and the TMs. Secondary
sources, intervening at low frequency, are the capacitive actuation noise on
the second TM, forces on the TMs coming from outside the SC, the $o_{12}$ and
$o_{1}$ sensing noises and the thruster actuation noise.

Figure 4.6: Theoretical noise projection of the $o_{1}$ data channel sensing
the relative motion between the SC and the first TM for the nominal dynamics
along $x$. The thruster actuation noise dominates the total noise budget
(dashed line) in the frequency band of interest. Secondary sources are the
infrared thermal emission of the SC external surface and the $o_{1}$ sensing
noise. Figure 4.7: Theoretical noise projection of the $o_{12}$ data channel
sensing the relative motion between the TMs for the nominal dynamics along
$x$. At high frequency, the $o_{12}$ sensing noise dominates the total noise
budget (dashed line). At low frequency, $\nicefrac{{2}}{{3}}$ of the total
noise budget is due to force couplings between the SC and the TMs. Secondary
sources are the capacitive actuation noise on the second TM, forces on the TMs
coming from outside the SC, the $o_{12}$, $o_{1}$ sensing noises and the
thruster actuation noise.

The noise shapes of the interferometric readouts (with their cross-
correlation) are also used for simulation purposes. From those models, noise
shaping filters are derived and integrated into a multi-channel cross-
correlated noise generator [49]. Figure 4.8 reports an example of a noise run
lasting 12 hours and obtained by coloring an input zero-mean
$\delta$-correlated (white) Gaussian noise with those filters. $o_{12}$ shows
a huge red component caused by the increase of the PSD at low frequency, due
to forces on the TMs, as predicted by Figure 4.7. While $o_{1}$ is dominated
by the thruster jitter, $o_{12}$ becomes much less noisy at high frequency,
being dominated by readout noise only. The red noise shape of $o_{12}$ is an
expected feature during the experiments of the LPF mission.

Figure 4.8: A simulated noise run of about 12 hours. $o_{1}$ and $o_{12}$ are
the two interferometer readings. Notice the behavior of $o_{12}$ at low
frequency – an expected feature during the LPF mission – showing a huge red
component caused by force couplings between the TMs and the SC. At high
frequency, $o_{12}$ becomes much less noisy than $o_{1}$, the former being
dominated by only interferometer readout noise and the latter by thruster
noise.

###  4.3 Identification experiments

Among the series of experiments characterizing the LPF mission, a few of
capital importance will tackle system identification. This thesis considers
two identification experiments allowing for a complete identification of the 7
most important system parameters introduced in Section 4.1. As said,
considering bias injections at the level of controller guidance signals is
completely equivalent to applying direct force stimuli through the equivalence
given by (3.5). In the nominal $x$-dynamics two experiments are defined:

  1. 1.

an injection into the controller guidance of the $o_{1}$ channel, namely
$o_{\text{i},1}$, producing forces on the SC through thruster actuation;

  2. 2.

an injection into the controller guidance of the $o_{12}$ channel, namely
$o_{\text{i},12}$, producing forces on the second TM through capacitive
actuation.

To naively understand how the parameters can be determined from the above
experiments and the model described in Section 4.1, it is useful to make a
projection of the differential operator, whose inverse enters into the
transfer matrix through (4.2).

Figure 4.9 contains the projection of the differential operator from the first
channel to equivalent acceleration in terms of: (i) dynamics and sensing; (ii)
control. Clearly, the control dominates the transfer for almost the entire
frequency band, in order to attenuate the SC jitter. For this reason,
injecting a signal into the first controller guidance (i.e., applying a
thruster actuation on the SC) allows for the identification, in turn, of: the
actuation gain, $A_{\text{df}}$, the first TM coupling to the SC,
$\omega_{1}^{2}$, as well as a possible delay in the application of the same
bias, $\Delta t_{1}$.

Figure 4.9: Frequency dependence of the differential operator for the transfer
from $o_{1}$ to equivalent acceleration. The control dominates the transfer
for almost the entire frequency band, in order to attenuate the SC jitter.

Analogously, Figure 4.10 contains the projection of the differential operator
from the differential channel to equivalent acceleration. Below
$1\,\mathrm{mHz}$, the control dominates the transfer in order to compensate
the differential force disturbances. Above $1\,\mathrm{mHz}$, dynamics and
sensing dominate the transfer. For this reason, injecting a signal into the
second controller guidance (i.e., applying a capacitive actuation on the
second TM) allows for the identification, in turn, of: the actuation gain,
$A_{\text{sus}}$, the differential coupling between the TMs,
$\omega_{12}^{2}$, as well as a possible delay in the application of the same
bias, $\Delta t_{2}$. Given the cross-talk elucidated in Figure 4.3 at low
frequency, the sensing cross-talk, $S_{21}$, can also be determined.

Figure 4.10: Frequency dependence of the differential operator for the
transfer from $o_{12}$ to equivalent acceleration. The control dominates the
transfer at low frequency, in order to compensate the differential force
disturbances.

To conclude the discussion on the projection of the differential operator, it
is worth noting that the off-diagonal terms contribute with a figure of
$1\\!\times\\!10^{-7}\,\mathrm{s^{-2}}$. In particular, even if not shown in
any figures, the transfer from the first channel to the equivalent
differential acceleration is dominated by dynamics and sensing; the other by
control at low frequency and dynamics and sensing at high frequency. As the SC
motion is common-mode and the first and differential channel are correlated,
the estimation of the differential acceleration noise can not be performed
independently of the first channel, which is the only means by which the SC
jitter can be measured and subtracted. The details of such an estimation will
be given in Section 4.5.

The next section is devoted to the estimation of the 7 system parameters by
means of a MIMO approach that maximizes the overall information. The
identification experiments defined at the beginning of this section are
simulated for a total duration of almost 3 hours each – a suitable timescale
for the mission – by injecting stimulating biases. The following facts are
assumed:

  1. 1.

the noise $\bm{o}_{\text{n}}$ is generated as in Section 4.2, independently
from the noise-only run which is used for noise characterization, and is
Gaussian and stationary;

  2. 2.

the signals $\bm{o}_{\text{s}}$ are simulated in time domain with a MIMO
approach by means of (3.8b), i.e., by anti-Fourier transforming 444The
numerical implementation of the direct and inverse Fourier transform are the
Fast Fourier Transform (FFT) and Inverse FFT (IFFT). Being circular
operations, the input time-series needs to be zero-padded to avoid systematic
errors caused by wrapped-around data [50]. A conservative default value of one
data length is assumed. with $\mathcal{F}^{-1}$ the deterministic input
signals

$\bm{o}_{\text{s}}(t,\bm{p}_{\text{true}})=\mathcal{F}^{-1}\left[\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p}_{\text{true}})\,\bm{o}_{\text{i}}(\omega)\right](t)\leavevmode\nobreak\
,$ (4.7)

where $\bm{p}_{\text{true}}$ is the set of assumed true system parameter
values to be estimated from the analysis and which the estimation of residual
equivalent acceleration noise depends on;

  3. 3.

the superposition principle of signals and noise holds true in the hypothesis
of small motion and in absence of non-linearities in the system, so that the
“experimental” data are simulated by

$\bm{o}_{\text{exp}}=\bm{o}_{\text{s}}+\bm{o}_{\text{n}}\leavevmode\nobreak\
.$ (4.8)

The underlying idea in parameter estimation is to excite the system with
proper high SNR signals so that the modeled parameters can be measured. A
typical injected bias is a series of sine waves of logarithmically increasing
frequency, with integer number of cycles, divided by gaps of $150\,\mathrm{s}$
to allow for system relaxation. The sine stretches last $1200\,\mathrm{s}$
each. The amplitudes are conservatively selected not to exceed $1\%$ of the
operating range of the interferometer, corresponding to a maximum sensed
displacement of $1\,\mathrm{\mu m}$, and $10\%$ of the maximum allowed force
authority, corresponding to $10\,\mathrm{\mu N}$ of thruster actuation and
$0.25\,\mathrm{nN}$ of capacitive actuation. The biases are parameterized in
Table 4.2 and referred to the standard input signals used for the rest of the
analysis. Instead, Chapter 5 will focus on the optimization of the same input
signals.

Table 4.2: Controller guidance signals injected as biases for system identification. The sine stretches last $1200\,\mathrm{s}$ each and are separated by gaps of $150\,\mathrm{s}$. The sine waves perform an integer number of cycles, from 1 to 64. The amplitudes are selected to not exceed $1\%$ of the operating range of the interferometer and $10\%$ of the maximum force authority. $o_{\text{i},1}$ for Exp. 1 |  $o_{\text{i},12}$ for Exp. 2
---|---
$f\,\mathrm{[mHz]}$ | $a\,\mathrm{[\mu m]}$ | $f\,\mathrm{[mHz]}$ | $a\,\mathrm{[\mu m]}$
$0.83$ | $1.0$ | $0.83$ | $0.80$
$1.7$ | $1.0$ | $1.7$ | $0.48$
$3.3$ | $1.0$ | $3.3$ | $0.19$
$6.6$ | $1.0$ | $6.6$ | $0.088$
$13$ | $0.59$ | $13$ | $0.096$
$27$ | $0.28$ | $27$ | $0.18$
$53$ | $0.14$ | $53$ | $0.46$

Data are simulated at $10\,\mathrm{Hz}$ and decimated to $1\,\mathrm{Hz}$ to
ease data processing. During the mission, data will be collected at a sample
rate between $1$ and $10\,\mathrm{Hz}$, depending on the experiment and
available down-link bandwidth. The simulation of the first experiment, with
injection of the $o_{i,1}$ signal of Table 4.2, is shown in Figure 4.11. The
response of the system in $o_{1}$ is approximately equal to $o_{i,1}$, except
at high frequency where there is a modest gain due to the particular shape of
the first diagonal element of the transfer function at that frequency. A
residual signal in $o_{12}$ of absolute peak $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}4\\!\times\\!10^{-8}\,\mathrm{m}$ is also visible
and due to dynamical cross-talk. As said before, the gaps allow for system
relaxation, particularly at high frequency.

Figure 4.11: Exp. 1 synthetic data. An injection of sine-wave signals lasting
for almost 3 hours into the first controller guidance $o_{\text{i},1}$
produces a different response in the two interferometer readings. The response
in $o_{1}$ is approximately equal to $o_{\text{i},1}$ (dashed line), except at
high frequency where there is a modest gain. A residual signal in $o_{12}$ of
absolute peak $\raise 0.58554pt\hbox{$$ \mbox{\scriptsize$\sim$
}$$}4\\!\times\\!10^{-8}\,\mathrm{m}$ is due to dynamical cross-talk (see
inset at the left bottom side). Gaps between two cycles of injection allow for
system relaxation (see inset at the right top side).

The simulation of the second experiment, with injection of the $o_{i,12}$
signal of Table 4.2, is shown in Figure 4.12. The response of the system in
$o_{12}$ is evidently phase delayed to $o_{\text{i},12}$. At high frequency,
the very low gain of the transfer function almost suppresses the signal. Since
the transfer from $o_{i,12}$ to $o_{1}$ is negligible, in this experiment
$o_{1}$ has signal contribution completely hidden by noise. For this reason,
during the mission the $o_{1}$ readout will serve as a useful sanity check for
a first understanding of the model.

Figure 4.12: Exp. 2 synthetic data. An injection of sine-wave signals lasting
for almost 3 hours into the second controller guidance $o_{\text{i},12}$
produces a different response in the two interferometer readings. The response
in $o_{12}$ is evidently phase delayed to $o_{\text{i},12}$. At high
frequency, the very low gain of the transfer function almost suppresses the
signal (see inset). The $o_{1}$ data channel has negligible contribution
hidden by the noise.

###  4.4 Parameter estimation

During the mission, noise runs will be used to characterize the noise itself
and estimate the total equivalent input acceleration. The estimation of the
total equivalent acceleration is possible if LPF is properly modeled. For this
in the various experiments, signals will be injected along different degrees
of freedom to study the response of the system. Along $x$, LPF will be
characterized giving, as a first approximation, the nominal dynamics. Instead,
along others degrees of freedom, LPF will be characterized in terms of the
many cross-talk contributions arising from the dynamical couplings, the
imperfections in the sensing conversion and the imperfections in the
actuation.

This section handles the general problem of estimating the LPF parameters
modeled as a MIMO dynamical system, where different inputs enters into the
system and produce a response in different outputs. For the sake of
simplicity, for the rest only the two experiments introduced above – the
characterization of the nominal dynamics along the optical axis – are
considered, bearing in mind that the method is general enough to handle more
complex experiments. An example would be the identification of the $xy$ cross-
talk, in where guidance or force bias signals are injected, in turn, along
$y_{1}$, $y_{2}$, $\phi_{1}$, $\phi_{2}$ and $\phi_{\text{SC}}$ to study the
response along the optical axis.

Finally, this section develops and validates the estimation procedures on the
two most important experiments described in the previous section. It also
shows the application to a couple of non-standard scenarios that may happen
during the real LPF mission.

####  4.4.1 Review of the problem

The experimental data (either simulated or from the mission) can be modeled
superimposing deterministic signals with noise

$\bm{o}_{\text{exp}}=\bm{o}_{\text{s}}+\bm{o}_{\text{n}}\leavevmode\nobreak\
,$ (4.9)

where $\bm{o}_{\text{n}}$ is the output noise with cross PSD matrix
$\bm{S}_{\text{n}}$ and

$\bm{o}_{\text{s}}(t,\bm{p})=\mathcal{F}^{-1}\left[\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p})\,\bm{o}_{\text{i}}(\omega)\right](t)\leavevmode\nobreak\ ,$
(4.10)

are the so-called template signals obtained by injecting bias guidance signals
$\bm{o}_{\text{i}}$ into the system modeled by the transfer matrix
$\bm{T}_{o_{\text{i}}\rightarrow o}$.

It is useful to think that the experimental data depends on the true parameter
values

$\bm{o}_{\text{exp}}=\bm{o}_{\text{exp}}(t,\bm{p}_{\text{true}})\leavevmode\nobreak\
,$ (4.11)

that need to be estimated from fitting procedures. In the case of simulated
experiments, the true values are exactly those used in data generation. In the
case of real mission experiments, the true values are actually those giving
the best possible description of the data, the one that perfectly subtracts
the deterministic signals, hence recovering the instrumental noise shapes.

In the same way, the observed noise (either simulated or from the mission)
depends on the parameter values

$\bm{o}_{\text{n}}=\bm{o}_{\text{n}}(t,\bm{p}_{\text{true}})\leavevmode\nobreak\
;$ (4.12)

but can be considered constant with respect to the parameter values for the
timescale of an identification experiment where only high SNR signals will be
injected.

The scope of parameter estimation is to recover the best possible description
of the experimental data. If the residuals between the experimental data and
the modeled template signals are defined by

$\bm{o}_{\text{r}}=\bm{o}_{\text{exp}}-\bm{o}_{\text{s}}\leavevmode\nobreak\
,$ (4.13)

the best possible description of the experimental data is given by

$\bm{o}_{\text{r}}(t,\bm{p}_{\text{est}})\simeq\bm{o}_{\text{n}}(t,\bm{p}_{\text{true}})\leavevmode\nobreak\
,$ (4.14)

implying that the residuals evaluated at the estimated parameter values
$\bm{p}_{\text{est}}$ recover the true instrumental noise.

####  4.4.2 Estimation method

LPF is a MIMO dynamical system for which each experiment has a unique set of
meaningful parameters. Hence, for two generic experiments two sets of
parameters can be independently determined. Sometimes a subset may be shared
between the two; sometimes there could be parameters that can be estimated by
only a particular experiment. Moreover, each experiment has multiple readouts
sensitive to different parameters. Section 4.3 has given an intuitive hint of
such an idea.

The first approach is to build an information-weighted average [40, 47] of
different parameter estimates coming from all readouts and experiments. If
$\bm{p}_{ij}$ are the parameter estimates of the $i$-th experiment and $j$-th
readout, the corresponding Fisher information matrix [51]

$\bm{\mathcal{I}}_{ij}=\int{\nabla_{\bm{p}}\bm{o}_{\text{r}}^{(ij)}(\omega,\bm{p}_{\text{est}})}^{*}\,\bm{S}_{\text{n}}^{(ij)}(\omega)^{-1}\,\nabla_{\bm{p}}\bm{o}_{\text{r}}^{(ij)}(\omega,\bm{p}_{\text{est}})\,\text{d}\omega\leavevmode\nobreak\
,$ (4.15)

where $\bm{S}_{\text{n}}^{(ij)}$ is the noise PSD of $i$-th experiment and
$j$-th readout, $\bm{o}_{\text{r}}^{(ij)}$ is the corresponding vector of
residuals, $\nabla_{\bm{p}}$ is the gradient with respect to the parameters
and ∗ is the conjugate transpose. The final combined parameter estimates are
given by

$\bm{p}=\bm{\mathcal{I}}^{-1}\sum_{i=1}^{N_{\text{exp}}}\sum_{j=1}^{N_{o}}\bm{\mathcal{I}}_{ij}\,\bm{p}_{ij}\leavevmode\nobreak\
,$ (4.16)

where $N_{\text{exp}}$ is the number of experiments and $N_{o}$ the number of
readouts per experiment assumed the same across the experiments. The combined
Fisher information matrix is

$\bm{\mathcal{I}}=\sum_{i=1}^{N_{\text{exp}}}\sum_{j=1}^{N_{o}}\bm{\mathcal{I}}_{ij}\leavevmode\nobreak\
.$ (4.17)

Notice that the estimates $\bm{p}_{ij}$ may have different dimension depending
of the $i$-th experiment and $j$-th interferometric readout; the same happens
for the corresponding information matrices. The issue can be easily solved by
inserting zeros where there is no information.

An example can readily show that the definition of (4.16) is not robust. In
fact, suppose that the estimation of the system parameters is performed
independently on each readout and one of those parameters has a biased value
for an inaccuracy of the transfer matrix model. Therefore, the information
matrix for that estimate is biased and the combined one in (4.17) as well. The
numerical inversion in (4.16) inexorably amplifies that bias to the combined
parameter estimates. To overcome the problem, one could try removing the
failing estimates (which is possible only if one has good indication of what
the real values are, for example, from ground measurements or previous
independent experiments), but in doing so information and precision would
definitely be lost.

The only solution is to attack the problem by a complete MIMO approach where
the poor information coming from the biased model of a readout is continuously
compensated by the others as the optimization goes on. One other advantage is
that a joint information can likely remove or, at least, reduce the effect of
parameter degeneracies.

The MIMO-Multi-Experiment joint log-likelihood of the system is a
generalization of the standard definition [51] and is given by

$\chi^{2}(\bm{p})=\int{\bm{o}_{\text{r}}(\omega,\bm{p})}^{*}\,\bm{S}_{\text{n}}(\omega)^{-1}\,\bm{o}_{\text{r}}(\omega,\bm{p})\,\text{d}\omega\leavevmode\nobreak\
,$ (4.18)

where

$\bm{o}_{\text{r}}(\omega,\bm{p})=\bm{o}_{\text{exp}}(\omega)-\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p})\,\bm{o}_{\text{i}}(\omega)\leavevmode\nobreak\ ,$ (4.19)

are the residuals between the experimental data $\bm{o}_{\text{exp}}$ and the
modeled system response. $\bm{o}_{\text{i}}$ are the controller biases,
$\bm{T}_{o_{\text{i}}\rightarrow o}$ the transfer matrix depending on all
system parameters $\bm{p}$ (stiffness constants, sensing cross-talk, etc.),
$\bm{S}_{\text{n}}$ the cross output noise PSD matrix assumed constant to the
system parameters. For two experiments and two interferometric readouts each,
$\bm{o}_{\text{i}}$ is a 4-vector, null in the second and third element, since
the injection is in $o_{\text{i},1}$ (first experiment) and $o_{\text{i},12}$
(second experiment); $\bm{T}_{o_{\text{i}}\rightarrow o}$ is a block diagonal
$4\times 4$-matrix replicating the same $2\times 2$-matrix;
$\bm{S}_{\text{n}}$ is a $4\times 4$-matrix of cross PSDs between different
readouts and experiments; $\bm{o}_{\text{exp}}$ is a 4-vector of all
experimental readouts.

Assuming that all readouts are sampled at the same rate and last for the same
duration, the overall number $\nu$ of degrees of freedom for the problem is
defined as

$\nu=N_{\text{exp}}\times N_{o}\times
N_{\text{data}}-N_{p}\leavevmode\nobreak\ ,$ (4.20)

where $N_{\text{exp}}$ is the number of experiments; $N_{o}$ is the number of
readouts per experiment (assumed the same across the experiments);
$N_{\text{data}}$ is the number of data points per readout; $N_{p}$ is the
dimension of the parameter space. For example, $\nu\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}4\\!\times\\!10^{4}$ for two experiments, two
readouts each, lasting for about 3 hours and sampled at $1\,\mathrm{Hz}$. For
the rest, if not otherwise stated, the reduced log-likelihood $\chi^{2}/\nu$
will be used in place of the standard definition, as its expectation value is
1.

Notice that system identification may be also implemented in the domain of
equivalent acceleration. If the $\bm{\Delta}$ operator is invertible, the two
approaches – identification in acceleration and displacement – are completely
equivalent. In fact,

$\begin{split}\chi^{2}&=\int{\left(\bm{f}_{\text{exp}}-\bm{f}_{\text{mdl}}\right)}^{*}\,\bm{S}_{\text{n},f}^{-1}\left(\bm{f}_{\text{exp}}-\bm{f}_{\text{mdl}}\right)\,\text{d}\omega\\\
&=\int{\left(\bm{o}_{\text{exp}}-\bm{o}_{\text{mdl}}\right)}^{*}\,{\bm{\Delta}}^{*}\left(\left.{\bm{\Delta}}^{*}\right.^{-1}\bm{S}_{\text{n},o}^{-1}\,\bm{\Delta}^{-1}\right)\bm{\Delta}\left(\bm{o}_{\text{exp}}-\bm{o}_{\text{mdl}}\right)\,\text{d}\omega\\\
&=\int{\left(\bm{o}_{\text{exp}}-\bm{o}_{\text{mdl}}\right)}^{*}\,\bm{S}_{\text{n},o}^{-1}\,\left(\bm{o}_{\text{exp}}-\bm{o}_{\text{mdl}}\right)\;\text{d}\omega\leavevmode\nobreak\
.\end{split}$ (4.21)

where $\bm{f}_{\text{mdl}}$ and $\bm{f}_{\text{exp}}$ are the modeled and
experimental equivalent accelerations; $\bm{o}_{\text{mdl}}$ and
$\bm{o}_{\text{exp}}$ are the modeled and experimental displacement readouts.
In the preceding equation, $\bm{\Delta}$ and $\bm{\Delta}^{-1}$ are used to
transform the sensed relative motion into equivalent acceleration (and vice-
versa) and contain the dependence on the modeled parameters. The main benefit
of working with accelerations is the automatic subtraction of system
transients as described in Section 3.2 and that is numerically demonstrated at
the end of this chapter. Despite the identification in displacement where the
parameters are explicit in the modeled template signal, in the identification
in acceleration the parameters are implicit in the estimated acceleration.
Even though there is no real experimental acceleration because this must be
estimated from the displacement readouts, system identification in
acceleration domain can be still carried out numerically with a non-standard
approach based upon a closed-loop optimization over the estimated acceleration
data, whereas the modeled forces are the injected bias signals. For the rest,
the following discussion employs the estimation in the domain of displacement
readouts, as the other approach is currently under investigation.

The MIMO-Multi-Experiment Fisher information matrix for the parameter
estimates $\bm{p}_{\text{est}}$ is the local curvature of the log-likelihood
surface around the minimum and is given by

$\bm{\mathcal{I}}=\int{\bm{o}_{\text{i}}(\omega)}^{*}\,{\nabla_{\bm{p}}\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p}_{\text{est}})}^{*}\,\bm{S}_{\text{n}}(\omega)^{-1}\,\nabla_{\bm{p}}\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p}_{\text{est}})\,\bm{o}_{\text{i}}(\omega)\,\text{d}\omega\leavevmode\nobreak\
,$ (4.22)

where $\nabla_{\bm{p}}$ is the gradient with respect to all 7 system
parameters. As above, if $\bm{T}_{o_{\text{i}}\rightarrow o}$ is a $4\times
4$-matrix, then $\nabla_{\bm{p}}\bm{T}_{o_{\text{i}}\rightarrow o}$ is a
$7\times 4\times 4$-tensor and the information is a $7\times 7$-matrix as
required. The very high SNR regime of the signals in Figure 4.11 and Figure
4.12 assures that the linear approximation of (4.22) holds true and no
corrective terms arise as pointed out by [52] and more recently by [53]. As
the inverse of the information matrix provides the estimated covariance
matrix, the validity of the linear approximation is checked a posteriori in
Section 4.4.6 by inspecting the statistics of a Monte Carlo simulation.

####  4.4.3 Whitening

The colored noise behavior of a typical LPF run makes mandatory to decorrelate
the data used for system identification in order for a generic statistical
estimator be unbiased. Consider for example a stationary noisy time-series
$o(t)$ with noise PSD $S_{\text{n}}(\omega)$. The SNR of the signal [51] can
be recast as

$\begin{split}\rho^{2}&=\int\frac{{o}^{*}(\omega)\,o(\omega)}{S_{\text{n}}(\omega)}\,\text{d}\omega\\\
&=\int
o^{*}_{\text{w}}(\omega)\,o_{\text{w}}(\omega)\,\text{d}\omega\leavevmode\nobreak\
,\end{split}$ (4.23)

which can be viewed as the acting of the whitening filter
$W(\omega)=1/\sqrt{S_{\text{n}}(\omega)}$ on $o(\omega)$ to produce the
whitened series

$o_{\text{w}}(\omega)=W(\omega)\,o(\omega)\leavevmode\nobreak\ .$ (4.24)

Here “whitened” is equivalent to saying that the noise PSD of the filtered
series is approximately frequency-independent. The discrete time-domain
version of the preceding involves the noise covariance matrix
$\bm{C}_{\text{n}}$

$\begin{split}\rho^{2}&={\bm{o}}^{\mathsf{T}}\,\bm{C}_{\text{n}}^{-1}\bm{o}\\\
&={\bm{o}}^{\mathsf{T}}_{\text{w}}\,\bm{\Lambda}_{\text{n}}^{-1}\bm{o}_{\text{w}}\leavevmode\nobreak\
,\end{split}$ (4.25)

which again can be viewed as the acting of the whitening filter $\bm{W}$, an
orthogonal matrix satisfying
$\bm{C}_{\text{n}}^{-1}={\bm{W}}^{\mathsf{T}}\bm{\Lambda}_{\text{n}}^{-1}\bm{W}$
555In fact, if
$\bm{C}_{\text{n}}={\bm{U}}^{\mathsf{T}}\bm{\Lambda}_{\text{n}}\bm{U}$ where
$\bm{U}$ is an orthogonal matrix and $\bm{\Lambda}_{\text{n}}$ is the eigen-
decomposition of $\bm{C}_{\text{n}}$, then it turns out that
$\bm{U}^{-1}={\bm{W}}^{\mathsf{T}}$., on $\bm{o}$ to produce the whitened
unit-variance series

$\bm{o}_{\text{w}}=\bm{W}\bm{o}\leavevmode\nobreak\ .$ (4.26)

As above, “whitened” means that the process diagonalizes the covariance
matrix, so that $\bm{\Lambda}_{\text{n}}$ effectively becomes an identity
matrix.

For simulation and analysis purposes, whitening a time-series is formally the
inverse process of noise generation. Whitening filters are obtained by
performing a fit in the $z$-domain to the inverse of the estimated PSD
666Throughout this thesis, if not otherwise stated, it is assumed that a PSD
is estimated by means of the Welch (modified periodogram) method [54] using a
4-sample 92-dB Blackman-Harris window [55], 16-segments averaged, $66\%$
overlap and mean detrended.. Figure 4.13 reports an example of whitening
777Data filtering can produce fake transients at the beginning of the filtered
time-series. To avoid this possibility, an initial segment of data is usually
cut away. a typical 28-hour run of interferometric noise. The effect of the
whitening filters, as required, is to flatten the noise shapes, i.e., to
decorrelate the time-series.

Figure 4.13: Whitening of a simulated noise run. $o_{1}$ and $o_{12}$ are the
two interferometer readings with PSD reported on the basis of the scale on the
left end side. $o_{1,\text{w}}$ and $o_{12,\text{w}}$ are the whitened
counterparts with PSD reported on the basis of the scale at the right end
side. They show how the whitening filters can flatten the noise shapes. The
convolution with a low-pass filter of the data resampling from 10 to
$1\,\mathrm{Hz}$ is the cause of the drop around $0.5\,\mathrm{Hz}$.

Despite the PSD shapes which seem reasonably good at first sight, a residual
red component still persists. Table 4.3 reports two higher-order moments
(skewness and excess kurtosis) of the empirical distribution together with
their uncertainties [50]. By inspecting the values, it turns out that the
sample mean of the differential channel $o_{12}$ is not compatible with zero,
as one would expect. Usually, a first or second order polynomial fit is
necessary to subtract that residual component. The result is not surprising:
the intrinsical difficulty is that the whitening process is performed on a
restricted frequency band (the one of the estimated PSD) and low-frequency
components may survive after the filtering.

Table 4.3: Sample mean $\mu$, standard deviation $\sigma$ and higher moments, the sample skewness $\gamma_{1}$ and the excess kurtosis $\gamma_{2}$, for the whitened data channels $o_{1}$ and $o_{12}$. Assuming Gaussian-distributed data, the approximate standard deviations are $\sigma_{\mu}\simeq\sigma/\sqrt{N}$, $\sigma_{\sigma}\simeq\sigma/\sqrt{2N}$, $\sigma_{\gamma_{1}}\simeq\sqrt{\nicefrac{{6}}{{N}}}$, $\sigma_{\gamma_{2}}\simeq\sqrt{\nicefrac{{24}}{{N}}}$, with $N$ the number of data samples. Data | $\mu$ | $\sigma$ | $\gamma_{1}$ | $\gamma_{2}$
---|---|---|---|---
$o_{1,\text{w}}$ | $0.008\,\pm\,0.003$ | $0.970\,\pm\,0.002$ | $(-5\,\pm\,8)\\!\times\\!10^{-3}$ | $(0\,\pm\,2)\\!\times\\!10^{-2}$
$o_{12,\text{w}}$ | $-0.254\,\pm\,0.003$ | $1.002\,\pm\,0.002$ | $(0\,\pm\,8)\\!\times\\!10^{-3}$ | $(3\,\pm\,2)\\!\times\\!10^{-2}$

The extent to which the idea of this section holds true depends on the
assumption of stationarity and Gaussianity. Even though for LPF the
interferometric noise is not explicitly dependent on the system parameters, it
may depend implicitly through the coupling between the external force noise
and the system response. Yet, as it will be discussed later in this chapter,
the estimated equivalent acceleration noise depends explicitly on the system
parameters through the transfer matrix given by the differential operator.

As a matter of fact, a non-stationarity in any of the system parameters
implies a non-stationarity in the noise. In fact, if $o=o\left(t,p(t)\right)$
is a generic interferometer readout depending, for simplicity, on just one
parameter fluctuating of $\delta p$ around the nominal value $p_{0}$, then to
first order $o\simeq o_{0}+o^{\prime}\,\delta p$, where $o_{0}=o(t,p_{0})$ and
$o^{\prime}=\left.\partial o(t,p)/\partial o\right|_{p_{0}}$. For a zero-mean
process the total variance is

$\text{Var}[o]\simeq\text{Var}[o_{0}]+\text{Var}^{\prime}[o_{0}]\,\delta
p+\text{Var}[o^{\prime}]\,\delta p^{2}\leavevmode\nobreak\ ,$ (4.27)

where the linear and quadratic terms come from the covariance between $o_{0}$
and $o^{\prime}$ and the variance of $o^{\prime}$ itself (see Appendix A.5 for
details). Therefore, if any of the system parameters fluctuates, noise is
likely to become non-stationary. In LPF all PSDs must be estimated piecewise
along data segments approximately stationary on a timescale given by the one
of the fluctuating parameter. The converse, i.e., a non-stationarity in the
noise implies a non-stationarity in any of the parameters is not assured,
since other effects, independent from those parameters, may still be relevant.
For example, Section 4.4.8 describes the possible existence of glitches, a
non-stationary behavior in the noise, and its impact to system identification.
Instead, Appendix A.6 introduces the time-frequency approach in the study of
non-stationarity noise.

####  4.4.4 Search algorithm

The joint log-likelihood (4.18) for two experiments, two readouts each, is
implemented in time domain by means of FFT/IFFT the whitened time-series. The
relevant iteration steps of the process taking to the final estimates of the
system parameters, in loop of increasing accuracy, are:

  1. 1.

the whitening filters are estimated on a long noise run, as in Section 4.4.3;

  2. 2.

the interferometric readouts of each experiments are whitened;

  3. 3.

the templates are generated according to (4.10) for the current parameter
values;

  4. 4.

the templates are whitened;

  5. 5.

the log-likelihood is evaluated, i.e., “models fit the data”, for the current
parameter values;

  6. 6.

the parameter values are updated according to the adopted optimization scheme.

From the optimization viewpoint, the log-likelihood is named the merit
function, i.e., the one being minimized as the parameter values are updated.
Figure 4.14 shows a sketch of the whole process of system identification. The
data production provides for the noise run and the experiments, with both
interferometric readouts and injected biases. Instead, the modeling provides
for the proper transfer matrix being used for simulating the template signals.
Finally, the data analysis concerns the estimation of the whitening filters
and the algorithm for the log-likelihood minimization.

Figure 4.14: Sketch of the system identification process for the two simulated
experiments along the optical axis. Noise run and experiments pertain to data
production. The modeling provides for the transfer matrix being used for
simulation and analysis. For system identification, data analysis comprises
the estimation of whitening filters and the log-likelihood optimization. The
estimated parameters are output together with their covariance matrix.

The algorithm performs a log-likelihood minimization by taking advantage of
the most recent developments in numerical non-linear optimization [50]. During
the work for this thesis, an investigation of different optimization
algorithms was carried out. Non-standard schemes like the simulated annealing,
genetic algorithms and the pattern search, with or without a multi start (an
initial Monte-Carlo-like exploration of the parameter space in which the
initial most likely points are taken into account for further processing),
were considered for the purposes of system identification. They were also
compared to a mixed strategy employing more standard and widely-used
optimization algorithms applied in sequence:

  1. 1.

the preconditioned conjugate gradient search (alternatively, the quasi-Newton
method) explores the parameter space to large scales;

  2. 2.

the derivative-free simplex allows to reach the required numerical accuracy.

The key advantage of mixing different approaches is that the global structure
of the parameter space can be explored while keeping the numerical accuracy.
Such an investigation proved that for the LPF system identification non-
standard schemes have comparable performances with respect to the one proposed
above which is assumed for the rest. The optimization is numerically
controlled and stopped until either the function tolerance or the average
parameter tolerance meets the requirement of $1\\!\times\\!10^{-4}$. The final
parameter estimates are output from the fitting tool, together with the
estimated covariance matrix, obtained by inverting the Fisher information
matrix (4.22) around the minimum.

Before showing the application to data simulated specifically for LPF, the
tool was checked against more simpler cases like the linear fit (which is
analytically solvable to any order), the chirped sine and the harmonic
oscillator [56, 57]. The result is that there are no systematic errors and
parameter uncertainties are in accordance to a Monte Carlo simulation. A
similar check for LPF is discussed in Section 4.4.6.

####  4.4.5 ESA simulator

A very important test-bench on both system modeling and validation of the
estimation techniques is the analysis of realistic data, closer to the actual
LPF mission than the ones simulated and shown in this work. A real LPF
simulator, named Off-line Simulation Environment (OSE), provided by ESA and
written by ASTRIUM has given the chance to promptly analyze the data as they
were realistically produced during the mission. The OSE is a state-space
representation of a 3-dimensional LPF model written under the MATLAB® and
Simulink® [58] environments. It contains the most relevant disturbances and
noise sources, the same actuation algorithms for drag-free, electrostatic
suspension and attitude controls (DFACS) embedded in LPF, all couplings within
the dynamics along the optical axis and between different degrees of freedom.
The OSE was written to mainly check all procedures, the mission timeline, the
experiments and validate the noise budget.

Several extended data analysis operational exercises were called in the past 2
years, where parameter estimation had a pivotal role, and therefore very
similar to a mock data challenge, where data production is strictly separated
from data analysis. The operational exercises culminated with the sixth one
targeted to parameter estimation using a linear fit with singular value
decomposition, a Markov-chain Monte Carlo method and the one described in this
thesis. The first application of system identification on that operational
exercise is contained in [59]. The final conclusion of the activity on the
same exercise is recently described in [60]. The three methods are apparently
in good agreement to each other, particularly the first and third approaches,
but an investigation of the fit residuals, like the one in Figure 4.4.7, shows
a mismatch in the first experiment between data and model at high frequency.
The fact is confirmed by a statistical comparison between the residual PSDs to
a noise-only measurement with a very general and model-independent method
based on the Kolmogorov-Smirnov test [61]. The explanation of such a mismatch
will be given in the near future with further operational exercises and much
more detailed knowledge of the simulator.

####  4.4.6 Monte Carlo validation

The aim of this section is to statistically validate the estimation method
presented so far. A Monte Carlo simulation of $1000$ different noise
realizations is used to check for consistency of the method. The estimation is
identically repeated at each step, enabling fine tuning and the study of the
statistics for every system parameter.

Table 4.4 reports on the comparison between the mean best-fit values and the
true values: the accordance is at the level of at most 2 standard deviations
and demonstrates that the estimation method is statistically unbiased.
Secondarily, it shows the best-fit standard deviations, i.e., the parameter
fluctuations due to noise, compared to the mean expected standard deviations
(the mean fit errors).

Table 4.4: Monte Carlo validation of 1000 independent noise realizations on which parameter estimation is repeated identically at each step. The mean best-fit values are compatible with the true values within 2 standard deviations. The terms in brackets are the error relative to the rightmost digit. The mean expected standard deviations (estimated from the fit) and the best-fit standard deviations are approximately the same order of magnitude. The mean log-likelihood is $\chi^{2}=0.96$ with $\nu=79993$. Parameter | True | Mean | Best-fit | Mean
---|---|---|---|---
| best-fit | st. dev. | exp. st. dev.
$\omega_{1}^{2}\,[10^{-6}\,\text{s}^{-2}]$ | $-1.303$ | $-1.303006(7)$ | $2\times 10^{-4}$ | $1\times 10^{-3}$
$\omega_{12}^{2}\,[10^{-6}\,\text{s}^{-2}]$ | $-0.698$ | $-0.697998(6)$ | $2\times 10^{-4}$ | $5\times 10^{-4}$
$S_{21}\,[10^{-4}]$ | $0.9$ | $0.90004(9)$ | $3\times 10^{-3}$ | $4\times 10^{-3}$
$A_{\text{df}}$ | $1.003$ | $1.00297(1)$ | $4\times 10^{-4}$ | $4\times 10^{-4}$
$A_{\text{sus}}$ | $0.9999$ | $0.9999001(1)$ | $4\times 10^{-6}$ | $2\times 10^{-5}$
$\Delta t_{1}\,[\text{s}]$ | $0.06$ | $0.059995(3)$ | $9\times 10^{-5}$ | $3\times 10^{-4}$
$\Delta t_{12}\,[\text{s}]$ | $0.05$ | $0.05000(3)$ | $8\times 10^{-4}$ | $1\times 10^{-3}$

Figure 4.15 shows a more in-depth analysis of all parameter statistics. The
accordance between the sample statistics of the Monte Carlo simulation and the
scaled theoretical Gaussian Probability Density Function (PDF) (evaluated at
the sample mean and standard deviation) is self-evident and demonstrates that:
(i) the estimation is statistically unbiased; (ii) the parameters are Gaussian
distributed.

|  |
---|---|---
(a) | (b) | (c)

|  |  |
---|---|---|---
(d) | (e) | (f) | (g)

Figure 4.15: Monte Carlo validation of 1000 independent noise realizations on
which parameter estimation is repeated identically at each step. The plots
show the statistics for all parameter estimates (a)-(g). The scaled Gaussian
PDF is evaluated at the sample mean (dashed vertical lines) and sample
standard deviation (half horizontal bars), which are compared to the true
values (solid vertical lines).

Analogously, Figure 4.16 shows the statistics for the estimated variances.
Theory prescribes that the variance must be $\chi^{2}$ distributed, but for
$\nu=79993$ the $\chi^{2}$ distribution tends to a Gaussian distribution with
very good approximation, as is clear from the plots.

|  |
---|---|---
(a) | (b) | (c)

|  |  |
---|---|---|---
(d) | (e) | (f) | (g)

Figure 4.16: Monte Carlo validation of 1000 independent noise realizations on
which parameter estimation is repeated identically at each step. The plots
show the statistics for all parameter variances (a)-(g). The scaled Gaussian
PDF is evaluated at the sample mean and standard deviation.

Appendix A.7 also discusses some other interesting features of the Monte Carlo
statistics, like the parameter correlation, related to the rotation of the
log-likelihood paraboloid principal axis around the minimum, and the scatter
of the estimation chains due to the noise fluctuation.

The final, and most remarkable check, is the comparison between the fit
$\chi^{2}$ log-likelihood and the one calculated on pure noise data contained
in Figure 4.17. It is worth stressing that both the fit and the noise
$\chi^{2}$ showed agreement between each other, but they were both positively
skewed in a preliminary Monte Carlo simulation. The following facts explain
why. Section 4.4.3 has discussed the practical method to implement the
diagonalization of the noise covariance matrix with its main limitation. This
consists in the impossibility of filtering out the lowest frequencies, due to
the finiteness of the data stretches from which whitening filters are derived
and which causes the skewness. Transparently, the application of a high-pass
filter to the data has solved the issue. The comparison in the plot provides
for an important twofold test: on one side, the parameter variances are
statistically distributed as the fit $\chi^{2}$ log-likelihood, as required;
on the other, the fit $\chi^{2}$ log-likelihood is in agreement with the noise
$\chi^{2}$ log-likelihood, showing that the estimation method has
statistically suppressed the deterministic signals and recovered the noise
statistic with no extra bias.

|
---|---
(a) |  (b)
Figure 4.17: Monte Carlo validation of 1000 independent noise realizations on
which parameter estimation is repeated identically at each step. The plots
show the statistic for (a) the fit $\chi^{2}$ log-likelihood and (b) the noise
$\chi^{2}$ log-likelihood. The agreement between the two demonstrates that the
deterministic signals are statistically suppressed out of the data.

####  4.4.7 Non-standard scenario: under-performing actuators and under-
estimated couplings

System identification has a key role in compensating the SC jitter and the TM
couplings. Even in the unlikely (but possible) situation of under-performing
actuators or under-estimated force couplings, it is still possible to retrieve
the actual parameter values and allow for a precise estimation of the total
equivalent acceleration noise without loosing sensibility and getting into
systematic errors. The impact on the estimation of the total equivalent
acceleration noise will be illustrated in Section 4.5.

To introduce the problem, suppose that the predicted TM couplings are
$\omega_{1}^{2}=-1.3\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$ and
$\omega_{12}^{2}=-0.7\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$ and during the LPF
mission:

  1. 1.

the actual TM couplings are about two times the predicted ones, due to
unexpected/unmodeled stronger forces, like
$\omega_{1}^{2}=-3\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$ and
$\omega_{12}^{2}=-2\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$;

  2. 2.

the thruster and capacitive actuators unfortunately misfunction, due to both a
breakdown of one or more thruster clusters and a loss of efficiency in the
capacitive actuators on the second TM; this situation can be described by
gains sensitively lower than one, like $A_{\text{df}}=0.62$ and
$A_{\text{sus}}=0.6$;

  3. 3.

the interferometer introduces an extra cross-talk,
$S_{21}=1.5\\!\times\\!10^{-3}$, ten times the expected one $S_{21}\raise
0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}1\\!\times\\!10^{-4}$.

In this very unfortunate situation, system identification, see Table 4.5,
allows for the estimation of the true values within 1 standard deviation from
the true values, so maintaining precision, even though the optimizations
starts from initial guesses which are typically $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}10^{3}$ standard deviations away, so guaranteeing
accuracy too.

Table 4.5: Robustness to a non-standard scenario: under-performing actuators / under-estimated couplings. Initial estimates (guess) at $\chi^{2}=1.3\\!\times\\!10^{5}$, $\nu=79193$; best-fit values at $\chi^{2}=0.99$. The term in brackets is the error relative to the rightmost digit. In curly brackets the bias (absolute deviation from the real value in units of standard deviation) for each estimate. Parameter | True | Best-fit | Guess
---|---|---|---
$\omega_{1}^{2}\,[10^{-6}\,\text{s}^{-2}]$ | $-3$ | $-2.9998(2)$ | $\\{1.1\\}$ | $-1.3$ | $\\{7.8\times 10^{3}\\}$
$\omega_{12}^{2}\,[10^{-6}\,\text{s}^{-2}]$ | $-2$ | $-2.0000(1)$ | $\\{0.32\\}$ | $-0.7$ | $\\{1.0\times 10^{4}\\}$
$S_{21}\,[10^{-3}]$ | $-1.5$ | $-1.4998(1)$ | $\\{0.55\\}$ | $0$ | $\\{4.7\times 10^{3}\\}$
$A_{\text{df}}$ | $0.62$ | $0.61994(8)$ | $\\{0.77\\}$ | $1$ | $\\{4.9\times 10^{3}\\}$
$A_{\text{sus}}$ | $0.6$ | $0.599990(8)$ | $\\{1.3\\}$ | $1$ | $\\{5.1\times 10^{4}\\}$
$\Delta t_{1}\,[\text{s}]$ | $0.6$ | $0.6013(7)$ | $\\{1.8\\}$ | $0$ | $\\{8.4\times 10^{2}\\}$
$\Delta t_{12}\,[\text{s}]$ | $0.4$ | $0.398(2)$ | $\\{0.95\\}$ | $0$ | $\\{2.3\times 10^{2}\\}$

Figure 4.18 elucidates much more the results, showing the overall performances
of the estimation. The $\chi^{2}$ is reduced from $1\\!\times\\!10^{5}$ to
$\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}1$ – the required
optimum – within the given tolerances (set to $1\\!\times\\!10^{-4}$ in both
log-likelihood and parameter values), while keeping both accuracy and
precision. The figure reports two examples of estimation chains (for
$\omega_{1}^{2}$ and $\omega_{12}^{2}$), showing the correlation with the big
jumps in the $\chi^{2}$ chain and how the parameters saturate to the optimum
values. The estimation, as already said, is divided into two phases: a
gradient-based search, spanning the global structure of the parameter space,
and a simplex search, improving the final accuracy.

---
(a)
|
(b) |  (c)
Figure 4.18: Robustness to a non-standard scenario: under-performing actuators
/ under-estimated couplings. The estimation performances relative to the log-
likelihood minimization (a) from $\raise 0.58554pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}1\\!\times\\!10^{5}$ to the optimum $\raise
0.58554pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}1$ and two examples of
estimation chains for (b) $\omega_{1}^{2}$ and (c) $\omega_{12}^{2}$ showing
the correlation with the big jumps in the $\chi^{2}$ chain. A preliminary
global gradient search is followed by a local simplex. The process lasts for
1636 iterations and stops when the required tolerance is met.

The final and most important discussion is the analysis of residuals
summarized in Figure 4.19 for both identification experiments and
interferometric readouts. The estimated PSDs of both initial and best-fit
residuals are compared to the PSDs of an independent noise run. It is clear
that the deterministic signals are completely subtracted from the data, hence
recovering the noise shapes for all experiments and readouts. The improvement
is mostly evident at low frequency: for $o_{12}$ the residuals are suppressed
by $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}3$ orders of
magnitude around $1\,\mathrm{mHz}$. The same happens for $o_{1}$ in the first
experiment where the improvement is $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}2$ orders of magnitude. Only $o_{1}$ in the second
experiment shows no improvement for the reason already discussed in Section
4.3 (the signal is negligible).

|
---|---
(a) |  (b)
|
(c) |  (d)
Figure 4.19: Robustness to a non-standard scenario: under-performing actuators
/ under-estimated couplings. Analysis of residuals for all simulated
identification experiments and interferometric readouts. Initial and best-fit
residuals are compared to the expected noise shapes estimated from an
independent run. For $o_{12}$ the improvement in both experiments (b) and (d)
is of $\raise 0.58554pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}3$ orders of
magnitude around 1 mHz; (a) for $o_{1}$ in the first experiment is $\raise
0.58554pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}2$ orders of magnitude; (c)
contains no signal.

####  4.4.8 Non-standard scenario: non-Gaussianities

This section is devoted to showing the impact of non-Gaussianities in the
noise to parameter estimation. The main realistic behavior of experimental
noise is the possible presence of outliers: consequently, the sampling
distribution of the data may show some prominent tails. An example of such
outliers is the manifestation of glitches, very short noise transients due to
anomalous response in the readout/circuitry.

Given the non-Gaussian components in the noise, the log-likelihood defined so
far is no longer well-behaved. Because of the intrinsical assumption of
Gaussianity, it usually overweighs the outliers, and a systematic error may
arise. A standard approach, named local L-estimate [50] 888“L” stands for
“likelihood”., requires the generalization of the definition of log-
likelihood. The idea is to properly take care of the outliers by regularizing
the usual square of whitened residuals with other similar definitions by means
of a weighting function $\rho$

$\chi^{2}=\sum_{i}\rho(r_{\text{w},i})\leavevmode\nobreak\ ,$ (4.28)

where, as an example, three possible choices, the squared, absolute and
logarithmic deviations, are considered

$\rho(r_{\text{w},i})=\begin{cases}r_{\text{w},i}^{2}&\text{mean squared
dev.}\\\ \left|r_{\text{w},i}\right|&\text{mean absolute dev.}\\\
\log(1+r_{\text{w},i}^{2})&\text{mean logarithmic
dev.}\end{cases}\leavevmode\nobreak\ ,$ (4.29)

corresponding to the cases of data distributed according to Gaussian, log-
normal and Lorentzian distribution, respectively. The subscript $i$ is a
generalized index counting the data available from all experiments and
interferometric readouts and $r_{\text{w},i}$ is the whitened time-series of
residuals. Figure 4.20 compares the three weighing functions for residuals out
to 5 standard deviations. As is clear, the squared deviation overweighs the
outliers. The absolute deviation gives a slightly better weight at high
deviations, but performs poorly at low deviations. The logarithmic deviation
has much more flexibility as it behaves like the squared deviation at low
deviations and performs better than the absolute deviation.

Figure 4.20: Comparison of the three weighing functions of (4.29) for the the
proper weighing of outliers in the data. The logarithmic deviation is the most
accurate as it behaves like the squared deviation at low deviations and
performs better than the absolute deviation.

The method can be successfully applied to data with glitches. Noise glitches
are unpredictable high-frequency noise transients mostly due to failures in
the circuitry. Such outliers usually fall well beyond 3 standard deviations
and produce an excess at the tails of the statistic. Since the output of the
interferometer might be subject to similar phenomena, this section presents
the results of the investigation of a realistic experiment containing
glitches. Such transients are modeled as sine-Gaussian functions

$o_{\text{gl}}(t)=a\sin\left[2\pi
f_{0}(t-t_{0})\right]\exp\left[-\frac{(t-t_{0})^{2}}{\tau^{2}}\right]\leavevmode\nobreak\
,$ (4.30)

where the glitch parameters span a wide (uniformly distributed) range of
values. In particular, the glitch frequency, $f_{0}$, covers the whole
bandwidth $(10^{-4}\text{--}0.45)\,\mathrm{Hz}$; the injection time, $t_{0}$,
is distributed all along the time-series; the characteristic time, $\tau$,
giving the typical duration of the pulse is $(1\text{--}2)\,\mathrm{s}$; the
amplitude, $a$, falls outside the Gaussian statistic by $(3\text{--}20)$ noise
standard deviations. Moreover, the number of glitch injections is fixed as a
fractional part of the whole data series, conventionally choosing
$f_{\text{gl}}=N_{\text{gl}}/N_{\text{data}}=1\%$, since higher values are
very unlikely. Notice that this value represents only the number of
injections: the actual fraction of corrupted data is the order of
$3\,\text{E}[\tau]\,f_{\text{gl}}\simeq 5\%$.

Glitchy noise is readily produced by coloring a white, zero-mean, unitary
standard deviation input time-series, as in Section 4.2, corrupted by random
injections of glitches. Figure 4.21 shows how glitches appear in the
interferometric differential readout and in the estimated PSDs, compared to
the original noise stretches. The effect of glitches is that the PSD of the
simulated noise scales linearly with the frequency, up to
$4\\!\times\\!10^{-9}\,\mathrm{m\,Hz^{-\nicefrac{{1}}{{2}}}}$ and
$6\\!\times\\!10^{-11}\,\mathrm{m\,Hz^{-\nicefrac{{1}}{{2}}}}$ around
$0.2\,\mathrm{Hz}$ for the first and differential readout, respectively. This
excess noise sums up to the original one and is shown as high-frequency
components. Obviously, the noise statistic contains an excess at the tails.
For example, $o_{1}$ has an excess kurtosis of $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}19$, compared to the original one of
$-9\\!\times\\!10^{-3}$. No significant difference in skewness is detected
since the statistic does not loose symmetry with the glitch injections.

|
---|---
(a) |  (b)
Figure 4.21: Robustness to a non-standard scenario: non-Gaussianities. (a)
simulated original and glitchy noise for $o_{12}$; (b) PSDs of the simulated
original and glitchy noise for $o_{1}$ and $o_{12}$. The level of data
corruption is evident and glitches appear as high-frequency bumps around
$0.2\,\mathrm{Hz}$.

Whitening filters are derived from the glitchy noise stretches with the same
procedure described in Section 4.4.3. However, since the whitening process
works assuming stationarity, glitches are not filtered out from the data.

Table 4.6 shows the results of three different parameter estimations with the
definitions of the weighting functions in (4.29).

Table 4.6: Robustness to a non-standard scenario: non-Gaussianities. The
comparison between three parameter estimations with the three definitions in
(4.29). $\nu=79193$. The term in brackets is the error relative to the
rightmost digit. In curly brackets the bias (absolute deviation from the real
value in units of standard deviation) for each estimate.

Parameter | Real | Best-fit | Best-fit | Best-fit | Guess
---|---|---|---|---|---
| (mean sq. dev.) | (mean abs. dev.) | (mean log. dev.)
| $\chi^{2}=10$ | $\chi^{2}=2.1$ | $\chi^{2}=0.95$
$\omega_{1}^{2}\,[10^{-6}\,\text{s}^{-2}]$ | $-1.32$ | $-1.320(1)$ | {0.061} | $-1.3188(6)$ | {2.0} | $-1.3192(4)$ | {2.0} | $-1.3$
$\omega_{12}^{2}\,[10^{-6}\,\text{s}^{-2}]$ | $-0.68$ | $-0.6798(7)$ | {0.29} | $-0.68000(3)$ | {0.011} | $-0.6804(2)$ | {1.8} | $-0.7$
$S_{21}\,[10^{-4}]$ | $1.1$ | $1.10(2)$ | {0.074} | $1.113(7)$ | {1.8} | $1.116(5)$ | {3.4} | $0$
$A_{\mathrm{df}}$ | $1.01$ | $1.011(3)$ | {0.29} | $1.010(1)$ | {0.23} | $1.0109(8)$ | {1.2} | $1$
$A_{\mathrm{sus}}$ | $0.99$ | $0.99000(5)$ | {0.035} | $0.98959(2)$ | {20} | $0.99001(1)$ | {0.99} | $1$
$\Delta t_{1}\,[\text{s}]$ | $0.1$ | $0.100(3)$ | {0.045} | $0.090(1)$ | {8.3} | $0.1007(8)$ | {0.90} | $0$
$\Delta t_{12}\,[\text{s}]$ | $0.1$ | $0.098(5)$ | {0.36} | $-0.0290(2)$ | {58} | $0.098(2)$ | {1.2} | $0$

The most conservative least square estimator provides overestimated errors
since they scale as $\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$
}$$}\sqrt{\chi^{2}}$. The absolute and logarithmic deviations provide better
statistics and lower errors, but the first gives biased estimates of
$A_{\text{sus}}$, $\Delta t_{1}$ and $\Delta t_{12}$ and the last one a
slightly biased estimate of $S_{21}$. The analysis of residuals demonstrates
that the three methods recover the noise shapes and are in agreement with each
other, so the systematic errors are only in the estimated parameters. These
estimators are also $30\%$ and $9\%$ faster than the Gaussian (mean squared
deviation), as the outliers have less influence on the estimation chains.

By inspecting the results, it turns out that there is no absolute rule that
can be applied when dealing with glitches. However, from the differences
between the estimates it is possible to infer the sensitivity of each single
parameter to glitches. For example, adopting the ratio between the biases as
the a-posteriori criterion for comparing two methods, it tends to one if that
parameter is not sensitive to glitches; otherwise, it tends to a very small or
very large number. In view of this consideration, the comparison between the
mean squared deviation and the mean logarithmic deviation gives that $S_{21}$
is the most sensitive parameter, whereas $\Delta t_{12}$ the least.

Starting from the fact that the three methods give the same results for purely
Gaussian noise, a proposed recipe is the following:

  1. 1.

apply the conservative approach (the ordinary mean squared deviation) directly
to corrupted time series and try with different estimators (mean absolute
deviation, mean logarithmic deviation, etc.);

  2. 2.

start removing some outliers giving them negligible weight;

  3. 3.

redo the analysis with all estimators;

  4. 4.

check for convergence and agreement between the estimators.

The overall process can be actually viewed as a reweighing analysis providing
for robust uncertainties and, at the same time, the removal of outliers in a
step-by-step smooth readjustment. Even though it would be possible in
principle to clean up the data just before the estimation, in that case the
results would likely be dependent on the statistical criterion used for such
cleaning. Even though it is beyond the scope of this thesis to implement the
idea, it is worth observing that the two main advantages of the preceding
recipe are its robustness in definition and the fact that data polishing is
smooth and model independent.

###  4.5 Estimation of total equivalent acceleration noise

This section justifies the efforts in developing the techniques introduced so
far with all tests and validation runs, showing the impact of system
identification on the estimation of the total equivalent acceleration noise.
As said throughout this thesis, the main objective of the LPF mission in view
of a real GW astronomy with spaced-based detectors is the characterization of
the Doppler link as the fundamental spacetime meter in terms of equivalent
differential acceleration. Even if LPF is different in design with respect to
LISA – no faraway optical measurement between two SCs is actually implemented
– yet the principle and, most of all, the performances in sensitivity can be
extrapolated and gather more confidence in the scientific scopes of any
spaced-based GW detector.

Assessing the performance in sensitivity as equivalent input acceleration
noise is a very effective way to put dynamics, sensing and control on the same
footing as described in Section 2.3. This can be achieved by means of the
$\bm{\Delta}$ operator of Section 3.1, connecting interferometric displacement
readouts to total equivalent acceleration and at the same time compensating
for TM couplings, SC jitter and sensing cross-talk.

Suppose that $\bm{S}_{\text{n},o}(\omega,\bm{p}_{\text{true}})$ is the
measured interferometric noise PSD. Then, the estimated total equivalent
acceleration noise PSD is given by

$\bm{S}_{\text{n},f}(\omega,\bm{p}_{\text{est}})=\bm{\Delta}(\omega,\bm{p}_{\text{est}})\,\bm{S}_{\text{n},o}(\omega,\bm{p}_{\text{true}})\,{\bm{\Delta}(\omega,\bm{p}_{\text{est}})}^{*}\leavevmode\nobreak\
,$ (4.31)

where $\bm{\Delta}(\omega,\bm{p}_{\text{est}})$ models the transfer from
interferometric displacement readouts to total equivalent acceleration and
$\bm{p}_{\text{est}}$ are the parameter estimates as obtained by system
identification. It is worth noting that if $\bm{S}_{\text{n},o}$ was assumed
constant to the parameter values in first approximation, the transfer to total
equivalent acceleration would anyhow couple the output noise with the dynamics
so that the estimated total equivalent acceleration noise becomes explicitly
dependent on the parameter values. This shows that parameter estimation serves
not only for system identification, but also for the actual identification of
the total equivalent acceleration noise.

Furthermore, suppose that
$\bm{p}_{\text{est}}\simeq\bm{p}_{\text{true}}+\delta\bm{p}$, with
$\delta\bm{p}$ the parameter biases being much larger than the statistical
uncertainties on $\bm{p}_{\text{est}}$. It is easy to show that the parameter
biases propagate to the differential operator
$\bm{\Delta}_{\text{est}}\simeq\bm{\Delta}_{\text{true}}+\delta\bm{\Delta}$,
where $\bm{\Delta}_{\text{true}}=\bm{\Delta}(\omega,\bm{p}_{\text{true}})$ and
$\bm{\Delta}_{\text{est}}=\bm{\Delta}(\omega,\bm{p}_{\text{est}})$. Systematic
errors found in the parameter values produce systematic errors in the
recovered total equivalent acceleration noise

$\delta\bm{S}_{\text{n},f}\simeq\delta\bm{\Delta}\,\bm{S}_{\text{n},o}\,{\bm{\Delta}}^{*}+\bm{\Delta}\,\bm{S}_{\text{n},o}\,{\delta\bm{\Delta}}^{*}\leavevmode\nobreak\
,$ (4.32)

where the subscript “true” is dropped out for clearness. As pointed out in
[62], the statistical uncertainty on the parameter values are masked by the
statistical uncertainty on the estimated spectrum. Despite this, systematic
errors in the estimated parameters can fall well outside the confidence levels
of the optimal spectrum and show themselves as not mere excess noise, but
producing really different noise shapes. Hence, it is expected that the
estimation of the total equivalent acceleration noise is biased if the
parameter values are not correctly assessed from system identification.

To demonstrate the impact of system identification on the estimation of the
total equivalent acceleration noise, a very long noise run, $\raise
0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}6$ days, is simulated with the
same procedures of Section 4.2, i.e., by coloring a sequence of white Gaussian
input time-series with cross-correlating noise shaping filters. The
interferometric displacement noise model is derived in a non-standard
configuration of LPF, as in Section 4.4.7, namely in the case of stronger-
than-expected TM couplings, malfunctioning actuators and a higher sensing
cross-talk. In this case, the estimation of the total equivalent acceleration
noise with naively guessed parameter values will surely contain systematic
errors.

The estimation of the total equivalent acceleration noise is readily performed
on the multi-channel interferometric run with a scheme described in details in
[62, 41], by applying a time-domain version of the $\bm{\Delta}$ operator of
Section 3.1. The issues connected to numerical derivatives in LPF are
extensively discussed and solved in [63]. As said, system identification
effectively helps in the calibration of the operator. In support of the
statement, the numerical estimation of the total equivalent acceleration noise
is performed assuming three different parameter sets that can be found in
Table 4.5:

  1. 1.

the initial guess values, as it was without a preliminary system
identification: $\omega_{1}^{2}=-1.3\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$,
$\omega_{12}^{2}=-0.7\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$, $S_{21}=0$,
$A_{\text{df}}=1$, $A_{\text{sus}}=1$ (typically $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}10^{4}$ standard deviations away from the real
values);

  2. 2.

the best-fit values, as it was with a preliminary system identification, i.e.,
after having calibrated the differential operator:
$\omega_{1}^{2}=-2.9998(2)\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$,
$\omega_{12}^{2}=-2.0000(1)\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$,
$S_{21}=-1.4998(1)\\!\times\\!10^{-3}$, $A_{\text{df}}=0.61994(8)$,
$A_{\text{sus}}=0.599990(8)$;

  3. 3.

the true values, used for consistency checks:
$\omega_{1}^{2}=-3\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$,
$\omega_{12}^{2}=-2\\!\times\\!10^{-6}\,\mathrm{s^{-2}}$,
$S_{21}=-1.5\\!\times\\!10^{-3}$, $A_{\text{df}}=0.62$, $A_{\text{sus}}=0.6$.

The result of the analysis is contained in Figure 4.22, showing the total
equivalent differential acceleration noise, both numerically estimated and
modeled, for the three different cases.

Figure 4.22: Total equivalent differential acceleration noise numerically
estimated on synthetic data and compared to theoretical noise models obtained
by a full projection of fundamental noise sources. The estimation of the total
out-of-loop equivalent acceleration can be performed either with a preliminary
system identification or without it. The PSD estimated with a preliminary
system identification completely overlaps the one of a hypothetical estimation
assuming the knowledge of the true parameter values. The observed difference
shows that a preliminary system identification is mandatory to avoid
systematic errors in the reconstructed total equivalent acceleration noise.
The solid thinner lines indicate the reasons of such a discrepancy. Around
$50\,\mathrm{mHz}$ the bump is due to unsuppressed thruster noise exceeding
the interferometric $o_{12}$ readout noise. At low frequency and around
$0.4\,\mathrm{mHz}$, the two major contributions are the unsuppressed force
couplings between the TMs and the SC and the capacitive actuation noise.
Thanks to system identification, an improvement in performance of a factor 4
at low frequency is evident.

First, the agreement between modeled and estimated total equivalent
acceleration noise PSDs states that: (i) the generation of the interferometric
noise is accurate to the assumed models at least to within the statistical
uncertainty of the spectra; (ii) the numerical estimation of the total
equivalent acceleration in time-domain is accurately explained by the
frequency-domain transfer matrix from interferometric readouts to the total
equivalent acceleration.

Second but more important, the total equivalent acceleration noise estimated
with a preliminary system identification completely overlaps the one of a
hypothetical estimation assuming the complete knowledge of the true values.
Therefore it demonstrates that it is still possible to meet the sensitivity
requirements during under-performing mission operations.

The observed systematic errors in the total equivalent acceleration noise
estimated without identification show that system identification is strictly
mandatory to avoid such problems and guarantee the scientific objectives. The
systematic errors can be explained by the fact that the naive initial guess
values are sensitively different from the true values. Since the operator is
not calibrated on fiducial parameter values, it is not effective in
compensating, in turn, the SC jitter due to the thruster actuation noise, the
TM couplings and the capacitive actuation noise. In particular, around
$50\,\mathrm{mHz}$ the bump is the unsuppressed thruster noise exceeding the
interferometric $o_{12}$ readout noise: the effect is due to the uncalibrated
drag-free gain $A_{\text{df}}$. At low frequency and around
$0.4\,\mathrm{mHz}$, the major contributions are the coupling forces between
the TMs and the SC (two contributions, accounting for
$1.8\\!\times\\!10^{-13}\,\mathrm{m\,s^{-2}}$, almost the whole noise budget)
and the capacitive actuation noise
($7\\!\times\\!10^{-14}\,\mathrm{m\,s^{-2}}$): the effect is due to the
uncalibrated stiffness constants $\omega_{1}^{2}$ and $\omega_{12}^{2}$ and
the suspension gain $A_{\text{sus}}$.

The final improvement in the estimated total equivalent acceleration noise
with system identification is a factor $4$ around $0.4\,\mathrm{mHz}$ and a
factor $2$ around $50\,\mathrm{mHz}$ in units of $\sqrt{\text{PSD}}$. The
conclusion is that without a preliminary system identification – robust to
non-standard parameter values – the performance of the mission and the
characterization of the total equivalent acceleration noise would seriously be
compromised.

###  4.6 Suppressing transients in the total equivalent acceleration noise

This final section discusses on the suppression of system transients for
realistic data produced by the OSE and provided by ESA. Section 3.2 and in
particular (3.16) demonstrate that system transients can be suppressed to
within the accuracy to which the differential operator $\bm{\Delta}$ has been
calibrated on parameter values representative of the system. A supporting
example is provided in what follows.

Figure 4.23 shows the first 3 hours of a typical noise run of the OSE. In
complete realism, just after the TM release the system is firstly turned into
accelerometer mode, then into science mode (around
$1\\!\times\\!10^{4}\,\mathrm{s}$) 999It is worth recalling that in
accelerometer mode the TMs are both electrostatically suspended, whereas in
the main science mode one of the two is in drag-free. The resulting noise is
at least one order of magnitude lower in the second case, especially in the
differential readout.: transients appears as a direct consequence of the non-
zero initial conditions. In fact, the initial positions are $0.24\,\mathrm{\mu
m}$ ($o_{1}$) and $0.36\,\mathrm{\mu m}$ ($o_{12}$), whereas the estimated
velocities 101010Since only an order of magnitude is needed, a two point
forward difference is applied together with a low-pass filter with frequency
cut at $100\,\mathrm{mHz}$. are about $-500\,\mathrm{pm\,s^{-1}}$ ($o_{1}$)
and $-4\,\mathrm{pm\,s^{-1}}$ ($o_{12}$). The transient in $o_{1}$ lasts for
half a hour and in $o_{12}$ for about 2 hours – the timescale of typical
transients as predicted by Section 3.2.

Figure 4.23: The first 3 hours of a typical noise run of the OSE. After the TM
release, the system is firstly turned into accelerometer mode, then into
science mode around $1\\!\times\\!10^{4}\,\mathrm{s}$. The transient due to
non-zero initial conditions lasts for half a hour in $o_{1}$ and about 2 hours
in $o_{12}$.

The estimation of the total equivalent differential acceleration noise is
performed twice on the same data, including the initial transitory, assuming
each time a different set of parameter values modeling the system. On one
side, a fair approximation of those parameters – the so-called initial guess –
reproduces the situation in which the estimation of the total equivalent
acceleration is performed without a preliminary system identification, as in
the previous section. On the other side, a fiducial approximation of those
parameters – the so-called best-fit – reproduces the situation in which the
estimation of the total equivalent acceleration is performed with a
preliminary system identification. Moreover, the estimation is performed
parallelly on two data segments lasting $3\\!\times\\!10^{4}\,\mathrm{s}$
each: the first one just after the system is turned into science mode and
containing the transient state; the second one follows it and is driven by the
steady state. Figure 4.24 shows the estimated total equivalent differential
acceleration noise for the two segments and for the two sets of parameter
values. The comparison shows that the transient is suppressed, and there is no
relevant difference between the two segments.

|
---|---
(a) |  (b)
Figure 4.24: Total equivalent differential acceleration noise numerically
estimated on synthetic data produced by the OSE and shown in Figure 4.23. The
data are split into two segments: the first one just after the system is
turned into science mode and containing the transient state; the second one
follows it and is driven by the steady state. The estimation of the total
equivalent out-of-loop acceleration can be performed either with a preliminary
system identification (fiducial parameter values) or without it (approximate
parameter values) (lines with different colors) and apparently there is no
difference. (a) the transient is suppressed, compared to (b) where the system
is dominated by the steady state.

Figure 4.25 reports the PSDs of the estimated total equivalent differential
acceleration noise for the above time-series, i.e., assuming the two sets of
parameter values for both segments. At low frequency, the noise level of the
segment containing the transient is higher then the subsequent segment, but
system identification helps in suppressing part of the noise around
$1\,\mathrm{mHz}$. Below $0.7\,\mathrm{mHz}$ there is an evidence that there
is an unsuppressed residual transitory in the data.

Figure 4.25: Total equivalent differential acceleration noise numerically
estimated on synthetic data produced by the OSE and shown in Figure 4.24. The
estimation of the total equivalent out-of-loop acceleration can be performed
either with a preliminary system identification (fiducial parameter values) or
without it (approximate parameter values) on both segments: the first one
dominated by the transient state and the second one dominated by the steady
state. System identification helps in suppressing the transient around
$1\,\mathrm{mHz}$.

The results of this section demonstrate how transients due to initial
conditions can be suppressed with reasonably good approximation in the total
equivalent acceleration time-series. The accuracy to which the suppression is
effective depends on the accuracy to which the parameter values of the system
are known. As shown, system identification helps in mitigating the effect due
to transients in the data.

##  5 Design of optimal experiments

The previous chapter introduced system identification and its relevance for
the unbiased estimation of the total equivalent acceleration noise of the LPF
mission. A standard series of sine injections spanning the frequency band was
utilized therein. Such bias injections make possible the estimation of the
modeled system parameters the total equivalent acceleration noise depends on.
Clearly, the goal of system identification is the parameter accuracy. A
possible approach is also to search for optimized stimuli to assess the system
parameters with better precision, with the final aim of a better estimate of
the equivalent acceleration noise. The relevance is worth that this chapter
addresses the question and provides for a solution.

###  5.1 Review of the problem

Aiming at discriminating among different designs of the same system
identification experiment, it is a rather natural consequence to enter into
the field of the optimal design of experiments [64, 65]. This matter tries to
answer to those physical problems characterized by a design matrix that shall
be maximized in order to perform a targeted measurement with an optimized
precision. This links to some very recent examples of practical applications
of the optimal design theory, multidisciplinary and covering very different
research fields: from dynamical systems [66], to geophysics [67], quantum
state tomography [68] and even magnetic resonance in medical engineering [69].
A general review can be found in [70, 71]. In what follows the same philosophy
is applied to the LPF mission with its peculiarities: the level of complexity
is worth as the very example of a MIMO multi-degree-of-freedom dynamical
system with coupled closed loops and subjected to various constraints.

As described in the previous chapter, system identification is targeted to
measuring the system parameters $\bm{p}$ appearing within the transfer matrix
$\bm{T}_{o_{\text{i}}\rightarrow o}(\omega,\bm{p})$ connecting applied
controller biases to interferometric readouts, for the case of the
investigation along the optical axis. If the inputs $\bm{o}_{\text{i}}$ are
parameterized by a set of parameters $\bm{\theta}$, then the Fisher
information matrix in (4.22) becomes

$\bm{\mathcal{I}}(\bm{\theta})=\int{\bm{o}_{\text{i}}(\omega,\bm{\theta})}^{*}\,{\nabla_{\bm{p}}\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p}_{\text{est}})}^{*}\,\bm{S}_{\text{n}}(\omega)^{-1}\,\nabla_{\bm{p}}\bm{T}_{o_{\text{i}}\rightarrow
o}(\omega,\bm{p}_{\text{est}})\,\bm{o}_{\text{i}}(\omega,\bm{\theta})\,\text{d}\omega\leavevmode\nobreak\
.$ (5.1)

Requiring that the estimates $\bm{p}_{\text{est}}$ should be given with the
optimal precision implies that the preceding matrix must be optimized with
respect to the design given by $\bm{\theta}$.

Theory provides for a solution of the problem. In fact, the optimal design is
attained by building up a scalar estimator on the information matrix,
$\phi[\bm{\mathcal{I}}]$, which is mathematically a functional over that
matrix. With the parametrization introduced above, the functional simply
becomes a scalar function of $\bm{\theta}$

$\phi[\bm{\mathcal{I}}]=\phi(\bm{\theta})\leavevmode\nobreak\ ,$ (5.2)

for given noise PSDs, interferometric readouts and estimated system
parameters. Hereafter, three different choices of the functional $\phi$ are
considered

$\phi(\bm{\theta})=\begin{cases}\det(\bm{\mathcal{I}}(\bm{\theta}))&\text{D
optimality}\\\ \min(\text{eig}(\bm{\mathcal{I}}(\bm{\theta})))&\text{E
optimality}\\\ \text{tr}(\bm{\mathcal{I}}(\bm{\theta}))&\text{T optimality}\\\
\end{cases}\leavevmode\nobreak\ ,$ (5.3)

and the corresponding for the covariance matrix, obtained directly inverting
(5.1), since maximum information is equivalent to minimum variance. The
interpretation of each single criterion is readily discussed. The D optimality
is the determinant of the information matrix and averages the information
along all terms, diagonal and off-diagonal. The E optimality takes the minimum
eigenvalue and tries to balance it with the others, hence regularizing the
conditional number of the matrix 111The conditional number is defined as the
ratio between the minimum and maximum eigenvalues. It expresses the
sensitivity of the matrix to numerical inversions. Round-off errors affect the
operation when the conditional number is either very small or very large. A
number of order 1 is considered stable to inversions.. The T optimality gives
to the diagonal the highest weight and corresponds to the averaged information
along all parameters. Even though quite different in definitions, the criteria
share the same philosophy: maximizing/minimizing the information/covariance
volume in the system parameter space around the minimum.

For LPF there is one more point adding much more complexity. The typical
constraints that must be met during all operations and especially for the
experiment design are:

  1. 1.

the general shape of the biases being injected;

  2. 2.

the dynamical range of the interferometer, $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}100\,\mathrm{\mu m}$;

  3. 3.

the force authority for thruster actuation, $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}100\,\mathrm{\mu N}$;

  4. 4.

the force authority for capacitive actuation, $\raise 0.73193pt\hbox{$$
\mbox{\scriptsize$\sim$ }$$}2.5\,\mathrm{nN}$.

Concerning the first one, the typical duration of an identification experiment
shall not exceed $T\raise 0.73193pt\hbox{$$ \mbox{\scriptsize$\sim$ }$$}3$
hours, mostly to ensure noise stationarity. The system can be stimulated with
a series of sine-waves of constant duration each $\delta t\simeq
1200\,\mathrm{s}$, as already described in Section 4.3. To simplify the
problem, the duration is kept fixed during the optimization. Furthermore, the
requirement of avoiding possible system transients at the beginning and the
end of each cycle, suggests to set null Dirichlet boundary conditions (i.e.,
null initial and final values of the signals) and leave gaps of $\delta
t_{\text{gap}}\simeq 150\,\mathrm{s}$. The general expression of a guidance
signal is a windowed series of sines

$o_{\text{i}}(t)=\sum_{n=1}^{N_{\text{inj}}}a_{n}\,\sin(2\pi
f_{n}\,t)\,\theta(t-t^{\prime}_{n})\,\theta(t^{\prime\prime}_{n}-t)\leavevmode\nobreak\
,$ (5.4)

where $\theta$ is the Heaviside unit-step, $f_{n}=n/\delta t$ is the injected
frequency of the $n$-th cycle and $a_{n}$ the corresponding amplitude, through
the maximum number of injections $N_{\text{inj}}=7$, and
$t^{\prime}_{n}=t_{0}+(n-1)(\delta t+\delta t_{\text{gap}})$ and
$t^{\prime\prime}_{n}=t_{0}+n\,(\delta t+\delta t_{\text{gap}})-\delta
t_{\text{gap}}$ the initial and final instants of the $n$-th injection cycle,
with $t_{0}$ the starting instant of the experiment. Clearly, the frequencies
$f_{n}$ are set by general requirements on the experiment duration, whereas