* Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schaye et al. (2023) Schaye J., et al., 2023, MNRAS,
  * Schroetter et al. (2016) Schroetter I., et al., 2016, ApJ, 833, 39
  * Smith et al. (2008) Smith B., Sigurdsson S., Abel T., 2008, MNRAS, 385, 1443
  * Smith et al. (2017) Smith B. D., et al., 2017, MNRAS, 466, 2217
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, Annual Review of Astronomy and Astrophysics, 53, 51
  * Sorini et al. (2022) Sorini D., Davé R., Cui W., Appleby S., 2022, MNRAS, 516, 883
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
  * Springel et al. (2005) Springel V., Di Matteo T., Hernquist L., 2005, ApJ, 620, L79
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Stinson et al. (2006) Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074
  * Stone et al. (2016) Stone M., Veilleux S., Meléndez M., Sturm E., Graciá-Carpio J., González-Alfonso E., 2016, ApJ, 826, 111
  * Sun et al. (2009) Sun M., Voit G. M., Donahue M., Jones C., Forman W., Vikhlinin A., 2009, ApJ, 693, 1142
  * Thomas et al. (2019) Thomas N., Davé R., Anglés-Alcázar D., Jarvis M., 2019, MNRAS, 487, 5764
  * Thorne et al. (2021) Thorne J. E., et al., 2021, MNRAS, 505, 540
  * Tillman et al. (2023a) Tillman M. T., et al., 2023a, AJ, 166, 228
  * Tillman et al. (2023b) Tillman M. T., Burkhart B., Tonnesen S., Bird S., Bryan G. L., Anglés-Alcázar D., Davé R., Genel S., 2023b, ApJ, 945, L17
  * Tollet et al. (2019) Tollet É., Cattaneo A., Macciò A. V., Dutton A. A., Kang X., 2019, MNRAS, 485, 2511
  * Tumlinson et al. (2013) Tumlinson J., et al., 2013, ApJ, 777, 59
  * Tumlinson et al. (2017) Tumlinson J., Peeples M. S., Werk J. K., 2017, ARA&A, 55, 389
  * Veilleux et al. (2005) Veilleux S., Cecil G., Bland-Hawthorn J., 2005, ARA&A, 43, 769
  * Veilleux et al. (2020) Veilleux S., Maiolino R., Bolatto A. D., Aalto S., 2020, A&ARv, 28, 2
  * Vikhlinin et al. (2006) Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Van Speybroeck L., 2006, ApJ, 640, 691
  * Villaescusa-Navarro et al. (2021) Villaescusa-Navarro F., et al., 2021, ApJ, 915, 71
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2014) Vogelsberger M., et al., 2014, MNRAS, 444, 1518
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Wadsley et al. (2004) Wadsley J. W., Stadel J., Quinn T., 2004, New Astron., 9, 137
  * Wang et al. (2015) Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83
  * Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
  * Werk et al. (2014) Werk J. K., et al., 2014, ApJ, 792, 8
  * Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Theuns T., Dalla Vecchia C., Tornatore L., 2009, MNRAS, 399, 574
  * Wright et al. (2020) Wright R. J., Lagos C. d. P., Power C., Mitchell P. D., 2020, MNRAS, 498, 1668
  * Wright et al. (2021) Wright R. J., Lagos C. d. P., Power C., Correa C. A., 2021, MNRAS, 504, 5702
  * Yang et al. (2022) Yang T., Cai Y.-C., Cui W., Davé R., Peacock J. A., Sorini D., 2022, MNRAS, 516, 4084
  * Zabl et al. (2019) Zabl J., et al., 2019, MNRAS, 485, 1961
  * Zahedy et al. (2019) Zahedy F. S., Chen H.-W., Johnson S. D., Pierce R. M., Rauch M., Huang Y.-H., Weiner B. J., Gauthier J.-R., 2019, MNRAS, 484, 2257
  * van Rossum & Drake Jr (1995) van Rossum G., Drake Jr F. L., 1995, Python reference manual. Centrum voor Wiskunde en Informatica Amsterdam
  * van de Voort (2017) van de Voort F., 2017, in Fox A., Davé R., eds, Astrophysics and Space Science Library Vol. 430, Gas Accretion onto Galaxies. p. 301, doi:10.1007/978-3-319-52512-9_13

##  Appendix A Comparison of gas flow measurements with previous literature

In this section, we compare our measurements of gas flow rates with those
measured from the same simulations (where applicable) in previous literature.
As outlined in §2.2, we elect to use an Eulerian or “instantaneous” approach
to calculating gas flow rates, rather than a Lagrangian particle tracking
method.

Figure 9 shows a comparison of $z\approx 2$ and $z\approx 0$ halo-scale mass
loading factors at the ISM scale (here, $0.20\times R_{\rm 200c}$) in EAGLE as
a function of $M_{\rm 200c}$ mass. While Mitchell et al. (2020) take a
Lagrangian particle-tracking approach to measure gas flow rates, we recover
good agreement with their results (to within $\approx 0.2\,$dex across the
halo mass range). Our instantaneously measured mass loading factors are
systematically (very) slightly higher than the Mitchell et al. (2020) results,
likely due to the latter not including the subset of gas that is very quickly
recycled within the $\Delta t$ used for the Lagrangian calculation. This
offset is highest at $M_{\rm 200c}\approx 10^{13}{\rm M}_{\odot}$, indicating
that some AGN-driven outflows in EAGLE may be recycled quite quickly at the
ISM scale.

Figure 10 shows a comparison of $z\approx 2$ $10{\rm kpc}$-scale mass loading
factors as a function of $M_{\star}$ in TNG for two different choices of
outflow velocity cut: $50\,{\rm km\,s^{-1}}$, and $150\,{\rm km\,s^{-1}}$.
Utilising a very similar method, our measurements agree very well, to within
$\approx 0.1\,{\rm dex}$ for both velocity thresholds, across the range of
stellar masses shown. As discussed in §2.2, in the main body of this paper we
do not enforce a minimum velocity threshold for the outflows. With
experimentation, we find that imposing a velocity threshold does not
qualitatively influence our results.

Figure 9: Comparison of $z\approx 2$ (blue) and $z\approx 0$ (black) halo-
scale mass loading factors ($\eta\equiv\langle\dot{M}_{\rm
out}\rangle/\langle\dot{M}_{\star}\rangle$) as a function of $M_{\rm 200c}$
mass measured in this work (solid lines, no shading), and the measurements
presented in Mitchell et al. (2020) (square markers connected with dot-filled
solid lines). In our calculations, hatched shaded regions correspond to the
bootstrap-generated $90\%$ confidence interval on the medians at a given mass.
At both redshifts, our measurements agree well, to within $\approx 0.2\,{\rm
dex}$. Figure 10: Comparison of $z\approx 2$ mass loading factors
($\eta\equiv\langle\dot{M}_{\rm out}\rangle/\langle\dot{M}_{\star}\rangle$) at
$10\,{\rm kpc}$ measured in this work (solid lines, no shading) and those
presented in Nelson et al. (2019) (square markers connected with dot-filled
solid lines). The data presented in orange correspond to outflow measurements
with a minimum outflow velocity of $150\,{\rm km\,s^{-1}}$, while the data
presented in blue correspond to outflow measurements with a minimum outflow
velocity of $50\,{\rm km\,s^{-1}}$. In our calculations, hatched shaded
regions correspond to the bootstrap-generated $90\%$ confidence interval on
the medians at a given mass. For both velocity thresholds, our measurements
agree very well, to within $\approx 0.1\,{\rm dex}$.

##  Appendix B Galaxy density profiles

In Fig. 11 and 12, we show the density profiles of galaxies in each of the
simulations at $z=2$ and $z=0$ respectively. Each row represents a different
bin in halo mass – top: $\log_{10}(M_{\rm 200c}/{\rm M}_{\odot})$
$\in[10.75,11.25]$; middle: $\log_{10}(M_{\rm 200c}/{\rm M}_{\odot})$
$\in[11.75,12.25]$; bottom: $\log_{10}(M_{\rm 200c}/{\rm M}_{\odot})$
$\in[12.75,13.25]$). Radial $x$-values are quoted as fractions of $R_{\rm
200c}$, with a broken scale above $1\times R_{\rm 200c}$ to include radial
bins up to $3\times R_{\rm 200c}$. For reference, in each panel we include a
dashed line indicating the expected outer slope of an isothermal $\beta$
profile with $\beta=2/3$; i.e. $\rho(r)\propto r^{-2}$.

In general, density profiles at $z=2$ agree fairly well between the
simulations. This corresponds to the universally higher baryon content of
galaxies and their CGM at this redshift, as demonstrated in Fig. 15.
Considering the lowest mass bin at $z=0$, we note that the density of the CGM
in EAGLE is quite low – $\approx 1\,{\rm dex}$ lower than the predictions from
TNG and SIMBA at the edge of the ISM. This low density allows for the
entrainment of gas with outflows in EAGLE, with the outflows being
overpressurised relative to the surrounding medium. We observe a similar
effect in SIMBA in the two higher mass bins, with CGM densities in the
intermediate halo mass bin being $\approx 1\,{\rm dex}$ lower than the
predictions from EAGLE and TNG. This signifies the efficient evacuation of the
CGM via AGN feedback, which also allows subsequent outflows to be
overpressurised relative to the CGM and entrain gas in the outflows towards
larger scales.

Figure 11: Radial density profiles of galaxies at $z=2$ in EAGLE (blue),
EAGLE-RECAL (teal), TNG (purple) and SIMBA (pink).Each column represents a
different bin in halo mass, increasing top to bottom. Error-bars correspond to
the $16^{\rm th}-84^{\rm th}$ percentile range in inflow rates at a given
mass, and hatched regions correspond to the bootstrap-generated $90\%$
confidence interval on the medians at a given mass. Gas density profiles are
relatively similar between the simulations at this redshift, corresponding to
the universally higher halo baryon fractions presented in Fig. 15.  Figure 12:
Radial density profiles of galaxies at $z=0$ in EAGLE (blue), EAGLE-RECAL
(teal), TNG (purple) and SIMBA (pink). Each row represents a different bin in
halo mass, increasing top to bottom. Error-bars correspond to the $16^{\rm
th}-84^{\rm th}$ percentile range in inflow rates at a given mass, and hatched
regions correspond to the bootstrap-generated $90\%$ confidence interval on
the medians at a given mass. Density profiles diverge between the simulations
at this redshift, with particularly low densities in EAGLE in the lowest halo
mass bin, and SIMBA for the two higher halo mass bins.

##  Appendix C Galaxy temperature profiles

In Fig. 13 and 14, we show the temperature profiles of galaxies in each of the
simulations at $z=2$ and $z=0$ respectively. Each row represents a different
bin in halo mass – top: $\log_{10}(M_{\rm 200c}/{\rm M}_{\odot})$
$\in[10.75,11.25]$; middle: $\log_{10}(M_{\rm 200c}/{\rm M}_{\odot})$
$\in[11.75,12.25]$; bottom: $\log_{10}(M_{\rm 200c}/{\rm M}_{\odot})$
$\in[12.75,13.25]$). Radial $x$-values are quoted as fractions of $R_{\rm
200c}$, with a broken scale above $1\times R_{\rm 200c}$ to include radial
bins up to $3\times R_{\rm 200c}$. For reference, in each panel we include a
horizontal line indicating the expected virial temperature of haloes in this
mass range, with $T_{\rm vir}\approx 1.1\times 10^{6}\,{\rm
K}\,\left(\mu/0.59\right)\left({M_{\rm halo}}/{10^{12}{\rm
M}_{\odot}}\right)^{2/3}(1+z)$, from van de Voort (2017).

As discussed in §4.1, considering haloes with mass $M_{\rm 200c}\lesssim
10^{13}{\rm M}_{\odot}$ at $z\approx 2$, predictions for temperatures in the
inner CGM are quite different between the simulations. In particular, in EAGLE
(and to some extent, SIMBA) there remains a significant presence of cool gas
between $0.3-0.6\times R_{\rm 200c}$. Comparatively, in TNG, temperatures in
this region typically exceed $10^{5}\,{\rm K}$. At lower redshift, the
predictions between the simulations are relatively similar in the two lower
mass bins (signifying a universal shift to a hotter CGM, and the dominance of
hot-mode accretion).

Figure 13: Radial temperature profiles of galaxies at $z=$ in EAGLE (blue),
EAGLE-RECAL (teal), TNG (purple) and SIMBA (pink). Each row represents a
different bin in halo mass, increasing top to bottom. Error-bars correspond to
the $16^{\rm th}-84^{\rm th}$ percentile range in inflow rates at a given
mass, and hatched regions correspond to the bootstrap-generated $90\%$
confidence interval on the medians at a given mass. Temperature predictions
diverge between the simulations at this redshift in the two lower mass bins,
particularly in the inner CGM. Figure 14: Radial temperature profiles of
galaxies at $z=0$ in EAGLE (blue), EAGLE-RECAL (teal), TNG (purple) and SIMBA
(pink). Each row represents a different bin in halo mass, increasing top to
bottom. Error-bars correspond to the $16^{\rm th}-84^{\rm th}$ percentile
range in inflow rates at a given mass, and hatched regions correspond to the
bootstrap-generated $90\%$ confidence interval on the medians at a given mass.
Temperature profiles are in slightly better agreement between the simulations
compared to $z\approx 2$, with generally hotter CGM.

##  Appendix D Halo baryon fractions at redshift 2

In this section, we add to the $z\approx 0$ results shown in §3 to compare the
total baryon content of haloes at as a function of halo mass, and how this
baryon content is distributed between stars, cool gas, and hot gas in each of
the simulations. As in the main body of this paper, in each simulation, we
define “cool gas” as the gas that is either considered star-forming, or below
$5\times 10^{4}{\rm K}$, and “hot gas” as the remaining gas within $R_{\rm
200c}$ that does not meet this criteria.

At this redshift, the baryon content of haloes is universally higher in the
simulations – for $M_{\rm 200c}\gtrsim 10^{11.5}{\rm M}_{\odot}$, all
simulations predict a baryon fraction above $50\%$ of the cosmic mean. Similar
to $z\approx 0$, the baryon content of EAGLEgalaxies at low halo mass –
$M_{\rm 200c}\lesssim 10^{12}{\rm M}_{\odot}$ – is lower than the predictions
from SIMBA and TNG, indicating that the ejective and preventative influence of
stellar feedback is already at play by $z\approx 2$. The transition masses
${\mathcal{M}}^{{\rm T}_{1}}$ and ${\mathcal{M}}^{{\rm T}_{2}}$ are similar to
those reported at $z\approx 0$, with the exception of SIMBA, where the
${\mathcal{M}}^{{\rm T}_{1}}$ mass at $z\approx 2$ increases to $\approx
10^{12}{\rm M}_{\odot}$ (as opposed to $\approx 10^{11.3}{\rm M}_{\odot}$ at
$z\approx 0$, with ${\mathcal{M}}^{{\rm T}_{2}}$ remaining at $\approx
10^{12.8}{\rm M}_{\odot}$). For SIMBA, this indicates that AGN feedback
becomes more efficient at lower masses towards $z\approx 0$.

Figure 15: The baryon content of haloes at $z\approx 2$ in EAGLE (blue), TNG
(purple) and SIMBA (pink). Top panel, $1^{\rm st}$ column: The median total
baryon fraction within $R_{\rm 200c}$ (stars, gas, black holes) of haloes as a
function of $M_{\rm 200c}$ mass. Top panels, $2^{\rm nd}-4^{\rm th}$ columns:
Total baryon content for each simulation within $R_{\rm 200c}$ broken down
into stars (hatched bars), “cool” gas (coloured shaded bars), and “hot gas”
(grey bars) in each simulation. In the top panels, this baryon content is
normalised by $f_{\rm b}\times M_{\rm 200c}$. Bottom panels: The bottom row of
panels shows the same breakdown, but re-normalised by the actual baryonic mass
in each halo, $M_{\rm bar}(r<R_{\rm 200c})$. In each of the simulations in the
left-hand panel, error-bars represent the $16^{\rm th}-84{\rm th}$ percentiles
in halo baryon fractions for a given mass bin; the hatched regions represent
the $90\%$ confidence interval on the median from 100 bootstrap re-samplings
in each bin. The difference between the simulations in terms of $z\approx 2$
baryon fractions are qualitatively similar compared to that shown at $z\approx
0$ in Fig. 2, but the baryon content of haloes is universally higher.


# Ab-initio study of magneto-ionic mechanisms in ferromagnet/oxide multilayers

Adriano Di Pietro1,2, Rohit Pachat3, Liza Herrera Diez3, Johannes W. van der
Jagt4,5, Dafiné Ravelosona3,4, Gianfranco Durin1
1 Istituto Nazionale di Ricerca Metrologica, Strada delle Cacce 91, 10135,
Torino, Italy
2Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129, Torino, Italy
3Centre de Nanosciences et de Nanotechnologies, CNRS, Université Paris-Saclay,
91120 Palaiseau, France
4 Spin-ion Technologies, Boulevard Thomas Gobert 10, 91120, Paliseau, France
5 Université Paris Saclay, 3 rue Juliot Curie, 91190, Gif-sur-Yvette, France

###### Abstract

The application of gate voltages in heavy metal/ferromagnet/Oxide multilayer
stacks has been identified as one possible candidate to manipulate their
anisotropy at will. However, this method has proven to show a wide variety of
behaviours in terms of reversibility, depending on the nature of the
metal/oxide interface and its degree of oxidation. In order to shed light on
the microscopic mechanism governing the complex magneto-ionic behaviour in
$\text{Ta/CoFeB/}\text{HfO}_{2}$, we perform ab-initio simulations on various
setups comprising $\text{Fe/O, Fe/HfO}_{2}$ interfaces with different oxygen
atom interfacial geometries. After the determination of the more stable
interfacial configurations, we calculate the magnetic anisotropy energy on the
different unit cell configurations and formulate a possible mechanism that
well describes the recent experimental observations in
$\text{Ta/CoFeB/}\text{HfO}_{2}$.

Electric field effects

††preprint: AAPM/123-QED

##  I  Introduction

The ever-increasing demand for memory storage in the modern IT industry has
made the need for new energy efficient storage alternatives all the more
important. Voltage control of magnetic anisotropy (VCMA) [1] has gained
scientific interest as one of the prime candidates to develop ultra-low energy
memory storage devices [2, 3, 4, 5] and is usually studied in 2 different
variants: the first one aims at modifying magnetic properties of thin films by
pure charge accumulation/depletion effects induced by voltage application [6,
7, 4, 8]. The second variant makes use of voltage induced ionic motion in
heavy metal (HM)/ferromagnet (FM)/Oxide (Ox) thin film multilayers to
carefully tune the oxygen/ferromagnet chemical and electrostatic interaction,
enabling the control of magnetic anisotropy [9, 10, 3]. The main advantage of
ionic manipulation in comparison to pure charge accumulation/depletion
techniques is the non-volatility of the magnetization switching, while the
trade off is a more complex reversibility mechanism combining ion mobility and
chemical composition at the FM/Ox interface. The work of Fassatoui et al. [11]
clearly shows how the application of gate voltages causes a reversible
magnetization switching in Pt/Co/AlOx as well as in Pt/Co/TbOx. At the same
time, applying a gate voltage in Pt/Co/MgOx has the effect of irreversibly
pushing the anisotropy easy axis out-of-plane. This discrepancy has been
attributed to the result of the different character of the ionic mobility of
the oxides: TbOx and AlOx have the common property of being oxides with a
predominantly oxygen based ionic mobility [11, 12], while MgOx is known to
have Mg as the principal ionic carrier under the application of gate voltages
[13]. CoFeB/Oxide multilayer structures are of great technological interest as
they have shown promise for the design of non-volatile, high-density memory
storage devices thanks to their high tunneling mangeto-resistance (TMR), low
damping and perpendicular magnetic anisotropy (PMA) [14, 15, 16]. A recent
work from Pachat et. al [17] highlighted a more complex magneto-ionic behavior
in Ta/CoFeB/HfO2 multilayers. The application of a gate voltage to the as-
grown material with in-plane anisotropy (IPA) initially causes a non-volatile,
irreversible spin-reorientation transition (SRT) to a perpendicular anisotropy
state (PMA). Further application of the gate voltage causes the transition to
a fully reversible regime. This 2-step process is in contrast with the picture
presented in [11] because ionic mobility in HfO2 is attributed to oxygen [18].
To formulate a hypothesis on the mechanism governing these different
reversibility behaviors in Ta/CoFeB/HfO2, we perform ab-initio simulations
using density functional theory (DFT) to determine the structural and magnetic
properties of two FM/Oxide interfaces. The work is structured as follows: In
Sec. II, we provide the computational details of our simulations and a brief
overview of the theoretical framework used to describe magnetic anisotropy in
FM/Oxide interfaces. In Sec. III we proceed and analyze the structural
properties of two different FM/Oxide interfaces; each one displaying either
interstitial or frontal oxygen positioning (see Fig. 1 and Fig. 3-(c)). After
having determined the optimal oxygen configurations of these setups, we
proceed and compute the magnetic anisotropy energy of the Fe/HfO2 unit cells.
We also explore the role of ionic mobility in determining the magnetic
anisotropy properties of Fe/HfO2 interfaces and highlight how the energy costs
involved in ionic mobility is different depending on the site occupied. In
Sec. IV, we discuss the results and compare them to the experimental data [17]
and theoretical predictions [7, 19] in order to formulate a hypothesis for the
appearance of different magneto-ionic regimes in CoFeB/HfO2 multilayers.
Finally, in Sec. V we provide a summary of the findings and outline new
possible systems to analyze to further probe our hypothesis.

##  II Methods and computational details

###  II.1 Structural relaxations

We perform structural relaxation using density functional theory and applying
the full-potential linearized augmented plane wave (FLAPW) method [20], as
implemented in the FLEUR code [21]. In particular, we rely on the generalized
gradient approximation for the exchange-correlation potential as implemented
by the Perdew–Burke-Ernzerhof (PBE) functional [22]. Since the simulation of
amorphous systems such as CoFeB is beyond the capabilities of ab-initio
methods, we reduced the ferromagnetic component of the system to the Fe atoms
only. This approximation is justified on the basis of the composition of the
Co20Fe60B20/HfO2 stacks studied in the literature, which are iron rich [17].
Furthermore, this approximation for CoFeB in ab-initio simulations is commonly
used in the literature [23]. We designed 5 different unit cells comprising: an
Fe/O interface (structures (I) and (II) in Fig. 1) composed of 5 magnetic
layers (ML) of Fe sandwiched between 2 mono-atomic layers of oxygen on each
side, an Fe/HfO2 system (structures (III) and (IV) in Fig. 1) composed of 5 ML
of Fe atoms sandwiched between 2 ML of HfO2 on each side. Finally, to account
for oxygen coming from atmospheric interaction with the sample, we designed an
Fe/HfO2 system displaying both frontal and interstitial oxygen atoms [24] .
This system is composed of 5 ML of Fe atoms sandwiched between 2 ML of HfO2 on
each side with an additional O-layer located in the interstitial site of 2 Fe
atoms (Fig. 3-(c)). We refer to this kind of system as ”mixed interface”
throughout the work. The in-plane lattice constant of the system is fixed to
the value of BCC Fe, i.e $a_{Fe}=2.87$ Å. Before performing structural
relaxation, we selected an energy cutoff value for the plane wave basis of
$275.5$ eV and a k-point mesh of dimensions $10\times 10\times 1$. These
parameters allow us to obtain self consistent energies converged to at least
$0.009$ eV/Atom. To determine the optimal interfacial configurations of oxygen
atoms, structural relaxation in shape and volume is performed up to $0.001$
eV/Å.

###  II.2  Magnetic anisotropy energy

The ab-initio calculation of the magnetic anisotropy energy [25] is performed
according to the following procedure [26, 24, 6]: after having determined the
more stable interfacial geometries via the structural relaxation outlined in
II.1, spin orbit coupling (SOC) is included in the Kohn-Sham Hamiltonian and
the magnetic anisotropy energy (MAE) is computed by comparing sums of one-
electron energies (”Ev-sum” in the MAE plots in Figs. 2,3,4-(b)) via the
magnetic force theorem [27, 28]. This method is widely used and has been
validated for transition metal interfaces [29]. The additional broken
symmetries emerging from the inclusion of magnetic moments in the unit cell
requires an increase in the k-point mesh to $12\times 12\times 1$. Both with
and without spin-orbit coupling, the energy cutoff is kept at $275.5$ eV. With
these parameters, we are able to obtain self consistent energies converged to
at least $0.009$ eV/Atom.

###  II.3 Theoretical background

The origins of oxygen enabled anisotropy manipulation in FM/Oxide interfaces
have been discussed extensively in the literature [19, 24]. The underlying
theoretical frameworks have been developed by Bruno [30] and Van Der Laan
[31], which successfully linked the appearance of a finite orbital magnetic
moment anisotropy to an anisotropy contribution in the energy by means of
spin-orbit coupling. These theoretical frameworks predict a heavy dependence
of the anisotropy energy on the exact shape and hybridization of the 3d
orbitals. Oxygen atoms, when posed at the right distance, can have dramatic
effects in the disruption of the otherwise almost anisotropic magnetic moment
distribution of the 3d orbitals of transition metal ferromagnets [24]. In
particular, the 3d-electrons tend to hybridize very effectively with the
oxygen $2p_{z}$ orbitals. Therefore, depending on the relative position of Fe
and O, different 3d-orbitals are going to be more or less involved in the
hybridization. If, for instance, we imagine an oxygen atom sitting on top of
an Fe atom, the atomic orbitals that are more likely to hybridize have OOP
symmetry [32] ($3d_{z^{2}},3d_{xz},3d_{yz}$). This results in an under
compensation of orbital magnetic moment coming from electrons in orbitals with
IP symmetry ($3d_{x^{2}-y^{2}},3d_{xy}$), which therefore shift the anisotropy
easy axis out of the plane [24] according to the relation

$\Delta E_{SO}=\xi_{SO}\frac{\Delta\mu}{4\mu_{B}}$ (1)

where $\Delta E_{SO}$ represents the anisotropy energy , $\Delta\mu$ the
orbital moment anisotropy and $\xi_{SO}$ the material dependent spin-orbit
coupling constant. If on the other hand, the oxygen atom is sitting in the
same plane of the Fe atom, hybridization is going to involve orbitals with IP
symmetry {$3d_{x^{2}-y^{2}},3d_{xy}$} and the orbital magnetic moment under
compensation of the {$3d_{z^{2}},3d_{xz},3d_{yz}$} orbitals will turn the
trend and shift the easy axis in-plane. Despite the thinness of the FM layers,
one should always consider that the contributions to magnetic anisotropy are
not limited to the first layer, but in fact often involve the second layer and
possibly beyond (as in the case of Fe/MgO MTJs [6]). In addition, the
appearance of magnetic anisotropy is still constrained by the symmetries of
the crystal field that is coupled to the spin of the electrons via SOC. This
implies that different lattice geometries have different angular dependencies
of the anisotropy energy [30]. BCC structures (which are the ones we are
concerned with) are predicted to have the usual relation

$\Delta E_{SO}=a+b\sin^{2}{\theta}$ (2)

where $\theta$ is the spin quantization axis angle with respect to the
$\hat{z}-$axis and $a,b$ represent material dependent constants of anisotropy.
An in-plane magnetic anisotropy (IPA) correspond to the minimum of $\Delta
E_{SO}$ for $\theta=\frac{\pi}{2}$, while perpendicular magnetic anisotropy
(PMA) correspond to the minimum of $\Delta E_{SO}$ for $\theta=0$. This is the
fitting function that we are going to use in all our MAE calculations
presented in Section III. The advantage of the perturbational approach is that
it does not require the comparison of self-consistent energies for the
determination of the fitting parameters $a$ and $b$ in (2), but rather allows
for the use of the magnetic force theorem [27, 28] outlined in II.2.

Figure 1: Comparison of the different interfacial configurations in different
samples. Unit cells (I) and (II) represent Fe/O interfaces while unit cells
(III) and (IV) represent the Fe/HfO2 interfaces. The bar plots represent the
ground state energies and highlight the more stable configurations. The yellow
highlighted number represents the absolute value of the energy difference
between the 2 relaxed structures of each type.

##  III Results

###  III.1 Pure interfaces

As a first step, we perform structural relaxation on 6 different unit cells
where we suppose no interaction with atmospheric oxygen has taken effect (we
refer to these as ”pure interfaces”). As can be seen by the plotted energies
(Fig. 1), the optimal oxygen configurations change depending on the system
considered: a pure Fe/O interface (structures (I) and (II) of (Fig. 1))
favours oxygen atoms to be located in the interstitial site with respect to Fe
atoms. This preference appears to be reverted in the case of HfO2, where a
frontal positioning of oxygen atoms with respect to to the Fe atoms seems to
be strongly favoured. After having determined the optimal configurations for
the different unit cells, we proceed and include spin-orbit coupling to
compute their magnetic anisotropy energy. In this case we focus specifically
on the Fe/HfO2 interface. By observing Fig. 2 we notice how the more stable
frontally aligned oxygen setup, (structure (IV) of (Fig. 1)) displays IP
magnetic anisotropy (Fig. 2-(a)) while an interstitial oxygen configuration
(structure (III) of (Fig. 1)) yields PMA (Fig. 2-(d)). At this point, we
introduce ionic mobility on the interfacial oxygen atoms and recalculate the
magnetic anisotropy. We notice how shifting the frontal oxygen atom in order
to reduce the Fe-O distance by $\approx 2$ Å, we are able to achieve PMA (Fig.
2-(b)). As a side note, we remark how shifting the interstitial oxygen atom of
structure (III) in Fig. 1 has the effect of recovering IP magnetic anisotropy
(Fig. 2-(e)).

Figure 2: Magnetic anisotropy energy (MAE) comparison of the pure Fe/HfO2
interface (structures (III) and (IV)). The $\theta$ angle of the spin
quantization axis (SQA) from eq.(2) is shown in the inset of (a). MAE for the
ground state of the Fe/HfO2 unit cell with frontal Fe-O distance of (a)
$5.763$ Å (equilibrium) and (b) $3.763$ Å (shifted). (c) Fe/HfO2 unit cell
with frontal oxygen positioning. d(O-Fe) represents the distance of the
interstitial oxygen species from the Fe surface respectively (marked with the
dashed line). Ground state of the Fe/HfO2 unit cell with interstitial Fe-O
distance of (d) $5.338$ Å (equilibrium) and (e) $3.338$ Å (shifted). (f)
Fe/HfO2 unit cell with interstitial oxygen positioning. d(O-Fe) represents the
distance of the interstitial oxygen species from the Fe surface respectively
(marked with the dashed line). Figure 3: Magnetic anisotropy energy (MAE)
($\theta$ from eq.(2)) in the Fe/HfO2 setup shown with shifted oxygen species
at the interface. (a) MAE of the mixed surface in its ground state. (b) MAE of
the mixed surface with a frontal Fe-O distance of $3.334$ Å . (c) MAE of the
mixed surface in with a frontal Fe-O distance of $3.246$ Å and an interstitial
Fe-O distance of $-0.077$ Å. (d) MAE of the mixed surface in with a frontal
Fe-O distance of $4.346$ Å and an interstitial Fe-O distance of $-0.068$ Å.
(e) MAE of the mixed surface in with a frontal Fe-O distance of $3.334$ Å and
an interstitial Fe-O distance of $-0.068$ Å. (f) Side view of the Fe/HfO2 unit
cell with the mixed setup. The d${}_{\text{(O-Fe)}1}$ and
d${}_{\text{(O-Fe)}2}$ denote the distance of the frontal and interstitial
oxygen species from the Fe surface respectively (marked with the dashed line).
Figure 4: (a) Energetic cost of shifting interstitial and frontal oxygen atoms
in the mixed interface setup displayed in panel (f) of Fig. 3. The starting
positions are $d_{(O-Fe)1}=1.123$ Å for the interstitial oxygen atom and
$d_{(O-Fe)2}=4.436$ Å for the frontal oxygen atoms (i.e. the setup of panel
(a) in Fig. 3) (b) Eigenvalue sum of the system with SQA point OOP -
Eigenvalue sum of the system with SQA point IP at different frontal oxygen
positions $d_{(O-Fe)2}$ represented in panel (f) of Fig. 3. The starting
positions are $d_{(O-Fe)1}=-0.0077$ Å for the interstitial oxygen atom and
$d_{(O-Fe)2}=3.236$ Å for the frontal oxygen atoms (i.e the setup of panel (c)
in Fig. 3). Positive values on the y-axis indicate that the system has IP
magnetic anisotropy.

###  III.2 Mixed Fe/HfO2 interface

As mentioned in Sec. II, we believe the actual experimental setup is not
presenting a pure interface. To account for the presence of additional mobile
oxygen impurities, we designed an Fe/HfO2 unit cell displaying both
interstitial and frontal oxygen alignment at the interface (see Fig. 3-(d)).
By analyzing the effect on the magnetic anisotropy of different oxygen species
mobility, we expect to be able to paint a plausible picture for the
microscopic mechanism governing magneto-ionic regimes in an experimental, more
disordered scenario. As can be seen in Fig. 3-(a), the relaxed structure with
the both frontal and interstitial oxygen displays IP magnetic anisotropy. If
we shift the frontal oxygens by a nominal distance of $1$ Å, we notice no
effect on magnetic anisotropy (Fig. 3-(b)). If on the other hand, we shift the
interstitial oxygen atom by the same amount, we notice how the system displays
PMA (Fig. 3-(d)). Despite having checked the effect of oxygen shifts on
anisotropy manipulation in this mixed surface setup, we expect the ionic
mobility behavior of these two oxygen species to be different given their
different environment. By observing Fig. 4-(a), we can see that the energy
cost of a frontal shift is lower than the energy cost of an interstitial
shift, in accordance with our expectations. The energy costs are obtained by
comparing the energy of identical unit cells that differ only by the position
of the frontal/interstitial oxygen atom [33]. If we proceed and further reduce
the distance of frontal oxygens from the Fe surface, we notice how the PMA of
Fig. 3-(d) is destroyed (Fig. 3-(c)). From Fig. 3-(e) we can see how reducing
the frontal Fe-O distance after PMA is achieved restores IP magnetic
anisotropy. We also point out how the disruption of PMA can also be obtained
by pushing the interstitial oxygens deeper in the ferromagnetic layer (Fig.
3-(c)) [24]. From Fig. 4-(b) we can however observe that once PMA is destroyed
by this type of oxygen incorporation, it cannot be restored by moving the
frontal oxygen atoms closer to the surface. This observation is also in
agreement with the experiments [17], where samples that were exposed to a
negative gate voltage for long times (for a reference of the field direction,
see Fig. 5) did not display any reversibility of the voltage induced SRT,
which was attributed to over-oxidation.

##  IV Discussion

Considering the results presented above, we propose the following hypothesis
for the appearance of different magneto-ionic regimes in HfO${}_{2}/$CoFeB. As
can be seen from Fig. 2-(a) and Fig. 1, the pure Fe/HfO2 surface with
frontally aligned oxygens appears to be both the more stable structure while
at the same time displaying IP magnetic anisotropy. This is in accordance with
the experimental data, where samples in the as-grown form displayed IP
magnetic anisotropy [17]. We also point out how the ground state of the
Fe/HfO2 surface with both frontal and interstitial oxygens also displays IP
magnetic anisotropy (Fig. 3-(a)), leading us to hypothesize that the actual
as-grown surface displays these 2 interfacial oxygen arrangements (see Fig.
5-(a)). As observed in [17], the application of a gate negative voltage across
the Ta/CoFeB/HfO2 favouring the diffusion of oxygen species into the CoFeB
layer, results in an irreversible SRT from IPA to PMA. We hypothesize that in
this first regime the frontal oxygen atoms, which are easier to mobilize (Fig.
4-(a)), are too far away to contribute to the anisotropy manipulation (Fig.
3-(b)). In contrast, the (harder to mobilize) interstitial oxygen atoms of the
mixed surface have a decisive impact on the anisotropy of the system and can
induce PMA with very small displacements (Fig. 3-(d)). We therefore conclude
that the initial irreversible anisotropy change in Ta/CoFeB/HfO2 is due to the
irreversibility of the shifts of oxygens located at the interstitial site of
mixed interface configurations (Fig. 5-(b)). Once PMA is achieved (Fig.
3-(d)), we hypothesize that the interstitial oxygen atoms become essentially
immobile under the effect of the electric field because of the screening
effects of the conductive material [34]. The experimental data [17] shows that
the application of the gate voltage beyond the PMA state has the effect of
pushing the magnetic anisotropy easy axis in-plane, albeit in a reversible
way. Our results suggest that this switch to a reversible behavior beyond PMA
could be explained by the following: the continued application of the gate
voltage has the effect of mobilizing the remaining frontal oxygen atoms at the
pure and mixed sites, which move closer to the surface and restore in-plane
magnetic anisotropy (Fig. 3-(e)). Mobilizing frontal oxygen species is easier
(Fig. 4) leading us to believe that the reversibility of the regime presented
in ref.[17] could be attributed to the reversible shifts of the frontal oxygen
species of both the pure and the mixed interface (Fig. 5-(c)). As shown in
Fig. 3-(c), the transition PMA $\rightarrow$ IPA can also be achieved by
shifting the interstitial oxygens deeper inside the ferromagnetic layer. We do
not attribute the manipulation of anisotropy beyond initial PMA to these
oxygen species because, once the anisotropy is shifted in-plane by means of
the interstitial oxygen species pushed deeper in the sample, it is impossible
to manipulate the magnetic anisotropy of the system by shifting the frontal
oxygen species closer to the Fe surface (as displayed in Fig. 4-(b)).

Figure 5: Hypothesis for the mechanism governing the different magneto-ionic
regimes in CoFeB/HfO2 multilayers. (a) Shows the ground state of the system
with both pure and mixed interfaces. (b) Irreversible magnetization switching
via interstitial oxygen shifts. (c) Reversible magnetization switching via
frontal oxygen shifts. The bottom cartoon in all three panels represents the
magnetization direction and the switching process.

##  V Conclusion

In this work we performed an ab-initio analysis of the interplay between
oxygen ionic mobility and anisotropy manipulation in $2$ FM/Oxide interfaces
in order to formulate a hypothesis for the appearance of magneto-ionic regimes
in Ta/CoFeB/HfO2 stacks [17]. We found out that the different nature of the
oxide at the interface plays an important role in determining the optimal
interfacial oxygen geometry. In particular, we discovered how the pure Fe/HfO2
interface displays a preferential frontal oxygen alignment which corresponds
to an IP magnetic anisotropy. We observed how frontal oxygen mobility can
induce PMA [4]. The inclusion of oxygen species in additional interstitial
sites was investigated in order to determine their role in the appearance of
magneto-ionic regimes [17]. We have shown how in these so called mixed
surfaces, the ionic mobility of frontal oxygen species is energetically more
favorable than the mobility of interstitial ones. We have also shown how the
interplay of mobility between these 2 different oxygen species can change the
magnetic anisotropy of the sample. We conclude that the irreversibility of the
transition between IPA (under-oxidised) and PMA (optimally oxidised) is mainly
due to the interaction of the interstitial oxygen of the mixed surface and
that the reversibility of the second regime is due to the mobility of the
frontal oxygen species of the pure and mixed surface (Fig. 5). We point out
how the Ta/CoFeB/HfO2 multilayer in [17] is composed by amorphous materials
that do not display the ordered structure of our unit cells. In spite of this
difference, we find qualitative agreement between our results and the ones
reported in [17] and are therefore led to believe that anisotropy manipulation
in these kinds of systems may be a consequence of relative FM/oxygen
positioning localized at the interface. To further probe the validity of the
hypothesis, one could analyze the relative range of magneto-ionic regimes in
Ta/CoFeB/HfO2 samples with different degrees of atmospheric oxygen
interaction. We predict that in the extreme case where there has been little
to no exposure to atmospheric oxygen should result in a purely reversible
system. In addition, a comparison with crystalline and policrystalline
structures could prove useful in understanding the role played by the
amorphous nature of the materials. Understanding the precise nature of the
interplay between ionic mobility and magnetic property tuning could prove very
useful for the optimization of highly energy efficient read/write mechanisms
for next-generation memory storage devices.

##  VI Acknowledgement

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 860060 “Magnetism and the effect of Electric Field” (MagnEFi).
Computational resources were provided by HPC@POLITO, a project of Academic
Computing within the Department of Control and Computer Engineering at the
Politecnico di Torino (http://hpc.polito.it)

## References

  * One  _et al._ [2021] R.-A. One, H. Béa, S. Mican, M. Joldos, P. B. Veiga, B. Dieny, L. D. Buda-Prejbeanu, and C. Tiusan, Route towards efficient magnetization reversal driven by voltage control of magnetic anisotropy, Sci. Rep. 11, 8801 (2021).
  * Suwardy  _et al._ [2019] J. Suwardy, M. Goto, Y. Suzuki, and S. Miwa, Voltage-controlled magnetic anisotropy and Dzyaloshinskii-Moriya interactions in CoNi/MgO and CoNi/Pd/MgO, Jpn. J. Appl. Phys. 58, 60917 (2019).
  * Dieny and Chshiev [2017] B. Dieny and M. Chshiev, Perpendicular magnetic anisotropy at transition metal/oxide interfaces and applications, Rev. Mod. Phys. 89, 025008 (2017).
  * Hallal  _et al._ [2013] A. Hallal, H. X. Yang, B. Dieny, and M. Chshiev, Anatomy of perpendicular magnetic anisotropy in Fe/MgO magnetic tunnel junctions: First-principles insight, Phys. Rev. B 88, 1 (2013).
  * Belabbes  _et al._ [2016] A. Belabbes, G. Bihlmayer, S. Blügel, and A. Manchon, Oxygen-enabled control of Dzyaloshinskii-Moriya Interaction in ultra-thin magnetic films, Sci. Rep. 6, 24634 (2016).
  * Ibrahim  _et al._ [2016] F. Ibrahim, H. X. Yang, A. Hallal, B. Dieny, and M. Chshiev, Anatomy of electric field control of perpendicular magnetic anisotropy at Fe/MgO interfaces, Phys. Rev. B 93, 184423 (2016).
  * Huang  _et al._ [2013] Z. Huang, I. Stolichnov, A. Bernand-Mantel, J. Borrel, S. Auffret, G. Gaudin, O. Boulle, S. Pizzini, L. Ranno, L. H. Diez, and N. Setter, Ferroelectric control of magnetic domains in ultra-thin cobalt layers, Appl. Phys. Lett. 103, 222902 (2013).
  * Niranjan  _et al._ [2010] M. K. Niranjan, C. G. Duan, S. S. Jaswal, and E. Y. Tsymbal, Electric field effect on magnetization at the Fe/MgO(001) interface, Appl. Phys. Lett. 96, 222504 (2010).
  * Monso  _et al._ [2002] S. Monso, B. Rodmacq, S. Auffret, G. Casali, F. Fettar, B. Gilles, B. Dieny, and P. Boyer, Crossover from in-plane to perpendicular anisotropy in Pt/CoFe/AlO x sandwiches as a function of Al oxidation: A very accurate control of the oxidation of tunnel barriers, Appl. Phys. Lett. 80, 4157 (2002).
  * Bauer  _et al._ [2015] U. Bauer, L. Yao, A. J. Tan, P. Agrawal, S. Emori, H. L. Tuller, S. Van Dijken, and G. S. Beach, Magneto-ionic control of interfacial magnetism, Nat. Mater. 14, 174 (2015).
  * Fassatoui  _et al._ [2020] A. Fassatoui, J. P. Garcia, L. Ranno, J. Vogel, A. Bernand-Mantel, H. Béa, S. Pizzini, and S. Pizzini, Reversible and Irreversible Voltage Manipulation of Interfacial Magnetic Anisotropy in Pt / Co /Oxide Multilayers, Phys. Rev. Applied 14, 064041 (2020).
  * Flynn [1962] H. G. Flynn, Transport Phenomena, Phys. Today 15, 96 (1962).
  * Wu  _et al._ [2017] N. Wu, W. Wang, Y. Wei, and T. Li, Studies on the effect of nano-sized MgO in magnesium-ion conducting gel polymer electrolyte for rechargeable magnesium batteries, Energies 10, 1215 (2017).
  * Liu  _et al._ [2014] T. Liu, Y. Zhang, J. W. Cai, and H. Y. Pan, Thermally robust Mo/CoFeB/MgO trilayers with strong perpendicular magnetic anisotropy, Sci. Rep. 4, 5895 (2014).
  * Ikeda  _et al._ [2010] S. Ikeda, K. Miura, H. Yamamoto, K. Mizunuma, H. D. Gan, M. Endo, S. Kanai, J. Hayakawa, F. Matsukura, and H. Ohno, A perpendicular-anisotropy CoFeB-MgO magnetic tunnel junction, Nat. Mater. 9, 721 (2010).
  * Wang  _et al._ [2020] L. Z. Wang, X. Li, T. Sasaki, K. Wong, G. Q. Yu, S. Z. Peng, C. Zhao, T. Ohkubo, K. Hono, W. S. Zhao, and K. L. Wang, High voltage-controlled magnetic anisotropy and interface magnetoelectric effect in sputtered multilayers annealed at high temperatures, Science China: Physics, Mechanics and Astronomy 63, 277512 (2020).
  * Pachat  _et al._ [2021] R. Pachat, D. Ourdani, J. van der Jagt, M.-A. Syskaki, A. D. Pietro, Y. Roussigné, S. Ono, M. Gabor, M. Chérif, G. Durin, J. Langer, M. Belmeguenai, D. Ravelosona, and L. H. Diez, Multiple magnetoionic regimes in ta/co20fe60b20/HfO2, Phys. Rev. Applied 15, 064055 (2021).
  * Schie  _et al._ [2017] M. Schie, M. P. Müller, M. Salinga, R. Waser, and R. A. De Souza, Ion migration in crystalline and amorphous HfOX, J. Chem. Phys. 146, 094508 (2017).
  * Liang  _et al._ [2014] S. H. Liang, T. T. Zhang, P. Barate, J. Frougier, M. Vidal, P. Renucci, B. Xu, H. Jaffrès, J. M. George, X. Devaux, M. Hehn, X. Marie, S. Mangin, H. X. Yang, A. Hallal, M. Chshiev, T. Amand, H. F. Liu, D. P. Liu, X. F. Han, Z. G. Wang, and Y. Lu, Large and robust electrical spin injection into GaAs at zero magnetic field using an ultrathin CoFeB/MgO injector, Phys. Rev. B 90, 085310 (2014).
  * Andersen [1975] O. K. Andersen, Linear methods in band theory, Phys. Rev. B 12, 3060 (1975).
  * [21] Fleur code : www.flapw.de.
  * Perdew  _et al._ [1996] J. P. Perdew, M. Ernzerhof, and K. Burke, Rationale for mixing exact exchange with density functional approximations, J. Chem. Phys. 105, 9982 (1996).
  * Lin  _et al._ [2020] W. Lin, B. Yang, A. P. Chen, X. Wu, R. Guo, S. Chen, L. Liu, Q. Xie, X. Shu, Y. Hui, G. M. Chow, Y. Feng, G. Carlotti, S. Tacchi, H. Yang, and J. Chen, Perpendicular Magnetic Anisotropy and Dzyaloshinskii-Moriya Interaction at an Oxide/Ferromagnetic Metal Interface, Phys. Rev. Lett. 124, 217202 (2020).
  * Yang  _et al._ [2011] H. X. Yang, M. Chshiev, B. Dieny, J. H. Lee, A. Manchon, and K. H. Shin, First-principles investigation of the very large perpendicular magnetic anisotropy at Fe$|$MgO and Co$|$MgO interfaces, Phys. Rev. B 84, 054401 (2011).
  * Strange  _et al._ [1991] P. Strange, J. B. Staunton, B. L. Györffy, and H. Ebert, First principles theory of magnetocrystalline anisotropy, J. Phys.: Condens. Matter 1, 3947 (1991).
  * Zimmermann  _et al._ [2019] B. Zimmermann, G. Bihlmayer, M. Böttcher, M. Bouhassoune, S. Lounis, J. Sinova, S. Heinze, S. Blügel, and B. Dupé, Comparison of first-principles methods to extract magnetic parameters in ultrathin films: Co/Pt(111), Phys. Rev. B 99, 214426 (2019).
  * Weinert  _et al._ [1985] M. Weinert, R. E. Watson, and J. W. Davenport, Total-energy differences and eigenvalue sums, Phys. Rev. B 32, 2115 (1985).
  * Wang  _et al._ [1996] X. Wang, D. S. Wang, W. Ruqian, and A. J. Freeman, Validity of the force theorem for magnetocrystalline anisotropy, J. Magn. Magn. Mater. 159, 337 (1996).
  * Błoński and Hafner [2009] P. Błoński and J. Hafner, Density-functional theory of the magnetic anisotropy of nanostructures: An assessment of different approximations, J. Phys.: Condens. Matter 21, 426001 (2009).
  * Bruno [1989] P. Bruno, Tight-binding approach to the orbital magnetic moment and magnetocrystalline anisotropy of transition-metal monolayers, Phys. Rev. B 39, 865 (1989).
  * Van Der Laan [1998] G. Van Der Laan, Microscopic origin of magnetocrystalline anisotropy in transition metal thin films, J. Phys.: Condens. Matter 10, 3239 (1998).
  * Butler [2008] W. H. Butler, Tunneling magnetoresistance from a symmetry filtering effect, Sci. Technol. Adv. Mat. 9, 014106 (2008).
  * Kyritsakis  _et al._ [2019] A. Kyritsakis, E. Baibuz, V. Jansson, and F. Djurabekova, Atomistic behavior of metal surfaces under high electric fields, Phys. Rev. B 99, 205418 (2019).
  * Zhang [1999] S. Zhang, Spin-dependent surface screening in ferromagnets and magnetic tunnel junctions, Phys. Rev. Lett. 83, 640 (1999).


# Windscreen Optical Quality for AI Algorithms: Refractive Power and MTF not
Sufficient

1st Dominik Werner Wolf  Glass Laboratory
Volkswagen Group
Wolfsburg, Germany
<EMAIL_ADDRESS>2nd Markus Ulrich  Inst. of Photogram. and
Remote Sens.
Karlsruhe Institute of Technology
Karlsruhe, Germany
<EMAIL_ADDRESS>3rd Alexander Braun  Faculty of Electr. Eng. and Info.
Tech.
University of Applied Sciences Duesseldorf
Duesseldorf, Germany
<EMAIL_ADDRESS>
###### Abstract

Windscreen optical quality is an important aspect of any advanced driver
assistance system, and also for future autonomous driving, as today at least
some cameras of the sensor suite are situated behind the windscreen.
Automotive mass production processes require measurement systems that
characterize the optical quality of the windscreens in a meaningful way, which
for modern perception stacks implies meaningful for artificial intelligence
(AI) algorithms. The measured optical quality needs to be linked to the
performance of these algorithms, such that performance limits – and thus
production tolerance limits – can be defined. In this article we demonstrate
that the main metric established in the industry – refractive power – is
fundamentally not capable of capturing relevant optical properties of
windscreens. Further, as the industry is moving towards the modulation
transfer function (MTF) as an alternative, we mathematically show that this
metric cannot be used on windscreens alone, but that the windscreen forms a
novel optical system together with the optics of the camera system. Hence, the
required goal of a qualification system that is installed at the windscreen
supplier and independently measures the optical quality cannot be achieved
using MTF. We propose a novel concept to determine the optical quality of
windscreens and to use simulation to link this optical quality to the
performance of AI algorithms, which can hopefully lead to novel inspection
systems.

###### Index Terms:

windscreen optical quality, AI algorithms, computer vision, refractive power,
MTF

##  I Introduction

Every car has a windscreen. The number of newly produced windscreens therefore
ranges in the millions every year. Following quality processes for automotive
mass production established since the 1960ies – like the outdated ISO/TS 16949
[18] or the more recent VDA6.3 [31] – these windscreens are tested end-of-line
(EOL) at the suppliers (Tier 1) production line using well-defined optical
measurements. Importantly, the windscreen quality is measured at the
production site alone, independent of any production tolerances that may arise
during assembly of the whole car at the site of the car manufacturer (original
equipment manufacturer, OEM). Economically, this is mandatory, as a thorough
testing of the whole windscreen after assembly by the OEM is prohibitively
expensive.

For several decades the optical quality of these windscreens has been judged
acceptable if humans could look through it with low impact on the perception
of the driver. With the rise of advanced driver assistance systems (ADAS) and
the future promise of autonomous driving (AD) many cars nowadays are equipped
with several camera systems, many of which are situated behind the windscreen.
A camera is not a human observer, and it is now not enough to qualify a
windscreen using human perception, especially as the quality and resolution of
the cameras are steadily increasing. The influence of the optical quality on
the image quality and further on the computer vision algorithms evaluating
these images has to be precisely determined.

In theory, the working limits of the computer vision algorithms are
determined, and production tolerance limits are derived from these algorithmic
working limits through a number of processes defined in the above mentioned
quality norms. Opto-mechanical tolerance calculations, numerical simulations
and test campaigns in the real world form three important pillars of these
studies accompanied by environmental stress tests and aging simulations [15,
7].

In practice, though, modern camera-based ADAS applications are based on
artificial intelligence (AI) and employ deep convolutional neural networks,
due to their superior performance in comparison to traditional, rule-based
computer vision algorithms. The difference in performance is such that
currently there is no alternative to using AI algorithms. As these AI
algorithms are ’black boxes’ in nature [13], i.e. the output cannot be
predicted, the link between optical quality and AI algorithm performance
cannot be easily established [25]. And due to the lack of quantitative working
limits w.r.t. the AI algorithms, production tolerance limits for the
windscreens can not be straightforwardly deduced [4].

In this article we are evaluating the two main measurement processes that are
currently used in the automotive industry to qualify windscreen optical
quality: refractive power and the modulation transfer function (MTF). While
refractive power is the established measurement method and has been
standardized already in the 1990ies [6, 9], the MTF – or equivalently the
spatial frequency response (SFR) – has gained recent attention as automotive
researchers [32] look for alternatives to refractive power because of the
increasing ADAS camera performances in terms of the number of pixels per field
angle. Novel startups are even forming around the promise of using MTF to
characterize windscreens.

We find and mathematically demonstrate in this work that both refractive power
and MTF are not sufficient to quantify windscreen quality for AI algorithm
performance. This is a fundamental finding in that our results are derived
from first principles of optics, and apply very generally. First, we
recapitulate the optical basics in Sec. II. Importantly, the optical quality
is described in terms of wavefront aberrations, using the Zernike formalism to
mathematically decompose the nature of the optical perturbations. Then in Sec.
III, using these basics we show how refractive power is fundamentally not
capable of accounting for a distinct number of wavefront aberrations, while at
the same time these aberrations have a demonstrable effect on AI algorithm
performance [26, 24].

In Sec. IV we then show how the windscreen and the camera system form a joint
optical system, that – again fundamentally – cannot be separated into two
distinct optical systems. This separation, though, is a necessary requirement
in linear system theory for the multiplicativity of the system MTF w.r.t. the
individual optical elements [12]. Therefore, this prohibits any MTF
measurement on the windscreen alone, and thus from using MTF as a qualifying
measurement at the production site of the Tier 1.

Optical quality has many different aspects. For this article, we concentrate
solely on the ’sharpness’ of the camera image, which is deteriorated by
optical path variations across the windshield plane and is typically
quantified by the MTF in optical linear system theory. In general, lens
distortions, which describe the failure of a lens to map lines into lines and
represent a curvilinear mapping [28], might also deteriorate the performance
of ADAS functionalities. Effects of optical distortions will not be considered
in the following.

In summary, we will show how the two only current measurement techniques in
the automotive industry are not sufficient to measure the sharpness of the
windscreen alone. These results have far reaching implications for the
automotive industry, which needs to focus more effort on finding alternatives.
We finally propose a concept on how to find a novel measurement process,
combining optical modeling, numerical simulation and AI algorithms to link the
optical quality of windscreens to the performance of AI algorithms.

##  II Optical Quality and Mathematical Models

Maxwell’s equations are the fundamental physical model of electromagnetic
radiation, and the wave equation forms the basis for the technological
application of light. If all elements in the optical system are large compared
to the wavelength of the light, geometrical optics may be used. It plays an
important role in the development of optical systems as well, in the form of
raytracing simulations. A windscreen is large in mechanical dimensions, both
laterally as well as axially, but previous work has shown that the aberrations
originating inside the windscreen cannot be neglected [20, 5]. Thus, it is not
sufficient to take only the geometry of the windscreen into account – which
would allow for a raytracing approach – but a comprehensive optical model
needs to be based on the wave description of light. This is why in the
following we use the fundamental Zernike approach [10] to model wavefront
aberrations, where the optical path difference mathematically models the
aberrations present in the windscreen.

###  II-A Wavefront Modelling with Zernike Polynomials

The optical path difference $W$, defined on the principle plane, is usually
expressed as a decomposition into Zernike polynomials $Z_{n}$ with
corresponding Zernike coefficients $c_{n}$ (in units of meters) as [3]:

$W(\rho,\;\phi)=\sum\limits_{n=0}^{\infty}c_{n}Z_{n}(\rho,\;\phi)\;\;\;\text{,}\;\;\;c_{n}:\overset{(\ref{eq:Zernike_scalar_product})}{=}\left<W,\;Z_{n}\right>\;.$
(1)

Here, the domain of the principle plane of the optical element is
parameterized by normalized polar coordinates with radius $\rho$ and polar
angle $\phi$. There are different numbering schemes for Zernike polynomials,
i.a. a linear numbering scheme according to the American National Standards
Institute (ANSI) which has been adopted within this work. The Zernike
polynomials reproduce the aberration pattern on the unit circle and correspond
to different, independent optical perturbations like defocus or astigmatism.

The independence of the perturbations is mathematically reflected by the
orthogonality relation of the scalar product:

$\left<Z_{i},\;Z_{j}\right>\coloneqq\int\limits_{0}^{2\pi}\int\limits_{0}^{1}Z_{i}(\rho,\;\phi)\cdot
Z_{j}(\rho,\;\phi)\cdot\rho\;\mathrm{d}\rho\;\mathrm{d}\phi=\pi\cdot\delta_{ij}\;.$
(2)

This is important, because we will demonstrate how certain Zernike polynomials
are simply not present in the refractive power measurement, and the
orthogonality fundamentally implies that this information can not be
recovered. Table I indicates the normalized Zernike polynomials defined by ISO
24157 [17] up to the third order.

$Z_{i}$ | Zernike polynomial | Harmonic
---|---|---
| Polar coordinates | Cartesian coordinates |
$Z_{0}$ | $1$ | $1$ | ✓
$Z_{1}$ | $2\rho\sin\phi$ | $2y$ | ✓
$Z_{2}$ | $2\rho\cos\phi$ | $2x$ | ✓
$Z_{3}$ | $\sqrt{6}\rho^{2}\sin 2\phi$ | $2\sqrt{6}xy$ | ✓
$Z_{4}$ | $\sqrt{3}(2\rho^{2}-1)$ | $\sqrt{3}(2x^{2}+2y^{2}-1)$ | $\mathbf{\times}$
$Z_{5}$ | $\sqrt{6}\rho^{2}\cos 2\phi$ | $\sqrt{6}(x^{2}-y^{2})$ | ✓
$Z_{6}$ | $\sqrt{8}\rho^{3}\sin 3\phi$ | $\sqrt{8}(3x^{2}y-y^{3})$ | ✓
$Z_{7}$ | $\sqrt{8}(3\rho^{3}-2\rho)\sin\phi$ | $\sqrt{8}(3x^{2}y+3y^{3}-2y)$ | $\mathbf{\times}$
$Z_{8}$ | $\sqrt{8}(3\rho^{3}-2\rho)\cos\phi$ | $\sqrt{8}(3x^{3}+3xy^{2}-2x)$ | $\mathbf{\times}$
$Z_{9}$ | $\sqrt{8}\rho^{3}\cos 3\phi$ | $\sqrt{8}(x^{3}+3xy^{2})$ | ✓

TABLE I: Zernike polynomials up to the third order.

###  II-B Refractive Power

Refractive power measures how much focusing power a lens has. It is given in
units of diopters, i.e. in inverse distance of the focal length of the lens. A
comprehensible way to visualize refractive power is two parallel light rays
entering the optical system – here: the windscreen – and upon exit are not
parallel anymore, but either divergent or convergent. In the convergent case,
the focal length is the distance from the refractive element to the
intersection of the two rays, and its inverse is the numerical value of the
refractive power. For concave lenses, the diverging rays are extended in the
negative direction until these two rays intersect, and the negative distance
now forms the focal length.

Figure 1: Refractive power measurement of the ADAS camera area of a VW series
production windshield under an inclination angle of $\epsilon=63$°. The
difference in magnitude between horizontal and vertical direction originates
from the inclination angle, which amplifies the refractive power according to
the Kerkhof model [33].

For windscreens, the refractive power is not a single number for the whole
glass, but the measurement has a spatial resolution as depicted by Fig. 1. In
the early days, two actual parallel laser beams were deflected, and the whole
setup was laterally moved to achieve a certain spatial resolution [6]. More
modern systems such as the one produced by ISRA use the Moiré effect to
spatially resolve the refractive power over a limited area by observing the
location dependency of the perturbed grid spacing between Moiré interferences
[23]. In addition, new refractive power measurement systems like the one
produced by LaVision [22] use the Background Oriented Schlieren (BOS) imaging
method to overcome the resolution limitation of the Moiré approach [33].

Importantly, the refractive power depends on the direction, as the two
parallel lines form a plane together with the principal plane of the optical
element. In principle, this direction can be rotated full circle by
$360\text{\,}\mathrm{\SIUnitSymbolDegree}$, but in practice, the refractive
power is determined and specified only in the horizontal and vertical
direction.

###  II-C Modulation Transfer Function

The modulation transfer function (MTF) – and its non-harmonic equivalent, the
spatial frequency response (SFR) – are established metrics to characterize
optical systems, based on linear system theory and scalar diffraction theory
[2, 10]. In image space, the transfer function of the system under test is
called the point spread function (PSF), and in frequency space, it is denoted
as the optical transfer function (OTF). The MTF is given by the absolute value
of the OTF and is of particular importance if the intensity distribution is
the matter of interest. The PSF and the MTF are highly non-linear functions
over the image field (radius, azimuth), and they also depend on the defocus
$\Delta z$, due to the refractions on different lens element surfaces. Hence,
the input space of the PSF is in general three dimensional.

The MTF is measured by using either harmonic input signals (MTF, e.g.
sinusoidal Siemens star) or a step function type input (SFR, e.g. slanted
edge). ISO12233 defines a norm to measure the MTF [16], and IEEE P2020 is
currently finalizing an automotive extension of this norm [14]. In this
article, we will use slanted edge measurements.

According to scalar diffraction theory, the MTF is proportional to the
absolute value of the Fourier transform of the wavefront in the aperture plane
of the lens (more general: the optical element). The wavefront is transformed,
normalized, and the absolute value is taken to yield the MTF. This allows for
an analytical relationship between the MTF and the wavefront aberrations,
which can be parameterized by Zernike coefficients $c_{n}$:

$\mathrm{MTF}(\vec{k}|\lambda)=\left|\dfrac{\leavevmode\resizebox{65.04034pt}{}{$\oiint\limits_{P_{+}\cap
P_{-}}$}\exp{\left(\dfrac{2\pi\mathrm{i}}{\lambda}\sum\limits_{n=0}^{\infty}c_{n}\left[Z_{n}(\vec{\xi}+\vec{\Delta})-Z_{n}(\vec{\xi}-\vec{\Delta})\right]\right)}\;\mathrm{d\xi^{2}}}{\leavevmode\resizebox{43.36464pt}{}{$\iint_{\mathbb{R}^{2}}$}|\;P(\vec{\xi})\;|^{2}\;\mathrm{d\xi^{2}}}\right|\;.$

(3)

Eq. (3) is motivated by Goodman [10] and the domain of integration is
determined by the aperture stop of the camera. In detail, $P$ describes the 2D
aperture stop function, which is given by a circular top hat function with
magnitude one and baseline zero. The displaced aperture stop function $P_{+}$
is shifted by:

$\vec{\Delta}\coloneqq\lambda z_{a\mapsto o}\dfrac{\vec{k}}{2}\approx\lambda
f\dfrac{\vec{k}}{2}\;.$ (4)

The intersection of $P_{+}$ and $P_{-}$ (shift by $-\vec{\Delta}$) determines
the domain of integration in Eq. (3). In addition, $\lambda$ characterizes the
wavelength under consideration and $z_{a\mapsto o}$ quantifies the distance on
the optical axis from the aperture stop to the observer plane, which roughly
equals the focal length $f$ of the camera lens if the Gaussian lens equation
is approximated. This simplification holds if the working distance is by
magnitudes larger than the image distance $z_{a\mapsto o}$. Finally, $\vec{k}$
denotes the spatial frequency vector.

As a side note, the $\mathrm{MTF}$ depends on the wavelength $\lambda$
according to Eq. (3). The polychromatic $\mathrm{MTF}$ can be retrieved by
integrating the conditioned, monochromatic $\mathrm{MTF}$ over the normalized
power spectral density $\mathrm{PSD}$ of the light source. The area under the
$\mathrm{PSD}$ curve quantifies the likelihood of emitting a photon in the
wavelength range $[\lambda,\;\lambda+\Delta\lambda]$ by the light source. As a
result, the polychromatic $\mathrm{MTF}$ is given by:

$\mathrm{MTF}(\vec{k})=\int\limits_{0}^{\infty}\mathrm{MTF}(\vec{k}|\lambda)\cdot\mathrm{PSD}(\lambda)\;\mathrm{d}\lambda\;.$
(5)

Consequently, chromatic aberrations will potentially also influence the
performance of AI-based algorithms for autonomous driving but they are not
discussed in detail in this paper.

At this point, it has to be emphasized, that the first three Zernike
coefficients for piston and tilt do not represent optical aberrations in the
classical sense. Even though there is a wavefront perturbation of the light
beam, the image quality is not influenced by those terms because the curvature
of the wavefront is not affected. Instead, tilts induce image distortions and
generate a non-conformal mapping. As we focus in this article on the sharpness
of the optical system we will not further investigate this influence.

If the Zernike polynomials of Table I are transformed into Cartesian
coordinates it becomes obvious that the difference in Eq. (3) vanishes for
Zernike polynomials of zeroth order. For the y-tilt $Z_{1}$ and the x-tilt
$Z_{2}$, the integrand evaluates to a constant phasor. Hence, it holds that:

$\displaystyle\begin{split}\mathrm{MTF}(\vec{k}|\lambda)&\overset{(\ref{eq:Zernike_parameterization})}{=}\left|\;\mathrm{e}^{\left(2\pi\mathrm{i}\cdot
c_{1,2}\cdot z_{a\mapsto
o}\cdot\vec{k}\right)}\cdot\dfrac{\oiint\limits_{P_{+}\cap
P_{-}}1\;\mathrm{d\xi^{2}}}{\iint_{\mathbb{R}^{2}}|\;P(\vec{\xi})\;|^{2}\;\mathrm{d\xi^{2}}}\;\right|\\\\[2.0pt]
\Leftrightarrow\mathrm{MTF}(\vec{k}|\lambda)&\;=\left|\;\mathrm{e}^{\left(2\pi\mathrm{i}\cdot
c_{1,2}\cdot z_{a\mapsto
o}\cdot\vec{k}\right)}\;\right|\cdot\dfrac{\oiint\limits_{P_{+}\cap
P_{-}}1\;\mathrm{d\xi^{2}}}{\iint_{\mathbb{R}^{2}}|\;P(\vec{\xi})\;|^{2}\;\mathrm{d\xi^{2}}}\\\\[5.0pt]
\Leftrightarrow\mathrm{MTF}(\vec{k}|\lambda)&\coloneqq\left|\;\mathrm{e}^{\left(2\pi\mathrm{i}\cdot
c_{1,2}\cdot z_{a\mapsto
o}\cdot\vec{k}\right)}\;\right|\cdot\mathrm{MTF_{diff}}(\vec{k}|\lambda)\\\\[5.0pt]
\Leftrightarrow\mathrm{MTF}(\vec{k}|\lambda)&\;=\mathrm{MTF_{diff}}(\vec{k}|\lambda)\;.\end{split}$
(6)

As a consequence, the diffraction limited $\mathrm{MTF}$ is not modulated by
the Zernike coefficients $c_{0}$ up to $c_{2}$. Therefore, the second order
Zernike coefficients are those of main interest for studying the effect of
optical aberrations in terms of sharpness degradation on convolutional neural
networks for autonomous driving. Nonetheless, the optical distortion does
influence the rectilinear mapping, but as stated before this is an independent
effect we do not study in this work.

Finally, note that the MTF is unfortunately currently not traceable to
fundamental physical quantities, and therefore a calibration chain to national
metrological institutes like the PTB in Germany or the NIST in the US can not
be established at the moment. For the automotive industry, this is a major
source of discussion, as the implementations of ISO12233 are very sensitive to
many diverse influences: stable lighting conditions (spectrum, intensity,
direction, homogeneity), a reproducible mechanical setup (target distance,
field-of-view) and well-defined camera settings (gain, exposure, HDR, ISP …)
are the goal, which are not met in practice. Instead, comparability of two
camera systems is only possible to a (relatively) good accuracy within
measurements from the very same experimental setup. Comparison between two
different measurement sites – even with the same nominal setup – is quite
difficult and error-prone [19, 27, 1].

##  III Refractive Power

In this section, we will demonstrate how relevant Zernike polynomials are not
captured by refractive power measurements. Our argument is based on the theory
that the refractive power is given by the second derivative of the wavefront
modulation of a plane wave passing a refractive element [29]. Using this
relationship, we demonstrate how several Zernike polynomials are simply not
covered by a refractive power measurement, as zeroth and first order
polynomials in $x$ and $y$ vanish if the second derivative is considered. We
start by introducing the measurement principle of a Shack-Hartmann sensor in
Sec. III-A to motivate how the Zernike coefficients $c_{i}$ are retrieved. In
Sec. III-C, we then mathematically derive which Zernike polynomials are not
present in the refractive power measurement. As the Zernike polynomials form
an orthogonal function basis, this proves that there are optical aberrations
that cannot be captured by the refractive power. Finally, we experimentally
demonstrate the validity of our assumption in Sec. III-D by using a Shack-
Hartmann sensor to measure the wavefront modulation induced by a high-quality
convex lens of well-known refractive power.

###  III-A From Shack-Hartmann Measurements to the Wavefront Aberration Map

Figure 2: Sketch of the measurement principle of a Shack-Hartmann sensor [34].

Fig. 2 demonstrates the method of operation of a Shack-Hartmann wavefront
sensor. If a collimated light beam is transmitting a refractive optical
element then the wavefront gets modulated. A Shack-Hartmann sensor consists of
a microlens array, which resolves the local wavefront perturbations by
focusing a wavefront snippet on a CCD or CMOS sensor. Without any aberrations,
the wavefront sensor will capture the light in the center of each pixel. If
aberrations are present, then the focusing spot will be displaced locally by
$d_{x}$ and $d_{y}$, respectively. The resulting local gradient of the optical
path difference ($W$) is given by:

$\displaystyle\begin{split}\vec{\beta}&\coloneqq\left[\leavevmode\resizebox{238.49231pt}{}{$\begin{array}[]{cccccc}\beta_{x}\Bigr{|}_{\vec{x}_{1}}&\ldots&\beta_{x}\Bigr{|}_{\vec{x}_{m}}&\beta_{y}\Bigr{|}_{\vec{x}_{1}}&\ldots&\beta_{y}\Bigr{|}_{\vec{x}_{m}}\end{array}$}\right]^{\mathrm{T}}\\\\[5.0pt]
\Leftrightarrow\vec{\beta}&\;=\left[\leavevmode\resizebox{336.05284pt}{}{$\begin{array}[]{cccccc}\dfrac{d_{x_{1}}}{\sqrt{f_{\mathrm{sh}}^{2}+d_{x_{1}}^{2}}}&\ldots&\dfrac{d_{x_{m}}}{\sqrt{f_{\mathrm{sh}}^{2}+d_{x_{m}}^{2}}}&\dfrac{d_{y_{1}}}{\sqrt{f_{\mathrm{sh}}^{2}+d_{y_{1}}^{2}}}&\ldots&\dfrac{d_{y_{m}}}{\sqrt{f_{\mathrm{sh}}^{2}+d_{y_{m}}^{2}}}\end{array}$}\right]^{\mathrm{T}}.\end{split}$
(7)

Here, $f_{\mathrm{sh}}$ denotes the focal length of the microlenses and $m$
specifies the number of microlenses within the array. With the Shack-Hartmann
measurement of the local wavefront gradients $\beta_{i}$, the Zernike
coefficients $c_{i}$ of the wavefront aberration map are determined by:

$\displaystyle\begin{split}\vec{\beta}=\dfrac{1}{\rho_{a}}\cdot\left[\leavevmode\resizebox{151.76964pt}{}{$\begin{array}[]{ccc}\dfrac{\partial
Z_{4}}{\partial\tilde{x}}\Bigr{|}_{\vec{x}_{1}}&\dots&\dfrac{\partial
Z_{n}}{\partial\tilde{x}}\Bigr{|}_{\vec{x}_{1}}\\\ \vdots&\ddots&\vdots\\\
\dfrac{\partial
Z_{4}}{\partial\tilde{x}}\Bigr{|}_{\vec{x}_{m}}&\dots&\dfrac{\partial
Z_{n}}{\partial\tilde{x}}\Bigr{|}_{\vec{x}_{m}}\\\\[8.0pt] \dfrac{\partial
Z_{4}}{\partial\tilde{y}}\Bigr{|}_{\vec{x}_{1}}&\dots&\dfrac{\partial
Z_{n}}{\partial\tilde{y}}\Bigr{|}_{\vec{x}_{1}}\\\ \vdots&\ddots&\vdots\\\
\dfrac{\partial
Z_{4}}{\partial\tilde{y}}\Bigr{|}_{\vec{x}_{m}}&\dots&\dfrac{\partial
Z_{n}}{\partial\tilde{y}}\Bigr{|}_{\vec{x}_{m}}\end{array}$}\right]\cdot\left[\begin{array}[]{c}c_{4}\\\
\vdots\\\
c_{n}\end{array}\right]\;\eqqcolon\;\dfrac{1}{\rho_{a}}\cdot\boldsymbol{\mathcal{M}}\cdot\vec{c}\;.\end{split}$
(8)

The Zernike decomposition coefficients $c_{i}$ are uniquely determined if
$|\boldsymbol{\mathcal{M}}^{T}\boldsymbol{\mathcal{M}}|\neq 0$. In other
words, the Gramian matrix
$\boldsymbol{\mathcal{M}}^{T}\boldsymbol{\mathcal{M}}$ has to be invertible,
wherefore $\boldsymbol{\mathcal{M}}^{T}\boldsymbol{\mathcal{M}}$ needs to have
full rank. If this condition is fulfilled, then the Zernike coefficient vector
$\vec{c}$ can be retrieved from the measured local wavefront gradient vector
$\vec{\beta}$ by:

$\vec{c}\overset{(\ref{eq:Zernike_numerical_analysis})}{=}\rho_{a}\cdot\left[\boldsymbol{\mathcal{M}}^{T}\boldsymbol{\mathcal{M}}\right]^{-1}\cdot\boldsymbol{\mathcal{M}}^{T}\cdot\vec{\beta}\;.$
(9)

###  III-B From Wavefront Aberration Maps to local Refractive Power

From Sec. III-A we know how to determine the Zernike coefficients $c_{i}$,
wherefore we can reconstruct the wavefront aberration map according to Eq.
(1). If the reference wavefront has been characterized by a plane wave, then
the local refractive power of an optical element is given by the second
derivative of the wavefront aberration map $W$ with respect to the axis of
interest [29, 30]. Hence, the refractive power $D_{x_{i}}$ along the axis
$x_{i}$ is given by:

$\displaystyle\begin{split}D_{x_{i}}(\vec{x}_{a})&=\dfrac{\partial^{2}}{\partial
x_{i}^{2}}W(\vec{x}_{a})\;.\end{split}$ (10)

Here, the input vector $\vec{x}_{a}\in\mathbb{R}^{2}$ is restricted to the
principal plane of the refractive element. The validity of Equation (10) can
be proven for the special case of a spherical thin lens:

$\displaystyle\begin{split}f_{x_{a_{1}}}^{2}&\;=x_{a_{1}}^{2}+\left(f_{x_{a_{1}}}-W(x_{a_{1}})\right)^{2}\;\text{,
w.l.o.g.:}\;\;x_{a_{2}}\overset{!}{=}0\\\\[5.0pt] \Rightarrow
W(x_{a_{1}})&\;=f_{x_{a_{1}}}\leavevmode\resizebox{151.76964pt}{}{$\left(1-\sqrt{1-\left(\dfrac{x_{a_{1}}}{f_{x_{a_{1}}}}\right)^{2}}\right)$}\\\\[5.0pt]
\Leftrightarrow
W(x_{a_{1}})&\;=f_{x_{a_{1}}}\leavevmode\resizebox{299.19885pt}{}{$\left(1-\left(1-\dfrac{1}{2}\left(\dfrac{x_{a_{1}}}{f_{x_{a_{1}}}}\right)^{2}+\mathcal{O}\left\\{\left(\dfrac{x_{a_{1}}}{f_{x_{a_{1}}}}\right)^{4}\right\\}\right)\right)$}\\\\[5.0pt]
\Rightarrow
W(x_{a_{1}})&\;\approx\dfrac{x_{a_{1}}^{2}}{2f_{x_{a_{1}}}}\eqqcolon\dfrac{D_{x_{a_{1}}}}{2}\cdot
x_{a_{1}}^{2}\\\\[5.0pt] \Rightarrow
D_{x_{a_{1}}}&\overset{(\ref{eq:refractive_power})}{=}\dfrac{\partial^{2}}{\partial
x_{a_{1}}^{2}}W(x_{a_{1}})=D_{x_{a_{1}}}\;.\;\;\;\;\;\;\blacksquare\end{split}$
(11)

###  III-C Information Content of Refractive Power Measurements

Eq. (10) determines the relationship between refractive power measurements $D$
and wavefront aberration measurements $W$ via the curvature of the optical
path difference map. In Sec. II-C we introduced the concept of the PSF as a
Fourier optical merit function, which serves as the impulse response function
or the Green’s function of an optical system [21]. In addition to the Fourier
optical approach there is also a ray optics approximation to describe the PSF
in terms of the area of a blurring ellipse, which encloses a certain amount of
light around the focusing spot in relation to the total amount of energy
entering the system through the aperture stop. The area of this blurring
ellipse is proportional to the Gaussian curvature [8] of the wavefront
aberration map or equivalently speaking, proportional to the determinant of
the Hessian matrix of the wavefront aberration function [29]:

$\displaystyle\begin{split}\oiint_{\mathcal{C}}\mathrm{PSF}(\vec{x}_{o})\Bigr{|}_{\hat{z}_{o}}\mathrm{d}^{2}x_{o}\propto\leavevmode\resizebox{216.81pt}{}{$\left|\left(\begin{array}[]{cc}\dfrac{\partial^{2}}{\partial
x_{1}^{2}}W(\vec{x}_{a})&\dfrac{\partial}{\partial
x_{1}}\dfrac{\partial}{\partial x_{2}}W(\vec{x}_{a})\\\\[10.0pt]
\dfrac{\partial}{\partial x_{1}}\dfrac{\partial}{\partial
x_{2}}W(\vec{x}_{a})&\dfrac{\partial^{2}}{\partial
x_{2}^{2}}W(\vec{x}_{a})\end{array}\right)\right|\;.$}\end{split}$ (12)

Here, $\mathcal{C}$ denotes the contour confining the domain of integration,
which is given by the blurring ellipse. Due to the relationship presented in
Eq. (10), this matrix is also known as the dioptric power matrix
$\boldsymbol{\mathcal{D}}$ [11]. The determinant can be rewritten in terms of
the traces of the dioptric power matrix:

$\displaystyle\begin{split}\oiint_{\mathcal{C}}\mathrm{PSF}(\vec{x}_{o})\Bigr{|}_{\hat{z}_{o}}\mathrm{d}^{2}x_{o}\overset{(\ref{eq:Hessian})}{\propto}\dfrac{1}{2}\left[\left(\mathrm{tr}\left(\boldsymbol{\mathcal{D}}\right)\right)^{2}-\mathrm{tr}\left(\boldsymbol{\mathcal{D}}^{2}\right)\right]\;.\end{split}$
(13)

So far, the automotive industry exclusively specifies requirements in terms of
the refractive power w.r.t. the horizontal and vertical directions.
Consequently, only the trace of $\boldsymbol{\mathcal{D}}$ is measured and
off-diagonal elements in the Hessian matrix are not investigated. This
demonstrates that there is a blind spot in the quality assurance chain at the
moment. This conclusion can be further underpinned by an mathematical
argument. The trace of $\boldsymbol{\mathcal{D}}$ is given by:

$\displaystyle\begin{split}\mathrm{tr}\left(\boldsymbol{\mathcal{D}}\right)=\sum\limits_{i=1}^{d}D_{x_{i}}(\vec{x}_{a})=\mathop{}\\!\mathbin{\bigtriangleup}W(\vec{x}_{a})\;.\end{split}$
(14)

Consequently, the trace of $\boldsymbol{\mathcal{D}}$ is unaffected by
wavefront aberration fields, which fulfill the Laplace equation:

$\mathop{}\\!\mathbin{\bigtriangleup}\Gamma(\vec{x}_{a})\overset{!}{=}0\;.$
(15)

As a result, the trace of $\boldsymbol{\mathcal{D}}$ is Gauge invariant under
aberration fields $\Gamma(\vec{x}_{a})$ that are composed of harmonic
functions. Hence, Zernike polynomials in Table I that are harmonic functions
(like astigmatism or trefoil) will not alter the trace of
$\boldsymbol{\mathcal{D}}$.

In a nutshell, refractive power measurements are not sensitive for optical
distortions quantified by $c_{1}$ and $c_{2}$. Furthermore, the refractive
power is invariant under oblique astigmatism given by $c_{3}$ if the
refractive power requirements are specified exclusively along the horizontal
and vertical axis, as it is the current governing standard in the automotive
industry. Finally, those quality standards are insufficient for extracting
more fundamental information about the optical system in terms of the
$\mathrm{PSF}$. Nonetheless, the aberrations associated with these polynomials
have been proven to have an influence on the performance of AI algorithms [26,
24].

###  III-D Experimental Verification

Since Eq. (10) is not well established in the automotive industry, we
experimentally demonstrate the validity of the relationship by a Shack-
Hartmann wavefront measurement of a calibration lens. The lens under test was
produced by Zeiss and is traced back to national standards by an accredited
calibration authority. The local wavefront gradients $\beta_{i}$ are measured
by a Shack-Hartmann sensor and the refractive power is retrieved by utilizing
Eq. (10). As demonstrated by Eq. 8 the Shack-Hartmanm measurement yields the
first derivative of the wavefront. Here, we measure the lens and numerically
determine the second derivative by a simple central difference scheme, which
should result in the specified refractive power. Fig. 3 illustrates the
outcome w.r.t. the refractive power map over the lens aperture. From the
frequency distribution of the local refractive power across the entire
principal plane, the expectation value for the global refractive power of the
optical element in the $x$\- and $y$-plane can be deduced. The expectation
values meet the certified refractive power values of the calibration lens
within the uncertainty intervals. Hence, the validity of Eq. (10) has also
been experimentally confirmed.

Figure 3: Wavefront measurement performed on a ${\left<D\right>=100.3}$ mdpt
reference lens. In order to cover the entire aperture of the lens, several
Shack-Hartmann measurements have been stitched together. This procedure has
introduced artifacts, which are visible in the measurement data by strongly
pronounced vertical and horizontal lines. In total, 15 measurements have been
performed over the calibration lens aperture of $d=10$ cm.

Summarizing, we have demonstrated that fundamentally several optical
aberrations are not captured by a refractive power measurement. The image
quality can be deteriorated even though the refractive power measurement
indicates a compliant windscreen sample. From previous studies on the effect
of oblique astigmatism ($c_{3}$) on road sign classification [26, 24] it
becomes evident, that refractive power measurements are insufficient for
specifying the quality of a windshield in order to ensure reliable computer
vision for autonomous driving vehicles.

##  IV Modulation Transfer Function

In this section, we will demonstrate why the windscreen and the camera form a
joint optical system that cannot be separated into two independent
constituents, such that the MTF cannot be determined for the two systems
separately. First, we argue how the refractive power of the windscreen
interacts with the focal length of the camera system. In a second step, this
is experimentally verified using a MTF measurement with and without a
windscreen. A discussion elaborates on several implications for the production
and testing process.

###  IV-A Field Curvature

The focal length of an imaging system varies over the field of view, with the
so-called _field curvature_ being a prominent optimization goal for any lens
designer. The semiconductor production processes produce completely flat image
sensors, which for the imaging optics is a challenge, as the field curvature
needs to be flat as well to minimize aberrations. This field curvature, as a
design property of the lens, is given by the offset $\Delta z_{\text{fc}}$
over field in units of length, typically on the micrometer range. A symbolic
field curvature is visualized in Fig. 4.

Figure 4: Windscreen and lens form a joint optical system. $H$ and
$H^{\prime}$ are the principle planes of the lens, $f$ is the nominal focal
length. The blue line visualizes the field curvature (not to scale). Normally,
parallel rays are focused onto the field curvature (yellow line). Windscreen
refractive power shortens or prolongs the effective focal length of the lens
(red line). There are two different focus offsets $\Delta z_{\text{fc}}$ and
$\Delta z_{\text{ws}}$ which may add or even cancel at different fields of
view.

As explained above, the refractive power of the windscreen leads to parallel
rays converging or diverging. Taking the two elements windscreen and lens
together yields a second focus offset $\Delta z_{\text{ws}}$ for the camera
system, as the converging (diverging) rays will shorten (prolong) the focal
length of the camera system. Fig. 4 depicts this situation. The two offsets
are added for the system offset, such that:

$\Delta z=\Delta z_{\text{ws}}+\Delta z_{\text{fc}}\;.$ (16)

Importantly, both $\Delta z_{\text{ws}}$ and $\Delta z_{\text{fc}}$ can have
positive or negative values, and thus the overall offset may vanish when these
terms cancel. A vanishing offset value implies a sharpening of the system.
Here, a MTF measurement of the camera alone would yield a certain number,
while putting a windscreen in front of the camera would act like glasses and
the image would become sharper. That this is indeed the case in practice is
presented in the following section.

###  IV-B Experimental Validation

(a) ADAS camera

(b) ADAS camera + windshield

Figure 5: MTF measurement for an ADAS system based on the Slanted Edge method
according to ISO 12233 [16].

Fig. 5 depicts two slanted edge measurements, one without a windscreen (5(a)),
and one with a windscreen placed in front of the camera system (5(b)). The
insets indicate the MTF values derived from the numerical evaluation for all
four edges, using an ISO12233-compliant algorithm [16]. There are two
horizontal and two vertical values. The two vertical values (top and bottom)
distinctly decrease from $52\pm 1.5\thinspace\%~{}[95\thinspace\%]$ to $39\pm
1.7\thinspace\%~{}[95\thinspace\%]$ when a windscreen is placed in front of
the camera. However, for the horizontal direction (left and right) the MTF
values both significantly increase from $45/47\pm
1.5\thinspace\%~{}[95\thinspace\%]$ to $52/54\pm
1.6\thinspace\%~{}[95\thinspace\%]$ when the windscreen is placed in front of
the camera. The results experimentally confirm that the defocus $\Delta
z_{\text{ws}}$ and $\Delta z_{\text{fc}}$ may cancel to a certain degree,
increasing the sharpness like glasses would do for a myopic person. This
conclusion is well established in physics [12] for decades but the
implications for the quality assurance testing procedure of ADAS systems in
the automotive industry are not well prevalent.

###  IV-C Discussion

Both the field curvature of the lens and the refractive power of the
windscreen are spatially variant. The field curvature not only varies over
field (radially) as the name implies, but due to production tolerances the
rotation symmetry of the lens is usually broken to a certain degree. The field
of view of the lens projected on the windscreen yields a trapezoidal cutout
(cf. Fig. 1). I.e., the (almost) rotational symmetry of the lens projected on
the windscreen combines with the local refractive power variation of the
windscreen in this cutout. Not only that, but the windscreen also has
distinctly different refractive power for the horizontal and the vertical
direction as given by the Kerkhof model [33].

Taken together, it is apparent that a windscreen cannot be qualified by a MTF
measurement if both the windscreen and the camera are measured separately. The
experimental sharpening unambiguously demonstrates a non-linear process,
proving that the two elements cannot be separated using linear system theory
(read: the individual MTFs cannot be multiplied). The way individual
production tolerances will add is not predictable. In brief, it is not
possible to determine individual MTF limits, as the combination of individual
tolerances may hold both good and bad surprises, either sending a good system
to scrap or a bad system to the field.

Therefore, any solution using MTF would have to measure the MTF on the
combined system of produced windscreen and camera system with their individual
production tolerances. This could be either at the production site of the Tier
1 or the OEM. But there are still several important open questions that make
this an unattractive proposal: if an assembly is non-compliant, is it worth
finding a compliant combination, does it make economically sense? How big are
the assembly tolerances fitting the windscreen into the car body? If the OEM
wants the measurement system at the site of the Tier 1, one should be aware
that the assembly of the windscreen into the car produces distinct mechanical
tolerances, changing the shape and internal tension of the windscreen. As we
are looking for subtle differences in optical quality, this may affect the
pre-assembled camera system as well. Finally, from an automotive process view
it is clear that an independent measurement of the camera and the windscreen
is much preferred.

Summarizing, the MTF is a measure of ’sharpness’ based on linear system
theory. The windscreen and the camera form a combined optical system that
cannot be separated, which prohibits its use for windscreen characterization
without the actual, produced camera system in place. Taken together with the
possibility of finding a better metric we are skeptical that the MTF should be
prioritized for windscreen characterization going forward.

##  V Simulating Optical Properties

Having shown that basically no current measurement system in the automotive
windscreen industry is capable of a meaningful characterization of the
windscreen optical quality for downstream AI algorithm consumption, what could
be a way forward? A comprehensive experimental study using thousands of actual
cameras and windscreens is out of the question. Therefore, the AI performance
needs to be linked to the windscreen optical quality by simulation, using
physical-realistic optic models. These simulations need to model the
production tolerances of both windscreens and cameras. Then, the performance
requirements of the AI-based ADAS functionalities can be translated to optical
quality specifications for windscreen production.

The waveform description is fundamental and includes all optical effects and
aberrations, and can be measured by a Shack-Hartmann sensor. Currently, this
is not a viable approach to windscreen characterization at the site of the
Tier 1, as it is too expensive and more importantly, too slow for a 100 % part
check. Nonetheless, we believe it is possible to use special laboratory-grade
equipment to create the physical-realistic optical models necessary for the
simulations, and then derive from these simulations an understanding of the
optical properties that are really necessary for the AI performance. Finally,
from this we can deduce a simplified form of measurement that captures this
newfound knowledge of the required optical properties. A first example of this
process is published in [20].

Therefore, the challenge is to understand those optical properties that are
really necessary for a robust AI algorithm performance. We believe that this
is a necessary step, and without it, the move from ADAS to AD will be
prohibitively difficult, as production tolerances combined with the complexity
of the world create an unmanageable number of combinations.

##  VI Summary

In automotive mass production, the inspection systems at the suppliers need to
measure the quality of windscreens in a meaningful way for the final device
performance. Modern ADAS and future AD camera systems are based on AI
algorithms, wherefore the windscreen quality needs to be linked to the
performance of these algorithms. Currently, there are two types of
measurements established in the industry to measure the optical quality of a
windscreen: refractive power and MTF. In this article, we demonstrated how
both these measurements are fundamentally not capable of capturing relevant
optical properties of the windscreen.

The refractive power measurement does not include several aberrations – given
by Zernike polynomials, e.g. oblique astigmatism – while these aberrations
obviously affect the performance of the AI algorithms: oblique astigmatism
causes a directional blurring of the scene, and blurring causes a degradation
of the performance. Because of the orthogonality of the Zernike polynomials,
it is clear that this information is simply lacking in refractive power
measurements.

MTF is based on linear system theory, where independent optical systems might
be multiplied in frequency space to yield the system MTF. This is, for
example, the case for the lens and the imager. Here, we demonstrated
mathematically and experimentally that the windscreen forms a novel optical
system together with the lens of the camera system, which cannot be separated
into individual components. Therefore, measuring the MTF on the windscreen
alone will not yield the performance of the combined system. Thus, the final
assembly of the windscreen and camera system in the car may be both better or
worse than the EOL measurement at the windscreen production site, either
sending good parts to scrap or bad parts into the field.

Every car has a windscreen. Using the knowledge presented in this article we
believe that the automotive industry needs to focus their efforts on finding
novel measurement methods that qualify the optical quality of windscreens in a
meaningful way for the downstream AI algorithms. We propose a concept using
fundamental wave (and Fourier) optics to characterize the windscreens and
combine wavefront measurements and physical-realistic simulations to reach an
understanding of what optical properties are really important for AI computer
vision algorithms. We believe that it cannot be said generally that the
optical quality of windscreens is too low – what is currently lacking is not
optical quality, but understanding about the robustness of the algorithms
against optical aberrations. We simply do not know what optical quality is
needed exactly. Taking these elements together we believe that a novel metric
can be found that contains the relevant information, while at the same time
the measurement is practical enough to be used stand-alone at the windscreen
production site. This is the great windscreen challenge the automotive
industry currently faces.

## References

  * [1] Uwe Artmann.  Quantify aliasing a new approach to make resolution measurement more robust.  Electronic Imaging, 31(10):320, 2019.
  * [2] Glenn D. Boreman.  Modulation Transfer Function in Optical and Electro-Optical Systems, Second Edition.  SPIE Press, 2001.
  * [3] Max Born, et al.  Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light.  Cambridge University Press, 7 edition, 1999.
  * [4] Alexander Braun.  Automotive mass production of camera systems: Linking image quality to AI performance.  tm - Technisches Messen, 2022.
  * [5] David Compertore.  Adas windshield measurements - white paper.  Technical report, Lumetrics, 2021.
  * [6] DIN:52305.  Determination of the optical deviation and refractive power of safety glass for vehicle glazing.  Technical report, German Institute for Standardization, 1995.
  * [7] DIN:52353.  Prüfung der abriebfestigkeit von fahrzeugverglasung mit dem wischer-test.  Technical report, Deutsches Institut für Normung (DIN), 2019.
  * [8] Manfredo Perdigão do Carmo.  Differential geometry of curves and surfaces.  Prentice-Hall, 1976.
  * [9] ECE:43.  Uniform provisions concerning the approval of safety glazing materials and their installation on vehicles.  Technical report, Regulation No. 43 of the Economic Commission for Europe of the United Nations (UN/ECE), 2014.
  * [10] Joseph W. Goodman.  Introduction to Fourier optics.  Stanford University / McGraw-Hill, 1968.
  * [11] William Harris.  The matrix representation of dioptric power. part 1: An introduction.  South African Optometrist, 47:19–23, 01 1988.
  * [12] Eugene Hecht.  Optik.  De Gruyter, 2018.
  * [13] Michael Heizmann, et al.  Implementing machine learning: chances and challenges.  at – Automatisierungstechnik, 70(1):90–101, 2022.
  * [14] IEEE:P2020.  Automotive imaging - white paper.  Technical report, IEEE Standards Association, 2018.
  * [15] M. J. Irland.  Windshield optics.  Applied Optics, 8(9):1787–1790, 1969.
  * [16] ISO:12233.  Photography - electronic still picture imaging - resolution and spatial frequency responses.  Technical report, International Organization for Standardization, 2023\.
  * [17] ISO:24157.  Ophthalmic optics and instruments - reporting aberrations of the human eye.  Technical report, International Organization for Standardization, 2008\.
  * [18] ISO/TS:16949.  Quality management systems — particular requirements for the application of iso 9001:2008 for automotive production and relevant service part organizations.  Technical report, International Organization for Standardization, 2009\.
  * [19] Norman Koren, et al.  Measuring mtf with wedges: pitfalls and best practices.  Electronic Imaging, 29(19):6, 2017.
  * [20] Christian Krebs, et al.  Impact of windshield optical aberrations on visual range camera based classification tasks performed by cnns.  Proc. London Imaging Meeting: Imaging for Deep Learning, pages 83–87, 2021.
  * [21] Prem K. Kythe.  Green’s Functions and Linear Differential Equations: Theory, Applications, and Computation.  CRC Press, 2011.
  * [22] LaVision.  Automotive imaging for safer driving.  Technical report, LaVision, 2023.
  * [23] Thomas Mitra.  Benefits of optical distortion measurement.  Technical report, ISRA Vision, 2022.
  * [24] Patrick Müller et al.  Simulating optical properties to access novel metrological parameter ranges and the impact of different model approximations.  In 2022 IEEE International Workshop on Metrology for Automotive (MetroAutomotive), pages 133–138, 2022.
  * [25] Patrick Müller et al.  Mtf as a performance indicator for ai algorithms?  Electronic Imaging, 35(AVM-125), 2023.
  * [26] Patrick Müller, et al.  Impact of realistic properties of the point spread function on classification tasks to reveal a possible distribution shift.  NeurIPS Workshop DistShift, 2022.
  * [27] J. Roland.  A study of slanted-edge mtf stability and repeatability.  In Electronic imaging, 2015.
  * [28] Carsten Steger, et al.  Machine Vision Algorithms and Applications.  Wiley-VCH, Weinheim, 2nd edition, 2018.
  * [29] Larry N. Thibos.  Calculation of the geometrical point-spread function from wavefront aberrations.  Ophthalmic and Physiological Optics, 39(4):232–244, June 2019.
  * [30] L. N. Thibos, et al.  Accuracy and precision of objective refraction from wavefront aberrations.  Journal of Vision, 4(4):9, April 2004.
  * [31] VDA:6.3.  Qualitätsmanagement in der automobilindustrie - prozessaudit.  Technical Report 4, Verband der Deutschen Automobilindustrie (VDA), 2023\.
  * [32] Korbinian Weik, et al.  Imaging through curved glass: windshield optical impact on automotive cameras.  Proceedings of SPIE, 12231(122310A), 2022.
  * [33] Dominik W. Wolf, et al.  Optical power measurement in the automotive world.  Currently under review at Metrologia, 2023.
  * [34] Z. Yang, et al.  Generalized hartmann-shack array of dielectric metalens sub-arrays for polarimetric beam profiling.  Nature Commun, 9(4607), 2018.