# Magneto-Optical Detection of the Orbital Hall Effect in Chromium

Igor Lyalin  Department of Physics, The Ohio State University, Columbus, Ohio
43210, United States  Sanaz Alikhah  Department of Physics and Astronomy,
Uppsala University, P. O. Box 516, SE-75120 Uppsala, Sweden  Marco Berritta
Department of Physics and Astronomy, Uppsala University, P. O. Box 516,
SE-75120 Uppsala, Sweden  Department of Physics and Astronomy, University of
Exeter, Stocker Road, Exeter EX4 4QL, United Kingdom  Peter M. Oppeneer
Department of Physics and Astronomy, Uppsala University, P. O. Box 516,
SE-75120 Uppsala, Sweden  Roland K. Kawakami  Department of Physics, The Ohio
State University, Columbus, Ohio 43210, United States

###### Abstract

The orbital Hall effect has been theoretically predicted but its direct
observation is a challenge. Here, we report the magneto-optical detection of
current-induced orbital accumulation at the surface of a light 3$d$ transition
metal, Cr. The orbital polarization is in-plane, transverse to the current
direction, and scales linearly with current density, consistent with the
orbital Hall effect. Comparing the thickness-dependent magneto-optical
measurements with ab initio calculations, we estimate an orbital diffusion
length in Cr of $6.6\pm 0.6$ nm.

††preprint: APS/123-QED

The spin Hall effect (SHE) allows for the conversion of an electric charge
current into a transverse spin current [1, 2, 3, 4]. It is one of the key
phenomena in the field of spintronics and has been extensively studied over
the past two decades [5, 6]. The SHE has become an indispensable tool to
generate and detect spin currents, and has been envisioned to play an
important role in many spintronic devices, including next-generation magnetic
memory based on spin-orbit torque (SOT). Recent theoretical studies have
indicated that in many materials the SHE arises due to the orbital Hall effect
(OHE) [7, 8, 9, 10, 11], arguing that the OHE is more fundamental than the SHE
[8, 9]. The OHE converts a charge current into a transverse flow of orbital
angular momentum (OAM). It is predicted theoretically that unlike the SHE, the
OHE does not require spin-orbit coupling (SOC), while in the presence of SOC,
a current of orbital angular momentum can be efficiently converted into a spin
current, thus converting the OHE into the SHE [7, 10].

A few experimental studies of the orbital Hall effect were conducted so far
[12, 13, 14, 15, 16, 17]. They focused on measuring SOT in non-
magnetic/ferromagnetic (NM/FM) bilayer systems. Unexpectedly large SOT signals
in seemingly trivial systems without a heavy metal, e.g. CuO/FM, Cr/FM, were
attributed to the OHE [12, 13, 14, 16]. In such systems, however,
contributions to the total measured torque acting on the FM magnetization from
Rashba-type SOC at the interface [18, 19] and anomalous SOT [20, 21] generated
by the FM itself complicate the analysis [22]. It has also been recently shown
that while experimentally measurable total torque is usually large and
negative for Ta/FM, it depends on the choice of FM and is positive in Ta/Ni
bilayers [14] as well as in Cr/Ni [13, 16]. This was attributed to the torque
due to OHE dominating over SHE-induced torque. The interplay between these
torques of different origin makes it difficult to disentangle a torque due to
the OHE, making experimental detection of the OHE less direct. Therefore, in
order to gain an insight into the OHE it would be advantageous to study it in
a single non-magnetic layer.

The magneto-optical Kerr effect (MOKE) has proved to be a valuable technique
to detect spin accumulation due to the SHE in both semiconductors [3] and
metals [23]. While MOKE cannot distinguish the signals from the orbital and
spin accumulation, experiments on materials with a negligible SOC have a great
advantage. In these materials, the OHE can be 1-2 order of magnitude larger
than the SHE, as a result, surface orbital accumulation is expected to be
significantly larger than spin accumulation [10, 24]. Furthermore, MOKE is
sensitive to the orbital magnetization because the electric dipole transition
couples light helicity directly to the orbital degree of freedom, whereas the
detection of spin polarization requires spin-orbit coupling. Therefore,
signals from SHE can be reduced by choosing materials with weak SOC, e.g.
light 3$d$ transition metals, which facilitates the separate detection of OHE
as pursued recently by Choi et al. [25], in contrast to materials with strong
SOC where both OHE and SHE coexist.

In this Letter, we demonstrate the detection of the orbital Hall effect in one
such material, namely, in single crystal bcc Cr(001), using longitudinal MOKE
combined with ac current modulation. MOKE microscopy with a sensitivity of
about $2$ nrad allows us to detect a current-induced magnetization at the
surface of Cr which, in combination with ab initio calculations of the OHE,
SHE, and MOKE in Cr, we demonstrate to be the orbital accumulation due to the
OHE. Comparison of the experimental data with the theoretical calculations
yields quantitative values for the orbital Hall angle and the orbital
diffusion length of Cr. Thus, our results demonstrate a possibility to detect
and quantify the OHE in a NM, independent of an adjacent FM. Along with Choi
et al. [25], this provides strong evidence for the OHE.

Experiments are performed on Cr(001) films grown by molecular beam epitaxy on
MgO(001) substrates with film thicknesses ranging between 7.5 and 50 nm.
Details on the growth and characterization of Cr films is given in
Supplemental Material, Sec. S1 (SM:S1) [26]. All Cr films are capped with 5 nm
CaF2 to protect them from oxidation, unless otherwise indicated. The samples
are lithographically patterned into $20\,\mu$m wide channels for the MOKE
measurements and Hall bars for resistivity measurements.

In the longitudinal MOKE measurements, a s-polarized laser beam (to exclude
contribution from transverse MOKE) with wavelength $\lambda=800$ nm is
incident at $45\degree$ from the sample normal. The laser beam is focused to a
$\sim$3 $\mu$m spot size on the sample using an objective lens. The reflected
light is collected using a second objective positioned at $90\degree$ to the
first one. The sample is mounted on an xy translation stage which allows
scanning across the sample surface. A sine-modulated current is applied to the
lithographically patterned device (patterned along Cr [100] direction),
inducing orbital accumulation on the top and bottom surfaces. The Kerr
rotation of the reflected light is measured using a combination of a half-wave
plate, a polarization-splitting Wollaston prism, a balanced photodetector, and
a lock-in amplifier at the frequency of the current modulation (details in
SM:S2,S3).

(a)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(a)

(b)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(b)

Figure 1: (a) Schematics of the orbital Hall effect. The charge current $j$
generates a transverse orbital current, leading to orbital accumulation on the
sample’s surfaces. (b) Measurement setup utilizing the longitudinal MOKE to
detect the in-plane orbital accumulation.

Figure 1b shows the longitudinal MOKE measurement geometry. The OHE generates
an orbital accumulation oriented along the in-plane y direction transverse to
the electric current flowing along x. The generated orbital accumulation has
opposite signs at the opposite surfaces of the Cr film. The longitudinal MOKE
detects the superposition of the polarization rotation caused by the
y-polarized OAM accumulated at the top and bottom surfaces. For quantitative
analysis, one needs to take into account the light attenuation in a metallic
film, with the characteristic penetration length of the light usually on the
order of 10-30 nm. At the same time, the orbital diffusion length determines
the spatial distribution and the amount of OAM accumulation, which is reduced
for films thinner than $l_{o}$. For sufficiently thick films these two effects
will lead to a saturation of the OAM accumulation detected by MOKE. Therefore,
a thickness-dependent study can estimate orbital diffusion length in a single
layer film. This is important as the characteristic length of orbital
transport is debated theoretically, with some works predicting very short
orbital accumulation distances on the order of several atomic layers [27], and
other works predicting orbital diffusion distances longer than the spin
diffusion length [25, 28, 29].

(a)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(a)

(b)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(b)

(c)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(c)

(d)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(d)

Figure 2: Magneto-optical Kerr rotation induced by the orbital and spin Hall
effect. (a) MOKE line scan across a 20 $\mu$m wide, 30 nm thick Cr channel at
rms current density $j=0.28\times 10^{11}$ A/m2. MOKE signal as a function of
current density: (b) for the laser beam incidence shown in Fig. 1b, (c) for
the reversed optical path of the laser beam (30 nm sample), (d) for device
rotated $90\degree$ about the surface normal (40 nm sample).

Stamm et al. [23] has shown that special care should be taken to detect MOKE
signals due to the SHE in Pt and W in order to avoid artificial signals due to
a change in reflectivity of the sample caused by Joule heating $\sim j^{2}$
(SM:S4). Similar to Stamm et al. [23], we measure the current-induced Kerr
rotation as the voltage output of the balanced photodetector at the frequency
of the sinusoidal driving current and subtraction of MOKE signals measured
with opposite current polarities (inverted by changing the wiring of the
device) is performed to exclude any remaining thermal artifacts. Figure 2a
shows the Kerr rotation measured on a 30 nm thick Cr film at current density
of $0.26\times 10^{11}$ A/m2. We observe a Kerr rotation signal from the
surface of the conducting Cr channel, which is of the order of a few
nanoradians. The MOKE signal is approximately constant over the Cr channel
surface (SM:S5) and scales linearly with the applied current, consistent with
the orbital/spin accumulation picture. The absence of antisymmetry in the y
spatial profile in Fig. 2a as well as measurements on a Cr/MgO/Cr sample show
that the MOKE signal does not originate from an in-plane Oersted field (see
SM:S6). The value of Kerr rotation per current density for a 30 nm Cr sample
is $12.2\pm 0.7$ nrad per $10^{11}$ A/m2 as found from a linear fit to the
data shown in Fig. 2b.

Further evidence that the measured Kerr rotation is consistent with the
y-polarized orbital accumulation in Cr comes from the following control
measurements.

(a)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(a)

(b)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(b)

(c)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(c)

(d)left-2pt+/2right-3pt-/2up+5pt-/2down3pt+/2[left,up][vc,hc]

(d)

Figure 3: Magneto-optical Kerr rotation induced by the orbital Hall effect.
(a) MOKE line scan across a 20 $\mu$m wide Cr/MgO (30/3 nm) channel at current
density $j=0.26\times 10^{11}$ A/m2 (b) MOKE line scan across a 40 nm thick Pt
channel at $j=0.10\times 10^{11}$ A/m2 (c) and (d) MOKE signal as a function
of current density for Cr and Pt, respectively.

First, the longitudinal MOKE geometry chosen in our work is sensitive to in-
plane y and out-of-plane z components of magnetization. Upon inversion of the
light incidence angle, the Kerr rotation due to the y component should change
sign, while the Kerr rotation due to the z component should stay constant.
Figure 2c shows that the Kerr rotation measured with the opposite angles of
incidence reverses sign. The linear fits to the data in Fig. 2b and 2c yield
slopes that, within error bars, are the same. The linear dependence with
current density demonstrates the absence of a thermally induced signal and the
same magnitudes of the slopes exclude the presence of a polar MOKE
contribution due to a magnetization along z.

Second, if the measured MOKE signal originates from an orbital accumulation
due to the OHE, then when the sample is rotated by $90\degree$ about the
normal, the Kerr rotation should decrease to zero, as the orbital accumulation
is polarized along x in this case. Figure 2d shows that upon rotating sample
by $90\degree$ Kerr signal indeed decreases and can be said to be zero within
error bars (SM:S7). This behavior is consistent with the OHE/SHE origin of the
signal.

Third, to investigate whether the observed signals are sensitive to the
interface, we insert 3 nm of MgO in between the Cr and CaF2 capping layer.
Fig. 3c shows that the different interface does not change the magnitude of
the signal. The measured Kerr rotation is the same, within error bars, for the
devices with Cr/MgO and Cr/CaF2 interfaces. This result demonstrates that it
is unlikely that the origin of the effect is interfacial (e.g., spin or
orbital Rashba-Edelstein effect), and points to the bulk origin of the effect,
as expected theoretically for OHE.

Fourth, we compare the signs of Kerr rotation in Cr and Pt films. Fig. 3b and
3d show the experimentally measured MOKE signal for a sputtered 40 nm Pt
sample. The value of Kerr rotation per current density has a value smaller in
magnitude than measured in Cr and of the opposite sign, $-3.2\pm 1.5$ nrad per
$10^{11}$ A/m2. We find that the opposite signs in Cr and Pt are consistent
with the opposite Kerr coefficients theoretically calculated for Cr and Pt at
800 nm wavelength (SM:S11). To follow the sign convention in Stamm et al.
[23], we choose the sign of current-induced Kerr rotation in Pt to be negative
(see SM:S8). To confirm the magnitude and sign of MOKE in Pt, we conducted
measurements on an MBE-grown Pt sample which, within error bars, yields the
same MOKE values (SM:S9).

Fifth, to examine the orbital diffusion in Cr we perform ab initio
calculations of the OHE, SHE, and wavelength-dependent longitudinal MOKE. To
compute the magneto-optical signal due to a spin or an orbital polarization in
Cr metal, we model the induced spin and orbital moment by including a Zeeman
term in the relativistic DFT Hamiltonian, $\mathcal{H}=\mu_{\rm
B}\bm{B}_{ext}\cdot(\bm{\ell}+2\bm{s})$, with $\bm{B}_{ext}$ the applied field
and $\bm{\ell}$ and $\bm{s}$ the orbital and spin angular momentum operators,
respectively. Coupling $\bm{B}_{ext}$ selectively to either spin or orbital
angular momentum gives us an induced spin or orbital polarization. We then use
the linear-response formalism to compute the magneto-optical spectrum [30] for
either a spin or an orbital polarization in bulk Cr (SM:S10). Note that when
we couple e.g. the Zeeman field only to the spin, there will nonetheless be a
very small orbital moment generated by the nonzero spin-orbit coupling. We
find that the magneto-optical conductivity can be approximated as
$\sigma_{zx}(\omega)=\sigma_{zx}^{s}(\omega)m_{s}+\sigma_{zx}^{o}(\omega)m_{\ell}$,
i.e., there are two distinct contributions that are linear in the spin or
orbital moment. In Fig. 4(a) we show the ab initio calculated magneto-optical
conductivity spectra for induced spin or orbital polarizations. It can be
recognized immediately that an induced orbital moment leads to a much larger
$\textrm{Im}\,\sigma_{zx}(\omega)$ than an induced spin moment. Next, we use
the calculated optical conductivity $\sigma_{xx}(\omega)$ and compute the
resulting longitudinal MOKE for a $45^{\circ}$-angle of incidence (see [30]
and SM:S10 for details). The calculated Kerr rotation spectra for bulk Cr are
shown in Fig. 4(b). Again, we observe that the Kerr rotation generated by
orbital polarization is substantially larger ($>10\times$) than that due to
spin polarization.

Figure 4: (a) The ab initio calculated magneto-optical conductivity
$\textrm{Im}\,\sigma_{zx}(\omega)$ of Cr for an induced spin moment $m_{s}$ or
orbital moment $m_{\ell}$ of 0.1 $\mu_{\rm B}$. (b) The calculated
longitudinal Kerr rotation $\theta_{\rm K}$ of Cr due to a spin or an orbital
moment. Note that $\theta_{\rm K}$ for the induced spin polarization is
multiplied by 10.

After verifying that detected MOKE signal in Cr (001) is not a thermal
artifact and fully consistent with the symmetries of the OHE, we investigate
the characteristic length of the orbital transport by studying the Kerr
rotation as a function of Cr film thickness, $t$. Figure 5 shows the
experimental results for the measured Kerr rotation normalized by a current
density $j=10^{11}$ A/m2. Intuitively, three effects can contribute to the
thickness dependence: the saturation of orbital accumulation with a length
scale of $l_{o}$ as discussed earlier, the limited penetration depth of the
light ($\sim 18$ nm in Cr, see SM:S10), and the opposite sign of the orbital
accumulation at the top and bottom interfaces.

To model the thickness-dependent MOKE we use the formalism of Ref. [23], but
suitably modified to describe separately the Kerr rotations due to orbital and
spin accumulation, respectively. Thus, we assume that orbital diffusion is
described by expressions equivalent to those for spin diffusion. The
thickness-dependent Kerr rotation due to orbital accumulation is given by

$\displaystyle\\!\\!\\!\\!\theta_{\rm K}^{o}$ $\displaystyle=$
$\displaystyle\frac{l_{o}\sigma_{zx}^{\rm OH}\,\rho(t)^{2}jD(E_{\rm
F})\,\mathrm{e}^{\frac{t}{2l_{o}}}}{\cosh(\frac{t}{2l_{o}})}\times$ (1)
$\displaystyle{\rm Re}$ $\displaystyle\left\\{\Phi_{\rm
K}^{bulk,o}\kappa\bigg{(}\frac{(\mathrm{e}^{-\kappa^{-}t}-1)\mathrm{e}^{-\frac{t}{l_{o}}}}{\kappa^{-}}-\frac{\mathrm{e}^{-\kappa^{+}t}-1}{\kappa^{+}}\bigg{)}\\!\right\\}\\!,$

where $\sigma_{zx}^{\rm OH}$ is the bulk OHE, $\rho(t)$ the resistivity,
$D(E_{\rm F})$ the density of states at the Fermi energy $E_{\rm F}$,
$\Phi_{\rm K}^{bulk,o}$ is the bulk complex Kerr effect due to orbital
polarization, and $\kappa=(4\pi i\bar{n}\cos\psi)/\lambda$, with wavelength
$\lambda$, $\cos\psi=(1-\sin^{2}\phi_{i}/\bar{n}^{2})^{1/2}$, with $\bar{n}$
the complex index of refraction, $\phi_{i}$ the angle of incidence, and
$\kappa^{\pm}=\kappa\pm 1/l_{o}$. A similar expression exists for $\theta_{\rm
K}^{s}$ [23], but containing the spin diffusion length $l_{s}$, the SHE
conductivity $\sigma_{zx}^{\rm SH}$, and bulk Kerr effect due to spin
polarization, $\Phi_{\rm K}^{bulk,s}$.

Our ab initio calculations provide values $\sigma_{zx}^{\rm SH}$ and
$\sigma_{zx}^{\rm OH}$ as well as $\Phi_{\rm K}^{bulk,s}$ and $\Phi_{\rm
K}^{bulk,o}$. $\bar{n}$ is evaluated from the computed optical conductivity
$\sigma_{xx}(\omega)$, which is in good agreement with measurements (SM:S10).
Using both ab initio calculated quantities and quantities that are given from
the experiment ($\rho(t)$, $t$, $j$, $\phi_{i}$) we can use the measured MOKE
as function of thickness together with Eq. (1) to obtain values for $l_{o}$
and $l_{s}$. At this point it is instructive to mention that the calculated
[24] $\sigma_{zx}^{\rm OH}$ is $7000\,(\hbar/e)$ $\Omega^{-1}\rm{cm}^{-1}$
whereas $\sigma_{zx}^{\rm SH}\approx-70\,(\hbar/e)$ $\Omega^{-1}\rm{cm}^{-1}$,
i.e., much smaller. As a consequence, we find that the measured thickness-
dependent $\theta_{\rm K}$ data in Fig. 5 cannot be described by spin
diffusion. This provides our fifth evidence that the observed MOKE signal is
due to orbital accumulation.

Figure 5: Measured Kerr rotation $\theta_{\rm K}$ as function of Cr film
thickness (symbols), for a current density of $10^{11}$ Am-2. The full curve
is a fit to the measured MOKE rotation based on Eq. (1) with ab initio
calculated quantities to obtain the orbital diffusion length $l_{o}$.

To extract the orbital diffusion length from the measured data, we adopt two
approaches. First, we can divide $\theta_{\rm K}^{o}$ in Eq. (1) by the
measured resistivity $\rho(t)^{2}$ (SM:S12). The data for $\theta_{\rm
K}^{o}/\rho^{2}$ as function of $t$ can then be fitted by using the
theoretical input values to obtain $l_{o}$. Second, we can model the
thickness-dependent resistivity by first fitting the resistivity data
$\rho(t)$ with the Fuchs-Sondheimer equation [31, 32] (SM:S13). After that, we
use the obtained resistivities as input values in Eq. (1) and fit $\theta_{\rm
K}^{o}$ as a function of $t$ to obtain $l_{o}$. In Fig. 5 we show the result
from the fit of $\theta_{\rm K}^{o}$ with Eq. (1) (second procedure) together
with the experimental data. The theoretical curve based on ab initio
calculated values for $\sigma_{zx}^{\rm OH}$, $\kappa$, and $\Phi_{\rm
K}^{bulk,o}$ captures the measured MOKE data very well. The fit value of the
orbital diffusion length is $l_{o}=6.66\pm 0.48$ nm. Applying the first
procedure, we obtain $l_{o}=6.48\pm 0.62$ nm (SM:S14). Both approaches give
thus consistently very similar values, from which we find an averaged orbital
diffusion length $l_{o}=6.6\pm 0.6$ nm.

Using the theoretically calculated value of the orbital Hall conductivity and
the average resistivity of the films (SM:S12), we can estimate an orbital Hall
angle in our Cr films, $\theta^{\rm OH}=\sigma_{zx}^{\rm
OH}\rho\,(2e/\hbar)\approx 0.28$. This value is comparable to values of the
effective spin Hall angle in Pt/FM bilayers [6, 33].

Orbital transport is currently only poorly understood. A key unsolved question
is the size of the orbital diffusion length. So far, a rather wide range of
orbital diffusion lengths (from 3 nm to 74 nm) have been reported [13, 25, 28,
34, 17]. Liu et al. [34] obtained $l_{o}\sim 3.1$ nm for Nb from torque
measurements on Co/Nb bilayers. Fukunaga et al. [17] extracted from torque
measurements on Ni/Zr an $l_{o}\approx 9.7$ nm for Zr. Using MOKE as well as
torque measurements, Choi et al. [25] determined an orbital diffusion length
of $61-74$ nm for Ti. However, they had to reduce the ab initio computed
$\sigma^{\rm OH}$ by a factor of hundred to match it with their MOKE
measurements. From torque measurements on a Ni/Ti bilayer, Hayashi et al. [28]
deduced $l_{o}\approx 50$ nm for Ti. Both these values are roughly ten times
larger than our $l_{o}$ for Cr. Interestingly, Lee et al. [13] performed
torque measurements on Ni/Cr bilayers and obtained $l_{o}=6.1\pm 1.7$ nm for
Cr. This is remarkably consistent with the $l_{o}$ value we obtained using
magneto-optics on a single Cr layer, which arguably should provide more
accurate, intrinsic $l_{o}$ values. It should be emphasized that the orbital
diffusion length of Cr is longer than the Cr spin diffusion length of 4.5 nm,
determined on Fe/Cr films [35]. Consequently, long-range orbital transport
carried by a huge orbital conductivity is expected for Cr.

In conclusion, we have used MOKE microscopy combined with ab initio
calculations of MOKE and orbital Hall conductivity to measure the OHE induced
orbital accumulation in Cr thin films. Our results conclusively demonstrate
the ability to detect and quantify the OHE in a single non-magnetic layer and
provide direct evidence of a large current-induced orbital accumulation. This
opens a way to study the orbital Hall effect in metals free of effects caused
by an adjacent ferromagnet, which could bring new insights into the OHE, the
SHE, and the connection between them.

We thank Shuyu Cheng, Katherine Robinson, Daniel Russell, and Fengyuan Yang
for the growth of Pt samples. We further thank Pavel Novák, Ján Rusz, and
Markus Weißenhofer for helpful comments. I.L. and R.K.K. acknowledge support
from the Center for Emergent Materials, an NSF MRSEC, under award number
DMR-2011876. This work was furthermore supported by the Swedish Research
Council (VR), the Swedish National Infrastructure for Computing (SNIC) (Grant
No. 2018-05973), and the K. and A. Wallenberg Foundation (Grant No.
2022.0079).

## References

  * D’Yakonov and Perel’ [1971] M. I. D’Yakonov and V. I. Perel’, Possibility of Orienting Electron Spins with Current, JETP Lett. 13, 467 (1971).
  * Hirsch [1999] J. E. Hirsch, Spin Hall Effect, Phys. Rev. Lett. 83, 1834 (1999).
  * Kato  _et al._ [2004] Y. K. Kato, R. C. Myers, A. C. Gossard, and D. D. Awschalom, Observation of the Spin Hall Effect in Semiconductors, Science 306, 1910 (2004).
  * Valenzuela and Tinkham [2006] S. O. Valenzuela and M. Tinkham, Direct electronic measurement of the spin Hall effect, Nature 442, 176 (2006).
  * Sinova  _et al._ [2015] J. Sinova, S. O. Valenzuela, J. Wunderlich, C. H. Back, and T. Jungwirth, Spin Hall effects, Rev. Mod. Phys. 87, 1213 (2015).
  * Manchon  _et al._ [2019] A. Manchon, J. Železný, I. M. Miron, T. Jungwirth, J. Sinova, A. Thiaville, K. Garello, and P. Gambardella, Current-induced spin-orbit torques in ferromagnetic and antiferromagnetic systems, Rev. Mod. Phys. 91, 035004 (2019).
  * Tanaka  _et al._ [2008] T. Tanaka, H. Kontani, M. Naito, T. Naito, D. S. Hirashima, K. Yamada, and J. Inoue, Intrinsic spin Hall effect and orbital Hall effect in $4d$ and $5d$ transition metals, Phys. Rev. B 77, 165117 (2008).
  * Kontani  _et al._ [2009] H. Kontani, T. Tanaka, D. S. Hirashima, K. Yamada, and J. Inoue, Giant Orbital Hall Effect in Transition Metals: Origin of Large Spin and Anomalous Hall Effects, Phys. Rev. Lett. 102, 016601 (2009).
  * Go  _et al._ [2018] D. Go, D. Jo, C. Kim, and H.-W. Lee, Intrinsic Spin and Orbital Hall Effects from Orbital Texture, Phys. Rev. Lett. 121, 086602 (2018).
  * Jo  _et al._ [2018] D. Jo, D. Go, and H.-W. Lee, Gigantic intrinsic orbital Hall effects in weakly spin-orbit coupled metals, Phys. Rev. B 98, 214405 (2018).
  * Bhowal and Satpathy [2020] S. Bhowal and S. Satpathy, Intrinsic orbital moment and prediction of a large orbital Hall effect in two-dimensional transition metal dichalcogenides, Phys. Rev. B 101, 121112 (2020).
  * Kim  _et al._ [2021] J. Kim, D. Go, H. Tsai, D. Jo, K. Kondou, H.-W. Lee, and Y. Otani, Nontrivial torque generation by orbital angular momentum injection in ferromagnetic-metal/Cu/Al2O3 trilayers, Phys. Rev. B 103, L020407 (2021).
  * Lee  _et al._ [2021a] S. Lee, M.-G. Kang, D. Go, D. Kim, J.-H. Kang, T. Lee, G.-H. Lee, J. Kang, N. J. Lee, Y. Mokrousov, S. Kim, K.-J. Kim, K.-J. Lee, and B.-G. Park, Efficient conversion of orbital Hall current to spin current for spin-orbit torque switching, Commun. Phys. 4, 234 (2021a).
  * Lee  _et al._ [2021b] D. Lee, D. Go, H.-J. Park, W. Jeong, H.-W. Ko, D. Yun, D. Jo, S. Lee, G. Go, J. H. Oh, K.-J. Kim, B.-G. Park, B.-C. Min, H. C. Koo, H.-W. Lee, O. Lee, and K.-J. Lee, Orbital torque in magnetic bilayers, Nat. Commun. 12, 6710 (2021b).
  * Ding  _et al._ [2022] S. Ding, Z. Liang, D. Go, C. Yun, M. Xue, Z. Liu, S. Becker, W. Yang, H. Du, C. Wang, Y. Yang, G. Jakob, M. Kläui, Y. Mokrousov, and J. Yang, Observation of the Orbital Rashba-Edelstein Magnetoresistance, Phys. Rev. Lett. 128, 067201 (2022).
  * Sala and Gambardella [2022] G. Sala and P. Gambardella, Giant orbital Hall effect and orbital-to-spin conversion in $3d$, $5d$, and $4f$ metallic heterostructures, Phys. Rev. Research 4, 033037 (2022).
  * Fukunaga  _et al._ [2023] R. Fukunaga, S. Haku, H. Hayashi, and K. Ando, Orbital torque originating from orbital Hall effect in Zr, Phys. Rev. Res. 5, 023054 (2023).
  * Baek  _et al._ [2018] S.-h. C. Baek, V. P. Amin, Y.-W. Oh, G. Go, S.-J. Lee, G.-H. Lee, K.-J. Kim, M. D. Stiles, B.-G. Park, and K.-J. Lee, Spin currents and spin–orbit torques in ferromagnetic trilayers, Nature Mater. 17, 509 (2018).
  * Amin  _et al._ [2018] V. Amin, J. Zemen, and M. Stiles, Interface-Generated Spin Currents, Phys. Rev. Lett. 121, 136805 (2018).
  * Wang  _et al._ [2019] W. Wang, T. Wang, V. P. Amin, Y. Wang, A. Radhakrishnan, A. Davidson, S. R. Allen, T. J. Silva, H. Ohldag, D. Balzar, B. L. Zink, P. M. Haney, J. Q. Xiao, D. G. Cahill, V. O. Lorenz, and X. Fan, Anomalous spin–orbit torques in magnetic single-layer films, Nat. Nanotechnol. 14, 819 (2019).
  * Amin  _et al._ [2019] V. P. Amin, J. Li, M. D. Stiles, and P. M. Haney, Intrinsic spin currents in ferromagnets, Phys. Rev. B 99, 220405 (2019).
  * Go  _et al._ [2020] D. Go, F. Freimuth, J.-P. Hanke, F. Xue, O. Gomonay, K.-J. Lee, S. Blügel, P. M. Haney, H.-W. Lee, and Y. Mokrousov, Theory of current-induced angular momentum transfer dynamics in spin-orbit coupled systems, Phys. Rev. Research 2, 033401 (2020).
  * Stamm  _et al._ [2017] C. Stamm, C. Murer, M. Berritta, J. Feng, M. Gabureac, P. M. Oppeneer, and P. Gambardella, Magneto-Optical Detection of the Spin Hall Effect in Pt and W Thin Films, Phys. Rev. Lett. 119, 087203 (2017).
  * Salemi and Oppeneer [2022] L. Salemi and P. M. Oppeneer, First-principles theory of intrinsic spin and orbital Hall and Nernst effects in metallic monoatomic crystals, Phys. Rev. Mater. 6, 095001 (2022).
  * Choi  _et al._ [2021] Y.-G. Choi, D. Jo, K.-H. Ko, D. Go, K.-H. Kim, H. G. Park, C. Kim, B.-C. Min, G.-M. Choi, and H.-W. Lee, Observation of the orbital Hall effect in a light metal Ti, arXiv:2109.14847 [cond-mat] (2021).
  * [26] See Supplemental Material for details on the sample growth, experimental setup, resistivity measurements, measurements of additional samples, theoretical calculations, discussion of Oersted field, and the fitting of the thickness dependence.
  * Salemi  _et al._ [2021] L. Salemi, M. Berritta, and P. M. Oppeneer, Quantitative comparison of electrically induced spin and orbital polarizations in heavy-metal/$3d$-metal bilayers, Phys. Rev. Mater. 5, 074407 (2021).
  * Hayashi  _et al._ [2023] H. Hayashi, D. Jo, D. Go, T. Gao, S. Haku, H.-W. Mokrousov, Yuriy Lee, and K. Ando, Observation of long-range orbital transport and giant orbital torque, Commun. Phys. 6, 32 (2023).
  * Go  _et al._ [2021] D. Go, D. Jo, K.-W. Kim, S. Lee, M.-G. Kang, B.-G. Park, S. Blügel, H.-W. Lee, and Y. Mokrousov, Long-Range Orbital Transport in Ferromagnets, arXiv:2106.07928 [cond-mat] (2021).
  * Oppeneer [2001] P. M. Oppeneer, Magneto-Optical Kerr Spectra, in  _Handbook of Magnetic Materials_ , Handbook of Magnetic Materials, Vol. 13, edited by K. H. J. Buschow (Elsevier, Amsterdam, 2001) pp. 229–422.
  * Fuchs [1938] K. Fuchs, The conductivity of thin metallic films according to the electron theory of metals, Proc. Cambridge Phil. Soc. 34, 100 (1938).
  * Sondheimer [1952] E. H. Sondheimer, The mean free path of electrons in metals, Adv. Phys. 1, 1 (1952).
  * Zhu  _et al._ [2021] L. Zhu, D. C. Ralph, and R. A. Buhrman, Maximizing spin-orbit torque generated by the spin Hall effect of Pt, Appl. Phys. Rev. 8, 10.1063/5.0059171 (2021).
  * Liu  _et al._ [2023] F. Liu, B. Liang, J. Xu, C. Jia, and C. Jiang, Giant efficiency of long-range orbital torque in Co/Nb bilayers, Phys. Rev. B 107, 054404 (2023).
  * Bass and Pratt [2007] J. Bass and W. P. Pratt, Spin-diffusion lengths in metals and alloys, and spin-flipping at metal/metal interfaces: an experimentalist’s critical review, J. Phys.: Condens. Matter 19, 183201 (2007).


# Multi-fidelity information fusion with concatenated neural networks

Suraj Pawar
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Omer San
School of Mechanical & Aerospace Engineering,
Oklahoma State University,
Stillwater, Oklahoma - 74078, USA.
<EMAIL_ADDRESS>
&Prakash Vedula
School of Aerospace & Mechanical Engineering,
The University of Oklahoma,
Norman, OK 73019, USA.
&Adil Rasheed
Department of Engineering Cybernetics,
Norwegian University of Science and Technology,
7465 Trondheim, Norway.
Department of Mathematics and Cybernetics,
SINTEF Digital,
7034 Trondheim, Norway.
Trond Kvamsdal
Department of Mathematics and Cybernetics,
SINTEF Digital,
7034 Trondheim, Norway.
Department of Mathematical Sciences,
Norwegian University of Science and Technology,
7491 Trondheim, Norway

###### Abstract

Recently, computational modeling has shifted towards the use of statistical
inference, deep learning, and other data-driven modeling frameworks. Although
this shift in modeling holds promise in many applications like design
optimization and real-time control by lowering the computational burden,
training deep learning models need a huge amount of data. This big data is not
always available for scientific problems and leads to poorly generalizable
data-driven models. This gap can be furnished by leveraging information from
physics-based models. Exploiting prior knowledge about the problem at hand,
this study puts forth a concatenated neural network approach to build more
tailored, effective, and efficient machine learning models. For our analysis,
without losing its generalizability and modularity, we focus on the
development of predictive models for laminar and turbulent boundary layer
flows. In particular, we combine the self-similarity solution and power-law
velocity profile (low-fidelity models) with the noisy data obtained either
from experiments or computational fluid dynamics simulations (high-fidelity
models) through a concatenated neural network. We illustrate how the knowledge
from these simplified models results in reducing uncertainties associated with
deep learning models applied to boundary layer flow prediction problems. The
proposed multi-fidelity information fusion framework produces physically
consistent models that attempt to achieve better generalization than data-
driven models obtained purely based on data. While we demonstrate our
framework for a problem relevant to fluid mechanics, its workflow and
principles can be adopted for many scientific problems where empirical,
analytical, or simplified models are prevalent. In line with grand demands in
novel physics-guided machine learning principles, this work builds a bridge
between extensive physics-based theories and data-driven modeling paradigms
and paves the way for using hybrid physics and machine learning modeling
approaches for next-generation digital twin technologies.

_K_ eywords Multi-fidelity modeling, Physics-guided machine learning, Reduced-
order modeling

##  1 Introduction

The modeling of spatiotemporal dynamics of multiscale and multiphysics systems
is an open problem relevant to many scientific and engineering applications.
For instance, wind energy is a highly complex system whose dynamics is
governed by global atmospheric processes to turbulent boundary layer formed
around the blades that span nine orders of magnitudes [1]. Over the past
several decades, we have improved our understanding of such multiphysics
systems by developing accurate numerical models for governing equations of the
system, such as Navier-Stokes equations for fluid flows. However, these
numerical models can be computationally prohibitive, especially for nonlinear
multiscale systems, and their use in real-time optimization and control is
scarce. The recent advancement in machine learning (ML) and deep learning (DL)
holds the great potential for tackling the challenge of modeling and analysis
of high-dimensional systems and has been successful in diverse applications,
such as fluid mechanics [2], earth science [3], and material science [4].
These advances have been driven by a vast amount of data generated from high-
fidelity numerical simulations, experimental and satellite measurements, and
computing power along with the emergence of effective and efficient algorithms
that can extract relevant patterns from the data. ML/DL techniques has been
successfully applied for turbulence closure modeling [5], super-resolution of
climate data [6], predicting clustered weather patterns [7], reduced-order
modeling [8], and many more.

ML/DL models are capable of providing insights from data, exploiting these
insights in building predictive tools, and continuously updating themselves as
the new streams of data get available. Despite these advantages, ML/DL
techniques lack interpretability and suffer from the curse of dimensionality.
The interpretability issue can be addressed by understanding the physical
implications of ML/DL models [9], and understanding the neural network
correlations discovered from the data [10, 11]. By the curse of
dimensionality, we mean that DL models are data-hungry in nature. For
instance, Bonavita and Laloyaux [12] showed that the amount of the training
data for nonlinear dynamical systems grows exponentially with the
dimensionality of the system. Furthermore, pure ML/DL models lead to poor
extrapolation/generalization, i.e., they fit the observations data very well,
but predictions may be poor and physically inconsistent for data beyond the
distribution of the training dataset. To this end, physics-informed learning
algorithms that leverage prior knowledge based on the physical and
mathematical understanding of the system are proposed in several studies [13].

Figure 1: The proposed multi-fidelity information fusion framework. The
prediction from the low-fidelity models (e.g., self-similarity solution or law
of the wall) is concatenated with the latent variables at the certain hidden
layer of the neural network. An ensemble of neural networks is trained using
the negative log-likelihood loss function to estimate the uncertainty
associated with the prediction. Here, $\mathbf{x}$ refers to independent
variables or design parameters, or $\mathbf{y}$ indicates the quantities of
interest, and $h(\mathbf{x})$ represents a low-fidelity or simplified model
that is fast-to-compute.

There are two main techniques based on inductive biases, and learning biases
to embed physics into ML/DL models in combination with the observations data
[13]. Inductive biases techniques relate to building tailored ML/DL model
architecture that exploits the prior knowledge about the problem at hand to
build physically consistent data-driven models. The representative examples
includes embedding invariance property into neural network architecture [14,
15], imposing conservation laws of physical quantities or analytical
constraints into neural network [16, 17, 18], using prediction from simplified
models as the bias [19], and equivariant transformer networks [20]. The other
approach based on learning biases imposes constraints such as governing
equations in a soft manner by penalizing the loss function of ML/DL models.
Some of the examples of this approach pertain to physics-informed neural
networks [21], statistically constrained generative adversarial networks [22],
and Bayesian framework with auto-regressive model [23]. There is also a class
of hybrid analysis and modeling approaches that utilizes pure data-driven and
physics-based models in tandem [24]. While many methods have been demonstrated
to be successful in enforcing physics into ML/DL models, they offer many
possibilities to fuse domain knowledge to improve the generalizability and
data efficiency of data-driven models [25].

In this work, we introduce a tailored neural network architecture based on a
concatenation layer for surrogate modeling of high-fidelity data given the
solution from a computationally inexpensive low-fidelity model. There is a
hierarchy of numerical models ranging from simple empirical relations to
highly accurate numerical discretization-based models across all scientific
disciplines. For example, flow around airfoils can be modeled using tools
extending from panel methods [26] to the fully resolved direct numerical
simulation (DNS) [27]. Another example is the wake modeling of wind turbines,
where fast but typically inaccurate analytical models are adopted for tasks
like layout optimization and wind farm control [28]. However these models are
insufficient to take the unsteady nature of interactions of turbine wake with
other wakes as well as atmospheric turbulence into account, and such effects
can be modeled with computationally demanding but accurate models like large-
eddy simulation (LES) [29]. Our work draws inspiration from these multi-
fidelity modeling approaches and exploits the real-time prediction from low-
fidelity models to inform a DL model of high-fidelity observations. The
information fusion from the multi-fidelity sources of data leads to a robust
and generalizable surrogate model in comparison to purely data-driven model
trained on high-fidelity data. Figure 1 graphically illustrates the proposed
data fusion approach from multimodal data streams in the process of generating
physics-guided machine learning (PGML) models.

In a nutshell, our report puts forth a novel data-driven framework to take
prior knowledge about the system into account when generating a black-box deep
learning predictive model. How should we inject physics and domain knowledge
into machine learning models? How deep learning can be constructed as a
trustworthy approach toward more accurate real-time prediction of nonlinear
complex systems such as turbulent flows? In this paper, these are the
fundamental research questions that we tackle, and provide our insights.

We demonstrate the application of our framework for boundary layer flows.
Boundary layer phenomenon is one of the most important flows and is of
engineering concern in many scientific and industrial applications [30]. The
behavior of flow in the boundary layer has implications on the drag force in
ship hulls and aircraft, the energy required to move oil through pipes, and
the distribution of heat in the atmosphere, and therefore boundary layer flows
are extensively studied in the literature [31, 32, 33]. Given that boundary
layer flows are prevalent in engineering applications, building a
computationally efficient and accurate surrogate model is of paramount
importance for online tasks like boundary layer control to achieve lift
enhancement, noise mitigation, drag reduction, and wall cooling. Additionally,
boundary layer flows can be described using models that have different levels
of fidelity spanning analytical models to direct numerical simulations and
hence represent an interesting test case for illustrating the effectiveness of
concatenated neural networks. As we will detail in Methods section, our multi-
fidelity data-fusion framework approach consists of a concatenated neural
network architecture where the self-similarity solution is a low-fidelity
model and the Reynolds-Averaged Navier-Stokes Equations (RANSE) solver is a
high-fidelity model. This framework can be easily scaled to large problems
such as wake behind bluff bodies with information from many models fused to
predict the high-fidelity data. While in this work we consider only two levels
of fidelity, the proposed framework can be applied for blending information
from various levels of fidelity. Moreover, the neural architecture search
tools can be utilized to discover more complex and optimal architectures
automatically [34].

##  2 Results

In this section, we demonstrate the capability of the proposed approach to
reconstruct laminar and turbulent boundary layer flows around the flat plate.

Figure 2: Boundary layer prediction for laminar flat plate flow at $x/L=0.5$
along with the observations used for training the ML and PGML model. The
amount of observations data used for training the model is 10% (left), 30%
(middle), and 50% (right). The shaded area corresponds to two standard
deviation (2-SD) band.

###  2.1 Laminar flow past a flat plate

We refer to a simple feed-forward neural network as the machine learning (ML)
model and the concatenated neural network is called the physics-guided machine
learning (PGML) model. The ML model is trained purely based on the data, while
the PGML model uses the prediction from physics-based low-fidelity model
(Blasius equations, see Equation 10 and 11) along with the data. Figure 2
shows the profile of horizontal component of velocity versus distance from the
wall, at the location $x/L=0.5$, for different amounts of data used for
training the ML and PGML models. The velocity is normalized with the
freestream velocity and the vertical distance is normalized using the boundary
layer thickness for laminar flow over the flat plate. The ML model fails to
capture the accurate velocity profile when the velocity field information at
10% of locations within the computational domain is utilized for training. The
mean velocity profile predicted by the PGML model is highly accurate with just
10% of the observations. The predicted mean velocity profile is also
accompanied by a confidence interval spanning one standard deviation (SD) on
either side of the mean velocity profile and is an outcome of the uncertainty
quantification mechanism built into deep ensembles. Deep ensembles achieve
uncertainty quantification by training an ensemble of neural networks with the
negative log-likelihood loss function. The uncertainty estimate associated
with the PGML model is lower than the ML model and this is particularly
notable near the wall within the boundary layer, i.e., for $y/\delta<0.8$.

###  2.2 Turbulent flow past a flat plate

Next, we evaluate the performance of the proposed PGML approach for
reconstruction of turbulent boundary layer flow over a (smooth) flat plate.
The prior knowledge we concatenate in this case is the one-seventh power law
velocity profile (i.e., see Equation 13). Figure 3 displays the variation of
the normalized velocity profile in the vertical direction at $x/L=0.5$ for
different amounts of data used for training the ML and PGML models. We can
observe that the ML model performs very poorly for the data-sparse regime
(i.e., 5% and 10% of the observations). Such situations are very common in
scientific applications where a collection of high-fidelity data either from
experiments or numerical simulations can be prohibitive. The PGML model on the
other hand leads to an accurate prediction by exploiting the correlation
between low- and high-fidelity data. The prediction from the ML model is also
not reliable as indicated by the high width of the confidence band. The slope
of the boundary layer profile $\partial u/\partial y$ at the wall determines
the skin friction drag along the wall. This quantity is not predicted
accurately with the ML model and this can lead to poor estimation of the
quantity of interests like the total drag. The PGML model is successful in
predicting the correct slope of the velocity profile at the wall and therefore
will lead to a more accurate estimation of total drag.

Figures 4\- 6 shows the spatial variation of the predicted mean of the
velocity field, confidence interval of two standard deviations, and the error
with respect to the true velocity field near the wall region for 5%, 10%, and
30% of the training data, respectively. As the training data increases, the
error decreases for both ML and PGML models. The confidence estimate
associated with the PGML model is substantially higher (i.e., lower
uncertainty) than the ML model for all three datasets. Moreover, the error of
the PGML model is greatly reduced compared to the ML model. One other benefit
of constructing a PGML approach is its modular nature that can provide an
opportunity of bridging the gap between domain-specific knowledge and physics-
agnostic models.

In our previous numerical experiments, we focused on the reconstruction task
within the interpolation region. Both ML and PGML models were trained using
the data sampled from the whole domain, i.e., up to $L=2.0$, where $L$ is the
length of the flat plate. In our next numerical experiments, we sample
observations only from the region till $L=1.5$. Therefore, the region between
$L=1.5$ to $L=2.0$ corresponds to the extrapolation region. We quantify the
performance of the ML and PGML model using the variation of root mean squared
error (RMSE) percentage along the streamwise direction as follows

$\text{RMSE}(x)=100\times\bigg{(}\frac{1}{N_{y}}\sum_{j=1}^{N_{y}}\bigg{(}\frac{u_{T}(y_{j})-u_{P}(y_{j})}{u_{T}(y_{j})}\bigg{)}^{2}\bigg{)}^{1/2}$
(1)

where $u_{T}$ is the velocity of the high-fidelity model, $u_{p}$ is the
velocity predicted from the data-driven model, $N_{y}$ is the spatial
resolution in the wall-normal direction. From Figure 7, we can see that the
RMSE increases substantially in the extrapolation region for the ML model,
especially when the observations are very sparse, i.e., 5% of the data. This
is a well-known limitation of DL models to extrapolate poorly in the absence
of dense data. The PGML model on the other hand has RMSE almost one order of
magnitude less than the ML model in the interpolation region. Additionally,
the increase in RMSE is not significant in the extrapolation region. This
demonstrates that the PGML model is robust for the unseen condition, and it
performs well for out-of-distribution examples.

Figure 3: Boundary layer prediction for turbulent flat plate flow at $x/L=0.5$
along with the observations used for training the ML and PGML model. The
amount of observations data used for training the model is 5% (left), 10%
(middle), and 30% (right). The shaded area corresponds to two standard
deviation (2-SD) band. Figure 4: Prediction of the turbulent flat plate
boundary layer with 5% of data used for training the ML and PGML model. The
error is calculated as the difference between the true flow field and the flow
field predicted by ML and PGML models. Figure 5: Prediction of the turbulent
flat plate boundary layer with 10% of data used for training the ML and PGML
model. The error is calculated as the difference between the true flow field
and the flow field predicted by ML and PGML models. Figure 6: Prediction of
the turbulent flat plate boundary layer with 30% of data used for training the
ML and PGML model. The error is calculated as the difference between the true
flow field and the flow field predicted by ML and PGML models. Figure 7:
Variation of the normalized RMSE (in percentage) along the streamwise
direction for ML and PGML models. The amount of observations data used for
training the model is 5% (left), 10% (middle), and 30% (right). The gray color
shaded area represents the extrapolation region.

##  3 Discussion

This study aims to develop a physics-guided machine learning (PGML) framework
to improve data-driven models using prior knowledge from low-fidelity models.
The PGML is a new deep neural network architecture that makes it possible to
inject known physics during the training and _deployment_ processes to reduce
uncertainty and consequently improve the trade-off between efficiency and
accuracy. Our design of a hierarchically sequential learning algorithm allows
us to embed simplified theories, low order models, or empirical laws directly
into deep learning models. These physics-based injections assist the neural
network models in constraining the output to a manifold of the physically
consistent solution and leads to improved reliability and generalizability.
The PGML model trained using the deep ensembles algorithm provides us an
estimate of the uncertainty associated with the prediction. This uncertainty
information can be used for several applications like active learning, sensor
placement, and optimization. Some of the questions addressed in this study are
as follows

  * •

How prior information on the physics of the problem can be used to improve
black-box machine learning models?

  * •

Can a concatenated neural network architecture augmented by a simplified or
empirical model outperform a pure data-driven reconstruction model?

  * •

How can these data-driven predictive models be quantified regarding their
uncertainties?

  * •

What is the generalizability of predictive performance across ML and PGML
architectures, when these are applied to unseen conditions?

Our analysis indicates that an injection of the empirical relations like _one-
seventh power law_ improves the predictive models significantly for estimating
canonical turbulent flat plate boundary layer flows. We found that the PGML
outperforms its ML counterpart by reducing the RMSE to nearly an order of
magnitude lower levels. We also demonstrated that the proposed PGML framework
substantially reduces the model uncertainty even when only sparse observations
are available. Furthermore, generalizability of results is also supported when
we integrate our predictive models for unseen conditions. Specifically, the
RMSE distributions for the ML model (with respect to true data) are around 20%
and 30% for interpolation and extrapolation regions, respectively. However,
the PGML model has superior performance in both interpolation and
extrapolation regions with the RMSE distribution in the range of 3%. For
laminar flows, we also showed that a PGML approach (via injecting the Blasius
solution) outperforms the ML approach. The Blasius approximation can also be
extended for heat transfer problems (based on Falkner-Skan solutions) and for
strongly nonlinear problems using non-similarity solutions [35]. The
concatenated neural networks are capable of discovering a correlation between
low- and high-fidelity data allowing for smaller training dataset sizes,
training time, and improved extrapolation performance. Therefore, the PGML
model has a great potential for a vast number of physical systems where a
hierarchy of models is commonly used. We also highlight that the PGML approach
could be useful for generating physics consistent initial conditions to
accelerate large-scale high-fidelity computations with eliminating non-
physical initial transient time.

##  4 Methods

###  4.1 Multi-fidelity concatenated neural network

A neural network is a computational graph composed of several layers
consisting of the predefined number of neurons. Each neuron is associated with
certain coefficients called weights and some bias. The input from the previous
layer is multiplied by a weight matrix as shown below

$S^{l}=\mathbf{W}^{l}\mathbf{x}^{l-1},$ (2)

where $\mathbf{x}^{l-1}$ is the output of the $(l-1){\text{th}}$ layer,
$\mathbf{W}^{l}$ is the matrix of weights for the $l{\text{th}}$ layer. The
summation of the above input-weight product and the bias is then passed
through a node’s activation function which is usually some nonlinear function.
The introduction of nonlinearity through activation function allows the neural
network to learn highly complex relations between the input and output. The
output of the $l{\text{th}}$ layer can be written as

$\mathbf{x}^{l}=\zeta(S^{l}+B^{l}),$ (3)

where $B^{l}$ is the vector of biasing parameters for the $l{\text{th}}$ layer
and $\zeta$ is the activation function. If there are $L$ layers between the
input and the output in a neural network, then the neural network mapping
$\mathcal{F}:\mathbf{x}\rightarrow\mathbf{y}(\mathbf{x})$ can be represented
as follows

${\mathbf{y}}=\zeta_{L}(\cdot;\boldsymbol{\Theta}_{L})\circ\cdots\circ\zeta_{2}(\cdot;\boldsymbol{\Theta}_{2})\circ\zeta_{1}(\mathbf{x};\boldsymbol{\Theta}_{1})$
(4)

where $\boldsymbol{\Theta}$ represents the weight and bias of the
corresponding layer of the neural network and $\mathbf{y}$ is the final output
of the neural network. For the concatenated neural network, the information
from the low-fidelity model is injected at a certain intermediate layer of the
neural network as follows

${\mathbf{y}}=\zeta_{L}(\cdot;\boldsymbol{\Theta}_{L})\circ\cdots\circ\underbrace{\mathcal{C}(\zeta_{i}(\cdot;\boldsymbol{\Theta}_{i}),h(\mathbf{x}))}_{\text{Concatenation
layer}}\circ\cdots\circ\zeta_{1}(\mathbf{x};\boldsymbol{\Theta}_{1}),$ (5)

where $\mathcal{C}(\cdot,\cdot)$ represents the concatenation operation and
the information from the low-fidelity model, i.e., $h(\mathbf{x})$ is injected
at $i$th layer, and $\boldsymbol{\Theta}$ are the trainable parameters of the
corresponding layer. The concatenation operator takes the latent variables at
a particular layer and combines them with information from the low-fidelity
model to return a vector. Specifically, if the $i$th layer has $D_{h}$ neurons
and $h(\mathbf{x})\in\mathbb{R}^{D_{l}}$, the output of the $i$th layer will
be in $\mathbb{R}^{D_{h}+D_{l}}$. For the boundary layer reconstruction
problem, we use a neural network with three hidden layers and twenty neurons
in each hidden layer. The prediction from the low-fidelity model is
concatenated at the second hidden layer. One important caveat in a
concatenated neural network is the selection of an appropriate intermediate
layer at which to inject the low-fidelity model information, and tools like
AutoML can be applied to automate this search. The concatenation operator
given by Equation 5 can also be applied to all hidden layers simultaneously as

${\mathbf{y}}=\underbrace{\mathcal{C}(\zeta_{L}(\cdot;\boldsymbol{\Theta}_{L}),h(\mathbf{x}))\circ\cdots\circ\mathcal{C}(\zeta_{i}(\cdot;\boldsymbol{\Theta}_{i}),h(\mathbf{x}))\circ\cdots\circ\mathcal{C}(\zeta_{1}(\cdot;\boldsymbol{\Theta}_{1}),h(\mathbf{x}))}_{\text{Concatenation
layers}}.$ (6)

We highlight here that the proposed PGML framework is modular and Equation 6
can be generalized for fusing information from multiple low-fidelity models.
However, we dedicate this study to investigate the feasibility of the proposed
PGML framework for only two levels of approximations. An additional
optimization problem can be constructed to search the best architecture in
terms of relevant hyperparameters such as the number of hidden layers, and the
location and sparsity of the concatenation structure. Such auto PGML
investigations will be a topic that we will pursue in our future works.

###  4.2 Deep Ensembles: Training and prediction

Deep learning algorithms like neural networks approximate the mapping from
inputs to outputs using trainable parameters called weights and biases. The
parameters of the neural network are determined through the minimization of
the loss function. The prediction from the neural network is usually a point
estimate, i.e., continuous outputs for regression tasks and discrete classes
for classification problems. However, the information about the confidence in
model’s prediction might be crucial for many scientific applications [36]. The
uncertainty estimates can also be useful for applications like sensor
placement and Bayesian optimization. In this study, we apply the deep
ensembles algorithm for estimating the probabilistic distribution function
(PDF) of output conditioned on the inputs [37]. While there are state-of-the-
art methods like Bayesian neural networks [38] that quantifies uncertainty by
learning the distribution of weights, deep ensembles is adopted due to their
simplicity and scalability.

Here, we briefly discuss the uncertainty quantification mechanism of deep
ensembles. We assume that our training dataset $\mathcal{D}$ consists of $N$
samples $\mathcal{D}=\\{\mathbf{x}_{i},\mathbf{y}_{i}\\}_{i=1}^{N}$, where
$\mathbf{x}\in\mathbb{R}^{P}$ represents the $P-$dimensional features and the
label is $Q-$dimensional, i.e., $\mathbf{y}\in\mathbb{R}^{Q}$. A neural
network is trained by minimizing the loss function
$L(\mathbf{y},\tilde{\mathbf{y}}(\mathbf{x};\theta))$, where
$\tilde{\mathbf{y}}$ is the predicted label from a neural network
parameterized by $\mathbf{\theta}$. The most common loss function for
regression tasks is the mean squared error (MSE) between true and predicted
labels averaged over all samples in the dataset. The MSE loss function does
not give an estimate of the probability distribution of
$p(\tilde{\mathbf{y}}|\mathbf{x})$ and hence the uncertainty estimate is
usually absent with the prediction from neural networks.

In order to quantify the predictive uncertainty, the neural network is trained
to output the mean and variance of the Gaussian distribution in the output
layer. The weights of the neural network are determined by minimizing the
negative log-likelihood $\mathcal{L}$ as follows

$\mathbf{\theta}=\operatorname{arg\,min}_{\mathbf{\theta}}[\mathcal{L}],\quad\text{where}\quad\mathcal{L}=\sum_{i=1}^{N}\frac{1}{2}\text{log}\sigma^{2}(\mathbf{x}_{i})+\frac{(\mathbf{y}_{i}-\mu(\mathbf{x}_{i}))^{2}}{2\sigma^{2}(\mathbf{x}_{i})},$
(7)

where the mean $\mu$ and the variance $\sigma^{2}$ are parameterized by the
neural network. The positivity constraint is enforced for the variance by
passing the output corresponding to variance of distribution though the
softplus function $\text{log}(1+\text{exp}(\cdot))$, and adding a minimum
variance (for example $10^{-6}$) for numerical stability. If we assume the
variance to be constant in Equation 7 (i.e., it does not depend on input
features), then the negative log-likelihood loss function becomes analogous to
the MSE loss function. Therefore, from a probabilistic point of view,
minimizing the MSE is equivalent to minimizing negative log-likelihood with an
assumption of Gaussian distribution with constant standard deviation [39, 40].

The ensemble of neural networks has been demonstrated to be successful in
improving the predictive performance of machine learning models [41]. There
are broadly two methods of generating ensembles, (i) randomization-based
approaches where the ensembles can be trained in parallel without any
interaction, and (ii) boosting-based approaches where the ensembles are
trained sequentially [42]. The randomization procedure for generating
ensembles of neural networks should be such that prediction from individual
models are de-correlated and each individual models are strong (i.e., high
accuracy). In this work, the random initialization of weights of the neural
network is used for generating ensembles. There are other schemes such as
bagging where the ensembles of neural networks are trained on a different
subset of the original training data. However, random initialization is better
than bagging for improving predictive accuracy and uncertainty [37, 43]. This
simple and yet robust randomization approach is highly scalable as it allows
for distributed training of neural networks and can be applied to many
scientific problems. For computing the predictive probability distribution, we
approximate the ensemble prediction as a Gaussian whose mean and variance are
computed as follows

$\displaystyle\mu_{*}(\mathbf{x})$
$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}\mu_{\mathbf{\theta}_{j}}(\mathbf{x}),$
(8) $\displaystyle\sigma_{*}^{2}(\mathbf{x})$
$\displaystyle=\frac{1}{M}\sum_{j=1}^{M}(\sigma_{\mathbf{\theta}_{j}}^{2}(\mathbf{x})+\mu_{\mathbf{\theta}_{j}}^{2}(\mathbf{x}))-\mu_{*}^{2}(\mathbf{x}),$
(9)

where $\mu_{\mathbf{\theta}_{j}}$ and $\sigma_{\mathbf{\theta}_{j}}$ is the
mean and standard deviation of predicted probability distribution by the $j$th
neural network. We employ an ensemble of five neural networks in this study
(i.e., $M=5$).

###  4.3 Multi-fidelity data fusion for laminar boundary layer

The first test case considered in this study is the laminar boundary layer
flow. Boundary layer flows can be characterized by dividing the flow into two
regions, one inside the boundary layer where the viscosity dominates and one
outside the boundary layer where the effect of viscosity can be neglected. The
low-fidelity model considered for laminar flow is the steady-state two-
dimensional laminar boundary layer described using Blasius equation [44]. The
core idea behind Blasius equation is transforming a partial differential
equation (PDE) comprised of the flat plate boundary layer equations, with zero
pressure gradient, into a single ordinary differential equation (ODE) by using
a similarity solution approach. The derivation of the Blasius equation can be
found in many texts on fluid mechanics and we describe only the final form.
The Blasius equation and its boundary conditions can be written as

$\displaystyle 2f{{}^{\prime\prime\prime}}+ff^{\prime\prime}=0,$ (10)
$\displaystyle f(0)=f^{\prime}(0)=0,\ f^{\prime}(\infty)=1,$ (11)

where $f(\eta)$ is a function of similarity variable $\eta$. The similarity
variable $\eta$ is defined as $\eta=y\sqrt{u_{\infty}/{(x\nu)}}$, where $y$ is
the direction normal to the plate, $x$ is the direction along its length with
zero being the leading edge, $u_{\infty}$ is the freestream velocity, and
$\nu$ is the kinematic viscosity of the fluid. The third-order ODE is first
split into a coupled system of three first-order ODEs. Then we apply the
shooting method to determine the initial value for $f^{\prime\prime}(0)$, and
the first-order ODEs are numerically integrated with the fourth-order Runge-
Kutta scheme [45]. The velocity profile from the Blasius solution can be
determined using the relation $\bar{u}=u_{\infty}f^{\prime}$, where overbar
symbol is used to indicate the low-fidelity model estimate. The high-fidelity
observations are generated by solving the incompressible Navier-Stokes
equations with the PISO algorithm [46] available in OpenFoam. We get the
velocity (components along streamwise and wall-normal directions) and the
pressure distribution from CFD simulation. The Reynolds number based on the
length of the flat plate used for generating data is $\text{Re}_{L}=5\times
10^{4}$, where $L$ is the length of the flat plate. The training data for the
concatenated neural network is sampled randomly from the whole domain and the
velocity field is contaminated by adding a white Gaussian noise with zero mean
and a standard deviation of 0.05. Advanced sampling methods like Latin
hypercube sampling, clustered sampling can be utilized to reduce the number of
samples required for training and we will consider this as part of our future
work.

For the laminar boundary layer reconstruction task, the input to the neural
network is the location of the sensor, i.e., $\mathbf{x}=[x,y]$, where $x$ and
$y$ are the positions of the sensors in streamwise and wall-normal directions.
The output of the neural network is the probability distribution of $u$, $v$,
$p$ represented by their mean and standard deviation, where $u$ is the
velocity in the streamwise direction, $v$ is the velocity in the wall-normal
direction, and $p$ is the pressure at the sensor’s location. Additionally, the
velocity profile obtained from the Blasius solution is used as the low-
fidelity model, i.e., $h(\mathbf{x})=[\bar{u}]$. Following our previous
discussion, the problem formulation can be written as

$\\{\mu(\mathbf{x}),\sigma(\mathbf{x})\\}=\mathcal{F}(\mathbf{x},h(\mathbf{x})),\quad
p(\mathbf{y}|\mathbf{x},h(\mathbf{x}))=\mathcal{N}(\mu_{*}(\mathbf{x}),\sigma_{*}(\mathbf{x})).$
(12)

In a nutshell, Figure 8 shows the proposed multi-fidelity data-fusion
framework for laminar boundary layer flows.

Figure 8: Illustration of the multi-fidelity data-fusion framework applied to
the flat plate boundary layer prediction task. The self-similarity Blasius
solution is replaced by the one-seventh power law when applying to the
turbulent boundary layer flows.

###  4.4 Multi-fidelity data fusion for turbulent boundary layer

The boundary layer around the flat plate transitions to turbulence at high
Reynolds number. Before the advent of supercomputing, it was not possible to
numerically solve the Navier-Stokes equations for turbulent flows and fluid
dynamicists had to resort to experimental studies to derive empirical
relations for high Reynolds number flows. In this study, the low-fidelity
approximation for the turbulent boundary layer is obtained using the one-
seventh power law. The one-seventh power law [47] for computing the mean (or
ensemble-averaged) velocity profile for flat-plate turbulent boundary layer is
given as follows

$\frac{\bar{u}}{u_{\infty}}\approx\begin{cases}\big{(}\frac{y}{\delta}\big{)}^{1/7}\quad&\text{for}\quad
y\leq\delta,\\\ 1\quad&\text{for}\quad y>\delta,\end{cases}$ (13)

where $u_{\infty}$ is the freestream velocity, $y$ is the direction normal to
the plate, and the turbulent boundary layer thickness $\delta$ is computed as
follows [47]

$\delta\approx\frac{0.38x}{(\text{Re}_{x})^{1/5}},$ (14)

where $\text{Re}_{x}$ is the Reynolds number at a given $x$-location. There
are many such empirical relations available to approximate turbulent boundary
layers such as the log law, and Spalding’s law of the wall [48]. The high-
fidelity data is generated for the flat-plate turbulent boundary layer by
solving the incompressible Navier-Stokes equations with the SIMPLE algorithm
and $k-\omega$-SST turbulent model implemented in OpenFoam. The Reynolds
number based on the length of the flat plate for turbulent boundary layer
simulation is $\text{Re}_{L}=1\times 10^{7}$. One important parameter in
turbulence modeling is the dimensionless distance in the normal direction
called wall $y^{+}$ and is defined as $y^{+}=\sqrt{yu_{\tau}/\nu}$, where
$u_{\tau}$ is the friction velocity. The friction velocity is calculated based
on the wall shear stress as $u_{\tau}=\sqrt{\tau_{w}/\rho}$. The mesh is
refined near the flat plate in such a way that the near-wall $y^{+}$ is below
5. The locations for collecting the data are sampled randomly and are
contaminated by adding a white Gaussian noise with zero mean and a standard
deviation of 1.0 to mimic the measurement error. The formulation of turbulent
boundary layer reconstruction is similar to the laminar boundary
reconstruction as given in Equation 12 except for the low-fidelity model. The
low-fidelity model prediction for turbulent boundary layer flow is calculated
using Equation 13.

## Acknowledgements

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Advanced Scientific Computing Research under
Award Number DE-SC0019290. O.S. gratefully acknowledges their support.
Disclaimer: This report was prepared as an account of work sponsored by an
agency of the United States Government. Neither the United States Government
nor any agency thereof, nor any of their employees, makes any warranty,
express or implied, or assumes any legal liability or responsibility for the
accuracy, completeness, or usefulness of any information, apparatus, product,
or process disclosed, or represents that its use would not infringe privately
owned rights. Reference herein to any specific commercial product, process, or
service by trade name, trademark, manufacturer, or otherwise does not
necessarily constitute or imply its endorsement, recommendation, or favoring
by the United States Government or any agency thereof. The views and opinions
of authors expressed herein do not necessarily state or reflect those of the
United States Government or any agency thereof.

## References

  * [1] Paul Veers, Katherine Dykes, Eric Lantz, Stephan Barth, Carlo L Bottasso, Ola Carlson, Andrew Clifton, Johney Green, Peter Green, Hannele Holttinen, et al.  Grand challenges in the science of wind energy.  Science, 366(6464):eaau2027, 2019.
  * [2] Steven L Brunton, Bernd R Noack, and Petros Koumoutsakos.  Machine learning for fluid mechanics.  Annual Review of Fluid Mechanics, 52:477–508, 2020.
  * [3] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, et al.  Deep learning and process understanding for data-driven earth system science.  Nature, 566(7743):195–204, 2019.
  * [4] Jonathan Schmidt, Mário RG Marques, Silvana Botti, and Miguel AL Marques.  Recent advances and applications of machine learning in solid-state materials science.  NPJ Computational Materials, 5(1):1–36, 2019.
  * [5] Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao.  Turbulence modeling in the age of data.  Annual Review of Fluid Mechanics, 51:357–377, 2019.
  * [6] Karen Stengel, Andrew Glaws, Dylan Hettinger, and Ryan N King.  Adversarial super-resolution of climatological wind and solar data.  Proceedings of the National Academy of Sciences, 117(29):16805–16815, 2020.
  * [7] Ashesh Chattopadhyay, Pedram Hassanzadeh, and Saba Pasha.  Predicting clustered weather patterns: A test case for applications of convolutional neural networks to spatio-temporal climate data.  Scientific Reports, 10(1):1–13, 2020.
  * [8] Stefania Fresca, Luca Dede, and Andrea Manzoni.  A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes.  Journal of Scientific Computing, 87(2):1–36, 2021.
  * [9] Amy McGovern, Ryan Lagerquist, David John Gagne, G Eli Jergensen, Kimberly L Elmore, Cameron R Homeyer, and Travis Smith.  Making the black box more transparent: Understanding the physical implications of machine learning.  Bulletin of the American Meteorological Society, 100(11):2175–2199, 2019.
  * [10] Grégoire Montavon, Alexander Binder, Sebastian Lapuschkin, Wojciech Samek, and Klaus-Robert Müller.  Layer-wise relevance propagation: an overview.  Explainable AI: interpreting, explaining and visualizing deep learning, pages 193–209, 2019.
  * [11] Imme Ebert-Uphoff and Kyle Hilburn.  Evaluation, tuning, and interpretation of neural networks for working with images in meteorological applications.  Bulletin of the American Meteorological Society, 101(12):E2149–E2170, 2020.
  * [12] Massimo Bonavita and Patrick Laloyaux.  Machine learning for model error inference and correction.  Journal of Advances in Modeling Earth Systems, 12(12):e2020MS002232, 2020.
  * [13] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [14] Julia Ling, Andrew Kurzawski, and Jeremy Templeton.  Reynolds averaged turbulence modelling using deep neural networks with embedded invariance.  Journal of Fluid Mechanics, 807:155–166, 2016.
  * [15] Laure Zanna and Thomas Bolton.  Data-driven equation discovery of ocean mesoscale closures.  Geophysical Research Letters, 47(17):e2020GL088376, 2020.
  * [16] Arvind T Mohan, Nicholas Lubbers, Daniel Livescu, and Michael Chertkov.  Embedding hard physical constraints in neural network coarse-graining of 3d turbulence.  arXiv preprint arXiv:2002.00021, 2020.
  * [17] Tom Beucler, Michael Pritchard, Stephan Rasp, Jordan Ott, Pierre Baldi, and Pierre Gentine.  Enforcing analytic constraints in neural networks emulating physical systems.  Physical Review Letters, 126(9):098302, 2021.
  * [18] Sam Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  arXiv preprint arXiv:1906.01563, 2019.
  * [19] Suraj Pawar, Omer San, Burak Aksoylu, Adil Rasheed, and Trond Kvamsdal.  Physics guided machine learning using simplified theories.  Physics of Fluids, 33(1):011701, 2021.
  * [20] Kai Sheng Tai, Peter Bailis, and Gregory Valiant.  Equivariant transformer networks.  In International Conference on Machine Learning, pages 6086–6095. PMLR, 2019.
  * [21] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [22] Jin-Long Wu, Karthik Kashinath, Adrian Albert, Dragos Chirila, Heng Xiao, et al.  Enforcing statistical constraints in generative adversarial networks for modeling chaotic dynamical systems.  Journal of Computational Physics, 406:109209, 2020.
  * [23] Nicholas Geneva and Nicholas Zabaras.  Modeling the dynamics of pde systems with physics-constrained deep auto-regressive networks.  Journal of Computational Physics, 403:109056, 2020.
  * [24] Omer San, Adil Rasheed, and Trond Kvamsdal.  Hybrid analysis and modeling, eclecticism, and multifidelity computing toward digital twin revolution.  GAMM-Mitteilungen, page e202100007, 2021.
  * [25] K Kashinath, M Mustafa, A Albert, JL Wu, C Jiang, S Esmaeilzadeh, K Azizzadenesheli, R Wang, A Chattopadhyay, A Singh, et al.  Physics-informed machine learning: case studies for weather and climate modelling.  Philosophical Transactions of the Royal Society A, 379(2194):20200093, 2021.
  * [26] John L Hess.  Panel methods in computational fluid dynamics.  Annual Review of Fluid Mechanics, 22(1):255–274, 1990.
  * [27] Parviz Moin and Krishnan Mahesh.  Direct numerical simulation: a tool in turbulence research.  Annual Review of Fluid Mechanics, 30(1):539–578, 1998.
  * [28] Cristina L Archer, Ahmadreza Vasel-Be-Hagh, Chi Yan, Sicheng Wu, Yang Pan, Joseph F Brodie, and A Eoghan Maguire.  Review and evaluation of wake loss models for wind energy applications.  Applied Energy, 226:1187–1207, 2018.
  * [29] S-P Breton, J Sumner, Jens Nørkær Sørensen, Kurt Schaldemose Hansen, Sasan Sarmast, and Stefan Ivanell.  A survey of modelling methods for high-fidelity wind farm simulations using large eddy simulation.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 375(2091):20160097, 2017.
  * [30] Alexander J Smits and Ivan Marusic.  Wall-bounded turbulence.  Phys. Today, 66(9):25–30, 2013.
  * [31] Javier Jiménez.  Cascades in wall-bounded turbulence.  Annual Review of Fluid Mechanics, 44:27–45, 2012.
  * [32] Kiran Bhaganagar, John Kim, and Gary Coleman.  Effect of roughness on wall-bounded turbulence.  Flow, turbulence and combustion, 72(2):463–492, 2004.
  * [33] Xiaohua Wu and Parviz Moin.  Direct numerical simulation of turbulence in a nominally zero-pressure-gradient flat-plate boundary layer.  Journal of Fluid Mechanics, 630:5–41, 2009.
  * [34] Thomas Elsken, Jan Hendrik Metzen, and Frank Hutter.  Neural architecture search: A survey.  The Journal of Machine Learning Research, 20(1):1997–2017, 2019\.
  * [35] Shi-Jun Liao.  A general approach to get series solution of non-similarity boundary-layer flows.  Communications in Nonlinear Science and Numerical Simulation, 14(5):2144–2159, 2009.
  * [36] Dario Amodei, Chris Olah, Jacob Steinhardt, Paul Christiano, John Schulman, and Dan Mané.  Concrete problems in ai safety.  arXiv preprint arXiv:1606.06565, 2016.
  * [37] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 6405–6416, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [38] Radford M Neal.  Bayesian learning for neural networks, volume 118.  Springer Science & Business Media, 2012.
  * [39] Anthony Christopher Davison.  Statistical models, volume 11.  Cambridge University Press, Cambridge, 2003.
  * [40] Michael A Nielsen.  Neural networks and deep learning, volume 25.  Determination Press, San Francisco, 2015.
  * [41] Thomas G Dietterich.  Ensemble methods in machine learning.  In International workshop on multiple classifier systems, pages 1–15. Springer, Berlin, 2000.
  * [42] Artur J Ferreira and Mário AT Figueiredo.  Boosting algorithms: A review of methods, theory, and applications.  In Ensemble machine learning, pages 35–85. Springer, Berlin, 2012\.
  * [43] Stefan Lee, Senthil Purushwalkam, Michael Cogswell, David Crandall, and Dhruv Batra.  Why m heads are better than one: Training a diverse ensemble of deep networks.  arXiv preprint arXiv:1511.06314, 2015.
  * [44] Frank M White and Joseph Majdalani.  Viscous fluid flow, volume 3.  McGraw-Hill, New York, 2006.
  * [45] Parviz Moin.  Fundamentals of engineering numerical analysis.  Cambridge University Press, Cambridge, 2010.
  * [46] Raad I Issa.  Solution of the implicitly discretised fluid flow equations by operator-splitting.  Journal of Computational Physics, 62(1):40–65, 1986.
  * [47] Yunus A.. Çengel and John M.. Cimbala.  Fluid Mechanics: Fundamentals and Applications.  McGraw-Hill, New York, 2010.
  * [48] D. B. Spalding.  A Single Formula for the “Law of the Wall”.  Journal of Applied Mechanics, 28(3):455–458, 1961.