# On the relation between connectivity, independence and generalized
caterpillars

M. Pedramfar<EMAIL_ADDRESS>M. Shokrian<EMAIL_ADDRESS>M. Tefagh<EMAIL_ADDRESS>Department of Mathematical Sciences, Sharif
University of Technology, P.O. Box 11155-9415,
Tehran, I.R. IRAN

###### Abstract

A spanning generalized caterpillar is a spanning tree in which all vertices of
degree more than two are on a path. In this note, we find a relation between
the existence of spanning generalized caterpillar and the independence and
connectivity number in a graph. We also point out to an error in a “theorem”
in the paper “Spanning spiders and light-splitting switches”, by L. Gargano et
al. in Discrete Math. (2004), and find out a relation between another
mentioned theorem and the existence of spanning generalized caterpillar.

###### keywords:

Spanning spider; Caterpillar; Spanning generalized caterpillar; Independence

Here, we consider only finite connected graphs without loops or multiple
edges. For standard graph-theoretic terminologies not explained in this note,
we refer to [1].

We denote the $\it{degree}$ of a vertex $x$ in a graph $G$ by $d_{G}(x)$, the
$\it{independence\ number}$ by $\alpha(G)$ and its $\it{connectivity\ number}$
by $\kappa(G)$. In a tree $T$, we call a vertex $v\in V(T)$ a $\it{branch\
vertex}$ when $d_{T}(v)>2$.

As defined in [2], a tree is called a $\it{generalized\ caterpillar}$ if all
vertices of degree more than two are on a path or equivalently a caterpillar
in which hairs edges incident to the spine, are replaced by paths.

A motivation to define such a tree comes from the definition of $\it{spider}$,
i.e. a star in which hairs are replaced by paths. Obviously a generalized
caterpillar is also a generalization of spider.

In earlier results, it is shown that:

######  Theorem A.

_([3])_ The bandwidth problem is NP-complete for generalized caterpillars of
hair length at most $3$.

But we are interested in the conditions under which a graph has a spanning
generalized caterpillar (or for simplicity an SGC). In order to find some
sufficient conditions, we concentrated on the special case of spiders. In this
way we looked at the paper [4], where it is discussed about the different
aspects of the existence of spanning spiders in a graph:

######  Theorem B.

_(Proposition 1 of[4])_ It is NP-complete to decide whether a given graph $G$
admits a spanning spider.

So finding necessary and sufficient conditions for the existence of SGC in a
graph does not seem to be an easy task! One of the results which was
particularly interesting for us is the following:

######  Theorem C.

_(Theorem 8 of[4])_ Let $G$ be a connected graph. Then $s(G)\leq
2\lceil\frac{\alpha(G)}{\kappa(G)}\rceil-2$, where $s(G)$ is the minimum
number of branch vertices in a spanning tree in $G$.

But we found that the following wrong theorem was used to prove Theorem C
which is:

######  Theorem D.

_([4])_ Vertices of any graph G can be covered by at most
$\lceil\frac{\alpha(G)}{\kappa(G)}\rceil$ vertex disjoint paths.

We checked the reference [5] to which the authors have referred, and we could
not find that used result in there. Indeed, we have a counterexample to
“Theorem D”:

######  Counterexample to “Theorem D”.

Let $G=K_{m,2m}$. We have $\alpha(G)=2\kappa(G)=2m$ but we can not cover the
vertices of $G$ by at most two vertex disjoint paths $P$ and $Q$.

###### Proof.

If $P=\emptyset$, it means that $G$ has a hamiltonian path which is obviously
wrong. So $P$ and $Q$ are not empty. Since $G$ is connected there is an edge
$e$ which is incident with a vertex in $P$ and also a vertex in $Q$. Adding
this edge leads to a spanning tree $T$ which has all vertices of degree at
most $2$, except for at most two vertices of degree $3$. $G$ is bipartite, so
the edges of $T$ have exactly one vertex in the part with $m$ vertices. But we
have:

$\displaystyle|E(T)|\leq 2+2+\cdots+2+3+3=2(m-2)+6=2m+4.$

On the other hand, since $T$ is a spanning tree we know that
$|E(T)|=|V(G)|-1=3m-1$ which is greater than $2m+4$ for $m>5$. So Theorem D is
wrong. ∎

In fact, we have:

######  Theorem E.

_([6])_ Vertices of any graph G can be covered by at most
$\lceil\frac{\alpha(G)}{\kappa(G)}\rceil$ cycles.

But what is the relation of Theorem C to the existence of SGC? At first
glance, the relation between connectivity, independence and existence of
certain kind of spanning trees might not be obvious. But there exist many
theorems discussing about this relation [7].

The following theorem answers the above question:

######  Theorem 1.

If $s(G)\leq\kappa(G)$, then $G$ has a spanning generalized caterpillar.

###### Proof.

We use the following well-known lemma to prove it.

######  Lemma.

_([1])_ In any graph $G$, any set with $\kappa(G)$ vertices can be covered by
a cycle.

The proof of the lemma given above can be found in [1], page 170. Now since
there exists a spanning tree $T$ with at most $\kappa(G)$ vertices, we can add
a cycle $C$ to the edges of $T$ so that in the new spanning subgraph
$T^{\prime}$ all vertices of degree more than $2$ are on that cycle. Now
delete one arbitrary edge from $C$, like $e$, and from each cycle in
$T^{\prime}$ delete one arbitrary edge which is not in $C-e$ (as $C-e$ is a
path we can do this). At last, we obtain a spanning tree in which all vertices
of degree more than $2$ are on a path : $C-e$. So we have an SGC. ∎

Now if we consider graphs for which
$\alpha(G)\leq\frac{(\kappa(G))^{2}+\kappa(G)}{2}$, then by Theorem C we have:

$\displaystyle s(G)\leq
2\lceil\frac{\frac{(\kappa(G))^{2}+\kappa(G)}{2}}{\kappa(G)}\rceil-2=\left\\{\begin{array}[]{ll}\kappa(G),&\mbox{
if }2|\kappa(G)\\\ \kappa(G)-1,&\mbox{ if }2|\kappa(G)+1.\end{array}\right.$

So as a corollary of Theorem 1, we can deduce that for all graphs with
$\alpha(G)\leq\frac{(\kappa(G))^{2}+\kappa(G)}{2}$ we have an SGC.

Now we try to construct graphs which do not have SGC:

######  Theorem 2.

There exists a graph G having no spanning generalized caterpillar as in the
following construction:

  1. 1.

Consider an independent set of vertices $S=\\{v_{1},\ldots,v_{m}\\}$.

  2. 2.

Consider $m+2$ graphs $G_{i}=K_{2m+1,m}\ ,\ 1\leq i\leq m+2$. From each
$G_{i}$ select $m$ arbitrary vertices from the part with $2m+1$ vertices like
$\\{t_{i,1},\ldots,t_{i,m}\\}$ and for each $1\leq i\leq m+2$ add the edge
$t_{i,j}v_{r}$ for all $1\leq r,j\leq m$.

Figure 1: Cases $m=1$ and $m=2$

###### Proof.

We have shown the cases $m=1$ and $m=2$ in Figure 1. We show by contradiction
that G does not have an SGC. Suppose that we have an SGC like $T$. Consider
$P$, the set that contains all branch vertices. Let $A_{i}$ be the part with
$m$ vertices in $G_{i}$ and $B_{i}$ be the other part which has $2m+1$
vertices.

######  Claim 1.

for each $i$, $1\leq i\leq m+2$ we have $G_{i}\cap P\neq\emptyset$.

$\it{Proof\ of\ the\ claim\ 1.}$ If for some $i$ we have $G_{i}\cap
P=\emptyset$, by considering the definition of GC, it means that all vertices
of $G_{i}$ are covered by, let us say, $n$ vertex disjoint paths. And also,
each of these paths has at least one of its ends, in the set of vertices
$\\{t_{i,1},\ldots,t_{i,m}\\}$, because these are the only vertices which are
connected to the vertices outside of $G_{i}$. Hence we get $n\leq m$.

If we denote the paths by $\\{P_{1},\ldots,P_{n}\\}$, then obviously as the
graph $G_{i}$ is bipartite we have $|P_{u}\cap B|\leq|P_{u}\cap A|+1$, for
each $u$, $1\leq u\leq n$. So we obtain:

$\displaystyle 2m+1=|B|=\sum_{u=1}^{n}|P_{u}\cap
B|\leq\sum_{u=1}^{n}(|P_{u}\cap A|+1)=|A|+n\implies m+1\leq n,$

which is a contradiction with $n\leq m$. Therefore for each $i$, $1\leq i\leq
m+2$ we have $\ G_{i}\cap P\neq\emptyset$.
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\
\ \ \square$

Define $Q$ to be the shortest path in $T$ that contains all branch vertices.
Now we claim the following and the contradiction is immediate:

######  Claim 2.

We have more than $2m$ edges which are in $Q$ and also are incident to a
vertex in $S$.

$\it{Proof\ of\ the\ claim\ 2.}$ In the last claim we proved that for each
$i$, $1\leq i\leq m+2$ we have $G_{i}\cap Q\neq\emptyset$ . Obviously
$G_{i}\cap Q$ is a set of vertex disjoint paths. Take one of them like $L$.
Now we look at the ends of $L$, namely $x$ and $y$. There exist three cases:

  * (i)

$x$ or $y$ (and not both) is at the end of $Q$.

  * (ii)

both of $x,y$ are vertices of degree $2$ in $Q$.

  * (iii)

$x,y$ are the ends of $Q$.

In the case (i), if $y$ is at the end of $Q$ then since $x$ has degree $2$ in
$Q$, one of two edges is in $E(G_{i})$ and the other is an edge between
$G_{i}$ and $S$. On the other hand as $Q$ has two ends, this case can happen
for at most two $G_{i}$’s.

In the case (ii), both of the vertices have the same condition of $x$ in the
last case. So there are two edges between $G_{i}$ and $S$ which are in $Q$.

As for the case (iii), it can not happen, because $Q$ is connected and has
vertices in at least two different $G_{i}$’s. So there must exist an edge in
$Q$ but not in $L$ and incident to one of the vertices of $L$ so that $L$ can
be connected to other parts of $Q$ in other $G_{i}$’s.

Hence by considering the fact that (i) happens at most two times there exist
at least:

$\displaystyle 2+2+\cdots+2+1+1=2(m+2-2)+2=2m+2>2m$

edges between each $G_{i}$ and $S$ which are in $Q$.
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ \ \ \ \square$

But $|S|=m$ and each vertex in $Q$ has degree at most $2$. As
$\frac{2m+2}{|S|}>2$, there is a vertex $v\in S\cap Q$ which has degree at
least $3$ in $Q$. Thus $G$ does not have any SGC. ∎

It is easy to see that $\alpha(G)=(2m+1)(m+2)$ and $\kappa(G)=m$. Hence, we
have proved that if $\alpha(G)\geq(2\kappa(G)+1)(\kappa(G)+2)$, then we may
not have SGC. We don’t know whether Theorem C is true or wrong but if we can
construct ‘better’ graphs which have no SGC, i.e. a graph for which the ratio
$\frac{\alpha(G)}{\kappa(G)}$ is less than $\frac{\kappa(G)}{2}$, then we may
be able to find a counterexample for $s(G)\leq
2\lceil\frac{\alpha(G)}{\kappa(G)}\rceil-2$.

At last, we provide a bound for $\alpha(G)$ in terms of $\kappa(G)$ which
guarantees the existence of SGC:

######  Theorem 3.

If $\alpha(G)\leq 2\kappa(G)+1$, then $G$ has an spanning generalized
caterpillar.

###### Proof.

Define $V_{3}^{T}=\\{v|v\in T,d_{T}(v)=3\\}$ and let a $k-$tree be a tree with
$\Delta(T)\leq k$. We use the following result which appears in [8]:

######  Theorem F.

_([8])_ If $\alpha(G)\leq\kappa(G)+n+1$, where $0\leq n\leq\kappa(G)$, then we
have a spanning $3-$tree $T$ with $|V_{3}^{T}|\leq n$.

By using Theorem F for the case $n=\kappa(G)$ and Theorem 1 we have an SGC.
Since before adding the cycle $C$ mentioned in Theorem 1, the degree of each
vertex was at most $3$, we get an SGC with maximum degree at most $5$. ∎

## Acknowledgements

The authors greatly appreciate Prof. E.S. Mahmoodian for his corrections and
many valuable comments. We also would like to acknowledge Prof. S. Akbari for
introducing us this problem.

## References

## References

  * [1] D. B. West, Introduction to graph theory, Prentice-Hall, Inc, United States of American, 2001.
  * [2] A. M. S. Shrestha, S. Tayu, S. Ueno, Bandwidth of convex bipartite graphs and related graphs, in: Proceedings of the 17th Annual International Conference on Computing and Combinatorics, COCOON’11, Springer-Verlag, Berlin, Heidelberg, 2011, pp. 307–318.
URL http://dl.acm.org/citation.cfm?id=2033094.2033122

  * [3] B. Monien, The bandwidth minimization problem for caterpillars with hair length 3 is np-complete, SIAM J. Algebraic Discrete Methods 7 (4) (1986) 505–512.  doi:10.1137/0607057.
URL http://dx.doi.org/10.1137/0607057

  * [4] L. Gargano, M. Hammar, P. Hell, L. Stacho, U. Vaccaro, Spanning spiders and light-splitting switches, Discrete Math. 285 (1-3) (2004) 83–95.  doi:10.1016/j.disc.2004.04.005.
URL http://dx.doi.org/10.1016/j.disc.2004.04.005

  * [5] V. Chvátal, P. Erdős, A note on Hamiltonian circuits, Discrete Math. 2 (1972) 111–113.
  * [6] M. Kouider, Cycles in graphs with prescribed stability number and connectivity, J. Combin. Theory Ser. B 60 (2) (1994) 315–318.  doi:10.1006/jctb.1994.1023.
URL http://dx.doi.org/10.1006/jctb.1994.1023

  * [7] K. Ozeki, T. Yamashita, Spanning trees: a survey, Graphs Combin. 27 (1) (2011) 1–26.  doi:10.1007/s00373-010-0973-2.
URL http://dx.doi.org/10.1007/s00373-010-0973-2

  * [8] M. Tsugaki, A note on a spanning 3-tree, Combinatorica 29 (1) (2009) 127–129.  doi:10.1007/s00493-009-2349-x.
URL http://dx.doi.org/10.1007/s00493-009-2349-x


# Giant Spin Lifetime Anisotropy and Spin-Valley Locking in Silicene and
Germanene from First-Principles Density-Matrix Dynamics

Junqing Xu  Department of Chemistry and Biochemistry, University of
California, Santa Cruz, CA 95064, USA  Hiroyuki Takenaka  Department of
Chemistry and Biochemistry, University of California, Santa Cruz, CA 95064,
USA  Adela Habib  Department of Physics, Applied Physics and Astronomy,
Rensselaer Polytechnic Institute, 110 8th Street, Troy, New York 12180, USA
Ravishankar Sundararaman<EMAIL_ADDRESS>Department of Materials Science and
Engineering, Rensselaer Polytechnic Institute, 110 8th Street, Troy, New York
12180, USA  Yuan Ping<EMAIL_ADDRESS>Department of Chemistry and
Biochemistry, University of California, Santa Cruz, CA 95064, USA

###### Abstract

Through First-Principles real-time Density-Matrix (FPDM) dynamics simulations,
we investigate spin relaxation due to electron-phonon and electron-impurity
scatterings with spin-orbit coupling in two-dimensional Dirac materials -
silicene and germanene, at finite temperatures and under external fields. We
discussed the applicability of conventional descriptions of spin relaxation
mechanisms by Elliott-Yafet (EY) and D’yakonov-Perel’ (DP) compared to our
FPDM method, which is determined by a complex interplay of intrinsic spin-
orbit coupling, external fields, and electron-phonon coupling strength, beyond
crystal symmetry. For example, the electric field dependence of spin
relaxation time is close to DP mechanism for silicene at room temperature, but
rather similar to EY mechanism for germanene. Due to its stronger spin-orbit
coupling strength and buckled structure in sharp contrast to graphene,
germanene has a giant spin lifetime anisotropy and spin valley locking effect
under nonzero $E_{z}$ and relatively low temperature. More importantly,
germanene has extremely long spin lifetime ($\sim$100 ns at 50 K) and
ultrahigh carrier mobility, which makes it advantageous for spin-valleytronic
applications.

## Introduction

Shortly after the discovery of graphene, significant advances have been made
in the field of spintronics, exploiting spin transport instead of charge
transport, with much less dissipation and unprecedented potentials for low-
power electronics. Several properties are key parameters for optimizing spin
transport, such as long spin lifetime and diffusion length, high electronic
mobility, spin-valley locking effect. Long spin lifetime and diffusion length
ensure robust spin state during propagation in a device. The spin-valley
locking effect is valuable for the emerging field of research -
”valleytronics”1, 2, 3, which utilizes the valley-pseudospin degree of freedom
as basic unit for quantum information technology.

Graphene is a very promising spintronic material4, in particular, with the
extremely high electronic mobility5 and the longest known spin diffusion
length at room temperature6. However, due to its weak spin-orbit coupling
(SOC), spin-valley locking effect may be realized only through external
effects, e.g., through proximity effect - by interfacing with large SOC
materials such as transition metal dichalcogenides (TMDs)7, 8. Other 2D
materials9 including TMDs10 have also shown exciting properties for
spin-/valley-tronics1, 2, 3. For instance, in Ref. 10, the ultralong
spin/valley lifetime is observed (e.g. 2 $\mu$s in p-type monolayer WSe2 at 5
K) and the spin lifetime is insensitive to in-plane magnetic fields, which is
a signature of the spin-valley locking effects. However, in general, TMDs have
much lower electric mobility than graphene11, 12, which is not ideal for
electronic device applications.

As the counterparts of graphene, silicene and germanene have attracted
significant attention throughout a decade due to their resemblance and
distinction from graphene13, 14, 15, 16, 17. They possess several remarkable
properties, including gate-tunable carrier concentration, high electric
mobility18, 19 (Fig. S8), quantum spin hall effects20, 21, etc. Moreover, due
to the buckled geometry, their SOC strength is highly enhanced and they will
also have the advantages of tunable band structures by applying perpendicular
electric fields $E_{z}$ and high spin polarizations, which are promising for
spintronic applications14. Their electronic structure under finite $E_{z}$ has
similarity to those of TMDs: the signs of spins in $K$ and $K^{\prime}$
valleys being opposite which is imposed by time-reversal symmetry, band
splittings induced by SOC and highly polarized states. Therefore, silicene and
germanene can combine the advantages of both graphene and TMDs but avoid some
shortcomings such as limited spin lifetimes ($\leq$ tens of ns) in graphene
samples22, 23, the absence of spin-valley locking in graphene, and low
mobilities in TMDs11, 12.

Figure 1: Electronic quantities of silicene and germanene under different
perpendicular electric field ($E_{z}$). (a) The schematic diagrams of band
structures of silicene and germanene under different $E_{z}$ (see calculated
band structures in Fig. S3-S5), where $E_{cr}$ is the critical electric field
leading to zero band gap. (b) The schematic diagrams of internal magnetic
fields ${\bf B}^{\mathrm{in}}$ (blue arrows) at a Fermi circle near one Dirac
cone of planar silicene, (buckled) silicene, (buckled) germanene induced by
breaking inversion symmetry under finite $E_{z}$. ${\bf
B}_{kn}^{\mathrm{in}}=2\Delta_{kn}{\bf
S}_{kn}^{\mathrm{exp}}/\left(g_{e}\mu_{B}\right)$, where $k$ and $n$ are
k-point and band indices, respectively. $g_{e}\mu_{B}$ is the electron spin
gyromagnetic ratio. $\Delta_{kn}$ is energy difference between two conduction
(valence) bands when $n$ is a conduction (valence) band. ${\bf
S}^{\mathrm{exp}}$ is a vector and ${\bf
S}^{\mathrm{exp}}\equiv\left(S_{x}^{\mathrm{exp}},S_{y}^{\mathrm{exp}},S_{z}^{\mathrm{exp}}\right)$,
where $S_{i}^{\mathrm{exp}}$ is spin expectation value along direction $i$ and
is the diagonal element of spin matrix $s_{i}$ in Bloch basis. (c) Band gaps
(black lines) and the averaged band splitting energies
$\left\langle\Delta\right\rangle$ (red lines) between two conduction/valence
bands. (d) The averaged out-of-plane ($B_{z}^{\mathrm{in}}$) and in-plane
internal magnetic fields ($B_{x}^{\mathrm{in}}$). Throughout this work,
$\left\langle A\right\rangle$ means the average24 of electronic quantity $A$
and $\left\langle
A\right\rangle=\sum_{kn}f^{\prime}\left(\varepsilon_{kn}\right)A_{kn}/\sum_{kn}f^{\prime}\left(\varepsilon_{kn}\right)$.
$f^{\prime}$ is the derivative of Fermi-Dirac distribution function. In this
figure, for averaging, $T=300$ K and chemical potential $\mu$ is set in the
middle of the band gap.

Despite of their promising properties, the potential for spin-based
information technologies by silicene and germanene has yet to be demonstrated.
Understanding spin dynamics and transport in materials is of key importance
for spintronics and spin-based quantum information science, and one key metric
of useful spin dynamics is spin lifetime $\tau_{s}$. Compared with graphene,
of which $\tau_{s}$ has been extensively studied experimentally and
theoretically6, 25, 26, $\tau_{s}$ of silicene and germanene have not been
measured and the existing few theoretical studies were done based on
relatively simple models27, 28, which do not include realistic interactions
with phonons and impurities. Recently we developed a First-Principles Density-
Matrix (FPDM) approach with quantum descriptions of scattering processes
between electron-phonon (e-ph), electron-impurities (e-i) and electron-
electron, to simulate spin-orbit mediated spin dynamics in general solid-state
systems with arbitrary symmetry29, 30. By employing this new method, we will
be able to predict $\tau_{s}$ of silicene and germanene at finite temperatures
with realistic interactions with environment without introducing any
simplified model or empirical parameters, for ps to $\mu$s timescale
simulation.

## Results and discussions

### Electronic structure.

We first show electronic quantities of silicene and germanene under different
$E_{z}$, which are closely related to spin dynamics and essential for
understanding spin relaxation mechanisms.

Fig. 1(a) describes a schematic picture of band structures under $E_{z}$. DFT
results of band gaps and the average of band splittings
$\left\langle\Delta\right\rangle$ between two conduction/valence bands (broken
Kramers’ degeneracy under $E_{z}$ field) are shown in Fig. 1(c).
’$\left\langle\right\rangle$’ represents taking average of electronic quantity
$A$ by24 $\left\langle
A\right\rangle=\sum_{kn}f^{\prime}\left(\epsilon_{kn}\right)A_{kn}/\sum_{kn}f^{\prime}\left(\epsilon_{kn}\right)$,
where $f^{\prime}$ is the derivative of Fermi-Dirac distribution function, and
$k$ and $n$ are k-point and band indices respectively. At $E_{z}=0$, due to
time-reversal and inversion symmetries, every two bands form a Kramers
degenerate pair.31 A finite $E_{z}$ splits a Kramers pair into spin-up and
spin-down bands due to broken inversion symmetry and the splitting increases
with $E_{z}$ (see Fig.1(c)). As the earlier model Hamiltonian study 13
presented, the phase transition from topological insulators to band insulators
happens at the critical electric field $E_{cr}$ with the schematic band
structures shown in Fig. 1(a). Based on our DFT calculations, the band gaps
close in silicene and germanene (black solid and dashed lines in Fig. 1(c)) at
0.2 and 2.5 V/nm respectively, and the gaps open again above $E_{cr}$.

The band splittings between spin-up and down under a finite $E_{z}$ are
effectively induced by $k$\- and band-dependent “internal” magnetic fields
${\bf B}^{\mathrm{in}}$, which are SOC fields induced by broken inversion
symmetry.31 ${\bf B}^{\mathrm{in}}$ is defined as $2\Delta_{kn}{\bf
S}_{kn}^{\mathrm{exp}}/\left(g_{e}\mu_{B}\right)$, $g_{e}\mu_{B}$ is the
electron spin gyromagnetic ratio. ${\bf
S}^{\mathrm{exp}}\equiv\left(S_{x}^{\mathrm{exp}},S_{y}^{\mathrm{exp}},S_{z}^{\mathrm{exp}}\right)$,
where $S_{i}^{\mathrm{exp}}$ is spin expectation value along direction $i$ and
is the diagonal element of spin matrix $s_{i}$ in Bloch basis.

In Fig. 1(b), we depict schematic diagrams of ${\bf B}^{\mathrm{in}}$ (blue
arrows) on a Fermi circle near one Dirac cone of planar silicene, buckled
silicene, and buckled germanene induced by finite $E_{z}$. In planar silicene,
${\bf B}^{\mathrm{in}}$ is purely in-plane, known as Rashba Spin-Orbit fields
${\bf B}_{\mathrm{R}}^{\mathrm{in}}$. ${\bf B}^{\mathrm{in}}$ in buckled
silicene is analogous to an admixture of ${\bf B}_{\mathrm{R}}^{\mathrm{in}}$
and out-of-plane field $B_{z}^{\mathrm{in}}$, because buckled geometry results
in the presence of $B_{z}^{\mathrm{in}}$. Different from silicene, ${\bf
B}^{\mathrm{in}}$ in buckled germanene is nearly fully out-of-plane. This
indicates that stronger intrinsic SOC in germanene significantly increases
out-of-plane $B_{z}^{\mathrm{in}}$ and, as a result, the proportion of in-
plane component diminishes. We examine the averaged $B^{\mathrm{in}}$ along z,
$\left\langle B_{z}^{\mathrm{in}}\right\rangle$, and along x, $\left\langle
B_{x}^{\mathrm{in}}\right\rangle$, under $E_{z}$ shown in Fig.1(d). While
$\left\langle B_{x}^{\mathrm{in}}\right\rangle$ for silicene and germanene
increase slowly with similar values, $\left\langle
B_{z}^{\mathrm{in}}\right\rangle$, for germanene due to stronger intrinsic
SOC, rises much more rapidly than that for silicene as a function of $E_{z}$.
${\bf B}^{\mathrm{in}}$ and its anisotropy between z and x are very important
for spin relaxation.

### Spin relaxation under finite $E_{z}$.

Spin lifetime, and even spin relaxation mechanism, can be tuned by applying
electric fields32, 26. Understanding the effect of electric field is critical
for the control and manipulation of spin relaxation. we start our theoretical
studies of spin relaxation from its electric field dependence.

As we perform FPDM calculations of spin lifetimes in this work, it is
important to understand the connection and distinction between FPDM method and
previous theoretical models. Previously, spin relaxation mechanisms are often
analyzed based on phenomenological models, such as Elliott–Yafet (EY) and
D’yakonov-Perel’ (DP) mechanisms31. EY represents the spin relaxation pathway
due to spin-flip scattering. DP is activated when inversion symmetry is broken
which results in random spin precession between adjacent scattering events. We
denote their corresponding spin lifetime with $\tau_{s}^{\mathrm{EY}}$ and
$\tau_{s}^{\mathrm{DP}}$ respectively. They are often approximated by some
simplified relations33, 34, 31: (i)
$\left(\tau_{s,i}^{\mathrm{EY}}\right)^{-1}\approx 4\left\langle
b_{i}^{2}\right\rangle\left\langle\tau_{p}^{-1}\right\rangle$ (EY relation)
along direction $i$, where $\tau_{p}$ is carrier lifetime.
$b_{i}^{2}=1-2S_{i}^{\mathrm{exp}}$ is the degree of mixture of spin-up and
spin-down states, so called ”spin mixing”31, 35, and is calculated at
$E_{z}=0$. (ii)
$\left(\tau_{s,i}^{\mathrm{DP}}\right)^{-1}\approx\left\langle\tau_{p}^{-1}\right\rangle^{-1}\left\langle\text{$\Omega^{2}$}-\Omega_{i}^{2}\right\rangle$
(DP relation), where $\Omega_{i}=g_{e}\mu_{B}B_{i}^{\mathrm{in}}$ is Larmor
precession frequency31 with $B_{i}^{\mathrm{in}}$ defined earlier. Another
qualitative estimation of the total spin relaxation rate by taking into
account both mechanisms is using36
$\left(\tau_{s}^{\mathrm{E+D}}\right)^{-1}=\left(\tau_{s}^{\mathrm{EY}}\right)^{-1}+\left(\tau_{s}^{\mathrm{DP}}\right)^{-1}$.
In the following, we will compare FPDM calculations and these phenomenological
models with first-principles input parameters including $\tau_{p}$,
$b_{i}^{2}$ and $\Omega_{i}$.

Figure 2: Spin lifetimes $\tau_{s}$ (left $y$-axis) and its anisotropy (right
$y$-axis, red lines) of (a) intrinsic silicene and (c) intrinsic germanene as
a function of $E_{z}$ at 300 K. (b) and (d) are $\tau_{s,z}$ obtained by
different methods of intrinsic silicene and germanene, respectively. “FPDM” -
first-principles real-time density-matrix calculations. “EY” and “DP”
correspond to $\tau_{s,z}^{\mathrm{EY}}$ and $\tau_{s,z}^{\mathrm{DP}}$
evaluated by EY and DP relations, respectively. “E+D” corresponds to
$\tau_{s}^{\mathrm{E+D}}=1/\left[1/\tau_{s}^{\mathrm{EY}}+1/\tau_{s}^{\mathrm{DP}}\right]$.

We first investigate out-of-plane and in-plane spin lifetime $\tau_{s,z}$ and
$\tau_{s,x}$, respectively, and their anisotropy ($\tau_{s,z}/\tau_{s,x}$) at
$E_{z}=0$ and 300K. From Fig.2 (a) and (b), we find that (i) $\tau_{s,z}$ and
$\tau_{s,x}$ of silicene are much longer than germanene; (ii) large anisotropy
(10-100) of $\tau_{s}$ is observed for both materials, e.g., 47 and 14 for
silicene and germanene at zero $E_{z}$ respectively, much greater than 0.5 for
graphene26, 7. Both phenomena may be qualitatively understood based on EY
relation (typically dominant in inversion symmetric systems) as discussed in
the following. The comparison between FPDM calculations and EY relation with
first-principles inputs is shown in Fig.2 (c) and (d), and they give the same
order of magnitude of spin lifetime. Roughly speaking, the larger spin
lifetime of silicene is mainly from the smaller spin mixing $b^{2}$ compared
with germanene based on their intrinsic SOC strength. While the large spin
anisotropy for both systems is a result of large ratio of
$b_{z}^{2}/b_{x}^{2}$ (Fig. S6(a)).

We then discuss $E_{z}$ dependence of spin relaxation. From Fig. 2(a) and
2(b), $\tau_{s,x}$ of silicene and germanene rapidly reduces with increasing
$E_{z}$. This trend can be qualitatively understood as follows. Finite $E_{z}$
breaks inversion symmetry and splits Kramers degeneracy. This induces ${\bf
B}^{\mathrm{in}}$ with a rapidly increased $z$ component as shown in Fig.
1(d), and thus leads to fast in-plane spin relaxation (perpendicular to
$B_{z}^{\mathrm{in}}$) and reduces in-plane spin lifetime $\tau_{s,x}$.

Unlike $\tau_{s,x}$, out-of-plane spin lifetime $\tau_{s,z}$ of germanene is
insensitive to $E_{z}$ in Fig. 2(c) although $\tau_{s,z}$ of silicene
decreases fast with $E_{z}$(Fig. 2(a)). To better understand $\tau_{s,z}$
dependence on $E_{z}$, we compare FPDM $\tau_{s,z}$ with the model
calculations based on EY ($\tau_{s,z}^{\mathrm{EY}}$), DP
($\tau_{s,z}^{\mathrm{DP}}$) and EY+DP ($\tau_{s,z}^{\mathrm{E+D}}$)
mechanisms as introduced earlier, in Fig. 2(b) and 2(d). From Fig. 2(b) for
silicene, we show that $\tau_{s,z}^{\mathrm{E+D}}$ and
$\tau_{s,z}^{\mathrm{DP}}$ approximately agree with FPDM $\tau_{s,z}$ in
trends. For germanene, however, from Fig. 2(d) we find that
$\tau_{s,z}^{\mathrm{EY}}$ is in good agreement with the FPDM $\tau_{s,z}$,
but neither $\tau_{s,z}^{\mathrm{E+D}}$ nor $\tau_{s,z}^{\mathrm{DP}}$ capture
the qualitative trend. Therefore, $z$-direction spin relaxation in germanene
should be mostly driven by EY mechanism, insensitive to $E_{z}$. The
suppression of DP mechanism in germanene under finite $E_{z}$ may be due to
the huge ${\bf B}^{\mathrm{in}}$ anisotropy
$B_{z}^{\mathrm{in}}/B_{x}^{\mathrm{in}}$ (see Fig. 1(d) and Fig. S6(b)): as
$B_{z}^{\mathrm{in}}$ is so strong, any in-plane spins will be quickly relaxed
and all spins are pinned along $z$; thus the total spin can only decay through
direct spin-flip processes but not through spin precession driven by ${\bf
B}^{\mathrm{in}}$ which changes spin direction gradually.

### Temperature dependence of $\tau_{s,z}$ and spin-valley locking.

Figure 3: Temperature dependent $\tau_{s,z}$ of intrinsic (a) silicene and (b)
germanene under $E_{z}=0$ and $E_{z}\neq 0$ (2 and 5 V/nm for silicene and
germanene, respectively). (c) $\tau_{s,z}$ of intrinsic silicene under
$E_{z}\neq 0$ obtained by different methods. The labels of the curves have the
same meanings as in Fig. 2. (d) Relative intervalley spin relaxation
contribution $\eta$ of germanene under $E_{z}=0$ and $E_{z}\neq 0$. $\eta$ is
defined as
$\eta=\frac{\left(\tau_{s,z}^{\mathrm{inter}}\right)^{-1}}{\left(\tau_{s,z}^{\mathrm{inter}}\right)^{-1}+w\left(\tau_{s,z}^{\mathrm{intra}}\right)^{-1}}$,
where $\tau_{s,z}^{\mathrm{inter}}$ and $\tau_{s,z}^{\mathrm{intra}}$ are
intervalley and intravalley spin lifetimes, corresponding to scattering
processes between $K$ and $K^{\prime}$ valleys and within a single $K$ or
$K^{\prime}$ valley, respectively. $\eta$ being close to 1 or 0 corresponds to
intervalley or intravalley scattering dominant spin relaxation, respectively.
See more details of $\eta$ in Ref. 37.

It is important to understand the sensitivity of spin lifetime to temperature
and determine the optimal operating temperature38, 39, 30. Therefore, we show
temperature dependence of $\tau_{s,z}$ without $E_{z}$ and with $E_{z}=$ 2 and
5 V/nm (higher than $E_{cr}$ by $\sim$2 V/nm) in intrinsic silicene and
germanene respectively in Fig. 3(a) and (b). Without $E_{z}$, $\tau_{s,z}$ of
both silicene and germanene increase fast on cooling. This is the usual
behavior of EY spin lifetime simply due to weaker e-ph scattering with
lowering temperature. Since phonon occupation is smaller at a lower
temperature, $\tau_{p}$ is longer (Fig. S8), so $\tau_{s}$ is longer
($\tau_{s}\propto\tau_{p}$ with EY mechanism). Under finite $E_{z}$,
$\tau_{s,z}$ of germanene are similar to the values under zero $E_{z}$ as
shown in Fig. 2(b), which are expected from the discussions on $E_{z}$
dependence of germanene above. In sharp contrast, finite $E_{z}$ significantly
reduces $\tau_{s,z}$ of silicene and modifies the temperature dependence. To
interpret such complex temperature dependence, in Fig. 3(c), we compare FPDM
$\tau_{s,z}$ and the model ones $\tau_{s,z}^{\mathrm{EY}}$,
$\tau_{s,z}^{\mathrm{DP}}$ and $\tau_{s,z}^{\mathrm{E+D}}$ for silicene. All
model relations fail to reproduce the temperature dependence under finite
$E_{z}$ for silicene. The failure of the DP relation in particular below 200 K
is probably because it is inapplicable for weak scattering regime (weak e-ph
scattering at low temperatures)31. Other relations have been proposed for weak
scattering31, 34, 40, e.g., $\tau_{s}\sim|\Omega|^{-1}$ and $\tau_{s}\sim
2\tau_{p}$, but none can capture the correct temperature dependence of
$\tau_{s,z}$ of silicene under finite $E_{z}$, e.g.,
$\tau_{s}\sim|\Omega|^{-1}$ predicts $\tau_{s,z}<20$ ps at all temperatures
investigated here. Our theoretical studies highlight the importance of
simulating spin lifetime using FPDM method for reliable prediction of spin
lifetimes for large variations of external conditions.

Figure 4: Spin relaxation in germanene at 50 K under 0 and 5 V/nm. (a)
$\tau_{s,z}$ and (b) relative intervalley spin relaxation contribution $\eta$
of germanene as a function of excess carrier density $n$, which is electron
density $n_{e}$ minus hole density $n_{h}$, and controlled by chemical
potential $\mu$. Positive carrier density corresponds to electron doping. (c)
$\tau_{s,z}$ of germanene as a function of in-plane magnetic field - $B_{x}$.
(d) $\tau_{s,z}$ of germanene as a function of impurity density $n_{i}$ of
neutral Ge vacancy (with $n_{e}$ and $n_{h}$ kept the same as intrinsic
germanene). $\eta$ being close to 1 or 0 correspond to spin relaxation being
dominated by intervalley or intravalley scattering, respectively. The blue and
green dotted vertical lines correspond to $\mu$ at the minima of first and
second conduction bands respectively. The black and red dotted lines are
$\tau_{s,z}$ in the zero $n_{i}$ limit under $E_{z}=0$ and 5 V/nm,
respectively.

Since the bands near the Fermi energy are composed of the Dirac cone electrons
around $K$ and $K^{\prime}$ valleys in silicene and germanene, spin relaxation
process mostly arises from intervalley and intravalley scatterings. To
scrutinize contributions of intervalley and intravalley scatterings to spin
lifetime, we examine relative intervalley spin relaxation contribution $\eta$
(see its definition in the caption of Fig. 3). $\eta$ being close to 1 or 0
correspond to spin relaxation being dominated by intervalley or intravalley
scattering, respectively. From Fig. 3(d), it is found that for germanene under
$E_{z}$=5 V/nm, the $\eta$ becomes close to 1 at $T\leq$70 K. This indicates
that under these conditions spin relaxation in germanene is dominated by
intervalley process, which is a signature of so-called spin-valley locking. In
addition, the long $\tau_{s,z}$ reaches around 100 ns at 50K and $E_{z}$=5
V/nm and indicates stability of spin states against spin relaxation. Different
from germanene, in silicene, since the relatively small SOC splitting even at
$E_{z}$=2 V/nm as shown in Fig. 1(b), the spin relaxation is mostly through
intravalley scattering (not shown), which indicates silicene will hardly show
the spin-valley locking property.

Formally the terminology “spin-valley locking” means that spin index (spin-up
and spin-down) becomes locked with the valley index (e.g., $K$ or
$K^{\prime}$)41. This implies that (i) two valleys (e.g., $K$ and
$K^{\prime}$) exist with opposite spin polarizations. (ii) Within one valley,
spins are all highly polarized along one direction, and carriers should have
the same sign of spin. This needs large SOC splitting, i.e., the spin-up and
down bands being largely separated.

Point (ii) will cause highly suppressed intravalley spin relaxation (i.e.,
relaxation through scattering processes within one valley). We explain the
reason using germanene as an example: Under 5 V/nm, the SOC splitting for
germanene is large, 23 meV at $K$. To have spins with the same sign, most of
carriers should be located around the band edges (which require low
temperatures and low carrier densities). Then the intravalley spin-flip
transition between an occupied state at band edges and an empty state at the
second conduction/valence band will be rather weak at relatively low
temperatures, since phonon occupation becomes negligible at the corresponding
phonon energy (comparable to SOC splitting).

On the other hand, for intervalley e-ph processes, the corresponding phonon
wavevectors are away from $\Gamma$ (e.g., around $K$) and the minimum phonon
frequency is finite (e.g., 7 meV for germanene, see Fig. S7). Therefore, with
spin-valley locking, spin relaxation through e-ph scattering (mostly
intervalley) will be highly suppressed leading to long spin lifetime, because
of small phonon occupations at relatively low temperatures.

T (K) |  $n_{i}$ (cm-2) |  $\overline{\mu_{c}}$ (cm2/V/s) |  $D$ (cm2/s) |  $\tau_{s,z}$ (ns) |  $l_{||,s_{z}}$ ($\mu$m)
---|---|---|---|---|---
300 | 0 | 3.2$\times$104 | 830 | 0.1 | 2.9
300 | 1011 | 2.5$\times$104 | 620 | 0.1 | 2.5
50 | 0 | 3.8$\times$106 | 16700 | 97 | 400
50 | 1011 | 4.5$\times$105 | 2000 | 76 | 120
50 | 1012 | 5.8$\times$104 | 250 | 26 | 25
Table 1: Spin dynamic and transport properties of intrinsic germanene under
$E_{z}$=5 V/nm without and with neutral impurities (with impurity density
$n_{i}$). $\overline{\mu_{c}}$ is the average of electron and hole mobility.
The method of calculating mobility is given in Supporting Information. The
theoretical results of electron and hole mobility are given in Fig. S8. $D$ is
diffusion coefficient. $l_{||,s_{z}}$ is spin diffusion length of
$z$-direction spin. The formula computing $D$ and $l_{||,s_{z}}$ are given in
the main text.

The existence of spin-valley locking not only leads to long spin lifetime but
also allows the utilization of previously developed valleytronic technologies
to design germanene-based devices1, 3. Our calculations provide guidance on
the necessary conditions to realize spin-valley locking in germanene.

### Carrier density and magnetic field dependence of $\tau_{s,z}$ for
Germanene at low T.

As intervalley scattering was shown being dominant in germanene at 50 K with a
very long spin lifetime, we further investigate its spin relaxation as a
function of carrier density, which can be easily tuned by electrical gate
experimentally. We will also determine the range of carrier density where
spin-valley locking happens.

As shown in Fig. 4(a), we plot $\tau_{s,z}$ of germanene as a function of
excess carrier density $n$ (which is electron density $n_{e}$ minus hole
density $n_{h}$, and controlled by chemical potential $\mu$). We found it is
very sensitive to $n$ at the degenerate doping range ($n\gtrsim$4$\times$1010
cm-2, corresponding to $\mu$ above the conduction band minimum). As the
carriers contributing to spin relaxation have higher energies when $n$
increases, the strong $n$ dependence of $\tau_{s,z}$ should be mainly a result
of the strength of the e-ph scattering being enhanced at higher energies (Fig.
S10).

Moreover, we find that the relative intervalley spin relaxation contribution
$\eta>$0.9 at $n\leqslant 4\times 10^{10}$ cm-2 under $E_{z}$=5 V/nm in Fig.
4(b). This indicates spin relaxation being dominated by intervalley processes
and the presence of spin-valley locking at the corresponding condition,
consistent with our above discussions about spin-valley locking.

We then investigate effects of spin-valley locking under in-plane magnetic-
field $B_{x}$. From Fig. 4(c) at 50K for germanene, $B_{x}$ has weak effects
on spin lifetime of germanene under finite $E_{z}$. This is because the
applied external magnetic field is too weak compared with internal B field
${\bf B}^{\mathrm{in}}$ under finite $E_{z}$. This weak $B_{x}$ dependence of
spin lifetime is often used as an experimental evidence of spin-valley
locking10.

### Impurity effects and spin diffusion length.

Finally, we investigate the effects of the electron-impurity scattering. As an
initial theoretical investigation, we will consider only one common neutral
defect in germanene in this work - single Ge atom vacancy42. From Fig. 4(d),
we observe that $\tau_{s,z}$ are reduced by introducing impurities and the
reduction becomes significant when impurity density $n_{i}$ approaches 1012
cm-2. Another interesting observation is that under finite $E_{z}=5$ V/nm,
$\tau_{s,z}$ reduces much less than the one with $E_{z}=0$. Our simulations
suggest that if $n_{i}$ can be controlled below 1012 cm-2, especially under a
finite $E_{z}$, germanene can exhibit long spin lifetime over 100 ns at or
below 50 K.

At the end, we compute in-plane spin diffusion length $l_{||,s_{z}}$ for
$z$-direction (out-of-plane) spin polarization of germanene using the
relation31 $l_{||,s_{z}}=\sqrt{D\tau_{s,z}}$, where $D$ is diffusion
coefficient. $D$ can be estimated using the general form of Einstein
relation43 $D=\mu_{c}(n_{e}+n_{h})/\frac{d(n_{e}+n_{h})}{d\mu}$, where
$\overline{\mu_{c}}$ is the average of electron and hole mobility, which are
obtained from solving Boltzmann equation as detailed in the method section.
From Table 1, we can see $l_{||,s_{z}}$ of germanene at 300 K is 2-3 $\mu$m,
shorter than the longest measured value of graphene samples, $\sim$12 $\mu$m6.
At 50 K, as mobilities are higher and $\tau_{s,z}$ are longer, $l_{||,s_{z}}$
become 400 $\mu$m without impurities and 120 $\mu$m with $n_{i}=$1011 cm-2,
which are much longer than experimental values of graphene at different
temperatures ranging from 1 to 40 $\mu$m9.

## Conclusions

By employing our newly developed first-principles density-matrix dynamics
approach, we computed spin lifetime of two Dirac 2D materials - silicene and
germanene as a function of temperature, electrical doping and neutral
impurities, as well as applied electric and magnetic fields. We find silicene
and germanene have qualitative different spin relaxation mechanisms under
finite $E_{z}$ fields. We did systematic comparisons between our FPDM
$\tau_{s}$ and those estimated by phenomenological models with first-
principles input parameters. We find that germanene out-of-plane spin
relaxation can be qualitatively understood by EY relation, regardless of the
one with and without E field. On the other hand, spin relaxation in silicene
is more complicated: although at room temperature, the trends of $E_{z}$
dependence of silicene $\tau_{s,z}$ can be captured by a combination of EY and
DP relation, the temperature dependence of $\tau_{s,z}$ under finite $E_{z}$
cannot be explained by any simplified model relations.

We demonstrated giant spin lifetime anisotropy (two order of magnitude higher
than graphene) and provided the condition for spin-valley locking with long
spin lifetime in germanene. Specifically, we show that at a low T - 50 K,
$\tau_{s}$ of germanene can reach 100 ns and $l_{s}$ can exceed 100 $\mu$m
(longer than graphene samples), if impurity density is controlled low
($\leq$1011 cm-2). This is very promising because the spin-valley locking
property has only been realized in either TMDs which usually have much lower
carrier mobility, or graphene on substrates which have complexity of
interfacial engineering. The realization of spin-valley locking in single
materials with long spin lifetime and ultrahigh mobility opens up highly
promising pathways for spin-valleytronic applications.

## Methods

To predict spin relaxation time from first principles, we employ our newly
developed _ab initio_ density-matrix dynamics approach, which includes quantum
descriptions of various scattering processes and is applicable to general
solid-state systems.29, 30 The density matrix master equation due to the e-ph
and e-i scattering in interaction picture reads:

$\displaystyle\frac{d\rho_{12}\left(t\right)}{dt}=$
$\displaystyle\frac{1}{2}\sum_{345}\left\\{\begin{array}[]{c}\left[I-\rho\left(t\right)\right]_{13}\rho_{45}\left(t\right)\times\\\
\left[P_{32,45}^{\mathrm{e-ph}}\left(t\right)+P_{32,45}^{\mathrm{e-i}}\left(t\right)\right]\\\
-\left[I-\rho\left(t\right)\right]_{45}\rho_{32}\left(t\right)\times\\\
\left[P_{45,13}^{\mathrm{e-ph}}\left(t\right)+P_{45,13}^{\mathrm{e-i}}\left(t\right)\right]^{*}\end{array}\right\\}$
(5) $\displaystyle+H.C.,$ (6)

where $\rho$ is density matrix. H.C. is Hermitian conjugate. The subindex,
e.g., “1” is the combined index of k-point and band. The weights of k points
must be considered when doing sum over k points. $P^{\mathrm{e-ph}}$ and
$P^{\mathrm{e-i}}$ are the generalized scattering-rate matrices for the e-ph
and e-i scattering respectively. Note that $P^{c}$ with $c$ being a scattering
channel is related to its value $P^{S,c}$ in the Schrodinger picture as
$P_{1234}^{c}\left(t\right)=P_{1234}^{S,c}\mathrm{exp}\left[i\left(\varepsilon_{1}-\varepsilon_{2}-\varepsilon_{3}+\varepsilon_{4}\right)t\right].$
$P^{S,c}$ is time independent and is computed from corresponding e-ph or e-i
matrix elements and electron and phonon energies.

All energies and matrix elements are calculated on coarse $k$ and $q$ meshes
using the DFT software JDFTx,44 and are then interpolated to extremely fine
meshes in a basis of maximally localized Wannier functions.45, 46, 47 Starting
from at an initial state with a net spin, we evolve the density matrix
$\rho\left(t\right)$ through the master equation Eq. 6 for a long enough
simulation time, typically from ns to $\mu$s. We then obtain the evolution of
spin observable $S_{i}\left(t\right)$ ($i=x,y,z$) from $\rho\left(t\right)$
(Eq. S1). At the end, spin lifetime $\tau_{s,i}$ is obtained by fitting
$S_{i}\left(t\right)$ to an exponential decay curve with decay constant
$\tau_{s,i}$.

More details are given in Supporting Information Sec. SI and SII and Ref. 29.

Using the same first-principles electron and phonon energies and matrix
elements on fine meshes, we calculate the carrier mobility by solving the
linearized Boltzmann equation using a full-band relaxation-time
approximation12 (Supporting Information Sec. SVI).

## Author contributions

J.X. and H.T. performed the _ab initio_ calculations and analyses. R.S. and
Y.P. designed and supervised all aspects of the study. All authors contribute
to the writing of the manuscript.

## Acknowledgements

We thank Mani Chandra for helpful discussions. This work is supported by the
Air Force Office of Scientific Research under AFOSR Award No. FA9550-YR-1-XYZQ
and National Science Foundation under grant No. DMR-1956015. A. H.
acknowledges support from the American Association of University Women(AAUW)
fellowship program. This research used resources of the Center for Functional
Nanomaterials, which is a US DOE Office of Science Facility, and the
Scientific Data and Computing center, a component of the Computational Science
Initiative, at Brookhaven National Laboratory under Contract No. DE-SC0012704,
the lux supercomputer at UC Santa Cruz, funded by NSF MRI grant AST 1828315,
the National Energy Research Scientific Computing Center (NERSC) a U.S.
Department of Energy Office of Science User Facility operated under Contract
No. DE-AC02-05CH11231, the Extreme Science and Engineering Discovery
Environment (XSEDE) which is supported by National Science Foundation Grant
No. ACI-1548562 48, and resources at the Center for Computational Innovations
at Rensselaer Polytechnic Institute.

## References

  * Schaibley et al. 2016 Schaibley, J. R.; Yu, H.; Clark, G.; Rivera, P.; Ross, J. S.; Seyler, K. L.; Yao, W.; Xu, X. Valleytronics in 2D materials. _Nat. Rev. Mater._ 2016, _1_ , 1–15
  * Zhang et al. 2014 Zhang, L.; Gong, K.; Chen, J.; Liu, L.; Zhu, Y.; Xiao, D.; Guo, H. Generation and transport of valley-polarized current in transition-metal dichalcogenides. _Phys. Rev. B_ 2014, _90_ , 195428
  * Tao et al. 2020 Tao, L.; Naeemi, A.; Tsymbal, E. Y. Valley-spin logic gates. _Phys. Rev. Appl._ 2020, _13_ , 054043
  * Han et al. 2014 Han, W.; Kawakami, R. K.; Gmitra, M.; Fabian, J. Graphene spintronics. _Nat. Nanotechnol._ 2014, _9_ , 794–807
  * Das Sarma et al. 2011 Das Sarma, S.; Adam, S.; Hwang, E. H.; Rossi, E. Electronic transport in two-dimensional graphene. _Rev. Mod. Phys._ 2011, _83_ , 407–470
  * Drogeler et al. 2016 Drogeler, M.; Franzen, C.; Volmer, F.; Pohlmann, T.; Banszerus, L.; Wolter, M.; Watanabe, K.; Taniguchi, T.; Stampfer, C.; Beschoten, B. Spin lifetimes exceeding 12 ns in graphene nonlocal spin valve devices. _Nano Lett._ 2016, _16_ , 3533–3539
  * Cummings et al. 2017 Cummings, A. W.; Garcia, J. H.; Fabian, J.; Roche, S. Giant Spin Lifetime Anisotropy in Graphene Induced by Proximity Effects. _Phys. Rev. Lett._ 2017, _119_ , 206601
  * Gmitra and Fabian 2017 Gmitra, M.; Fabian, J. Proximity Effects in Bilayer Graphene on Monolayer ${\mathrm{WSe}}_{2}$: Field-Effect Spin Valley Locking, Spin-Orbit Valve, and Spin Transistor. _Phys. Rev. Lett._ 2017, _119_ , 146401
  * Avsar et al. 2020 Avsar, A.; Ochoa, H.; Guinea, F.; Özyilmaz, B.; van Wees, B. J.; Vera-Marun, I. J. Colloquium: Spintronics in Graphene and Other Two-Dimensional Materials. _Rev. Mod. Phys._ 2020, _92_ , 021003
  * Dey et al. 2017 Dey, P.; Yang, L.; Robert, C.; Wang, G.; Urbaszek, B.; Marie, X.; Crooker, S. A. Gate-Controlled Spin-Valley Locking of Resident Carriers in WSe2 Monolayers. _Phys. Rev. Lett._ 2017, _119_ , 137401
  * Cheng and Liu 2018 Cheng, L.; Liu, Y. What limits the intrinsic mobility of electrons and holes in two dimensional metal dichalcogenides? _J. Am. Chem. Soc._ 2018, _140_ , 17895–17900
  * Ciccarino et al. 2018 Ciccarino, C. J.; Christensen, T.; Sundararaman, R.; Narang, P. Dynamics and spin-valley locking effects in monolayer transition metal dichalcogenides. _Nano Lett._ 2018, _18_ , 5709–5715
  * Ezawa 2012 Ezawa, M. Valley-polarized metals and quantum anomalous Hall effect in silicene. _Phys. Rev. Lett._ 2012, _109_ , 055502
  * Tsai et al. 2013 Tsai, W.-F.; Huang, C.-Y.; Chang, T.-R.; Lin, H.; Jeng, H.-T.; Bansil, A. Gated silicene as a tunable source of nearly 100% spin-polarized electrons. _Nat. Commun._ 2013, _4_ , 1–6
  * Dávila et al. 2014 Dávila, M. E.; Xian, L.; Cahangirov, S.; Rubio, A.; Le Lay, G. Germanene: a novel two-dimensional germanium allotrope akin to graphene and silicene. _New J. Phys._ 2014, _16_ , 095002
  * Zhang et al. 2016 Zhang, L.; Bampoulis, P.; Rudenko, A.; Yao, Q. v.; Van Houselt, A.; Poelsema, B.; Katsnelson, M.; Zandvliet, H. Structural and electronic properties of germanene on MoS2. _Phys. Rev. Lett._ 2016, _116_ , 256804
  * Zhao et al. 2016 Zhao, J.; Liu, H.; Yu, Z.; Quhe, R.; Zhou, S.; Wang, Y.; Liu, C. C.; Zhong, H.; Han, N.; Lu, J.; et al., Rise of silicene: A competitive 2D material. _Prog. Mater. Sci._ 2016, _83_ , 24–151
  * Shao et al. 2013 Shao, Z.-G.; Ye, X.-S.; Yang, L.; Wang, C.-L. First-principles calculation of intrinsic carrier mobility of silicene. _J. Appl. Phys._ 2013, _114_ , 093712
  * Ye et al. 2014 Ye, X.-S.; Shao, Z.-G.; Zhao, H.; Yang, L.; Wang, C.-L. Intrinsic carrier mobility of germanene is larger than graphene’s: first-principle calculations. _RSC Adv._ 2014, _4_ , 21216–21220
  * Liu et al. 2011 Liu, C.-C.; Feng, W.; Yao, Y. Quantum Spin Hall Effect in Silicene and Two-Dimensional Germanium. _Phys. Rev. Lett._ 2011, _107_ , 076802
  * Acun et al. 2015 Acun, A.; Zhang, L.; Bampoulis, P.; Farmanbar, M.; van Houselt, A.; Rudenko, A. N.; Lingenfelder, M.; Brocks, G.; Poelsema, B.; Katsnelson, M. I.; Zandvliet, H. J. W. Germanene: the germanium analogue of graphene. _J. Phys. Condens. Matter_ 2015, _27_ , 443002
  * Han and Kawakami 2011 Han, W.; Kawakami, R. K. Spin relaxation in single-layer and bilayer graphene. _Phys. Rev. Lett._ 2011, _107_ , 047207
  * Kamalakar et al. 2015 Kamalakar, M. V.; Groenveld, C.; Dankert, A.; Dash, S. P. Long distance spin communication in chemical vapour deposited graphene. _Nat. Commun._ 2015, _6_ , 1–8
  * Restrepo and Windl 2012 Restrepo, O. D.; Windl, W. Full first-principles theory of spin relaxation in group-IV materials. _Phys. Rev. Lett._ 2012, _109_ , 166604
  * Cummings and Roche 2016 Cummings, A. W.; Roche, S. Effects of dephasing on spin lifetime in ballistic spin-orbit materials. _Phys. Rev. Lett._ 2016, _116_ , 086602
  * Habib et al. 2020 Habib, A.; Xu, J.; Ping, Y.; Sundararaman, R. Electric Field and Substrate Effects Dominate Spin-Phonon Relaxation in Graphene. _arXiv preprint arXiv:2012.11550_ 2020,
  * Bishnoi and Ghosh 2013 Bishnoi, B.; Ghosh, B. Spin transport in silicene and germanene. _RSC Adv._ 2013, _3_ , 26153–26159
  * Babaee Touski 2019 Babaee Touski, S. Spin transport in armchair silicene nanoribbon on the substrate: The role of charged impurity. _Phys. Status Solidi B_ 2019, _256_ , 1900082
  * Xu et al. 2020 Xu, J.; Habib, A.; Sundararamany, R.; Ping, Y. Ab initio Ultrafast Spin Dynamics in Solids. _arXiv preprint arXiv:2012.08711_ 2020,
  * Xu et al. 2020 Xu, J.; Habib, A.; Kumar, S.; Wu, F.; Sundararaman, R.; Ping, Y. Spin-Phonon Relaxation from a Universal Ab Initio Density-Matrix Approach. _Nat. Commun._ 2020, _11_ , 2780
  * Žutić et al. 2004 Žutić, I.; Fabian, J.; Sarma, S. D. Spintronics: Fundamentals and Applications. _Rev. Mod. Phys._ 2004, _76_ , 323
  * Guimarães et al. 2014 Guimarães, M. H.; Zomer, P. J.; Ingla-Aynés, J.; Brant, J. C.; Tombros, N.; van Wees, B. J. Controlling spin relaxation in hexagonal BN-encapsulated graphene with a transverse electric field. _Phys. Rev. Lett._ 2014, _113_ , 086602
  * Fabian and Sarma 1998 Fabian, J.; Sarma, S. D. Spin relaxation of conduction electrons in polyvalent metals: Theory and a realistic calculation. _Phys. Rev. Lett._ 1998, _81_ , 5624
  * Leyland et al. 2007 Leyland, W.; Harley, R.; Henini, M.; Shields, A.; Farrer, I.; Ritchie, D. Oscillatory Dyakonov-Perel spin dynamics in two-dimensional electron gases. _Phys. Rev. B_ 2007, _76_ , 195305
  * Kurpas et al. 2019 Kurpas, M.; Junior, P. E. F.; Gmitra, M.; Fabian, J. Spin-orbit coupling in elemental two-dimensional materials. _Phys. Rev. B_ 2019, _100_ , 125422
  * Bronold et al. 2002 Bronold, F. X.; Martin, I.; Saxena, A.; Smith, D. L. Magnetic-field dependence of electron spin relaxation in n-type semiconductors. _Phys. Rev. B_ 2002, _66_ , 233206
  * 37 $w$ is a weight factor related to how many percent of total $S_{z}$ can be relaxed out by intravalley scattering itself. $w$ being close to 0 and 1 correspond to intravalley scattering can only relax a small part and most of excess spin respectively. In Supporting Information Sec. SVIII and Fig. S11, we give more details about $w$.
  * Yang et al. 2015 Yang, L.; Sinitsyn, N. A.; Chen, W.; Yuan, J.; Zhang, J.; Lou, J.; Crooker, S. A. Long-Lived Nanosecond Spin Relaxation and Spin Coherence of Electrons in Monolayer MoS2 and WS2. _Nat. Phys._ 2015, _11_ , 830–834
  * Kikkawa and Awschalom 1998 Kikkawa, J. M.; Awschalom, D. D. Resonant Spin Amplification in n-Type GaAs. _Phys. Rev. Lett._ 1998, _80_ , 4313
  * Wu et al. 2010 Wu, M.; Jiang, J.; Weng, M. Spin dynamics in semiconductors. _Phys. Rep._ 2010, _493_ , 61–236
  * Xu et al. 2014 Xu, X.; Yao, W.; Xiao, D.; Heinz, T. F. Spin and pseudospins in layered transition metal dichalcogenides. _Nat. Phys._ 2014, _10_ , 343–350
  * Hastuti et al. 2019 Hastuti, D. P.; Nurwantoro, P., et al. Stability study of germanene vacancies: The first-principles calculations. _Mater. Today Commun._ 2019, _19_ , 459–463
  * Kubo 1966 Kubo, R. The fluctuation-dissipation theorem. _Rep. Prog. Phys._ 1966, _29_ , 255
  * Sundararaman et al. 2017 Sundararaman, R.; Letchworth-Weaver, K.; Schwarz, K. A.; Gunceler, D.; Ozhabes, Y.; Arias, T. A. JDFTx: Software for Joint Density-Functional Theory. _SoftwareX_ 2017, _6_ , 278–284
  * Marzari and Vanderbilt 1997 Marzari, N.; Vanderbilt, D. Maximally Localized Generalized Wannier Functions for Composite Energy Bands. _Phys. Rev. B_ 1997, _56_ , 12847
  * Brown et al. 2016 Brown, A. M.; Sundararaman, R.; Narang, P.; Goddard, W. A.; Atwater, H. A. Nonradiative Plasmon Decay and Hot Carrier Dynamics: Effects of Phonons, Surfaces, and Geometry. _ACS Nano_ 2016, _10_ , 957–966
  * Habib et al. 2018 Habib, A.; Florio, R.; Sundararaman, R. Hot Carrier Dynamics in Plasmonic Transition Metal Nitrides. _J. Opt._ 2018, _20_ , 064001
  * Towns et al. 2014 Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lathrop, S.; Lifka, D.; Peterson, G. D.; Roskies, R.; Scott, J. R.; Wilkins-Diehr, N. XSEDE: Accelerating Scientific Discovery. _Comput. Sci. Eng._ 2014, _16_ , 62–74


Table SII: Example of candidate terms for the rhs of the Euler equation $v_{t}=G(\cdot)$, Eq. (S91). Selected terms (see last column) have the following signature with respect to $P$\- and $T$-inversion: $(P,T)=(-,+)$. Terms | $P$ | $T$ | Select
---|---|---|---
$const$, $\rho$, $\rho^{2}$, $\rho^{3}$, $\ldots$, $\rho^{5}$, $v^{2}$ | + | + | $\times$
$\rho_{x}$, $\rho\rho_{x}$, $\ldots$ $\rho^{5}\rho_{x}$ | - | + | ✓
$\rho_{xx}$, $\rho\rho_{xx}$, $\rho^{2}\rho_{xx}$, $\ldots$, $v^{2}\rho_{xx}$ | + | + | $\times$
$v^{2}\rho$, $v^{2}\rho_{xx}$ | + | + | $\times$
$\rho_{x}v_{x}$, $v_{xx}$, $\rho v_{xx}$, $\rho^{2}v_{xx}$ | - | - | $\times$
$\rho_{x}^{2}$, $v_{x}^{2}$ | + | + | $\times$
$v$, $\rho v$, $\rho^{2}v$, $\rho^{3}v$ | - | - | $\times$
$v_{x}$, $\rho v_{x}$, $\rho^{2}v_{x}$ | + | - | $\times$
$vv_{x}$, $\rho vv_{x}$, $\rho^{2}vv_{x}$, $\ldots$, $\rho^{5}vv_{x}$ | - | + | ✓
$v^{2}\rho_{x}$, $v^{2}\rho\rho_{x}$, $\ldots$, $v^{2}\rho^{5}\rho_{x}$, $(\log{\rho})_{x}$, $v^{2}(\log{\rho})_{x}$ | - | + | ✓

###  III.6 Learning of hydrodynamic PDEs from partial observations

In an experimental setting, it is quite common that only some physical
observables can be directly measured. In this subsection, we propose two
approaches for PDE-learning of hydrodynamic equations from partial
observations, when only data for the evolution of density (but not velocity)
is available. Such a situtation is common in ultracold-atom experiments, since
it is relatively easy to measure the density of the atomic cloud via optical
absorption, but it is hard to directly measure the velocity field of the
atomic cloud in situ.

The evolution of density between different time snapshots could be considered
as a “movie” that contains information about particle velocity at each spatial
point. The velocity field can then be extracted by integrating the continuity
equation, Eq. (S57):

$v(t,x)=-\frac{1}{\rho(t,x)}\partial_{t}\left[\int_{-\infty}^{x}dy\,\rho(t,y)\right].$
(S98)

The right-hand side of Eq. (S98) can be directly evaluated from the data at
spatiotemporal points of interest. After the extraction of the velocity field
$v(t,x)$, we can proceed with the standard PDE-reconstruction procedure of the
hydrodynamic equation for the velocity, described in Section III.1. Applying
the method described above to the $\rho(t,x)$ data shown in Fig. S9, we
reconstruct the Euler equation from the library of candidate terms
$v_{t}=G(\cdot)$, Eq. (S66), using the BruteForce, CrossEntropy, and STRidge
algorithms (all three algorithms leading to the same result),

$v_{t}+0.97vv_{x}+9.42\rho\rho_{x}=0,$ (S99)

which is in good agreement with the theoretically expected equation (S58) and
nearly identical to Eq. (S68), where both density and velocity data were used
for PDE learning.

The method presented above efficiently solves the problem of reconstructing
the Euler equation for the velocity $v_{t}=G(\cdot)$ from partial observations
(only from the density data $\rho(t,x)$). However, this method has a few
drawbacks. (i) The velocity reconstruction procedure via Eq. (S98) introduces
additional numerical errors due to the finite-difference computation of the
time-derivative $\partial_{t}[\ldots]$. (ii) The situation worsens in the
regions when the density approaches zero, $\rho\to 0$, resulting in a
vanishing denominator in Eq. (S98), that amplifies numerical errors. (iii) The
velocity reconstruction trick (S98) works only if the continuity equation is
valid. In problems where the total number of particles is not conserved (e.g.
in the presence of three-body loss in cold-atom experiments), the continuity
equation is no longer exact and has to be modified with appropriate loss
terms.

Problem (ii) can be partially alleviated by the considering hydrodynamic
equation for the particle current,

$(\rho v)_{t}+(\rho v^{2})_{x}=-\partial_{x}P(\rho),$ (S100)

so that the we search for an unknown equation of the form $(\rho
v)_{t}=G(\cdot)$. We utilized this method in the main text for extracting PDEs
from experimental data (boson gas expansion on an atom chip).

In addition, we propose a modified PDE learning method to address problems
(ii) and (iii). We assume that

$\displaystyle\rho_{t}$
$\displaystyle=F(\xi_{1};\rho,v,\rho_{x},\rho_{xx},v_{x},v_{xx},\ldots),$
(S101) $\displaystyle v_{t}$
$\displaystyle=G(\xi_{2};\rho,v,\rho_{x},\rho_{xx},v_{x},v_{xx},\ldots),$
(S102)

where $\xi_{1,2}$ are the coefficients that parametrize functions $F$ and $G$.
We define the objective function as

$\mathcal{L}(\xi_{1,2};\rho^{*},\lambda_{0})=\sum_{t_{k},x_{i}}\left|\rho^{*}(t_{k},x_{i})-\textrm{PDESolve}(t_{k},x_{i},\rho_{0},v_{0},F(\xi_{1};\cdot),G(\xi_{2};\cdot),BC\right|+\lambda_{0}||\xi_{2}||_{0},$
(S103)

where $\rho_{0}$ and $v_{0}$ are the initial conditions at $t=0$, $BC$
represents a set of boundary conditions, and $\rho^{*}(t_{k},x_{i})$ are the
data points for the density evaluated at the spatiotemporal grid
$\\{t_{k},x_{i}\\}$. $\textrm{PDESolve}(\ldots)$ denotes a PDE solver that
takes initial and boundary conditions and the coefficients $\xi_{1,2}$
parametrizing the unknown function $f(\xi_{1,2};\cdot)$ and outputs a solution
$\rho_{pde}$ at the grid points $\\{t_{k},x_{i}\\}$. We assume that we know
the initial conditions for the velocity [in quench experiments with ultracold
atoms, one usually has $v(t=0,x)=0$]. The goal is to find a set of
coefficients $\xi_{1,2}$ that minimize the objective function:

$\xi_{1,2}=\operatorname*{arg\,min}_{\xi_{1,2}}(\mathcal{L}(\xi_{1,2};\rho^{*},\lambda_{0})).$
(S104)

Optimization of Eq. (S104) is more computationally costly compared to the
sparse-regression methods discussed in Sec. I. The former approach requires
(i) discrete optimization to find optimal combinations of non-zero terms, (ii)
continuous optimization to find optimal values of PDE coefficients $\xi_{1,2}$
(e.g. by gradient-descent algorithm), (iii) solution of the PDE forward in
time for each optimization step in order to evaluate the current value of the
objective function (S103). Computational cost of this algorithm can be
significantly improved by combining PDE solvers and reverse mode automatic
differentiation [89].

###  III.7 Could the evolution of the density of non-interacting fermions be
described by a single PDE?

In this subsection, we address the following question: is it possible to
rewrite the coupled system of hydrodynamics equations describing fermion
dynamics,

$\begin{cases}\rho_{t}+(\rho v)_{x}=0,\\\
v_{t}+vv_{x}+\pi^{2}\rho\rho_{x}=0,\end{cases}$ (S105)

as a single closed form PDE for the fermion density? The answer is positive,
although surprisingly we did not find such an equation in the literature. We
introduce an auxiliary variable $w$ as:

$\displaystyle\rho v=-w_{t},$ (S106) $\displaystyle\rho=w_{x},$ (S107)

so that the continuity equation is automatically satisfied. Solving Eq.
(S107), we obtain $w(t,x)=\int_{0}^{x}dx^{\prime}\,\rho(t,x^{\prime})+g(t)$,
where it is easy to show that $g(t)=const$ due to fixed boundary conditions at
$x\to\pm\infty$. The physical meaning of $w(t,x)$ is the number of fermions to
the left of coordinate $x$. Expressing velocity from Eqs. (S106, S107) as
$v=-w_{t}/w_{x}$ and substituting the result into the Euler equation (S105),
we obtain

$-\partial_{t}\left(\frac{w_{t}}{w_{x}}\right)+\frac{w_{t}}{w_{x}}\left(\frac{w_{t}}{w_{x}}\right)_{x}+\pi^{2}w_{x}w_{xx}=0.$
(S108)

The resulting PDE (S108) is second-order in time and depends only on the
fermion density $\rho(t,x)$ via $w$. Unfortunately, Eq. (S108) has no
transparent physical interpretation, to our knowledge. Although PDE-learning
methodology allows one to reconstruct a second-order-in-time equation
$w_{tt}=F(w_{t},w_{x},w_{tx},\ldots)$, such an equation would be hard to
interpret, compared to a conventional hydrodynamic system of equations for the
density and the velocity.

###  III.8 Hydrodynamics of interacting fermions: emergent Euler and Navier-
Stokes equations

While in previous subsections we considered hydrodynamics in a non-interacting
fermion gas, here we focus on the case of interacting fermions described by
the 1D spinless Fermi-Hubbard model:

$H=-J\sum_{i}(c^{\dagger}_{i}c_{i+1}+c^{\dagger}_{i+1}c_{i})+U\sum_{i}n_{i}n_{i+1}-\mu\sum_{i}c^{\dagger}_{i}c_{i},$
(S109)

where $U$ is the interaction constant, and $\mu$ is the chemical potential of
the fermion gas.

We perform a search of hydrodynamic-type equations for the fermion density and
velocity using the following form of hydrodynamic equations:

$\begin{cases}\rho_{t}+(\rho v)_{x}=0,\\\
v_{t}+vv_{x}=G(\rho_{x},\rho_{xx},v_{x},v_{xx},\ldots),\end{cases}$ (S110)

where the function $G(\cdot)$ is unknown. We use a truncated library of terms
corresponding to Table SII, where the terms are preselected based on the
$P$-parity transformation, $P=-1$. Using the BruteForce search algorithm (with
$\lambda_{0}=10^{-2}$), we found the following Euler equation,

$\displaystyle v_{t}+vv_{x}+\kappa(U)\rho\rho_{x}=0,$ (S111)

describing the dynamics of interacting fermions. We now provide the summary of
our analysis, consolidate some of the statements from the main text, and
expand on the connection between the discovered Euler/Navier-Stokes fluid
models and the Tomonaga-Luttinger theory:

  * •

In the limit of vanishing interaction strength $U\to 0$, the coefficient
$\kappa$ approaches the value $\kappa\approx\pi^{2}/m^{2}=4J^{2}\pi^{2}$, in
agreement with the free-fermion theory for quadratic dispersion. At the same
time, minor deviations from the predicted theoretical value
$\kappa(U=0)\approx 4J^{2}\pi^{2}$ are primarily due to the lattice dispersion
corrections, as discussed in Sec. III.2.

  * •

In the limit of small density perturbations, $\delta\rho\ll\rho_{0}$, the
viscosity term in Eq. (S111) can be neglected, and the density dynamics takes
the form $\Bigl{[}\partial_{t}^{2}-v_{\rm
eff}^{2}(U)\partial_{x}^{2}\Bigl{]}\rho(t,x)=0$, where $v_{\rm
eff}(U)=\sqrt{\kappa(U)}$ is a renormalized quasiparticle velocity. This
equation can be derived independently in the framework of the Tomonaga-
Luttinger (T-L) theory[90], which predicts the renormalized velocity of the
form

$v_{\rm eff}(U)=v_{F0}\sqrt{1+\frac{U}{2\pi v_{F0}}},$ (S112)

where $v_{F0}=\pi\rho_{0}$ is the Fermi velocity in the non-interacting limit.
Hence, given the relation between effective velocity and $\kappa(U)$, the
value of $\kappa(U)$ in this regime must be close to
$\kappa(U)=\kappa(0)\left[1+U/(2\pi^{2}\rho_{0})\right]$. Indeed, we find a
quantitative agreement between extracted values of $\kappa(U)$ and the T-L
prediction for interaction strengths in the region $U/J\lesssim 1$. However,
in the strongly interacting regime, $U/J\gg 1$, the observed values of
$\kappa$ significantly deviate from the T-L theory (see Fig. 5 in the main
text). Notably, the effective hydrodynamic model (S111) works even beyond the
linearized regime for relatively large density perturbations in both weakly-
and strongly-interacting limits.

  * •

Although the linearization of the discovered Euler equation (S111) is
equivalent to the T-L theory, to the best of our knowledge, the non-linear
term $\kappa(U)\rho\rho_{x}$ cannot be explicitly derived from the T-L theory.

  * •

Larger perturbations and longer-time evolution would have notable features
that cannot be captured by a renormalized $\kappa(U)$. Lowering the penalty
constant $\lambda_{0}$ and extending the library of candidate terms, we found
that the deviations are well-captured by the equation

$v_{t}+vv_{x}+\kappa\rho\rho_{x}=\nu v_{xx}.$ (S113)

Comparison between TEBD data and the solution of the Navier-Stokes model
(S113) is shown in Fig. S12. The viscosity term on the right-hand side
significantly improves agreement with TEBD data at long evolution times.
Despite the fact that the viscosity term violates time-reversal invariance of
the effective Navier-Stokes model, it does not conflict with unitary dynamics,
due to the presence of relaxation processes. Short-range interactions result
in the entropy flow from small-scale fluid cells to large-scale structures, by
analogy with classical hydrodynamics. More accurate analysis shows that the
universal value of viscosity emerges only in late-time dynamics, see Sec.
III.9. Notably, viscosity effects are beyond linear T-L theory.

  * •

We verified that both Euler and Navier-Stokes PDEs (S111, S113) remain
accurate for various initial states, see Fig. S13.

Figure S12: (a, b) TEBD data for the evolution of fermion density and velocity
in the spinless Fermi-Hubbard model (S109). (c) Difference between the
solution of the reconstructed hydrodynamic PDE (“Navier-Stokes model”, Eq.
(S113)) and the TEBD data ($U/J=4$). In our simulations, we fixed the value of
the chemical potential at $\mu/J=-2$. Figure S13: Hydrodynamics in the Fermi-
Hubbard model (S109) in the strongly interacting regime, $U/J=4$. (a, b)
Evolution of fermion density $\rho(t,x)$ and velocity $v(t,x)$ for an
asymmetric double peak initial density profile. This plot demonstrates that
the discovered Euler and Navier-Stokes equations (S111, S113) remain valid for
various initial conditions.

###  III.9 Analysis of an emergent Navier-Stokes equation and the viscosity
coefficient

In this subsection, we discuss an emergent viscosity term in the hydrodynamics
of the 1D Fermi-Hubbard model in more detail and point out limitations of our
analysis.

In order to break integrability of the spinless Fermi-Hubbard model and
consider a more generic case, we introduce next-nearest-neighbor fermion-
fermion interactions: $H=H_{\rm f}+V_{\rm int}$. Here $H_{\rm f}$ is a tight-
binding Hamiltonian and fermion-fermion interaction has the form $V_{\rm
int}=U\sum_{i}n_{i}n_{i+1}+U_{2}\sum_{i}n_{i}n_{i+2}$. Parameters $U$ and
$U_{2}$ are the nearest-neighbor and the next-nearest-neighbor couplings,
respectively. Remarkably, the hydrodynamic equation (S111) remains valid in
the presence of next-nearest-neighbor interactions: see comparison between
TEBD data and the solution of the PDE in Fig. S14. In this case, the pressure
renormalization coefficient becomes a function of both coupling parameters:
$\kappa(U,U_{2})$. Hence the hydrodynamic model (S111) is applicable to a wide
range of models of interacting fermions (e.g. the Fermi-Hubbard model with
nearest-neighbor interactions, next-nearest-neighbor interactions, etc..).
Discovered Eq. (S111) for an ideal Euler fluid works well for short evolution
times. However, we found that for longer evolution times, the viscous Navier-
Stokes model (S113) becomes more accurate.

Figure S14: Fermion hydrodynamics in the Fermi-Hubbard model with nearest-
neighbor interaction strength $U$ and next-nearest-neighbor interaction
strength $U_{2}$. Panels (a, c) show data for fermion density, and panels (b,
d) show data for the velocity at the final evolution time, $t_{f}=60$.
Parameters of the Fermi-Hubbard Hamiltonian are (a, b) $U/J=4$, $U_{2}/J=0$,
(c, d) $U/J=2$, $U_{2}/J=-1$. The blue and red solid lines correspond to the
discovered Euler-like equation (S111) and Navier-Stokes-like equation (S113),
respectively.

The viscosity term in Eq. (S113) discovered by our algorithm has an important
qualitative role: it prevents the formation of the gradient catastrophe. We
found that the viscosity coefficient obtained from TEBD data $\nu$ could drift
with time, see Fig. S15. We extract the time dependence of the viscosity
coefficient by performing our PDE-reconstruction within a sliding temporal
window $[t,\,t+0.2T]$, where $T$ is the total evolution time. At the beginning
of the evolution, the effect of the viscosity term is negligible, and the
effective dynamics is well-described by the inviscid Euler equation. The
viscosity coefficient extracted from TEBD data is close to zero at the start
of the evolution and grows with evolution time approaching a fixed value
$\nu\sim J$ at long evolution times, $t\to\infty$. We checked that the
asymptotic value for the viscosity coefficient remains stable to variations of
the initial state (changing the amplitude of the external potential), see Fig.
S15. Generally, the value of the viscosity coefficient will depend on the
values of interaction constants $U$ and $U_{2}$. In the non-interacting limit
$U,U_{2}\to 0$, the viscosity coefficient must vanish, $\nu\to 0$. In our
simulations, the magnitude of the effective viscosity term $\nu v_{xx}$ is
much smaller compared to the dominant pressure term $\kappa\rho\rho_{x}$,
which affects the precision of the extracted value of $\nu$. The saturation of
the viscosity coefficient $\nu(t\to\infty,U,U_{2})\to const$ can be
interpreted as an onset of local equilibration in the interacting fermion gas.

Figure S15: Time dependence of the extracted viscosity coefficient
$\nu(U,U_{2},t)$ in the extended Fermi-Hubbard model with nearest-neighbor
interactions ($U/J=U_{2}/J=2$). Viscosity is extracted within the sliding time
window $[t_{start},t_{start}+0.2T]$, where $T$ is the total evolution time. We
present data for five initial conditions corresponding to diferent values of
the amplitude of the Gaussian potential $V_{0}$. At large evolution times, the
viscosity coefficient saturates to a constant value
$\nu(t\to\infty,U,U_{2})\sim J$.

##  IV Details of numerical simulations

In the present Section, we provide additional details of numerical simulations
for data generation. In Section IV.1, we analyze the errors in the
coefficients of recovered PDEs and the dependence of the errors on the
spatiotemporal resolution of the input data.

Simulations of dynamics in the 1D XX spin chain and the non-interacting
fermion chain were performed by exact diagonalization of the single-particle
density matrix, $\mathcal{G}_{ij}(t)=\langle
c_{i}^{\dagger}(t)c_{j}(t)\rangle$. In the cases where single-magnon dynamics
in the XXZ spin chain was considered, we exactly solved the Schrödinger
equation in the single-magnon sector of the full Hamiltonian and computed the
observables.

TEBD simulations for the dynamics of the domain-wall initial state in the XXZ
model and for the dynamics in the interacting Fermi-Hubbard model were
performed with the TenPy package [91]. The matrix-product-state (MPS) bond
dimension was set to $\chi=200$. We checked that an increase of the MPS bond
dimension to $\chi=300$ resulted in a small change in the values of the
recovered coefficients of the hydrodynamic PDEs ($<0.1\%$).

###  IV.1 Error analysis

In this subsection, we would like to comment on the sources of error in our
PDE-reconstruction procedure.

The primary sources of error encountered in PDE-learning from numerical
simulations are (i) errors in numerical schemes for the evaluation of high-
order derivatives from data, (ii) numerical errors in the dataset (e.g.
truncation errors in TEBD simulations), (iii) physical-model errors
originating from higher-order corrections to the approximate PDE: corrections
beyond the hydrodynamic approximation, high-order terms in gradient expansion,
etc… A substantial amount of noise in the data can confuse the sparse
regression algorithm, thereby introducing spurious terms and/or shifting the
values of the extracted coefficients. Of course, experimental data is usually
more noisy as compared to numerical simulations, thus affecting the
reliability of the recovered equations.

Below we discuss the role of the spatiotemporal resolution on the quality of
PDE reconstruction, see Table SIII. We found that leading semiclassical terms
are robustly identified with our method, even for very “pixelated” data, see
Fig. S16.

Table SIII: Dependence of PDE reconstruction performance on the spatiotemporal
resolution $(N_{t},N_{x})$ of the dataset for the quench problem in the non-
interacting fermion gas [see Sections III.1 and III.2]. While changing the
resolution of the dataset, we keep the spatiotemporal extent $(T,L)$ fixed.
The candidate terms in the Euler equation $v_{t}=G(\cdot)$ are preselected
based on $(P,T)$ symmetry, resulting in $M=20$ terms, see Table SII. The
leading WKB terms $v_{t}+vv_{x}+\pi^{2}\rho\rho_{x}+\ldots=0$ were correctly
identified by the BruteForce algorithm for the entire range of
$(N_{t},\,N_{x})$ presented. When decreasing the number of spatiotemporal
points some of the correction terms originating from the tigh-binding
dispersion were misidentified, the number of misidentified terms
$(b_{1}\rho^{3}\rho_{x},b_{2}v^{2}\rho\rho_{x},b_{3}\rho^{2}vv_{x})$ is shown
in the entries of the table.

Statistical uncertainty in the values of regression coefficients could be
estimated if the covariance matrix of the error term $\epsilon=y-A\xi$ is
known, $\Sigma=\mathbb{E}[\epsilon\,\epsilon^{T}]$. However, for a given
dataset, the residual term $\epsilon$ is a fixed vector, and therefore its
covariance is not known. As an alternative approach to estimating the
statistical error in $\xi$, we randomly select a subset of rows of the
regression vector and the regression matrix: we split data in 10 batches of
equal size each containing $10\%$ of the data points. Next we find the
regression vector $\xi$ for each batch and estimate the uncertainty in the
regression coefficients as an element-wise standard deviation of the values of
$\xi$ across batches. The resulting statistical error, as well as the
empirical error $|\xi-\xi_{true}|$ (deviation of recovered coefficients from
the exact theoretical values), is shown in Fig. S17, where $\xi_{true}$ are
the theoretically expected coefficients. Comparing left and right panels of
Fig. S17, we see that the statistical error remains significantly lower
compared to the empirical error. Therefore, the uncertainties in the values of
the regression coefficients are primarily systematic and generated by noise in
the numerically calculated derivatives, as well as by the mutual bias from
higher-order nonlinear terms.

Additionally, we analyze the dependence of the reconstruction error on the
choice of the upper cutoff of the evolution time $t_{f}$ in the dataset, see
Fig. S18.

Figure S16: Fermion hydrodynamics (non-interacting fermions): PDE-learning
from low resolution data for (a) density and (b) velocity. When significantly
decreasing data resolution to $(N_{t},N_{x})=(10,10)$, our algorithm was still
able to identify the correct form of the hydrodynamic PDE, although the error
in the values of the coefficients became significant: $\rho_{t}+1.39\,\rho
v_{x}+1.14\,v\rho_{x}=0$ ($\lambda_{0}=10^{-5}$), and
$v_{t}+1.94\,vv_{x}+10.3\rho\rho_{x}=0$ ($\lambda_{0}=10^{-4}$).  Figure S17:
Dependence of the empirical error $|\xi-\xi_{true}|$ on the dataset temporal
resolution $N_{t}$ in recovered coefficients for the $vv_{x}$ and
$\rho\rho_{x}$ terms in the PDE
$v_{t}=\xi_{vv_{x}}vv_{x}+\xi_{\rho\rho_{x}}\rho\rho_{x}$. Here
$\xi_{true}=(\xi_{vv_{x}},\xi_{\rho\rho_{x}})=-(1,\pi^{2})$ corresponds to the
theoretical values of the coefficients for the leading WKB terms. Scattered
points show empirical error in the coefficients (a) $\xi_{vv_{x}}$ and (b)
$\xi_{\rho\rho_{x}}$ for individual random batches of subsampled data. This
plot illustrates that the uncertainty in the reconstruction of coefficients in
nonlinear PDEs is primarily systematic, since the spread of the error
distribution in batches is much smaller than the mean error for a fixed value
of $N_{t}\in[25,1000]$. The number of spatial points was fixed at
$N_{x}=1000$. The input dataset corresponds to Fig. S9.  Figure S18: (a)
Dependence of the reconstruction error of the coefficients of the leading WKB
terms, $v_{t}=\xi_{vv_{x}}vv_{x}+\xi_{\rho\rho_{x}}\rho\rho_{x}$, on the
choice of the training time window. The reconstruction error is defined as
$|\xi-\xi_{true}|$. (b) Input data for the evolution of fermion density,
$\rho(t,x)$. The solid horizontal black lines in (b) show the upper cutoff
$t_{f}$ of the training window $t\in[0,t_{f}]$. At large values $t_{f}\sim T$,
the reconstruction error $\xi_{vv_{x}}$ and $\xi_{\rho\rho_{x}}$ grows due to
high spatial gradients of the density and the velocity (the “gradient
catastrophe”). The gradient catastrophe is a well-known feature of the
semiclassical description of fermionic dynamics [69].

##  V Details of PDE-learning from experimental data: interacting bosons on an
atom chip

In this Section, we provide details of PDE-learning from experimental data,
including details of data post-processing/interpolation.

We process the experimental data as follows. The data from Ref. 74 contains
density profiles of the atomic cloud. The spatiotemporal resolution of the
original data is $[565\times 7]$. The seven experimental snapshots at
different evolution times correspond to the time range $t\in[0,85]$ ms. The
total length of the 1D atomic cloud is $L\sim 10^{3}\mu m$. The original post-
processed experimental data contains high-frequency spatial noise resulting in
negative values of the measured density $\rho(t,x)$. In order to reconstruct
continuous equations, we remove the spatial noise and increase the time
resolution. We first remove the high-frequency component of the spatial noise
by applying a Gaussian filter with the variance parameter
$\sigma_{x}/L=2.5\times 10^{-2}$ for each of the seven snapshots. Next, to
suppress the remaining low-frequency noise, we apply the Savitsky-Gollay [92]
filter with a sliding window of length 41 and a polynomial of order 2.
Finally, we perform cubic 2D interpolation of the resulting data in order to
increase resolution along the temporal dimension, which results in the final
dataset with spatiotemporal dimensions $[200\times 200]$. The data for the
particle velocity was reconstructed by leveraging the continuity equation, Eq.
(S98)—the result is shown in Fig. S19(b).

The interpolation process may impact the precision of the learning process.
For example, the interpolated particle velocities at the start of the quench
(i.e. $t=0$) are non-zero, see Fig. S19 (middle panel). In contrast, based on
the experimental quench protocol, we expect zero particle velocity immediately
after the quench of the confining double-well potential. This mismatch is a
byproduct of insufficient temporal resolution of the original data. This
issue, however, has a limited adverse impact since the inferred velocities at
the initial times are not too large. Also, while solving the inferred PDE
forward, Fig. S19 (c), we plugged $\rho(t=0,x)=\rho_{data}(0,x)$ and
$v(t=0,x)=0$ as the initial conditions.

Figure S19: (a, b) Post-processed experimental data [from Ref. 74] for (a)
density $\rho(t,x)$ and (b) particle current $j(t,x)=\rho v$ corresponding to
boson cloud expansion on an atom chip from a double well potential. Here
$t_{max}=85$ms is the maximum evolution time in the experimental dataset. The
particle-current data was reconstructed from the density data $\rho(t,x)$ by
utilizing the continuity equation. (c) Atom cloud density (in arbitrary
units): original unprocessed data (thick fading line), post-processed data
(dashed line), and the solution of the inferred PDE (solid line, Eq. (20) in
the main text). The individual density profiles at different times were
shifted vertically relative to each for visualization purposes (black
horizontal lines correspond to the origin of the vertical axis). The seven
density profiles in (c) correspond to evolution times labeled in (a).