∎

<EMAIL_ADDRESS><EMAIL_ADDRESS>
11institutetext: Department of Physics and Astronomy, Stony Brook University,
Stony Brook, New York 11794-3800, USA 22institutetext: Institute of
Theoretical Physics, Jagiellonian University, S. Łojasiewicza 11, 30-348
Kraków, Poland

# Scaling properties of direct photon yields
in heavy ion collisions

Vladimir Khachatryane1,addr1
Michał Praszałowicz e2,addr2

(Received: date / Accepted: date)

###### Abstract

A recent analysis from the PHENIX collaboration of available direct photon
measurement results in collisions of various systems such as Au+Au, Cu+Cu, and
Pb+Pb, at different beam energies ranging from 39 to 2760 GeV, has shown a
universal, within experimental uncertainties, multiplicity scaling, in which
direct photon $p_{T}$-spectra for transverse momenta up to 2 GeV/$c$ are
scaled with charged hadron pseudorapidity density at midrapidity raised to
power $\alpha=1.25$. On the other hand, those direct photon $p_{T}$-spectra
also exhibit geometrical scaling in the similar $p_{T}$ range. Assuming power-
law dependence of the scaled photon spectra for both scaling laws, we
formulate two independent conditions for the power $\alpha$, which overshoot
experimental data by $\sim 10\%$ on average. We discuss possible sources that
might improve this estimate.

###### Keywords:

Direct photons, heavy ion collisions, saturation, geometrical scaling.

††journal: Eur. Phys. J. C

##  1 Introduction

Measurements of direct photons provide unique opportunities in probing and
studying the properties and evolution of the matter produced in heavy ion
collisions (HIC). These photons are defined to be produced from all the
sources except for hadronic decays. Since they hardly interact with the
“fireball” of quarks and gluons due to a small interaction cross section with
the medium, the information they carry from the time of their production is
not washed out by final state interactions. Experimentally measured low
momentum direct photon $p_{T}$-spectra by PHENIX (in Au+Au at
$\sqrt{s_{{}_{NN}}}$ = 200 GeV and 62.4 GeV, in Cu+Cu at $\sqrt{s_{{}_{NN}}}$
= 200 GeV) Adare:2008ab ; Afanasiev:2012dg ; Adare:2014fwh ; Adare:2018jsz ;
Adare:2018wgc ; Khachatryan:2018evz ; Khachatryan:2018uzc  and ALICE (in Pb+Pb
at $\sqrt{s_{{}_{NN}}}$ = 2760 GeV) Adam:2015lda  collaborations in HIC are
enhanced with respect to $N_{\rm coll}$ (number of binary nucleon collisions)
scaled reference yield (measured or calculated) in p+p collisions. Low
momentum direct photon measurements by STAR collaboration show less
enhancement STAR:2016use . Earlier low energy WA98 data Aggarwal:2000th  have
mostly upper bounds in the relevant $p_{T}\leq 2$ GeV$/c$ region. Direct
photons in HIC also show large anisotropy (elliptic flow) Khachatryan:2018uzc
; Adare:2015lcd ; Acharya:2018bdy .

Figure 1: Direct photon invariant yield as a function of $p_{T}$ for various
colliding nuclei and collision centrality selections at three center-of-mass
energies. The plot labels give the references of the shown data. The vertical
lines of the data points describe statistical uncertainties, the boxes
describe systematic uncertainties.

Fig. 1 shows several data sets of direct photon $p_{T}$-spectra at low- ($<$ 1
GeV$/c$) and intermediate- (from 1 GeV$/c$ up to $\sim$ 5 GeV$/c$) $p_{T}$
regions. There have been many theoretical attempts to reproduce the photon
yields shown in Fig. 1 and flow coefficients (that are not discussed here)
with, however, mixed success. Hydrodynamical simulations of the fireball
evolution DPSYJG:2011 ; Shen:2013cca ; SHPG:2014 ; Paquet:2015lta ,
calculations in the framework of the elliptic-fireball expansion scenario
vanHees:2011 ; Rapp:2014 ; vanHees:2015 , Parton-Hadron-String Dynamics
transport approach Bratkovskaya:2008 ; Bratkovskaya:2014 ; Linnyk:2014 ;
Linnyk:2015 , as well as the spectral function approach Dusling:2010 ;
Duslingthesis ; Lee:2014 ; Kim:2017 , have encountered difficulties in the
simultaneous description of the observed large yields and large anisotropies,
which have given rise to a challenging problem, commonly referred to as direct
(thermal) photon flow puzzle Shen:2016odt  (for recent review see
David:2019wpt ).

The early attempts to describe the photon yields by initial state models can
be found in Refs. Chiu:2012ij ; Klein-Bosing:2014uaa ; McLerran:2014oea ;
McLerran:2014hza ; McLerran:2015mda . More recently initial state models have
been used both as a part of hydrodynamic evolution Paquet:2015lta  and in a
bottom-up thermalization scenario Berges:2017eom ; Khachatryan:2018ori . In
the latter case Khachatryan:2018ori  good fits for the photon yields have been
obtained for PHENIX and ALICE data (the anisotropies are not discussed there).

Much work still needs to be done to include the anisotropies in the initial
state model calculations, in order to directly address the aforementioned
puzzle. Ref. McLerran:2014hza  discusses the possibility of late stage
elliptic flow effects to be explained in the combined Glasma and thermal
photon production scenarios, by which the photon mean emission time is shifted
toward a later time scale as compared to the mean emission time obtained only
from the thermal ansatz. Nonetheless, as it has been pointed out in
Khachatryan:2018ori , a more realistic description of the photon production
might come out within the combined framework of the initial and late stage
models. In this regard a promising step is already undertaken in
Paquet:2015lta , where the event-by-event hydrodynamical model uses IP-Glasma
initial state conditions.

However, the most important argument from Chiu:2012ij ; Klein-Bosing:2014uaa ;
Khachatryan:2018ori , directly applicable to the studies of this paper,
indicates that the initial state models potentially allow the photon
$p_{T}$-spectra to exhibit geometrical scaling (GS)111We will perform detailed
studies of GS of photon $p_{T}$-spectra in HIC elsewhere KhaPrasz .. In fact,
two types of scalings are observed in direct photon $p_{T}$ distributions in
HIC and hadron collisions:

### Direct photon multiplicity scaling (MS).

It is an experimental observation in HIC, which shows that direct photon
invariant yields follow one universal curve within experimental statistical
and systematic uncertainties for various colliding species and system sizes at
different center-of-mass energies, when scaled by charged hadron
multiplicities at midrapidity raised to the power of $\alpha=1.25$
Adare:2018wgc :

$\frac{1}{\left(dN_{\mathrm{ch}}/d\eta|_{\eta\approx
0}\right)^{\alpha}}\frac{dN_{\gamma}}{d^{2}p_{T}dy}=\frac{1}{Q_{0}^{2}}G(p_{T}),$
(1)

where $G$ is a universal energy- and multiplicity-independent function of
$p_{T}$ (see Fig. 2)222We will henceforth simply use $dN_{\mathrm{ch}}/d\eta$,
skipping the midrapidity notation $\eta\approx 0$ in the text, and keeping it
only in the figures., and $Q_{0}\sim 1$ GeV/$c$. MS holds for small and
intermediate $p_{T}$ up to $\sim 2\div 2.5$ GeV$/c$, precisely in the
enhancement region discussed above.

Figure 2: Direct photon $p_{T}$-spectra from Fig. 1 scaled by ${(dN_{\rm
ch}/d\eta)}^{1.25}$. The black curve corresponds to a power-law fit, the
details of which are shown in the plot. See the next sections for more
details.

### Geometrical scaling (GS).

Here the invariant yields of direct photons (and also of charged hadrons) can
be related to a universal energy-independent function of scaling variable
$\tau$ Klein-Bosing:2014uaa ; KhaPrasz ; McLerran:2010ex ; McLerran:2011 ;
Praszalowicz:2011rm ; Praszalowicz:2018vfw :

$\frac{1}{S_{T}}\frac{dN_{\gamma,\text{ch}}}{d^{2}p_{T}d\eta}=\,F_{\gamma,\text{ch}}(\tau),$
(2)

where $S_{T}$ is a parameter characterizing geometrical overlap area of
colliding nuclei. The variable $\tau$ is given by

$\tau=p_{T}/Q_{\mathrm{s}}(x).$ (3)

GS has been for the first time observed in deep inelastic scattering (DIS) in
e+p collisions Stasto:2000er . The quantity $Q_{\mathrm{s}}(x)$ in (3) is the
saturation scale, i.e., the transverse momentum at which the number of gluons
with longitudinal momentum fraction $x$, called Bjorken variable, cannot grow
any more due to the nonlinearity of QCD evolution. This happens at gluon
occupancy numbers being of the order of $1/\alpha_{s}$ ($\alpha_{s}$ – strong
coupling constant). The phenomenon of gluon saturation is generated from such
dynamics Gribov:1984tu ; Mueller:1985wy ; BK ; Kovchegov1 ; Kovchegov2 ;
Mueller:2002zm ; Munier:2003vc ; MPesch  and can be described by the Color
Glass Condensate (CGC) effective theory of high energy scattering
McLerran:1993ni ; McLRV2 ; McLRV2 ; elis:2010nm; McLerran:2010ub .

One still should remember that GS is more general. Indeed, it has been shown
in Ref. Caola:2008xr  that the approximate GS is in fact a general property of
solutions of the QCD collinear evolution with general boundary conditions.
Furthermore, the results of Ref. Iancu:2002tr  indicate that the rapidity
evolution preserves GS beyond saturation when starting from the boundary
conditions in the saturation region. It means that the scaling properties of
an initial state are preserved and even built up by the QCD evolution. This
may explain the fact that the geometrical scaling is observed in the hadronic
$p_{T}$-spectra and even more so in the case of photons.

Saturation models predict specific dependence of $Q_{\mathrm{s}}$ on the
collision energy, transverse momentum and collision centrality class,
quantified by the number of participating nucleons

$Q_{\rm{s}}=N_{\rm{part}}^{\delta/2}\,Q_{0}\left(\frac{p_{T}}{\sqrt{s_{{}_{\\!NN}}}\times
10^{-3}}\right)^{-\lambda/2},$ (4)

with $\delta\sim 2/3$ and $\lambda\sim 0.2\div 0.35$ Praszalowicz:2018vfw .
Here we have chosen $x_{0}=10^{-3}$ for a typical value of Bjorken $x$ where
the saturation effects become important. The direct photon GS is illustrated
in Fig. 3.

Figure 3: Direct photon $p_{T}$-spectra from Fig. 1 exhibit geometrical
scaling when plotted in terms of the scaling variable $\tau$ (at
$\lambda\approx 0.2$) and divided by $S_{T}\sim N_{\rm part}^{2/3}$. The black
curve corresponds to another power-law fit, the details of which are shown in
the plot. See the next sections for more details.

In the present paper we formulate two conditions that allow one to calculate
the power $\alpha$ of MS in terms of the parameters that enter the theoretical
parametrization of GS, assuming power-law dependence of the scaled photon
spectra. These two conditions result in two predictions for $\alpha$: 1.34 and
1.40, i.e. $\sim$ 10% above the experimental value. Given the simplicity of
the present analysis and approximate nature of both scaling laws, this result
may be considered as a satisfactory one. Nevertheless, we also discuss
possible corrections that might improve this theoretical prediction of
$\alpha$.

##  2 Basics of Geometrical Scaling

Geometrical Scaling is based upon the assumption that bulk properties of
charged particles and photons share scaling properties encoded in the initial
state wave function of colliding hadrons/nuclei. While detailed models of
direct photon production in the early and intermediate stages of HIC can be
found in Refs. Chiu:2012ij ; McLerran:2014oea ; McLerran:2014hza ;
Berges:2017eom ; Khachatryan:2018ori  as it was noted in the introduction,
here we discuss general assumptions and necessary conditions leading to GS.

Differential gluon production cross-section can be written in terms of the
unintegrated gluon distributions Gribov:1981kg :

$\frac{d\sigma}{dyd^{2}p_{\mathrm{T}}}=\frac{C}{p_{{T}}^{2}}{\displaystyle\int}d^{2}\vec{k}_{T}\,\alpha_{\mathrm{{s}}}\varphi^{(1)}(x_{1},\vec{k}_{T}^{2})\,\varphi^{(2)}(x_{2},(\vec{k}-\vec{p})_{T}^{2}),$
(5)

where $C$ contains color factors and numerical constants. The Bjorken $x$’s of
colliding partons read

$x_{1,2}=e^{\pm y}\,{p_{T}}/{\sqrt{s}}\,.$ (6)

In the midrapidity region $y\approx 0$, hence $x_{1}=x_{2}=x$.

There exist many models of unintegrated gluon distributions
$\varphi(x,\vec{k}_{T}^{2})$ (see _e.g._ GolecBiernat:1998js ;
GolecBiernat:1999 ; Kharzeev:2002ei ; Szczurek:2003vq ) that enter (5). Most
of them share two important features: geometrical scaling and dependence on
the transverse area parameter $A_{T}$:

$\varphi(x,\vec{k}_{\mathrm{{T}}}^{2})=A_{T}\,\phi(k_{T}^{2}/Q_{\mathrm{{s}}}^{2}(x)),$
(7)

where $\phi$ is a dimensionless function of the scaling variable
$\tau^{2}=k_{T}^{2}/Q_{\mathrm{{s}}}^{2}(x)$. The precise meaning of $A_{T}$
is best understood in a picture where the impact parameter is also taken into
account Kowalski:2003hm ; Tribedy:2010ab .

Ignoring momentum dependence of the strong coupling constant, on dimensional
analysis grounds, we arrive at

$\frac{d\sigma}{dyd^{2}p_{{T}}}=A^{(1)}_{T}A^{(2)}_{T}{F}(\tau),$ (8)

where ${F}$ is a universal, energy-independent function of the scaling
variable $\tau$ in (3).

If one assumes that $A_{T}$ is an energy-independent constant, which is true
in the case of the GBW model GolecBiernat:1998js ; GolecBiernat:1999  for DIS,
then it is the differential cross-section that should exhibit GS. Indeed, it
has been shown in Ref. Praszalowicz:2015dta  that for charged particles the
differential cross-section in p+p scattering scales better (over larger
$p_{T}$ interval) than multiplicity, and the exponent $\lambda\approx 0.3$ is
compatible with the DIS scaling Praszalowicz:2012zh . In contrast, for
multiplicity scaling $\lambda\approx 0.2$.

In order to obtain the multiplicity distribution, we should divide both sides
of (8) by an appropriate inelastic cross-section:

$\frac{dN}{dyd^{2}p_{{T}}}=\frac{A^{(1)}_{T}A^{(2)}_{T}}{\sigma_{\mathrm{inel}}}{F}(\tau).$
(9)

This distribution would scale if

$S_{T}={A^{(1)}_{T}A^{(2)}_{T}}/{\sigma_{\mathrm{inel}}}$ (10)

were energy-independent. This assumption leads to (2). It is the energy
dependence of $S_{T}$ that is responsible for different scaling properties of
the cross-section and multiplicity of charged hadrons in p+p collisions
Praszalowicz:2012zh . In the case of HIC for a fixed centrality, $S_{T}$ has
geometrical interpretation as an overlap area of two colliding nuclei
Kharzeev:2002ei ; Kharzeev:2000ph  that scales like $N_{\mathrm{part}}^{2/3}$.
Throughout this paper we shall assume $Q_{\rm s}\sim S_{T}\sim
N_{\mathrm{part}}^{\delta}$. Possible weak energy dependence of $S_{T}$ would
lead to the violation of geometrical scaling. For small systems (like $e.g.$
dA or pA) one has two saturation scales that differently scale with $N_{\rm
part}$ and the effective saturation scale has to be taken as
$\sqrt{Q_{\rm{sat}}[{\rm large}]Q_{\rm{sat}}\left[{\rm{small}}\right]}$ Klein-
Bosing:2014uaa ; Dumitru:2001ux . We address this issue in KhaPrasz . Here we
concentrate on large systems only.

The functional form of the saturation scale $Q_{\mathrm{{s}}}(x)$ in (4)
follows from the nonlinear QCD evolution Munier:2003vc , and in the case of
heavy ions, also from the collision geometry Kharzeev:2002ei ; Kharzeev:2000ph
, resulting in the following form of the scaling variable $\tau$:

$\tau=\frac{1}{N_{\text{part}}^{\delta/4}}\frac{p_{T}}{Q_{0}}\left(\frac{p_{T}}{W}\right)^{\lambda/2},$
(11)

where $W=\sqrt{s_{{}_{\\!NN}}}\times 10^{-3}$ is an effective energy scale
shown in (4).

Should there be any energy dependence of $S_{T}$ (which enters the definition
of $Q_{\rm s}$), it is already included in $\lambda$, as we fix $\lambda$ from
the fits to data that do not distinguish different sources (QCD nonlinear
evolution, $S_{T}$ energy dependence) of the energy increase of
$Q_{\mathrm{s}}(W)$. These fits KhaPrasz , as well as a direct inspection of
Fig. 3, indicate that for direct photons $\lambda\approx 0.2$, rather than
0.3. The same conclusion has been reached in an earlier study of Ref. Klein-
Bosing:2014uaa .

##  3 Charged particle multiplicity

Figure 4: Scaled multiplicity density (14) for charged particles at
$\lambda=0.2$ and $\lambda=0.35$ as a function of $N_{\rm part}$.

Let us first calculate the charged particle multiplicity for large collision
systems from the scaling formula (2). To this end we will need a Jacobian to
change the integration variable from $p_{T}$ to $\tau$:

$p_{T}dp_{T}=\frac{2Q_{0}^{2}}{2+\lambda}N_{\text{part}}^{\frac{\delta}{2+\lambda}}\left(\frac{W}{Q_{0}}\right)^{\frac{2\lambda}{2+\lambda}}\tau^{\frac{2-\lambda}{2+\lambda}}d\tau,$
(12)

which yields

$\displaystyle\frac{dN_{\text{ch}}}{d\eta}$
$\displaystyle={\displaystyle\int}d^{2}p_{T}\frac{dN_{\text{ch}}}{d^{2}p_{T}d\eta}$
$\displaystyle=S_{T}\,Q_{0}^{2}\,\frac{2}{2+\lambda}N_{\text{part}}^{\frac{\delta}{2+\lambda}}\left(\frac{W}{Q_{0}}\right)^{\frac{2\lambda}{2+\lambda}}{\displaystyle\int}F_{\text{ch}}(\tau)\,\tau^{\frac{2-\lambda}{2+\lambda}}d\tau$
$\displaystyle=N_{\text{part}}^{\frac{3+\lambda}{2+\lambda}\delta}\left(\frac{W}{Q_{0}}\right)^{\frac{2\lambda}{2+\lambda}}\kappa,$
(13)

where the constant $\kappa$ includes the integral of $F_{\text{ch}}$, $S_{T}$
and other irrelevant constants. Note that the energy dependence $\sim
W^{2\lambda/(2+\lambda)}$ is compatible with multiplicity growth in p+p as
measured at the LHC Aamodt:2009aa . For small $\lambda$ and $\delta=2/3$ we
have the quantity $N_{\text{part}}^{\frac{3+\lambda}{2+\lambda}\delta}\approx
N_{\text{part}}$, as has been assumed e.g. in Ref Kharzeev:2000ph , though it
has been now established that better scaling of charged multiplicities is
obtained when $N_{\rm part}$ is replaced by the number of participating quarks
(see below).

In order to assess the quality of the scaling formula (13), we plot the scaled
charged hadron multiplicity in Fig. 4:

$N_{\text{part}}^{-\frac{3+\lambda}{2+\lambda}\delta}\left(\frac{Q_{0}}{W}\right)^{\frac{2\lambda}{2+\lambda}}\,\frac{dN_{\text{ch}}}{d\eta}$
(14)

for the data from Refs. Adare:2015bua ; Aamodt:2010cz ; Adam:2015ptt ;
Acharya:2018hhy  as a function of $N_{\rm part}$ for $\delta=2/3$ and for two
different choices of $\lambda$. We see that, in contrast to the photon GS-
scaled yields, the best alignment of all data is achieved for $\lambda\sim
0.35$. Here $Q_{0}\sim 1$ GeV/$c$ as before. Some systematic shift is seen in
the case of Cu+Cu data. We would expect that the scaled multiplicity in (14)
should not depend on $N_{\rm part}$, while we see that it rises with $N_{\rm
part}$. This is a well-known problem in HIC, and it has been argued that
better scaling is obtained with the number of participating quarks
Adare:2015bua . Here it is enough to slightly change the value of $\delta$
from 2/3 to 3/4 to ensure $N_{\rm part}$ independence of (14). However, as we
shall see in the following, the value of $\delta$ does not enter the
conditions for the MS exponent $\alpha$.

The fact that $\lambda$ for charged particle spectra is different than the one
for photons is – to our best knowledge – not understood theoretically. It may
be due to the interactions of charged particles in medium. The other reason
might be that for charged particles we extract $\lambda$ from integrated
yields that include the scaling part up to $2\div 3$ GeV/$c$
Praszalowicz:2018vfw  and the non-scaling tail, while in obtaining Eq. (13) we
have assumed that GS is present in the whole experimentally measured range. In
the following we will accept this difference as a phenomenological
observation, and distinguish $\lambda$ for photons from $\lambda_{\rm ch}$ for
charged particles.

##  4 Relating the scaling laws

Let us rewrite the scaling laws (2) and (1) in the following way

$S_{T}\,F_{\gamma}(\tau(p_{T}))=\frac{dN_{\gamma}}{d^{2}p_{T}dy}=\left(\frac{dN_{\text{ch}}}{d\eta}\right)^{\alpha}\frac{1}{Q_{0}^{2}}G(p_{T}),$
(15)

which can be satisfied only when the left hand side and the right hand side of
Eq. (15) have the same $p_{T}$, $W$ and $N_{\text{part}}$ dependence.

It is clear that one cannot proceed further without explicit knowledge of the
functions $F_{\gamma}$ and $G$. To this end we shall assume a power-law
dependence

$F_{\gamma}(\tau)\sim\frac{1}{\tau^{n}}~{}~{}~{}\mathrm{and}~{}~{}~{}G(p_{T})\sim\left(\frac{Q_{0}}{p_{T}}\right)^{m}.$
(16)

Indeed, as shown in Figs. 2 and 3, where we plot $F_{\gamma}$ and $G$ for
Au+Au Adare:2008ab ; Afanasiev:2012dg ; Adare:2014fwh , Cu+Cu Adare:2018jsz
and Pb+Pb Adam:2015lda  direct photon data at different energies and
centralities for $\lambda=0.2$ and $\delta=2/3$, the power-law fall-off Klein-
Bosing:2014uaa  with $n,\,m\sim 4,\,4.5$ works pretty well at small and
intermediate $\tau$ (or $p_{T}$). Note that $m=4$ is a generic prediction for
radiation from the CGC Dumitru:2001ux  at intermediate $p_{T}$, while – as we
have checked by looking at prompt photons up to transverse momenta $\sim 2$
GeV$/c$ – fits to the data prefer $m=5$. Altogether

$\displaystyle\\!\\!\\!\frac{dN_{\gamma}}{d^{2}p_{T}dy}\sim
N_{\text{part}}^{\delta}\left(N_{\text{part}}^{\delta/4}\left(\frac{W}{Q_{0}}\right)^{\frac{\lambda}{2}}\left(\frac{Q_{0}}{p_{T}}\right)^{\frac{2+\lambda}{2}}\right)^{n}$
$\displaystyle~{}~{}~{}\sim\left(N_{\text{part}}^{\delta}N_{\text{part}}^{\frac{1}{2+\lambda_{\rm
ch}}\delta}\left(\frac{W}{Q_{0}}\right)^{\frac{2\lambda_{\rm
ch}}{2+\lambda_{\rm
ch}}}\right)^{\alpha}\left(\frac{Q_{0}}{p_{T}}\right)^{m},$ (17)

where we have used $S_{T}\sim N_{\mathrm{part}}^{\delta}$. The functions in
(17) are proportional to each other if

$m=\frac{2+\lambda}{2}\,n,~{}~{}\frac{4+n}{4}\,\delta=\frac{3+\lambda_{\rm
ch}}{2+\lambda_{\rm
ch}}\,\alpha\delta,~{}~{}\frac{\lambda}{2}\,n=\frac{2\lambda_{\rm
ch}}{2+\lambda_{\rm ch}}\,\alpha.$ (18)

From the first equation we see that the power-like fall-off should be faster
for MS than for GS, and this prediction works quite well, as can be seen from
Figs. 2 and 3.

Thereby, from (18) we obtain the following:

$\displaystyle\alpha=\left.\frac{4+n}{4}\,\frac{2+\lambda_{\rm
ch}}{3+\lambda_{\rm ch}}\right|_{\begin{subarray}{c}n=4\\\ \lambda_{\rm
ch}=0.35\end{subarray}}$ $\displaystyle=1.40,$ (19)
$\displaystyle\alpha=\left.n\,\frac{\lambda}{\lambda_{\rm
ch}}\frac{2+\lambda_{\rm ch}}{4}\right|_{\begin{subarray}{c}n=4\\\
\lambda=0.2\\\ \lambda_{\rm ch}=0.35\end{subarray}}$ $\displaystyle=1.34.$
(20)

Both estimates give the power $\alpha$ compatible with the experimental value
of 1.25 (given our crude assumptions and approximations). Note that the first
estimate follows from the assumptions concerning the geometry of A+A
collisions (although it does not depend on the actual value of $\delta$),
whereas the second one follows from the energy dependence.

##  5 Discussion and conclusions

Let us first observe that the equality of two estimates for $\alpha$, (19) and
(20), within an acceptable range of $n$, $\lambda$ and $\lambda_{\rm ch}$,
cannot be obtained if $\lambda=\lambda_{\rm ch}$. As mentioned in Sec. 3 there
is no theoretical understanding of this difference. One possibility might be
the apparent violation of GS that has been mentioned in Sec. 2, namely the
energy dependence of $S_{T}$. We do not have any experimental handle on this
dependence, however recent Glauber model calculations Loizides:2017ack  show
slight increase of $S_{T}$ with energy, depicted in Fig. 5, which can be
effectively parametrized as a power-law:

$S_{T}=\left(\frac{W}{Q_{0}}\right)^{\lambda^{\prime}}\frac{N_{\text{part}}^{\delta}}{Q_{0}^{2}}$
(21)

with $\lambda^{\prime}\sim 0.06$. While (21) may explain a difference in
quality of GS for multiplicity vs. cross-section, it is difficult to
disentangle this dependence from the genuine energy dependence of $Q_{s}$ on
$x$. This is because we effectively include the energy dependence of $S_{T}$
in the process of tuning $\lambda$ in search for GS. The naive inclusion of
(21) in Eqs. (15) and (13) would only modify the second equation for $\alpha$
in (20) shifting it to 1.29. However, a proper way of including the energy
dependence of $S_{T}$ would require a global fit with $\lambda$,
$\lambda^{\prime}$ and $\delta$ as free parameters, rather than adding
$\lambda^{\prime}$ in (20).

Figure 5:  Glauber model predictions Loizides:2017ack  for $S_{T}$ scaled by
$N_{\mathrm{part}}^{\delta}$, with $\delta=2/3$ for ALICE Pb+Pb 2760 GeV data
as upper (green) squares and PHENIX Au+Au 200 GeV data as lower (red) circles.
The shaded band corresponds to the PHENIX data scaled by $(2760/200)^{0.06}$,
including systematic uncertainties of $N_{\rm part}$.

There are many other sources that lead to the violation of GS in HIC, which we
have ignored. For example, the factorization of the unintegrated gluon
densities in (7) into $A_{T}$ and $\phi$, neglecting the momentum dependence
of $\alpha_{s}$ in (5), as well as simplifying the assumptions about $N_{\rm
part}$ dependence and the above-mentioned energy increase of $S_{T}$. It is
obvious from Fig. 5 that these types of dependence do not factorize, moreover,
much better scaling is obtained with the power $\delta=3/4$. The latter,
however, is of no importance for our analysis, as $\delta$ drops out from the
second equation of (18). Finally, simple power-law dependence of the photon
spectra is also a simplification that affects the accuracy of our analysis. A
choice of more accurate data-driven functions could give an improved value of
$\alpha$, closer to that extracted from experiment. This would require,
however, a numerical analysis that we wish to avoid for the clarity of
argumentation. One should not forget that both scaling laws observed
experimentally are only approximate. This concerns also the charged particle
spectra used here do calculate $dN_{\rm ch}/d\eta$ in (13).

Even though one might argue that $\sim$ 10% accuracy is quite satisfactory,
given a simplicity of the present analysis, our results may indicate that
there might be other components of direct photons that exhibit GS but have
different functional dependences than the initial stage component. This might
resolve the difference between $\lambda$ and $\lambda_{\rm ch}$.
Identification of such components would be of great interest.

Our conclusion is that we have linked the multiplicity scaling and geometrical
scaling of direct photon $p_{T}$-spectra, estimating the scaling power
$\alpha$ from simple assumptions based on the functional forms of the scaling
functions $G$ and $F_{\gamma}$, as well as on the energy dependence and
$N_{\rm part}$ dependence of the saturation scale. We have obtained two
independent, albeit consistent, estimates of $\alpha$ that overshoot
experimental value by $\sim$ 10%. This result relies on the fact that
$\lambda\neq\lambda_{\rm ch}$. We have argued that possible explanation of
this difference, apart from trivial facts that charged particles interact with
the medium while photons do not and that integrated yields of charged
particles include non-scaling tails, might also indicate that other
contributions of direct photon production exhibit GS but possibly have
different functional forms than the one of the early stage production.

## Acknowledgments

We are grateful to Axel Drees and Larry McLerran for many important and
fruitful discussions. The research of VK is supported under Department of
Energy Contract No. DE-FG02-96-ER40988. The research of MP is supported by the
NAWA (Polish National Agency for Academic Exchange) Bekker program and
partially by the COST Action CA15213 THOR. MP wishes to thank the Institute
for Nuclear Theory at the University of Washington for its kind hospitality,
stimulating research environment and for partial support by INT’s U.S.
Department of Energy Grant No. DE-FG02-00ER41132.

## References

  * (1) A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 104, 132301 (2010).
  * (2) S. Afanasiev et al. [PHENIX Collaboration], Phys. Rev. Lett. 109, 152302 (2012).
  * (3) A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 91, 064904 (2015).
  * (4) A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 98, 054902 (2018).
  * (5) A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 123, 022301 (2019).
  * (6) V. Khachatryan [PHENIX Collaboration], Nucl. Phys. A 982,763 (2019).
  * (7) V. Khachatryan [PHENIX Collaboration], Acta Phys. Pol. B Proc. Suppl. 12, 457 (2019).
  * (8) J. Adam et al. [ALICE Collaboration], Phys. Lett. B 754, 235 (2016).
  * (9) L. Adamczyk et al. [STAR Collaboration], Phys. Lett. B 770, 451 (2017).
  * (10) M. M. Aggarwal et al. [WA98 Collaboration], Phys. Rev. Lett. 85, 3595 (2000).
  * (11) A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 94, 064901 (2016).
  * (12) S. Acharya et al. [ALICE Collaboration], Phys. Lett. B 789, 308 (2019).
  * (13) M. Dion, J.-F. Paquet, B. Schenke, C. Young, S. Jeon and C. Gale, Phys. Rev. C 84, 064901 (2011).
  * (14) C. Shen, U. W. Heinz, J. F. Paquet, I. Kozlov and C. Gale, Phys. Rev. C 91, 024908 (2015).
  * (15) C. Shen, U. Heinz, J.-F. Paquet, and C. Gale, Phys. Rev. C 89, 044910 (2014).
  * (16) J. F. Paquet, C. Shen, G. S. Denicol, M. Luzum, B. Schenke, S. Jeon and C. Gale, Phys. Rev. C 93, no. 4, 044906 (2016).
  * (17) H. van Hees, C. Gale, and R. Rapp, Phys. Rev. C 84, 054906 (2011).
  * (18) R. Rapp, H. van Hees, M. He, Nucl. Phys. A 931, 696 (2014).
  * (19) H. van Hees, M. He, R. Rapp, Nucl. Phys. A 933, 256 (2015).
  * (20) E. L. Bratkovskaya, S. M. Kiselev and G. B. Sharkov, Phys. Rev. C 78, 034905 (2008).
  * (21) E. L. Bratkovskaya, Nucl. Phys. A 931, 194 (2014).
  * (22) O. Linnyk, W. Cassing, E. Bratkovskaya, Phys. Rev. C 89, 034908 (2014).
  * (23) O. Linnyk, V. Konchakovski, T. Steinert, W. Cassing, E. L. Bratkovskaya, Phys. Rev. C 92, 054914 (2015).
  * (24) K. Dusling and I. Zahed, Phys. Rev. C 82, 054909 (2010).
  * (25) K. Dusling, “Hydrodynamic Description of Dilepton Production”, PhD thesis, Stony Brook University (2008).
  * (26) C.-H. Lee and I. Zahed, Phys. Rev. C 90, 025204 (2014).
  * (27) Y.-M. Kim, C.-H. Lee, D. Teaney and I. Zahed, Phys. Rev. C 96, 015201 (2017).
  * (28) C. Shen, Nucl. Phys. A 956, 184 (2016).
  * (29) G. David, Rept. Prog. Phys. 83, 046301 (2020).
  * (30) M. Chiu, T. K. Hemmick, V. Khachatryan, A. Leonidov, J. Liao and L. McLerran, Nucl. Phys. A 900, 16 (2013).
  * (31) C. Klein-Bösing and L. McLerran, Phys. Lett. B 734 (2014) 282
  * (32) L. McLerran, Acta Phys. Polon. B 45, no. 12, 2307 (2014).
  * (33) L. McLerran and B. Schenke, Nucl. Phys. A 929, 71 (2014).
  * (34) L. McLerran and B. Schenke, Nucl. Phys. A 946, 158 (2016).
  * (35) J. Berges, K. Reygers, N. Tanji and R. Venugopalan, Phys. Rev. C 95, 054904 (2017).
  * (36) V. Khachatryan, B. Schenke, M. Chiu, A. Drees, T. K. Hemmick and N. Novitzky, Nucl. Phys. A 978, 123 (2018).
  * (37) V. Khachatryan and M. Praszalowicz, in preparation.
  * (38) L. McLerran and M. Praszalowicz, Acta Phys. Polon. B 41, 1917 (2010).
  * (39) L. McLerran and M. Praszalowicz, Acta Phys. Polon. B 42, 99 (2011).
  * (40) M. Praszalowicz, Acta Phys. Polon. B 42, 1557 (2011).
  * (41) M. Praszalowicz, EPJ Web Conf. 206, 02002 (2019).
  * (42) A. M. Stasto, K. J. Golec-Biernat and J. Kwiecinski, Phys. Rev. Lett. 86, 596 (2001).
  * (43) L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100, 1 (1983).
  * (44) A. H. Mueller and J. W. Qiu, Nucl. Phys. B 268, 427 (1986).
  * (45) I. Balitsky, Nucl. Phys. B 463, 99 (1996).
  * (46) Y. V. Kovchegov, Phys. Rev. D 60, 034008 (1999).
  * (47) Y. V. Kovchegov, Phys. Rev. D 61, 074018 (2000).
  * (48) A. H. Mueller and D. N. Triantafyllopoulos, Nucl. Phys. B 640, 331 (2002).
  * (49) S. Munier, R. B. Peschanski, Phys. Rev. Lett. 91, 232001 (2003).
  * (50) S. Munier, R. B. Peschanski, Phys. Rev. D 69, 034008 (2004).
  * (51) L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 2233 (1994).
  * (52) L. D. McLerran and R. Venugopalan, D 49, 3352 (1994).
  * (53) L. D. McLerran and R. Venugopalan, D 50, 2225 (1994).
  * (54) F. Gelis, E. Iancu, J. Jalilian-Marian and R. Venugopalan, Ann. Rev. Nucl. Part. Sci. 60, 463 (2010).
  * (55) L. McLerran, Acta Phys. Pol. B 41, 2799 (2010).
  * (56) F. Caola and S. Forte, Phys. Rev. Lett. 101, 022001 (2008).
  * (57) E. Iancu, K. Itakura and L. McLerran, Nucl. Phys. A 708, 327 (2002).
  * (58) L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Lett. B 100, 173 (1981).
  * (59) K.J. Golec-Biernat and M. Wüsthoff, Phys. Rev. D 59, 014017 (1998).
  * (60) K.J. Golec-Biernat and M. Wüsthoff, Phys. Rev. D 60, 114023 (1999).
  * (61) D. Kharzeev, E. Levin and M. Nardi, Nucl. Phys. A 730, 448 (2004). Erratum: [Nucl. Phys. A 743, 329 (2004)].
  * (62) A. Szczurek, Acta Phys. Pol. B 35, 161 (2004).
  * (63) H. Kowalski and D. Teaney, Phys. Rev. D 68, 114005 (2003).
  * (64) P. Tribedy and R. Venugopalan, Nucl. Phys. A 850, 136 (2011) Erratum: [Nucl. Phys. A 859, 185 (2011)].
  * (65) M. Praszalowicz and A. Francuz, Phys. Rev. D 92, 074036 (2015).
  * (66) M. Praszalowicz and T. Stebel, JHEP 1303, 090 (2013).
  * (67) A. Dumitru and L. D. McLerran, Nucl. Phys. A 700, 492 (2002).
  * (68) D. Kharzeev and M. Nardi, Phys. Lett. B 507, 121 (2001).
  * (69) K. Aamodt et al. [ALICE Collaboration], Eur. Phys. J. C 65, 111 (2010).
  * (70) A. Adare et al. [PHENIX Collaboration], Phys. Rev. C 93, 024901 (2016).
  * (71) K. Aamodt et al. [ALICE Collaboration], Phys. Rev. Lett. 106, 032301 (2011).
  * (72) J. Adam et al. [ALICE Collaboration], Phys. Rev. Lett. 116, 222302 (2016).
  * (73) S. Acharya et al. [ALICE Collaboration], Phys. Lett. B 790, 35 (2019).
  * (74) C. Loizides, J. Kamin and D. d’Enterria, Phys. Rev. C 97, no. 5, 054910 (2018) Erratum: [Phys. Rev. C 99, no. 1, 019901 (2019)].


# Formal and Practical Elements
for the Certification of Machine Learning Systems

Jean-Guillaume Durand1 and Arthur Dubois2  1Head of Perception<EMAIL_ADDRESS>
2Director of Engineering<EMAIL_ADDRESS>
Xwing, San Francisco, CA 94102  Robert J. Moss  PhD Student, Department of
Computer Science
<EMAIL_ADDRESS>
Stanford University, Stanford, CA 94305

###### Abstract

Over the past decade, machine learning has demonstrated impressive results,
often surpassing human capabilities in sensing tasks relevant to autonomous
flight. Unlike traditional aerospace software, the parameters of machine
learning models are not hand-coded nor derived from physics but learned from
data. They are automatically adjusted during a training phase, and their
values do not usually correspond to physical requirements. As a result,
requirements cannot be directly traced to lines of code, hindering the current
bottom-up aerospace certification paradigm. This paper attempts to address
this gap by 1) demystifying the inner workings and processes to build machine
learning models, 2) formally establishing theoretical guarantees given by
those processes, and 3) complementing these formal elements with practical
considerations to develop a complete certification argument for safety-
critical machine learning systems. Based on a scalable statistical verifier,
our proposed framework is model-agnostic and tool-independent, making it
adaptable to many use cases in the industry. We demonstrate results on a
widespread application in autonomous flight: vision-based landing.

††©2023 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

Figure 1: Overview of the proposed framework.

##  I Introduction

Machine learning (ML) is a long-established field that builds input-output
models by automatically learning from example data. Over the past decade, a
succession of seminal research breakthroughs [1] rekindled the field and
unlocked a wide range of new applications. ††That progress was largely
possible due to increased amounts of available data and computational power.
Deep learning is a subset of machine learning, with a subtle difference. In
machine learning, the model learns from data features which are handcrafted by
subject matter experts (i.e., object color, edges, texture, velocity, etc.).
In deep learning, the model learns these features on its own. For many
specialized tasks, such models have outperformed traditional methods, and even
human specialists. These tasks range from medicine diagnosis to satellite
imagery analysis, disease classification, board and video games strategies,
and many others. Researchers and practitioners regard deep learning as a key
enabler for autonomy in robotics, allowing unprecedented levels of safety. In
particular, deep learning has been behind the autonomous driving successes of
the past decade, and has the potential to unlock the true benefits of
autonomous flight [2]. In spite of these achievements, the widespread adoption
of machine learning has been slow for safety-critical systems such as cars and
aircraft. This is largely because of a lag in regulations and the absence of
standardized certification policies, especially in aerospace.

Figure 2: ARP4754A framework [3]

In aerospace, certification aims at ensuring that hardware, software, and
processes follow a defined set of design and performance requirements, and
ultimately have a very low probability of failure. The key to obtaining
operational approval lies in quantitatively demonstrating reliability over an
entire Operational Domain (OD). For aerospace, government agencies grant
certification, such as the Federal Aviation Administration (FAA) in the US,
and the European Union Aviation Safety Agency (EASA) in Europe. They do not
set exact practices to follow but rather issue advisory circulars to recognize
acceptable means for developing and certifying aircraft, aerospace systems, or
components.

These agencies are supported by professional associations such as SAE
International, RTCA Inc., or the European Organisation for Civil Aviation
Equipment (EUROCAE). They provide a forum to devise technical standards and
recommended practices for the design of aircraft systems. These documents do
not carry any legal force but are often referenced as acceptable means of
compliance. The regulatory organization then works with each applicant to
approve the use of each standard to define the formal Means of Compliance
(MOC).

Figure 3: ARP4754A V-lifecycle [3]

The widely accepted paradigm for aerospace is to follow guidance from SAE
ARP4761 [4] and ARP4754A / ED-79A [3] (Figure 2). This framework follows a
classic systems engineering V-lifecycle (Figure 3) to determine the level of
certification required by system components. Using Functional Hazard Analysis
(FHA), the system architecture is decomposed into different levels; from the
aircraft level, down to subsystems. The functions of each level are then
listed out and mapped to respective subsystems. The failure conditions and
severity of each function then determine the Development Assurance Level (DAL)
and establish the safety requirements to meet (Figure 4). The applicant then
follows a Means of Compliance (MOC) to verify that the failure condition
requirements are met.

Figure 4: Assurance levels based on failure severity and probability [3], p.14

DO-254 / ED-80 [5] and DO-178C / ED-12C [6] are companion documents which
provide guidance for the certification of hardware and software components,
respectively. Similarly, DO-200B [7] provides guidelines for handling
aeronautical data and could be extended to the management of machine learning
data. For the scope of this paper, we focus on the software certification
aspect. Following DO-178C, each requirement must be traceable to multiple
artifacts: the lines of code implementing it, the test cases verifying the
correctness of that code, and the results of such tests. This ensures that
each requirement is fulfilled by the source code, and that there is no
extraneous code.

With machine learning models, mapping lines of code to requirements is not
straightforward, hence creating a gap in the current certification process.
First, such models often contain millions of parameters, and subsets of those
are activated depending on the properties of the input. This makes it
impossible to trace the impact of a single parameter on the resulting output.
Second, those parameters are not directly derived from physics or requirements
as is the case with traditional software (e.g., a Kalman filter, a yaw
controller). The parameters are automatically learned from example data which
is fed to the model during its training phase. Thus in order to trace this
parameter to a requirement, that data itself should be mapped to the
requirements. This shifts the focus of current certification processes from
not only the code itself, but also towards data traceability. Third, some
machine learning models build intermediate representations of the data that
are increasingly more complex, powerful, but also not interpretable. By
looking at those intermediate layers in the code, it is not possible to map
their parameters to requirements.

The previous observations prevent the bidirectional mapping that is required
from requirements to code and parameters: additional means of compliance and
standards are needed. The next section describes attempts from key research
and industry actors to fill this gap.

##  II Related Work

###  II-A Certification frameworks

In their 2016 report on the verification of adaptive systems [8], the FAA
first established some useful terminology, and drew the boundary between
adaptive and non-adaptive systems. For the scope of this paper, we are
concerned only with non-adaptive systems which do not change behavior in
operation.

EASA then proposed an incremental Artificial Intelligence (AI) roadmap [9, 10]
focused on trustworthiness and the need for social acceptance. They propose
ramping up capabilities from simpler crew assistance (2022-2025), to
human/machine collaboration (2025-2030), and full autonomy (2035+). They
identify a paradigm shift to focus the assurance process on the correctness
and completeness of datasets, bias mitigation, model accuracy and performance,
as well as novel verification methods.

Figure 5: W-shaped cycle for learning assurance [11]

EASA and Daedalean AG published the joint CoDANN report [11] to study the
learning assurance concept introduced in the EASA roadmap. They present a
novel W (double V) assurance cycle (Figure 5). In the first V, the model is
trained and verified iteratively until a satisfactory performance is met. In
the second V, the traditional verification process continues with the model
implementation. For dataset validation, they present a distribution
discriminator to evaluate the completeness of datasets. The report establishes
a distinction between robustness to training data (algorithm robustness), and
to input perturbations (model robustness).

Shortly after, CoDANN II [12] addressed additional topics such as
implementation, inference, and model explainability. CoDANN I and II are
demonstrated on vision-based landing and Detect-and-Avoid (DAA) applications
[13], attesting to the modularity of the proposed framework. Following this
work, EASA then published usable guidance for AI applications of level 1 & 2
(i.e. human assistance), establishing clear objectives for applicants [14,
15].

Xwing Inc. proposed a modified V-diagram [16] in three stages: 1) data
requirements and processes to verify data used to build and test the model, 2)
model performance, and 3) reliability of operation at runtime (Figure 6).

Figure 6: Modified V-diagram from [16]

Xwing presented results on a DAA application and compared a traditional radar-
based approach with a computer vision method using machine learning. Their
paper identifies the clear difference between the traditional physics-based
bottom up approach, and the top-down approach required for machine learning
applications. The authors emphasize that model parameters are just as
important, if not more, than the code processing the network. They also
propose to treat training data identically to aeronautical design data in
traditional certification. However, the paper does not detail how to validate
a machine learning model against performance requirements. We propose to
extend this framework in our solution.

We can also note the contribution of working groups on the topic such as SAE
G-34 and AFE 87. These groups are currently a work-in-progress and their
initial findings echo the research aforementioned.

###  II-B Machine learning verification

####  II-B1 Verification of Neural Networks (VNN)

A major failure mode of ML models is local instability: a small input
perturbation can result in large output variations. Adversarial attacks
exploit this vulnerability to trick models into giving incorrect outputs. VNN
gives formal guarantees that a model is not susceptible to any attack within a
given input range [17]. While it is unlikely to adversely affect image
formation at an airport or airspace level, robustness to small input
variations should still be considered as they can occur in nominal operation.

Formal verifiers such as VENUS [18, 19], AI2 [20] or VeriNet [21, 22] solve
verification problems using a branch-and-bound approach. Counterexamples of
unexpected behavior are obtained as a byproduct. However formal verifiers are
slow and limited to small feed-forward architectures with ReLU activation
functions, which fail to scale up to real-world aerospace models.

Reluplex [23], Marabou [24], DeepTest [25], DeepZ [26], DeepSafe [27], and
DeepHunter [28] all require knowledge of the model activations or architecture
and are thus not scalable. Some of these approaches were successfully applied
to the ACAS Xu aircraft collision avoidance algorithm [23, 24]. Most methods
focus on pixel-space noise which can model realistic sensor errors, but fails
to capture semantic errors (e.g., brightness, saturation, scaling, etc.).
Semantify-NN [29] and Tiler [30] provide robustness to families of semantic
perturbations but modify the training process or require a full image
formation model, which may prove impractical.

####  II-B2 Explainability

A large body of work focuses on making neural networks interpretable by
peeking into the hidden layers [31]. While core to the human-centric EASA
roadmap, we consider it intractable due to the large and ever-changing number
of model architectures, and possible variants. Model-agnostic interpretation
techniques exist but still face many challenges [32]. They do not provide the
level of interpretability expected in traditional software certification but
could still be useful for development purposes or post-mortem analyses: we
will thus favor black-box verification.

####  II-B3 Black-box verification

Black-box systems are opaque to the user, only the input and output can be
observed, discarding the inner workings of the model. Corso et al. [33]
standardizes the problem formulation and studies applications to Autonomous
Vehicles (AV) and aircraft collision avoidance. Adaptive stress testing (AST)
[34] uses reinforcement learning to find likely system failure paths in
simulation. The approach is limited to sequential decision systems and has not
been generalized to AI models which process single images. By ranking failure
scenarios, the simulation testing framework from [35] performed the first
independent black-box evaluation of a commercial AV system.

The FAA considers the characteristics of training data as necessary conditions
for verification [36]. Focusing on a black-box approach to neural networks,
they propose the well-known $R^{2}$ indicator to determine the goodness-of-fit
of a given model, and introduce the $S^{2}$ indicator to determine if the
training set adequately covers the model operating range. While limited to a
single neuron with no activation function, additional efforts are underway to
extend the approach.

Failures can be rare and challenging to estimate in real flights (e.g.,
prohibitive number of tests [37] and dangerous scenarios), or in simulation
(computationally expensive). NASA thus introduced SysAI [38] to efficiently
generate simulated test cases using active learning and DynaTrees. This
approach effectively uses less samples than Monte Carlo but requires
parameterized boundary shapes for model failures, and is tailored to
sequential data. Moss et al. [39] then proposed a statistical verifier which
scales better than formal verifiers. It uses Bayesian optimization and
importance sampling to efficiently sample the Operational Domain (OD) as a
parametric space. Using this minimal set of design points, a surrogate model
is fitted for an inexpensive but robust estimate of the system failure
probability. Results are successfully demonstrated on a simulated vision-based
landing application.

##  III Proposed approach

We choose not to open the AI black box. Indeed, certifying individual
architectures would be ill-fated as their sheer number is intractable;
architectures evolve rapidly in this very dynamic research field. Even the
handful of well-established architectures are customized to optimize
performance in use cases that differ for each applicant. Alternatively, a
restriction to interpretable-only models would limit the capacity and
performance of models that could be used, ultimately preventing the safety
improvements offered by machine learning. As for explainability, a deeper
understanding of the inner layers of a model might be useful for development
or debugging, but is unnecessary for the certification argument where the
holistic performance of the system is evaluated. Similar to hardware
qualification, this black-box approach is a key differentiator and shifts the
paradigm from requirement-based code, to performance-based code.

Whether it takes the form of a W-diagram [11] or a modified V-diagram [16],
the research community has converged on three-stage verification: dataset,
model, and inference. We follow this decomposition and start by establishing
performance requirements for the ML subsystem.

###  III-A Requirements development

We follow the traditional ARP4761 and ARP4754 standards in place. The
applicant should first define a complete and correct Operational Design Domain
(ODD) of the system: an exhaustive list of conditions in which the model will
operate, and their expected proportions. For autonomous aircraft, examples can
be flight location, weather conditions, time of day, runway type, etc.

Figure 7: Notional functional hazard analysis with ML function highlighted
(yellow)

The traditional FHA and Fault Tree Analysis (FTA) are then carried out. As a
result, the ML function is assigned a DAL level, and an associated probability
of failure $P_{\text{fail}}$ (Figure 7). EASA recognizes that estimating
$P_{\text{fail}}$ is necessary but not sufficient on its own, and additional
model performance artifacts are needed for certification (see [15]). For the
scope of this paper, we present a methodology to substantiate the
$P_{\text{fail}}$ computation. The functional decomposition above helps
isolating the behavior of the ML component and sets its context. As can be
seen in Figure 11, it is a small function of a more comprehensive system and
does not operate the controls of the aircraft directly. The output of the
model is bounded and monitored by safety layers and other non-ML subsystems
which follow the traditional certification process.

The problem is then reduced to two things: 1) following software development
processes consistent with the DAL level, and 2) computing the probability of
failure for the ML system over the ODD (Figure 1).

###  III-B Dataset validation and verification

The data is split between independent training, validation, and test sets
which identically represent the ODD (i.e. drawn from the same distribution).
The training and validation sets are used by developers to optimize
respectively weights and hyperparameters of the model. The test set or
certification set is kept separate by a quality or certification team. Never
seen by developers, this set ensures that good model performance actually
stems from unbiased generalization on unseen data, rather than overfitting the
model to the training/validation data.

These datasets evolve continuously: applicants might add additional data
focusing on model failures, correct ground truth labels, or remove outdated
data (e.g., new sensor). For complete traceability, the applicant should track
these changes as Engineering Change Orders (ECO) in a Product Lifecycle
Management (PLM) system. This version control applies to the data itself, but
also to the associated metadata: ground truth labels, timestamp,
geolocalization, aircraft attitude, and other relevant descriptions (Figure
8). Those attributes can then be used to analyze the completeness of the
operational domain. The traceability of the certification set is especially
important: the applicant should validate that no examples from that set were
used in developing the model. By versioning the changes, it is also always
possible to reproduce previous models. We recommend the same traceability for
open-source pre-trained models, usually developed by large corporations on
very expensive cloud infrastructure. In theory, this would enable reproducing
a given model. Synthetic data should follow the same process and any simulator
should be validated through flight test in order to minimize any simulation-
to-reality gap (i.e., the sim2real problem).

Figure 8: Data version control. Data and metadata (here images and labels) are
hashed into a repository to keep track of their integrity. Each datum can be
assigned multiple labels (or no label). Versioned data and metadata can then
be used to assemble datasets for training or testing algorithms.

Finally, dataset verification also ensures that acceptable performance on the
test set translates to acceptable performance in the operational domain. As
such, the test set should capture the variations observed in the ODD, and span
the parameter ranges that the applicant defined in their parametric
decomposition. Rare corner cases should be also included to test the
robustness of the model and weighted according to their operational
likelihood. To verify the coverage of the test dataset, we propose to use the
$R^{2}$ and $S^{2}$ indicators [36], and the distribution discriminator [11].

###  III-C Training verification

We focus on parametric supervised learning models which are deterministic: for
a given input, the model will always yield the same output. We do not consider
online learning because the model would keep learning after being deployed
into operation, and inference would run on an unverified model.

Model iterations are often required to keep on improving a baseline
performance. Each experiment might use different data, hyperparameters, model
architecture, or training process. When a model is trained, a number of
artifacts need to be tracked in order to ensure experiments are reproducible
(Figure 1):

  * •

The data and labels used to train the model (including training/validation
split)

  * •

Model architecture, hyperparameters, and weights

  * •

Training hardware and software (including versioned dependencies)

  * •

Model initialization weights (before training)

  * •

Seeds used for pseudo-random number generators

  * •

Training logs and metrics

For certification, the applicant should implement traceability of at least:
the artifacts necessary to instantiate the model (architecture,
hyperparameters, weights), and the certification dataset. Akin to hardware
certification, we argue that these two artifacts alone are sufficient to build
the model and verify that it passes the robustness and performance tests on
the ODD. Many different models could be designed from the same
training/validation data. Following our black-box approach, we argue that
traceability of other artifacts and training/validation datasets is a healthy
and encouraged internal practice when developing the models, but might not be
as critical for certification substantiation.

###  III-D Inference verification

Once a suitable model is obtained in the development phase, it has to be
ported as executable code on the target hardware. This porting exercise is
usually necessary because aircraft computers differ in requirements (e.g.,
real-time constraints, memory, etc.) from the development computers used to
train the models. This model conversion can slightly modify the original
research model (layers grouping, weight quantization, pruning, hardware
threading determinism, etc.) and impact detection performance. Inference
software and hardware can be developed with traditional means of compliance
such as DO-178C and DO-254. The ported model can now be considered as a black-
box system to verify. Inference is verified in terms of performance and
operational runtime guarantees.

For performance, we need to compute the probability of failure
$P_{\text{fail}}$ (see section III-A). To this end, we center our framework
around the Bayesian safety validation (BSV) algorithm [39] to efficiently
sample the ODD and obtain an accurate estimate of $P_{\text{fail}}$. The
applicant should define metrics that evaluate the performance of the model
over its expected tasks and operational domain. For machine learning models,
some common metrics include: mean average precision (mAP), precision and
recall, false positive rate (FPR), mean squared error (MSE), mean absolute
error (MAE), etc. These metrics depend on the type of task (e.g.,
classification, regression, tracking) and should be defined by the applicant.

As for runtime requirements, they can be verified and validated using
traditional methods. For machine learning models, this step ensures that the
model obtained in the research phase is correctly ported to the target
hardware. Again, the traceability of the model architecture, weights, and
hyperparameters is essential to verify this step. Other aspects to be verified
include: throughput, latency, execution errors, determinism. It is important
to not that these requirements do not apply to the development/training code.

##  IV Results

This section introduces the use case of vision-based landing, before showing
how we apply the proposed approach.

###  IV-A Vision-based landing

Figure 9: Vision-based landing using a machine learning component to detect a
runway. Runway corners and sidelines are detected in blue, defining a coarse
bounding box in red.

A common application in autonomous flight, Vision-Based Landing (VBL) (Figure
9) uses only camera images to augment the landing guidance obtained from GPS
and/or Instrument Landing System (ILS). This is a typical example of computer
vision application that was enabled by the resurgence of deep learning.
Looking at the high-level architecture of a notional autonomous flight system
(Figure 10), we find interconnected components from the high-level mission
management system, down to the flight controls. The low-level flight control
components help make the system independently safe [40]: they maintain the
flying state of the aircraft and its safety, before taking into account other
inputs from high-level components.

Figure 10: Notional autonomous flight components. The vision-based subsystem
is a small part of the comprehensive system.

One of the systems is in charge of vertical and lateral guidance to the
runway. In that subsystem, vision-based landing provides one solution, along
with ILS and GPS solutions. In the VBL block, an ML component is used to
detect the runway in the image, while other traditional computer vision
algorithms perform tracking across frames, and pose estimation (Figure 11).
This separation of concerns helps understand where the ML model fits in the
overall system, and visualize how it constitutes only a small portion of the
complete solution.

Figure 11: Notional VBL architecture. The machine learning model is used only
for one step of the overall algorithm.

The input to the VBL system is an image with the associated sensor calibration
(intrinsics and extrinsics), the main output is the position and orientation
of the aircraft with respect to the runways detected in the image. The
different components can be described as follows:

  * •

Camera driver: Acquires images from the sensor and feeds them to the rest of
the pipeline.

  * •

Runway detector: The ML component. A model trained to detect runways in
images. For each runway detected in the image, the model outputs the four
runway corners position in the image, as well as uncertainties on the
prediction. Many convolutional neural network (CNN) architectures can
implement this task: single-shot detectors, region-proposal networks,
anchorless networks, or more recently transformers. Note that machine learning
algorithms other than CNNs can also be used such as Histogram of Gradients
(HoG) combined with Support Vector Machines (SVM) or Random Forests.

  * •

Pose estimator: Uses the corners of the runway and camera sensor calibration
to compute the position and orientation (i.e., pose) of the camera (and
aircraft) with respect to the runway.

  * •

Tracker: Uses traditional a computer vision approach (i.e., non ML-based) to
track the corners of the runway across frames and provide robust position
estimates.

  * •

Health monitor: Runs consistency checks on the output of the previous steps:
Is the prediction confidence high enough? Is this the expected runway? Are the
runway corners geometrically consistent? Is there a large jump in position
estimate?

Note that the output of the machine learning model is gated and validated in
many ways: by its own confidence score, by the geometric consistency of its
output, by the VBL health monitor, by comparison with other estimators (ILS,
GPS+INS), by the mission management system, and finally the vehicle management
system. Additionally, each system has redundancies and consistency checks that
can be run on the output of several VBL algorithms running in parallel.

###  IV-B Verification

####  IV-B1 Requirements development

We carry out the FHA and FTA of the system and its components (Figure 11), and
arrive at a DAL level C with an associated failure probability of $10^{-4}$
per approach over the operational domain (Figure 7), for the runway detector
model.

Given that the runway detector is only active during the approach phase of the
flight, we condition our probability distribution on the approach:
$P(\text{fail}\mid\text{approach})$. We then parameterize our operational
domain along known dimensions of variability. For pedagogic and plotting
purposes, we decompose our ODD along only two continuous dimensions: glide
slope angle $\alpha\sim\mathcal{N}\left(3,0.3,\left[1,7\right]\right)$
degrees, and distance to the runway
$d\sim\mathcal{N}\left(0,1.5,\left[0,4\right]\right)$ nautical miles. Note
that $\mathcal{N}\left(\mu,\sigma,\left[a,b\right]\right)$ represents the
truncated normal distribution with mean $\mu$, standard deviation $\sigma$,
and truncation interval $\left[a,b\right]$. Additional dimensions are
presented on Figure 12 for example.

Figure 12: Some example ODD dimensions after a draw of 10,000 samples (roll
angle, elevation, location, time of day)

####  IV-B2 Dataset verification

We use data version control to track changes to our data and metadata.
Datasets are assembled from versioned data and metadata and used to train and
test the models. The integrity of the data is guaranteed by computing a hash,
a unique digital signature that depends on the content of the data itself. If
the data is edited, or corrupted in any way, the hash also changes.

Runway images are labeled in batches through third-party providers. Annotators
are given precise instructions on how to label a runway and each batch is
subject to quality analysis on a random sample. A runway label consists of
annotating the full extent of visible runway edges, from which we can derive
runway corners for example (Figure 9).

We augment our dataset by using synthetic images generated by a physics-based
rendering engine (Figure 13).

Figure 13: Synthetic images covering different conditions of the parametric
space (sun in the field of view, weather, dangerous aircraft position and
orientation, etc.)

For independence of the different sets, we use a special splitting policy so
that no two images from the same landing sequence end up in both the training
and the certification set. In addition to version control, we also propose a
strict data access policy. We propose that every data collection campaign be
gated by an internal certification team or in collaboration with Designated
Engineering Representatives (DER). A given amount of data is kept to
contribute to the certification test set, while the rest is released to the
engineering team to develop the models. Keeping this separation ensures that
no practitioner can cheat and boost model performance by accessing the
certification test data while developing the model. As part of the development
cycle, the model is then provided to the certification team or DERs to be run
on the certification data, and the evaluation results are returned to the
engineering team for approval or further developments.

We propose healthy operational overfitting by training on all the runways
where the aircraft is expected to be deployed. This includes 65 airports where
Xwing Air111Xwing Air is the cargo airline of Xwing, Inc. Under a Part 135
certificate, their fleet of 35 Cessna Caravans operates more than 400 weekly
flights to shuttle packages across 34 UPS delivery routes. operates (Figure
12). This set is extended to also include alternate airports in case of
diversion. This healty overfitting ensures that any runway encountered in
operation has been seen by the model in training, and should result in good
detection performance. This set represents 80% of the total development data
and is completed by additional airports and runways to promote performance
generalization further. Note that the test set or certification set consist
only of real-world data from the operational domain and is not accessed by
developers during the training phase. Using a common 90/10% training/test
split, a breakdown is presented in Table I.

####  IV-B3 Model verification

We use the PyTorch deep learning framework in deterministic mode to develop
our models. For every training session, we take a snapshot of the code, its
software dependencies, and the hardware configuration. This ensures that the
training is fully deterministic and repeatable.

Split | Flights | Images
---|---|---
Development |  |
Training (real) | $45$ |  $23{,}442$ ($65\%$)
Training (synthetic) | N/A |  $5{,}000$ ($14\%$)
Validation (real) | $10$ |  $4{,}067$ ($11\%$)
Certification |  |
Test (real) | $9$ |  $3{,}488$ ($10\%$)
Total | $65$ |  $35{,}997$ ($100\%$)
TABLE I: Dataset composition

We use transfer learning (i.e. pre-trained weights) to leverage proven state-
of-the-art models for single-class (i.e. runway) object detection. Our model
stems from a Keypoint R-CNN architecture with a ResNet-50-FPN backbone [41].
The ResNet backbone acts as a task-agnostic feature extractor and is pre-
trained on the very large corpus of the ImageNet dataset [42]. The keypoint
detector is pre-trained on the COCO 2017 Keypoint dataset [43]. The model
outputs a bounding box around the runway, a confidence score, as well as
keypoints of interest (the corners of the runway). The loss minimized during
training is a combination of bounding box regression loss (position, size),
and Object Keypoint Similarity (OKS).

Figure 14: Learning curves plotting the loss function per optimizer iteration
(mini-batch). The checkpoint around iteration $50{,}000$ achieves best
generalization performance on the validation set.

We minimize the loss function with a Stochastic Gradient Descent (SGD)
optimizer with a momentum of $0.9$ and a maximum of $90{,}000$ iterations
(Figure 14). We feed $4$ images per mini-batch with a base learning rate
$lr=0.02$. For warmup and scheduling, we decay $lr$ by $\gamma=0.1$ after
$30{,}000$, $60{,}000$, and $80{,}000$ iterations. The model is trained with
on-premise deep learning infrastructure. During training, we apply random
image augmentations (e.g., blur, color jitter, rotation, scale, fog, rain,
snow, sunflare, etc.) to prevent overfitting and improve the robustness of the
network. At the end of model development, we choose the model around iteration
$50{,}000$ which minimizes the validation loss; using early stopping as a form
of regularization to minimize overfitting (Figure 14).

####  IV-B4 Inference verification

Confidence scores below $0.7$ are discarded when computing predictions. We
track two main metrics for the runtime performance of the model:

  * •

$AP_{BB}$: The average bounding box precision. Computed for Intersection Over
Union $IoU\in\left[0.50,0.95\right]$ for all object sizes. We limit the number
of model detections to a maximum of $100$.

  * •

$AP_{KP}$: The average keypoints precision. Computed for Object Keypoint
Similarity [43] $OKS\in\left[0.50,0.95\right]$ for all object sizes. We limit
the number of model detections to a maximum of $20$.

After a few development iterations and hyperparameter adjustments, the best
and final model achieves $AP_{BB}=0.89$ and $AP_{KP}=0.99$ on our validation
set. We convert the model from the developmental PyTorch format to an
inference executable using runtime libraries such as NVIDIA TensorRT, Apache
TVM, or ONNX. The inference engine compiles the model into an optimized set of
instructions compatible with GPUs. By tracing model architecture, weights, and
hyperparameters (Figure 1), we can easily verify that the executable
corresponds to the ported model. Note that the porting exercise might perform
some weight quantization. This is acceptable since the conversion is
deterministic, and the final certification metrics and performance are
computed on the final inference model.

####  IV-B5 Requirements verification

Using our inference model, we run the probabilistic safety validation
algorithm and obtain a failure boundary similar to that of Figure 15. For the
sake of the example, we stop the BSV algorithm early at only 33 iterations,
yielding an estimated failure probability of the ML subsystem is
$P_{\text{fail}}=5.8\times 10^{-3}$. By placing ourselves in a situation where
performance does not meet the initial $10^{-4}$ requirement, we can now
illustrate example follow-up steps for applicants:

  * •

Refining the ODD: we can decide to revise the subsystem requirements or
concept of operations (CONOPS) to restrict the model to regions where it
performs well. In our example, the applicant could decide to only activate the
runway detector below 2 nautical miles from the runway to now meet the
$10^{-4}$ $P_{\text{fail}}$ requirement.

  * •

Improving the model: BSV highlights failure regions which can be used to
update the model to perform better in these areas. To prevent simply
overfitting the certification dataset, the data access policy should prevent
the counter-examples themselves from being used directly to re-train the
model. Rather, general characteristics of those failure modes and regions
could be communicated to developers.

  * •

Validating the surrogate: the BSV generates a minimal set of test points which
can be collected with real-world flight test data to validate the simulator,
and the BSV failure boundary.

This example highlights the iterative nature of the machine learning model
development. We repeat all the steps of the process, keeping track of the
different experiments and artifacts, until a satisfactory performance is met.

Figure 15: Probabilistic failure surrogate from an example iteration of the
BSV algorithm on the runway detector. The algorithm efficiently samples (the
squares) the ODD (distance, glideslope) to find the failure boundary (white).
The runway detector failures are represented in red, successes in green.

##  V Conclusion

Safety is an emergent system property, not a component property: we considered
a machine learning model as a small component of a more comprehensive
autonomous flight stack. This hierarchical decomposition bounds the output of
the model with many redundancies and safety layers. The stack follows the
conventional certification process for traditional software subsystems. For
the machine learning subsystem, we propose a model-agnostic and tool-
independent framework applicable to a wide variety of use cases, data types,
and machine learning architectures.

Throughout a modified V-diagram, we emphasize the traceability of data,
metadata, and model artifacts. After defining a complete operational domain,
we sample simulated test cases to run through the model. We present the first
use of the Bayesian safety validation (BSV) algorithm [39] in a complete
certification argument. This statistical verifier scales to the large
architectures currently used in research and industry. As opposed to other
verification approaches, we sample the operational domain (e.g., distance to
runway, weather), not the sensor space (e.g., pixel color, brightness). The
verifier converges on an accurate estimate of the probability of failure
$P_{\text{fail}}$ of the ML component, which we can tie back to requirements,
closing the ML certification gap. This work is currently being used to
supplement the FAA certification process of an autonomous cargo aircraft.

## References

  * [1] J. Dean, “A Golden Decade of Deep Learning: Computing Systems & Applications,” _Daedalus_ , vol. 151, no. 2, pp. 58–74, 05 2022.
  * [2] J.-G. Durand, “Journey to Autonomy: The Potential of Vision-Based Systems,” Oct. 2022. [Online]. Available: https://www.linkedin.com/pulse/journey-autonomy-potential-vision-based-systems-xwing
  * [3] S-18 Aircraft and Sys Dev and Safety Assessment Committee, “Guidelines for Development of Civil Aircraft and Systems,” SAE International, Tech. Rep. ARP4754A / ED-79A, Dec. 2010.
  * [4] SAE International, “Guidelines and Methods for Conducting the Safety Assessment Process on Civil Airborne Systems and Equipment,” SAE International, Tech. Rep. ARP4761, Jan. 1996.
  * [5] RTCA SC-180 / EUROCAE WG-46, “Design Assurance Guidance for Airborne Electronic Hardware,” RTCA, Inc. / EUROCAE, Tech. Rep. DO-254 / ED-80, Apr. 2000\.
  * [6] RTCA SC-205 / EUROCAE WG-12, “Software Considerations in Airborne Systems and Equipment Certification,” RTCA, Inc. / EUROCAE, Tech. Rep. DO-178C / ED-12C, Dec. 2011.
  * [7] RTCA, “Standards for Processing Aeronautical Data,” RTCA, Inc., Tech. Rep. DO-200B, Jun. 2015.
  * [8] FAA, “Verification of Adaptive Systems,” Federal Aviation Administration, Tech. Rep. Final report, Apr. 2016.
  * [9] EASA, “Artificial Intelligence Roadmap: A human-centric approach to AI in aviation,” European Union Aviation Safety Agency, Tech. Rep. 1.0, Feb. 2020\.
  * [10] ——, “Artificial Intelligence Roadmap: A human-centric approach to AI in aviation,” European Union Aviation Safety Agency, Tech. Rep. 2.0, May 2023.
  * [11] Innovation Network and EASA AI Task Force and Daedalean AG, “Concepts of Design Assurance for Neural Networks (CoDANN),” European Union Aviation Safety Agency, Tech. Rep. 1.0, Mar. 2020.
  * [12] ——, “Concepts of Design Assurance for Neural Networks (CoDANN) II,” European Union Aviation Safety Agency, Tech. Rep. 1.0, May 2021.
  * [13] G. Balduzzi, M. Ferrari Bravo, A. Chernova, C. Cruceru, L. van Dijk, P. de Lange, J. Jerez, N. Koehler, M. Koerner, C. Perret-Gentil, Z. Pillio, R. Polak, H. Silva, R. Valentin, I. Whittington, and G. Yakushev, “Neural network based runway landing guidance for general aviation autoland,” Federal Aviation Administration, Tech. Rep., Nov 2021, DOT/FAA/TC-21/48.
  * [14] EASA, “EASA Concept Paper: First usable guidance for Level 1 machine learning applications,” European Union Aviation Safety Agency, Tech. Rep. 01, Apr. 2021\.
  * [15] ——, “EASA Concept Paper: First usable guidance for Level 1&2 machine learning applications,” European Union Aviation Safety Agency, Tech. Rep. 02, Feb. 2023.
  * [16] M. Gariel, B. Shimanuki, R. Timpe, and E. Wilson, “Framework for Certification of AI-Based Systems,” Mar. 2021, xwing Inc. internal document.
  * [17] C. Liu, T. Arnon, C. Lazarus, C. Strong, C. Barrett, M. J. Kochenderfer _et al._ , “Algorithms for verifying deep neural networks,” _Foundations and Trends® in Optimization_ , vol. 4, no. 3-4, pp. 244–404, 2021.
  * [18] E. Botoeva, P. Kouvaros, J. Kronqvist, A. Lomuscio, and R. Misener, “Efficient verification of relu-based neural networks via dependency analysis,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, 2020, pp. 3291–3299.
  * [19] P. Kouvaros, T. Kyono, F. Leofante, A. Lomuscio, D. Margineantu, D. Osipychev, and Y. Zheng, “Formal analysis of neural network-based systems in the aircraft domain,” in _International Symposium on Formal Methods_. Springer, 2021, pp. 730–740.
  * [20] T. Gehr, M. Mirman, D. Drachsler-Cohen, P. Tsankov, S. Chaudhuri, and M. Vechev, “AI2: Safety and robustness certification of neural networks with abstract interpretation,” in _2018 IEEE symposium on security and privacy (SP)_. IEEE, 2018, pp. 3–18.
  * [21] P. Henriksen and A. Lomuscio, “Efficient neural network verification via adaptive refinement and adversarial search,” in _ECAI 2020_. IOS Press, 2020, pp. 2513–2520.
  * [22] P. Henriksen, K. Hammernik, D. Rueckert, and A. Lomuscio, “Bias field robustness verification of large neural image classifiers,” in _British Machine Vision Conference (BMVC21)_ , 2021.
  * [23] G. Katz, C. Barrett, D. L. Dill, K. Julian, and M. J. Kochenderfer, “Reluplex: An efficient SMT solver for verifying deep neural networks,” in _International conference on computer aided verification_. Springer, 2017, pp. 97–117.
  * [24] G. Katz, D. A. Huang, D. Ibeling, K. Julian, C. Lazarus, R. Lim, P. Shah, S. Thakoor, H. Wu, A. Zeljić _et al._ , “The Marabou Framework for Verification and Analysis of Deep Neural Networks,” in _International Conference on Computer Aided Verification_. Springer, 2019, pp. 443–452.
  * [25] Y. Tian, K. Pei, S. Jana, and B. Ray, “Deeptest: Automated testing of deep-neural-network-driven autonomous cars,” in _International Conference on Software Engineering_ , 2018, pp. 303–314.
  * [26] G. Singh, T. Gehr, M. Mirman, M. Püschel, and M. Vechev, “Fast and effective robustness certification,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [27] D. Gopinath, G. Katz, C. S. Păsăreanu, and C. Barrett, “Deepsafe: A data-driven approach for assessing robustness of neural networks,” in _International symposium on automated technology for verification and analysis_. Springer, 2018, pp. 3–19.
  * [28] X. Xie, L. Ma, F. Juefei-Xu, M. Xue, H. Chen, Y. Liu, J. Zhao, B. Li, J. Yin, and S. See, “Deephunter: a coverage-guided fuzz testing framework for deep neural networks,” in _Proceedings of the 28th ACM SIGSOFT International Symposium on Software Testing and Analysis_ , 2019, pp. 146–157.
  * [29] J. Mohapatra, T.-W. Weng, P.-Y. Chen, S. Liu, and L. Daniel, “Towards verifying robustness of neural networks against a family of semantic perturbations,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 244–252.
  * [30] Y. Yang and M. Rinard, “Correctness verification of neural networks,” _arXiv preprint arXiv:1906.01030_ , 2019.
  * [31] P. Linardatos, V. Papastefanopoulos, and S. Kotsiantis, “Explainable AI: A review of machine learning interpretability methods,” _Entropy_ , vol. 23, no. 1, p. 18, 2020.
  * [32] M. T. Ribeiro, S. Singh, and C. Guestrin, “”Why should i trust you?” Explaining the predictions of any classifier,” in _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_ , 2016, pp. 1135–1144.
  * [33] A. Corso, R. J. Moss, M. Koren, R. Lee, and M. J. Kochenderfer, “A Survey of Algorithms for Black-Box Safety Validation of Cyber-Physical Systems,” _Journal of Artificial Intelligence Research_ , 2021.
  * [34] R. Lee, O. Mengshoel, A. Saksena, R. Gardner, D. Genin, J. Silbermann, M. Owen, and M. Kochenderfer, “Adaptive Stress Testing: Finding Likely Failure Events with Reinforcement Learning,” _Journal of Artificial Intelligence Research_ , vol. 69, pp. 1165–1201, 12 2020.
  * [35] J. Norden, M. O’Kelly, and A. Sinha, “Efficient black-box assessment of autonomous vehicle safety,” _Computing Research Repository (CoRR) in arXiv_ , vol. abs/1912.03618, 2019.
  * [36] T. T. Pham, “Identifying the Characteristics of Training Data as Necessary Conditions for Machine Learning Verification,” Oct. 2022, Federal Aviation Administration internal document.
  * [37] N. Kalra and S. M. Paddock, “Driving to safety: How many miles of driving would it take to demonstrate autonomous vehicle reliability?” _Transportation Research Part A: Policy and Practice_ , vol. 94, pp. 182–193, 2016.
  * [38] Y. He and J. Schumann, “A framework for the analysis of deep neural networks in aerospace applications using Bayesian statistics,” in _2020 International Joint Conference on Neural Networks (IJCNN)_. IEEE, 2020, pp. 1–9.
  * [39] R. J. Moss, M. J. Kochenderfer, M. Gariel, and A. Dubois, “Bayesian Safety Validation for Black-Box Systems,” in _AIAA Aviation_ , 2023.
  * [40] G. Kimchi, “Moving Forward Safely: Amazon’s Approach to Drone Delivery,” https://www.youtube.com/watch?v=yHY-ZWpC8T4, Jun. 2019.
  * [41] K. He, G. Gkioxari, P. Dollár, and R. Girshick, “Mask R-CNN,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 2961–2969.
  * [42] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.
  * [43] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft COCO: Common Objects in Context,” in _European conference on computer vision_. Springer, 2014, pp. 740–755.