# Galaxy cluster hydrostatic bias in Kottler spacetime

Shantanu Desai<EMAIL_ADDRESS>Department of Physics, Indian Institute of
Technology, Hyderabad, Telangana-502285, India

(September 2021)

###### Abstract

We calculate the relativistic corrections to hydrostatic X-ray masses for
galaxy clusters in Kottler spacetime, which is the spherically symmetric
solution to Einstein’s equations in General relativity endowed with a
cosmological constant. The hydrostatic masses for clusters (calculated
assuming Newtonian gravity) have been found to be underestimated compared to
lensing masses, and this discrepancy is known as hydrostatic mass bias. Since
the relativistic hydrostatic X-ray masses are automatically lower than lensing
masses, under the edifice of Kottler metric, we check if the hydrostatic mass
bias problem gets alleviated using this ansatz. We consider a sample of 18
galaxy clusters for this pilot test. We find that the ratio of X-ray to
lensing mass is close to unity even in Kottler spacetime. Therefore, the
effect of relativistic corrections to hydrostatic X-ray masses for galaxy
clusters is negligible.

##  I Introduction

Galaxy clusters are the most massive gravitationally collapsed objects in the
universe Vikhlinin et al. (2014); Kravtsov and Borgani (2012); Allen et al.
(2011). Galaxy clusters have proved to be wonderful laboratories for
cosmology, galaxy evolution, modified gravity theories, and fundamental
Physics Allen et al. (2011); Kravtsov and Borgani (2012); Vikhlinin et al.
(2014); Desai (2018); Bora and Desai (2021); Mendonça et al. (2021). In the
past three decades a large number of galaxy clusters have been discovered
through a whole suite of optical, infrared, X-Ray and mm-wave surveys.
Historically, galaxy clusters have also been one of the key probes, which
helped corroborate current concordance $\Lambda$CDM model of Cosmology
consisting of about 70% Dark energy and 25% Dark matter Vikhlinin et al.
(2009); Rozo et al. (2010); Mantz et al. (2014); Bocquet et al. (2015).

However, with the availability of more precise data, it was found that some of
the cosmological parameters obtained using cluster counts (eg. $\sigma_{8}$)
do not agree with those from primary CMB Ade et al. (2014). Secondly, the
hydrostatic masses of galaxy clusters estimated using X-ray measurements were
found to be underestimated with respect to the Weak lensing masses Ettori et
al. (2013); Sereno and Ettori (2015); Lovisari et al. (2020). This discrepancy
has been characterized by positing a hydrostatic mass bias parameter, which
quantifies the difference between Weak lensing and hydrostatic masses. Such a
mass bias can also resolve the tensions in Cosmology between cluster counts
and primary CMB. The main cause of the hydrostatic mass bias is due to the
non-thermal pressure support in clusters Shi et al. (2016); Ota et al. (2018);
Gianfagna et al. (2021). An uptodate compilation of this bias parameter from
both observational data and simulations can be found in Gianfagna et al.
(2021).

X-ray masses have also been used in conjunction with lensing masses to test
modified theories of gravity, which dispense with dark matter and dark energy
Tian et al. (2020); Pradyumna et al. (2021); Gopika and Desai (2020);
Nieuwenhuizen (2017), and also to test non-standard dark matter scenarios such
as fermionic dark matter. Some examples of these studies tests using X-ray and
lensing data for A1689 and A1835 can be found in Nieuwenhuizen et al. (2021);
Nieuwenhuizen (2020).

However, all these comparisons of X-Ray masses with weak lensing masses have
done by calculating the X-ray hydrostatic masses using Newtonian gravity.
However, the mass of an object also depends on the theory of gravity and the
background space-time assumed. Since lensing is a purely General Relativity
(GR) based effect and does not occur in Newtonian gravity, it is important to
also calculate the X-ray masses using GR, which is the main goal of this work.
For this purpose, we consider the Kottler metric for the background spacetime,
and evaluate its impact in ameliorating some of the discrepancies in the mass
estimates of galaxy clusters. The Kottler metric (also known as Schwarzschild-
Desitter metric) is the unique spherically symmetric solution of Einstein’s
General Relativity in the presence of the Cosmological constant $\Lambda$. We
know for more than two decades that we live in a dark energy dominated
universe, whose equation of state is close to the Cosmological constant
Huterer and Shafer (2018). Therefore, it is natural to consider such a metric
while calculating relativistic hydrostatic masses within GR. The line element
for the Kottler metric can be written as Kottler (1918):

$ds^{2}=A(r)dt^{2}-\frac{dr^{2}}{A(r)}-r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})$ (1)

where $A(r)=1-\frac{2GM}{r}-\frac{\Lambda}{3}r^{2}$, and $\Lambda$ is the
Cosmological constant.

In this work, we carry out a pilot study to study how the X-ray masses for
galaxy clusters and galaxy groups change under the aegis of Kottler metric,
and whether it alleviates the galaxy cluster mass bias problem. Previously,
relativistic corrections to Newtonian hydrostatic masses have been computed
for galaxy clusters using Tolman-Oppenheimer-Volkoff equation and shown to be
negligible Gupta and Desai (2020).

The outline of this manuscript is as follows. In Sec. II, we recap some
previous works in literature related to galaxy clusters in Kottler spacetime.
The impact of Kottler space time on the hydrostatic masses for our sample of
clusters is considered in Sec. III. We conclude in Sec. IV.

##  II Cluster dynamics in Kottler Spacetime

There are multiple methods to estimate galaxy cluster mass Ettori et al.
(2013); Lovisari et al. (2020). The classic methods involve galaxy kinematics
using line of sight galaxy velocity dispersions and caustic methods, X-ray and
Sunyaev-Zeldovich (SZ) observations involving the assumption that the gas is
in hydrostatic equilibrium.

The first seminal study of the impact of Kottler metric on galaxy cluster mass
was carried out by Bambi Bambi (2007) who pointed out the effective Newtonian
mass of a galaxy cluster ($M_{eff}(R)$) at a distance ($R$) estimated using
X-ray/SZ or velocity dispersion measurements in the weak field limit is given
by Bambi (2007):

$M_{eff}(R)=M(R)-\frac{8}{3}\pi r^{3}\rho_{\Lambda},$ (2)

where $M(R)$ is the true mass at radius $R$ and $\rho_{\Lambda}$ is the energy
density of the cosmological constant. However, the hydrostatic mass that is
usually determined. (in literature) from temperature, density or pressure
profiles is assumed to be the same as $M(R)$ Sarazin (1986). However, as
emphasized by Bambi, what these methods really measure is $M_{eff}(R)$ and not
$M(R)$.

The other widely used technique to measure cluster mass is gravitational
lensing Hoekstra et al. (2013). There has been a long debate in literature on
whether the bending of light is affected by the cosmological constant or not.
One school of thought has argued that since $\Lambda$ does not appear in the
null geodesic equation, it does not affect the lensing results Kagramanova et
al. (2006); Khriplovich and Pomeransky (2008); Simpson et al. (2010); Bambi
(2007). However, this viewpoint has been disputed by other groups, who argue
that $\Lambda$ affects the bending of light, because of the background space-
time, which has $\Lambda$ in-built into it Rindler and Ishak (2007); Sereno
(2008); Ishak et al. (2010); Ishak and Rindler (2010). Upper limits on the
cosmological constant were also set, using the Einstein radii measured for
galaxies and clusters Ishak et al. (2008). The most recent exposition on this
issue can be found in Hu et al. (2021), who found no dependence of light
bending on $\Lambda$ using numerical integration of the geodesic equation of
motion for a Swiss cheese model made up of a point mass and a vacuole. Bambi
(2007) has pointed out that the lensing mass does not get affected by
$\Lambda$. Therefore, assuming $\rho_{\Lambda}\approx 6\times 10^{-30}g/cc$,
corresponding to 70% composition by dark energy, one gets the following
expression for the ratio of $M_{eff}$ (which is what X-ray hydrostatic masses
correspond to) and $M$, which one gets using lensing Bambi (2007)

$\frac{M_{eff}}{M}=1-0.007\left(\frac{10^{14}M_{\odot}}{M}\right)[r(Mpc)]^{3}$
(3)

We evaluate the impact of this equation on a sample of few galaxy clusters in
the next section.

We also briefly discuss some other works in the literature related to Kottler
metric and clusters. Bisnovatyi-Kogan and Chernin (2012) studied the effect of
Kottler metric on VIRGO-like galaxy clusters and showed that that the radial
extent as well as the average density of dark matter haloes is determined
$\Lambda$. In a followup paper, Chernin et al. (2013) then studied the impact
of Kottler metric on the mass of Coma cluster and showed that the effective
gravitating mass of the Coma cluster could be three times smaller than the
regular Newtonian mass at 14 Mpc. Teerikorpi et al. (2015) constructed a
$\Lambda$-significance graph which characterizes the parameter space in Mass-
radius plane for bound objects, where $\Lambda$ dominates the dynamics. They
also calculated the zero-gravity radius where the repulsive force due to the
$\Lambda$ term is equal to the attractive Newtonian gravitational force for
various clusters and superclusters. Other results in literature related to
clusters include studies of dynamical stability Gurzadyan et al. (2020) of
Kottler spacetime, and estimation of $\Lambda$ using velocity dispersions in
groups and clusters Gurzadyan and Stepanian (2019).

##  III Hydrostatic masses in Kottler spacetime

For our pilot study, we now test the impact of Kottler metric on the
hydrostatic masses for galaxy clusters. For this purpose, we calculate the
hydrostatic masses for 18 clusters compiled in Mahdavi et al. (2008), where
Newtonian hydrostatic and Weak lensing masses were calculated. The Weak
lensing and X-ray data was obtained with observations obtained from CFHT
telescope and Chandra X-ray observatory, respectively. This aforementioned
work had found a decreasing trend of the ratio of X-ray ($M_{X}$) to WL mass
($M_{L}$) with increasing radii, with $\frac{M_{X}}{M_{L}}=0.78\pm 0.09$ at
$R_{500}$ Mahdavi et al. (2008). We now use Eq. 3 to compute this ratio for
all the 18 clusters in this sample. Our results are tabulated in Table 1. We
find that this ratio is close to 1 (99.8%) for all clusters. Therefore, the
effective mass in Kottler spacetime is almost the same as the Newtonian mass.
We also calculated the same ratio for a sample of galaxy groups from Sun et
al. (2009), whose masses are about 10 times smaller than clusters. Even for
this group sample, the ratio of $\frac{M_{X}}{M_{L}}$ is about 99.7%.
Therefore, we conclude that even though the hydrostatic masses are less than
Weak lensing masses in Kottler space-time, their ratio is close to one, and
this does not alleviate the hydrostatic mass bias problem found.

Cluster | $r_{500}$ | $M_{Newt}$ | $M_{eff}/M]_{r=r_{500}}$
---|---|---|---
| (kpc) | $(M_{\bigodot})$ | (%)
A68 | 1.22 | 6.64 | 99.8
A209 | 1.27 | 7.14 | 99.8
A267 | 1.21 | 6.29 | 99.8
A370 | 1.47 | 13.27 | 99.8
A383 | 1.09 | 4.49 | 99.8
A963 | 1.06 | 4.16 | 99.8
A1689 | 1.61 | 14.29 | 99.8
A1763 | 1.43 | 10.47 | 99.8
A2218 | 1.22 | 6.1 | 99.8
A2219 | 1.42 | 10.27 | 99.8
A2390 | 1.35 | 8.79 | 99.8
CL0024.0+1652 | 1.32 | 9.87 | 99.8
MS0015.9+1609 | 1.56 | 19.51 | 99.8
MS0906.5+1110 | 1.41 | 9.46 | 99.8
MS1358.1+6245 | 1.18 | 6.64 | 99.8
MS1455.0+2232 | 1.09 | 4.83 | 99.8
MS1512.4+3647 | 0.89 | 2.94 | 99.8
MS 1621.5+2640 | 1.19 | 7.64 | 99.8
Table 1: Newtonian hydrostatic mass (third column) as well as the ratio of
effective mass to the total mass in Kottler spacetime (fourth column)
evaluated at $r_{500}$, for 18 galaxy clusters tabulated in Mahdavi et al.
(2008) using Eq. 3. In all cases, the difference between effective and true
mass is less than 0.2%

##  IV Conclusions

A large number of studies have found that the X-ray masses of galaxy clusters
are underestimated with respect to lensing masses. This problem is known as
the hydrostatic mass bias problem and is attributed to astrophysical
mechanisms such as a non-thermal pressure support.

Here, we point out all estimates of hydrostatic masses of galaxy clusters in
literature have been calculated using Newtonian gravity. Since we know for
more than two decades that we live in an accelerating universe, which could be
driven by a cosmological constant $\Lambda$, we look at the effect on galaxy
masses in GR using the Kottler metric for the background spacetime, which is
the spherically symmetric solution to Einstein’s field equations endowed with
a cosmological constant. Although, there have been a few previous works
studying the dynamics of clusters in Kottler space-time (cf. Sect. II), no one
has looked the effect of hydrostatic mass bias using the Kottler metric for
observed cluster samples.

It has been pointed out more than a decade ago by Bambi Bambi (2007) that in
Kottler spacetime, hydrostatic and velocity dispersion based masses are
smaller than the lensing mass. Here, we try to assess the change in the
hydrostatic masses using this ansatz for an observed galaxy cluster sample, to
see if the hydrostatic mass bias problem would get alleviated in Kottler
spacetime. For this pilot study, we considered a sample of 18 galaxy clusters
with both X-ray and Weak lensing masses for which $\frac{M_{X}}{M_{L}}$ was
found to be about 70% Mahdavi et al. (2008).

We then calculated the ratio of X-ray to Weak lensing masses in Kottler space
time using Eq. 3. Our results can be found in Table 1. We find that the
$\frac{M_{X}}{M_{L}}$ is equal to 99.8% for all the clusters in our sample.
Therefore, the impact of general relativistic corrections to hydrostatic
masses of galaxy clusters using the Kottler metric is negligible.

###### Acknowledgements.

We are grateful to Cosimo Bambi for useful discussions on this subject.

## References

  * Vikhlinin et al. (2014) A. A. Vikhlinin, A. V. Kravtsov, M. L. Markevich, R. A. Sunyaev, and E. M. Churazov, Physics Uspekhi 57, 317-341 (2014).
  * Kravtsov and Borgani (2012) A. V. Kravtsov and S. Borgani, Ann. Rev. Astron. Astrophys. 50, 353 (2012), eprint 1205.5556.
  * Allen et al. (2011) S. W. Allen, A. E. Evrard, and A. B. Mantz, Ann. Rev. Astron. Astrophys. 49, 409 (2011), eprint 1103.4829.
  * Desai (2018) S. Desai, Physics Letters B 778, 325 (2018), eprint 1708.06502.
  * Bora and Desai (2021) K. Bora and S. Desai, JCAP 2021, 012 (2021), eprint 2008.10541.
  * Mendonça et al. (2021) I. E. C. R. Mendonça, K. Bora, R. F. L. Holanda, S. Desai, and S. H. Pereira, JCAP 2021, 034 (2021), eprint 2109.14512.
  * Vikhlinin et al. (2009) A. Vikhlinin, A. V. Kravtsov, R. A. Burenin, H. Ebeling, W. R. Forman, A. Hornstrup, C. Jones, S. S. Murray, D. Nagai, H. Quintana, et al., Astrophys. J. 692, 1060 (2009), eprint 0812.2720.
  * Rozo et al. (2010) E. Rozo, R. H. Wechsler, E. S. Rykoff, J. T. Annis, M. R. Becker, A. E. Evrard, J. A. Frieman, S. M. Hansen, J. Hao, D. E. Johnston, et al., Astrophys. J. 708, 645 (2010), eprint 0902.3702.
  * Mantz et al. (2014) A. B. Mantz, S. W. Allen, R. G. Morris, D. A. Rapetti, D. E. Applegate, P. L. Kelly, A. von der Linden, and R. W. Schmidt, Mon. Not. R. Astron. Soc. 440, 2077 (2014), eprint 1402.6212.
  * Bocquet et al. (2015) S. Bocquet, A. Saro, J. J. Mohr, K. A. Aird, M. L. N. Ashby, M. Bautz, M. Bayliss, G. Bazin, B. A. Benson, L. E. Bleem, et al., Astrophys. J. 799, 214 (2015), eprint 1407.2942.
  * Ade et al. (2014) P. A. R. Ade et al. (Planck), Astron. Astrophys. 571, A20 (2014), eprint 1303.5080.
  * Ettori et al. (2013) S. Ettori, A. Donnarumma, E. Pointecouteau, T. H. Reiprich, S. Giodini, L. Lovisari, and R. W. Schmidt, Space Science Reviews 177, 119 (2013), eprint 1303.3530.
  * Sereno and Ettori (2015) M. Sereno and S. Ettori, Mon. Not. R. Astron. Soc. 450, 3633 (2015), eprint 1407.7868.
  * Lovisari et al. (2020) L. Lovisari, S. Ettori, M. Sereno, G. Schellenberger, W. R. Forman, F. Andrade-Santos, and C. Jones, Astron. & Astrophys. 644, A78 (2020), eprint 2010.03582.
  * Shi et al. (2016) X. Shi, E. Komatsu, D. Nagai, and E. T. Lau, Mon. Not. R. Astron. Soc. 455, 2936 (2016), eprint 1507.04338.
  * Ota et al. (2018) N. Ota, D. Nagai, and E. T. Lau, PASJ 70, 51 (2018), eprint 1507.02730.
  * Gianfagna et al. (2021) G. Gianfagna, M. De Petris, G. Yepes, F. De Luca, F. Sembolini, W. Cui, V. Biffi, F. Kéruzoré, J. Macías-Pérez, F. Mayet, et al., Mon. Not. R. Astron. Soc. 502, 5115 (2021), eprint 2010.03634.
  * Tian et al. (2020) Y. Tian, K. Umetsu, C.-M. Ko, M. Donahue, and I. N. Chiu, Astrophys. J. 896, 70 (2020), eprint 2001.08340.
  * Pradyumna et al. (2021) S. Pradyumna, S. Gupta, S. Seeram, and S. Desai, Physics of the Dark Universe 31, 100765 (2021), eprint 2011.06421.
  * Gopika and Desai (2020) K. Gopika and S. Desai, Physics of the Dark Universe 30, 100707 (2020), eprint 2006.12320.
  * Nieuwenhuizen (2017) T. M. Nieuwenhuizen, Fortschritte der Physik 65, 201600050 (2017), eprint 1610.01543.
  * Nieuwenhuizen et al. (2021) T. M. Nieuwenhuizen, M. Limousin, and A. Morandi, European Physical Journal Special Topics 230, 1137 (2021).
  * Nieuwenhuizen (2020) T. M. Nieuwenhuizen, Fluctuation and Noise Letters 19, 2050016-301 (2020).
  * Huterer and Shafer (2018) D. Huterer and D. L. Shafer, Reports on Progress in Physics 81, 016901 (2018), eprint 1709.01091.
  * Kottler (1918) F. Kottler, Annalen der Physik 361, 401 (1918).
  * Gupta and Desai (2020) S. Gupta and S. Desai, Physics of the Dark Universe 28, 100499 (2020), eprint 1909.07408.
  * Bambi (2007) C. Bambi, Phys. Rev. D 75, 083003 (2007), eprint astro-ph/0703645.
  * Sarazin (1986) C. L. Sarazin, Reviews of Modern Physics 58, 1 (1986).
  * Hoekstra et al. (2013) H. Hoekstra, M. Bartelmann, H. Dahle, H. Israel, M. Limousin, and M. Meneghetti, Space Science Reviews 177, 75 (2013), eprint 1303.3274.
  * Kagramanova et al. (2006) V. Kagramanova, J. Kunz, and C. Lämmerzahl, Physics Letters B 634, 465 (2006), eprint gr-qc/0602002.
  * Khriplovich and Pomeransky (2008) I. B. Khriplovich and A. A. Pomeransky, International Journal of Modern Physics D 17, 2255 (2008), eprint 0801.1764.
  * Simpson et al. (2010) F. Simpson, J. A. Peacock, and A. F. Heavens, Mon. Not. R. Astron. Soc. 402, 2009 (2010), eprint 0809.1819.
  * Rindler and Ishak (2007) W. Rindler and M. Ishak, Phys. Rev. D 76, 043006 (2007), eprint 0709.2948.
  * Sereno (2008) M. Sereno, Phys. Rev. D 77, 043004 (2008), eprint 0711.1802.
  * Ishak et al. (2010) M. Ishak, W. Rindler, and J. Dossett, Mon. Not. R. Astron. Soc. 403, 2152 (2010), eprint 0810.4956.
  * Ishak and Rindler (2010) M. Ishak and W. Rindler, General Relativity and Gravitation 42, 2247 (2010), eprint 1006.0014.
  * Ishak et al. (2008) M. Ishak, W. Rindler, J. Dossett, J. Moldenhauer, and C. Allison, Mon. Not. R. Astron. Soc. 388, 1279 (2008), eprint 0710.4726.
  * Hu et al. (2021) L. Hu, A. Heavens, and D. Bacon, arXiv e-prints arXiv:2109.09785 (2021), eprint 2109.09785.
  * Bisnovatyi-Kogan and Chernin (2012) G. S. Bisnovatyi-Kogan and A. D. Chernin, Astrophys. and Space Science 338, 337 (2012), eprint 1206.1433.
  * Chernin et al. (2013) A. D. Chernin, G. S. Bisnovatyi-Kogan, P. Teerikorpi, M. J. Valtonen, G. G. Byrd, and M. Merafina, Astron. & Astrophys. 553, A101 (2013), eprint 1303.3800.
  * Teerikorpi et al. (2015) P. Teerikorpi, P. Heinämäki, P. Nurmi, A. D. Chernin, M. Einasto, M. Valtonen, and G. Byrd, Astron. & Astrophys. 577, A144 (2015), eprint 1503.02805.
  * Gurzadyan et al. (2020) V. G. Gurzadyan, A. A. Kocharyan, and A. Stepanian, European Physical Journal C 80, 24 (2020), eprint 2001.02634.
  * Gurzadyan and Stepanian (2019) V. G. Gurzadyan and A. Stepanian, European Physical Journal C 79, 169 (2019), eprint 1902.07171.
  * Mahdavi et al. (2008) A. Mahdavi, H. Hoekstra, A. Babul, and J. P. Henry, Mon. Not. R. Astron. Soc. 384, 1567 (2008), eprint 0710.4132.
  * Sun et al. (2009) M. Sun, G. M. Voit, M. Donahue, C. Jones, W. Forman, and A. Vikhlinin, Astrophys. J. 693, 1142 (2009), eprint 0805.2320.


# Ionization efficiency for nuclear recoils in silicon from $\sim 50$ eV to
$3$ MeV

Y. Sarkis, A. Aguilar-Arevalo and J.C. D’Olivo  Instituto de Ciencias
Nucleares, Universidad Nacional Autónoma de México, 04510 CDMX, Mexico

(August 28, 2024)

###### Abstract

We present a model for the nuclear recoil ionization efficiency in silicon
based on an extension of Lindhard’s theory where atomic bond disruption is
modeled as a function of the initial ion energy, the interatomic potential,
and the average ion-vacancy production energy. A better description of the
electronic stopping than the one assumed by Lindhard, the effect of electronic
straggling, as well as charge screening and Coulomb repulsion effects of ions
are also considered. The model describes the available data over nearly four
orders of magnitude in nuclear recoil energy.

quenching factor, dark matter, nuclear recoils

Introduction. Silicon is a common detection medium in experiments searching
for the low-energy depositions expected from coherent elastic neutrino-nucleus
scattering (CE$\nu$NS) or dark matter (DM) interactions. In these
applications, Si detectors typically observe the ionization signal produced by
a nuclear recoil after interaction with the incident particle. The ionization
efficiency, or quenching factor (QF), defined as the ratio of the amount of
electronic excitation produced by a nuclear recoil to that produced by a
recoiling electron of the same energy, is a key feature of these detectors. As
early as 1963, Lindhard [1] crafted the basic integral equation describing the
partitioning of the energy deposited by a nuclear recoil between electronic
and atomic motion. Since then, his approximate solution for the QF has been
known to describe the data [2, 3, 4] at high energies ($>10$ keV), where
atomic binding energy effects can be safely ignored. It has been shown [5]
that when properly introduced into the model, a constant average binding
energy of the order of 150 eV, consistent with that associated with $2s2p$
shell electrons in Si, can describe the data available at the lowest energies
[6, 7]. Although a constant binding energy model works well for the existing
data, it predicts a cut-off at $\sim$300 eV, significantly greater than the
$36$ eV of energy required in average to create a stable ion-vacancy pair
(Frenkel [8] pair), which can be argued to represent an upper bound on the
minimum energy required to be deposited by a nuclear recoil to produce an
observable ionization signal in a Si detector.

In order to obtain a model that is valid at lower nuclear recoil energies, in
this work we implement several improvements to the integro-differential
equation with binding energy presented in [5]. We assume that ionization can
only occur when an ion can move freely through the lattice but, instead of a
constant threshold, the energy lost due atomic bond breaking is described as a
“binding energy” function dependent on the energy. In the low energy regime it
corresponds to the Frenkel pair creation energy, while at higher energies it
describes the inner atomic excitation energies. We also make use of improved
modeling of the electronic stopping power at low energies, and add a term
related to electronic straggling [9, 10] to the equation.

Consideration of the aforementioned effects results in a first-principles
model of the QF capable of describing the available published data in Si from
0.68 keV to 3 MeV nuclear recoil energy, and provides a prediction for this
quantity down to a few tens of eV.

The nuclear recoil QF and Lindhard’s basic integral equation. When a neutrino
or a DM particle elastically scatters off a nucleus in a Si detector, the
struck ion recoils with an energy $E_{R}$. If the timescale of the collision
is much smaller than that of atomic processes, the ion will loose some energy
to the atomic degrees of freedom emerging from the collision with kinetic
energy $E=E_{R}-U$, where $U$ is the energy used to disrupt the atomic
binding. In general, $U$ is not limited to the energy needed to remove the ion
from its site, but it can also include contributions to excitation or
ionization of bound atomic electrons, and therefore effectively incorporates
the Migdal effect [11, 12].

Lindhard’s theory [1] assumes that $E_{R}$, the total energy deposited by the
incoming particle in the material, is separated between the energy given to
electrons, $H$ and the energy given to atomic motion, $N$ so that $E_{R}=H+N$.
For binary collisions, the Ziegler scaling length [13] $a=(0.8853)a_{0}/Z$,
with $Z=Z_{1}^{0.23}+Z_{2}^{0.23}$, and $a_{0}$ the Bohr radius, can be used
with suitable values of the effective atomic numbers of the incident ion and
the target atom $Z_{1}$ and $Z_{2}$, respectively, to define the reduced
dimensionless quantities $\varepsilon=CE$, $\eta=CH$, and $\nu=CN$ where
$C=16.2616/ZZ_{1}Z_{2}$ keV-1. Assuming that this separation holds on average
over the large number of collisions occurring in the process, the average
quantities also satisfy $\varepsilon_{R}=\bar{\eta}+\bar{\nu}$. The nuclear
recoil QF is given by
$f_{n}=\bar{\eta}/\varepsilon_{R}=(\varepsilon+u-\bar{\nu})/(\varepsilon+u)$,where
$\varepsilon_{R}=\varepsilon+u$, and $u=CU$.

Consider an ion moving in the material with energy $E$ and colliding with an
atom, transferring an energy $T_{n}$ to its center of mass and an energy
$\sum_{i}T_{ei}$ to a certain number of ionized electrons. The ion will
scatter with an energy $E-T_{n}-\sum_{i}T_{ei}$, and the struck atom will
recoil with an energy $T_{n}-U$. Lindhard’s basic integral equation for
$\bar{\nu}(\varepsilon)$,

$\int
d\sigma_{n,e}\left[\bar{\nu}(E-T_{n}-\Sigma_{i}T_{ei})+\bar{\nu}(T_{n}-U)-\bar{\nu}(E)\right]=0,$
(1)

states that the average energy given to atomic motion by the initial ion with
energy $E$ equals the sum of the contributions of the scattered ion and the
struck recoiling ion, where the contribution due to ejected electrons has been
neglected (approximation I). Integration over the total nuclear and electronic
cross sections $\int d\sigma_{n,e}$ represents the sum over all possible
impact parameters for nuclear and electronic collisions.

Lindhard made use of four additional approximations to construct an
approximate integro differential equation for $\bar{\nu}(\varepsilon)$: II)
the energy transferred to ionized electrons is small; III) electronic and
atomic collisions can be treated separately; IV) $T_{n}$ is also small
compared to $E$; V) neglect the binding energy $U$. As has been shown in [5],
dropping V results in a higher order approximation to the integro-differential
equation. In what follows, we adopt the same approach.

Improved integro-differential equation Using approximations I-IV, Eq.(1) can
be written in a form suitable for numerical solution, where the electronic
stopping and straggling appear naturally. If $\Delta E$ is the energy lost by
an ionizing projectile, straggling is defined as the mean-square fluctuation
$\Omega^{2}=\langle(\Delta E-\langle\Delta E\rangle)^{2}\rangle$ [10].
Lindhard’s partitioning of the deposited energy into electronic and atomic
motion, guarantees that the electronic and nuclear contributions to straggling
are equal [1]. The electronic straggling per unit length can be expressed as
$N_{e}^{-1}(d\Omega^{2}/dR)_{e}=\int d\sigma_{e}(\Sigma_{i}T_{ei})^{2}$ [10],
where $N_{e}$ is the electron number density and $R$ is the distance traveled
by the projectile. In terms of the reduced quantities $\varepsilon$,
$\omega^{2}=C^{2}\Omega^{2}$ and $\rho=\pi a^{2}N_{e}R$, we have
$W(\varepsilon)=d\omega^{2}/d\rho$. Similarly, for the electronic stopping
power $N_{e}^{-1}(dE/dR)_{e}=\int d\sigma_{e}(\Sigma_{i}T_{ei})$, and
$S_{e}(\varepsilon)=d\varepsilon/d\rho$.

Expanding the first term in Eq.(1) up to second order in
$\Sigma_{i}T_{ei}/(E-T_{n})$ (approximations II-IV), as
$\bar{\nu}\left(E-T_{n}-\Sigma_{i}T_{ei}\right)\approx\bar{\nu}\left(E-T_{n}\right)-\bar{\nu}^{\prime}(E)\left(\Sigma_{i}T_{ei}\right)+\bar{\nu}^{\prime\prime}(E)T_{n}\left(\Sigma_{i}T_{ei}\right)+\frac{1}{2}\bar{\nu}^{\prime\prime}(E){\left(\Sigma_{i}T_{ei}\right)^{2}}$,
and following the arguments in [5], the electronic stopping power and
straggling can be identified in the terms containing $\Sigma_{i}T_{ei}$ to
first and second order, respectively. With these considerations, we write the
modified simplified integral equation with straggling as

$\displaystyle-\frac{1}{2}$ $\displaystyle\varepsilon
S_{e}(\varepsilon)\left(1+\frac{W(\varepsilon)}{S_{e}(\varepsilon)\varepsilon}\right)\bar{\nu}^{\prime\prime}(\varepsilon)+S_{e}(\varepsilon)\bar{\nu}^{\prime}(\varepsilon)=$
(2) $\displaystyle\int_{\varepsilon
u}^{\varepsilon^{2}}dt\frac{f\left(t^{1/2}\right)}{2t^{3/2}}\times[\bar{\nu}(\varepsilon-t/\varepsilon)+\bar{\nu}(t/\varepsilon-u)-\bar{\nu}(\varepsilon)],$

where $u(\varepsilon)$ is the binding energy,
$t=\varepsilon^{2}\sin^{2}(\theta/2)$, with $\theta$ the scattering angle in
the center of mass frame, and $f(t^{1/2})$ is related to the nuclear
differential cross section in the Lindhard-Scharff-Schiott (LSS)
approximation, $d\sigma_{n}=\pi a^{2}dtf(t^{1/2})/2t^{3/2}$ [1].

The contribution of $S_{e}(\varepsilon)$ in Eq. (2) is dominant compared to
$W(\varepsilon)$. For the latter we use the parametrization given by Wilson et
al. [14], $W(\varepsilon)=1/[4+A\varepsilon^{-B}+C\varepsilon^{-D}]$, where
the constants $A,B,C$, and $D$ depend on the choice of the interatomic
potential, characterized by the choice of the screening function $\phi(x)$. In
this work, we consider four such potentials: Tomas Fermi, Molière, and Average
from Ref. [14], and Ziegler from [13], each giving a model for
$u(\varepsilon)$ and $f(t^{1/2})$. For $S_{e}(\varepsilon)$ we consider three
detailed models, which discussed below.

Electronic stopping power ($S_{e}$). Lindhard’s theory of the electronic
stopping power [15] is appropriate for atomic collisions down to energies of
the order of a few keV. It assumes point-like interactions between an incident
atom with $Z$ electrons, whose velocity is not affected by interatomic
potentials, and the electron cloud of the target atom, which is approximated
by a Fermi gas. An important assumption in the theory is that all the
electronic states of the degenerate gas with energies up to the Fermi energy
$E_{F}$ of the system, participate in the ionization process. Lindhard also
disregards Coulomb repulsion effects between colliding atoms, which at low
energies ($<1$ keV) prevent electron clouds from penetrating each other
completely, allowing only interactions in regions of low electron density
[16]. When taken into account, the Coulomb repulsion effects lead to an
appreciable departure from the proportionality with velocity of the electronic
stopping power assumed by Lindhard.

In general, a detailed model of electronic stopping requires a non-
perturbative analysis of electron dynamics [17]. When the moving ion collides
with an atom in the lattice, multiple inner electronic transitions and
electron promotion will occur, contributing to electronic stopping even for
low sub-keV energies [18]. Two general approaches are considered in the
literature to study electronic stopping: kinetic theory [19, 20] and the
dielectric function approach, pioneered by Lindhard [21]. From the first
category, in this work we consider the models due to Tilinin [22, 23] and
Kishinevsky [24, 25, 26, 27, 28], and from the second, the model due to Arista
[29, 30]. All three are theoretically well founded semi-classical models,
which assume that the kinetic energy of the incoming ion is always sufficient
to effectively excite an atom of the material by electron promotion effects.
They also rely on the determination of the distance of closest approach of the
ion in the interatomic potential. The three models are compared in Fig. 1 for
values of $\varepsilon$ between $0.0006$ and $80$ (about 20 eV to 3 MeV for
Si). Details on our implementation of these models can be found in the
supplemental materials.

Scaling length. A general feature of electronic stopping studies is the use of
the Thomas Fermi (TF) model to define the scaling length $a$, introduced
above. All the atomic electrons in the free gas, with average kinetic energy
$3/5E_{F}$, are assumed to participate in the momentum transfer, an
approximation valid only for sufficiently high energies. However, as noted in
[22], at low energies, the excitation energies of the electrons are typically
much smaller than $E_{F}$ and, since transitions between occupied states are
forbidden, the only electrons participating in the momentum transfer are those
close to the Fermi level, with kinetic energy $\sim E_{F}$. This results in a
correction to the TF scaling length by a factor of $5/3$ at low energies. For
models like those of Refs.[22], [24], and [29], where
$S_{e}\propto\varepsilon^{1/2}\tau(\varepsilon,Z_{1},Z_{2})$, with $\tau$
being a dimensionless function incorporating the effects of Coulomb repulsion,
it can be shown that changing the scaling length introduces a change in
$S_{e}$ dominated by a multiplicative factor of $\xi_{e}=(5/3)^{1.5}\approx
2.152$. For high energies, in principle, $\xi_{e}=1$. Averaging the low and
high energy effects on the scaling $a$, an effective value of $\xi_{e}$
between 2.152 and 1 may be used for all energies. Notice, however, that
Lindhard introduced by hand [31] a semi-empirical factor of $\approx Z^{1/6}$
(1.55 for Si) to match the stopping power data available at his time, which
could be interpreted as the effective value of $\xi_{e}$ suitable for high
energies.

A change in $a$ will also change the interatomic potential, by a factor of
$\xi_{e}^{-2/3}$, and, therefore, introduce modifications in any quantity
depending on it, such as the nuclear and electronic stopping, and the variable
binding energy in the model discussed below. In solving Eq.(2) we will take
into account the uncertainty in the value of $\xi_{e}$, considering it as an
additional parameter of the electronic stopping, i.e. changing
$S_{e}(\varepsilon)$ to $S_{e}(\varepsilon,\xi_{e})$.

Figure 1: Electronic stopping power $S_{e}(\varepsilon,\xi_{e})$, with
$\xi_{e}=1.26$, from Tilinin, Kishinievsky, and Arista, for the various
interatomic potentials: Ziegler (solid), TF (dash-dot), Moliere (dotted),
Average (dash-dot-dot). The gray band shows the range of $\xi_{e}$ consistent
with the spread of QF measurements in Fig. 3. Available data for Si-on-Si [32,
33, 34] are shown as points with error bars . The red line in the right panel
is the Arista-Ziegler combination with $\xi_{e}=1.46$, which best fits the
data. The Lindhard stopping power ($\propto\leavevmode\nobreak\
\varepsilon^{1/2}$) is shown as the constant dashed line.

Bohr stripping and Z oscillations. When a projectile of velocity $v$ collides
with a target atom, accoring to Bohr’s criterion [35], the projectile is
stripped of those electrons with orbital velocities smaller than $v$. This
leads to ions having an effective number of electrons $Z^{\dagger}<Z_{\alpha}$
($\alpha=1,2$), which can be computed using the condition $v\leq\sqrt{3}v_{F}$
[36]. In terms of the screening function $\phi(x)$,
$Z^{\dagger}=Z_{\alpha}[1-\phi(x_{0})+x_{0}\phi^{\prime}(x_{0})]$, where
$x_{0}$ is the distance at which $v=\sqrt{3}v_{F}(x_{0})$. Friedel [37]
predicted that for low energies ($\ll$ 30 MeV), the change in electron density
around an ion or impurity induces an oscillatory behavior of the electron
transport cross section as a function of $Z^{\dagger}$, which is inherited by
the electronic stopping power of the material. For Si we used the oscillatory
factor computed in [38] normalized to its value at $Z_{1}=12$ and applied it
as an energy dependent multiplicative correction to the electronic stopping.
For each ion energy an effective $Z^{\dagger}$ is extracted from Bohr’s
criterion. We notice that for Si, including this effect results in good
agreement with measurements of the electronic stopping power in the interval
from 680 eV to 3 MeV, with effective values of $\xi_{e}$ ranging between 1.3
and 1.7, consistent with Lindhard’s observed value, see Fig. 1.

Binding energy model. Part of the energy transferred to a target atom in the
material by an incident ion is used to excite atomic degrees of freedom and
extract the atom from its site in the lattice. The energy lost in these
processes, $U$ in Eq.(1), is in general dependent on the incident ion energy,
$E$. At low energies it corresponds to the energy required to remove the atom
from a lattice site creating a Frenkel-pair (FP) [8], while at higher energies
it is dominated by the energy spent in disrupting the atomic binding in the
outer valence orbitals. To account for both effects, we define an energy
dependent binding energy function $U(E)=U_{\rm{FP}}+U_{\rm{EDFT}}(E)$, shown
in Fig. (2). The first constant term is the FP creation energy, for which we
consider the interval $U_{\rm{FP}}=23.6^{+15.2}_{-12.1}$ eV, where the central
value is consistent with the weighted average of the predictions from
Hölstrom, Stallinger-Weber, and Tersov quoted in [39].The lower and upper
limits cover the wide variability of both, experimental measurements and
predictions, including the highest range reported from recent simulations
based on electron density functional theory (EDFT). For the second term of
$U(E)$, the energy stored in the electron cloud of an atom in the EDFT model,
from the distance of closest approach, $x(E,\xi_{e})$, to infinity, with a
given choice for the interatomic potential, is considered a measure of the
energy absorbed by the struck atom prior to its recoil in the sudden
approximation. This energy is computed in terms of the screening function
$\phi$ as

$U_{\rm{EDFT}}(E,\xi_{e})={U_{0}}\int^{\infty}_{x}dy\>\phi^{5/2}(y,\xi_{e})/y^{1/2}$
(3)

where ${U_{0}}=18.6/\xi^{4/3}_{e}\;Z^{1.23}$ eV.

Solving the equation with variable $u(\varepsilon)$. Generalizing the ideas in
[5] for an energy dependent $u(\varepsilon)$, it can be seen that Eq. (2) is
only applicable for $\varepsilon\geq u(\varepsilon)$. When the equality holds,
the rhs of Eq.(2) is equal to zero, since the upper and lower limits of
integration are equal. This will occur at the threshold energy
$\varepsilon^{*}=u(\varepsilon^{*})$, which for silicon is only slightly
greater than the FP creation energy. The appropriate generalization of the
parameterized solution for $\bar{\nu}(\varepsilon$) studied in [5] is the
following:

$\bar{\nu}[s(\varepsilon)]=\left\\{\begin{array}[]{ll}s(\varepsilon)+u(\varepsilon),&s(\varepsilon)\leq
u(\varepsilon),\\\
s(\varepsilon)+u(\varepsilon)-\lambda[s(\varepsilon)],&s(\varepsilon)\geq
u(\varepsilon),\\\ \end{array}\right.$ (4)

where $s(\varepsilon)$ is any well behaved function of $\varepsilon$, such as
those appearing as arguments of $\bar{\nu}$ in the rhs of Eq.(2) and $\lambda$
is a continuous function satisfying $\lambda(\varepsilon^{*})=0$, but with
discontinuous derivatives at this point. This form will guarantee that at
$\varepsilon^{*}$ the integrand in the rhs of Eq.(2) goes to zero, as is
required by energy conservation. In the case $s(\varepsilon)=\varepsilon$, the
discussion in [5] leading to the implementation of the shooting method for
solving the equation for $\bar{\nu}(\varepsilon)$, applies for varying
$u(\varepsilon)$ as well, with minor modifications. For values below
$\varepsilon^{*}$ it suffices that $\bar{\nu}(\varepsilon)$ take the form
$\bar{\nu}(\varepsilon)=\varepsilon+u(\varepsilon)=\varepsilon_{R}$. For
$\varepsilon\geq\varepsilon^{*}$ the solution is found numerically. To take
into account the effect that the moving ion will lose some of its valence
electrons as a consequence of its motion through the lattice, following [40]
we consider that the effective number of valence electrons in the incident ion
is $Z_{1}=12$, while $Z_{2}=14$. Given that the Bohr Stripping effect
describes the experimental data of electronic stopping power up to energies of
the order of 3 MeV [41], we evaluate the solution from $U_{FP}$ to 3 MeV, in
order to compare to QF measurements available at those energies (Sattler [3]).

Figure 2: Variable binding energy model with $\xi_{e}=1.26$, based on EDFT,
see Eq.(3) , for the four interatomic potencials. Includes both, the Frenkel-
pair creation energy in Si and the inner electron excitations.

Results. We find the numerical solution to Eq.(2) for the three different
models of the electronic stopping power and the four interatomic potentials
introduced. For each combination we evaluate the solution in a grid of values
of $U_{FP}$ and $\xi_{e}$ constrained to lie within the physically motivated
intervals discussed earlier. Each curve was compared to the available QF data
in Si summarized in Table II of [5], to which we have added the Sattler [3],
Gerbier [42] and Agnese [43] data points. Overall, the data ($N_{data}=102$
points) span a range of energies from 680 eV to 3 MeV. We calculated the
$\chi^{2}$ for each curve as
$\chi^{2}=\sum_{i=1}^{N_{data}}(D_{i}-C_{i})^{2}/\sigma_{i}^{2}$, where
$D_{i}$ and $\sigma_{i}$ are the QF value and its error for data point $i$,
and $C_{i}$ is the value predicted by the curve at the same energy. Fixing the
Frenkel pair energy at $U_{FP}=23.6$ eV, and varying $\xi_{e}$, we obtain that
the model that best fits the data, shown as the black solid line in Fig. 3
corresponds to the combination of the Arista stopping power with Ziegler
potential (Arista-Ziegler), with $\xi_{e}=1.26$ ($\chi^{2}/d.o.f=1870/101$).
The high $\chi^{2}$ values are expected from the tension among the different
data sets. We set an uncertainty on $\xi_{e}$ of $\pm 0.25$
($\Delta\chi^{2}\gtrsim 3100$) so as to approximately span the spread in the
QF measurements. Notice that the electronic stopping data in Si at higher
energies prefer a larger value of $\xi_{e}$, around 1.46, but are consistent
within the uncertainty interval determined for the QF measurements (see Fig.
1). The $\chi^{2}$ has a weak dependence on $U_{FP}$, and the solution at the
center of the physically motivated interval has nearly the same goodness of
fit than that at the edges ($\Delta\chi^{2}\approx 40$). When $U_{FP}$ is
allowed to vary freely, the fit prefers the lower edge of the interval.
Although the main effect of $U_{FP}$ is to set the threshold, it also has a
small effect in the detailed shape of the curve at intermediate energies
(1-100 keV).

Figure 3:  Published Si QF data (points with error bars) compared to the
Arista-Ziegler model with $\xi_{e}=1.26\pm 0.25$ and
$U_{FP}=23.6^{+15.2}_{-12.1}$ eV (black solid line with gray band). Also shown
are the Arista-Ziegler model with $\xi_{e}=1.46$ (black dashed line), which
best fits the $S_{e}$ data in Fig.1, and the Lindhard model (dotted line).

In Fig. 3 we show the Arista-Ziegler model with $\xi_{e}=1.26$ and
$U_{FP}=23.6$ eV (black solid line). The gray band is drawn varying $\xi_{e}$
between 1.0 and 1.5, and varying $U_{FP}$ in the interval 11.5-38.8 eV. For
comparison, we also show the model with $\xi_{e}=1.46$ (dashed black line),
which gives the best fit to the $S_{e}$ for Si-on-Si data [32, 33, 34] in Fig.
1, for the Arista-Ziegler combination. The Sarkis 2020 [5] model (constant
binding energy) is shown in the red line.

In the supplemental materials we provide tabulated versions of the curves in
Fig. 3, as well as a comparison of the curves obtained with fixed $\xi_{e}$
and $U_{FP}$ for the twelve combinations of models of the stopping power and
interatomic potential.

Conclusions. In this letter, we present new results on the study of the QF for
low-energy nuclear recoils in Si that introduce significant improvements over
previous work based on the Lindhard integral equation for the energy given to
atomic motion [5]. We use more accurate descriptions of the electronic
stopping power with several interatomic potentials. Instead of a constant
binding energy, we treat it as a function of the initial ion energy, dependent
also on the choice of the interatomic potential. At low energies, this
function is taken to correspond with the energy required to create a stable
Frenkel pair. The high energy effects of Bohr stripping and Friedel
Z-Oscillations were also incorporated and provide, simultaneously, a good
match to the QF data at higher energies, also constraining the low energy
behavior. The main uncertainty on the models comes from the scaling parameter
chosen to define the dimensionless energy scale $\xi_{e}$, affecting both, the
low and high energy regimes. The value of the Frenkel pair creation energy
plays a significant role at the lowest energies and might be more relevant
when forthcoming data becomes available. Slight tension between the sparse
measurements of the Si-on-Si electronic stopping power at high energies and
the QF data might be resolved with more measurements performed in both the low
and high energy regimes. Our model could be improved allowing the scaling
factor to change with the ion energy. Phonon excitations, not included here,
could also have an effect near the threshold. With these improvements, and if
more data of the electronic stopping becomes available it could be used in a
combined fit with the QF data. The found solution is able to describe
reasonably well all the available published data for the nuclear recoil QF in
Si from 50 eV to 3 MeV of nuclear recoil energy, and provides a unified
theoretical picture of the effects of electronic and nuclear stopping from
high energies down to the Frenkel pair creation threshold on the ionization
efficiency in Si.

Acknowledgements This research was supported in part by DGAPA-UNAM grants
number PAPIIT-IN106322 and PAPIIT-IT100420, and Consejo Nacional de Ciencia y
Tecnología (CONACYT) through grant CB-2014/240666.

## References

  * Lindhard  _et al._ [1963] J. Lindhard, V. Nielsen, M. Scharff, and P. Thomsen, Kong.Dan.Vid.Sel.Mat.Fys.Med. 33, 10 (1963).
  * Dougherty [1992] B. L. Dougherty, Phys. Rev. A 45, 2104 (1992).
  * Sattler [1965] A. R. Sattler, Phys. Rev 138, 1815 (1965).
  * Zecher  _et al._ [1990] P. Zecher, D. Wang, J. Rapaport, C. J. Martoff, and B. A. Young, Phys. Rev. A 41, 4058 (1990).
  * Sarkis  _et al._ [2020] Y. Sarkis, A. Aguilar-Arevalo, and J. C. D’Olivo, Phys. Rev. D 101, 102001 (2020).
  * Chavarria  _et al._ [2016] A. E. Chavarria _et al._ , Phys. Rev. D94, 082007 (2016), arXiv:1608.00957 [astro-ph.IM] .
  * Izraelevitch  _et al._ [2017] F. Izraelevitch _et al._ , JINST 12 (06), P06014, arXiv:1702.00873 [physics.ins-det] .
  * Frenkel [1926] J. Frenkel, Zeitschrift für Physik 35, 652 (1926).
  * J.Lindhard [1954] J.Lindhard, Matematisk-fysiske Meddelelser Vol: 28, No. 8 (1954).
  * Sigmund [2006] P. Sigmund,  _Particle Penetration and Radiation Effects_ (Springer-Verlag Berlin Heidelberg, 2006).
  * Migdal [1941] A. B. Migdal, J. Phys 4, 449 (1941).
  * Baur  _et al._ [1983] G. Baur, F. Rosel, and D. Trautmann, Journal of Physics B: Atomic and Molecular Physics 16, L419 (1983).
  * F. Ziegler, J.P. Biersack, and U. Littmark [1985] F. Ziegler, J.P. Biersack, and U. Littmark, Pergamon Press New York (1985).
  * Wilson  _et al._ [1977] W. D. Wilson, L. G. Haggmark, and J. P. Biersack, Phys. Rev. B 15, 2458 (1977).
  * Lindhard and Scharff [1961] J. Lindhard and M. Scharff, Phys. Rev. 124, 128 (1961).
  * Semrad [1986] D. Semrad, Phys. Rev. A 33, 1646 (1986).
  * Ovchinnikov  _et al._ [2004] S. Ovchinnikov, G. Ogurtsov, J. Macek, and Y. Gordeev, Physics Reports 389, 119 (2004).
  * Fano and Lichten [1965] U. Fano and W. Lichten, Phys. Rev. Lett. 14, 627 (1965).
  * Sigmund [1982a] P. Sigmund, Phys. Rev. A 26, 2497 (1982a).
  * [20] B. Trubnikov and Yavlinskii, Zh. Eksp. Teor.Fiz. 48, 253.
  * J.Lindhard [1961] J.Lindhard, Kong.Dan.Vid.Sel.Mat.Fys.Med. 28 (1961).
  * Tilinin [1988] I. S. Tilinin, Zh. Eksp. Teor. Fiz. 94, 96 (1988).
  * Sigmund [1982b] P. Sigmund, Phys. Rev. A 26, 2497 (1982b).
  * Kishinevsky [1962] L. Kishinevsky, Izv. Akad. Nauk SSSR 26, 1410 (1962).
  * Firsov [1959] Z. E. Firsov, O.B., Teor. Fiz. 36, 1517 (1959).
  * Robinson [1974] M. Robinson, Third National Conference on Atomic Collisions with Solids, Kiev USRR (1974).
  * S. Oen and T. Robinson [1976] O. S. Oen and M. T. Robinson, Nuclear Instruments and Methods 132, 647 (1976).
  * Littmark [1976] U. Littmark, JPP 9/21 (Max-Planck Institut fur Plasmaphysik, 1976).
  * Fernández-Varea and Arista [2014] J. M. Fernández-Varea and N. R. Arista, Radiation Physics and Chemistry 96, 88 (2014).
  * Sigmund [2014] P. Sigmund,  _Particle Penetration and Radiation Effects Volume 2_ (Springer International Publishing, 2014).
  * Sigmund [2008] P. Sigmund, Bulletin of the Russian Academy of Sciences: Physics 72, 569 (2008).
  * SeD [2000] Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms https://doi.org/10.1016/S0168-583X(00)00050-1 (2000).
  * [33] N.Hecking, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 59, 619.
  * Zhang and Weber [2003] Y. Zhang and W. J. Weber, Phys. Rev. B 68, 235317 (2003).
  * Gervasoni and Cruz-Jiménez [1996] J. L. Gervasoni and S. Cruz-Jiménez, Radiation Physics and Chemistry 48, 433 (1996).
  * [36] H. Gümüş and Önder Kabaday, Journal of Physical Science and Application 3.
  * Friedel [1952] J. Friedel, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 43, 153 (1952), https://doi.org/10.1080/14786440208561086 .
  * Peñalba  _et al._ [1992] M. Peñalba, A. Arnau, and P. Echenique, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 67, 66 (1992).
  * Holmström  _et al._ [2008] E. Holmström, A. Kuronen, and K. Nordlund, Phys. Rev. B 78, 045202 (2008).
  * Ying-Tai  _et al._ [2002] L. Ying-Tai, Z. Qi-Ren, and G. Chun-Yuan, Communications in Theoretical Physics 38, 361 (2002).
  * Sarkis  _et al._ [2021] Y. Sarkis _et al._ , Magnificent CEVNS 2021, https://indi.to/DqCzz (6 to 7 Oct. 2021).
  * Gerbier  _et al._ [1990] Gerbier _et al._ , Phys. Rev. D 42, 3211 (1990).
  * Agnese  _et al._ [2018] R. Agnese _et al._ , Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 905, 71 (2018).


# Physics Informed Machine Learning (PIML) methods for estimating the
remaining useful lifetime (RUL) of aircraft engines
††thanks: Work performed in Southwest Research Institute (SwRI)

Sriram Nagaraj  Southwest Research Institute (SwRI)
San Antonio, Texas, USA
<EMAIL_ADDRESS>Truman Hickok  Southwest Research Institute (SwRI)
San Antonio, Texas, USA
<EMAIL_ADDRESS>
###### Abstract

This paper is aimed at using the newly developing field of physics informed
machine learning (PIML) to develop models for predicting the remaining useful
lifetime (RUL) aircraft engines. We consider the well-known benchmark NASA
Commercial Modular Aero-Propulsion System Simulation (C-MAPSS) data as the
main data for this paper, which consists of sensor outputs in a variety of
different operating modes. C-MAPSS is a well-studied dataset with much
existing work in the literature that address RUL prediction with classical and
deep learning methods. In the absence of published empirical physical laws
governing the C-MAPSS data, our approach first uses stochastic methods to
estimate the governing physics models from the noisy time series data. In our
approach, we model the various sensor readings as being governed by stochastic
differential equations, and we estimate the corresponding transition density
mean and variance functions of the underlying processes. We then augment LSTM
(long-short term memory) models with the learned mean and variance functions
during training and inferencing. Our PIML based approach is different from
previous methods, and we use the data to first learn the physics. Our results
indicate that PIML discovery and solutions methods are well suited for this
problem and outperform previous data-only deep learning methods for this data
set and task. Moreover, the framework developed herein is flexible, and can be
adapted to other situations (other sensor modalities or combined multi-physics
environments), including cases where the underlying physics is only partially
observed or known.

###### Index Terms:

Predictive maintenance, Physics Informed Machine Learning (PIML), Deep
Learning

##  I Introduction

Aircraft maintenance is a a critical component in ensuring the safety,
dependability, and efficiency of aviation systems and their operational modes
([1, 2, 3, 4]). Classical maintenance methods, such as scheduled maintenance,
aimed at corrective repair/upkeep may still be relevant but are often plagued
by limitations such as operational (including logistic) disruptions, periods
of unavailability (downtime) and increased costs, for even routine issues. On
the other hand, predictive maintenance (PM), which refers to maintenance
approaches relying on data-driven or data-augmented techniques, is a proactive
set of approaches that leverage past data (and modes) to foresee system (or
subsystem) failures prior to occurrence ([5, 6]). Quantities of interest that
may be part of a PM diagnosis are probability of engine failure/stalls,
estimating remaining useful lifetime (RUL) of engines, probability of
catastrophic multi-system failure etc. A characteristic of PM systems are
real-time (or near real-time) information about the system and dynamic
updating of data to provide estimates when the mission or flight is in
progress. Moreover, these systems can be incorporated into the flight system
itself and not only in separate command centers.

###  I-A Salient Features:

Other important points to note about PM are:

  * •

Since PM systems are (near) real-time, they allow for early detection and
identification of aircraft degradation or faulty functioning of aircraft
subsystems. In the case of commercial aircraft, this allows to potentially
save passenger lives by reducing the risk of dangerous malfunction being
unnoticed. PM allows for corrective, risk-mitigation maneuvers on the part of
pilots ([5]).

  * •

When early onset of faulty system performance is picked by PM diagnosis, this
allows for action to be taken that would help avoid further catastrophic
outcomes such as complete engine failure. This would help not only saving
human life, but also reduce total maintenance costs of the mission/flight ([3,
4]).

  * •

PM allows for addressing maintenance needs during both scheduled and
unscheduled periods. This in turn minimizes any unplanned downtime of aircraft
systems and helps to optimize aircraft utilization.

  * •

Continuing the earlier remarks, early detection of system malfunction allows
for the use of maintenance resources for mission-critical tasks. Thus, PM
ensures the focusing of in-depth maintenance efforts on aircraft/avionics
systems or components that are known to require attention, thereby optimizing
the allotment of maintenance resources. Finally, this helps to reduce costs
overall in maintenance.

  * •

As another means of reducing overall costs, PM reduces the chances of
unexpected maintenance issues, and thereby reduces operational delays (such as
flight delays, cancellation or diversions) which are known to impact
commercial airlines significantly. For military aircraft, operational delays
may impact the efficacy of the mission, and may even lead to a critical
mission being compromised ([5]).

  * •

PM and ML/AI models: PM relies on innovative use of sensor data from various
aircraft systems, such as engines, avionics, and airframes. These sensor
readings are fused via ML/AI techniques and processed to predict anomalous
behavior, potential known or unknown failures; further, in non-critical
scenarios, the ML/AI methods are used to monitor equipment health. Simply put,
the ML/AI techniques are trained on a large swathe of past sensor data in
order to detect/identify trends and risk factors that are indicative of
equipment failures or system degradation. Once detected or identified,
information about these issues can be communicated to a subject matter expert
(SME) to determine the severity of the issues and what corrective actions need
to be taken ([2]).

  * •

Regulatory agencies such as the Federal Aviation Administration (FAA) and the
European Union Aviation Safety Agency (EASA), have developed standards and
regulations that need to be met by commercial aircraft in order for them to be
flight worthy. PM techniques can be used to augment existing risk management
methods to further ensure compliance with these regulatory standards.

There are several challenges to full-scale incorporation of PM methods. First,
data quality/integrity and interoperability challenges often prevent the
implementation of PM methods. Scarce data or uneven data between different
sensors limit the applicability of PM techniques. Integrating real-time
systems on platform is a separate challenge, as the existing hardware and
infrastructure on aircraft needs to be sufficient for running PM diagnostics.
Finally, reliability of PM methods requires extensive testing, and both
commercial and military operators need to allocate sufficient resources to
adequate testing. However, these challenges can be surmounted with careful
analysis and a well-defined systems engineering approach during design and
development. The future directions of a PM based fleet are to include internet
of things (IoT) devices into PM data sources and the design of digital twins
to further improve the efficacy of PM.

PM is therefore an emerging, and continually evolving paradigm in aircraft
maintenance practices ([1]). PM techniques can be incorporated with other
evolutions in the aviation industry with newer hardware/sensing capabilities.
PM can then co-evolve to ensure the safety and reliability of ever increasing
complex aircraft systems. Moreover, the financial savings afforded by PM can
be used to further increase aircraft capabilities.

###  I-B Physics Informed Machine Learning

Currently, a large amount of computational modeling is based on deterministic
models and/or multi-fidelity approximate solutions that do not directly
account for inherent variability (stochasticity) in the material properties.
When predicting real-world expected performances, these models are laboriously
recomputed many times and deliver solutions that do not generalize efficiently
to situations where the data needed for simulations is partially observed.
Moreover, real experimental test data is a scarce commodity, so there is a
need to utilize such data in the most efficient way possible ([7]).

####  I-B1 PIML in a Nutshell

Physics Informed Machine Learning (PIML) models offer a promising way forward
for balancing computational fidelity and data-driven generalization. PIML is a
data driven ML modeling framework that is augmented by the underlying
scientific domain (mechanics, electromagnetics, acoustics, thermal or
multiphysics). Thus, scientific knowledge captured by some relevant set of
physical principle/equation (conservation laws, Maxwell equations, mechanical
properties, heat diffusion etc.) is directly integrated with the ML modeling
to predict quantities of interest (failure of materials, energy transmission
efficiency, stress, drag etc.) which are often needed for mission-critical
decision making ([7, 8]). PIML models achieve better generalization than
traditional “data only” ML models since they can be viewed as incorporating
real-world prior knowledge (i.e. a Bayesian reformulation). This also allows
PIML models to operate in low/partially observed data regimes where the
physics-informed component can act as data surrogate as needed. Due to
inherent regularization, PIML models are less prone to spurious solutions
while dealing with noisy data as well. Finally, PIML models yield more
succinct models characterized by fewer ML parameters while achieving (often
exceeding) the accuracy/efficiency of complex “data only” deep learning
models. This also renders then more explainable than “black box” ML models and
is significant when downstream decision making requires explainable
predictions ([9, 10]).

####  I-B2 Methodologies in Physics-Informed Machine Learning

There are several ways in which PIML performed.

  * •

Differential Equation Constraints: PIML frameworks utilize differential
equations to incorporate dynamic system behavior into machine learning models,
facilitating time-series forecasting, and dynamical system identification
([10]).

  * •

Generative Data: Systems governed by stochastic differential equation (SDE)
models can be used to perform physics-informed data generation. This generated
data would capture well the underlying physics, and being a data augmentation
method, can directly be incorporated into the ML modeling to ensure model
robustness ([11]).

  * •

Hybrid Physics-Data Models: Hybrid approaches are also possible which
integrate physics-based simulations with data-driven (ML) techniques,
balancing model accuracy and computational efficiency ([9]).

####  I-B3 Modes of PIML Operation

Roughly, PIML methods operate in two modes. The first is the “discovery mode”
wherein the physical parameters (coefficients of differential equations,
structural parameters etc.) are learning from the underlying data. The second
is the “solution mode” which uses a fully specified physical model along with
a data-driven component to drive predictions. The exact means by which the two
modes operate are largely problem/domain dependent, and specific
implementations of either mode of operation may change widely across different
areas ([11, 8]).

####  I-B4 Physics Informed Neural Networks (PINNs)

PINNs are an exemplary class of PIML models. PINNs are deep neural networks
which are trained on a loss functions composed of two components ([8]). The
first is the usual data driven (often mean square error MSE), while the other
component is a physics-informed one. The physics informed loss includes all
physical constraints of the problem include modeling choices,
differential/integral equations, boundary conditions, uncertain coefficients
etc. Both components of the overall loss work in tandem to deliver a robust
model. While PINNs are one example of a class of PIML models, it is possible
to have other PIML approaches as well. In this regard, _PIML is a modeling
framework, rather than a specific set of model types_.

####  I-B5 Applications of PIML

  * •

Aerospace Engineering: PIML methods are very useful for turbulence modeling,
wing design and shape optimization, where traditional computational methods
such as finite element/finite difference methods can be augmented by data-
driven models ([12, 9]).

  * •

Mechanical Engineering: Applications in mechanical engineering includes
modeling of fatigue of materials, CFD, elasticity etc. ([13, 14, 9])

  * •

Material Science: PIML accelerates materials discovery by integrating physics-
based simulations with high-throughput experimentation, guiding the design of
novel materials with desired properties ([9, 15]).

  * •

Biomedical Engineering: PIML aids in medical image analysis, patient
monitoring, and drug discovery, leveraging physics-based models and clinical
data to improve diagnosis and treatment outcomes ([16, 17, 18]).

Scalability of PIML methods to high-dimensional data remains a challenge.
Lower order modeling, model order reduction methods, and surrogate models have
proven to be promising ways of mitigating this challenge ([19]). Another area
requiring further investigation is uncertainty quantification (UQ) of PIML
predictions. UQ is needed for meaningful decision-making and risk assessment,
and robust methods for uncertainty quantification and propagation are
essential. Finally, software packages for easy implementation of differential
equation solvers using deep learning methods are available as well ([20]).

###  I-C Prior Work: Data-Driven Deep Learning for PM using the C-MAPSS
Dataset

In recent years, machine learning techniques, particularly long short-term
memory (LSTM) and convolutional neural networks (CNNs), have shown promise in
PM applications. Towards this end, the Commercial Modular Aero-Propulsion
System Simulation (C-MAPSS) dataset (a widely used benchmark dataset for
prognostics and health management research) is a standard dataset used for
studying the efficiency of ML models for PM ([21, 22]). C-MAPSS consists of
multiple simulated engine run-to-failure trajectories, each containing sensor
measurements and corresponding remaining useful life (RUL) labels. The data
consists of multiple noisy time series that correspond to sensor readings from
different parts of an aircraft engine. Thus far, the LSTM Networks and CNNs
have been well-studied architectures for PM usinng the C-MAPSS dataset ([23,
24, 25, 26]). LSTMs are a type of recurrent neural network (RNN) designed to
capture long-term dependencies in sequential data. LSTMs consist of memory
cells, input and forget gates, and output gates, allowing them to retain
information over long sequences and adapt to varying temporal patterns. CNNs
on the other hand, are deep learning architectures designed for feature
extraction from structured data, such as images or time-series signals. CNNs
consist of convolutional layers, pooling layers, and fully connected layers,
enabling hierarchical feature learning and pattern recognition. CNNs are
applied to sensor data for fault detection, anomaly detection, and condition
monitoring, leveraging their ability to automatically learn relevant features
from raw data. LSTMs have been successfully applied to time-series data for
equipment health monitoring, remaining useful life estimation, and fault
diagnosis. CNNs have been applied to sensor data for fault detection, anomaly
detection, and condition monitoring, leveraging their ability to automatically
learn relevant features from raw data. However, these methods are purely data-
driven and do not incorporate the physics underlying the data.

(a)

(b)

Figure 1: (a,b) Sample sensor outputs for two different modes of operations
(FD001 and FD003 respectively). (a) is well modeled as geometric brownian
motion, while (b) shows multi-modal behavior of sensor output.

###  I-D Our Contributions:

We consider the C-MAPSS dataset as the main data for this paper. As we have
seen, this is a well-studied dataset, and much work exists in the literature
that address RUL prediction on C-MAPSS data with classical and data-driven
(deep) ML methods. Our PIML based approach is different from previous methods,
and we use the data to first learn the physics. In the absence of published
empirical physical laws governing the C-MAPSS data, our approach first uses
stochastic PIML methods to estimate the governing stochastic physics models
from the noisy time series data (this is referred to as the “discovery” mode
of PIML). In this phase, we model the sensor data as being governed by a
system of first order stochastic differential equation (SDE) and estimate the
mean and variance functions. When the underlying statistical distribution at a
given timestep is known to be unimodal, we can use the sample mean and
variance as estimates of the mean and variance functions at that timestep, or
equivalently, the estimate of the $K=1-$mean cluster median and variance. If
the distribution is multimodal (with $K-$modes) we use $K$-means clustering
estimation to define the estimated mean and variance functions at the timestep
as the $K-$medians of the collection of sample paths at the timestep in
question. We are thereby able to obtain a surrogate for the underlying
physics. $K$-means estimation is a well-established nonparametric method, and
using it, we are able to model multi-modal sensor output wherein the
statistical distribution of the same sensor output is a multi-modal
distribution (see Figure 1(b) and Figure 2(b)). We then use existing ML models
(LSTM) with the learned physics to generate predictions of RUL (the “solution”
mode of PIML). Unlike traditional PIML methods, in the solution mode, we
directly incorporate the mean and variance functions as data augmentation
rather than indirectly through a loss function. The proposed set of methods
delivers reliable, robust solutions to the challenging problems of physics
discovery and physics informed prediction, relying on the superior
generalization capabilities of modern deep machine learning. Our results
indicate that PIML discovery and solutions methods are well suited for this
problem and outperform previous data-only deep learning methods for this data
set and task. Moreover, the framework developed herein is flexible, and can be
adapted to other situations (other sensor modalities or combined multi-physics
environments), including cases where the underlying physics is only partially
observed or known.

##  II Stochastic Physics Model

We now describe our physics model for the different datasets. Recall the
absence of existing physics based models for the data. We therefore appeal to
the approach of _learning_ the physics from the raw data. To this end, we
model the output $S(t)$ of each (viable, information bearing) sensor via a
stochastic differential equation (SDE). Note that, inspecting the sample paths
of the sensors from FD001-FD004 (see section III for more details on the
dataset) shows that this is a reasonable modeling choice. Indeed, several
sensors are well-modeled as geometric brownian motion. We assume the following
form for each sensor output $S:$

$dS(t)=a(t)dt+b(t)dW(t),$ (1)

where $W(t)$ is a one-dimensional martingale stochastic process (not
necessarily Gaussian, but with independent increments) and $a(t),b(t)$ are the
_drift and diffusion_ coefficients respectively. The stochastic integrals are
assumed to be in the sense of Ito. Since $S(t)$ is a stochastic process, we
will, as needed, explicitly denote the dependence on a random outcome in the
event space $\omega$ as $S(t,\omega)$. In reality, the drift and diffusion
functions are themselves dependent on both $t$ as well as $S(t),$ i.e.,
$a(t)=\tilde{a}(t,S(t))$ and $b(t)=\tilde{b}(t,S(t))$, although for our
purposes, we shall denote them as functions of $t$ alone, with the dependence
on $S(t)$ being implicit. Thus, each data set consists of multiple cycles of
several sensor values, $S_{i}(t)$ with $i=1,\ldots,N.$ One possible approach
would be to estimate the corresponding $a_{i}(t),b_{i}(t),$ from the data
directly to fully specified the SDE model. However, we take a slightly
different approach. We instead estimate the mean and variance process defined
by:

$\displaystyle\mu(t)$ $\displaystyle=\mathbb{E}[S(t)],$ (2)
$\displaystyle\rho(t)$ $\displaystyle=\mathbb{E}[|S(t)-\mu(t)|^{2}].$ (3)

These are the instantaneous mean and variance functions of $S(t)$. Note that
if $p(t,x)$ denotes the density of the process $S(t,\omega),$ i.e.,
$S(t,\cdot)\sim p(t,x)$ we have that $\mu(t)=\int xp(t,x)dx$ and
$\rho(t)=\int(x-\mu(t))^{2}p(t,x)dx.$ We take this approach for the following
reasons. First, an inspection of sensor data as a function of time (over
multiple cycles) indicates that several sensors have non-Gaussian distribution
at different time instances (see for e.g. Figure 1(b)). Typical use cases of
the standard diffusive Ito SDE assume Gaussian, independent increment
processes, which would not be valid in the current case. Second, estimating
the coefficients $a(t),b(t)$ directly from data would need further assumptions
on the underlying process $S(t),$ which, again, may not be valid. Instead, we
stay with the most general assumptions possible, and estimate the mean and
variance functions instead. Here, the only assumption being made is that the
process has finite first and second moments, which is evidently true.

###  II-A Estimating $\mu(t),\rho(t)$

We now discuss the estimation of $\mu(t),\rho(t)$. We first assume that the
distribution of $S(t,\omega)$ is unimodal for all $t$. This is true for
several sensors in the dataset. For a given sample path $S_{i}(t,\omega_{k}),$
we define the support of the path as

$\rho_{S_{i}}(k):=\\{t>0:|S_{i}(t,\omega_{k})|>0\\}.$ (4)

We assume that $\rho_{S_{i}}(k)=[0,T_{S_{i}}(k)]$ for some $T_{S_{i}}(k)>0.$
We set $T_{max}=\text{sup}_{\omega_{k},S_{i}}(T_{S_{i}}(k))$ to be the maximum
length of the vector of any sensor’s readings over all cycles. If a given
sample path $S_{i}(t,\omega_{k})$ ends before $T_{max}$, we extend the sample
path by zero, i.e., we set $S_{i}(t,\omega_{k})=0,T_{S_{i}}(k)<t<T_{max}.$

Now, given a fixed partition $\Delta=\\{0<t_{1},\ldots,t_{n}=T_{max}\\},$ for
each $t_{k}$ we define a random variable as
$\tilde{S_{i}}(k,\omega)=S_{i}(t_{k},\omega)$. That is,
$\tilde{S_{i}}(k,\cdot)$ is obtained by collecting the sample paths
$S_{i}(t,\cdot)$ and evaluating them at $t=t_{k},$ thereby resulting in the
random variable $\tilde{S_{i}}(k,\cdot):\Omega\rightarrow\mathbb{R}.$ We then
define $\hat{\mu}_{i}(k)=\mathbb{E}(\tilde{S_{i}}(k,\cdot)),$ and
$\hat{\rho}_{i}(k)=\mathbb{E}(\tilde{S_{i}}(k,\cdot)^{2})-\hat{\mu}_{i}(k)^{2}$
to be the mean and variance of the random variable $\tilde{S_{i}}(k,\cdot).$
Finally, we can define the estimates of the drift and diffusion functions via
an interpolation operator $\Pi:\mathbb{R}^{n}\rightarrow C([0,T_{max}])$ as
$\hat{\mu}_{i}(t)=\Pi(\\{\hat{\mu}_{i}(k),k=1,\ldots,n\\})$ and likewise
$\hat{\rho}_{i}(t)=\Pi(\\{\hat{\rho}_{i}(k),k=1,\ldots,n\\}).$

In practise, we shall be in need of only the values
$\hat{\mu}_{i}(k),k=1,\ldots,n$ (resp. $\hat{\rho}_{i}(k)$), and hence will
not specify the exact form of the interpolation operator $\Pi(\cdot).$ Indeed,
we preprocess the sensor data and compute the drift and diffusion functions
using an appropriately defined histogram (or density estimate) of the random
variable $\tilde{S_{i}}(k,\omega)$ at each time instant $t_{k}\in\Delta.$
Computing $\hat{\mu}_{i}(k)$ is done using $K-$ means clustering. Indeed,
since we know the distribution is unimodal, we compute the single centroid of
the distribution of $\tilde{S_{i}}(k,\cdot)$ and define it to be the mean
$\hat{\mu}_{i}(k)$. This way, we avoid explicit computation of the
expectation.

####  II-A1 Estimating $\mu(t),\rho(t)$: Multimodal Case

Thus far, we have tacitly assumed that the histogram (density estimate) of
$\tilde{S_{i}}(k,\omega)$ is unimodal when computing the drift and diffusion
functions. In some cases, this is not true, and we need to modify our approach
to include multi-modal densities. In such cases, we resort to modeling the
multi-modal density of $\tilde{S_{i}}(k,\omega)$ using a $K-$means estimate of
the mean. Typically, there are not more than two components in the mixure.
Accordingly, the same approach as earlier can be carried forward (using
$K=2$), with the caveat that we now have two mean function estimates
$\hat{\mu}^{1}(t),\hat{\mu}^{2}(t)$ (resp.
$\hat{\rho}^{1}(t),\hat{\rho}^{2}(t)$). The particular sensors that exhibit
this behavior can be ascertained beforehand, and this multi-modal analysis can
be done only in such cases.

###  II-B Generative Modeling: Synthetic Data Generation

Note that the above approach also yields a _generative_ model for information-
bearning sensors as well. Indeed, in the standard diffusive SDE case, assuming
that $W(t)$ is a Brownian motion, we know from standard SDE theory that the
density $p(x,t)$ of $S(t)$ is a solution of the Fokker-Planck equation:

$\frac{\partial}{\partial t}p(x,t)=-\frac{\partial}{\partial
x}\mu(x,t)p(x,t)+\frac{\partial^{2}}{\partial x^{2}}\rho(x,t)p(x,t).$ (5)

In our current case, we do not have such a clear mathematical form for the
evolution of the density. However, we can sample from the $K$-means density
estimate that we used to compute the mean process at each time instant
$t_{k}\in\Delta$ to obtain generated sample paths. This is a unique feature of
our approach, and yields a physics-informed synthetic data generation
mechanism that can be used to augment the raw (real) sensor data. Figures 3(a)
and 3(b) show how we can generate samples from the time-varying densities
which may possibly be multi-modal. The samples that are generated evidently
have similar mean/variance properties as the “real” data. Notice that, since
$dW(t)$ is assumed to be a zero-mean process, we trivially have

$\displaystyle\frac{d\mu(t)}{dt}$ $\displaystyle=\bar{a}(t),$ (6)
$\displaystyle\frac{d\rho(t)}{dt}$
$\displaystyle=\frac{dR(t)}{dt}-2\mu(t)\bar{a}(t),$ (7)

where $\bar{a}(t)=\mathbb{E}[a(t)],\,R(t)=\mathbb{E}[S(t)^{2}].$

(a)

(b)

Figure 2: (Evolution of the histogram over time for the same sensor. Notice
the shift from unimodal to bimodal behavior.

(a)

(b)

Figure 3: Sample generated (synthetic) data from the estimated distributions.
In bold is the mean of the actual (real) data.

##  III Experiments

###  III-A Dataset

We focus our experiments on the aforementioned C-MAPSS dataset produced by
NASA, which is organized into four subsets corresponding to distinct aircraft
operating conditions. Each subset contains synthetic sensor readings meant to
simulate run-to-failure trajectories of turbofan engines. Each trajectory is
usually a few hundred timesteps. We drop 8 of the 22 available sensors
(sensors 1, 5, 6, 10, 16, 18, 19, and 22) and split each trajectory into
20-timestep windows, which are shuffled during training.

###  III-B Deep Learning Model

Consistent with previous work [23], we use a small LSTM neural network to
predict the RUL at the end of each 20-timestep window. Specifically, we use a
single-layer LSTM with 12 hidden units and apply a two-layer MLP to the LSTM’s
hidden state after all 20 timesteps have been processed.

We train on each of the four operating conditions in isolation, which aligns
with real-world applications. We train the model to minimize mean squared
error (MSE) of the RUL at the end of each 20-timestep window:

$J=\frac{1}{B}\sum_{B}\left\lVert RUL^{pred}_{t}-RUL^{gt}_{t}\right\rVert^{2}$
(8)

where $B$ is the batch size (64 in our case), $RUL^{pred}_{t}$ is the RUL
predicted by the model at timestep $t$ (end of the 20-timestep window), and
$RUL^{gt}_{t}$ is ground truth. During testing, we predict the RUL at the end
of each trajectory in the test set and report MSE across the entire test set
(same as Equation 8) along with the average L1 error. We report the average of
these metrics across 3 random seeds.

We use an Adam optimizer and a learning rate of 0.001. We train the model for
100 epochs on the FD001 and FD003 subsets and for 1000 epochs on the FD002 and
FD004 subsets. We also used early stopping, testing the best-performing model
from any epoch.

###  III-C Results

Table I presents the MSE and L1 error of the model under each operating
condition. The leftmost columns indicate the operating condition and whether
$\mu$ and/or $\rho$ were included as input during each training run.

Under all four operating conditions, including $\mu$ and $\rho$ significantly
outperforms the models where neither $\mu$ or $\rho$ are included. $\rho$,
however, appears to contain less information than $\mu$, as shown by our
$\rho$-only ablations; these runs nonetheless outperformed the baseline under
all operating conditions. FD002 and FD004 are clearly of much greater
complexity than FD001 and FD003, which is consistent with previous work [23].

Condition | Mu | Rho | Test MSE | Test L1
---|---|---|---|---
FD001 | ✗ | ✗ | 4.73 | 1.21
FD001 | ✓ | ✗ | 1.24 | 0.55
FD001 | ✗ | ✓ | 2.13 | 0.91
FD001 | ✓ | ✓ | 1.02 | 0.59
FD002 | ✗ | ✗ | 430.44 | 15.41
FD002 | ✓ | ✗ | 309.23 | 13.15
FD002 | ✗ | ✓ | 282.44 | 12.62
FD002 | ✓ | ✓ | 277.66 | 12.52
FD003 | ✗ | ✗ | 3.47 | 1.40
FD003 | ✓ | ✗ | 1.69 | 0.42
FD003 | ✗ | ✓ | 3.32 | 0.57
FD003 | ✓ | ✓ | 0.15 | 0.29
FD004 | ✗ | ✗ | 944.21 | 22.50
FD004 | ✓ | ✗ | 670.47 | 19.05
FD004 | ✗ | ✓ | 720.76 | 19.57
FD004 | ✓ | ✓ | 659.21 | 18.79
TABLE I: Test set mean squared error and L1 error under different operating
conditions (FD001, FD002, FD003, FD004), according to whether mu and rho were
added as inputs. Results are averaged across 3 seeds.

##  IV Summary and Conclusions

We have presented a physics informed machine learning (PIML) approach to
solving the problem of remaining useful lifetime (RUL) prediction using the
C-MAPSS dataset. Our approach, based on stochastic models of the underlying
data, is novel in that we do not use the traditional PIML architectures such
as PINNs, but rather augment the training data using estimated underlying
quantities. Finally, our apprach was shown to allow for generative synthesis
of data (i.e. physics aware synthetic data generation). The statistical
distribution of the sensor data was seen to exhibit multi-modal behavior, and
our approach is valid in this situation as well. Our experiments clearly
indicate the increased accuracy (measured both using MSE and absolute error)
over traditional deep learning (LSTM) approaches. It is expected that future
work in this area would involve increasing the scale and scope of physics
informed ML methods for predictive diagnostics.

## References

  * [1] Y. Huan, Q. Hu, Y. Jin, and G. Hu, “Fault detection and prognostics for avionics systems: A review,” _Aerospace Science and Technology_ , vol. 99, p. 106012, 2020.
  * [2] S. Singh, K. Iqbal, and R. Pathak, “A review on machine learning applications in aircraft maintenance,” _Aerospace Science and Technology_ , vol. 101, p. 106052, 2020.
  * [3] Z. Tian, S.-C. Tso, Y. Li, D. Xie, and K.-L. Tsui, “Health monitoring systems and prognostics for aircraft structures: A review,” _Aerospace Science and Technology_ , vol. 101, p. 106038, 2020.
  * [4] X. Tong, Y. Zhang, and X. Chen, “A review of prognostics and health management for aircraft systems,” _Aerospace Science and Technology_ , vol. 78, pp. 469–476, 2018.
  * [5] Z. Wang, D. Zhou, H. Jiang, and X. Shao, “Data-driven prognostics for aircraft components: A review,” _Aerospace Science and Technology_ , vol. 89, pp. 654–662, 2019.
  * [6] D. Yang, Y. Xu, and D. Wang, “Condition-based maintenance for aircraft engines based on big data analytics,” _Aerospace Science and Technology_ , vol. 94, pp. 74–82, 2019.
  * [7]
  * [8] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” _Journal of Computational Physics_ , vol. 378, pp. 686–707, 2019.
  * [9] J. Xu, H. Xiao, and G. E. Karniadakis, “Physics-informed neural networks: A comprehensive review,” _Journal of Computational Physics_ , vol. 404, p. 109180, 2020.
  * [10] W. Zhang, Z. Qi, J. Chen, L. Wang, and L. Tang, “Physics-informed machine learning: A survey on physics-constrained neural networks,” _arXiv preprint arXiv:2003.10811_ , 2020.
  * [11] Q. Yang, H. J. Kang, Y. Zhang, and D. N. Metaxas, “Physics-informed neural networks for inverse problems in imaging: A review,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 372, p. 113391, 2020.
  * [12] A. Ghavamian, E. Arabi, M. Gholami, and A. Borji, “Physics-informed machine learning for augmented reality flight simulation,” _Aerospace Science and Technology_ , vol. 101, p. 106042, 2020.
  * [13] J. Ling, A. Kurzawski, and J. Templeton, “Reynolds averaged turbulence modelling using deep neural networks with embedded invariance,” _Journal of Fluid Mechanics_ , vol. 807, pp. 155–166, 2016.
  * [14] J. Wang and P. Perdikaris, “Physics-informed deep learning for inverse problems in fluid dynamics,” _arXiv preprint arXiv:2001.04271_ , 2020.
  * [15] W. Chen, L. Zhang, and D. Lee, “Physics-informed machine learning for material science,” _Advanced Materials_ , vol. 25, no. 4, pp. 512–525, 2023.
  * [16] Y. Zhu, J. Qiao, and N. Zabaras, “Physics-informed deep learning for structural health monitoring and non-destructive evaluation: A review,” _Mechanical Systems and Signal Processing_ , vol. 148, p. 107194, 2021.
  * [17] J. Hamilton, I. Frosio, S. Tyree, and M. Prasad, “Learning physics from medical imaging using self-supervised deep neural networks,” _Nature Machine Intelligence_ , vol. 3, pp. 780–791, 2021.
  * [18] J. Onieva, F. Santini, S. Nikolenko, and S. Bakas, “Physics-informed deep learning for dynamic contrast-enhanced mri reconstruction,” _arXiv preprint arXiv:2002.05700_ , 2020.
  * [19] Y. Zhu and N. Zabaras, “Physics-informed neural networks for high-dimensional surrogate modeling and uncertainty quantification without labeled data,” _Journal of Computational Physics_ , vol. 403, p. 109151, 2020.
  * [20] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis, “Deepxde: A deep learning library for solving differential equations,” _arXiv preprint arXiv:1910.07291_ , 2019.
  * [21] A. Saxena, J. R. Celaya, B. Saha, and K. Goebel, “Damage propagation modeling for aircraft engine run-to-failure simulation,” _Aerospace Science and Technology_ , vol. 12, no. 4, pp. 277–288, 2008.
  * [22] X.-S. Si, W. Wang, and C. Hu, “Prognostics of equipment health: A review of data-driven methods and hybrid approaches with implementation guidelines,” in _2011 IEEE International Conference on Prognostics and Health Management_. IEEE, 2011, pp. 1–10.
  * [23] S. Zheng, K. Ristovski, A. Farahat, and C. Gupta, “Long short-term memory network for remaining useful life estimation,” in _2017 IEEE International Conference on Prognostics and Health Management (ICPHM)_ , 2017, pp. 88–95.
  * [24] X. Jiang, Y. Yang, F. Dong, H. Jiang, and Y. Miao, “Hybrid lstm–rnn model with adversarial training for remaining useful life prediction,” _Journal of Mechanical Science and Technology_ , vol. 34, no. 4, pp. 1491–1502, 2020.
  * [25] Y. Yang, F. Dong, Y. Miao, J. Yan, and J. Lee, “Deep learning approach for remaining useful life prediction based on lstm networks,” in _2018 Prognostics and System Health Management Conference (PHM-Asia)_. IEEE, 2018, pp. 1–6.
  * [26] K. Zhou, Z. Zhao, K. Mao, B. Hu, and Y. Zhang, “Remaining useful life estimation for turbofan engines based on a convolutional long short-term memory recurrent neural network,” _Journal of Mechanical Science and Technology_ , vol. 32, no. 9, pp. 4287–4297, 2018.