# Full-scale modal testing of a Hawk T1A aircraft for benchmarking vibration-
based methods

Marcus Haywood-Alexander111The work published here was undertaken at the
University of Sheffield, but this author is now affiliated at ETH Zurich.
<EMAIL_ADDRESS>Robin S. Mills<EMAIL_ADDRESS>Max D. Champneys
<EMAIL_ADDRESS>Matthew R. Jones
<EMAIL_ADDRESS>Matthew S. Bonney<EMAIL_ADDRESS>David Wagg<EMAIL_ADDRESS>Timothy J. Rogers
<EMAIL_ADDRESS>
###### Abstract

Research developments for structural dynamics in the fields of design, system
identification and structural health monitoring (SHM) have dramatically
expanded the bounds of what can be learned from measured vibration data.
However, significant challenges remain in the tasks of identification,
prediction and evaluation of full-scale structures. A significant aid in the
roadmap to the application of cutting-edge methods to the demands of in-
service engineering structures, is the development of comprehensive benchmark
datasets. With the aim of developing a useful and worthwhile benchmark dataset
for structural dynamics, an extensive testing campaign is presented here. This
recent campaign was performed on a decommissioned BAE system Hawk T1A aircraft
at the Laboratory for Verification and Validation (LVV) in Sheffield. The aim
of this paper is to present the dataset, providing details on the structure,
experimental design, and data acquired. The collected data is made freely and
openly available with the intention that it serve as a benchmark dataset for
challenges in full-scale structural dynamics. Here, the details pertaining to
two test phases (frequency and time domain) are presented. So as to ensure
that the presented dataset is able to function as a benchmark, some baseline-
level results are additionally presented for the tasks of identification and
prediction, using standard approaches. It is envisaged that advanced
methodologies will demonstrate superiority by favourable comparison with the
results presented here. Finally, some dataset-specific challenges are
described, with a view to form a hierarchy of tasks and frame discussion over
their relative difficulty.

###### keywords:

Modal testing , Hawk T1A , Full-scale dynamic testing , Shaker test , SIMO
testing , Aircraft dynamics

††journal: Journal of Sound Vibration

[uos_add]organization=Dynamics Research Group, University of
Sheffield,addressline=Mappin Street, city=Sheffield, postcode=S1 4ET,
country=UK

[eth_add]organization=Chair of Structural Mechanics and Monitoring, ETH
Zurich,addressline=Wolfgang-Pauli-Strasse, city=Zurich, postcode=8049,
country=Switzerland

[ati_add]organization=The Alan Turing Institute,city=London, postcode=NW1 2DB,
country=UK

##  1 Introduction

How a structure behaves under dynamic loading is an imperative part of
engineering design, operation, and maintenance. In both experimental and
operational contexts, the measuring and monitoring of modal characteristics
can provide vital information on the structure [1, 2]. These modal
characteristics can be used to, either directly or indirectly, infer the
physical properties of the system. One such use of this information is to
determine the existence of damage in important locations on the structure.
When damage, such as cracks or disbonds, occur on the structure, this
corresponds to a localised loss of stiffness, which manifests as a change in
the eigenfrequencies of the system. However, structural nonlinearities can
also cause apparent changes to modal properties therefore causing false-
positive or false-negative detection [3].

Using a well-defined model, it is possible to generate estimates of physical
parameters using analytical, or purely-physics modelling (derived from first-
principles analysis of elastic deformation theory). The presence of
nonlinearity significantly complicates this task. The modelling of _known_
nonlinearities is in itself a difficult task [4, 5], however, when the
nonlinearities are unknown, or ill-prescribed, in the model, this presents an
even more difficult challenge. This, along with the challenge of learning
model parameters, has led to the increased development of data-driven
techniques to estimate or learn the models and parameters [6, 7, 8]. As these
algorithms are developed, it is important that they are consistently assessed
on meaningful examples, or benchmarks. The meaningfulness, as we define it
here, of such benchmarks can be described by two characteristics; they present
a reasonable challenge in that they that cannot be trivially solved by simple,
well-established approaches, and that they provide data from a variety of
loading or operating conditions.

In this work, we present a dataset with the aim to fulfil the requirements of
a meaningful benchmark for structural modal analysis. In order to do so,
experimental acceleration data was collected from a BAE system Hawk T1A
aircraft undergoing a variety of dynamic input excitation, with simulated
damage scenarios. The aim of the dataset is to provide a challenging benchmark
for structural dynamics, in fields such as system identification and
structural health monitoring. As well as presenting the data, an ancillary
contribution of this paper is to also present a baseline standard in the tasks
of system identification and structural health monitoring. The techniques
presented here do not necessarily represent the state of the art in these two
tasks but instead offer a first analysis; indicative of approaches that might
be employed before more advanced techniques are sought. The rationale behind
these choices is to establish how much insight can be gained from basic
techniques and to act as a baseline against which future approaches can be
assessed. A further motivation is to qualitatively establish the difficulty of
selected system identification and structural health monitoring tasks on the
Hawk data. This is an important undertaking as it highlights the issues
presented when advanced machine learning (ML) methods are applied to case
studies that can be solved by more simplistic, robust, and interpretable
methods.

The Hawk dataset is freely provided with the motivation that it might be used
to develop new algorithms for structural dynamics. The data is available in
open _.hd5_ format, and access is provided via a simple Python API, which is
also provided.222All data and documentation for the API is available at
https://doi.org/10.15131/shef.data.22710040.v1.

The remainder of this paper is arranged as follows. In Section 2, details of
the aircraft structure are given, showing the boundary conditions and general
information. This is followed by Section 3, where details on the experiment
and data collection are provided, including information on the full
experimental regime, data collected, control system, hardware, and
acquisition. At this point the paper moves on to the baseline tasks applied to
the dataset. Section 4 presents results of a series of system identification
tasks applied to the dataset, both in the frequency and time domain. Following
this, results from a structural health monitoring perspective are given in
Section 5. Finally, in Section 6 some observations are provided on the
challenges of the dataset, in order to showcase the use and meaningfulness of
the data.

##  2 Structure

The structure used in this dataset is a BAE systems Hawk T1A aircraft [9].
This is a real, decommissioned aircraft, that was used for advanced training
in the RAF. The aircraft was tested at the Laboratory for Verification and
Validation (LVV) in Sheffield, England. For this dataset, only the starboard
wing was measured, as well as supplementary measurements at the top of the
aircraft and the exhaust, more details on the exact locations of the sensors
will be given later in Section 3.2, Figure 3. The aircraft was positioned
resting on its wheels on the ground. Therefore, between the main body and the
ground are the wheels, supports and hydraulic shock absorbers. All these
substructures are joined with a variety of fixing methods and types, which are
not detailed here. Figure 1 shows both an image of the real aircraft, as well
as a schematic, indicating the boundary conditions of the wing.

(a) Image of aircraft.

(b) Diagram showing front view of aircraft.

Figure 1: Images of BAE systems Hawk T1A aircraft.

##  3 Experiment

The aim of this work was to generate modal analysis data for the starboard
wing of the Hawk aircraft. To do so, the wing was excited with a shaker, and
the acceleration response over the upper and lower surfaces of the wing were
measured using accelerometers, and processed to return the response in the
frequency domain. This section will detail all the experiment details,
beginning with the testing regime, the experiment setup and hardware, and then
acquisition and processing details.

###  3.1 Regime outline

The total experimental regime was split into two phases; open-loop control
using off-the-shelf hardware, and a custom-built, closed-loop control system.
Specific details on the control system and hardware will be given in Section
3.2. Each of these experimental regimes followed a hierarchical structure
based on the aims of the experiment; this is shown in Figure 2. The lowest
level of this hierarchy is herein referred to as a ‘test’, and each test was
performed with 10 repetitions. For clarification; a repetition is where the
run is finished, allowed to reset, and then repeated. This repeat style was
done to allow for quantification and mitigation of experimental error.

Figure 2: Testing regime for dataset. For the damaged tests, the acronyms
represent the locations of the added mass, where TLE, RLE, and CTE stand for
_tip-leading-edge_ , _root-leading-edge_ , and _central-trailing-edge_ ,
respectively. The weight numbers are ID-representations, the value of mass is
provided in the experiment details.

####  3.1.1 Phase 1: Open-loop control

In the first phase, there were two main sub-regimes; the nominally-undamaged
state, and the pseudo-damaged state. In the first case, the aircraft was left
in its assumed healthy (free from significant structural damage) condition,
and in the second case, damage was ‘simulated’ by adding masses of different
weights in specified locations. For the undamaged tests, the method of forcing
input was set to either burst random, or forward/reverse sine sweeps in the
linear/logarithmic scale. The burst random uses a white noise forcing input,
the amplitude of which is ramped at the beginning and end of the burst. For
the sine sweeps, the forcing input is a sine wave signal of which the
frequency increases during the time period of the measurement. Forward/reverse
refers to the direction of the frequency sweep, and linear/logarithmic refers
to the pattern of the sweep. This means that for forward logarithmic sweeps,
more time is spent exciting the lower frequencies, and for reverse logarithmic
sweeps, more time is spent exciting the higher frequencies. For each of the
excitation types, the amplitude of the forcing was increased by increasing the
output voltage in increments of 0.2V.

For the damaged cases in Phase 1, only burst random excitation was used, and
to simulate a range of damage cases, the weight and location of the added
masses was changed. In Figure 2, TLE, RLE, and CTE refer to these locations,
and stand for _tip-leading-edge_ , _root-leading-edge_ , and _central-
trailing-edge_ , respectively. The locations on the wing are shown in Figure
3, and exact details are given in the accompanying technical report [10]. Much
like the undamaged tests, the amplitude of the forcing was increased by
increasing the voltage, but this time in increments of 0.4V.

####  3.1.2 Phase 2: Closed-loop control

The second phase of testing is split into sub-regimes based on the target
signal types, by iteratively generating the signal passed to the shaker. More
details on the control algorithm are given in Section 3.2.1, along with the
definition of each of the actuation signal types. To determine information on
the effect of randomness, each of these signal types are tested with different
initial seeds. To increase the amplitude of the forcing, as this is a closed-
loop control which uses the measured amplitude, changing the voltage output
will simply cause the control algorithm to decrease the signal level.
Therefore, instead the amplitude of the target spectra were increased to
corresponding levels. It is important to note, these do not necessary
correspond exactly to the levels used in Phase 1, they are instead values
which return amplitudes of a similar spectral density.

###  3.2 Setup and hardware

To measure the structural response, PCB Piezotronics accelerometers were used
with a nominal sensitivity of 10 mV/g. When converting voltage readings to
acceleration values, the sensor-specific sensitivities, accurate to 0.01 mV/g,
were used. A total of 54 sensors were placed on the wing, which are placed in
three lines on the upper surface, and two lines on the lower surface of the
wing, where a line runs from the root to the tip. As well as these 5 strips,
additional sensors were placed to provide information at boundary condition
locations of the root and the open panel for the landing wheel. The force was
applied to the lower surface of the wing using a Tira TV 51140-MOSP shaker,
attached directly to a PCB Piezotronics 208C02 force transducer with a
sensitivity of 11.11 mV/N.

A diagram showing the locations of the sensors, forcing, and masses is given
in Figure 3. The locations of the sensors in this diagram are the nominal
locations, but may have been slightly shifted in order to avoid components of
the wing, exact coordinates are provided in the data and the technical report
[10]. The black sensor locations indicate those placed for the purposes of
experimental modal analysis of the general structure. The blue sensor
locations indicate those placed for the collection of data relating to the
boundary conditions. Also shown are the shaker attachment location, indicating
where the excitation force is applied, and the locations of the added masses
and their corresponding ID tags.

Figure 3: Sensor layout on the wing, also indicating the locations of the
shaker and additional masses.

For the open-loop control, a Siemens PLM LMS SCADAS unit was used to generate
the forcing signal, and simultaneously measure the acceleration responses. For
the closed-loop control system, National Instruments hardware was used, along
with LabView, to develop a frequency-domain control system, targeting the
force signal at the location of the shaker. A table of the full hardware is
shown in Table 1.

Name | Quantity | Use
---|---|---
NI-cDAQ-9179 | 2 | Chassis for total system
NI-9232 | 20 | Measures voltage signals from accelerometers
NI-9260 | 1 | Outputs voltage signal
Table 1: Hardware components for closed-loop control system.

####  3.2.1 Control algorithm

In the second phase of testing, a custom algorithm was implemented to control
the actuation signal in the frequency domain. This allows the use of pre-
prescribed forcing types, for both amplitude and phase of the frequencies. The
overall control flow is shown in Figure 4, showing the process which covers
one repeat of a test. The general aim of the algorithm is to control the
output of the shaker to the structure in the frequency domain. For this work,
the frequency domain information is defined and calculated as the power-
spectral-density (PSD), and contains both magnitude and phase information. The
target PSD of this force is denoted $S_{r}(\omega)$, and the measured PSD from
the transducer is denoted as $S_{f}(\omega)$. The algorithm estimates a model
of the mapping between the driving signal, in terms of voltage, to the applied
force. This is achieved by assuming a linear transfer function,

$|F(\omega)|=|H(\omega)||S_{d}(\omega)|$ (1)

where $F(\omega)$ is the force applied to the structure, $H(\omega)$ is the
transfer function, and $S_{d}(\omega)$ is the drive signal in volts, all in
the frequency domain as functions of the circular frequency $\omega$.

The target actuation $S_{r}(\omega)$ is defined in terms of the PSD of the
force applied. Details of how these target signals are generated are provided
in Section 3.2.2. At the first step of the control algorithm, the drive PSD
$S_{d}(\omega)$ is simply set to the reference (i.e. $|H|=1$). At each
iteration, the drive signal $s_{x}(t)$, is calculated using an inverse Fourier
transform,

$s_{d}(t)=\mathcal{F}^{-1}(|S_{d}(\omega)|)$ (2)

When the actuation signals were generated, a linear amplitude ramp was applied
to window the signal in the time domain, to reduce shock to the structure and
equipment. This signal is sent out from the control system, and the PSD of the
force, $S_{f}(\omega)$ is measured. This measurement is then used to update
the estimate of the linear transfer function using,

$|H(\omega)|=\alpha\sqrt{\frac{S_{r}(\omega)}{S_{f}(\omega)}}$ (3)

where $0.0\geq\alpha\leq 1.0$ is a proportional weighting factor which is
included to slowly iterate the estimated transfer function to mitigate
potential out-of-control issues, or overloading of the system. The value of
this is chosen by the control system designer, and in this case was chosen to
be 0.8. Then, the drive PSD is updated, and subsequently the drive signal
using Equation 2,

$|S_{d}(\omega)|=\frac{|S_{r}(\omega)|}{|H(\omega)|}$ (4)

During the drive updating, the target PSD error is defined as the mean squared
error,

$\varepsilon=\frac{1}{N}\sum_{i=1}^{N}\left(|S_{f}(\omega^{(i)})|-|S_{r}(\omega^{(i)})|\right)^{2}$
(5)

where $N$ is the number of frequency points in the discretised frequency
domain. Once the error is below a specified target, $\eta$, the latest drive
PSD is used for the acquisition stage. For this work, the target error was set
to 0.01.

Figure 4: Control algorithm process flow, indicating the total workflow for
one _repeat_ of a test.

####  3.2.2 Frequency Domain Actuation Types

This subsection aims to explain the difference between the actuation types in
Phase 2, which are defined in the frequency domain. Both actuation types
generate the time domain signal as a sum of the discrete frequencies
represented by their definitions in terms of power-spectral-density. The
actuation is defined using a vector of frequencies, $\boldsymbol{\omega}$,
with corresponding vectors of amplitude, $\mathbf{A}$, and phase,
$\boldsymbol{\phi}$. Where $N$ is the length of these vectors, and is the same
as the PSD definitions in the above subsection, the discrete time signal,
$s(t)$ is defined as as scalar function of the discrete time points,

$s(t)=\sum_{i=1}^{N}A_{i}\sin(\omega_{i}t+\phi_{i})$ (6)

The vector of frequencies is chosen based on the bandwidth of modes of
interest, and its discretisation is aligned with the acquisition parameters.
Generally, the vector $\mathbf{A}$ is designed to ‘ramp up and ramp down’ the
amplitude with respect to the frequencies. This is done in order to maintain
periodicity in the time domain signal to avoid shock in the physical system.
The phase is generated by randomly sampling from a uniform distribution
between 0 and $2\pi$. Using these generation methods with no modification
forms the random phase actuation type. However, for the odd phase actuation,
every _other_ element in the vectors $\mathbf{A}$ and $\boldsymbol{\phi}$ are
set to zero,

$A_{i}=0,~{}~{}~{}~{}\phi_{i}=0,~{}~{}\text{for}~{}~{}i=2n-1~{}~{}\text{for}~{}~{}0\leq
n\leq N/2$ (7)

Figure 5: Examples of the target driver PSDs for (left) random phase and
(right) odd phase multisine types. The top plots show the target PSD power in
terms of dB, and the bottom plots show the corresponding phase.

###  3.3 Acquisition and Processing

For all the tests, a number of acquisition parameters (listed in Table 2) were
kept constant. The number of samples indicates the number of sampling points
taken in the time domain. The burst/sweep time represents the portion of the
measured signal in which the random burst or the sine sweep takes place; i.e.
for 90% burst time, the first and last 5% of the acquisition time has no
actuation.

Phase 1 | Phase 2
---|---
Parameter | Value | Parameter | Value
Sample Rate | 512 Hz | Sample Rate | 2048 Hz
Number of samples | 16,384 | Number of samples | 32,768
Burst/sweep time | 90% | Burst time | 90%
Amplifier level | 50% | Amplifier level | 50%
Repetitions | 10 | Repetitions | 10
Table 2: Acquisition parameters for each test phase.

The data is then saved with identifiers which relate to the test series as
shown in Figure 2. The sub-regime identifiers are shown in Table 3, where the
placeholder ‘X’ in the identifiers are the test numbers identifiers. These are
integer values which correspond to the voltage level, added mass weight,
and/or seed used. Further details of the exact labels for each individual test
can be found in the technical report [10].

Phase 1 | Phase 2
---|---
Undamaged | Damaged | Random Phase | Odd Phase Multisine
BR_AR_X | DS_TLE_X | RPH_AR_X | ORP_VS_X
FSS_AR_X | DS_RLE_X | FLS_AR_X |
RSS_AR_X | DS_CTE_X | RLS_AR_X |
FLS_AR_X |  |  |
RLS_AR_X |  |  |
Table 3: Test identifiers used for naming data in acquisition and processing.

###  3.4 Dataset structure

In order to keep this paper self-contained, a brief description of the dataset
structure is included here. All data is stored in the open-source
_hierarchical data format_ (hd5). The hd5 format can be accessed using several
open-source tools, but the authors also provide a python interface.333A python
API for the hawk data is available at https://github.com/MDCHAMP/hawk-data The
hd5 format essentially mimics the tree-like file structure of modern
computers. Data can be accessed by specifying the path to the desired sensor
channel and data series. The general form of the path for the Hawk SBW test
data is:

`/SBW/{test phase}/{test series}/{repetition}/{sensor}/{signal}`

where SBW (starboard wing test) is the overall test campaign the bracketed
quantities refer to the following optional fields:

  * 1.

`test phase` $\in$ `{LMS, NI}`: Refers to the frequency domain and time domain
data respectively, as per Table 2 and Figure 2.

  * 2.

`test series` $\in$ `{BR_AR, DS_TLE, ...}` (Phase 1) and $\in$ `{RPH_AR,
ORP_VS}` (Phase 2): Test series identifier (see Figure 2 for definitions).

  * 3.

`repetition` $\in$ `{01, 02, ...}`: Test repetition number. For each test
series there are nominally 10 repetitions, although for some tests fewer
repetitions are available.

  * 4.

`sensor` $\in$ `{FRC, LLE-01, ..., LLC-07}`: Sensor location (see Figure 3 for
definitions).

  * 5.

`signal` $\in$ `{(coherenceSpectrum, spectra, frf}` (Phase 1) and $\in$ `{acc,
(force)}` (Phase 2): Desired signal. Note that frequency-response-functions
(FRFs) are only available for output channels in the phase 1 tests.

The self-describing nature of the .hd5 format allows metadata to be stored
alongside the collected series. Test data related to sensor position,
sensitivities, units and test parameters (added mass, excitation level etc.)
can all be accessed by querying the metadata fields at each level of the .hd5
hierarchy. Additional detail regarding the dataset structure and interface can
be found at https://doi.org/10.15131/shef.data.22710040.v1.

##  4 System identification

With the description of the experimental campaign complete, attention can now
be turned to establishing a benchmark for the Hawk data in two common tasks in
structural dynamics analysis _system identification_ and _structural health
monitoring_. It is envisaged that the results here will act as a baseline
level against which more advanced methods might be compared in future work.

One of the key challenges in structural dynamics is obtaining accurate models
of system dynamics. For structures that can be considered predominantly
linear, the default framework for modelling dynamics is that of _modal
analysis_ [11]. Modal analysis is familiar to many engineers as a method of
decomposing system dynamics into a number of eigenfrequencies and
eigenfunctions that encode the response of the system to any input. These
properties are invariant to any excitation and can be conveniently expressed
in terms of a discrete number of interpretable properties; the natural
frequencies $\omega_{r}$, damping ratios $\zeta_{r}$ and modal matrix
$\Phi$.444In this section, the $r$th vector of the modal matrix (commonly
known as the $r$th modeshape) will be denoted by $\bm{\phi}^{(r)}$.

Experimental modal analysis (EMA) [11] is a general term for methods that
attempt to specify the modal properties from measured structural dynamics
data. Broadly, these methods can be divided into methods that operate in the
time and frequency domains. In this paper, two EMA methods are presented on
the proposed Hawk dataset. In the frequency domain, the method of rational
fractional polynomials (RFP) [12] is employed. In the time domain, an output-
only stochastic subspace identification (SSI) is shown.

###  4.1 Modal analysis: Frequency domain

Frequency domain methods for EMA rely on approximating the form of the
frequency response function (FRF), often employing techniques from curve
fitting. Assuming that the dynamics are linear, the FRF is a linear transfer
function from the Fourier-transformed input force (N) to the output
acceleration (g) that has a generalised form that depends only on the modal
properties. In order to keep the notation in this section general (and in
keeping with the system identification literature), the input forcing spectra
will be denoted $X(\bm{\omega})$ and the output acceleration $Y(\bm{\omega})$,
where $\bm{\omega}$ is the vector of frequencies for which data were recorded.
The accelerance FRF can be expressed as,

$H(\bm{\omega})=\frac{Y(\bm{\omega})}{X(\bm{\omega})}$ (8)

The expression for a general accelerance FRF in terms of the modal parameters
is,

$H_{il}(\bm{\omega})=\sum_{r=1}^{n}\frac{-\bm{\omega}^{2}\bm{\phi}_{i}^{(r)}\bm{\phi}_{j}^{(r)}}{-\bm{\omega}^{2}+2j\zeta_{r}\omega_{r}\bm{\omega}+\omega_{r}^{2}}$
(9)

Where $\omega_{r}$, $\zeta_{r}$ and $\phi^{(r)}$ are the natural frequency,
damping ratio and modeshape corresponding to the $r$th vibration mode and $i$
and $l$ are the indices of the input and output location respectively.

The method of RFP works by approximating (9) as a ratio of two polynomials.
Considering here a single-input, single output (SISO) system comprised of $n$
modes, one has,

$H(\bm{\omega})=\frac{\sum_{k=0}^{m}a_{k}(j\bm{\omega})^{k}}{\sum_{k=0}^{n}b_{k}(j\bm{\omega})^{k}}$
(10)

where the coefficients $a_{k}$ and $b_{k}$ depend on the modal properties.
Note that the number of terms in the numerator polynomial does not necessarily
have to be the same as the denominator. Additional terms in the numerator
polynomial can be used to compensate for out-of-band effects (i.e modes from
outside the frequency range of interest).

RFP learns the $a_{k}$ and $b_{k}$ coefficients from the measured $H$ and then
infers the modal properties. Rather than learn the coefficients by
optimisation, it is preferable to construct a least squares scheme. In order
for the coefficients to be identifiable we must include a constraint that
fixes one of the polynomial coefficients (or else there would be a constant
factor invariance). Here, after [12], the lead coefficient of the denominator
is arbitrarily set to unity $b_{n}=1$. In order to simplify notation, the
vectors $\bm{a},\bm{b}$ are defined as the unknown coefficients of the
numerator (up to order $m$) and the denominator (up to order $n-1$)
polynomials respectively.

Representing the rational FRF form above in matrix form,

$H=(\Phi_{\bm{a}}\bm{a})(b_{n}(j\bm{\omega})^{n}+\Phi_{\bm{b}}\bm{b})^{-1}$
(11)

where the matrices $\Phi_{\bm{a}}$ and $\Phi_{\bm{b}}$ are polynomial bases of
order $n$ and $m$ respectively for the vector of frequency lines
$\bm{\omega}$, and $(j\bm{\omega})^{n}$ is the vector $j\bm{\omega}$ raised to
the power $n$ element-wise. In order to construct the least squares scheme,
the error between the measured $\hat{H}(\bm{\omega})$ and predicted FRFs can
be expressed as,

$\bm{e}=\Phi_{\bm{a}}\bm{a}-((j\bm{\omega})^{n}+\bm{b}\Phi_{\bm{b}})\hat{H}(\bm{\omega})$
(12)

In order to simplify the forthcoming notation (and for compatibility with the
notation of [12]), the following notation is introduced here,

$P=\Phi_{\bm{a}}$ (13)

$T=\Phi_{\bm{b}}\hat{H}(\bm{\omega})$ (14)

$\bm{w}=b_{n}(j\bm{\omega})^{n}\hat{H}(\bm{\omega})$ (15)

Thus the objective function becomes,

$J(\bm{a},\bm{b}|\hat{H})=\bm{e}^{\dagger}\bm{e}=(P\bm{a}-(\bm{w}+T\bm{b}))^{\dagger}(P\bm{a}-(\bm{w}+T\bm{b}))$
(16)

where ${\dagger}$ represents the Hermitian transpose. The above is a positive
quadratic in both $\bm{a}$ and $\bm{b}$ and therefore has a single optimum
($\bm{a}^{*}$, $\bm{b}^{*}$), obtained exactly by setting the partial
derivatives in $\bm{a}$ and $\bm{b}$ to zero,

$\frac{\partial
J(\bm{a},\bm{b}|\hat{H})}{\partial\bm{a}}\bigg{|}_{\bm{a}^{*},\ \bm{b}^{*}}=0$
(17)

and,

$\frac{\partial
J(\bm{a},\bm{b}|\hat{H})}{\partial\bm{b}}\bigg{|}_{\bm{a}^{*},\ \bm{b}^{*}}=0$
(18)

Now evaluating these derivatives starting with the former condition,

$\frac{\partial
J}{\partial\bm{a}}=2\mathbb{Re}(P^{\dagger}P)\bm{a}-2\mathbb{Re}(P^{\dagger}T)\bm{b}-2\mathbb{Re}(P^{\dagger}\bm{w})$
(19)

where $\mathbb{Re}$ denotes the real part of the complex vector. Next, for the
second condition,

$\frac{\partial
J}{\partial\bm{b}}=-2\mathbb{Re}(T^{\dagger}P)\bm{a}+2\mathbb{Re}(T^{\dagger}T)\bm{b}+2\mathbb{Re}(T^{\dagger}\bm{w})$
(20)

Finally, applying the condition of zero gradient at the optimum, one has the
block-linear system that can be solved for the polynomial coefficients.

$\mathbb{Re}\left(\begin{bmatrix}P^{\dagger}P&-P^{\dagger}T\\\
-T^{\dagger}P&T^{\dagger}T\end{bmatrix}\right)\begin{bmatrix}\bm{a}\\\
\bm{b}\end{bmatrix}=\mathbb{Re}\left(\begin{bmatrix}P^{\dagger}\bm{w}\\\
-T^{\dagger}\bm{w}\end{bmatrix}\right)$ (21)

It is worth noting that this expression differs in several respects from the
expression in [12] due to the presence of several errors in the original
manuscript. Finally the modal properties (natural frequencies and damping
ratios) can be recovered from the coefficients of the denominator polynomial
as,

$\omega_{r}=|R_{r}|$ (22)

$\zeta_{r}=\frac{-\text{Re}(R_{r})}{|R_{r}|}$ (23)

where $R_{r}$ is the $r$ th complex-conjugate root of the denominator
polynomial in (10). Extraction of the modeshapes is also possible, requiring
both the coefficients of the numerator and the denominator polynomials.
However it is not considered here.

In Figure 6, the FRF of the hawk at sensor LTC-05 (Lower trailing edge centre,
burst-random excitation, 0.4V LMS input voltage, averaged over 10 repeats) is
plotted. In order to demonstrate the RFP approach, some modal properties will
be extracted from this signal.

Figure 6: FRF from sensor LTC-05 subject to burst-random excitation at 0.4V.

As a first analysis, all frequency lines are included in the computation of
the modal properties with the number of modes set to 10 (corresponding to a
20th order polynomial). The resulting FRF is plotted in Figure 7 with dashed
lines corresponding to the natural frequencies. As can be seen, the higher
frequencies are well approximated but there is significant errors in the lower
frequencies and many of the dominant modes have been missed entirely. This is
due to an inherent bias toward higher frequencies in the RFP method that can
be attributed to the choice of the error criterion $\bm{e}$.555This bias can
be seen by noting that the choice of $\bm{e}$ does not correspond directly to
the difference between the measured and predicted FRF, but a difference that
is scaled by the frequency.

Figure 7: Measured vs. identified RFP FRF for the sensor at LTC-05, vertical
lines are the identified natural frequencies.

A more practical approach (indicative of how frequency-domain RFP is often
implemented) is to instead perform the curve fit on a number of pre-selected
bounds. In order to compensate for the effects of ‘out-of-bound’ modes,
additional terms can be added to the numerator polynomial. Here, for each
range, 6 additional terms are added to the numerator polynomial. The predicted
FRF and natural frequencies are plotted in Figure 8 for identification bounds
given in Figure 4. In the figure, the quality of the FRF fit is visibly
superior to the global fit and the the peak locations are well identified
within the selected bands. Note that the total FRF can be reconstructed by
first estimating the damping ratios and then using the general form of the
modal transfer function although this is not considered here.

It is clear that such an approach cannot be considered to be ‘automatic’ and
that the performance comes to depend strongly on the choice of the
identification ranges. However it is important to consider such approaches as
a baseline to which more advanced methods might be compared.

Table 4: Bands and identified modal parameters from the banded RFP identification. Lower bound (Hz) | Upper bound (Hz) | Modes in band |  $\omega_{n}$ (Hz) | $\zeta_{n}$
---|---|---|---|---
5 | 9 | 1 | 7.094 | 0.02434
12 | 14 | 1 | 13.05 | 0.01552
15 | 19 | 1 | 17.1 | 0.005918
22 | 24 | 2 | 22.66, 23.56 | 0.008452, 0.009352
26 | 30 | 1 | 27.35 | 0.01584
30.5 | 31 | 1 | 30.77 | 0.0004578
35 | 37 | 1 | 36.06 | 0.00005626
40.5 | 44 | 2 | 41.16, 42.85 | 0.008827, 0.007161
48.5 | 54 | 2 | 49.61, 52.07 | 0.00995, 0.00469
86 | 90 | 1 | 87.41 | 0.01138
92 | 100 | 1 | 92.54 | 0.02798
112 | 118.5 | 1 | 116.9 | 0.01307
119.5 | 122 | 1 | 120 | 0.003427
122 | 125 | 1 | 124.6 | 0.008256
135 | 138 | 1 | 124.6 | 0.008256
154 | 162 | 1 | 157.1 | 0.01505
Figure 8: Measured vs. identified banded RFP FRF for the sensor at LTC-05,
vertical lines are the identified natural frequencies.

###  4.2 Modal analysis: Time domain

Several approaches to EMA consider the identification problem in the time-
domain. Of these, a common approach is stochastic subspace identification
(SSI) [13]. The SSI approach outlined here has the added advantage that it
operates in an output-only setting (subject to a broadband excitation
assumption). This is especially attractive for applications whereby
measurement of system input in-service (such as in an aircraft during flight)
is challenging.

The collection of time-series data from the Hawk allows the effectiveness of
SSI to be tested. Here, a classical approach of “covariance-driven” SSI (Cov-
SSI) is used. In this work the terms SSI and Cov-SSI will now be used
interchangeably. A full background to SSI will not be presented in this paper
in the name of brevity, instead the reader is referred to [13, 14] for general
introductions to the family of methods or towards [15] in reference to the
output-only _operational modal analysis_ (OMA) setting or more recently with
respect to uncertainty quantification [16] (although the subject of
uncertainty is not addressed in this work, see section 6).

Briefly, the application of SSI requires formulation of the (block) Hankel
matrices of the “past” and “future”. These are constructed through stacking
time-delayed (lagged) versions of the measured time-series responses of the
system. The number of lags chosen in the Hankel matrices defines the maximum
order of the system which can be identified, where this maximum order is twice
the maximum number of modes. After constructing these matrices, $Y_{p}$ and
$Y_{f}$, for the past and future respectively, the cross covariance can be
computed for $N$ datapoints as,

$\Sigma=\begin{bmatrix}\Sigma_{pp}&\Sigma_{pf}\\\ \Sigma_{fp}&\Sigma_{ff}\\\
\end{bmatrix}=\frac{1}{N}\begin{bmatrix}Y_{p}\\\
Y_{f}\end{bmatrix}\begin{bmatrix}Y_{p}^{\mkern-1.5mu\mathsf{T}}&Y_{f}^{\mkern-1.5mu\mathsf{T}}\end{bmatrix}$
(24)

This cross covariance is partitioned to consider the auto-covariance within
the $Y_{p}$ and $Y_{f}$, and the cross covariance between them. To perform
SSI, a canonical correlation analysis between these two sets of variables is
performed based on the covariances. Mathematically, a singular value
decomposition (SVD) is taken of the normalised covariance,

$USV^{\mkern-1.5mu\mathsf{T}}=\Sigma_{pp}^{1/2}\Sigma_{pf}\Sigma_{ff}^{1/2}$

The singular values (along the diagonal of $S$) are sorted in descending order
and, to select an appropriate model order, a number of the smallest of these
values is discarded. If one wishes to construct the consistency
(stabilisation) diagram, sequentially more and more values are discarded such
that the model order is reducing and identification is performed at each of
these orders. Since the computationally intensive stage of SSI is the
formation of the Hankel matrices, covariances and SVD, this procedure for
generating the consistency diagram is quite efficient. The truncated set of
singular values can be denoted as $\hat{S}$, the associated columns of $U$ and
$V$ are also removed with these truncated matrices as $\hat{U}$, $\hat{V}$.
Given that it can be shown that the cross covariance is the product of the
observability and controllability matrices,
$\Sigma_{fp}=\mathcal{O}\mathcal{C}$, these can be recovered by:

$\mathcal{O}=\Sigma_{ff}^{1/2}\hat{U}\hat{S}^{1/2}$ (25a)
$\mathcal{C}^{\mkern-1.5mu\mathsf{T}}=\Sigma_{pp}^{1/2}\hat{V}\hat{S}^{1/2}$
(25b)

From $\mathcal{O}$ the system matrices $A$ and $C$ of a linear dynamic state
space model with states $\vec{x}$:

$\displaystyle\vec{x}_{t+1}=A\vec{x}_{t}+\vec{v}_{t}$ (26)
$\displaystyle\vec{y}_{t}=C\vec{x}_{t}+\vec{w}_{t}$

where $\vec{v}_{t}$ and $\vec{w}_{t}$ are samples from multivariate Gaussian
distributions with covariances $Q$ and $R$ respectively. Given that the
observability matrix is constructed as,

$\mathcal{O}=\begin{bmatrix}C\\\ CA\\\ CA^{2}\\\ \vdots\\\
CA^{K}\end{bmatrix}$

where $K$ is the maximum model order considered, the $C$ matrix of the system
can be recovered as the first block row of this matrix and then the $A$ matrix
can be recovered via a least squares solution. The $A$ and $C$ matricies fully
define the dynamics of the system and it is therefore possible to recover the
modal properties. A common step in the SSI approach is to vary the size of the
$A$ and $C$ matrices of the system and hence the model order by reducing the
level of truncation of the singular values (usually in steps of two since the
underlying system is assumed to be of a set of second order modes). Then,
comparing the recovered modal properties a consistency666Here the term
_consistency_ is preferred to _stabilisation_ since there is not any link
between the pole stability in terms of a requirement for the system to possess
positive damping and the similarity of the modal properties between each model
order. Therefore, the terminology could be ambiguous and give some implication
about the system which may be misinterpreted. (also “stabilisation”) diagram
is constructed based on the changes in these modal properties as the model
order varies.

In Figure 9, the time series from all accelerometers in the random phase
amplitude test (excitation level 0.4V) are plotted. Each burst signal
represents 32,768 datapoints captured over a 16s second interval. Unlike the
RFP approach above, SSI is able to simultaneously incorporate data from all of
the sensors at once. In order to speed up the computation of the SSI, the time
series are decimated by a factor of 5 prior to the calculation of the Hankel
matrices.

Figure 9: Time series collected at each of the 55 sensor locations subject to
closed-loop random phase excitation at 0.4V.

In order to evaluate the performance of the SSI, the singular-valued spectrum
(SVS) [17] is calculated. It is useful to consider the SVS as a convenient
visualisation as it is expected that all modes should be evident in only the
first few singular values, whereas any given sensor may not have all modes
visible (for example if a sensor is close to a nodal point in one of the
modeshapes). In Figure 10, the modal consistency diagram is plotted over the
SVS. This diagram plots the predicted natural frequencies at each model order
up to a maximum of 50. A standard approach in EMA is to interpret the ‘true’
modes as those that are ‘consistent’ in the diagram (that is they are present
at each model order).

As can be seen in the figure, the SSI approach has identified the major
resonances well although there are a significant number of non-physical modes
especially in the higher order-models. Once again it useful to interpret this
result as a baseline against which more advanced output-only time domain
methods might be considered.

Figure 10: SSI consistency diagram (black dots) compared to the first 5
components of the SVS (lines), subject to closed-loop random phase excitation
at 0.4V.

##  5 Structural health monitoring

To equip engineering asset management teams with the information to make more
optimal operation, maintenance and repair decisions for their high value
infrastructure, the last several decades have seen a rise in developments of
structural health monitoring strategies to provide insight into the condition
of structures using remote sensing [18]. A health monitoring system is defined
as the implementation of a strategy that facilities the real time assessment
of the condition of a structure - often for civil infrastructure or high value
assets such as aircraft - for ideally the span of the structure’s operational
life. Broadly speaking, installing an SHM strategy first involves acquiring
data relating to the dynamic behaviour of the structure, such as accelerometer
measurements. Features sensitive to damage are then extracted, such as the
natural frequencies of the structure, which are subsequently used as inputs to
a statistical model that relate the features to some information pertaining to
the damage state of the system.

The necessary complexity and sophistication of the model used can vary
significantly, and is generally dependent on which level of the SHM hierarchy
one is interested, of which there are five stages [19]:

  1. 1.

Detection - is damage present?

  2. 2.

Localisation - where is the damage?

  3. 3.

Classification - what type of damage?

  4. 4.

Assessment - how severe is the damage?

  5. 5.

Prognosis - what is the remaining useful life of the structure?

As progress is made up the hierarchy, the complexity of the task increases,
with classification, assessment and prognosis being particularly challenging.
Regardless of the complexity of the SHM task that is being attempted, to
transition monitoring algorithms from research to operational implementation,
it is vital that the algorithms are validated on real world data to ensure
that they are robust and reliable such that the risk of malfunction is inline
with the guidelines of the relevant industry. Datasets that can be used to
validate such methods are therefore critical in progressing structural health
monitoring into engineering practice.

In the following section, an overview of some basic SHM approaches that have
been applied to the dataset introduced in this paper will be provided,
including the use of both unsupervised and supervised machine learning
techniques. In the analysis that follows, data are taken from the damage
simulation tests in the frequency domain. In total, 27 test runs are taken at
amplitude levels between 0.4V and 2.0v for three levels of damage severity
(added mass) and at three different locations along the wing. For undamaged
condition data, the tests are taken from the burst-random amplitude ramp test
series with excitation values between 0.4V and 2.0v. In total this corresponds
to 32 test series and 320 FRFs at each sensor.

In this paper, two SHM approaches are presented. As a supervised learning
benchmark, a Gaussian mixture model (GMM) is fitted to the principal
components of the natural frequencies. In the unsupervised case, a novelty
detector based on the Mahalanobis distance is employed with damage sensitive
spectral lines used as features.

###  5.1 Damage detection: Supervised learning

In order to demonstrate a first analysis in supervised SHM, a GMM is fitted to
features arising from a modal identification as per the previous section.
Natural frequencies are extracted using an RFP approach with ranges defined in
Table 4. Overall, 18 natural frequencies are derived for each sensor for each
of the test runs. In order to reduce the variance in the identification
process, the natural frequencies are averaged over each sensor location.
Figure 11 plots a boxplot of the resulting features. As can be seen in the
figure, there are several outliers in the identification of the 10th peak,
indicating that this peak has not been well identified. To avoid the variance
in this mode dominating the features, the 10th natural frequency is removed
from the feature set for an overall feature dimension of 18.

Figure 11: Boxplot of the natural frequencies from the RFP analysis for each
of the testing scenarios.

In order to further reduce the feature dimension, a principal component
analysis (PCA) decomposition is taken on the natural frequencies. The first
two principal components of the resulting features are plotted in Figure 12.
In the plot the samples are coloured by damage extent and the relative sizes
of the markers represent the excitation level. Considering the data in Figure
12, it is clear that the classes are not linearly separable and therefore the
GMM classifier might not be expected to perform well. However it is still a
useful result to consider the GMM as baseline for clustering performance such
that advanced approaches might be compared to first analyses.

Figure 12: First two principal components of the natural frequency features.
In the plot colour denotes damage severity (added mass) and marker size
denotes excitation level (larger markers for higher excitation levels).

In order to classify the data, a GMM is trained. For training data, a subset
of half of the training examples is used with balanced membership of each
class in the training and testing sets. The training data are then used to
empirically fit Gaussian distributions to each of the classes by evaluating
the sample mean and covariances. The overall probability that an unseen
example $x^{*}$ has a label $y^{*}$ corresponding to the $k$th class can be
written,

$p(y^{*}|X)=\sum_{k}\pi_{k}\mathcal{N}(\bm{\mu}_{k},\bm{\Sigma}_{k})$ (27)

where $X$ are the training data and,

$\pi_{k}=\frac{n_{k}}{N}$ (28)

the proportion of training samples in the $k$th class and,

$\bm{\mu}_{k}=\mathbb{E}[\bm{x}_{i}]\quad\forall i\ \ \text{s.t.}\ \
\bm{y}_{i}=k$ (29)

$\bm{\Sigma}_{k}=\operatorname*{cov}[\bm{x}_{i}]\quad\forall i\ \ \text{s.t.}\
\ \bm{y}_{i}=k$ (30)

are the maximum likelihood estimators of the class mean and covariances. The
classification performance of the GMM in several classification tasks is
plotted in Figure 13, and confusion matrices are depicted in Figure 14. In the
figures, datum points used during training are marked with a grey cross. The
predicted label is denoted by the centre of each marker and the true label is
denoted by the colour of the border.

(a) Presence of damage.

(b) Added mass amount.

(c) Damage location.

(d) Excitation level.

Figure 13: Supervised GMM clustering results on the first two principal
components of the natural frequency features.

As can be seen in the figures, the performance of the classification is strong
in the simplest damage detection case (i.e. binary classification) but the
performance in other scenarios is poor. It appears that the simplistic model
and feature set used here is insufficiently rich to separate the classes.

(a) Presence of damage.

(b) Added mass amount.

(c) Damage location.

(d) Excitation level.

Figure 14: Confusion matrices from the supervised GMM clustering.

###  5.2 Damage detection: Unsupervised learning

In SHM tasks, it is common for data to lack labels pertaining to the damage
condition of the structure that the data are representative of, often leading
to SHM practitioners undertaking unsupervised learning approaches. In this
setting, it is common to perform outlier analysis, where incoming data from
the structure in an unknown condition are compared by some distance metric to
data generated by a baseline "normal" condition. If outliers are observed in
the incoming data, then it is believed that damage may have developed in the
structure 777Whether an outlier is truly representative of damage or not is
strongly dependent on the quality of the feature extraction, and how sensitive
to damage whilst insensitive to external influences such as operational and
environmental fluctuations the chosen feature is..

The principle of outlier detection is the computation of some distance metric
between a test data point $\mathbf{x}_{i}$ and the centre of the distribution
that we believe the normal condition data to be drawn from. Assuming that the
data in the baseline condition follow a normal distribution, a standard
distance or discordancy measure is that of the Mahalanobis distance, defined
as,

$D=(\mathbf{x}_{i}-\bm{\mu})^{T}\Sigma^{-1}(\mathbf{x}_{i}-\bm{\mu})$ (31)

where $\mu$ and $\Sigma$ are the mean and covariance of the normal condition
observations. Monte Carlo sampling can then be used to determine a threshold
for $D$ over which a test point $\mathbf{x}_{i}$ is believed to be an outlier.

Following the definition of a discordancy measure, all that remains to
construct the outlier analysis is the selection of a feature sensitive to
damage. As mentioned previously, natural frequencies are a common damage
sensitive feature to use given that they will shift inline with local
stiffness and mass changes that will be caused as a consequence of structural
damage. To demonstrate, Figure 15 plots the FRF extracted from the undamaged
and damaged (898g added mass at the tip) condition from sensor ULC-03 (sensor
positions are in reference to Figure 3) subject to burst random excitation at
2.0V. It can be seen that there is a shift in both the peak frequencies and
amplitudes in the damage case relative to the undamaged, confirming a suitable
damage sensitive feature.

Figure 15: Comparison between undamaged and damaged (898.0g added mass at the
tip) FRFs for sensor ULC-03, subject to burst-random excitation at 2.0V.

To mitigate the curse of dimensionality and avoiding taking the entire FRF as
the feature vector, it is sensible to focus on just one peak of the FRF (e.g.
around a natural frequency). Zooming in to the FRF around the natural
frequency at 156Hz, shown in Figure 16, it can be seen that there is still a
clear shift between the two damage cases. The feature vector is then taken to
be 8 spectral lines between 155 and 157Hz.

Figure 16: Zoom into damage sensitive region of the FRFs. Vertical lines
indicate the frequency range used as features for the novelty detection
algorithm.

Having identified a suitable input feature for the outlier analysis, it is
possible to use data that is known to have been collected from the structure
in of an undamaged case to form a baseline for the novelty index, which sets
the threshold for whether a datum point is believed to be an outlier, and thus
indicative of damage, or not. Beginning with just a comparison between data
from an undamaged and damaged condition, Figure 17(a) shows that the novelty
detection is easily able to distinguish between the two damage classes, which
includes varying levels of mass (pseudo-damage level), mass location and
excitation, as shown in Figures 17(b) \- 17(d). With outlier analysis,
however, one is limited to a binary state classification, which in this case
is whether the data is representative of a healthy or damaged(anomalous)
structural condition. Attempts to attempt to classify on the basis of the
novelty measure online will generally be futile, as shown in Figures 17(b) and
17(c), where there exists no discernible separation in the novelty measure
space for data across differences in mass amount and location.

(a) Presence of damage.

(b) Added mass amount.

(c) Damage location.

(d) Excitation level.

Figure 17: Results from the novelty detection algorithm, trained on half of
the normal condition data. Features are spectral lines between 155 and 157Hz.

##  6 Dataset challenges

The aim of this publication has been to layout a benchmark dataset for
practitioners interested in structural dynamics, SHM and digital twins. To
this end, a set of baseline analysis has been conducted to establish what is
possible with well known techniques within the field. As such, it seems
fitting to describe here what the remaining open challenges are when
considering this dataset. It is the hope of the authors that this (admittedly
non-exhaustive) list may prompt other members of the research community to
apply their state-of-the-art methods.

###  6.1 Identification

The first category of challenges in this dataset is related to the system
identification of a large and complex structure. From the perspective of the
authors, the purpose of this system identification task is to provide a model
which can:

  1. 1.

Replicate the dynamic behaviour of the physical system through a mathematical
model.

  2. 2.

Provide some level of insight into the physical mechanisms which have led to
the observed response of the system.

  3. 3.

Quantify the uncertainty/confidence of the modeller, perhaps considering if
the observed uncertainty is both epistemic and aleatoric.

Any given approach may target one or more of these objectives and, ideally,
all three. It is also worth noting that the analysis carried out should be
usable in such a way that it informs further decisions, e.g. the results of a
linear modal system identification may be used as features within an SHM
scheme.

Considering the methodologies shown so far in this publication which are
intended as a baseline, it is worth noting some of the major assumptions made
thus far which (the authors would argue) should and can be challenged in
future work. The analysis thus far has:

  * 1.

Considered only the _low_ frequency regime, this has neglected _local modes_ ,
e.g. the response of individual panels on the wing which may behave as thin
plates.

  * 2.

Relied on there being a _modal_ decomposition, such that the observed transfer
function of the system is a superposition of linear modes.

  * 3.

Extending this modal assumption, all analysis has considered the system
linear, despite some observed sensitivity in the FRF to excitation amplitude.

  * 4.

Considered the system to be completely time-invariant across a testing
campaign of several weeks.

  * 5.

The aircraft is treated directly in a reduced order form as a relatively small
number of ordinary differential equations (in the decoupled modal space),
compared to the overall complexity of the structure.

  * 6.

All analysis shown in this paper has avoided a formal treatment or
quantification of uncertainty, instead only point estimates are given.

Removing or replacing these assumptions and simplifications would be valuable
topics for further research. Some potentially interesting avenues of
investigation include: consideration of the best manner in which to build the
low order model of the system; verification of the identification results by
other algorithms; deeper investigation into the potential nonlinearity in the
system and its impact.

###  6.2 Prediction

In addition to the tasks listed above which consider the identification of the
system, it is also useful to imagine scenarios in which the practitioner may
wish to predict certain aspects of the system behaviour. The extension of the
identification procedures to prediction has not been shown in this work but
should be an area of investigation in the future. The value of such
predictions to the engineer would be to allow replication of (as yet) unseen
behaviours; to explore the response of the system in extreme events or
critical scenarios; the ability to infer the response of the system at
locations which are not measurable in operation (i.e. virtual sensing); to
allow validation of the identified models against other collected dataset; to
provide an objective means for comparison between models (especially if
looking at the relative performance of physics-based and data-driven models).

It is instructive to further consider two of these tasks in greater detail.
For the problem of virtual sensing which can be considered to reproducing data
which would have been measured at some dense (by some measure) set of sensors
on the basis of observations at a sparse set of sensors. There exist several
methodologies in the literature which provide promising approaches for the
virtual sensing problem. Applying these to the Hawk data presented here in
future work would form a further valuable benchmark on this system. In
addition to this there are several outstanding questions which are relevant
when considering the problem of virtual sensing. Firstly, one could explore
what would be the minimal set of sensors required to provide some bounded
uncertainty/error in the locations which are being reconstructed. Secondly, it
is important to explore which sensors are (in some sense) optimal to perform
the reconstruction and how the user may choose this subset. Thirdly, the
verification and validation of virtual sensing systems especially in terms of
their reconstruction under extrapolation, e.g. when operating in a larger
amplitude excitation regime.

Another prediction task of interest, for which this data could serve as a
benchmark, is that of load estimation. In many settings of interest, unlike
the scenario encountered in this laboratory test, it may be impractical to
acquire measurements of the loads applied to the structure. In aerospace
applications, it is very challenging to acquire the aerodynamic loads on the
wings during flight. Similar scenarios are encountered in wind turbine or
bridge applications. In all such cases, the fluid structure interaction forces
are typically unavailable. Owing to this difficulty, in recent years, there
has been some attention given to estimation of the loads alongside other
properties of interest in the system. This is referred to as either load
estimation/reconstruction or as joint input-state(-parameter) estimation [20].
In the context of using in service measurements for asset management,
knowledge regarding the loads is highly valuable as the loads are
intrinsically linked to the fatigue life of the structure. In practice the
estimation of these loads is complicated by noise on the measurements,
uncertainty regarding the dynamic system itself and a lack of full state
information/observation, i.e. only sparse measurements are available.

###  6.3 Evaluation

In addition to the understanding and reconstruction of the dynamic system, it
is of interest to consider the follow-on tasks which this form of
testing/monitoring will enable. One example of a task for which dynamics might
be studied (and demonstrated in this paper) is that of SHM which enables
automated asset management through monitoring. This desire for more tightly
integrated modelling, incorporation of data and asset management may also be
characterised, more recently, using the language of digital twins [21, 22]. In
this paper tasks related to the lower levels of Rytter’s hierarchy have been
shown. For problems of damage detection, it has been seen in this work how
(given insight about the form of the damage or some prior knowledge of
sensitive features) it is possible to build a good quality model. However,
doing this in the more representative context of only having access to the
outputs or without good prior knowledge of the features may yet pose a
challenge. What has been clear is that the more challenging tasks in the
hierarchy of localisation, classification and extent remain very challenging
based on the features extracted. This difficulty is present even when
considering the high level of prior knowledge used in this study.
Additionally, the interaction of the SHM procedure with, for example, the
nonlinearity or other confounding influences remains to be explored. The topic
prognosis has not been approached in this work and, although different extents
of (proxy) damage are used, this would not easily be accomplished with this
benchmark dataset. Some additional tasks which may add high value to the SHM
community would be to attempt to develop a generative model of damage which
would allow classifiers to be trained without examples of the previous damage
cases. For example, this could be to reproduce the effect of damage at the
trailing edge of the wing from the example data collected at the tip.

This section has attempted to provide a non-exhaustive list of future avenues
of investigation for which this dataset could prove a useful benchmark. It is
the belief of the authors that the presented data represents a challenging
full-scale example of some of the open questions within the field of
structural dynamics and SHM. Hopefully, this discussion will stimulate further
investigation of this dataset and subsequent research.

##  7 Conclusion

In this work, a new, large-scale structural dynamics benchmark dataset is
introduced. Details of the structure, experimental campaign, acquisition and
control systems are presented and the data is made freely available. As well
as the supplied data, this work also presents a first analysis of several
dataset challenges using standard approaches in structural dynamics. For the
task of system identification, a _rational fraction decomposition_ is applied
in the frequency domain, and this is compared to an output only time domain
identification in _stochastic subspace identification_.

For the task of damage identification, a _supervised_ machine learning
approach based on _Gaussian mixture models_ on the principal coordinates of
the natural frequencies is considered, with an outlier analysis presented as a
first benchmark as an _unsupervised_ method.

The analyses and results presented here are not groundbreaking and are
intended to serve as a baseline, against which, more advanced methods might be
compared. Some challenges which the dataset presents have been discussed, with
the view to inspire readers with motivation for new approaches to structural
dynamics analysis.

## Data availability

All data and documentation as well as a python API for accessing and
interacting with the data are openly available at
https://doi.org/10.15131/shef.data.22710040.v1.

## Author contributions

MHA: Conceptualisation (Supporting), Methodology (Lead (Experimental)),
Software (Lead (Experimental)), Investigation (Lead), Data Curation
(Supporting), Writing - Original Draft (Lead), Writing - Review & Editing
(Supporting)

RM: Conceptualisation (Supporting), Methodology (Supporting), Software
(Supporting), Investigation (Supporting), Resources (Lead), Data Curation
(Supporting)

MDC: Methodology (Lead (System ID)), Software (Lead (System ID)),
Investigation (Supporting), Data Curation (Lead), Writing - Original Draft
(Supporting), Writing - Review & Editing (Supporting)

MRJ: Methodology (Lead (SHM)), Software (Lead (SHM)), Investigation
(Supporting), Data Curation (Supporting), Writing - Original Draft
(Supporting), Writing - Review & Editing (Supporting)

MSB: Data Curation (Supporting)

DW: Resources (Supporting), Writing - Review & Editing (Lead), Supervision
(Supporting), Project administration (Lead), Funding acquisition (Lead)

TJR: Conceptualisation (Lead), Methodology (Supporting), Software
(Supporting), Resources (Supporting), Data Curation (Supporting), Writing -
Original Draft (Supporting), Writing - Review & Editing (Supporting),
Supervision (Lead), Project administration (Supporting), Funding acquisition
(Supporting)

## Acknowledgements

The authors of the paper are thankful for the support of The Alan Turing
Institute under the Data-Centric Engineering programme, specifically through
the project Digital twins for high-value engineering applications (DTHIVE), as
well as the EPSRC through grants EP/R006768/1, EP/W002140/1 and EP/S001565/1.
The experimental work made use of The Laboratory for Verification and
Validation (LVV), which was funded by the EPSRC (grant numbers EP/J013714/1
and EP/N010884/1),

The authors are also extremely grateful for the stimulating and insightful
discussions during an Isaac Newton Institute for Mathematical Sciences data-
driven deep-dive "From physics-based to data-driven assessment of structures"
in February 2023, which was supported through grant no EP/R014604/1.

Finally, the authors extend their thanks to Aidan Hughes, Brandon O’Connell,
Daniel Brennan, Daniel Pitchforth, George Tsialiamanis, Jack Poole, Jacques
Mclean, Joe Longbottom, Tristan Gowdridge for assistance in collecting the
data.

## References

  * [1] R. R. Craig Jr, A. J. Kurdila, Fundamentals of structural dynamics, John Wiley & Sons, 2006.
  * [2] D. H. Hodges, G. A. Pierce, Introduction to structural dynamics and aeroelasticity, Vol. 15, cambridge university press, 2011.
  * [3] F. M. Hemez, S. W. Doebling, Review and assessment of model updating for non-linear, transient dynamics, Mechanical Systems and Signal Processing 15 (1) (2001) 45–74.
  * [4] G. Kerschen, K. Worden, A. F. Vakakis, J.-C. Golinval, Past, present and future of nonlinear system identification in structural dynamics, Mechanical systems and signal processing 20 (3) (2006) 505–592.
  * [5] K. Worden, Nonlinearity in structural dynamics: detection, identification and modelling, CRC Press, 2019.
  * [6] T. J. Rogers, G. R. Holmes, E. Cross, K. Worden, On a grey box modelling framework for nonlinear system identification, in: Special Topics in Structural Dynamics, Volume 6: Proceedings of the 35th IMAC, A Conference and Exposition on Structural Dynamics 2017, Springer, 2017, pp. 167–178.
  * [7] K. Worden, R. Barthorpe, E. Cross, N. Dervilis, G. Holmes, G. Manson, T. Rogers, On evolutionary system identification with applications to nonlinear benchmarks, Mechanical Systems and Signal Processing 112 (2018) 194–232.
  * [8] R. Fuentes, R. Nayek, P. Gardner, N. Dervilis, T. Rogers, K. Worden, E. Cross, Equation discovery for nonlinear dynamical systems: A bayesian viewpoint, Mechanical Systems and Signal Processing 154 (2021) 107528.
  * [9] H. Fraser-Mitchell, The Hawk story, J. Aeronautical Hist (2011).
  * [10] M. Haywood-Alexander, R. Mills, M. Champneys, M. Bonney, D. Wagg, T. Rogers, BAE T1A Hawk Starboard Wing Modal Test (5 2023).  doi:10.15131/shef.data.22710040.v1.
URL
https://figshare.shef.ac.uk/articles/dataset/BAE_T1A_Hawk_Starboard_Wing_Modal_Test/22710040

  * [11] D. J. Ewins, Modal testing: theory, practice and application, John Wiley & Sons, 2009.
  * [12] M. Richardson, D. Formenti, Global curve fitting of frequency response measurements using the rational fraction polynomial method, in: Proceeding of 3rd IMAC, 1985, pp. 390–397.
  * [13] P. Van Overschee, B. De Moor, N4sid: Subspace algorithms for the identification of combined deterministic-stochastic systems, Automatica 30 (1) (1994) 75–93.  doi:10.1016/0005-1098(94)90230-5.
  * [14] T. Katayama, Subspace methods for system identification, Vol. 1, Springer, London, 2005.
  * [15] B. Peeters, G. De Roeck, Stochastic system identification for operational modal analysis: a review, J. Dyn. Sys., Meas., Control 123 (4) (2001) 659–667.
  * [16] E. Reynders, K. Maes, G. Lombaert, G. De Roeck, Uncertainty quantification in operational modal analysis with stochastic subspace identification: Validation and applications, Mechanical Systems and Signal Processing 66 (2016) 13–30.
  * [17] H. Hassani, Singular spectrum analysis: Methodology and comparison, Journal of Data Science 5 (2) (2007) 239–257.  doi:10.6339/jds.2007.05(2).396.
  * [18] C. R. Farrar, K. Worden, Structural health monitoring: a machine learning perspective, John Wiley & Sons, 2012.
  * [19] K. Worden, J. M. Dulieu-Barton, An overview of intelligent fault detection in systems and structures, Structural Health Monitoring 3 (1) (2004) 85–98.
  * [20] S. Vettori, Development of virtual sensing techniques for dynamic virtualization, Doctoral thesis, ETH Zurich, Zurich (2023).  doi:10.3929/ethz-b-000623754.
  * [21] D. J. Wagg, K. Worden, R. J. Barthorpe, P. Gardner, Digital Twins: State-of-the-Art and Future Directions for Modeling and Simulation in Engineering Dynamics Applications, ASCE-ASME J Risk and Uncert in Engrg Sys Part B Mech Engrg 6 (3), 030901 (05 2020).  arXiv:https://asmedigitalcollection.asme.org/risk/article-pdf/6/3/030901/6534130/risk\\_006\\_03\\_030901.pdf, doi:10.1115/1.4046739.
URL https://doi.org/10.1115/1.4046739

  * [22] P. Gardner, M. Dal Borgo, V. Ruffini, A. J. Hughes, Y. Zhu, D. J. Wagg, Towards the development of an operational digital twin, Vibration 3 (3) (2020) 235–265.


$\displaystyle\tilde{\Psi}_{{\rm
QM}}(\overline{p},\sigma)=\sum_{\sigma^{\prime}=1}^{4}\;\delta(p^{0}-E_{p}^{(\sigma^{\prime})})\alpha_{\sigma^{\prime}}(\vec{p})\;u_{\sigma,\sigma^{\prime}}(\vec{p})\;,$
(184)

with $\alpha_{\sigma^{\prime}}(\vec{p})$ arbitrary functions, i.e.

$\displaystyle\Psi_{{\rm QM}}(\vec{x},\sigma|t)$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}p}{(2\pi)^{\frac{3}{2}}}e^{i\vec{p}\cdot\vec{x}}\sum_{\sigma^{\prime}=1}^{4}\;\frac{u_{\sigma,\sigma^{\prime}}(\vec{p})}{\sqrt{2\pi}}$
(185) $\displaystyle\times
e^{-iE_{p}^{(\sigma^{\prime})}t}\;{\alpha_{\sigma^{\prime}}(\vec{p})}\;,$

at the level of the 3D+1 spinor wave-function.

Expressed as in Eq. (173) it is easy to verify that, at the level of the GEB
distribution $|\Psi_{{\rm QM}}\rangle=\sum_{\sigma=1}^{4}\int
d^{4}x\;\Psi_{{\rm QM}}(\vec{x},\sigma|t)|\overline{x},\sigma\rangle$, the
Dirac equation (47) corresponds to the identity $J_{\rm D}|\Psi_{{\rm
QM}}\rangle=0$ with $J_{\rm D}$ as in Eq. (III.1). Indeed, to show this, we
need the fact that thanks to (5) $\tilde{\Psi}_{{\rm
QM}}(\overline{p},\sigma)$ provides the 4D-momentum spinor wave-functions
expansion of $|\Psi_{{\rm QM}}\rangle$, i.e. $|\Psi_{{\rm
QM}}\rangle=\sum_{\sigma=1}^{4}\int d^{4}p\;\tilde{\Psi}_{{\rm
QM}}(\overline{p},\sigma)|\overline{p},\sigma\rangle$. As mentioned in the
main text the operator $J_{\rm D}$ is not a self-adjoint: this is a direct
consequence of the fact that for all $i=1,2,3$ the matrices $\gamma^{i}$ are
anti-Hermitian (indeed $(\gamma^{i})^{\dagger}=-\gamma^{i}=\gamma_{i}$), while
$\gamma^{0}$ is Hermitian, so that $J_{\rm
D}^{\dagger}=\sum_{\mu=1}^{4}(\gamma^{\mu})^{\dagger}P_{\mu}-m=\sum_{\mu=1}^{4}\gamma_{\mu}P_{\mu}-m\neq
J_{\rm D}$. Notice however that exploiting the fact that
$\gamma^{0}\gamma^{0}=\openone$, and
$\gamma^{0}\gamma^{i}=\left(\begin{matrix}0&\sigma_{i}\cr\sigma_{i}&0\end{matrix}\right)$,
we can write

$\displaystyle J_{\rm D}=\gamma^{0}J_{\rm D}^{(H)}\;,$ (186)

where given the matrix elements
$M_{\sigma^{\prime},\sigma^{\prime\prime}}(\overline{p})$ of Eq. (176) $J_{\rm
D}^{(H)}$ is the self-adjoint operator

$\displaystyle J_{\rm
D}^{(H)}:=\sum_{\sigma^{\prime},\sigma^{\prime\prime}}\int
d^{4}p\;M_{\sigma^{\prime},\sigma^{\prime\prime}}(\overline{p})|\overline{p},\sigma^{\prime}\rangle\langle\overline{p},\sigma^{\prime\prime}|\;.$
(187)

Equation (51) finally follows by using Eq. (183) observing that the vectors

$\displaystyle|\phi_{\sigma}(\overline{p})\rangle:=\sum_{\sigma^{\prime}=1}^{4}u_{\sigma^{\prime},\sigma}(\vec{p})|\overline{p},\sigma^{\prime}\rangle\;,$
(188)

obey generalized orthonormal conditions (52) thanks to the unitary properties
of the matrix elements $u_{\sigma^{\prime},\sigma}(\vec{p})$: indeed with this
choice Eq. (187) becomes

$\displaystyle J_{\rm D}^{(H)}=\sum_{\sigma=1}^{4}\int
d^{4}p\;\lambda_{\sigma}(\overline{p})\;|\phi_{\sigma}(\overline{p})\rangle\langle\phi_{\sigma}(\overline{p})|\;,$
(189)

and hence

$\displaystyle K_{\rm D}$ $\displaystyle=$ $\displaystyle J_{\rm
D}^{\dagger}J_{\rm D}=\left(J^{(H)}_{\rm D}\right)^{2}$ (190) $\displaystyle=$
$\displaystyle\sum_{\sigma=1}^{4}\int
d^{4}p\;\lambda^{2}_{\sigma}(\overline{p})\;|\phi_{\sigma}(\overline{p})\rangle\langle\phi_{\sigma}(\overline{p})|\;.$

We conclude reporting explicit expressions for the
$|\phi_{\sigma}(\overline{p})\rangle$:

$\displaystyle|\phi_{1}(\overline{p})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}\sqrt{1-\frac{m}{E_{p}}}|\overline{p},s_{1}(\hat{n})\rangle-\frac{|\vec{p}|}{\sqrt{1-\frac{m}{E_{p}}}}|\overline{p},s_{2}(\hat{n})\rangle\Big{)}\;,$
(191) $\displaystyle|\phi_{2}(\overline{p})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}\sqrt{1+\frac{m}{E_{p}}}|\overline{p},s_{1}(\hat{n})\rangle+\frac{|\vec{p}|}{\sqrt{1+\frac{m}{E_{p}}}}|\overline{p},s_{2}(\hat{n})\rangle\Big{)}\;,$
(192) $\displaystyle|\phi_{3}(\overline{p})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}\sqrt{1-\frac{m}{E_{p}}}|\overline{p},s_{3}(\hat{n})\rangle+\frac{|\vec{p}|}{\sqrt{1-\frac{m}{E_{p}}}}|\overline{p},s_{4}(\hat{n})\rangle\Big{)}\;,$
(193) $\displaystyle|\phi_{4}(\overline{p})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\Big{(}\sqrt{1+\frac{m}{E_{p}}}|\overline{p},s_{3}(\hat{n})\rangle-\frac{|\vec{p}|}{\sqrt{1+\frac{m}{E_{p}}}}|\overline{p},s_{4}(\hat{n})\rangle\Big{)}\;,$
(194)

where for $\hat{n}:=\vec{p}/|\vec{p}|$,
$|\overline{p},s_{\sigma}(\hat{n})\rangle$ are the orthonormal vectors

$\displaystyle|\overline{p},s_{1}(\hat{n})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\sqrt{1+n^{3}}\;|\overline{p},1\rangle+\frac{n^{1}+in^{2}}{\sqrt{1+n^{3}}}\;|\overline{p},2\rangle\right)\;,$
(195) $\displaystyle|\overline{p},s_{2}(\hat{n})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\sqrt{1+n^{3}}\;|\overline{p},3\rangle+\frac{n^{1}+in^{2}}{\sqrt{1+n^{3}}}\;|\overline{p},4\rangle\right)\;,$
(196) $\displaystyle|\overline{p},s_{3}(\hat{n})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\sqrt{1-n^{3}}\;|\overline{p},1\rangle-\frac{n^{1}+in^{2}}{\sqrt{1-n^{3}}}\;|\overline{p},2\rangle\right)\;,$
(197) $\displaystyle|\overline{p},s_{4}(\hat{n})\rangle$ $\displaystyle:=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\sqrt{1-n^{3}}\;|\overline{p},3\rangle-\frac{n^{1}+in^{2}}{\sqrt{1-n^{3}}}\;|\overline{p},4\rangle\right)\;.$
(198)

Observe that via Eq. (188) these identities implicitly define the matrix
elements $u_{\sigma,\sigma^{\prime}}(\vec{p})$: for instance we get

$\displaystyle u_{1,1}(\vec{p})$ $\displaystyle=$
$\displaystyle\frac{1}{{2}}\sqrt{1-\frac{m}{E_{p}}}\sqrt{1+n^{3}}\;,\qquad
u_{2,1}(\vec{p})=\frac{1}{{2}}\sqrt{1-\frac{m}{E_{p}}}\frac{n^{1}+in^{2}}{\sqrt{1+n^{3}}}\;,$
(199) $\displaystyle u_{3,1}(\vec{p})$ $\displaystyle=$
$\displaystyle-\frac{1}{{2}}\frac{|\vec{p}|}{\sqrt{1-\frac{m}{E_{p}}}}\sqrt{1+n^{3}}\;,\qquad\cdots$
(200)

If one appropriately normalizes the state, one can recover, just as for
Bosons, a scalar probability density for each spinor component, since the
Dirac equation implies the Klein-Gordon one:
$(\gamma^{\mu}p_{\mu}+m)(\gamma^{\nu}p_{\nu}-m)=(p^{\mu}p_{\mu}-m^{2})\openone_{4}$.
In addition, one can, as usual, also introduce a conserved 4-current
$j^{\mu}\equiv\Psi^{\dagger}(\overline{x})\gamma^{0}\gamma^{\mu}\Psi(\overline{x})=\bar{\Psi}(\overline{x})\gamma^{\mu}\Psi(\overline{x})$,
where $\Psi$ is the column vector of the conditioned spinors in the position
representation, namely the column of position-representation amplitudes. Since
the zeroth component $j^{0}=\Psi^{\dagger}(\overline{x})\Psi(\overline{x})$ is
positive definite, it can be given a probabilistic interpretation as the
conditional probability density of finding a particle at position $\vec{x}$,
given that time is $t$, where $\overline{x}=({t},\vec{x})$.

##  Appendix H Extra observations on the multi-event QM/GEB correspondence

This section is dedicated to making explicit some technical aspects of the
QM/GEB correspondence in multi-event scenario discussed in Sec. III.2. We
start by showing that the vector (57) is uniquely defined; then we verify that
Eq. (64) gives the right prescription to compute the evolution of a QM 3D+1
spinor wave-function under Lorentz transformations.

###  H.1 Uniqueness of the the multi-event QM/GEB correspondence

Here we prove that the vector (57) is uniquely defined.

To begin with recall that, for all $t$, the QM spinor 3D wave-function
$\Psi^{[n]}_{\rm QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)$
of $n$ particles can be expressed as

$\Psi^{[n]}_{\rm
QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)=\langle\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|\psi^{[n]}_{\rm
QM}(t)\rangle\;,$ (201)

where

$\displaystyle|\psi^{[n]}_{\rm
QM}(t)\rangle=\sum_{\sigma_{1},\cdots,\sigma_{n}}\int d^{3}{x}_{1}\cdots\int
d^{3}{x}_{n}$ (202) $\displaystyle\quad\quad\Psi^{[n]}_{\rm
QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)|\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}\rangle\;,$

is the associated wave-vector, and $|\vec{x}_{j},\sigma_{j}\rangle$ the
generalized 3D position eigenvectors of the $j$-th particle. Recalling then
that we are dealing with non-interacting systems, we can now write

$\displaystyle|\psi^{[n]}_{\rm QM}(t)\rangle=U_{1}(t)\otimes\cdots\otimes
U_{n}(t)|\psi^{[n]}_{\rm QM}(0)\rangle\;,$ (203)

where for $j=1,\cdots,n$, $U_{j}(t)$ stands for the QM unitary transformation
that rules the free evolution of the $j$-th particle. Equation (203) leads to
Eq. (55) by expressing $|\psi^{[n]}_{\rm QM}(0)\rangle$ in terms of an
arbitrary local basis for the $n$ particles, i.e.

$|\psi^{[n]}_{\rm
QM}(0)\rangle=\sum_{\vec{\ell}}\alpha_{\vec{\ell}}\;|\psi^{(\ell_{1})}_{\rm
QM}\rangle\otimes\cdots\otimes|\psi^{(\ell_{n})}_{\rm QM}\rangle\;,$ (204)

and using the identities

$\displaystyle\Psi^{(\ell_{j})}_{\rm QM}(\vec{x}_{j},\sigma_{j}|t)$
$\displaystyle:=$
$\displaystyle\langle\vec{x}_{j},\sigma_{j}|U_{j}(t)|\psi^{(\ell_{j})}_{\rm
QM}\rangle\;.$ (205)

Replacing this into (56) and (57) we finally arrive to

$\displaystyle|\Psi^{[n]}_{\rm
QM}\rangle=\sum_{\sigma_{1},\cdots,\sigma_{2}}\int d^{4}x_{1}\cdots\int
d^{4}x_{n}|\overline{x}_{1},\sigma_{1};\cdots;\overline{x}_{n},\sigma_{n}\rangle$
$\displaystyle\;\times\langle\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|U_{1}(t_{1})\otimes\cdots\otimes
U_{n}(t_{n})|\psi^{[n]}_{\rm QM}(0)\rangle,$

which explicitly shows that $|\Psi_{\rm QM}\rangle$ carries no functional
dependence upon the specific choice of the local basis $\\{|\psi^{(\ell)}_{\rm
QM}\rangle\\}_{\ell}$ used in (204).

As an application of the above identities we report here the special cases of
particles obeying to the positive energy KG equation and the Dirac equation.
For the KG equation, setting $\tilde{\psi}_{\rm
QM}^{[n]}(\vec{p}_{1},\cdots,\vec{p}_{n}):=\sum_{\vec{\ell}}\alpha_{\vec{\ell}}{\psi}^{(\ell_{1})}(\vec{p}_{1})\cdots{\psi}^{(\ell_{2})}(\vec{p}_{n})$
and replacing (153) into Eq. (55) we get

$\displaystyle{\Psi}_{\rm QM}^{[n]}(\vec{x}_{1},\cdots,\vec{x}_{n}|t)$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{(2\pi)^{3/2}}\cdots\frac{d^{3}p_{n}}{(2\pi)^{3/2}}e^{i(\vec{p}_{1}\cdot\vec{x}_{1}+\cdots+\vec{p}_{n}\cdot\vec{x}_{n})}e^{-i(E_{p_{1}}+\cdots+E_{p_{n}})t}\;\tilde{\psi}_{\rm
QM}^{[n]}(\vec{p}_{1},\cdots,\vec{p}_{n})\;,$ (206)

with associated 4D GEB spinor wave-function

$\displaystyle{\Psi}_{\rm QM}^{[n]}(\overline{x}_{1},\cdots,\overline{x}_{n})$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{(2\pi)^{3/2}}\cdots\frac{d^{3}p_{n}}{(2\pi)^{3/2}}e^{i(\vec{p}_{1}\cdot\vec{x}_{1}+\cdots+\vec{p}_{n}\cdot\vec{x}_{n})}e^{-i(E_{p_{1}}t_{1}+\cdots+E_{p_{n}}t_{n})}\;\tilde{\psi}_{\rm
QM}^{[n]}(\vec{p}_{1},\cdots,\vec{p}_{n})\;.$ (207)

Similarly for the Dirac equation from Eq. (185), setting

$\displaystyle\phi^{[n]}(\vec{p}_{1},\sigma^{\prime}_{1};\cdots;\vec{p}_{n},\sigma^{\prime}_{n}):=\sum_{\vec{\ell}}\alpha_{\vec{\ell}}\;{\alpha_{\sigma^{\prime}_{1}}(\vec{p}_{1})}\cdots{\alpha_{\sigma^{\prime}_{n}}(\vec{p}_{n})}\;,$
(208)

we get

$\displaystyle\Psi^{[n]}_{{\rm
QM}}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)$ $\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{(2\pi)^{3/2}}\cdots\frac{d^{3}p_{n}}{(2\pi)^{3/2}}\sum_{\sigma^{\prime}_{1}=1}^{4}\frac{u_{\sigma_{1}^{\prime},\sigma_{1}}(\vec{p}_{1})}{\sqrt{2\pi}}\cdots\sum_{\sigma^{\prime}_{n}=1}^{4}\frac{u_{\sigma_{n}^{\prime},\sigma_{n}}(\vec{p}_{n})}{\sqrt{2\pi}}\;e^{i(\vec{p}_{1}\cdot\vec{x}_{1}+\cdots+\vec{p}_{n}\cdot\vec{x}_{n})}$
(209) $\displaystyle\qquad\times
e^{-i(E^{(\sigma^{\prime}_{1})}_{p_{1}}+\cdots+E^{(\sigma^{\prime}_{n})}_{p_{n}})t}\;\phi^{[n]}(\vec{p}_{1},\sigma^{\prime}_{1};\cdots;\vec{p}_{n},\sigma^{\prime}_{n})\;.$

which at the level of GEB corresponds to

$\displaystyle\Psi^{[n]}_{{\rm
QM}}(\overline{x}_{1},\sigma_{1};\cdots;\overline{x}_{n},\sigma_{n})$
$\displaystyle=$
$\displaystyle\int\frac{d^{3}p_{1}}{(2\pi)^{3/2}}\cdots\frac{d^{3}p_{n}}{(2\pi)^{3/2}}\sum_{\sigma^{\prime}_{1}=1}^{4}\frac{u_{\sigma_{1}^{\prime},\sigma_{1}}(\vec{p}_{1})}{\sqrt{2\pi}}\cdots\sum_{\sigma^{\prime}_{n}=1}^{4}\frac{u_{\sigma_{n}^{\prime},\sigma_{n}}(\vec{p}_{n})}{\sqrt{2\pi}}\;e^{i(\vec{p}_{1}\cdot\vec{x}_{1}+\cdots+\vec{p}_{n}\cdot\vec{x}_{n})}$
(210) $\displaystyle\qquad\times
e^{-i(E^{(\sigma^{\prime}_{1})}_{p_{1}}t_{1}+\cdots+E^{(\sigma^{\prime}_{n})}_{p_{n}}t_{n})}\;\phi^{[n]}(\vec{p}_{1},\sigma^{\prime}_{1};\cdots;\vec{p}_{n},\sigma^{\prime}_{n})\;.$

###  H.2 Lorentz transformations

Equation (201) represents the (time-dependent) 3D+1 spinor wave-function that
an observer $O$ would assign to describe the state of the $n$ particles on his
reference frame $R$. Assuming the dynamical evolution is relativistic
consistent (e.g. the particles obey KG or Dirac dynamical equations), we are
now interested in determining the spinor 3D wave-function a second observer
$O^{\prime}$ sitting in the reference frame $R^{\prime}$ with 4D coordinates
$\overline{x}^{\prime}$ that are linked with those of $R$ as in Eq. (11) will
assign to such a state. Since particles are independent (i.e. no interactions
are present in the model), this can be done using the decomposition (55) for
$\Psi^{[n]}_{\rm QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t$
and applying the single-particle transformation (54) to each individual term
$\Psi^{(\ell)}_{\rm QM}(\vec{x},\sigma|t)$. Accordingly we can write

$\displaystyle\\!\\!\\!\\!\\!\\!\Psi^{[n]\prime}_{\rm
QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)=\\!\\!\\!\\!\\!\\!\sum_{\sigma_{1},\cdots,\sigma_{n}}S_{\sigma^{\prime}_{1},\sigma_{1}}^{-1}(\Lambda)\cdots
S_{\sigma^{\prime}_{n},\sigma_{n}}^{-1}(\Lambda)$
$\displaystyle\\!\\!\\!\\!\\!\\!\times\sum_{\vec{\ell}}\alpha_{\vec{\ell}}\;\Psi^{(\ell_{1})}_{\rm
QM}(\Lambda^{-1}\overline{x}_{1},\sigma^{\prime}_{1})|_{t_{1}=t}\cdots\Psi^{(\ell_{n})}_{\rm
QM}(\Lambda^{-1}\overline{x}_{n},\sigma^{\prime}_{n})|_{t_{n}=t}\;,$ (211)

whose r.h.s. exactly matches with the one of Eq. (64) of the main text. It is
important to stress that while not immediately evident from the resulting
expression $\Psi^{[n]\prime}_{\rm
QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)$ does not depends
upon the specific choice of the local decomposition used in (204). One easy
way to verify this is e.g. to use the fact that (57) does not depends on such
a choice (see previous section) and the fact that thanks to Eq. (64) we can
write

$\displaystyle\Psi^{[n]\prime}_{\rm
QM}(\vec{x}_{1},\sigma_{1};\cdots;\vec{x}_{n},\sigma_{n}|t)$ (212)
$\displaystyle\qquad=\langle\overline{x}_{1},\sigma_{1};\cdots;\overline{x}_{n},\sigma_{n}|\Psi^{[n]\prime}_{\rm
QM}\rangle\Big{|}_{t_{1}=\cdots=t_{n}=t}$
$\displaystyle\qquad=\langle\overline{x}_{1},\sigma_{1};\cdots;\overline{x}_{n},\sigma_{n}|U^{\otimes
n}_{\Lambda}|\Psi^{[n]}_{\rm QM}\rangle\Big{|}_{t_{1}=\cdots=t_{n}=t}\;,$

where we made use of the invertion formula (58) and of Eq. (63).