# High throughput screening for spin-gapless semiconductors in quaternary
Heusler compounds

Qiang Gao  Ingo Opahle  Hongbin Zhang<EMAIL_ADDRESS>Institute of
Materials Science, Technische Universit$\ddot{a}$t Darmstadt, 64287 Darmstadt,
Germany

###### Abstract

Based on high throughput density functional theory calculations, we performed
systematic screening for spin-gapless semiconductors (SGSs) in quaternary
Heusler alloys XX′YZ (X, X′, and Y are transition metal elements without Tc,
and Z is one of B, Al, Ga, In, Si, Ge, Sn, Pb, P, As, Sb, and Bi). Following
the empirical rule, we focused on compounds with 21, 26, or 28 valence
electrons, resulting in 12, 000 possible chemical compositions. After
systematically evaluating the thermodynamic, mechanical, and dynamical
stabilities, we successfully identified 70 stable SGSs, confirmed by explicit
electronic structure calculations with proper magnetic ground states. It is
demonstrated that all four types of SGSs can be realized, defined based on the
spin characters of the bands around the Fermi energy, and the type-II SGSs
show promising transport properties for spintronic applications. The effect of
spin-orbit coupling is investigated, resulting in large anisotropic
magnetoresistance and anomalous Nernst effects.

##  I Introduction

In recent years, spin-gapless semiconductors (SGSs) have drawn intensive
attention to the spintronics community. Basically, SGSs are half metals with
the majority spin channel being semimetallic, i.e., the gap is zero, while
there is a finite band gap in the minority spin channel. Based on how the spin
characters of bands touching the Fermi level, there are four types of SGSs as
sketched in Fig. 1(a). wang-2008-prl In the type-I SGSs, the valence band
maximum (VBM) and conduction band minimum (CBM) are in the same spin channel
while there is a gap in the opposite spin channel. This is the conventional
SGS as mentioned. Moreover, the CBM and VBM can hold different spin
characters, hereafter dubbed as type-II SGS. Additionally, if the VBM (CBM) is
of one spin character while the CBMs (VBMs) are originated from both spin
channels, type-III (type-IV) SGSs will be defined. In principle, the VBM and
CBM can touch each other at the same or different $k$-points, corresponding to
the direct or indirect zero band gap. In comparison to the usual half metals,
the 100% spin polarized carriers can be excited from the valence to conduction
bands with no energy cost, leading to new functionalities and potential
applications in logic gates. For instance, the spin polarized transport
properties of SGSs can be tuned by shifting the Fermi energy with finite gate
voltages, wang-2008-prl ; c-felsher-2013-prl which are promising for future
spintronic applications.

Based on first principles calculations, it was originally proposed that Co-
doped PbPdO2 can host the SGS state. wang-2008-prl However, its Curie
temperature (TC) is just about 180K, pbpdo2apl well below the room
temperature. The first above room temperature SGS was experimentally observed
in the inverse Heusler Mn2CoAl (TC=720 K). c-felsher-2013-prl Later on, the
Heusler compounds have been considered as outstanding candidates for SGSs. For
example, ternary Heusler Ti2MnAl, quaternary Heusler CoFeMnSi, and DO3 type
Heusler V3Al have been predicted theoretically to be SGSs skaftouros-2013-apl
; xu-2013-epl ; gao-2013-apl and also confirmed by experimental measurements.
feng-2015-pssppl ; bainsla-2015-prb ; jamer-2015-prb Interestingly, during the
explorations of SGSs in the Heusler compounds, an empirical rule has been
discovered. That is, the Heusler compounds with 18, 21, 26, or 28 valence
electrons are more probable to realize the SGS phase. skaftouros-2013-apl ;
qiang_2015 ; sg9 However, there has been no systematic study to design novel
SGS Heusler systems. Particularly, there are still a few questions about SGSs
to be understood. For instance, all four types of SGSs should in principle
exist but most experimentally studied systems are of type-I and type-II.
skaftouros-2013-apl ; xu-2013-epl ; qiang_2015 A particular intriguing
question is the effect of spin-orbit coupling (SOC) on the transport
properties of SGSs, i.e., whether a band gap can be opened with nontrivial
topological properties. Wang has proposed recently that SGSs are promising for
massless and dissipationless spintronics and quantum anomalous Hall effects.
nsr-wang In this regard, SGSs with direct band touching will be very
interesting, since they may host nontrivial topological properties after
considering SOC.

On the other hand, high throughput (HTP) screening based on density functional
theory (DFT) calculations has been proven to be an efficient way to search for
materials with desired properties. Stefano-Curtarolo-2012-comp-mat-sci ; Ingo-
Opahle-2013-njp Using the AFLOWLIB database, Carrete et al. have done HTP
calculations on approximately 79, 000 half-Heusler compounds and found 75
systems which are thermodynamically stable, where the thermal conductivities
and thermoelectric performance have also been evaluated. j-Carrete-2014-prx
The Heusler compounds with ten valence electrons (X2YZ , X = Ca, Sr, and Ba;
Y= Au; Z = Sn, Pb, As, Sb, and Bi) are demonstrated to have ultra-low lattice
thermal conductivities according to He’s HTP calculations. jiangang-
he-2016-prl Whereas, in a more recent HTP study, He et al. have identified 99
new nonmagnetic semiconductors following the 18 valence electron rule with
promising thermoelectric properties. he-18 Furthermore, for spintronics
applications, Ma et al. have performed a systemic HTP study on 405 inverse
Heusler alloys resulting in 14 stable semiconductors and 10 half metals. ma-
arxive Focusing on the magnetic properties, Sanvito et al. did HTP
calculations on 36, 540 Heusler alloys, leading to 248 thermodynamically
stable compounds with 20 magnetic cases. Sanvito Moreover, 21
antiferromagnetic Heusler compounds with high N$\acute{e}$el temperature have
been proposed for spintronic applications. Jan-2017 Last but not least, among
286 Heusler compounds, HTP screening calculations suggest 62% have a
tetragonal structure due to the peak-and-valley character in the density of
states. Parkin-phyrev-appl

In this work, we have carried out a systematic HTP screening for SGSs in
Heusler compounds (including DO3 binary, ternary, and quaternary Heusler
systems). Based on the empirical rule, we considered 12, 000 systems with 21,
26, or 28 valence electrons and identified 80 novel SGSs, which are
thermodynamically stable based on the formation energies. Among them, 70 are
both mechanically and dynamically stable. It is noted that the Heusler alloys
with 18 valence electrons are also promising to realize SGSs, which will be
investigated in the future. We have identified all four types of SGSs in the
quaternary Heusler compounds, together with one case showing direct band
touching at the Fermi energy. The longitudinal and transversal transport
properties were also evaluated based on the semi-classical transport theory,
revealing that SGSs are promising materials for spintronic applications. It is
demonstrated that the magnetization direction can be used to tailor the
electronic structure and hence the physical properties for SGSs with heavy
elements, due to the anisotropy caused by SOC.

##  II Computational Details

We considered quaternary Heusler compounds with a general chemical formula
XX′YZ, where X, X′, and Y are the transition metal elements except for the
radioactive Tc, and Z is one of the main group elements among B, Al, Ga, In,
Si, Ge, Sn, Pb, P, As, Sb, and Bi. For convenience, the ternary and binary
(DO3-type) Heusler systems are considered as quaternary Heusler by allowing X,
X′, or Y to be the same element. As shown in Fig. 1(b), quaternary Heusler
XX′YZ has the so-called LiMgPdSn-type structure with space group F$\bar{4}$3m
(space group 216), consisting of 4 Wyckoff positions 4a(0,0,0),
4c($\frac{1}{4}$,$\frac{1}{4}$,$\frac{1}{4}$),
4b($\frac{1}{2}$,$\frac{1}{2}$,$\frac{1}{2}$), and
4d($\frac{3}{4}$,$\frac{3}{4}$,$\frac{3}{4}$).heusler-1 ; heusler-2 According
to the empirical rule for the number of valence electrons (N${}_{\text{V}}$),
all the possible chemical composition with 21, 26, and 28 valence electrons
are generated, leading to about 12,000 possible compounds. Moreover, three
site occupations are considered for each chemical coposition, as shown in Fig.
1(c). heusler-stru Lastly, we consider that all the transition metal elements
(X, X′, and Y) are magnetic while the main group element (Z) is non-magnetic
(NM). For each chemical composition in each site-occupation, we consider five
spin configurations, namely, the NM, FM, AF1, AF2, and AF3 phases (Fig. 1(d)).

Figure 1: (Color online) (a) Sketches of the density of states for the four
types of SGSs, defined based on the touching schemes of the majority (marked
in blue) and minority (marked in red) bands. (b) The crystal structure of
quaternary Heusler XX′YZ, where the solid (dashed) lines indicate the
conventional cubic (primitive rhombohedral) cell. (c) The three possible site
occupations for a quaternary Heusler alloy with a specific chemical
composition. (d) The possible spin configurations within the primitive cell.
(e) The work flow for the present HTP screening.

The HTP screening has been carried out in an automated way following the work
flow shown in Fig. 1(e), managed with our in-house developed high-throughput
environment (HTE). Ingo-Opahle-2013-njp ; ingo-highthrough The DFT
calculations are performed using the Vienna ab initio Simulation Package
(VASP). vasp1 ; vasp2 For each composition-occupation case, the structural
relaxation is done in a two-step manner to save computational time. In the
first step, ultrasoft pseudopotentials (US-PP) chaoyan are used in combination
with the PW91 vasp3 exchange correlation functional, where the cutoff energy
for the plane wave basis is set to 250 eV and and a k-space density of 40 Å-1.
The follow-up finer relaxation is done using the projector augmented plane
wave (PAW) method with the exchange-correlation functional under the
generalized gradient approximation (GGA) parameterized by Perdew, Burke, and
Ernzerhof (PBE). vasp4 The cutoff energy for the plane wave expansion is
increased to 350 eV and the $k$-mesh density is increased to 50 Å-1 to achieve
good convergence. The structural relaxations are done for each magnetic
configuration mentioned above.

After obtaining the magnetic ground state together with the optimized
crystalline structures, the formation energy ($\bigtriangleup H$) is evaluated
to verify the thermodynamic stability, i.e., the stability with respect to
decomposing into constituting elements. For a general quaternary Heusler
XX′YZ, the formation energy is expressed as

$\bigtriangleup
H^{XX^{\prime}YZ}=\frac{1}{4}[E^{XX^{\prime}YZ}-(E^{X}+E^{X^{\prime}}+E^{Y}+E^{Z})],$
(1)

where E${}^{\text{X}}$, E${}^{\text{X'}}$, E${}^{\text{Y}}$ and
E${}^{\text{Z}}$ are the energies of elements X, X′, Y, and Z in their bulk
forms, while E${}^{XX^{\prime}YZ}$ is the ground state energy of XX′YZ. The
electronic structure together with the magnetic moments of the compounds with
negative formation energies are calculated with a denser $k$-mesh of $21\times
21\times 21$ using the full-potential local-orbital minimum-basis band
structure scheme (FPLO). fplo1 ; fplo2 The SGS phase can be identified by
examining the value of magnetic moments (i.e., being an integer following the
Slater-Pauling rule as discussed below) and the band structure directly.

For the candidate SGSs, we further checked the mechanical and dynamical
stability. The mechanical stability describes the stability of the crystal
against deformations or distortions in terms of strain, which can be obtained
based on the elastic constants (Cij). Basically, the elastic constants are
associated with the second-order change of the internal energy for a crystal
under an arbitrary deformation of strain as

$C_{ij}=\frac{1}{V_{0}}(\frac{\partial^{2}E}{\partial\varepsilon_{i}\partial\varepsilon_{j}}),$
(2)

where E is the internal energy, V0 is the equilibrium volume of the crystal,
and $\varepsilon_{i}$ or $\varepsilon_{j}$ denote applied strains. For a cubic
crystal system (such as Heusler compounds in this work), the elastic constant
matrix has only three independent elements as

$C_{cubic=}{\left[\begin{array}[]{cccccc}C_{11}&C_{12}&C_{12}&0&0&0\\\
C_{12}&C_{11}&C_{12}&0&0&0\\\ C_{12}&C_{12}&C_{11}&0&0&0\\\
0&0&0&C_{44}&0&0\\\ 0&0&0&0&C_{44}&0\\\ 0&0&0&0&0&C_{44}\end{array}\right]}.$
(3)

Correspondingly, the Born stability conditions elstic-huang suggest

$C_{11}+2C_{12}>0,C_{11}-C_{12}>0,C_{44}>0,$ (4)

which are related the bulk, tetragonal, and shear moduli, respectively.

On the other hand, the dynamical stability describes the change of the total
energy with respect to the internal degrees of freedom, i.e., the atomic
displacements. In the harmonic approximation, the total energy of a crystal
can be expressed as in terms of displacements $D_{\mathbf{R}\sigma}$

$E=E_{0}+\frac{1}{2}\sum_{\mathbf{R},\sigma}\sum_{\mathbf{R}^{\prime},\sigma^{\prime}}D_{\mathbf{R}\sigma}\Phi^{\sigma\sigma^{\prime}}_{\mathbf{R}\mathbf{R^{\prime}}}D_{\mathbf{R}^{\prime}\sigma^{\prime}},$
(5)

where $\mathbf{R}$ is the position, $\sigma$ is the Cartesian index, and
$\Phi^{\sigma\sigma^{\prime}}_{\mathbf{R}\mathbf{R^{\prime}}}$ is the
interatomic force constant matrix. The dynamical stability is determined by
the dynamical matrix $D(\mathbf{q})$, which can be obtained from Fourier
transformation of $\Phi(\mathbf{R})$ as following

$D(\mathbf{q})=\frac{1}{M}\sum_{R}\Phi(\mathbf{R})e^{-i\mathbf{qR}}$ (6)

where $\mathbf{q}$ is the wave vector of phonon. Dynamical stability indicates
that $D(\mathbf{q})$ is positive-definite, meaning all the phonons have real
and positive frequencies $\omega(\mathbf{q})$. The phonon dispersion
calculations are carried out using the Phonopy phonopy package with force
constants obtained from VASP.

Finally, the transport properties are studied for a few representative SGS
candidates, including the anomalous Hall conductivity (AHC) and the
longitudinal conductivity. The AHC is calculated by integrating the Berry
curvature ($\Omega(\mathbf{k})$) over the whole Brillouin zone (BZ) as
$\sigma_{xy}=\frac{e^{2}}{\hbar}\int_{BZ}\Omega(\mathbf{k})d^{3}\mathbf{k}$,
with the Berry curvature given by

$\Omega_{xy}(\mathbf{k})=2\text{Im}\sum_{n<E_{F}}\sum_{m\neq
n}\frac{<\psi_{n{\mathbf{k}}}|\nu_{x}|\psi_{m{\mathbf{k}}}><\psi_{m{\mathbf{k}}}|\nu_{y}|\psi_{n{\mathbf{k}}}>}{(\epsilon_{m{\bf
k}}-\epsilon_{n{\mathbf{k}}})^{2}},$ (7)

where $\psi_{\alpha{\mathbf{k}}}$ is the spinor Bloch wave function
corresponding to the eigenenergy $\epsilon_{\alpha{\bf k}}$, and $\nu_{i}$ is
the $i$-th Cartesian component of the velocity operator. In our calculations,
in order to achieve numerical convergence, the AHC is obtained using the
Wannier interpolation technique based on the Maximally Localised Wannier
Functions (MLWFs). wannier90 Furthermore, the longitudinal conductivities at
finite temperature (300 K) for SGSs are calculate based on the semiclassical
theory with the BoltzTrap boltztrap code. Here the energy independent
relaxation time ($\tau$) is used to approximate the distribution function as

$(\frac{\partial f}{\partial t})_{s}=-\frac{f-f_{0}}{\tau}$ (8)

where $f_{0}$ and $f$ are the equilibrium and nonequilibrium distribution
functions, respectively. The conductivity is expressed by

$\sigma_{\alpha\beta}(T,\mu)=\frac{1}{V}\int\bar{\sigma}_{\alpha\beta}(\epsilon)[-\frac{\partial
f_{0}(T,\epsilon,\mu)}{\partial\epsilon}]d\epsilon$ (9)

where $\alpha$ and $\beta$ are the Cartesian indices, $V$ and $\mu$ indicate
the unit cell volume and the chemical potential, respectively. The transport
distribution function $\bar{\sigma}_{\alpha\beta}(\epsilon)$ can be evaluated
by

$\bar{\sigma}_{\alpha\beta}(\epsilon)=\frac{e^{2}}{N}\sum_{i,\mathbf{k}}\tau\cdot\nu_{\alpha}(i,\mathbf{k})\cdot\nu_{\beta}(i,\mathbf{k})\cdot\frac{\delta(\epsilon-\epsilon_{i,\mathbf{k}})}{d\epsilon}$
(10)
$\nu_{\alpha}(i,\mathbf{k})=\frac{1}{\hbar}\bigtriangledown_{k}\epsilon_{i,\mathbf{k}}$
(11)

where $\mathbf{k}$, i, and N are the wave vector, band index, and the number
of the sampled $\mathbf{k}$ points. For bonding analysis, the crystal orbital
Hamilton population (COHP) was evaluated using the LOBSTER code. cohp

##  III results and discussion

###  III.1 Validation

Table 1: Comparisons between our HTP calculations and previous reported SGSs. For the mechanical and dynamical stabilities, “1” (“0”) indicates the system is stable (unstable). “Ref. Exp.” and “Ref. Cal.” denote experimental and computational references. Compound | latt. |  Mtot | $\bigtriangleup H$ | Mechanical | Dynamical |
---|---|---|---|---|---|---
| (Å) |  ($\mu_{B}$) | (eV/at.) | stability | stability |
$N_{V}=21$ |  |  |  |  |  |
Ti2CoSi | 6.081 | 3.00 | -0.3718 | 0 | 1 |
Ref. Cal. skaftouros-2013-apl | 6.030 | 3.03 |  |  |  |
MnCrTiSi | 5.855 | 3.02 | -0.4103 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 5.860 | 2.98 |  |  |  |
MnCrVAl | 5.897 | 3.00 | -0.2110 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 5.900 | 2.99 |  |  |  |
MnVTiAs | 5.978 | 2.90 | -0.2353 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 5.990 | 2.87 |  |  |  |
CoVTiAl | 5.978 | 3.00 | -0.3248 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 6.040 | 3.00 |  |  |  |
FeVTiSi | 5.978 | 3.02 | -0.4351 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 5.910 | 2.99 |  |  |  |
FeCrTiAl | 5.964 | 3.02 | -0.2920 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 5.970 | 3.00 |  |  |  |
CoVHfGa | 6.193 | 2.95 | -0.2434 | 1 | 1 |
Ref. Cal. yang_2016 | 6.260 | 3.00 |  |  |  |
CrFeHfGa | 6.127 | 3.00 | -0.1858 | 1 | 1 |
Ref. Cal. yang_2016 | 6.261 | 3.02 |  |  |  |
ZrCoVIn | 6.445 | 2.97 | -0.0632 | 0 | 1 |
Ref. Cal. qiang_2015 | 6.468 | 3.00 | -0.3500 |  |  |
ZrFeCrIn | 6.408 | 3.02 | 0.0279 | 1 | 0 |
Ref. Cal. qiang_2015 | 6.419 | 3.00 | -0.0325 |  |  |
ZrFeCrGa | 6.177 | 3.00 | -0.1690 | 1 | 1 |
Ref. Cal. qiang_2015 | 6.184 | 3.00 | -0.2400 |  |  |
ZrFeVGe | 6.199 | 3.06 | -0.2069 | 1 | 0 |
Ref. Cal. qiang_2015 | 6.210 | 3.00 | -0.2500 |  |  |
$N_{V}=26$ |  |  |  |  |  |
Mn2CoAl | 5.729 | 2.01 | -0.2666 | 1 | 1 |
Ref. Exp. c-felsher-2013-prl | 5.798 | 2.00 |  |  |  |
Ref. Cal. xing_2008 | 5.760 | 2.00 |  |  |  |
CoFeCrAl | 5.692 | 2.00 | -0.1931 | 1 | 1 |
Ref. Exp. Nehra | 5.736 | 2.00 |  |  |  |
Ref. Cal. gao_jac | 5.710 | 2.00 | -0.2500 |  |  |
CoFeCrGa | 5.717 | 2.00 | -0.0686 | 1 | 1 |
Ref. Exp. Bainsla_prb_2015_1 | 5.736 | 2.00 |  |  |  |
Ref. Cal. gao_jac | 5.730 | 2.00 |  |  |  |
CoFeTiAs | 5.835 | 2.00 | -0.3615 | 1 | 1 |
Ref. Cal. kog_jap_2013 | 5.850 | 1.99 |  |  |  |
CoMnCrSi | 5.669 | 2.00 | -0.3280 | 1 | 1 |
Ref. Cal. xu-2013-epl | 5.630 | 2.00 | -0.3750 |  |  |
FeMnCrSb | 6.059 | 2.00 | 0.0996 | 1 | 1 |
Ref. Cal. xu-2013-epl | 5.980 | 2.00 |  |  |  |
ZrCoFeP | 5.941 | 2.00 | -0.3491 | 0 | 0 |
Ref. Cal. qiang_2015 | 5.944 | 2.00 | -0.6500 |  |  |
$N_{V}=28$ |  |  |  |  |  |
CoFeMnSi | 5.597 | 4.00 | -0.3833 | 1 | 1 |
Ref. Exp. bainsla-2015-prb | 5.658 | 4.00 |  |  |  |
Ref. Cal.dai_jap_2009 | 5.609 | 4.00 |  |  |  |
Mn2CuAl | 5.710 | 0.00 | -0.1066 | 0 | 0 |
Ref. Cal. luo_2016_jac | 5.650 | 0.00 |  |  |  |
Cr2ZnSi | 5.972 | 0.00 | 0.08745 | 1 | 0 |
Ref. Cal. zhang_epl | 5.850 | 0.00 |  |  |  |
Cr2ZnGe | 6.123 | 0.00 | 0.1898 | 1 | 1 |
Ref. Cal. zhang_epl | 6.140 | 0.22 |  |  |  |
Cr2ZnSn | 6.413 | 0.00 | 0.3079 | 1 | 0 |
Ref. Cal. zhang_epl | 6.530 | 0.14 |  |  |  |

To validate the stability criteria implemented in the HTE, we collected all
the previously reported Heusler SGSs, and compared with our DFT results (TABLE
LABEL:table:literature). The lattice constants, total magnetic moments (TABLE
LABEL:table:literature), and the electronic structure (not shown) are in good
agreement with the literature. However, even though the formation energies for
most of the reported Heusler SGSs are negative, ZrFeCrIn and Cr2ZnX (X = Si,
Ge, and Sn) turn out to be thermodynamically unstable in our HTP calculations.
For ZrFeCrIn, in the previous calculations, qiang_2015 the energies of
composite elements with the fcc structure are considered, which leads to an
underestimation of the formation energy. This explains also the big difference
for the formation energy of ZrCoVIn. For the Cr2ZnX compounds, only the
inverse Heusler structure is considered in Ref. zhang_epl, . According to our
calculations, all three compounds end up with the non-magnetic (NM) full
Heusler structure, with positive formation energies. That is, even though the
electronic structure might be interesting with the hypothetical crystal
structures, the stability should be checked before making valid predictions.
It is noted that for quaternary compounds, the thermodynamical stability with
respect to other competing binary, ternary, and quaternary phases, i.e., the
distance to the convex hull should also be evaluated. We note that 55
previously unknown, thermodynamically stable (low convex hull) quaternary
Heusler compounds are discovered among 2,000,000 compounds by using machine
learning method. apsmeeting

Furthermore, it is observed that the mechanical stability or the dynamical
stability criteria are also critical for some previously predicted compounds.
For instance, according to our calculations, Ti2CoSi and ZrCoVIn are
mechanically unstable, ZrFeVGe is dynamically unstable, and ZrCoFeP and
Mn2CuAl are both mechanically and dynamically unstable. We note that such
compounds may still be synthesized experimentally using molecular beam
epitaxy, which is known to be efficient in obtaining metastable crystalline
phases. For all the systems which have been experimentally synthesized, such
as CoFeCrAl, CoFeCrGa, CoFeMnSi, and Mn2CoAl, we observed that they fulfill
all three stability criteria based on our calculations. This confirms the
reliability of our theoretical framework to do HTP screening for novel SGSs.

Table 2: Basic information of the newly predicted SGS candidates with negative formation energies in our research. The red compounds are either mechanically or dynamically unstable.  XX′YZ |  aopt |  Mtot | $\bigtriangleup H$ | SGS |  |  XX′YZ |  aopt |  Mtot | $\bigtriangleup H$ | SGS |
---|---|---|---|---|---|---|---|---|---|---|---
(4a,4b,4c,4d) | (Å) |  ($\mu_{B}$) | (eV/atom) | type |  | (4a,4b,4c,4d) | (Å) |  ($\mu_{B}$) | (eV/atom) | type |
N${}_{\text{V}}$=21 |  |  N${}_{\text{V}}$=26 |
IrVYSn | 6.720 | 3.00 | -0.0942 | SOC-I |  | CoOsTiSb | 6.255 | 2.00 | -0.1635 | I |
CoVYSn | 6.620 | 3.00 | -0.0862 | II |  | CoFeHfSb | 6.232 | 2.00 | -0.2847 | I |
CoVScSn | 6.402 | 3.00 | -0.2049 | III |  | CoOsZrSb | 6.453 | 2.00 | -0.1075 | I |
IrVScSn | 6.518 | 3.00 | -0.2488 | SOC-II |  | RhFeTiSb | 6.259 | 1.95 | -0.3896 | I |
RhVScSn | 6.518 | 3.00 | -0.2773 | I |  | CoFeTiSb | 6.074 | 2.00 | -0.2948 | I |
CoVYGe | 6.377 | 3.00 | -0.0763 | II |  | IrFeTiSb | 6.287 | 1.99 | -0.2932 | III |
CoVScGe | 6.145 | 3.00 | -0.2749 | II |  | CoRuTiSb | 6.228 | 2.00 | -0.3261 | I |
IrVScGe | 6.300 | 3.00 | -0.3025 | II |  | CoFeNbGe | 5.961 | 2.00 | -0.2374 | I |
RhVScGe | 6.290 | 3.00 | -0.3318 | II |  |  CoOsNbSn | 6.352 | 2.00 | -0.0609 | I |
RhVYGe | 6.512 | 3.00 | -0.1377 | III |  | CoRuTaSn | 6.303 | 2.00 | -0.1268 | I |
CoVYSi | 6.297 | 3.00 | -0.1077 | II |  | IrFeTaSn | 6.354 | 1.98 | -0.1782 | I |
CoVScSi | 6.058 | 3.00 | -0.355 | II |  | CoOsTaGe | 6.143 | 2.00 | -0.0702 | I |
IrVScSi | 6.215 | 3.00 | -0.4254 | SOC-II |  | CoOsTaSi | 6.064 | 1.99 | -0.2546 | I |
RhVScSi | 6.210 | 3.00 | -0.4242 | II |  | CoOsTaSn | 6.332 | 2.00 | -0.007 | I |
RhVYSi | 6.438 | 3.00 | -0.1862 | III |  | CoFeTaGe | 5.938 | 2.00 | -0.2475 | I |
PtVScAl | 6.369 | 3.00 | -0.4431 | SOC-I |  | CoFeTaSi | 5.856 | 2.00 | -0.4222 | I |
PtVYAl | 6.608 | 3.00 | -0.2477 | I |  | CoFeTaSn | 6.154 | 2.00 | -0.1353 | I |
PtVYGa | 6.600 | 3.00 | -0.1867 | I |  | IrCoNbAl | 6.162 | 1.99 | -0.5563 | I |
FeCrHfAl | 6.142 | 3.00 | -0.2456 | II |  | IrCoNbGa | 6.173 | 2.00 | -0.4043 | I |
OsCrHfAl | 6.299 | 3.00 | -0.403 | II |  | IrCoNbIn | 6.360 | 2.00 | -0.1326 | I |
RuCrHfAl | 6.284 | 3.00 | -0.4544 | II |  | IrCoTaAl | 6.140 | 2.00 | -0.5579 | I |
FeCrTiAl | 5.964 | 3.00 | -0.292 | II |  | IrCoTaGa | 6.150 | 2.00 | -0.388 | I |
FeCrZrAl | 6.194 | 3.00 | -0.2156 | III |  | IrCoTaIn | 6.336 | 2.00 | -0.1622 | I |
OsCrZrAl | 6.347 | 3.00 | -0.3543 | SOC-II |  | CoCoNbAl | 5.970 | 2.00 | -0.4312 | I |
RuCrZrAl | 6.335 | 3.00 | -0.4154 | III |  | CoCoNbGa | 5.968 | 2.00 | -0.3299 | I |
FeCrScSi | 5.992 | 3.00 | -0.279 | II |  | CoCoNbIn | 6.179 | 2.00 | -0.0869 | I |
FeCrScSn | 6.364 | 3.00 | -0.0891 | II |  | IrCoTiPb | 6.380 | 2.00 | -0.0571 | I |
FeCrYSi | 6.236 | 3.00 | -0.0081 | III |  | IrCoTiSn | 6.276 | 2.00 | -0.3789 | I |
OsCrYSi | 6.386 | 3.00 | -0.0246 | SOC-III |  | IrCoTiSi | 5.965 | 2.00 | -0.6805 | I |
CoVHfAl | 6.211 | 3.00 | -0.2896 | I |  | CoRuCrAl | 5.848 | 2.01 | -0.2802 | II |
IrVHfAl | 6.346 | 3.00 | -0.4634 | II |  | NiCrMnAl | 5.809 | 2.00 | -0.2127 | III |
RhVHfAl | 6.342 | 3.00 | -0.3855 | II |  | NiReCrAl | 5.920 | 1.97 | -0.1633 | II |
CoVZrAl | 6.258 | 3.00 | -0.2662 | I |  | CoOsCrAl | 5.866 | 2.00 | -0.2412 | II |
CoVZrGa | 6.238 | 3.00 | -0.2317 | I |  |  N${}_{\text{V}}$=28 |
IrTiZrSn | 6.651 | 2.98 | -0.3335 | II |  | NiFeMnAl | 5.731 | 4.00 | -0.2773 | IV |
IrTiZrSi | 6.385 | 2.96 | -0.4232 | II |  |  Continue with $N_{V}=21$ |
FeVNbAl | 6.117 | 2.99 | -0.2012 | II |  | MnCrNbAl | 6.077 | 3.00 | -0.1912 | II |
FeVTaAl | 6.097 | 2.99 | -0.2202 | II |  | MnCrTaAl | 6.053 | 2.99 | -0.2124 | II |
MnCrZrGe | 6.157 | 2.99 | -0.1473 | II |  | FeVHfGe | 6.158 | 3.00 | -0.2094 | II |
MnCrZrSi | 6.076 | 3.00 | -0.2569 | II |  | FeVHfSi | 6.079 | 3.00 | -0.3187 | II |
MnCrZrSn | 6.393 | 3.00 | -0.0593 | II |  | FeVHfSn | 6.386 | 3.00 | -0.129 | II |

###  III.2 New SGS candidates

Having confirmed the SGSs reported in the literature, HTP calculations are
done following the work flow (Fig. 1(e)). In total, we have identified 80 new
SGS candidates with negative formation energies, as summarized in TABLE 2.
Among such candidates, 70 systems are also mechanically and dynamically stable
(TABLE 2). More detailed information about the mechanical stabilities
(including elastic constants), dynamical stabilities, and, magnetic moments
(including partial magnetic moments) are summarized in the big table of
Section S1 in the Supplemental Material, and the band structures shown in
Section S3 and S4 of the Supplemental Material.

In the following, we will focus on the 70 cases which fulfill all the three
stability criteria (TABLE 2). It is observed that all four types of SGSs as
sketched in Fig. 1 (a) can be found, namely, there are 28 (32, 9, and 1)
type-I (type II, type-III, and type-IV) SGSs, respectively.

According to TABLE 2, most newly predicted SGSs are quaternary Heusler
compounds. We have only found three new ternary SGSs (e.g., Co2NbX (X = Al,
Ga, and In), but they are either dynamically or mechanically unstable. As
shown in TABLE LABEL:table:literature, there are six ternary Heusler (e.g.,
Ti2CoSi, Mn2CoAl, Mn2CuAl, Cr2ZnX (X=Si, Ge, and Sn)) proposed to be SGSs. It
is mentioned above that the Cr2ZnX (X=Si, Ge, and Sn) would be more stable in
the full Heusler structure, in contrast to previous calculations. zhang_epl
The Ti2CoSi and Mn2CuAl are both mechanically unstable, thus dedicated
experimental techniques such as molecular beam epitaxy should be used to
synthesize such compounds. In this regard, Mn2CoAl is a special case. Our
analysis on systems with the same number of valence electrons as Mn2CoAl, such
as Mn2FeSi, reveals that the CBM and VBM have very strong overlap, destroying
the SGS behavior. Generally speaking, the occurrence of the SGS phase depends
significantly on detailed hybridization of the atomic orbitals, as discussed
below. Therefore, the empirical rule on the number of electrons serves only as
a qualitative guide, and explicit DFT calculations on the electronic structure
are required to identify such phases.

####  III.2.1 Magnetization

Essentially, SGSs are half metals thus the total magnetic moments are expected
to be integers, and they should obey the Slater-Pauling rule. kog_jap_2013 ;
qiang_2015 According to TABLE 2, it is obvious that when N${}_{\text{V}}$ is
26 and 28, the resulting magnetic moments are 2.0 $\mu_{B}$ and 4.0 $\mu_{B}$
following $M_{tot}=(N_{V}-24)\mu_{B}$, where $M_{tot}$ and $N_{V}$ are the
total magnetic moment and number of valence electrons per unit cell,
respectively. For the cases with N${}_{\text{V}}$ being 21, the total magnetic
moments are 3.0 $\mu_{B}$ following $M_{tot}=(N_{V}-18)$. This is consistent
with the expected values based on the Slater-Pauling rule.

Such behaviors of the magnetization for Heusler compounds can be understood
based on the atomic models, as demonstrated in previous studies. kog_jap_2013
; qiang_2015 Generally, the magnitude of the magnetic moments is caused by the
competition between the crystal field splitting (between t2g and eg states)
and the exchange splitting (between the majority and minority spin channels).
apl-material ; heusler4 ; parking_progress In Ref. kog_jap_2013, , a picture
with bonding and antibonding t2g and eg bands is applied to interpret the
quaternary Heusler compounds with one magnetic ion, due to significant
hybridization between the $d$-orbitals. It is realized that such a picture has
to be generalized in order to understand the magnetization of the quaternary
Heusler SGSs, especially for cases with more than one types of magnetic atoms.

The $t_{2g}$-$e_{g}$ picture is valid for compounds with one magnetic ion. For
instance, as shown by the density of states for PtVYAl (Fig. S1(a) in the
Supplementary), the $t_{2g}$ states in the majority spin channels are
occupied, resulting in a total magnetization of 3.0 $\mu_{B}$ per formula
unit. This is generally true for other cases with NV=21, such as XVScSn (X =
Co, Ir, and Rh), PtVYAl, and FeCrScSi. For the NV=26 cases, the t2g shells in
both spin channels are filled, while the eg state is only occupied in the
majority spin channel, leading to a total magnetization of 2.0 $\mu_{B}$ per
formula unit, as demonstrated by IrFeTiSb (Fig. S1(b) in the Supplementary).

Figure 2: (Color online) (a) Orbital splitting in octahedral crystal field
environment and the corresponding anisotropic bonding environment. Although
the energy level of original eg is higher than that of original $t_{2g}$. The
resulting $a_{1g}^{*}$ and $b_{1g}^{*}$ may be not higher than the $b_{1g}$ in
principle due to the energy level is shifted. (b) Local density of states
(LDOS) for two magnetic ions in FeVNbAl. The black and red curves are the
original $t_{2g}$ and $e_{g}$ shells. Up: LDOS for Fe with one magnetic
moment. Down: LDOS for V with 2 magnetic moments.

The crystal field splittings changes greatly for systems with two or more
magnetic atoms. It is well known that for full Heusler with chemical formula
X2YZ, the site symmetry for both Y and Z are m$\bar{3}$m (Oh), while that for
X is $\bar{4}$3m (Td). The $d$-shell will split into $t_{2g}$ and $e_{g}$
subshells in both Oh and Td crystallographic symmetries, where the difference
is that the $t_{2g}$ is lower (higher) in energy for the Oh (Td) case,
respectively. Galanakis_2002 For quaternary Heuslers, the site symmetry for X,
X′, Y, and Z sites is the same, i.e., of the Td type. However, it is observed
that the $e_{g}$-$t_{2g}$ picture is not applicable in Heusler alloys with two
magnetic ions. Taking FeVNbAl as an example, for both Fe and V atoms, as
indicated by the partial density of states shown in Fig. 2(b), the crystal
fields behave like the case with D4h symmetry, where the $t_{2g}$ shell splits
into a two-fold $e_{g}$ and an one-fold $b_{1g}$ subshells, and the $e_{g}$
shell into $a_{1g}^{*}$ and $b_{1g}^{*}$ (Fig. 2(a)). In the octahedral
crystal field environment, it is already clear that the $t_{2g}$ orbitals are
lower in energy than the $e_{g}$ orbitals. For the quaternary Heusler
compounds, it is observed that the $b_{1g}$ states originated from the
$t_{2g}$ shell can be either higher or lower in energy than the
$e_{g}$-derived $a_{1g}^{*}$ subshells (Fig. 2(a)), as discussed in detail
below.

Such a pattern of crystal field splitting can be attributed to the bonding
strength of different atomic pairs. For FeVNbAl, although the nearest neighbor
V-Nb bond length (2.64 Å) is the same as that of the nearest neighbor V-Al
bonds, the integrated COHP for the V-Nb bonds is about -4.04 eV, which is much
stronger than the V-Al bonding with an integrated COHP of -1.33 eV. Moreover,
the next nearest neighbor V-Fe bond length is about 3.05 Å, but the integrated
COHP is about -1.62 eV, which is comparable to that of the V-Al bonds. Such
features can be clearly observed from the DOS (cf. Fig. S2(a) in the
Supplemantary), where the hybridization between the $d$-orbitals of V and Nb
is obviously strong. Such splittings of the original $t_{2g}$ and $e_{g}$
shells are not resulting into the separation of the {dxy, $d_{yz}$, $d_{zx}$}
(due to local tetragonal crystal fields) or {$d_{x^{2}-y^{2}}$,
$d_{3z^{2}-r^{2}}$} (due to Jahn-Teller like distortions) orbitals (cf. Fig.
S6 in the Supplemantary), like the local tetragonal distortions on the
$d$-orbitals in the octahedral environment. In contrast, the {$d_{xy}$,
$d_{yz}$, $d_{zx}$} orbitals are still degenerated in the subshells $e_{g}$
and $b_{1g}$, which is the same for the {$d_{x^{2}-y^{2}}$,
$d_{3z^{2}-r^{2}}$} orbitals in the $a_{1g}^{*}$ and $b_{1g}^{*}$ subshells
(cf. Fig. S6 in the Supplemantary).

Following such a splitting scheme, the resulting magnetic moments for
compounds with two magnetic elements can be easily understood. For the NV=21
cases such as FeVNbAl (cf. Fig. S2 (a) and Fig. S4 in the Supplementary), the
magnetic moment of 2.0 $\mu_{B}$ on the V atoms is as a result of the occupied
$e_{g}$ subshell (originated form t2g shell) in the majority spin channel
(Fig. 2(b)). Moreover, for the Fe atoms, the $e_{g}$, $b_{1g}$, and
$a_{1g}^{*}$ subshells in the majority spin channel are occupied, while only
the $e_{g}$ and $a_{1g}^{*}$ subshells in the minority spin channel are
occupied, resulting in a magnetic moments of 1.0 $\mu_{B}$. It is noted that
in this case the $b_{1g}$ subshells can be higher in energy than the
$a_{1g}^{*}$ subshells. The magnetization of other two magnetic ion compounds
with NV=21 can also be understood in the similar way (not shown).

The similar picture can also be applied to the NV=26 cases with two magnetic
ions, where the total magnetic moments of 2.0 $\mu_{B}$ can be attributed to
1.0 $\mu_{B}$ atomic moments from two atoms. Here we take CoFeTaGe (cf. Fig.
S2 (b) and Fig. S5 in the Supplementary) as an example. The bond lengths of
nearest neighbor Co-Ge and Fe-Ge almost have the same value as 2.57 Å. However
the integrated COHP of Co-Ge is -1.62 eV, which is larger than that of Fe-Ge
(-1.03 eV). Moreover, the next nearest neighbor Co-Fe bonds have a comparable
bond length to the Co-Ge and Fe-Ge bonds (about 2.96 Å), but a much weaker
bonding with integrated COHP as -0.50 eV. The resulting crystal field
splittings are very comparable to those in the cases with NV=21 case (cf. Fig.
S5 in the Supplementary). For the Co atoms, the only unoccupied state is the
$b_{1g}^{*}$ sub-shell in the minority spin channel, whereas the majority spin
channel is fully occupied, resulting in 1.0 $\mu_{B}$ magnetic moment. For the
Fe atoms, the $t_{2g}$ is not split in both spin channels and lies below Fermi
level. The $e_{g}$ state is weakly split into $a_{1g}^{*}$ and $b_{1g}^{*}$
subshells below Fermi level in the majority spin channel. On the other hand,
in the minority spin channel the $e_{g}$ state is also split into a wide
spread $a_{1g}^{*}$ subshell below Fermi level and a localized $b_{1g}^{*}$
above Fermi level. So the majority spin channel also has one more state than
the minority spin channel, resulting in one $\mu_{B}$ magnetic moment (cf.
Fig. S2(b) and Fig. S5 in the Supplementary).

In short, it is observed that the magnetization of quaternary Heusler
compounds with two magnetic ions can be understood based on the crystal
splittings of the D4h like picture. Such splittings are originated from the
anisotropic bonding between the ions, which breaks the local Oh (Td) symmetry.
In this regard, the required band filling to achieve SGSs is more flexible for
quaternary Heusler compounds than the ternary cases. This explains also why we
found more candidate SGSs in the quaternary Heusler systems, as mentioned
above. On the other hand, as to the the only new SGS (NiFeMnAl) with 28
electrons (cf. Fig. S3 in the Supplementary), the hybridization between
$d$-orbitals from the Ni, Fe, and Mn atoms is so strong that the atomic
picture is not applicable. This is also true for the other cases with three
magnetic ions.

####  III.2.2 Properties of four types of SGSs

Figure 3: (Color online) (a), (b), (c), and (d) are the band structure (left),
anomalous hall conductivity (middle), and spin resolved longitudinal
conductivity (right) of PtVYAl (type I), FeVHfSi (type II), RuCrZrAl (type
III), and NiFeMnAl (type IV), respectively. The insert in the middle panel of
(b) displays the zoom-in of AHC $\pm$0.17 eV around the Fermi energy. The blue
and red lines denote the majority and minority spin channels, respectively.
The horizontal dashed lines indicate the Fermi energy.

For all four types of SGSs, the electronic structures for one representative
case in each class are shown in Fig. 3, together with the AHCs and spin-
resolved longitudinal conductivities. For PtVYAl, which is a type-I SGS, VBM
(at the $\Gamma$ point) and CBM (at the X point) touch with each other
indirectly in the majority spin channel, while there is a gap about 0.6 eV in
the minority spin channel. Thus, the system is expected to show typical
behavior of half metals, i.e., 100% spin polarized transport properties. For
type-II SGS as exemplified by FeVHfSi, the VBM and CBM have the opposite spin
characters and touch with each other indirectly at the Fermi level (Fig.
3(b)). In this case, the spin polarization of the resulting current can be
tuned by tailoring the Fermi energy. For RuCrZrAl which represents type-III
SGSs, the valence bands near the Fermi energy are mostly of the majority spin
character, while the conduction bands constitute both majority and minority
spin character carriers. This is in contrast to the case of NiFeMnAl (a type-
IV SGS), where the conduction bands are originated from one spin channel while
the valence bands have both majority and minority spin characters.

Such specific electronic structures for four types of SGSs can be reflected in
the transport properties in terms of the AHC and longitudinal conductivity,
shown as well in Fig. 3. Due to vanishing DOS at the Fermi energy, a common
phenomenon for the four representative SGSs is that the AHC vanishes at Fermi
level. That is, the indirect band gaps for such compounds are topologically
trivial, i.e., there exist no quantum anomalous Hall effect. This is
comparable with the experimental AHC of Mn2CoAl (also a type II SGS).
c-felsher-2013-prl Moreover, for type-II SGSs such as HfVFeSi, there is a sign
change for AHC around the Fermi energy, due to the fact that the spin
character of the carriers changes when they are excited from VBM to CBM. The
resulting derivative of AHC is as high as 1597 S/(cm$\cdot$eV), corresponding
to a large anomalous Nernst effect (ANE). In this sense, such type-II SGSs are
likely promising candidates for engineering spintronic field-effect
transistors.

The right panels of Fig. 3 display the spin resolved longitudinal
conductivities at 300 K for four types of SGSs. Like the AHC, the longitudinal
conductivities of all four SGSs are quite low due to the vanishing DOS at the
Fermi energy. For type-I SGSs as exemplified by PtVYAl, the conductivity is
mostly originated from the majority spin channel, showing typical behavior of
half-metals. For type-II SGSs (FeVHfSi), due to the VBM and CBM with opposite
spin characters, the spin polarization of the longitudinal conductivity can be
conveniently tuned by controlling the chemical potential. Such compounds may
be used to fabricate spin valves which are switchable via electrostatic
gating. In case of type-III SGSs (RuCrZrAl), above the Fermi energy the
conductivity has finite values for both spin channels, while the conductivity
is nonzero only for one spin channel below the Fermi energy (Fig. 3). Such
transport property is opposite to that of the type-IV SGSs (Fig. 3(d)). It is
an interesting question how such two types of SGSs can be utilized for future
spintronic devices.

####  III.2.3 Effects of spin-orbit coupling

Figure 4: (Color online) (a) and (b) are the band structures of IrVScSn
without and with SOC. Without SOC, the blue and the red curves are majority
and minority spin channels, respectively. (c) The calculated gap as a function
of azimuthal angle $\varphi$ (the angle between the magnetization direction
and the [100] axis) in the (001) plain. The horizontal dashed lines indicate
the Fermi energy.

It is observed that SOC can induce significant changes in the electronic
structure of SGSs, since the band gaps of SGSs are on average of small
magnitude (cf. Section S4 in Supplemental Materials). For instance, for
IrVScSn, the indirect band gap is about 58.4 meV without SOC (Fig. 4(a)). When
SOC is turned on with magnetization direction along the [001] direction, the
band gap is reduced to only 0.6 meV. Such a large change in the magnitude of
the band gap can be attributed to the fact that the CBM is mainly derived from
the Ir-$d$ orbitals, where the atomic SOC strength is about 0.5 eV. Such SOC
effect on electronic structure is particularly associated with compounds
constituted of heavy elements such as Os, Ir, and Pt, due to the strong atomic
SOC strength. Similarly, it is expected that SOC has significant influence on
the electronic structure for compounds with heavy elements such as Os, Ir, and
Pt. This is indeed confirmed by our explicit DFT calculations for IrVScSn,
IrVScSi, IrVYSn, PtVScAl, and OsCrZrAl (cf. Section S4 in the Supplemental
Materials), where the band gap size can be fine tuned by about 15 meV on
average. As to OsCrYSi, the gap is even closed and the CBM and VBM are
overlapping. Here we select IrVScSn as an example for SOC effect on SGSs.

As the SGSs are magnetic, the combination of magnetic ordering with SOC lowers
the symmetry of the systems, leading to magnetization direction dependent
physical properties. Fig. 4(c) shows the magnetization direction dependence of
the band gap for IrVScSn, as the magnetization direction rotates in the (001)
plane. Obviously, the magnitude of the band gap shows a monotonous behavior of
the sinusoidal type as a function of the azimuthal angle $\varphi$ (the angle
between the magnetization direction and the [100] axis). A maximal band gap of
14.8 meV is achieved for $\varphi=\frac{\pi}{4}$. Such changes in the fine
structure of electronic structure can be manifested by the anisotropic
magnetoresistance (AMR) effect. Using the constant relaxation time ($\tau$)
approximation, we estimated the AMR ratio at 300K following the semiclassical
transport theory, given by

$\frac{\rho(0)-\rho(\frac{\pi}{4})}{\rho(0)}=\frac{(\frac{1}{\sigma(0)}-\frac{1}{\sigma(\frac{\pi}{4})})}{\frac{1}{\sigma(0)}}=\frac{(\frac{1}{\sigma(0)/\tau}-\frac{1}{\sigma(\frac{\pi}{4})/\tau})}{\frac{1}{\sigma(0)/\tau}},$
(12)

where $\sigma(\varphi)$ ($\rho(\varphi)$) is the longitudinal conductivity
(resistivity) with the azimuthal angle $\varphi$ for the magnetization
direction in the (001) plane. This results in an AMR ratio as large as 33%. On
the other hand, the magnetocrystalline anisotropy energy between such two
cases with azimuthal angle $\varphi$=0 and $\frac{\pi}{4}$ is only about
$10^{-6}$ eV per formula unit, due to the underlying cubic symmetry.
Therefore, we suspect that such materials with large AMR ratio and easily
tunable magnetization directions can be applied for future spintronic
applications.

Figure 5: (Color online)(a) and (b) are the band structures of NiCrMnAl
without and with SOC. Without SOC, the blue and the red curves are majority
and minority spin channels, respectively. (c) The red and blue curves are AHC
and ANE results $\pm$ 0.15 eV around Fermi level, respectively. The horizontal
(a and b) and vertical (c) dashed lines indicate the Fermi energy.

####  III.2.4 SGS with direct band touching

As shown in Fig. 5, we find NiCrMnAl is a special SGS, where a direct band
touching occurs at the $\Gamma$ point. Without considering SOC, the CBM from
the minority spin channel touches the VBM with the opposite spin character.
That is, it is a type-II SGS following the classification discussed above.
Unfortunately, due to the presence of a conduction band which goes slightly
below the Fermi energy at the X point, the direct touching point is hidden.
When SOC is turned on, a band gap of 24 meV is opened locally at the $\Gamma$
point. However, the resulting band gap is topologically trivial according to
the AHC shown in Fig. 5(c), since the magnitude of the AHC changes its sign
around the Fermi energy, similar to the above discussions of FeVHfSi.
Moreover, the AHC shows a singularity for an energy about 50 meV above the
Fermi energy. This indicates there is band anti-crossing in the electronic
structure. Particularly, due to the drastic variation of the AHC with respect
to the chemical potential around Fermi level, the resulting derivative of AHC
is as large as -6,000 S/(cm$\cdot$eV) at the Fermi level. That is, gigantic
anomalous Nernst effect is expected in NiCrMnAl. Such a derivative of AHC is
much larger than that of the recent experimental realized large anomalous
Nernst effect in Mn3Sn with a value of -845 S/(cm$\cdot$eV).mn3ga1 ; mn3ga2 ;
mn3ga3 In this sense, type-II SGSs may be promising materials for anomalous
Nernst applications.

##  IV Conclusion

To summarize, we have carried out a systematic high-throughput screening for
spin-gapless semiconductors (SGSs) in quaternary Heusler compounds with 21,
26, and 28 valence electrons. After validating our calculations with the
previously reported cases, we predicted 80 new stable (based on the formation
energy) compounds as promising candidates of spin-gapless semiconductors,
where 70 cases are stable based on further evaluation of the mechanical and
dynamical stabilities. The magnetization of SGSs obeys the Slater-Pauling
rule, which can be interpreted based on a new scheme of crystal field
splitting of the D4h type. Interestingly, all four types of SGSs have been
identified among our candidate systems, where both the longitudinal
conductivity and transversal anomalous Hall conductivity are calculated. We
find the type-II SGSs are particularly interesting for spintronic applications
as the spin polarization of the longitudinal conductivity is very sensitive to
the chemical potential, while the anomalous hall conductivity changes its sign
across the Fermi level, leading to possible significant anomalous Nernst
effect. This is also true for the SGS candidate NiCrMnAl with direct touching.
Additionally, it is also demonstrated that spin orbit coupling can have
significant effect on the electronic structure of SGSs with heavy elements,
where the band gap can be tuned by the magnetization direction, resulting in
large anisotropic magnetoresistance in cubic crystals. Therefore, we suspect
that SGSs are promising materials for future spintronic applications, awaiting
further experimental and theoretical explorations.

## ACKNOWLEDGMENTS

Qiang Gao thanks the financial support from the China Scholarship Council. The
authors gratefully acknowledge computational time on the Lichtenberg High
Performance Supercomputer.

## References

  * (1)
  * (2) X.L. Wang, Phys. Rev. Lett. 100, 156404 (2008).
  * (3) S. Ouardi, G. H. Fecher, C. Felser, and J. K$\ddot{u}$bler, Phys. Rev. Lett. 110, 100401 (2013).
  * (4) M. Wang, R.P. Campion, A.W. Rushforth, K.W. Edmonds, C.T. Foxon, and B.L. Gallagher, Appl. Phys. Lett. 93, 132103 (2008).
  * (5) S. Skaftouros, K. Özdoğan, E. Şaşioğlu, and I. Galanakis, Appl. Phys. Lett. 102, 022402 (2013).
  * (6) G.Z. Xu, E.K. Liu, Y. Du, G.J. Li, G.D. Liu, W.H. Wang, and G.H. Wu, Europhys. Lett. 102, 17007 (2013).
  * (7) G. Y. Gao and K. L. Yao, Appl. Phys. Lett. 103, 232409 (2013).
  * (8) W. Feng, X. Fu, C. Wan, Z. Yuan, X. Han, N. V. Quang, and S. Cho, Phys. Status Solidi RRL 9, 641 (2015).
  * (9) L. Bainsla, A.I. Mallick, M.M. Raja, A.K. Nigam, B.S.D.C.S. Varaprasad, Y. K. Takahashi, A. Alam, K. G. Suresh and K. Hono, Phys. Rev. B 91, 104408 (2015).
  * (10) M.E. Jamer, B.A. Assaf, G.E. Sterbinsky, D.Arena, L.H. Lewis, A.A. Saul, G. Radtke, and D. Heiman, Phys. Rev. B 91 , 094409 (2015).
  * (11) J. Zhou, B. Sa, Z. Sun, C. Si, and R. Ahuja, RSC Adv. 5, 73814-73919 (2015).
  * (12) Q. Gao, H.H. Xie, L. Li, G. Lei, J.B. Deng, and X.R. Hu, Superlattice. Microstruct. 85, 536 (2015).
  * (13) X.T. Wang, et al., J. Mater. Chem. C 4, 7176-7192 (2016).
  * (14) X.L. Wang, Natl. Sci. Rev. 4, 252 (2017).
  * (15) S. Curtarolo, G.L.W. Hart, M. Buongiorno Nardelli, N. Mingo, S. Sanvito, and O. Levy, Nat. Mater. 12, 191 (2013).
  * (16) I. Opahle, A. Parma, E.J. McEniry, R. Drautz, and G.K.H. Madsen,  New J. Phys. 15, 105010 (2013).
  * (17) S. Sanvito, et al., Science Advances 14, 1602241 (2017).
  * (18) J. Carrete, W. Li, N. Mingo, S. Wang, and S. Curtarolo, Phys. Rev. X 4, 011019 (2014).
  * (19) J. Balluff, K. Diekmann, G. Reiss, and M. Meinert,  Phys. Rev. Materials 1, 034404 (2017) .
  * (20) J.G. He, M. Amsler, Y. Xia, S. S. Naghavi, V. I. Hegde, S. Hao, S. Goedecker, V. Ozolins, and C. Wolverton, Phys. Rev. Lett. 117, 046602 (2016).
  * (21) J.G. He, et al., arXiv:1802.04875.
  * (22) J.H. Ma, et al.,  e-print arXiv:1712.02278v1 .
  * (23) S. V. Faleev, Y. Ferrante, J. Jeong, M. G. Samant, B. Jones, and S. S. P. Parkin,  Phys. Rev. Appl. 7, 034022 (2017).
  * (24) S.V. Faleev, Y. Ferrante, J. Jeong, M.G. Samant, B. Jones, and S.S.P. Parkin,  Phys. Rev. Mater. 1 , 024402 (2017).
  * (25) G. Kresse and J. Furthmuller, Phys. Rev. B 54, 11169 (1996).
  * (26) G. Kresse and D. Joubert,  Phys. Rev. B 59, 1758 (1999).
  * (27) I. Opahle, G.K.H. Madsen, and R. Drautz, Phys. Chem. Chem. Phys. 14, 16197 (2012).
  * (28) D. Vanderbilt, Phys. Rev. B 41, 7892(R) (1990).
  * (29) P.E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * (30) J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * (31) P. Klaer, B. Balke, V. Alijani, J. Winterlik, G.H. Fecher, C. Felser, and H. J. Elmers,  Phys. Rev. B 84, 144413 (2011).
  * (32) X. P. Wei, Y. L. Zhang, Y. D. Chu, X. W. Sun, T. Sun, P. Guo, J. B. Deng, J. Phys. Chem. Solids 82, 28 (2015).
  * (33) V. Alijani, et al.,  Phys. Rev. B 84, 224416 (2011).
  * (34) M. Born and K. Huang, Dynamics Theory of Crystal Lattices (Oxford University Press, 1954).
  * (35) K. Koepernik and H. Eschrig, Phys. Rev. B 59, 1743(1999).
  * (36) I. Opahle, K. Koepernik, and H. Eschrig, Phys. Rev. B 60 14035 (1999).
  * (37) A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * (38) A.A. Mostofi, J.R. Yates, Y.S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 178, 685 (2008).
  * (39) G.K.H. Madsen and D.J. Singh, Comput. Phys. Commun. 175, 67 (2006).
  * (40) V.L. Deringer†, A.L. Tchougreeff, and R. Dronskowski, J. Phys. Chem. A 116, 21 (2011).
  * (41) K. Özdoğan, E. Şaşioğlu, and I. Galanakis, J. Appl. Phys. 113, 193903 (2013).
  * (42) X. Yang, X. Wu, B. Wu, Y. Feng, P. Li and H. Huang, Mater. Sci. Eng. B 209, 46 (2016).
  * (43) N. Xing, H. Li, J. Dong, R. Long, and C. Zhang, Comput. Mater. Sci. 42, 600 (2008).
  * (44) J. Nehra, V. D. Sudheesh, N. Lakshmi, and K. Venugopalan, Phys. Status Solidi RRL 7, 289 (2013).
  * (45) G.Y. Gao, L. Hu, K.L.Yao, B. Luo, N. Liu, J Alloy. Compd. 551, 539-543 (2013).
  * (46) L. Bainsla, A.I. Mallick, M.M. Raja, A.A. Coelho, A.K. Nigam, D.D. Johnson, A. Alam, K.G. Suresh,  Phys. Rev. B 92, 045201 (2015).
  * (47) X. Dai, G. Liu, G. H. Fecher, C. Felser, Y. Li, and H. Liu, J. Appl. Phys. 105, 07E901 (2009).
  * (48) H.Z. Luo, Y.P. Xin, B.H. Liu, F.B. Meng, H.Y. Liu, E.K. Liu, G.H. Wu, J. Alloy. Compd. 665, 180–185 (2016).
  * (49) Y. J. Zhang, Z. H. Liu, E. K. Liu, G. D. Liu, X. Q. Ma and G. H. Wu, Europhys. Lett. 111, 37009 (2015).
  * (50) K. Kim, L. Ward, J.G. He, A. Krishna, A. Agrawal, P. Voorhees and C. Wolverton. Accelerated Discovery of Quaternary Heusler with High-Throughput Density Functional Theory and Machine Learning. Paper presennted at: APS March Meeting 2018. March 5–9, 2018; Los Angeles, California.
  * (51) C. Felser, L. Wollmann, S. Chadov, G.H. Fecher, and S.S.P. Parkin, APL Mater. 3, 041518 (2015).
  * (52) L. Wollmann, A.K. Nayak, S.S. Parkin, C. Felser, Ann. Rev. Mater. Res. 47, 247 (2017).
  * (53) T. Graf, C. Felser, and S.S.P. Parkin, Prog. Solid State Chem. 39, 1 (2011).
  * (54) I. Galanakis, P. H. Dederichs, and N. Papanikolaou, Phys. Rev. B 66, 174429 (2002)
  * (55) G.Y. Guo and T.C. Wang, Phys. Rev. B 96, 224415 (2017).
  * (56) X.K. Li, L.C. Xu, L.C. Ding, J.H. Wang, M.S. Shen, X.F. Lu, Z.W. Zhu, and K. Behnia,  Phys.Rev.Lett. 119, 056601 (2017).
  * (57) M. Ikhlas, T. Tomita, T. Koretsune, M.T. Suzuki, D. Nishio- Hamane, R. Arita, Y. Otani, and S. Nakatsuji, Nat. Phys. 13, 1085 (2017).
  * (58)


# Detecting Stochastic Governing Laws with Observation on Stationary
Distributions 111This work is supported by the National Natural Science
Foundation of China (NSFC) (Grant No.11901536). Xiaoli Chen is supported by
the Ministry of Education, Singapore, under its Research Centre of Excellence
award to the Institute for Functional Intelligent Materials (I-FIM, project
No. EDUNC-33-18-279-V12)

Xiaoli<EMAIL_ADDRESS>Hui Wangc,333Corresponding author:
<EMAIL_ADDRESS>Jinqiao<EMAIL_ADDRESS>
aDepartment of Mathematics
National University of Singapore, 119077, Singapore
bInstitute for Functional Intelligent Materials
National University of Singapore, 117544, Singapore
c School of Mathematics and Statistics
Zhengzhou University, Zhengzhou, China
d Department of Applied Mathematics & Department of Physics,
Illinois Institute of Technology, Chicago, IL 60616, USA

###### Abstract

Mathematical models for complex systems are often accompanied with
uncertainties. The goal of this paper is to extract a stochastic differential
equation governing model with observation on stationary probability
distributions. We develop a neural network method to learn the drift and
diffusion terms of the stochastic differential equation. We introduce a new
loss function containing the Hellinger distance between the observation data
and the learned stationary probability density function. We discover that the
learnt stochastic differential equation provides a fair approximation of the
data-driven dynamical system after minimizing this loss function during the
training method. The effectiveness of our method is demonstrated in numerical
experiments.

Key words: Stochastic dynamical systems; Fokker-Planck equations; Machine
learning; Neural network.

##  1 Introduction

Mathematical models for scientific and engineering systems often involve
uncertainties and thus are often in the form of stochastic differential
equations (SDE). These stochastic dynamical systems are ubiquitous in biology,
physics, geosciences and other fields. Stochastic dynamical systems provide an
appropriate framework to investigate random phenomena [1, 2, 3, 4, 5]. Hence,
the determination of SDE models is crucial for quantifying and predicting
dynamical behaviors of the nonlinear system under random fluctuations. An SDE
is characterized by the drift term and diffusion terms. In this paper, we aim
to detect appropriate drift and diffusion terms with stationary probability
distribution data.

The stationary probability distribution of a stochastic dynamical system does
not change with time and it is the stationary solution of the corresponding
Fokker-Planck equation [6, 7, 8, 9, 10, 11, 12]. The stationary probability
distribution carries the information of the underlying stochastic system [13,
14]. The Hellinger distance is the distance between probability distributions
which characterizes how close two different distributions are. In this paper,
we use Hellinger distance to identify whether the constructed SDE is an
appropriate approximation of a data-driven stochastic dynamical system.

Neural networks can be represented as compositions of simple functions with
parameters, and such functional representations can be used for parameter
estimation of time-series data and kernel estimation [15]. There has been some
progress in learning stochastic differential equation models from noisy data.
A variation estimation method was used to learn the drift term with the
observation trajectory data [16, 17, 18, 19]. There was also an RNN-based
variational method [20], a sparse learning method [21], and a Kramers-Moyal
formulae [22] for learning stochastic dynamical systems. A stochastic adjoint
sensitivity method was proposed to learn stochastic differential equations
[23] or stochastic differential equations with jumps [24]. In [25], they used
small samples from just a few snapshots of unpaired data to infer the drift
and diffusion terms of stochastic differential equations. Moreover, in [26,
27], they learned Lévy noise parameters by deep neural networks. In [28], they
solved the steady-state Fokker-Planck equation with a small amount of data
through combining the deep KD-tree.

We have recently developed a data-driven approach [29, 30] to discover
stochastic differential equations with non-Gaussian Lévy noise using the
nonlocal Kramers-Moyal formulas, and further learned the stochastic
differential equations from discrete particle samples at different time
snapshots using the Fokker-Planck equation and physics-informed neural
networks [31].

However, in addition to sample path observation data, there are recent
advances in observing or measuring stationary probability distributions [7, 9,
11, 32]. To take advantage of these new types of data, we devise a neural
network method to extract stochastic dynamical system models with stationary
probability distribution or a long time trajectory as observation data. This
motivates our research reported in this paper. Specifically, we develop a
neural network method to extract stochastic governing laws based on
probability measures. Given observation data, we learn the drift and diffusion
terms which are approximated by two neural networks. Since if we learn the
drift and diffusion together, the results would not be unique. So in this work
we proposed two approaches. The first approach entails simply learning the
drift or diffusion terms. The second technique involves learning the drift and
diffusion terms simultaneously with one drift term observational data. We
compare our learned results in three-dimensional settings with Hellinger
distance substituted by Jensen-Shannon divergence and mean-square distance
which demonstrate the efficacy of our proposed approaches.

This paper is organized as follows. In Section 2, we present two methods for
learning stochastic governing laws based on physics informed neural networks
and Hellinger distance of probability distributions. In Section 3, we present
examples to learn the drift terms and the diffusion terms. Finally, we end
with some discussions in Section 4.

##  2 Methodology

###  2.1 Problem setup

Consider the following stochastic differential equation (SDE)

$dX_{t}=b(X_{t})dt+\sigma(X_{t})dB_{t},\;\;X_{0}=x_{0},$ (1)

where the n-dimensional vector function $b(\cdot)$ is the drift term, the
$n\times n$ matrix function $\sigma(\cdot)$ is the diffusion term, and $B_{t}$
is an n-dimensional Brownian motion.

The generator of the SDE (1) is [33]:

$\displaystyle Au=$ $\displaystyle\sum_{i=1}^{n}b_{i}\frac{\partial
u}{\partial
x_{i}}+\frac{1}{2}\sum_{i,j=1}^{n}(\sigma\sigma^{T})_{i,j}\frac{\partial^{2}u}{\partial
x_{i}\partial x_{j}}.$

The probability density function (PDF) is a quantity that carries information
of the stochastic system. The time evolving probability density function of
the solution process $X_{t}$ is governed by the Fokker-Planck equation, which
is written as follows:

$\displaystyle\partial_{t}p(x,t)=A^{*}p(x,t),~{}~{}~{}~{}~{}x\in\mathbb{R}^{n},t>0,$
(2) $\displaystyle p(x,0)=p_{0}(x),$

where $p_{0}(x)$ is the initial probability density function, $A^{*}$ is the
adjoint operator of the generator $A$ and has the following form:

$\displaystyle A^{*}p=$ $\displaystyle-\sum_{i=1}^{n}\frac{\partial}{\partial
x_{i}}(b_{i}p)+\frac{1}{2}\sum_{i,j=1}^{n}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}((\sigma\sigma^{T})_{i,j}p).$ (3)

Note that the Fokker-Planck equation is a deterministic linear partial
differential equation with an initial condition.

The stationary Fokker-Planck equation is:

$\displaystyle A^{*}p(x)=$ $\displaystyle 0,$ (4)

with a condition $\int_{{\mathbb{R}}^{n}}p(x)dx=1$. We assume that there
exists a unique stationary probability density function (still denote it by
$p(x)$) in this paper.

We consider the scenario when available data is time series observation data
or probability density function. Our objective is to infer the drift and
diffusion terms. Because of the stochasticity of the dynamical system, we
could not use the mean square error to get the loss function of the SDE (1).
The main issue is how to quantify the stochastic dynamics with use the
deterministic indexes. For example, we can use the maximal likelihood
estimation [34] or the most transition pathway [35] to extract or learn the
SDE model. Here we will use the stationary probability density function as the
deterministic index to learn the SDE. If the available data is long time
trajectory data of $X(t)$, we may first use kernel density estimation to learn
the probability density function. We will propose a machine learning method to
learn the drift and the diffusion terms of the SDE, with different measures
for the distance of the observed probability density function.

###  2.2 Machine Learning

As the drift $b$ and diffusion $\sigma$ characterize the uncertainty of the
SDE, we will estimate them based on observations of probability distributions
(i.e., probability measures) of the system paths $X_{t}$. Now we introduce the
Hellinger distance [36, 37] between two probability distributions. It is used
to quantify the distance between two probability distributions in the space of
probability measures. For our purpose here, the Hellinger distance $H(p,q)$
between two probability density functions $p(x)$ and $q(x)$ is defined as
follows,

$H(p,q)\triangleq\sqrt{\frac{1}{2}\int(\sqrt{p(x)}-\sqrt{q(x)}\;\;)^{2}dx},$
(5)

which satisfies the property: $0\leq H(p,q)\leq 1$.

With the observed stationary probability density $q(x)$, we determine or
estimate the drift term $b(x)$ by minimizing the Hellinger distance between
the true stationary probability density $p(x)$ for the solution process $X(t)$
and the observation probability density $q(x)$.

Note that Hellinger distance is a measure to describe the distance of two
probability density functions. Other distance also can describe the distance.
Such as, given the probability density function $p(x)$ and $q(x)$,
respectively, the Kullback-Leibler (KL) divergence is defined as

$\displaystyle H_{KL}(p||q)=\int p(x)\log(\frac{p(x)}{q(x)})dx.$ (6)

While the Kullback-Leibler divergence is asymmetry, there also exists a
symmetric measure between two probability density, which is Jensen-Shannon
divergence, introduced as follows:

$\displaystyle
H_{JS}(p||q)=\frac{1}{2}H_{KL}(p||\frac{q+p}{2})+\frac{1}{2}H_{KL}(\frac{q+p}{2}||p).$
(7)

Later we will also use the Jensen-Shannon divergence to measure the distance.

Given the noise intensity $\sigma(x)$ and observation of the stationary
probability density $q(x)$, we will learn the drift term $b$. We devise two
neural networks to approximate the drift term and stationary probability
density $p(x)$, where the input is the space domain $x$ and the output is the
$b_{NN}(x)$ and $p_{NN}(x)$.

On the one hand, the output of $p_{NN}(x)$ should satisfy the functional (21).
We define the loss function as:

$Loss_{H}=\frac{1}{2N_{H}}\sum_{i=1}^{N_{H}}(\sqrt{p_{NN}(x_{i})}-\sqrt{q(x_{i})})^{2},$
(8)

where $\\{x_{i}\\}_{i=1}^{N_{H}}$ are the points in the spatial domain to
compute the integral and $N_{H}$ is the number of the observation data.

On the other hand, the neural networks of the drift term $b_{NN}(x)$ and the
stationary probability density $p_{NN}(x)$ should satisfy the steady Fokker-
Planck equation (4).

Similar to the physics informed neural network [39, 40], we define the
residual neural network as

$f(x)=-\sum_{i=1}^{n}\frac{\partial}{\partial
x_{i}}(b_{iNN}p_{NN})+\frac{1}{2}\sum_{i,j=1}^{n}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}((\sigma\sigma^{T})_{i,j}p_{NN}).$ (9)

Then the loss function of the residual neural network is defined as:

$Loss_{f}=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}(f(x_{i}))^{2},$ (10)

where $\\{x_{i}\\}_{i=1}^{N_{f}}$ is the residual points in the spatial domain
and $N_{f}$ is the number of the residual points. Here we randomly choose the
residual points at each iteration step. The sketch of the method is shown in
Figure 1.

Figure 1: Schematic of the neural network for solving PDEs: two neural
networks for approximate the probability p and the drift term b, where the
input is x and $\varphi$ is the activation function.

The total loss function is

$Loss=Loss_{H}+Loss_{f}.$ (11)

For the unknown drift term and diffusion term, because $b(x)=0$ and
$\sigma(x)=0$ is also the minimization solution of loss function, thus we
could not learn the terms uniquely if we train the loss function (11). The
observation data [40] of drift term at some points need to know. To avoid the
zeros solution, the observation data of drift term at few points is given,
i.e. $\\{x_{i},b(x_{i})\\}_{i=1}^{N_{b}}$. The loss function of the drift term
is:

$Loss_{b}=\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}(b_{NN}(x_{i})-b(x_{i}))^{2},$ (12)

so the loss function is defined as:

$Loss_{2}=Loss_{H}+Loss_{f}+Loss_{b}.$ (13)

If we use the Jensen-Shannon divergence to measure the distance, we just
replace the $loss_{H}$ to

$Loss_{JS}=H_{JS}(p_{NN}||q).$ (14)

We also compare our method (11) with the traditional physics informed neural
network (PINN) method [39] for solving the inverse problem of the Fokker-
Planck equation. With the observation data of the stationary probability
density function $q(x)$, the loss function is written as mean square error:

$\displaystyle
Loss_{ob}=\frac{1}{N_{ob}}\sum_{i=1}^{N_{ob}}(p_{NN}(x_{i})-q(x_{i}))^{2}.$
(15)

The total loss of PINN method is defined as:

$\displaystyle Loss_{PINN}=Loss_{ob}+Loss_{f}.$ (16)

If the observed data is the trajectory of stochastic differential equation,
i.e., $X=(X_{t_{0}},X_{t_{1}},\cdots,X_{t_{N}})$, we can use the kernel
density estimation to obtain the probability density function, denoting as
$q_{KD}$. Similarly, we replace the probability density function loss
$loss_{H}$ in Eq. 8 or $loss_{JS}$ in Eq. 14 with another loss from the
estimated density, as $q_{KD}$.

Remark: For an SDE with non-Gaussian Lévy case

$\displaystyle dX_{t}$ $\displaystyle=b(X_{t})dt+\varepsilon
dL_{t}^{\alpha},~{}~{}X_{t}\in\mathbb{R}^{n},$ (17)

where $b(\cdot)$ is the vector drift term, and $L_{t}^{\alpha}$ is a symmetric
$\alpha$-stable Lévy process in $\mathbb{R}^{n}$. The generating triplet of
the Lévy process is $(0,0,\nu_{\alpha})$.
The corresponding nonlocal Fokker-Planck operator is [33]

$\displaystyle A^{*}p=$ $\displaystyle-\sum_{i=1}^{n}\frac{\partial}{\partial
x_{i}}(a_{i}p)+\varepsilon^{\alpha}\int_{\mathbb{R}^{n}\setminus\\{0\\}}[p(x+y)-p(x)]\nu_{\alpha}(dy),$
(18)

where $\nu_{\alpha}(dy)$ is the $\alpha$-stable Lévy measure and
$\nu_{\alpha}(dy)=C_{n,\alpha}||y||^{-n-\alpha}dy$,
$C_{n,\alpha}=\frac{\alpha\Gamma((n+\alpha)/2)}{2^{1-\alpha}\pi^{n/2}\Gamma(1-\alpha/2)}.$
The stationary probability density function is the solution for the nonlocal
equation $A^{*}p=0$. To use our method to learn the SDE (17) driven by Lévy
noise, we only need to change the loss function of residual neural network
(10) to (18). As for the nonlocal integral term, we discretize it with a
scheme in our earlier work [45], and while for the first order derivative
term, we evaluate with automatic differentiation [44].

##  3 Numerical Experiments

We first present an analytical example to learn a simple stochastic system,
with quite involved calculations. For more complex stochastic systems, we will
have to use our proposed machine learning method as demonstrated in the
following numerical experiments.

###  3.1 Analytical method for learning stochastic dynamical systems

Consider a scalar stochastic differential equation

$dX_{t}=b(X_{t})dt+\sigma dB_{t},$ (19)

with appropriate conditions on drift $b$ and diffusion $\sigma$ (see [38,
p.170]), such as, $b\leq 0$ and $\sigma\neq 0$ as well as some smoothness
requirements, there exists a unique stationary probability density $p(x)$ for
the SDE (19), as a solution of the steady Fokker-Planck equation,

$p(x)=\frac{C}{\sigma^{2}(x)}e^{\int_{x^{*}}^{x}\frac{2b(y)}{\sigma^{2}(y)}dy},$
(20)

where the positive normalization constant $C$ is chosen so that $p>0$ and
$\int_{{\mathbb{R}}}p(x)dx=1$, i.e.,

$C\triangleq
1/\int_{-\infty}^{\infty}\frac{e^{\int_{x^{*}}^{x}\frac{2b(y)}{\sigma^{2}(y)}dy}}{\sigma^{2}(x)}\;dx.$

Note that $x^{*}$ here may be an arbitrary reference point so that the
integral $\int_{x^{*}}^{x}\frac{2b(y)}{\sigma^{2}(y)}dy$ exists. Different
choice of $x^{*}$ only affects the normalization constant $C$. (Say, take
$x^{*}=0$ if that is possible).

Given the observed stationary probability density $q$, we like to find out the
true stationary probability density $p$. Consider Hellinger distance $H$
between probability densities $p(x)$ and $q(x)$.

$H(b,\sigma)\triangleq\sqrt{\frac{1}{2}}\int_{{\mathbb{R}}}(\sqrt{p(x)}-\sqrt{q(x)}\;\;)^{2}dx,$
(21)

where $p(x)$ is in (20). The corresponding Euler-Lagrange equation for $H^{2}$
is

$\frac{d}{dt}H^{2}_{\sigma}=H^{2}_{b}.$ (22)

Since the Euler-Lagrange equation is the necessary condition for the
functional to obtain the minimum. So we can solve the corresponding Euler-
Lagrange equation to get the minimum value of the functional.

######  Example 1.

Consider a specific scalar stochastic model

$dX=b(X)dt+dB_{t},$

with unknown drift $b(x)$, and given diffusion $\sigma=1$. Given an
“observation” of the stationary probability density
$q(x)=\frac{1}{\sqrt{2\pi}}\;e^{-\frac{1}{2}x^{2}}$ (the Gaussian
distribution). Find a function $b(x)$ so that the Hellinger distance
$H^{2}(b(x))=\frac{1}{2}\int_{{\mathbb{R}}}[\sqrt{p(x)}-\sqrt{q(x)}]^{2}dx$ is
minimized.

The true stationary probability density for the solution process $X_{t}$ is

$\displaystyle
p(x)=\frac{e^{2\int_{0}^{x}b(y)dy}}{\int_{-\infty}^{\infty}e^{2\int_{0}^{x}b(y)dy}dx}.$
(23)

The Euler-Lagrange equation is a necessary condition for functional minima

$I(b)=\frac{1}{2}\int_{{\mathbb{R}}}(p(x)+q(x)-2\sqrt{p(x)}\sqrt{q(x)}\;)dx.$
(24)

Submitting Eq. (23) and $q(x)$ into Eq. (24), we get

$\begin{split}I(b)&=\frac{1}{2\sqrt{2\pi}}\int_{{\mathbb{R}}}e^{-\frac{1}{2}x^{2}}dx+\frac{1}{2}\int_{{\mathbb{R}}}\frac{e^{2\int_{0}^{x}b(y)dy}}{\int_{{\mathbb{R}}}e^{2\int_{0}^{x}b(y)dy}dx}dx-\int_{{\mathbb{R}}}\sqrt{\frac{e^{-\frac{1}{2}x^{2}+2\int_{0}^{x}b(y)dy}}{\sqrt{2\pi}\int_{{\mathbb{R}}}e^{2\int_{0}^{x}b(y)dy}dx}}dx.\\\
\end{split}$ (25)

In order to get the minima of $I(b)$, we can obtain $I^{\prime}(b)=0$ and
$b(x)=-kx,k\geq 0$.
Submitting $b(x)$ into $p(x)$, then $p(x)=\sqrt{\frac{k}{\pi}}e^{-kx^{2}}$,
which satisfies $\int_{-\infty}^{\infty}p(x)dx=1$.
The error $Err=\|p(x)-q(x)\|_{H}$, submitting $p(x)$ and $q(x)$ into $Err$:

$\displaystyle Err=$
$\displaystyle\|p(x)-q(x)\|_{H}=\frac{1}{2}\int_{R}(\sqrt{p(x)}-\sqrt{q(x)})^{2}dx$
(26) $\displaystyle=$
$\displaystyle\frac{1}{4\sqrt{\pi}}+\sqrt{\frac{k}{2\pi}}-\sqrt{\frac{4k}{2\pi(k+\frac{1}{2})}}.$

This is a function $f(k)$ about variable $k$. $f_{min}$ attains when
$k=\frac{1}{2}$. So, $b(x)=-\frac{1}{2}x.$

In this example, we can luckily find the optimal drift term $b$ analytical.
While it is exceedingly difficult to compute the true drift term by hand in
many problems. So, in the cases below, we use our proposed machine learning
method to learn the drift and diffusion terms.

###  3.2 Machine learning for learning stochastic dynamical systems

The neural networks in our numerical experiments below have 4 hidden layers
and 20 neurons per layer, with $\tanh$ activation function. The weights are
initialized with truncated normal distributions. The biases are initialized as
zero. We use the Adam optimizer with a learning rate $10^{-4}$ to train the
loss function.

######  Example 2.

Consider a scalar stochastic model

$dX=b(X)dt+\sigma dB_{t},$

with drift function $b(x)=x-x^{3}$. Given an “observation” of the stationary
probability density
$q(x)=\frac{1}{A}e^{\frac{1}{\sigma^{2}}x^{2}-\frac{x^{4}}{2}}$, where
$A=\int_{{\mathbb{R}}}e^{\frac{1}{\sigma^{2}}x^{2}-\frac{x^{4}}{2}}dx$. Find a
drift function $b(x)$ so that the Hellinger distance
$I(b(x))=\frac{1}{2}\int_{{\mathbb{R}}}[\sqrt{p(x)}-\sqrt{q(x)}]^{2}dx$ is
minimized.

Given the noise intensity (diffusion) $\sigma=1$, we use two fully connected
neural networks to approximate the drift term and stationary probability
density respectively. We choose $N_{H}=1001$, $N_{f}=10000$ to train the loss
function (11). The results we learned are shown in Figure 2. In Figure 2 (a),
we plot the true drift term (black line) and the learned drift term (red
line). The neural network can approximate the drift very well. In Figure 2
(b), the given $q(x)$ and neural network result of $p_{NN}(x)$ can approximate
well too. We also plot the loss function evolves with the number of iterative
steps. The loss is less than $10^{-4}$.

(a)

(b)

(c)

Figure 2: Unknown drift of Example 2: (a) learned drift term; (b) learned
probability; (c) loss function.

For the unknown drift term and diffusion term, we train the loss function (13)
to get the optimal results. Only one observation data of the drift term at
$x=-2$ is given. The results are shown in Figure 3. We present learned drift
result in Figure 3 (a), and the diffusion term evolution predictions as the
iteration of the optimiser progresses in Figure 3 (b). The neural network of
probability density is shown in Figure 3 (c). We can see for one observation
data of the drift term, the drift and diffusion term can be learned well. When
$|x|>2$, the error of the learned drift term becomes larger than that for the
other $x$. The fundamental reason for this is that the probability in bigger
$x$ is almost zero, providing very little information there.

(a)

(b)

(c)

Figure 3:  Unknown drift and diffusion of Example 2: (a) the learned drift
term with 1 observation data at $x=-2$; (b) the learned diffusion term; (c)the
learned probability density.

We also use our method to learn the SDE model, with only one trajectory
observation data $X(t)$. The observation data is the long time trajectory of
$X(t)$ and is shown in Figure 4 (a). We use kernel density estimation to get
the probability density and then use our method to learn the drift and
diffusion terms. The results of the drift term are shown in Figure 4 (b) and
the probability density is shown in Figure 4 (c). The results validate that
our method also works with long time trajectory observation data, while the
error of the probability density function is larger than using the stationary
probability density function observation data. If we have more trajectory data
of the $X(t)$, we can learn the probability density function better.

Figure 4:  Unknown drift of Example 2 with one trajectory observation data:
(a) trajectory observation data $X(t)$; (b) the learned drift term; (c)the
learned probability density.

######  Example 3.

Consider a scalar stochastic model

$dX=b(X)dt+\sigma dB_{t},$

with drift function $b(x)$. Given an “observation” of the stationary
probability density $q(x)=\frac{1}{\pi}\frac{1}{1+x^{2}}$ . Find a function
$b(x)$ so that the Hellinger distance
$I(b(x))=\frac{1}{2}\int_{{\mathbb{R}}}[\sqrt{p(x)}-\sqrt{q(x)}]^{2}dx$ is
minimized.

Similar to the last example, we fix the noise intensity (diffusion)
$\sigma=1$, and use two neural networks to approximate the drift term and
stationary probability density respectively. Here we also choose $N_{H}=1001$,
$N_{f}=10000$.

$I(b(x))=\frac{1}{2}\int_{{\mathbb{R}}}(\sqrt{p(x)}-\sqrt{q(x)})^{2}dx,$ (27)

where $q(x)=\frac{1}{\pi}\frac{1}{1+x^{2}}$.

The results are shown in Figure 5. In this case, the true drift term $b(x)$ is
unknown, so we could not compare the drift term. By comparing the learned
stationary distribution and the observation data, the result shows that they
match well. What is more, the loss function is sufficiently small, and the
probability density distribution is concentrated around zero. So zero could be
the stable point of this system. Our learned drift term $b(x)$ has one stable
point zero, as in Figure 5. This is in line with our expectations.

(a)

(b)

(c)

Figure 5: Unknown drift of Example 3: (a) learned drift term; (b) learned
probability; (c) loss function of OM part $loss_{H}$ and equation part
$loss_{f}$.

Molecular and cell biology is playing an increasingly important role in life
sciences. For example, many research findings are about stochastic
fluctuations inducing phenotypic diversity in gene expression. Here we
consider a stochastic gene regulation model [41, 42].

######  Example 4.

This is a stochastic model for a transcription factor (i.e., a protein)
concentration evolution in a certain gene regulation network

$dX=b(X)dt+\sigma dB_{t},$

with drift function $b(x)=\frac{k_{f}x^{2}}{x^{2}+K_{d}}-k_{d}x+R_{bas}$. Here
the parameters are $K_{d}=10$, $k_{d}=1min^{-1}$, $k_{f}=6min^{-1}$, and
$R_{bas}=0.4min^{-1}$. We will find a drift function $b(x)$ and diffusion
$\sigma$ so that the Hellinger distance
$I(b(x),\sigma)=\frac{1}{2}\int_{{\mathbb{R}}}[\sqrt{p(x)}-\sqrt{q(x)}]^{2}dx$
is minimized.

Given the noise intensity $\sigma=1$, two neural networks are used to
approximate the drift term and stationary probability density respectively.
Here we still choose $N_{H}=1001$, $N_{f}=10000$. The result of the drift is
shown in Figure 6 (a), where the black line is the true drift term and the red
line is the neural network result. We find that he learned drift term can fit
the true result very well. And the learned probability density function is
shown in Figure 6 (b), which fits very well with the observation probability
$q(x)$. We also show the loss function of $loss_{H}$ and $loss_{f}$ in Figure
6 (c). We see that the loss function decreases fast with the iteration steps
increasing.

(a)

(b)

(c)

Figure 6:  Unknown drift of Example 4: (a) The learned drift term; (b) The
learned probability; (c) The loss function.

For the unknown drift term and diffusion term case, we train the loss function
(13) to get the optimal result. Only one observation data of drift term at
$x=5$ is given. The results are shown in Figure 7. We present learned drift
result in Figure 7 (a), and the diffusion term evolution predictions as the
iteration of the optimiser progresses in Figure 7 (b). The neural network
result of probability density is shown in Figure 7 (c). From the figures, we
see that the drift and diffusion terms can be learned well.

(a)

(b)

(c)

Figure 7:  Unknown drift of Example 4: (a) learned drift term with 1
observation data at $x=5$; (b) learned diffusion term; (c) learned probability
density function.

######  Example 5.

We now consider the following three dimensional stochastic dynamical systems
with non-polynomial drift [31]:

$d\left(\begin{array}[]{ccc}X_{t}\\\ Y_{t}\\\
Z_{t}\end{array}\right)=\left(\begin{array}[]{c}-\partial_{X_{t}}\Phi(X_{t},Y_{t},Z_{t})\\\
-\partial_{Y_{t}}\Phi(X_{t},Y_{t},Z_{t})\\\
-\partial_{Z_{t}}\Phi(X_{t},Y_{t},Z_{t})\\\
\end{array}\right)dt+\left[\begin{array}[]{ccc}\sigma_{1}&0&0\\\
0&\sigma_{2}&0\\\ 0&0&\sigma_{3}\\\
\end{array}\right]d\left(\begin{array}[]{c}B_{1,t}\\\ B_{2,t}\\\
B_{3,t}\end{array}\right),$ (28)

where the potential
$\Phi(x,y,z)=-\frac{1}{2}\log[(2\exp(\lambda_{01}(x-\lambda_{11})^{2}+\lambda_{02}(y-\lambda_{12})^{2}+\lambda_{03}(z-\lambda_{13})^{2})+\exp(\lambda_{04}(x-\lambda_{14})^{2}+\lambda_{05}(y-\lambda_{15})^{2}+\lambda_{06}(z-\lambda_{16})^{2})]$,
$\lambda_{0i}=-5,-2.5,-5,-1,-1,-1$, $\lambda_{1i}=1,1,1,-2,-1,-1$ and
$\sigma_{j}=1$, where $i=1,2,...,6$ and $j=1,2,3$. The “observation” of the
stationary probability density is $q(x,y,z)=1/Z\exp(-2\Phi(x,y,z))$, where $Z$
is the normalization parameter such that the integral of $q(x,y,z)$ on domain
$\mathrm{R}^{3}$ is equal to 1. Find the parameters in drift term and
diffusion term so that the Hellinger distance
$I=\frac{1}{2}\int_{{\mathbb{R}}^{3}}[\sqrt{p(x,y,z)}-\sqrt{q(x,y,z)}]^{2}dxdydz$
is minimized.

We use neural network to approximate the stationary probability density. And
here we choose $N_{H}=50000$, $N_{f}=5000$.

First, we learn all the parameters $\lambda_{0i}$, $\lambda_{1i}$ and
$\sigma_{j}$, for $i=1,2,...6$ and $j=1,2,3$. The results are shown in Figure
8. In Figure 8(a) and (c), the parameters of $\lambda_{0i}$ and drift term are
learned not well. While the parameters $\lambda_{2i}$ can be learned well, see
Figure 8 (b). So we will learn the parameter in the drift term and diffusion
term respectively.

(a)

(b)

(c)

Figure 8: Three dimensional results of Example 5: (a) learned $\lambda_{0i}$;
(b) learned $\lambda_{1i}$; (c) learned $\sigma_{j}$, where $i=1,2,...,6$ and
$j=1,2,3$.

On the one hand, we just learn the diffusion term given the drift term. The
results are shown in Figure 9. The parameters in the diffusion term approach
to the true parameter as the number of iterations increases.

Figure 9:  3D result the learned drift term of Example 5.

On the other hand, we learn the parameters $\lambda_{0i}$ and $\lambda_{1i}$
in the drift term given the diffusion term. The results are shown in Figure
10. We learn the results for three cases. For case I: the observation data is
clean, i.e. $q(x,y,z)$. For case II: the observation data is given with $5\%$
noise, i.e. $q(x,y,z)*(1+0.05\mathrm{N}(0,I))$, and for case III: the
observation data is given with $10\%$ noise, i.e.
$q(x,y,z)*(1+0.1\mathrm{N}(0,I))$. Here $\mathrm{N}(0,I)$ mean the standard
normal distribution. The parameters we learned are well even the observation
data has $10\%$ noise.

(Case I)

(Case II)

(Case III)

Figure 10: Three dimensional results of Example 5. Learn all parameters in the
drift terms with perturbation. Case I: clean observation data of the PDF; Case
II: $5\%$ noise observation data of the PDF; Case III: $10\%$ noise
observation data of the PDF. Left: learned $\lambda_{0i}$; right: learned
$\lambda_{1i}$, where $i=1,2,...,6$ and $j=1,2,3$.

In the following, We change the Hellinger distance to the Jensen-Shannon
divergence in the loss function. The results of learned parameters in the
drift term are shown in Figure (11). The unknown parameters can be learned
well. While compared with the Hellinger distance, this method needs more
iteration steps to train.

Figure 11: Three dimensional results of Example 5. Learn parameters in the
drift terms using Jensen-Shannon distance. Left: learned $\lambda_{0i}$;
right: learned $\lambda_{1i}$, where $i=1,2,...,6$ and $j=1,2,3$.

Here we also compare our results with the traditional physics informed neural
network (PINN) with the case of learning the parameter in the drift term. The
results are shown in Figure 12. Only several parameters in the drift term can
be learned well using PINN method. Compared with PINN method using mean square
error (16), our loss with Hellinger distance (11) would get better results. We
use the neural network to approximate the probability density function and
plot the learned probability density function when $z=-1,0.5,1$ using
different methods. The results are shown in Figure 13. Compared with the true
probability density function, the proposed method with Hillinger distance and
Jensen-Shannon divergence works better than the mean square distance.

Figure 12: Three dimensional results with PINN loss of Example 5. Learn all
the parameters in the drift term with clean observation data of the PDF. Left:
learned $\lambda_{0i}$; right: learned $\lambda_{1i}$, where $i=1,2,...,6$ and
$j=1,2,3$.

(a1)

(b2)

(c1)

(d1)

(a2)

(b2)

(c2)

(d2)

(a3)

(b3)

(c3)

(d3)

Figure 13: The probability density function of Example 5. (a1)-(a3): the true
PDF with $z=-1,0.5,1$; (b1)-(b3): the learned PDF using Hellinger distance
(11); (c1)-(c3): the learned PDF using Jensen-Shannon divergence loss (14 );
(d1)-(d3): the learned PDF using PINN loss (16).

Figure 14: The learned potential of Example 5. (a1)-(a3): the true potential;
(b1)-(b3): the learned potential using Hellinger distance (11); (c1)-(c3): the
learned potential using Jensen-Shannon divergence loss (14); (d1)-(d3): the
learned potential using PINN loss (16).

For the unknown drift term, we use our proposed method to recover the drift
term. The results are shown in Figure 14. The first row of the Figure 14
(a1,b1,c1,d1) are the the projection of $\Phi(x,y,z)$ var minimization of $z$,
i.e. $\min_{z}\Phi(x,y,z)$. The second and third rows are projected on $(x,z)$
domain and $(y,z)$ domain separately. The true projections are shown in Figure
14(a). The results with Hellinger distance are shown in Figure 14(b).
Comparing with the true potential, this example illustrates our method with
Hellinger distance works well. However, it is difficult to recover the SDE
with the Jensen-Shannon divergence (see Figure 14(c)) and PINN losses (see
Figure 14(d)). This indicates the efficiency of our proposed method using the
Hellinger distance in this example.

######  Example 6.

Finally, We now consider the following five dimensional stochastic dynamical
systems with non-polynomial drift:

$d\left(\begin{array}[]{ccccc}X_{t}\\\ Y_{t}\\\ Z_{t}\\\ V_{t}\\\
W_{t}\end{array}\right)=\left(\begin{array}[]{c}-\partial_{X_{t}}\Phi(X_{t},Y_{t},Z_{t},V_{t},W_{t})\\\
-\partial_{Y_{t}}\Phi(X_{t},Y_{t},Z_{t},V_{t},W_{t})\\\
-\partial_{Z_{t}}\Phi(X_{t},Y_{t},Z_{t},V_{t},W_{t})\\\
-\partial_{V_{t}}\Phi(X_{t},Y_{t},Z_{t},V_{t},W_{t})\\\
-\partial_{W_{t}}\Phi(X_{t},Y_{t},Z_{t},V_{t},W_{t})\\\
\end{array}\right)dt+\left[\begin{array}[]{ccccc}\sigma_{1}&0&0&0&0\\\
0&\sigma_{2}&0&0&0\\\ 0&0&\sigma_{3}&0&0\\\ 0&0&0&\sigma_{4}&0\\\
0&0&0&0&\sigma_{5}\\\ \end{array}\right]d\left(\begin{array}[]{c}B_{1,t}\\\
B_{2,t}\\\ B_{3,t}\\\ B_{4,t}\\\ B_{5,t}\end{array}\right),$

where the potential
$\Phi(x,y,z,v,w)=-\frac{1}{2}\log[exp((\lambda_{01}(x-\lambda_{11})+\lambda_{02}(y-\lambda_{12})+\lambda_{03}(z-\lambda_{13})+\lambda_{04}(v-\lambda_{14})+\lambda_{05}(w-\lambda_{15})))+exp((\lambda_{06}(x-\lambda_{16})+\lambda_{07}(y-\lambda_{17}))+\lambda_{08}(z-\lambda_{18})+\lambda_{09}(v-\lambda_{19})+\lambda_{10}(w-\lambda_{110}))]$,
$\lambda_{0i}=-1$, $\lambda_{1i}=(1,1,1,1.5,1.5,-2,-1,-1,-1,-2)$ where
$i=1,2,\cdots,10$, and $\sigma_{j}=1$, $j=1,2,3,4,5$. The “observation” of the
stationary probability density is $q(x,y,z,v,w)=1/Z\exp(-2\Phi(x,y,z,v,w))$,
where $Z$ is the normalization parameter such that the integral of
$q(x,y,z,v,w)$ on domain ${\mathbb{R}}^{5}$ is equal to 1. Find the parameters
in drift term so that the Hellinger distance
$I=\frac{1}{2}\int_{{\mathbb{R}}^{5}}[\sqrt{p(x,y,z,v,w)}-\sqrt{q(x,y,z,v,w)}]^{2}dxdydz$
is minimized.

We use a neural network to approximate the stationary probability density. And
here we choose $N_{H}=50000$, $N_{f}=5000$. We learn the parameters
$\lambda_{0i}$ and $\lambda_{1i}$ with $i=1,2,\cdots,10$ in the drift term
given the diffusion term. The results are shown in Figure 15, indicating that
our method also works for five dimensional case.

Figure 15: Five dimensional results of Example 6. Learn all the parameters in
the drift term with clean observation data of the PDF. Left: learned
$\lambda_{0i}$; right: learned $\lambda_{1i}$, where $i=1,2,...,10$.

Remark: With only one observation trajectory data, we first use kernel density
estimation to approximate the stationary probability density function, and
then learn the SDE model. This works well for one dimension, while for high
dimensional cases, we need data on multiple trajectories.

##  4 Discussion

Based on minimizing Hellinger distance between two probability distributions,
we have devised a data-driven method to extract stochastic dynamical systems
models from observation data of either long time trajectories or stationary
probability distributions. Our numerical results in one, three dimensionaland
five dimensional examples have verified that this method is feasible. We may
also take other distances, in the space of probability distributions, in our
method. Indeed, we have also tried our method using Jensen-Shannon divergence
and mean square distance.

In principle, we may extend our method to learn high dimensional stochastic
dynamical systems. But when dealing with a high dimensional case, the larger
search space for the Fokker-Planck equation makes it difficult to train the
neural network. We will try to use parallel computing [43] or active sampling
[46] to train the neural network. Moreover, our method is more stringent on
data requirements for higher dimensions cases. Therefore, in the future, we
are going to explore other method to recover stochastic differential equation
models. We will also try to use our method to learn the stochastic
differential equations with non-Gaussian Lévy noise.

This approach leads to a data-driven stochastic dynamical systems study of
random phenomena, as we can further examine dynamical behaviors of the leaned
stochastic governing models [33].

## Acknowledgements

We would like to thank Xi Chen for helpful discussions. This work is supported
by the National Natural Science Foundation of China (NSFC) (Grant No.11901536,
12141107). Xiaoli Chen is supported by the Ministry of Education, Singapore,
under its Research Centre of Excellence award to the Institute for Functional
Intelligent Materials (I-FIM, project No. EDUNC-33-18-279-V12)

## CRediT authorship contribution statement

X. Chen: Conceptualization, Software, Formal analysis, Writing-original draft.
H. Wang: Conceptualization, Formal analysis, Writing-original draft. J. Duan:
Discussion and Suggestion.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * [1] L. Arnold, Random Dynamical Systems. New York, Springer,(2003) Corrected 2nd printing.
  * [2] C. Pasquero, E. Tziperman, Statistical parameterization of heterogeneous oceanic convection, _Journal of Physical Oceanography_. (2007) 37: 214-229.
  * [3] C. Penland, P. Sura, Sensitivity of an ocean model to“details” of stochastic forcing. In _Proc. ECMWF Workshop on Represenation of Subscale Processes using Stochastic-Dynamic Models_. Reading, England, (2005)6-8 June.
  * [4] T. Gao, J. Duan, Quantifying model uncertainty in dynamical systems driven by non-gaussian Lévy stable noise with observations on mean exit time or escape probability. Communications in Nonlinear Science and Numerical Simulation. 39(2016) 1-6.
  * [5] D. Wu, M. Fu, J. Duan, Discovering mean residence time and escape probability from data of stochastic dynamical systems. Chaos. 29(9)(2019)093122.
  * [6] C. L. Hung, X. Zhang, N. Gemelke, C. Chin, Observation of scale invariance and universality in two-dimensional Bose gases. Nature. 470(7333)(2011) 236-239.
  * [7] G. Hairapetian, R. Stenzel, Observation of a stationary, current-free double layer in a plasma. Physical Review Letters. 65(2)(1990)175.
  * [8] E. J. Yarmchuk, M. J. V. Gordon, R. E. Packard, Observation of Stationary Vortex Arrays in Rotating Superfluid Helium. Physical Review Letters. 43(3)(1979)214-217.
  * [9] O. Gefen, O. Fridman, I. Ronin, N. Q. Balaban, Direct observation of single stationary-phase bacteria reveals a surprisingly long period of constant protein production activity. Proceedings of the National Academy of Sciences. 111(1)(2014)556-561.
  * [10] L. Arnold, V. Wishtutz, Stationary solutions of linear systems with additive and multiplicative noise. Stochastics-an International Journal of Probability & Stochastic Processes. 7(1-2)(1982)133-155.
  * [11] D. Liberzon, R. W. Brockett, Nonlinear feedback systems perturbed by noise: steady-state probability distributions and optimal control. Automatic Control IEEE Transactions on. 45(6)(2000)1116-1130.
  * [12] A. H. Gray, Uniqueness of Steady-State Solutions to the Fokker-Planck Equation. Journal of Mathematical Physics. 6(4)(1965)644-647.
  * [13] R. Khasminskii, Stochastic stability of differential equations,(2012) Springer.
  * [14] B. Schmalfuss, Lyapunov functions and non-trivial stationary solutions of stochastic differential equations. Dynamical Systems: An International Journal. 16(4)(2001) 303-317.
  * [15] F. Gerber, D. W. Nychka, Fast covariance parameter estimation of spatial Gaussian process models using neural networks. _Statistics and Probability_.(2021) 10(1).
  * [16] P. Batz, A. Ruttor, M. Opper, Variational estimation of the drift for stochastic differential equations from the empirical density. Journal of Statistical Mechanics- Theory and Experiment. (8)(2016) 083404.
  * [17] P. Batz, A. Ruttor, M. Opper, Approximate Bayes learning of stochastic differential equations. Physical Review E. 98(2)(2018) 022109.
  * [18] M. Opper, An estimator for the relative entropy rate of path measures for stochastic differential equations. Journal of Computational Physics. 330(2017) 127-133.
  * [19] M. Opper, Variational inference for stochastic differential equations. Annals of Physics. 531(3)(2019) 1800233.
  * [20] T. Ryder, A. Golightly, A. S. McGough, D. Prangle, Black-box variational inference for stochastic differential equations. ICML.(2018) 4423-4432.
  * [21] L. Boninsegna, F. Nüske, C. Clementi, Sparse learning of stochastic dynamical equations. Journal of chemical physics. 148(24)(2018) 241723.
  * [22] R. Tabar, Analysis and data-based reconstruction of complex nonlinear dynamical systems, Berlin/Heidelberg, (2019)Germany: Springer.
  * [23] X. Li, T. K. L. Wong, R. T. Chen, D. Duvenaud, Scalable gradients for stochastic differential equations. AISTATS. (2020) 3870-3882.
  * [24] J. Jia, A. R. Benson, Neural jump stochastic differential equations, NIPS. (2019)32.
  * [25] L. Yang, C. Daskalakis, G. E. Karniadakis, Generative Ensemble Regression: Learning Particle Dynamics from Observations of Ensembles with Physics-Informed Deep Generative Models. SIAM Journal on Scientific Computing. 44(1)(2022) B80-B99.
  * [26] H. Zhang, Y. Xu, Y. Li, J. Kurths, Statistical solution to SDEs with $\alpha$-stable Lévy noise via deep neural network. International Journal of Dynamics and Control. 8(4)(2020) 1129-1140.
  * [27] Y. Xu, H. Zhang, Y. Li, K. Zhou, Q. Liu, J. Kurths, Solving Fokker-Planck equation using deep learning. Chaos. 30(1)(2020) 013133.
  * [28] H. Zhang, Y. Xu, Q. Liu, X. Wang, Y. Li, Solving Fokker-Planck equations using deep kd-tree with a small amount of data. Nonlinear Dynamics. (2022) https://doi.org/10.1007/s11071-022-07361-2.
  * [29] Y. Li, J. Duan, A data-driven approach for discovering stochastic dynamical systems with non-Gaussian Lévy noise. Physica D. 417(2021) 132830.
  * [30] Y.Lu, Y.Li, J. Duan, Extracting stochastic governing laws by non-local Kramers¨CMoyal formulaePhil. Trans. R. Soc. A.380, (2022)20210195.
  * [31] X. Chen, L. Yang, J. Duan, G. E. Karniadakis, Solving Inverse Stochastic Problems from Discrete Particle Observations Using the Fokker-Planck Equation and Physics-Informed Neural Networks. SIAM Journal on Scientific Computing. 43(3)(2021) B811-30.
  * [32] Y. Yang, L. Nurbekyan, E. Negrini, R. Martin, M. Pasha. Optimal transport for parameter identification of chaotic dynamics via invariant measures. arXiv preprint arXiv:2104.15138 (2021).
  * [33] J. Duan, _An Introduction to Stochastic Dynamics_ , New York(2015) Cambridge University Press.
  * [34] F. Dietrich, A. Makeev, G. Kevrekidis, N. Evangelou, T. Bertalan, S. Reich, I. Kevrekidis. Learning effective stochastic differential equations from microscopic simulations: combining stochastic numerics and deep learning. arXiv preprint arXiv:2106.09004 (2021).
  * [35] X. Chen, J. Duan, J. Hu, D. Li. Data-driven method to learn the most probable transition pathway and stochastic differential equation. Physica D: Nonlinear Phenomena 443 (2023) 133559.
  * [36] S. Cha, Comprehensive Survey on Distance/Similarity Measures between Probability Density Functions. _International Journal of Mathematical Models and Methods in Applied Sciences_. 1(4)(2007) 300-307.
  * [37] R. Beran, Minimum hellinger distance estimates for parametric models. Annals of Statistics. 5(3)(1977)445-463.
  * [38] F. C. Klebaner, _Introduction to Stochastic Calculus with Applications_. Imperial College Press, London, (2005)2nd edition.
  * [39] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics. 378(2019) 686-707.
  * [40] X. Chen, J. Duan, G. E. Karniadakis, Learning and meta-learning of stochastic advection-diffusion-reaction systems from sparse measurements. European Journal of Applied Mathematics. 32(3) (2020)397-420.
  * [41] H. Wang, X. Cheng, J. Duan, J. Kurths, X. Li, Likelihood for transcriptions in a genetic regulatory system under asymmetric stable Lévy noise. Chaos. 28 (2018) 013121.
  * [42] X. Cheng, H. Wang, X. Wang, J. Duan, X. Li, Most probable transition pathways and maximal likely trajectories in a genetic regulatory system. Physica A. 531 (2019) 121779.
  * [43] K. Shukla, D. Ameya, G. Karniadakis. Parallel physics-informed neural networks via domain decomposition. Journal of Computational Physics 447 (2021): 110683.
  * [44] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind. Automatic differentiation in machine learning: a survey. Journal of Marchine Learning Research, 18 (2018): 1-43.
  * [45] X. Chen, F. Wu, J. Duan,J. Kurths, X. Li. Most probable dynamics of a genetic regulatory network under stable Lévy noise. Applied Mathematics and Computation, 348 (2019): 425-436.
  * [46] X. Yang, Y. Liu, C. Mi, X. Wang. Active learning Kriging model combining with kernel-density-estimation-based importance sampling method for the estimation of low failure probability. Journal of Mechanical Design, 140 (2018): 051402.