# Electrically switchable giant Berry curvature dipole in silicene, germanene
and stanene

Arka Bandyopadhyay  Solid State and Structural Chemistry Unit, Indian
Institute of Science, Bangalore 560012, India.  Nesta Benno Joseph  Solid
State and Structural Chemistry Unit, Indian Institute of Science, Bangalore
560012, India.  Awadhesh Narayan<EMAIL_ADDRESS>Solid State and
Structural Chemistry Unit, Indian Institute of Science, Bangalore 560012,
India.

###### Abstract

The anomalous Hall effect in time-reversal symmetry broken systems is
underpinned by the concept of Berry curvature in band theory. However, recent
experiments reveal that the nonlinear Hall effect can be observed in non-
magnetic systems without applying an external magnetic field. The emergence of
nonlinear Hall effect under time-reversal symmetric conditions can be
explained in terms of non-vanishing Berry curvature dipole arising from
inversion symmetry breaking. In this work, we availed realistic tight-binding
models, first-principles calculations, and symmetry analyses to explore the
combined effect of transverse electric field and strain, which leads to a
giant Berry curvature dipole in the elemental buckled honeycomb lattices –
silicene, germanene, and stanene. The external electric field breaks the
inversion symmetry of these systems, while strain helps to attain an
asymmetrical distribution of Berry curvature of a single valley. Furthermore,
the topology of the electronic wavefunction switches from the band inverted
quantum spin Hall state to normal insulating one at the gapless point. This
band gap closing at the critical electric field strength is accompanied by an
enhanced Berry curvature and concomitantly a giant Berry curvature dipole at
the Fermi level. Our results predict the occurrence of an electrically
switchable nonlinear electrical and thermal Hall effect in a new class of
elemental systems that can be experimentally verified.

##  I Introduction

The appearance of Hall current is invariably contingent on the breaking of
time-reversal symmetry (TRS) in the linear response regime Klitzing  _et al._
(1980). TRS in a material can be broken by an external magnetic field, or
suitable magnetic dopants. However, in recent experiments, _nonlinear_ Hall
effects (NHE) Du  _et al._ (2021); Ortix (2021) have been detected in non-
magnetic transition-metal dichalcogenides (TMDs) under time-reversal-symmetric
conditions Ma  _et al._ (2019); Kang  _et al._ (2019); Huang  _et al._ (2020).
In their seminal work, Sodemann and Fu Sodemann and Fu (2015) have explored
the quantum origin of this nonlinear response by introducing an intrinsic
effect of the dipole moment of the Berry curvature. This Berry curvature
dipole (BCD) can be observed in time-reversal invariant systems, but its non-
zero value is strictly protected by the breaking of inversion symmetry of the
crystal. Moreover, symmetry-based indicators are also crucial in determining
the strength of BCD in a noncentrosymmetric system. For example, a uniaxial
strain reduces the symmetry of TMDs and gives rise to enhanced BCD You  _et
al._ (2018); Zhou  _et al._ (2020); Son  _et al._ (2019). Similar enhancements
are also observed for few-layer TMDs Ma  _et al._ (2019); Kang  _et al._
(2019); Huang  _et al._ (2020); Joseph and Narayan (2021), where the lowering
of symmetry is the result of the stacking of monolayers. It has further been
observed that the electric field can efficiently tune the anomalous NHE in
low-symmetry TMDs Zhang  _et al._ (2018a). Moreover, a pressure-driven
topological phase transition in three-dimensional bismuth tellurium iodine
(BiTeI) with a strong Rashba effect is also predicted to lead to a large BCD
Facio  _et al._ (2018). In principle, the BCD-induced NHE can be considered as
a second-order response to the electric field in the system’s plane. The
combined effect of this in-plane electric field and BCD is responsible for
several exotic physical properties, such as giant magneto-optical effects Liu
and Dai (2020), orbital valley magnetization Son  _et al._ (2019), non-linear
Nernst effects Yu  _et al._ (2019); Zeng  _et al._ (2019), and thermal Hall
effect Zeng  _et al._ (2020).

In general, large Berry curvature segregation occurs at the Brillouin zone
(BZ) points, where two bands nearly touch each other. The shape of the Bloch
states rapidly modifies near such narrow-gap points of the BZ. Therefore,
massive tilted Dirac cones Sodemann and Fu (2015); Nandy and Sodemann (2019);
Du  _et al._ (2019) or Weyl cones Kumar  _et al._ (2021); Zeng  _et al._
(2021); Matsyshyn and Sodemann (2019); Singh  _et al._ (2020); Roy and Narayan
(2021) are the natural choices for realizing sizeable BCD. In these systems,
the BCD and corresponding NHE systematically provide the geometrical
information of Bloch wavefunctions even under TRS. Moreover, Battilomo et al.
have revealed that the Fermi surface warping triggers appreciable BCD in
uniaxially strained monolayer and bilayer graphene Battilomo  _et al._ (2019).
It is worth mentioning that the magnitude of BCD in the warped graphene
systems is comparable with that of the TMDs. Further, the merging of Dirac
points near the Fermi level can also lead to non-zero BCD even in the absence
of any tilt or warping term Samal  _et al._ (2021). The BCD and the
unconventional NHE drive the understanding of topological physics and quantum
transport phenomena to the nonlinear domain. This generalization opens up many
exciting prospects for direct applications, such as nonlinear photocurrents Xu
_et al._ (2018) and terahertz radiation detection Zhang and Fu (2021).

In a pioneering work, Kane and Mele Kane and Mele (2005) first explored the
fascinating quantum spin Hall effect (QSHE) in graphene. However, in reality,
the QSHE in graphene is not experimentally accessible because of the
negligible strength of spin-orbit coupling (SOC). The strength of SOC largely
determines the occurrence of helical edge states with well-defined spin
texture in topological insulators. In real systems, the requirement of large
SOC is considerably fulfilled by the experimental realization of ‘graphene
counterparts’ – silicene, germanene, and stanene Balendhran  _et al._ (2015);
Vogt  _et al._ (2012); Derivaz  _et al._ (2015); Zhu  _et al._ (2015).
Silicene, germanene, and stanene are two-dimensional (2D) Dirac materials
having a buckled honeycomb geometry. The buckling can be exploited by
employing a transverse electric field that can tune the electronic band
structure, particularly, the band gap Mak  _et al._ (2009). Furthermore, an
electric field driven topological phase transition from QSH to a normal
insulating (NI) state is a primary characteristic of these systems Ezawa
(2015); Deng  _et al._ (2018); Qin  _et al._ (2021).

In this work, we discover a large and electrically switchable NHE in these
elemental buckled honeycomb lattices silicene, germanene, and stanene. Using
tight-binding calculations, in conjunction with symmetry arguments, we explore
the tunability of the BCD in these systems, particularly near the topological
phase transitions. We demonstrate a giant enhancement of the BCD near the
electric field tuned topological critical points and connect it to the
underlying variations of the Berry curvature. Our findings put forward a new
class of systems to explore nonlinear topological phenomena, and highlight an
as-yet-unexplored aspect of elemental buckled honeycomb lattices. We hope that
our work motivates experimental as well theoretical work along this front in
the near future.

##  II Methodology

Figure 1: Buckled honeycomb lattices and their topological nature. (a) The top
and side view of buckled strained honeycomb structures – silicene, germanene,
and stanene. One of the three nearest bonds is stretched (up to $\approx
2\%$). Therefore, the angle between the translation vectors increases from the
usual 60o. A transverse electric field ($\vec{E}$) adequately exploits the
buckling and breaks the inversion symmetry of the lattice. Different colors
indicate two sub-lattices for clarity. The elongated bond destroys the $C_{3}$
rotational symmetry but protects the mirror axis ($M_{x}$) in the direction
perpendicular to the stretching. (b) The Brillouin zone of the symmetry
reduced buckled honeycomb structures. The coordinates of the symmetry point
$\mathcal{K}$ are distinct from that of the original systems. (c) The phase
diagram of the perturbed system ($\approx 2\%$ stretch) resembles the pristine
structure. The blue line represents the variation of band gap with external
electric field. Appreciable spin-orbit interaction gives rise to the quantum
spin Hall (QSH) phase with a finite band gap ($\Delta\mathcal{E}$). The band
gap decreases with increasing electric field, and above a critical value
($\lvert E_{c}\rvert$), a topological phase transition from topological
insulator to normal insulator (NI) occurs. The bottom panel (d) Band
structures at different electric fields. Two spin bands touch the Fermi level
for the electric field strength $\lvert E_{c}\rvert$. At this point, the
system behaves as a semimetal with gapless bands. The systems exhibit a band
gap below and above the critical field. The spin degeneracy of the bands is
broken due to the presence of a finite electric field.

The nonlinear current in response to an oscillating electric field
$\vec{E}(t)=Re\\{\vec{E}_{0}e^{i\omega t}\\}$, with a magnitude $E_{0}$ and
frequency $\omega$ can be expressed as
$J_{\alpha}=Re\\{J_{\alpha}^{(0)}+J_{\alpha}^{(2\omega)}e^{i2\omega t}$}
Sodemann and Fu (2015). Therefore, the response current has been clearly
decomposed into a static part (rectified current)
$J_{\alpha}^{(0)}=\chi_{\alpha\beta\gamma}\vec{E}_{\beta}\vec{E}_{\gamma}^{*}$
and a double frequency oscillating part (second harmonic)
$J_{\alpha}^{(2\omega)}=\chi_{\alpha\beta\gamma}\vec{E}_{\beta}\vec{E}_{\gamma}$.
Under time-reversal symmetric conditions, the nonlinear conductivity tensor
($\chi_{\alpha\beta\gamma}$) depends on the momentum derivative of the Berry
curvature over the occupied states as follows

$\chi_{\alpha\beta\gamma}=-\epsilon_{\alpha\delta\gamma}\frac{e^{3}\tau}{2(1+i\omega\tau)}\int_{k}[d\vec{k}]f_{0}\left(\frac{\partial\Omega_{\delta}}{\partial
k_{\beta}}\right).$ (1)

Here, $e$, $\tau$, $\epsilon_{\alpha\delta\gamma}$, $f_{0}$ and
$\Omega_{\delta}$ represent the electron charge, scattering time, Levi-Civita
symbol, equilibrium Fermi-Dirac distribution, and Berry curvature component
along $\delta$, respectively with $\alpha,\beta,\gamma,\delta\in\\{x,y,z\\}$.
Here the integration is performed with respect to $[d\vec{k}]$, which has the
expression $d^{d}k/(2\pi)^{d}$ in $d$ dimensions. BCD can be defined in
reciprocal space as
$D_{\alpha\beta}=\int_{k}[d\vec{k}]f_{0}\left(\partial\Omega_{\beta}/\partial
k_{\alpha}\right)$. Particularly, in 2D materials, only the out-of-plane ($z$)
component of Berry curvature is non-vanishing, i.e.,
$\Omega_{\beta}\equiv\Omega_{z}$. In the framework of the well-known Kubo
formalism, this $\Omega_{z}(\vec{k)}$ has the following form Xiao  _et al._
(2010)

$\Omega_{z}(\vec{k})=2i\sum_{i\neq j}\frac{\langle i|{\partial\hat{H}/\partial
k_{x}}|j\rangle\langle j|{\partial\hat{H}/\partial
k_{y}}|i\rangle}{\left(\varepsilon_{i}-\varepsilon_{j}\right)^{2}},$ (2)

where $\varepsilon_{i}$ and $\varepsilon_{j}$ are the eigenenergies of the
Hamiltonian $\hat{H}$ with eigenstates $|i\rangle$ and $|j\rangle$,
respectively. The methodology discussed above for calculating BCD is
implemented in the Wannier-Berri package Tsirkin (2021), which is compatible
with the PythTB module Pyt . It is worth mentioning that the Berry curvature
of time-reversal invariant systems is an odd function of momentum, i.e.,
$\mathcal{T}^{\dagger}\Omega_{z}(-\vec{k})\mathcal{T}=-\Omega_{z}(\vec{k})$,
where $\mathcal{T}$ is the time-reversal operator. In contrast, BCD is even
under the above situation, as it satisfies
$\mathcal{T}^{\dagger}D_{\alpha\beta}(-\vec{k})\mathcal{T}=D_{\alpha\beta}(\vec{k})$.
Therefore, it is clear that BCD can manifest a significant anomalous
electronic response, NHE, even in the presence of TRS.

In order to support the tight-binding results, first-principles calculations
are carried out based on the density functional theory (DFT) framework as
implemented in the quantum espresso code Giannozzi  _et al._ (2017, 2009). A
kinetic energy cut-off of $40$ Ry is considered, using the ultrasoft
pseudopotentials Vanderbilt (1990) to describe the core electrons, including
spin-orbit coupling interactions. We used the Perdew-Burke-Ernzerhof (PBE)
form for the exchange-correlation functional Perdew  _et al._ (1996). The
Brillouin zone is sampled over a uniform $\Gamma$-centered $k$-mesh of
$8\times 8\times 1$, and the monolayers were modeled with a 15 Å vacuum along
the $z$-direction to avoid any spurious interaction between the periodic
images. To study the topological properties, maximally localized Wannier
functions (MLWFs) were computed to derive a tight-binding model from the ab-
initio calculations, with complete $s,p$ orbitals as the basis, using the
wannier90 code Mostofi  _et al._ (2014). Further calculation of
$\mathbb{Z}_{2}$ topological invariants and analysis of the edge spectra is
performed using the WannierTools code Wu  _et al._ (2018).

Figure 2: The electronic properties obtained by density functional theory
methods. (a) The band structures of 2% strained buckled honeycomb lattices
exhibit a band splitting at $\mathcal{K}$ point. (b) Wannier charge centers
for strained silicene demonstrate the nontrivial $\mathbb{Z}{{}_{2}}$ index of
the system. (c) The protected edge states have been calculated in the
nontrivial topological region of the phase space.

##  III Results and discussion

###  III.1 Model and symmetry analysis

Next, we introduce the tight-binding (TB) model Hamiltonian for buckled
honeycomb lattices used in this work. In spite of the similarity in their
basic geometry, the buckled honeycomb lattices primarily differ in their bond
lengths ($d$), hopping integrals ($t$), and lattice parameters
($a=\sqrt{3}d$). Further, the strength of SOC ($\lambda_{SO}$) and buckling
height ($l$) between the two sublattices also depend on the atomic number of
the constituent atom and vary from system to system. All the above mentioned
parameters of buckled honeycomb lattices are compared with that of graphene in
Table 1. In our study, we shall extensively use these values in Table 1 to
extract the system specific information, while also using general symmetry
arguments. In the presence of buckling, the transverse electric field
($E_{z}$) assigns different mass terms to the two sublattices. This difference
in mass terms, in turn, breaks the inversion symmetry of the system. The
generalized TB Hamiltonian Ezawa (2015) of buckled honeycomb lattices in the
presence of an external transverse electric field can be expressed as

$\displaystyle\hat{H}=$ $\displaystyle-t\sum_{\langle
ij\rangle,\sigma}c_{i,\sigma}^{\dagger}c_{j,\sigma}+i\frac{\lambda_{SO}}{3\sqrt{3}}\sum_{\langle\langle
ij\rangle\rangle,\sigma}\sigma\zeta_{ij}c_{i,\sigma}^{\dagger}c_{j,\sigma}$
$\displaystyle-l\sum_{i,\sigma}\nu_{i}E_{z}c_{i,\sigma}^{\dagger}c_{i,\sigma}.$
(3)

Here $\langle ij\rangle$ and $\langle\langle ij\rangle\rangle$ indicate
hopping between $i$ and $j$ sites up to nearest and next-nearest neighbours.
Further, $\sigma$ represents spin degrees of freedom and denotes either
$\uparrow$ ($+1$), or $\downarrow$ ($-1$) spin. The SOC and staggered
sublattice potential, $\zeta_{ij}$ $(=\pm 1$) and $\nu$ $(=\pm 1$) explicitly
depend on the direction (clockwise or anticlockwise) of hopping and type of
sublattice, respectively. Furthermore, Eq. 3 reveals that buckled honeycomb
lattices exhibit a topologically nontrivial band gap,
$\Delta\mathcal{E}=2\lambda_{SO}$, in the absence of an external electric
field Ezawa (2015). Application of the electric field reduces this value of
the band gap to zero at a topological critical point $E_{C}=\lambda_{SO}/l$,
where the system behaves like a semimetal. Beyond this point,
$\Delta\mathcal{E}$ again increases and gives rise to a topologically trivial
phase. The tight-binding parameters for Eq. 3 vary across different systems we
considered. The parameters for silicene, germanene, and stanene are compared
with graphene in Table 1. Because of the negligible spin-orbit coupling and
planar geometry, graphene does not show an electrically tunable quantum spin
Hall state.

Table 1: The structural and tight binding model parameters of silicene, germanene, and stanene are compared to that of graphene. Here, $a$, $t$, $\lambda_{SO}$, and $l$ represent the lattice parameter, hopping integral, SOC strength, and buckling heights, respectively. The values are taken from Ref. Liu  _et al._ (2011). Systems |  $a$ (Å) |  $t$ (eV) |  $\lambda_{SO}$ (meV) |  $l$ (Å)
---|---|---|---|---
Graphene | 2.46 | 2.8 | 10-3 | 0.00
Silicene | 3.86 | 1.6 | 3.9 | 0.23
Germanene | 4.02 | 1.3 | 43.0 | 0.33
Stanene | 4.70 | 1.3 | 100.0 | 0.40

To trace this topological transition, we initiate our calculations in the
presence of a non-zero electric field. The electric field helps us to attain
the essential criterion of inversion symmetry breaking for BCD. Regardless of
the breaking of inversion symmetry in buckled honeycomb lattices, we obtain a
vanishing value of BCD at different electric field strengths. This can be
understood from the crystallographic symmetry of the buckled honeycomb
lattices. The buckling in buckled honeycomb lattices primarily eliminates the
$C_{6}$, $\sigma_{h}$, and three $\sigma_{d}$ symmetry elements of the planar
honeycomb lattices ($D_{6h}$). Consequently, the system possesses only three
$C_{2}$ rotational axes perpendicular to the principal axis of symmetry.
Further, three mirror planes of the buckled system bisect the angle between
each neighboring pair of these $C_{2}$ rotational axes. Therefore, the groups
of wavevectors in buckled honeycomb lattices are $D_{3d}$ (point group of
buckled honeycomb lattices) and $D_{3}$ at symmetry points $\Gamma$ and $K$
respectively. The order of the finite group $D_{3}$ is $3\times 2=6$ with 3
rotational and 3 reflection symmetry elements. In particular, the rotational
symmetry elements are $2\pi$, $2\pi/3$, and $4\pi/3$ rotation about the
$C_{3}$ axis, while the reflection symmetry elements represent the symmetry
planes ($\sigma_{d}$) passing through the 3 medians of the equilateral
triangle. The character tables are presented in Table S1 and Table S2 of ESI
Bandyopadhyay  _et al._ . The presence of two or more mirror axes in two-
dimensional buckled honeycomb lattices (here, three) relates non-linear Hall
conductance to a null pseudovector field Ortix (2021), so that the BCD is
zero. Fundamentally, the maximum permitted symmetry in two dimensions for the
occurrence of non-vanishing BCD is a single mirror line. This symmetry
analysis motivates us to reduce the symmetry of buckled honeycomb lattices
down to a single mirror axis by minimal operations. We have achieved this by
applying a small uniaxial strain that elongates one nearest neighbour bond
($d^{\prime}$) compared to the other two ($d$) as depicted in Fig 1(a). As a
representative value we have chosen to $d^{\prime}$ to be different from $d$
by 2%. Application of such a strain does not break the inversion symmetry of
these systems but essentially reduces its rotational symmetry. In the above
case, $C_{3}$ rotational symmetry is destroyed along with two $\sigma_{d}$
reflection symmetry planes. The strained buckled honeycomb lattices possess
only one $\sigma_{d}$ symmetry element, which corresponds to $M_{x}$ mirror
symmetry plane in our case. We expect that these strained buckled honeycomb
lattices will be suitable candidates for obtaining sizable BCD in the presence
of an external transverse electric field.

Figure 3: Berry curvature and electrical field. The Berry curvatures of the
four bands near $\mathcal{K}$ point are shown for stanene under (a) 0.15 V/Å
(non-trivial phase) and (b) 0.35 V/Å (trivial phase) electric field strength.
The spin degeneracy of the bands is lifted under the electric field. A clear
flipping of the Berry curvature of two bands indicates the topological phase
transition. The color scale is shown on the right. Figure 4: Wannier charge
centers. The Wannier charge centers (WCCs) are plotted for stanene under (a)
0.15 V/Å (non-trivial phase) and (b) 0.35 V/Å (trivial phase) electric field
strengths. In the first case, a single line parallel to $k_{y}$ crosses an odd
number of charge centers in the half of the BZ, giving rise to odd
$\mathbb{Z}_{2}$ invariant. This odd $\mathbb{Z}_{2}$ confirms the non-trivial
topological state of the material. (b) In the second case, due to an even
number of crossing points, the WCC reveals that the $\mathbb{Z}_{2}$ invariant
is now even, and the phase is trivial.

The unit cell of the strained buckled honeycomb lattices is defined by a new
set of lattice vectors: $\vec{a}_{1}=\sqrt{3}d\hat{x}$ and
$\vec{a}_{2}=\sqrt{3}d/2\hat{x}+(3d/2+2d/100)\hat{y}$. As a result, the angle,
$\theta$, between lattice parameters slightly increases than the usual
$60^{\circ}$ of the pristine case. In the reciprocal space, the strain changes
the high symmetry point $K$ $(0.3333,0.6667)$ to an equivalent point
$\mathcal{K}$ $(0.3404,0.6702)$ of the BZ as shown in Fig. 1(b). Moreover, the
band minimum or maximum of strained buckled honeycomb lattices occur at this
new point $\mathcal{K}$, near which the band dispersion is linear. Similar to
the unperturbed case, the band gap of the slightly strained buckled honeycomb
lattices is $2\lambda_{SO}$ in the absence of an electric field. Furthermore,
the applied electric field controls the band gap as $\Delta\mathcal{E}\approx
2l|E_{z}-\lambda_{SO}/l|$ as indicated in Fig. 1(c). Note that the strained
systems are semimetallic ($\Delta\mathcal{E}$ is zero) at a critical electric
field $E_{C}\approx\lambda_{SO}/l$. The low-energy effective Hamiltonian
$\hat{H}_{\mathrm{eff}}$ around $\mathcal{K}$ can be written as follows

$\hat{H}_{\mathrm{eff}}\approx\begin{bmatrix}\mu+h_{11}&v_{x}k_{x}+iv_{y}k_{y}\\\
v_{x}k_{x}-iv_{y}k_{y}&\mu-h_{11}\end{bmatrix},$ (4)

where $v_{x}=v_{y}=v_{F}$ is the Fermi velocity. The matrix elements $h_{11}$
and $\mu$ are related to the Pauli spin matrix $\sigma_{z}$ by the relations
$h_{11}=-\lambda_{SO}\sigma_{z}$ and $\mu=lE_{z}\sigma_{z}$. The three-
dimensional band structures obtained using $\hat{H}_{\mathrm{eff}}$ around
$\mathcal{K}$ are illustrated in Fig. 1(d). As expected, the band structures
of strained buckled honeycomb lattices exhibit band opening on both sides of
the critical field $E_{C}$ similar to the pristine case (see also, Fig. S1 of
ESI Bandyopadhyay  _et al._ ).

The DFT band structures for all the three strained systems, given in Fig.
2(a), also indicate that the strained buckled honeycomb lattices invariably
exhibit a band gap at the $\mathcal{K}$ point. Further, the nontrivial
topological nature of the band gap has been confirmed by the Wannier charge
centers plots ($\mathbb{Z}{{}_{2}}=1$) and symmetry protected edge states
given in Fig. 2(b) and Fig. 2(c), respectively. These results are in excellent
agreement with our tight-binding model findings.

###  III.2 Berry curvature

We then calculated the Berry curvature using Eq. 2, to understand the
topological aspects of the above-mentioned gapped states of our strained
buckled honeycomb lattices. It has been observed that the sign of Berry
curvature of two spin states near the Fermi level, i.e., conduction band (CB)
and valence band (VB) flips when electric field strength crosses the critical
value $E_{C}$ Zhang  _et al._ (2018a). A clear example of this flipping of
Berry curvature for strained stanene with $E_{C}\approx 0.2499$ V/Å is shown
in Fig. 3. It is evident from Fig. 3(a) and (b) that the Berry curvature
distribution of the CB and VB are distinct for the electric field strength
0.15 V/Å and 0.35 V/Å. In particular, the Berry curvature of the VB (CB) is
negative (positive) under 0.15 V/Å electric field, while the same is positive
(negative) for the field strength 0.35 V/Å. In a similar vein, we have also
observed the Berry curvature exchange between the VB and CB for strained
silicene ($E_{C}\approx 0.0170$ V/Å) and strained germanene ($E_{C}\approx
0.1303$ V/Å). Note that the Berry curvature of all the systems tends to
diverge while reaching the critical point from both sides. All the results
mentioned above are certainly the signatures of a topological phase transition
of the buckled buckled honeycomb lattices.

To confirm the presence of distinct topological phases in strained buckled
honeycomb lattices, the $\mathbb{Z}_{2}$ invariant has been evaluated from the
Wannier charge centers (WCC) Yu  _et al._ (2011); Soluyanov and Vanderbilt
(2011). Fig. 4(a) indicates that a straight line parallel to $k_{y}$
intersects WCC an odd number of times ($\mathbb{Z}_{2}=1$) in half of the BZ
of strained stanene under 0.15 V/Å electric field. On the other hand, the same
has an even number of intersections ($\mathbb{Z}_{2}=0$) for the case of 0.35
V/Å electric field [Fig. 4(b)]. The $\mathbb{Z}_{2}$ invariant establishes
that the external electric field drives the strained stanene from a non-
trivial topological state to a trivial one. Similar results have also been
observed for strained silicene and germanene structures. Therefore, the
electric field dependent topological phase transition in BHLs is robust to
this applied strain. On that account, the required criteria for tracing the
topological phase transition in BHLs by BCD are entirely fulfilled.

Figure 5: BCD and its variation with electric field. (a) The variation of BCD
with the energy exhibits a peak near the Fermi level for all the systems.
Here, BCD for the silicene under 0.0170 V/Å electric field, germanene under
0.1303 V/Å electric field, and stanene under 0.2499 V/Å electric field are
given as an example. Only the $D_{xz}$ component of the BCD tensor is non-zero
from symmetry-based indicators. (b) The variation of BCD at the Fermi level
($D_{xz0}$), with the electric field for silicene, germanene, and stanene,
shows a peak around the corresponding topological critical point. The sign of
BCD also changes on either side of the topological transition.

###  III.3 Berry curvature dipole

This rapid change with electric field in sign and magnitude of Berry curvature
in momentum space strongly suggests the possibility of large and tunable BCD
in strained buckled honeycomb lattices. We next explore this aspect in these
systems by calculating the BCD in the presence of an external electric field.
As we noted, the strained buckled honeycomb lattices possess only one mirror
symmetry $M_{x}$. In the presence of $M_{x}$ symmetry, $k_{x}$ and $k_{y}$
transform as odd and even parameters, respectively. $\Omega_{z}$ is the only
non-zero component of Berry curvature in these 2D systems, which is odd in
momentum. The above observations immediately show that the gradient of
$\Omega_{z}$ along $k_{x}$, or, $d_{xz}=\partial\Omega_{z}/\partial k_{x}$
will be the only even term. In other words, the momentum integrated $d_{xz}$,
i.e., $D_{xz}$ will be the only non-zero BCD tensor component in strained
buckled honeycomb lattices. We find that all these strained buckled honeycomb
lattices exhibit a substantial BCD in the presence of electric field induced
inversion symmetry breaking. Moreover, we have discovered a generic feature of
$D_{xz}$ in these systems – the contribution of $D_{xz}$ is enormous near the
Fermi level for all the strained buckled honeycomb lattices as shown in Fig.
5(a). This giant BCD value can be well-explained by the high concentration of
Berry curvature near the Fermi level. Furthermore, the sign of BCD is found to
be reversed when the direction of electric field, $E_{z}$, is flipped. A Fermi
smearing of 40 K is considered for quantitative discussions. The band gap of
these systems gradually decreases to zero at $E_{c}$ and then again increases
linearly. This band gap variation results in a maximum BCD response near the
topological transition point, as we found. In particular, the maximum BCD of
strained germanene at the Fermi level is $D_{xz0}=239.35$ Å, which is larger
than that of silicene ($D_{xz0}=217.16$ Å). In contrast, the maximum BCD of
stanene $D_{xz0}=179.73$ Å does not follow the above trend of with increasing
SOC. Similar results have been reported for Weyl semimetals TaAs, TaP, NbAs,
and NbP Zhang  _et al._ (2018b). In comparison, the maximum value of BCD for
other 2D materials, such as strained monolayer graphene Battilomo  _et al._
(2019), bilayer graphene Battilomo  _et al._ (2019), WTe2 Zhang  _et al._
(2018a), and MoTe2 Zhang  _et al._ (2018a) are approximately 0.01 Å, 10 Å, 60
Å, and 80 Å, respectively. Therefore, strained buckled honeycomb lattices can
provide an intriguing platform to achieve large and tunable BCD. In
particular, dual-gated, encapsulated devices can be fabricated based on the
strained ($\sim$ 2%) buckled honeycomb lattices with controllable chemical
potential and transverse electric field setup for realizing BCD at low
temperature ($\sim$ 100 K). The alternating electric field applied to the
device will result in a non-linear voltage, with doubled frequency, that can
be measured using a sensitive lock-in amplifier.

The diverging nature of BCD at topological transition point can be well
understood from the following discussion. It is evident from Eq. 4 that spin-
orbit coupling gives rise to massive Dirac cones by opening a gap
$\Delta\mathcal{E}\equiv\Delta\approx 2l|E_{z}-\lambda_{SO}/l|$ in the energy
spectra. In the case of these strained buckled honeycomb lattices, another
term proportional $k_{y}$ ($\lambda k_{y}$) can be included to address the
non-isotropic dispersion arising from the strain. The term proportional to
$\lambda$ gives the anisotropic velocity, depending on the applied strain. It
is worth noting that this $\lambda$ term is responsible for the non-zero value
of the Berry curvature dipole. Further, two non-equivalent massive Dirac cones
are related by the time-reversal symmetry. Hence, only the out-of-plane
component of Berry curvature is non-zero, which is segregated in these two
valleys with a different sign. Furthermore, the small external strain is
responsible for a perturbed Berry curvature distribution and assigns different
weights to it in the Brillouin zone. We can write a low energy model
Hamiltonian for the system considering the states near the Fermi level as

$H=\left(v_{x}k_{x}\sigma_{y}-\tau
v_{y}k_{y}\sigma_{x}\right)+\Delta\sigma_{z}+\tau\lambda k_{y}.$ (5)

Here valley index $\tau=\pm$, $\sigma_{x}$ and $\sigma_{y}$ are Pauli spin
matrices. The Hamiltonian given above allows only the $M_{x}$ crystal
symmetry, where $k_{x}\rightarrow-k_{x}$ symmetry is preserved. The dispersion
relation obtained using Eq. 5 is obtained as

$\varepsilon_{\pm}(k_{x},k_{y})=\tau\lambda
k_{y}\pm\left(\Delta^{2}+v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2}\right)^{1/2}.$
(6)

Now, the Berry curvature $\Omega_{z}(\vec{k})$ can be calculated using Eq. 2.
For buckled honeycomb lattices $\Omega_{z}(\vec{k})$ can be expressed as

$\Omega_{z}(k_{x},k_{y})=\pm\frac{1}{2}\frac{\tau
v_{x}v_{y}\Delta}{\left(\Delta^{2}+v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2}\right)^{3/2}}.$
(7)

The BCD is related to the Berry curvature by
$D_{\alpha\beta}=\int_{k}[d\vec{k}]f_{0}\left(\partial\Omega_{\beta}/\partial
k_{\alpha}\right)$. It is then straightforward to write down the expression of
Berry curvature dipole mentioned above in terms of Berry curvature dipole
density $d_{\alpha\beta}$ as follows

$D_{\alpha\beta}=-\sum_{n}\int\frac{\partial
f_{n}(\vec{k})}{\partial\varepsilon_{n}(\vec{k})}v_{\alpha}\Omega_{n\beta}(\vec{k})d[\vec{k}]=\int\frac{\partial
f_{n}(\vec{k})}{\partial\varepsilon_{n}(\vec{k})}d_{\alpha\beta}d[\vec{k}].$
(8)

In the above, $f_{n}(\vec{k})$ is the distribution function and the velocity
$v_{\alpha}$ can be obtained by $\partial\varepsilon_{n}(\vec{k})/\partial
k_{\alpha}$. The BCD density then has the following expression

$d_{\alpha\beta}=-\sum_{n}v_{\alpha}\Omega_{n\beta}(\vec{k}).$ (9)

The partial differentiation gives the delta function in the low energy limit,
which indicates that Berry curvature dipole is a Fermi surface property. We
note that the band gap of the system is
$E_{g}=2\left(\Delta^{2}+v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2}\right)^{1/2}$.
From the above discussion it is evident that

$\lim_{E_{g}\to 0}d_{xz}=\frac{1}{2}\lim_{E_{g}\to
0}\sum_{n}\frac{v_{\alpha}\tau
v_{x}v_{y}\Delta}{\left(\Delta^{2}+v_{x}^{2}k_{x}^{2}+v_{y}^{2}k_{y}^{2}\right)^{3/2}}\rightarrow\infty.$
(10)

Therefore, we find that the BCD density predominantly varies as $\sim$
$1/\Delta^{2}$ and diverges near the topological critical point where the band
gap closes.

###  III.4 Nonlinear thermal Hall effect

Similar to the electrical Hall effect, the thermal Hall current also vanishes
under time reversal symmetric condition in the linear response regime.
However, in case of the strained honeycomb lattices the perturbed distribution
function gives rise to nonlinear thermal Hall effect. The thermal Hall
current, $j_{a}^{T}$, in the nonlinear regime can be obtained using Boltzmann
equation as given below

$j_{a}^{T}=-\kappa_{abd}\Delta_{b}T\Delta_{d}T.$ (11)

Here $\Delta T$ represents temperature gradient and the subscripts {$a,b,d$}
$\in$ {$x,y,z$} stand for the components. The coefficient of nonlinear thermal
Hall effect is denoted by $\kappa_{abd}$. From the symmetry analysis it is
clear that the temperature gradient along a direction normal to $M_{x}$ gives
rise to nonlinear thermal Hall effect in the perpendicular direction. The
intrinsic contribution of the nonlinear thermal Hall coefficient,
$\kappa^{BCD}_{abd}$, comes from the BCD Zeng  _et al._ (2020) given as
follows

$\kappa^{BCD}_{abd}=\zeta T^{2}G^{1}_{0}(\mu)+\mathcal{O}[T^{4}],$ (12)

where $\zeta=\frac{7\tau_{0}\pi^{4}k_{B}^{4}}{14\hbar^{2}}$ is a constant
under the constant relaxation time ($\tau_{0}$) approximation. The chemical
potential dependent parameter $G^{1}_{0}(\mu)=\partial G_{0}(\mu)/\partial\mu$
is related to the BCD by

$G_{\alpha\beta}(\mu)=\int[dk]\delta\left(\varepsilon-\varepsilon(k)\right)\Omega_{\alpha}(k)\frac{\partial\varepsilon(k)}{\partial
k_{\beta}}.$ (13)

In presence of disorder, relaxation time $\tau_{0}$ can be written as
$\tau_{0}^{-1}=nV_{0}^{2}(\mu^{2}+3\Delta^{2})/4\hbar\mu v_{x}v_{y}$, where,
$n$ and $V_{0}$ are the concentration and strength of disorder, respectively
Zeng  _et al._ (2020). In the case of our buckled honeycomb lattices, the
$yxx$ component of nonlinear thermal Hall coefficient $\kappa^{BCD}$ can be
written as

$\kappa^{BCD}=\frac{7\pi^{3}k_{B}^{2}\lambda\Delta
v_{x}v_{y}(\mu^{2}-2\Delta^{2})}{5n\hbar^{2}V_{0}^{2}\mu^{4}(\mu^{2}+3\Delta^{2})}k_{B}T^{2}.$
(14)

We have studied the variation of the nonlinear Hall coefficient for the
strained buckled honeycomb lattices (here, shown for stanene) with chemical
potential, for different temperatures, presented in Fig. 6. We have considered
representative values of the different parameters for our buckled honeycomb
systems: $v_{x}=v_{y}=v_{F}=at=6.11$ eV Å (velocity for stanene given in Table
1), $\Delta=0.1$ eV (band gap of stanene given in Table 1), $\lambda=0.1v$,
$nV_{0}=100$ eV2 Å2. Therefore, starting from the situation where the Fermi
level touches the bottom of conduction band ($\mu\approx 0.1$ eV), the
chemical potential can be tuned using a gate voltage. We find that for all
temperatures, the nonlinear thermal Hall coefficient changes its sign and
exhibits a peak at $\mu>0.1$ eV for stanene. Further, the coefficient
invariably goes to zero for large value of $\mu$. The peak value of nonlinear
thermal Hall coefficient increases with increasing temperature. In a nutshell,
BCD can give rise to the nonlinear thermal Hall conductivity, which varies
quadratically in the temperature difference.

Figure 6: Nonlinear thermal Hall conductivity. The variation of BCD induced
nonlinear thermal Hall conductivity, $\kappa^{BCD}$, with chemical potential,
$\mu$, for different temperature values. The $\kappa^{BCD}$ values are
represented in the units of $10^{-2}k_{B}^{2}/\hbar$ Å.

We note that a stronger BCD response can be attained by approaching the
critical point with more precision. Finally, we compare the BCD response in
the two distinct topological phases. All the strained buckled honeycomb
lattices exhibit universal behaviour around the topological critical point.
Similar to the observation during reversal of applied electric field, the sign
of BCD changes when going from topologically non-trivial to trivial state.
This is presented in Fig. 5(b). This can be understood to be the result of the
exchange of Berry curvature between VB and CB around the topological phase
transition. Overall, our results put forward a new platform to explore large
and tunable BCD in buckled honeycomb lattices. This electrically switchable
BCD will facilitate the exploration of various exotic quantum mechanical
phenomena, such as NHE Sodemann and Fu (2015), chiral polaritonic effects
Basov  _et al._ (2016), nonlinear Nernst effect Yu  _et al._ (2019); Zeng  _et
al._ (2019), nonlinear thermal Hall effect Zeng  _et al._ (2020), and orbital-
Edelstein Yoda  _et al._ (2018) effects.

##  IV Conclusions

In summary, we introduced a class of elemental systems that exhibit
electrically switchable giant BCD at the Fermi level. In particular, the
elemental buckled honeycomb lattices – silicene, germanene, and stanene –
exhibit an electric field-driven topological phase transition. The transverse
electric field breaks the inversion symmetry of the systems and opens the
possibility of obtaining a large BCD. However, the non-zero value of BCD is
still restricted by the point group symmetry of the crystals. Therefore, we
proposed that the symmetry of the crystals is reduced down to a single mirror
plane using appropriate strain. The strain essentially perturbs the
distribution of Berry curvature and induces asymmetry in a valley.
Consequently, a sizable BCD is obtained for all the systems mentioned above.
Moreover, a vanishing band gap near the critical band gap closing point
triggers a giant BCD at the Fermi level. The value of BCD switches when the
electric field strength crosses a critical value. This flipping can be
explained in terms of the change in sign of Berry curvature across the
critical point. Our study paves the way for exploring field tunable electrical
and thermal nonlinear effects in a class of elemental systems.

##  V Acknowledgments

A.B. sincerely acknowledges Indian Institute of Science for providing
financial support. A.B. also wants to thank S. Roy, A. Bose, and A. Banerjee
for illuminating discussions. N.B.J. acknowledges the support from the Prime
Minister’s Research Fellowship. A.N. acknowledges support from the start-up
grant (SG/MHRD-19-0001) of the Indian Institute of Science and DST-SERB
(project number SRG/2020/000153).

## References

  * Klitzing  _et al._ (1980) K. v. Klitzing, G. Dorda, and M. Pepper, Physical review letters 45, 494 (1980).
  * Du  _et al._ (2021) Z. Du, H.-Z. Lu, and X. Xie, Nature Reviews Physics , 1 (2021).
  * Ortix (2021) C. Ortix, arXiv preprint arXiv:2104.06690 (2021).
  * Ma  _et al._ (2019) Q. Ma, S.-Y. Xu, H. Shen, D. MacNeill, V. Fatemi, T.-R. Chang, A. M. M. Valdivia, S. Wu, Z. Du, C.-H. Hsu, _et al._ , Nature 565, 337 (2019).
  * Kang  _et al._ (2019) K. Kang, T. Li, E. Sohn, J. Shan, and K. F. Mak, Nature materials 18, 324 (2019).
  * Huang  _et al._ (2020) M. Huang, Z. Wu, J. Hu, X. Cai, E. Li, L. An, X. Feng, Z. Ye, N. Lin, K. T. Law, _et al._ , arXiv preprint arXiv:2006.05615 (2020).
  * Sodemann and Fu (2015) I. Sodemann and L. Fu, Physical review letters 115, 216806 (2015).
  * You  _et al._ (2018) J.-S. You, S. Fang, S.-Y. Xu, E. Kaxiras, and T. Low, Physical Review B 98, 121109 (2018).
  * Zhou  _et al._ (2020) B. T. Zhou, C.-P. Zhang, and K. T. Law, Physical Review Applied 13, 024053 (2020).
  * Son  _et al._ (2019) J. Son, K.-H. Kim, Y. Ahn, H.-W. Lee, and J. Lee, Physical review letters 123, 036806 (2019).
  * Joseph and Narayan (2021) N. B. Joseph and A. Narayan, Journal of Physics: Condensed Matter 33, 465001 (2021).
  * Zhang  _et al._ (2018a) Y. Zhang, J. van den Brink, C. Felser, and B. Yan, 2D Materials 5, 044001 (2018a).
  * Facio  _et al._ (2018) J. I. Facio, D. Efremov, K. Koepernik, J.-S. You, I. Sodemann, and J. Van Den Brink, Physical review letters 121, 246403 (2018).
  * Liu and Dai (2020) J. Liu and X. Dai, npj Computational Materials 6, 1 (2020).
  * Yu  _et al._ (2019) X.-Q. Yu, Z.-G. Zhu, J.-S. You, T. Low, and G. Su, Physical Review B 99, 201410 (2019).
  * Zeng  _et al._ (2019) C. Zeng, S. Nandy, A. Taraphder, and S. Tewari, Physical Review B 100, 245102 (2019).
  * Zeng  _et al._ (2020) C. Zeng, S. Nandy, and S. Tewari, Physical Review Research 2, 032066 (2020).
  * Nandy and Sodemann (2019) S. Nandy and I. Sodemann, Physical Review B 100, 195117 (2019).
  * Du  _et al._ (2019) Z. Du, C. Wang, S. Li, H.-Z. Lu, and X. Xie, Nature communications 10, 1 (2019).
  * Kumar  _et al._ (2021) D. Kumar, C.-H. Hsu, R. Sharma, T.-R. Chang, P. Yu, J. Wang, G. Eda, G. Liang, and H. Yang, Nature Nanotechnology 16, 421 (2021).
  * Zeng  _et al._ (2021) C. Zeng, S. Nandy, and S. Tewari, Physical Review B 103, 245119 (2021).
  * Matsyshyn and Sodemann (2019) O. Matsyshyn and I. Sodemann, Physical review letters 123, 246602 (2019).
  * Singh  _et al._ (2020) S. Singh, J. Kim, K. M. Rabe, and D. Vanderbilt, Physical review letters 125, 046402 (2020).
  * Roy and Narayan (2021) S. Roy and A. Narayan, arXiv preprint arXiv:2110.03166 (2021).
  * Battilomo  _et al._ (2019) R. Battilomo, N. Scopigno, and C. Ortix, Physical review letters 123, 196403 (2019).
  * Samal  _et al._ (2021) S. S. Samal, S. Nandy, and K. Saha, Physical Review B 103, L201202 (2021).
  * Xu  _et al._ (2018) S.-Y. Xu, Q. Ma, H. Shen, V. Fatemi, S. Wu, T.-R. Chang, G. Chang, A. M. M. Valdivia, C.-K. Chan, Q. D. Gibson, _et al._ , Nature Physics 14, 900 (2018).
  * Zhang and Fu (2021) Y. Zhang and L. Fu, Proceedings of the National Academy of Sciences 118 (2021).
  * Kane and Mele (2005) C. L. Kane and E. J. Mele, Physical review letters 95, 226801 (2005).
  * Balendhran  _et al._ (2015) S. Balendhran, S. Walia, H. Nili, S. Sriram, and M. Bhaskaran, small 11, 640 (2015).
  * Vogt  _et al._ (2012) P. Vogt, P. De Padova, C. Quaresima, J. Avila, E. Frantzeskakis, M. C. Asensio, A. Resta, B. Ealet, and G. Le Lay, Physical review letters 108, 155501 (2012).
  * Derivaz  _et al._ (2015) M. Derivaz, D. Dentel, R. Stephan, M.-C. Hanf, A. Mehdaoui, P. Sonnet, and C. Pirri, Nano letters 15, 2510 (2015).
  * Zhu  _et al._ (2015) F.-f. Zhu, W.-j. Chen, Y. Xu, C.-l. Gao, D.-d. Guan, C.-h. Liu, D. Qian, S.-C. Zhang, and J.-f. Jia, Nature materials 14, 1020 (2015).
  * Mak  _et al._ (2009) K. F. Mak, C. H. Lui, J. Shan, and T. F. Heinz, Physical review letters 102, 256405 (2009).
  * Ezawa (2015) M. Ezawa, Journal of the Physical Society of Japan 84, 121003 (2015).
  * Deng  _et al._ (2018) J. Deng, B. Xia, X. Ma, H. Chen, H. Shan, X. Zhai, B. Li, A. Zhao, Y. Xu, W. Duan, _et al._ , Nature materials 17, 1081 (2018).
  * Qin  _et al._ (2021) M.-S. Qin, P.-F. Zhu, X.-G. Ye, W.-Z. Xu, Z.-H. Song, J. Liang, K. Liu, and Z.-M. Liao, Chinese Physics Letters 38, 017301 (2021).
  * Xiao  _et al._ (2010) D. Xiao, M.-C. Chang, and Q. Niu, Reviews of modern physics 82, 1959 (2010).
  * Tsirkin (2021) S. S. Tsirkin, npj Computational Materials 7, 1 (2021).
  * (40) “Pythtb: see link,” https://wannier-berri.org/.
  * Giannozzi  _et al._ (2017) P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. D. Jr, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Giannozzi  _et al._ (2009) P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (19pp) (2009).
  * Vanderbilt (1990) D. Vanderbilt, Phys. Rev. B 41, 7892 (1990).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Mostofi  _et al._ (2014) A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Computer Physics Communications 185, 2309 (2014).
  * Wu  _et al._ (2018) Q. Wu, S. Zhang, H.-F. Song, M. Troyer, and A. A. Soluyanov, Computer Physics Communications 224, 405 (2018).
  * Liu  _et al._ (2011) C.-C. Liu, H. Jiang, and Y. Yao, Physical Review B 84, 195430 (2011).
  * (48) A. Bandyopadhyay, N. B. Joseph, and A. Narayan, “Supplementary information: Electrically switchable large berry curvature dipole in silicene, germanene and stanene,” .
  * Yu  _et al._ (2011) R. Yu, X. L. Qi, A. Bernevig, Z. Fang, and X. Dai, Physical Review B 84, 075119 (2011).
  * Soluyanov and Vanderbilt (2011) A. A. Soluyanov and D. Vanderbilt, Physical Review B 83, 235401 (2011).
  * Zhang  _et al._ (2018b) Y. Zhang, Y. Sun, and B. Yan, Physical Review B 97, 041101 (2018b).
  * Basov  _et al._ (2016) D. Basov, M. Fogler, and F. G. De Abajo, Science 354 (2016).
  * Yoda  _et al._ (2018) T. Yoda, T. Yokoyama, and S. Murakami, Nano letters 18, 916 (2018).


# Hierarchical Higher-Order Dynamic Mode Decomposition for Clustering and
Feature Selection

Adrián Corrochano<EMAIL_ADDRESS>School of Aerospace Engineering,
Universidad Politécnica de Madrid, 28040 Madrid, Spain  Giuseppe D’Alessio
Université Libre de Bruxelles, École polytechnique de Bruxelles, Aero-Thermo-
Mechanics Laboratory, Bruxelles, Belgium  CRECK Modeling Lab, Department of
Chemistry, Materials and Chemical Engineering, Politecnico di Milano, Piazza
Leonardo da Vinci 32, 20133 Milano, Italy  Alessandro Parente  Université
Libre de Bruxelles, École polytechnique de Bruxelles, Aero-Thermo-Mechanics
Laboratory, Bruxelles, Belgium  Soledad Le Clainche  School of Aerospace
Engineering, Universidad Politécnica de Madrid, 28040 Madrid, Spain

###### Abstract

In this work, a new algorithm based on the application of higher-order dynamic
mode decomposition (HODMD) is proposed for feature selection and variables
clustering in reacting flow simulations. The hierarchical HODMD (h-HODMD)
performs a reduction of the model order, followed by the iterative selection
of the best reconstructed variables thus creating clusters of features which
can eventually be associated with distinct dynamical phenomena. Firstly,
h-HODMD is combined with different data pre-processing techniques to assess
their influence on the algorithm in terms of reconstruction error. Afterwards,
the proposed algorithm is applied to analyze three different databases
obtained from numerical simulations of a non-premixed co-flow methane flame,
and its performance are compared with the standard HODMD in terms of the
achievable degree of reduction as well as in terms of reconstruction error.
Two detailed kinetic mechanisms are used for the simulations, and different
time horizons are considered to assess the influence of the number of
considered cycles on the accuracy of the reduced model. Results show that
h-HODMD improves the reconstruction for all the variables when compared to the
standard HODMD algorithm. This condition is achieved thanks to the iterative
variables’ clustering: finding dedicated modes for a specific group of
features does in fact lead to a better reconstruction of the dynamics with
respect to the case when the same (global) modes are used to reconstruct the
entire set of variables. Finally, the clusters of variables found by means of
h-HODMD are analyzed, and it is observed that the algorithm can group, in an
unsupervised fashion, chemical species whose behavior is also consistent from
a kinetic point of view. In fact, radicals responsible for fast chemistry
processes (such as propagation and branching reactions) are grouped together,
and the same happens for nitrogen oxides (whose chemical timescales are
instead order of magnitudes slower). The possibility to identify and to
describe correctly the main dynamics with reduced order models for these two
classes of chemical species is extremely important for practical applications
in combustion. In fact, it allows for the possibility to formulate inexpensive
reduced dynamical models for predicting flames liftoff, as well as for
identifying the formation of local extinction and blowout conditions, to
formulate accurate reduced models to describe the formation of pollutants in
aviation, and for control purposes.

Reacting flows; higher-order dynamic mode decomposition; reduced order models

##  I Introduction

Understanding the complex physical and chemical phenomena occurring in
reacting systems is crucial because of the urgent need for a reduction of the
environmental impact of the next generation combustion systems. However,
performing detailed experimental and/or numerical investigations of combustion
systems (especially when large-scale systems are taken into account) is not
always straightforward because of the associated high computational cost [1].

Thus, the adoption of reduced order models (ROMs) represents a viable option
in this context, as the latter are capable of describing complex reacting
systems with low errors and minimal computational costs [2, 3, 4].

One of the most employed algorithms to construct ROMs for reacting flows is
represented by the principal component analysis (PCA) [5], also known as
proper orthogonal decomposition (POD) [6]. The latter is a dimensionality
reduction algorithm based on the projection of the original system onto a
lower-dimensional manifold whose directions (the principal components, PCs)
are orthogonal and obtained as linear combination of the original variables.
In literature, it is possible to find many applications exploiting PCA to
speed-up computational fluid dynamics (CFD) simulations of reacting systems by
reducing the number of the differential equations to be solved for combustion
or atmospheric re-entry applications [7, 8, 9, 10]. In addition, the algorithm
and its local formulation were recently validated for clustering in chemical
kinetics applications [11, 12], for data analysis of Direct Numerical
Simulation (DNS) of reacting jets [13, 14], and for feature extraction and
selection to enhance the classification and regression accuracy of Artificial
Neural Networks (ANN) [15, 16].

Recently, dynamic mode decomposition (DMD) [17] was proposed for the analysis
and the development of ROMs of unsteady systems. The main advantage for this
algorithm over alternative modal decomposition techniques, such as PCA, is
represented by the possibility to retrieve the frequency associated to the
modes to analyze the system. In addition, DMD finds patterns that can be
related to flow instabilities [18, 19, 20, 21], and it can also be used for
dimensionality reduction, as fewer DMD modes are usually capable of
reproducing the main flow dynamics with a reduced computational cost.
Nevertheless, its applicability is usually limited to specific conditions
[22]. Thus, higher order dynamic mode decomposition (HODMD) [23] was found to
be a more robust alternative to standard DMD for non-reactive complex flows
and noisy experimental data, and it has been widely used in the analysis of
several complex flows in industrial applications, either to create ROMs or to
study flow patterns [24, 25, 26, 27]. Recently, HODMD was also proposed and
validated for the development of reduced order models for reacting flows [28].

In the current work, a modified formulation for HODMD is proposed to
accomplish feature selection and clustering of variables when high-dimensional
reacting systems are considered. More specifically, the hierarchical HODMD
(hereinafter referred to as h-HODMD), iteratively identifies the best
reconstructed variables and associates them with the corresponding DMD modes
and frequencies. The performance of the proposed algorithm are assessed using
three databases obtained from CFD simulations of a laminar co-flow flame
taking into account different kinetic mechanisms, as well as different
boundary conditions. The two kinetic mechanisms used to describe methane
oxidation consist of 84 and 53 chemical species, respectively. A wide range of
chemical time scales is also considered since the nitrogen oxides’ chemistry
is also included in one of the two mechanisms. The considered flame, despite
being laminar, exhibits complex physics because of local extinction phenomena.
Finally, the results are compared with the ones obtained by using a standard
HODMD in terms of the achievable degree of reduction as well as in terms of
reconstruction error.

The current work is organized as follows. In Section II, the original HODMD is
explained and the modifications to obtain the h-HODMD are discussed. In
Section III, the numerical setup for the CFD simulations that were considered
for the training datasets is shown and discussed. Finally, the main results
and conclusions are presented in Sections IV and V, respectively.

##  II Higher Order Dynamic Mode Decomposition

Higher Order Dynamic Mode Decomposition (HODMD) [23] is a data driven method
used in fluid dynamics to identify flow patterns. The latter represents an
extension of Dynamic Mode Decomposition (DMD) [17] used for the analysis of
complex flow [29], noisy experiments [30] or turbulent flows [31]. HODMD
decomposes data equispaced in time as an expansion of Fourier-like DMD modes
as

$\boldsymbol{v}(x,y,t_{k})\simeq\sum_{m=1}^{M}a_{m}\boldsymbol{u}_{m}(x,y)e^{(\delta_{m}+i\omega_{m})t_{k}},\quad
k=1,\ldots,K,$ (1)

where $\boldsymbol{u}_{m}$ are the DMD modes and $a_{m}$, $\omega_{m}$,
$\delta_{m}$ are their associated amplitudes, frequencies and growth rates.
The number of DMD modes ($M$) is referred to as spectral complexity, whereas
$K$ in Eq. 1 represents the number of snapshots available for the analyzed
data, which is referred to as the temporal dimension. A thorough description
of the algorithm can be found in Ref. [23], while the associated code is
available in Ref. [32].

###  II.1 Main algorithm

Firstly, equidistant snapshots in time for the considered dynamical system
(with time interval $\Delta t$) are organized as the following snapshot matrix

$\boldsymbol{V}_{1}^{K}=[\boldsymbol{v}_{1},\boldsymbol{v}_{2},\ldots,\boldsymbol{v}_{k},\boldsymbol{v}_{k+1},\ldots,\boldsymbol{v}_{K-1},\boldsymbol{v}_{K}].$
(2)

From the definition in Eq. 2, it follows that the dimensions of the input
matrix are thus $J\times K$, where $K$ is the number of snapshots analyzed and
$J=N_{x}\times N_{y}$ is the number of grid points considered for the CFD
simulations. The HODMD algorithm can be summarized into two main steps: first,
a spatial reduction via singular value decomposition (SVD) is accomplished,
then the DMD-d algorithm is applied.

Within the first step, SVD is applied to the snapshot matrix to reduce noise
and spatial redundancies to only retain the large scales:

$\boldsymbol{V}_{1}^{K}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{T}^{T}=\boldsymbol{U}\hat{\boldsymbol{V}}_{1}^{K},$
(3)

where $\hat{\boldsymbol{V}}_{1}^{K}$ is the reduced snapshot matrix whose
dimensions are $N\times K$. The number of SVD modes to be retained ($N$) can
be determined setting a tolerance $\varepsilon_{1}$ as

$\dfrac{\sigma_{N+1}}{\sigma_{1}}\leq\varepsilon_{1}.$ (4)

By means of the operations discussed above, the spatial complexity can be
reduced from $J$ grid points to $N$ linearly independent vectors. When high-
dimensional data are analyzed, such as the ones obtained from complex reacting
systems, SVD is replaced with Higher Order Singular Value Decomposition
(HOSVD) [33]. To perform this decomposition, the snapshot matrix
$\boldsymbol{V}_{1}^{K}$ is organized in tensor form $\boldsymbol{V}_{ijlk}$
having dimensions $N_{v}\times N_{x}\times N_{y}\times K$. HOSVD can then be
used to reduce the dimensionality in an efficient way. The algorithm applies
standard SVD to the four matrices whose columns are formed by each one of the
four dimensions of the tensor. In particular, if a four-dimensional tensor is
considered, HOSVD is applied as follows

$\boldsymbol{V}_{ijlk}\simeq\sum_{p_{1}=1}^{N_{v}}\sum_{p_{2}=1}^{N_{x}}\sum_{p_{3}=1}^{N_{y}}\sum_{n=1}^{K}\boldsymbol{S}_{p_{1}p_{2}p_{3}n}\boldsymbol{W}^{(v)}_{ip_{1}}\boldsymbol{W}^{(x)}_{jp_{2}}\boldsymbol{W}^{(y)}_{lp_{3}}\boldsymbol{T}_{kn},$
(5)

where $S_{p_{1}p_{2}p_{3}n}$ is the core tensor and $\boldsymbol{W}^{(v)}$,
$\boldsymbol{W}^{(x)}$, $\boldsymbol{W}^{(y)}$ and $\boldsymbol{T}$ are the
modes of the decomposition in each dimension. As well as in SVD, it is
possible to retain a number of modes en each direction, in order to remove the
noise and retain the large scales. The temporal modes matrix $\boldsymbol{T}$
is used as the reduced snapshot matrix.

In the second step, i.e., the application of the DMD-d algorithm, the high
order Koopman assumption is applied to the reduced snapshot matrix as follows

$\hat{\boldsymbol{V}}_{d+1}^{K}\simeq\hat{\boldsymbol{R}}_{1}\hat{\boldsymbol{V}}_{1}^{K-d}+\hat{\boldsymbol{R}}_{2}\hat{\boldsymbol{V}}_{2}^{K-d+1}+\ldots+\hat{\boldsymbol{R}}_{d}\hat{\boldsymbol{V}}_{d}^{K-1}.$
(6)

The Koopman operators $\hat{\boldsymbol{R}}_{1}$, $\hat{\boldsymbol{R}}_{2}$,
$\cdots$, $\hat{\boldsymbol{R}}_{d}$ are encompassed into a single matrix,
which contains the dynamics of the system. When $d=1$ is considered, the
algorithm provides the same results to the standard DMD.

The eigenvalues provide the frequencies $\omega_{m}$ and the growth rates
$\delta_{m}$ of the eigenvectors, which are used to construct the DMD modes
$\boldsymbol{u}_{m}$. Finally, the amplitude of the DMD modes $a_{m}$ are
obtained via least squares fitting of the expansion which was previously
expressed in Eq. 1. Using a second tunable tolerance $\varepsilon_{2}$, it is
possible to then select $M$ DMD modes that are representative for the main
dynamics of the system as follows

$a_{M+1}/a_{1}\leq\varepsilon_{2}.$ (7)

When the considered database is related to the transient state of a numerical
simulation, it is necessary to perform another step to identify the leading
modes driving the flow dynamics. This step consists in computing the
contribution to the transient state $I_{m}$ defined as

$I_{m}=\sum_{k=1}^{K}\lvert
a_{m}e^{(\delta_{m}+i\omega_{m})t_{k}}\rvert\;\|u_{m}\|_{F}^{2}\times\Delta
t,$ (8)

where $\|\bullet\|_{F}^{2}$ is the square of the Frobenius norm. This index
$I_{m}$ takes into account also the damping or growth of the modes in time,
which is reflected in $\delta_{m}$. A much more detailed description about
this criterion can be found in [34].

Finally, the accuracy of the algorithm is measured in terms of the relative
root mean square error (RRMSE) as

$\text{RRMSE}=\sqrt{\frac{\|\boldsymbol{V}^{\text{DMD}}-\boldsymbol{V}\|_{F}}{\|\boldsymbol{V}\|_{F}}}.$
(9)

###  II.2 Pre-processing techniques for applying HODMD to combustion databases

When HODMD is applied to combustion databases, a pre-processing step of the
data in terms of centering and scaling must be applied because of their
multivariate nature [35, 28]. Centering is accomplished by subtracting the
temporal mean of the considered variable to each statistical observation of
the dataset, so that each sample can be seen as a fluctuation from the mean.
Scaling the variables, instead, is motivated by the need of having them lying
in the same order of magnitude. The centered and scaled $i$-th statistical
observation for the $j$-th variable is finally obtained as

$\tilde{x}_{i,j}=\dfrac{x_{i,j}-\bar{x}_{i,j}}{c_{j}},$ (10)

where $\bar{x}_{i,j}$ represents the temporal mean of the $i$-th point
considering the variable’s vector $x_{j}$ and $c_{j}$ is the associated
scaling factor. Four scaling techniques were considered for the current
analysis:

  1. 1.

Range scaling: The difference between the maximum and minimum value of each
variable is used as scaling factor.

  2. 2.

Auto scaling: the standard deviation of the $j$-th variable ($\sigma_{j}$) is
used as scaling factor.

  3. 3.

Pareto scaling: the square root of the standard deviation of the $j$-th
variable ($\sqrt{\sigma_{j}}$) is used as scaling factor.

  4. 4.

Vast scaling: the product of the standard deviation $\sigma_{j}$ and the
coefficient of variation $\sigma_{j}/\bar{x}_{j}$ is used as scaling factor.

The effects of the different scaling factors on modal decomposition techniques
such as PCA and HODMD have been throughly analyzed for combustion applications
in literature [35, 13, 28].

###  II.3 Hierarchical Higher Order Dynamic Mode Decomposition

The hierarchical higher order dynamic mode decomposition (h-HODMD) is a new
algorithm which aims at iteratively identifying the main dynamics of the
considered system, creating clusters of variables having the same dynamic
behavior.

After applying the chosen pre-processing technique to the considered database,
HODMD is applied to the scaled tensor. Once the model’s order is reduced, each
variable can be reconstructed back to the original space by multiplying by the
scaling factor and adding the mean value, as

$\boldsymbol{x}_{j}^{DMD}=\tilde{\boldsymbol{x}}_{j}^{DMD}*c_{j}+\bar{\boldsymbol{x}}_{j},$
(11)

The accuracy of the reconstruction of each variable is measured afterwards by
using the relative root mean square error. A tolerance $\varepsilon_{v}$ is
selected by the user to retain the variables verifying the criterion

$\sqrt{\frac{\|\boldsymbol{x}_{j}^{\text{DMD}}-\boldsymbol{x}_{j}\|_{F}}{\|\boldsymbol{x}_{j}\|_{F}}}<\varepsilon_{v}.$
(12)

Afterwards, HODMD is applied again to the remaining variables so they can be
grouped according to their reconstruction error. The variables of the
considered system are progressively better reconstructed as HODMD is applied
at each iteration, and new modes and frequencies are obtained. A summary of
the iterative process is shown in Figure 1: a group of variables (cluster) is
retained at each iteration and it is removed from the dataset for the next
iteration of HODMD. In this work, the number of clusters created in each
analysis (i.e., the number of iterations for the algorithm) is set to $4$.

Figure 1: Sketch for the h-HODMD algorithm: the variables of the considered
dynamical system are progressively grouped together according to their
reconstruction error, and clusters are obtained in an unsupervised fashion.

##  III Numerical simulations

The CFD simulations carried out to obtain the databases analyzed by h-HODMD in
the following sections are representative for the oxidation of a nitrogen-
diluted methane jet ($65\%$ of methane and $35\%$ nitrogen on a molar basis)
in air. The oxidizer is injected into the domain with a constant velocity of
35 $cm/s$, while the fuel is injected according to the following sinusoidal
velocity profile:

$v(r,t)=v_{max}\left(1-\frac{r^{2}}{R^{2}}\right)[1+A\ sin(2\pi f\ t)],$ (13)

where $v_{max}$ is the maximum velocity ($v_{max}=70cm/s$), $r$ is the radial
coordinate, $R$ represents the internal radius of the nozzle and $A$ is the
amplitude of the perturbation ($A=0.25$).

Three numerical simulations were carried out changing the frequency of the
perturbation $f$ and/or the kinetic mechanism, as summarized in Table 1. The
two detailed kinetic mechanisms employed for the CFD simulations were the
POLIMI_C1_C3_HT_1412 [36], which takes into account 84 species and 1698
reactions, and the GRI-Mech 3.0 [37], which instead consists of 53 species and
325 reactions. The simulations carried out using the GRI-Mech 3.0 represent a
particularly interesting test case for the proposed algorithm. In fact, the
kinetic mechanism takes into account the chemistry associated to the formation
and the consumption of nitrogen oxides (NOx), and both phenomena are well-
known to take place on completely different timescales with respect to other
chemical species and radicals. Thus, this represents an additional degree of
complexity for the formulation of an accurate reduced-order model for the
dynamical systems being considered in this work.

The CFD simulations were carried out by using LaminarSMOKE, an OpenFOAM-based
operator-splitting solver proposed by Cuoci _et al._ [38]. Additional
information regarding the numerical settings for the simulations can be found
in Refs. [11, 12].

| frequency [Hz] | Number of species | Number of reactions
---|---|---|---
Simulation 1 | 10 | 84 | 1698
Simulation 2 | 10 | 53 | 325
Simulation 3 | 20 | 53 | 325
Table 1: Frequency of the perturbation (in Hertz) and number of species and
reactions being used for the CFD simulations.

For each simulation reported in Table 1, a database is extracted in the form
of a fourth-order tensor having dimensions $n_{v}\times n_{x}\times
n_{y}\times n_{t}$, where

  * •

$n_{v}$ is the number of considered variables (temperature and $n_{v}-1$
chemical species, respectively);

  * •

$n_{x}$ is the number of grid points along the streamwise direction;

  * •

$n_{y}$ is the number of grid points along the spanwise direction;

  * •

$n_{t}$ is the number of snapshots equispaced in time, considering $\Delta
t=10^{-3}$.

For all considered simulations it was imposed $n_{x}=n_{y}=150$, while $n_{v}$
was changing depending on the number of species originally taken into account
by the kinetic mechanism, and $n_{t}$ was also different for each database.
The selected values for $n_{v}$ and $n_{t}$ for the three simulations are
reported in Table 2 with the associated number of completed periods (a period
is defined as: $T=2\pi/f$). The influence of this parameter on the performance
of h-HODMD will also be assessed later in Sec. IV.

| $n_{v}$ | $n_{t}$ | Completed periods
---|---|---|---
Database 1 | 83 | 250 | 1
Database 2 | 53 | 250 | 3
Database 3 | 53 | 250 | 5
Table 2: Number of variables, snapshots, and completed periods for each
database considered for the h-HODMD analysis.

##  IV Results

The algorithm discussed in Sec. II was finally applied to the data obtained
from the simulations summarized in Tab. 1. In particular, the databases were
firstly preprocessed by means of the centering and scaling operations reported
in Sec. II.2; afterwards, the algorithm discussed in Sec. II.1 was applied to
extract different groups of variables according to their reconstruction error.
Finally, the dynamics for each group was analyzed.

To quantitatively evaluate the performance of the algorithm proposed in this
work, three different metrics to compute the reconstruction error were
employed at each iteration:

  1. 1.

The relative root mean square error (RRMSE) just after the application of
HODMD:

$RRMSE_{sc}=\sqrt{\frac{\|\tilde{\boldsymbol{V}}^{\text{DMD}}-\tilde{\boldsymbol{V}}\|_{F}}{\|\tilde{\boldsymbol{V}}\|_{F}}},$
(14)

where $\tilde{\boldsymbol{V}}$ is the centered and scaled snapshot matrix and
$\tilde{\boldsymbol{V}}^{\text{DMD}}$ is the reconstruction of this snapshot
matrix after the application of HODMD.

  2. 2.

The RRMSE with the original tensor:

$RRMSE_{or}=\sqrt{\frac{\|\boldsymbol{V}^{\text{DMD}}-\boldsymbol{V}\|_{F}}{\|\boldsymbol{V}\|_{F}}}.$
(15)

In this metric, each variable of the reconstructed tensor is multiplied by its
scaling factor and the mean value is added, as in eq. 11.

  3. 3.

The average RRMSE value for the variables:

$RRMSE_{mean}=\dfrac{\sum_{j=1}^{n_{v}}RRMSE_{or}(x_{j})}{n_{v}}$ (16)

The algorithm was calibrated with the aim to ensure its robustness as well as
to obtain the lowest reconstruction error, according to the procedure
presented in [32, 28]. The values for $d$ and
$\varepsilon=\varepsilon_{1}=\varepsilon_{2}$ adopted for the current work are
reported in Tab. 3.

| Scaling Method
---|---
| Auto | Range | Pareto | Vast
| $d$ | $\varepsilon$ | $d$ | $\varepsilon$ | $d$ | $\varepsilon$ | $d$ | $\varepsilon$
Database 1 | $6$ | $4\times 10^{-2}$ | $6$ | $6\times 10^{-4}$ | $6$ | $6\times 10^{-4}$ | $6$ | $5\times 10^{-3}$
Database 2 | $15$ | $4\times 10^{-2}$ | $12$ | $6\times 10^{-4}$ | $12$ | $1\times 10^{-3}$ | $12$ | $4\times 10^{-3}$
Database 3 | $12$ | $4\times 10^{-2}$ | $10$ | $6\times 10^{-4}$ | $12$ | $2\times 10^{-3}$ | $12$ | $4\times 10^{-3}$
Table 3: Values for $d$ and $\varepsilon$ obtained in the HODMD calibration
process for the cases considered in the analysis.

For the iterative process in h-HODMD, the tolerance $\varepsilon_{v}$ (eq. 12)
was set equal to $\varepsilon_{v,iter}=0.02\times iter$, where $iter$ is the
number of iterations accomplished by the algorithm. Thus, the tolerance
becomes progressively higher with the number of iterations. The number of
variables in each cluster $n_{iter}$ is not kept constant, as this parameter
is only controlled by the tolerance $\varepsilon_{v}$. Other tolerances, i.e.
$\varepsilon_{v,iter}=0.01\times iter$ and $\varepsilon_{v,iter}=0.05\times
iter$, have been used for the calibration. However, a high tolerance means
that each iteration retains a high number of variables which can have
different dynamics among them. On the other hand, a small tolerance may not
include all the variables with the same dynamics in a single cluster.

The first database, analyzed in Figure 2, only considers $\sim$ 1 cycle, so
the simulation can be considered fully transient. Auto scaling appears to
provide similar main dynamics at each iteration, and this is in agreement with
what is expected from theory (as this scaling method gives the same importance
to all the variables). Range scaling leads instead to different dynamics for
each iteration, while Pareto and Vast scaling initially emphasize temperature
and most stable species, respectively, to then give to the rest of variables
the same importance as the number of iterations increases. Because of that,
after removing the best reconstructed variables, the algorithm give back
similar dynamics between two consecutive iterations.

Figure 2: From left to right: Three main frequencies calculated as in Equation
8, $RRMSE_{sc}$, $RRMSE_{or}$ and $RRMSE_{mean}$. From top to bottom: auto
scaling, range scaling, pareto scaling and vast scaling. The analysis
corresponds to Database 1, i.e., to the simulation obtained by imposing a
sinusoidal perturbation having frequency $f=10$ Hz and amplitude A = 0.25, and
carried out with a detailed kinetic mechanism accounting for 82 species.

The reconstruction error for the scaled ($RRMSE_{sc}$) and original tensor
($RRMSE_{or}$) increases with the number of iterations, as the algorithm
progressively removes the variables which are better reconstructed. The first
iteration has a low $RRMSE_{sc}$ for all the scaling methods except auto
scaling, as it does not give more relevance to main variables and/or
temperature. After the first iteration, the $RRMSE_{sc}$ increases for Pareto
and Vast scaling, while it reaches a plateau for Range scaling. The
qualitative behavior for $RRMSE_{or}$ is similar to the one observed for
$RRMSE_{sc}$, with Range scaling providing the best reconstruction error.
Finally, the mean RRMSE ($RRMSE_{mean}$) is stable for auto scaling, and this
appears to be justified for it gives to all variables the same importance. In
range scaling, instead, the $RRMSE_{mean}$ decreases as the number of
iterations increases. When Pareto scaling is applied, $RRMSE_{mean}$ on the
first iteration is much higher than the rest of the iterations. This happens
because the latter initially emphasizes temperature, thus ignoring the
reconstruction of the rest of the variables. After the first iteration, the
rest of the variables become more important, and $RRMSE_{mean}$ decreases. The
$RRMSE_{mean}$ when using vast scaling remains almost constant with the
iterations. Usually, Vast scaling focuses on extremely stable variables but
this effect is mitigated as the transient behavior of the considered database
for the analysis increases.

Figure 3: From left to right: Three main frequencies calculated as in Equation
8, $RRMSE_{sc}$, $RRMSE_{or}$ and $RRMSE_{mean}$. From top to bottom: auto
scaling, range scaling, pareto scaling and vast scaling. The analysis
corresponds to Database 2, i.e., to the simulation obtained by imposing a
sinusoidal perturbation having frequency $f=10$ Hz and amplitude A = 0.25, and
carried out with a detailed kinetic mechanism accounting for 53 species.

When databases 2 and 3 are analyzed (Figures 3 and 4, respectively), similar
results are obtained. For the main dynamics, Auto scaling provides similar
frequencies for all the iterations, Pareto and Vast scaling have different
dynamics for the first iteration and find same frequencies for the rest; while
Range finds different frequencies for each iteration. With regard to the
different reconstruction error metrics, the behavior of each scaling method
follows the same behavior observed for the first database. The only metric
that changes is the $RRMSE_{mean}$ for Vast scaling, which behaves similarly
to Pareto scaling. This is mainly due to the fact transient behavior for the
numerical simulations 2 and 3 does not dominate as more cycles are taken into
account.

Figure 4: From left to right: Three main frequencies calculated as in Equation
8, $RRMSE_{sc}$, $RRMSE_{or}$ and $RRMSE_{mean}$. From top to bottom: auto
scaling, range scaling, pareto scaling and vast scaling. The analysis
corresponds to Database 3, i.e., to the simulation obtained by imposing a
sinusoidal perturbation having frequency $f=20$ Hz and amplitude A = 0.25, and
carried out with a detailed kinetic mechanism accounting for 53 species.

Summarizing the results for the three databases, Auto scaling provides the
worst reconstruction error, although the clusters have similar dynamics. Range
scaling gives the lowest $RRMSE_{mean}$ on the last iteration, while each
cluster of variables on this scaling method have different main dynamics.
Pareto behaves similarly to Auto scaling from the second iteration onward,
although it allows for a better reconstruction of the first cluster. When Vast
scaling is adopted, the first cluster appears to find different dynamics with
respect to the rest of iterations and a great improvement for $RRMSE_{mean}$
is observed from second iteration onwards.

After assessing the contribution for the different scaling methods, the new
iterative formulation of the algorithm was compared with the standard HODMD.
To accomplish this task, the reconstruction error relative improvement
($\eta$) between the two different algorithms is analyzed, which is defined as
follows

$\eta=\dfrac{RRMSE_{HODMD}-RRMSE_{h-HODMD}}{RRMSE_{HODMD}},$ (17)

Equation 17 indicates that if a given variable is better reconstructed by
h-HODMD than by the standard HODMD, a positive value for the improvement is
obtained, wherease the opposite occurs when the aforementioned quantity is
negative.

Figure 5: Reconstruction error improvement for all the variables of the
chemical system computed as reported in Equation 17 for Database 1 (simulation
obtained by imposing a sinusoidal perturbation having frequency $f=10$ Hz and
amplitude A = 0.25, and carried out with a detailed kinetic mechanism
accounting for 82 species), Database 2 (simulation obtained by imposing a
sinusoidal perturbation having frequency $f=10$ Hz and amplitude A = 0.25, and
carried out with a detailed kinetic mechanism accounting for 53 species),
Database 3 (simulation obtained by imposing a sinusoidal perturbation having
frequency $f=20$ Hz and amplitude A = 0.25, and carried out with a detailed
kinetic mechanism accounting for 53 species), and considering different
scaling criteria (from top to bottom: Auto, Range, Pareto, Vast).

The improvement factor is shown in Figure 5 for all the databases and scaling
methods. By adopting Auto scaling, the worst results are achieved as the
reconstruction error gets worse for most of the species. The use of Range
scaling improves instead the reconstruction error for most of the variables,
and the same happens for Pareto scaling. Thus, these two scaling methods
appear to be a better choice when h-HODMD is adopted. Lastly, Vast scaling
appears to favor the iterative algorithm as the transient behavior decreases
(e.g., as when moving from database 2 to database 3).

###  IV.1 HODMD feature selection

As already discussed in Section II, the new iterative formulation for the
HODMD is capable of achieving a lower reconstruction error by grouping
different variables at each iteration, i.e., by applying a model reduction
considering clusters of variables instead of the full input tensor.

Iteration 1 | Iteration 2 | Iteration 3 | Iteration 5 | Iteration 7
---|---|---|---|---
$CH_{4}$ | $CO$ | $HO_{2}$ | $CH_{3}$ | $H_{2}O_{2}$
$CO_{2}$ | $N_{2}O$ | $O$ |  |
$H_{2}O$ | $NO$ | $OH$ |  |
$N_{2}$ | $NO_{2}$ |  |  |
$O_{2}$ |  |  |  |
$T$ |  |  |  |
Table 4: Clusters of variables obtained using h-HODMD for Database 3 (i.e.,
for the simulation obtained by imposing a sinusoidal perturbation having
frequency $f=20$ Hz and amplitude A = 0.25, and carried out with a detailed
kinetic mechanism accounting for 53 species) and using range scaling for the
preprocessing step.

The groups of variables that are found by the iterative algorithm for the
third database considering Range scaling are listed in Table 4, while all the
other cases (in terms of databases and scaling criteria) are made available to
the reader in the supplemental material. In this work, only the most relevant
species from a physical point of view were taken into account to perform the
variables’ clusters analysis, therefore in Table 4 only temperature, fuel,
oxidizer, main oxidation products, main initiation and propagation radicals,
and main nitrogen oxides are reported. The tolerance $\epsilon_{v}$ was set
equal to $\epsilon_{v}=\textit{tol}\times iter$, where the value for tol was
set to be equal to the RRMSE of the sixth best reconstructed variable, so as
the first cluster should always contain the main variables for the reacting
system (i.e., temperature, fuel, oxidizer, and main oxidation products such as
water and carbon dioxide). If an iteration of the algorithm does not report
any variable is then removed and the algorithm runs the next iteration until
there are no variables left, so in this analysis the number of iterations is
not fixed.

Pareto scaling always fails at identifying groups of variables which are
coherent from a chemical point of view, as seen in Table 4 and in Tables SI to
SVII of supplemental material. In fact, O and OH are placed in different
clusters, and they are grouped instead with nitrogen oxides which have a much
slower dynamics. In addition, it also requires the maximum number of
iterations to identify and to cluster all the selected variables. On the other
hand, the application of Vast, and Range scaling yields to better results not
only from a reconstruction error point of view, but also because of the
possibility to have physics-based clusters consisting of species whose kinetic
behavior is similar. In fact, the aforementioned scaling criteria result in
two different clusters separating the species having slow dynamics (such as
the nitrogen oxides) from fast dynamics (radicals responsible for propagation
and branching reactions such as the oxygen radicals, as well as the ones
involved in initiation reactions such as the methane radical). As shown in
Table 4, the three main nitrogen oxides (NO, NO2, and N2O) are grouped
together when Range scaling is utilized, while another cluster is identified
for oxy- and peroxy radicals (O, OH, and HO2). Also, this scaling only needs
seven iterations to get this result, while nine iterations are needed when
Vast scaling is employed.

Figure 6: Reconstruction error improvement for the main variables of the
chemical system computed as reported in Equation 17 for Database 2 (simulation
obtained by imposing a sinusoidal perturbation having frequency $f=10$ Hz and
amplitude A = 0.25, and carried out with a detailed kinetic mechanism
accounting for 53 species), and Database 3 (simulation obtained by imposing a
sinusoidal perturbation having frequency $f=20$ Hz and amplitude A = 0.25, and
carried out with a detailed kinetic mechanism accounting for 53 species),
considering different scaling criteria (from top to bottom: Auto, Range,
Pareto, Vast).

As can be observed from Figure 6, Range scaling allows for a significant
improvement ($\sim$40%) in the reduced model’s accuracy for the description of
radicals that are fundamental for the reacting layer dynamics and liftoff
height prediction such as OH, O, H2O2, and CH3. Also from a qualitative
comparison, such as the one reported in Figure 7 where the original (left side
of each figure) and the reconstructed (right side of each figure) variables’
profiles are compared, it is clear that the error observed when the proposed
algorithm is adopted is negligible.

Finally, the two most relevant modes for temperature and nitrogen oxides are
reported in Figure 8. In particular, the left side of each figure is
representative for $\omega_{1}\sim 126$, which is the frequency of
perturbation of the simulation ($f_{1}=\omega_{1}/2\pi\simeq 20$). This mode
has, for all the species, a dominant component along the streamwise direction
due to the effects of the sinusoidal velocity perturbation. Also, the modes
for the nitrogen oxides are non-negative in the regions where these chemical
species are formed and have the maximum concentration (i.e., in the regions
where the reaction terms are stronger) while negative values are observed
where the transport term is dominant and no reaction is occurring. On the
right side of each contour in Figure 8, the frequency with highest energy
$\omega_{2}=25$ is shown. This particular frequency was selected as it is one
sub-harmonic of the main frequency ($\omega_{2}=\omega_{1}/5$) thus it is
expected to have larger structures. For all the considered variables, it is
possible to observe that the mode is focused on the large flow structures
upstream and downstream, and it also has a peak in the liftoff and the flame
border regions which are due to the strong gradients (both in the streamwise
and in the spanwise directions) they undergo.

Figure 7: Reconstruction of a representative snapshot of the flame for
different variables using h-HODMD (Range scaling). The right side of each
contour represents the original snapshot, while the left side of each contour
represents the reconstruction with h-HODMD. Figures from left to right:
contours of temperature, contours of NO, contours of NO2, contours of N2O.

Figure 8: Real part of the two most important DMD modes (left: $\omega\sim
126$,right: $\omega\sim 25$) for some of the most relevant species using
h-HODMD with range scaling criteria. The modes have been normalised between
$-1$ and $1$.

##  V Conclusions

The article introduces a new iterative algorithm based on HODMD, the
hierarchical HODMD (h-HODMD), to formulate reduced order models for reacting
flows. The novelty in the proposed algorithm is represented by the possibility
to perform model reduction by selecting groups of variables (clusters) thus
achieving a reduction in the error observed for the reduced model. In
addition, for each cluster it is possible to identify a main dynamics which is
characteristic for the considered group of variables. This technique results
to be particularly effective for the analysis of reacting flows, as a high
number of variables is usually considered in the corresponding CFD
simulations. As the algorithm relies on modes with their frequencies, it can
be adapted to decouple diff dynamics in combustion systems. The h-HODMD
iteratively identifies and isolates the best reconstructed variables by
computing the RRMSE of each variable and by imposing a new tolerance settled
by the user. Consequently, in addition to reduced order modeling, it can also
be used for feature selection and variables’ clustering tasks.

The new algorithm has been successfully applied in all the simulations and
scaling methods analysed, obtaining a cluster at each iteration with an
associated main dynamics. Auto scaling makes all the clusters to have the same
dynamics, while it is not able to improve the reconstruction with standard
HODMD. The first cluster using Pareto scaling has different frequencies with
respect to the rest, while greatly improving the reconstruction error,
compared with HODMD. Vast scaling has similar behavior on main dynamics as
Pareto, while only making good improvements on the second and third database.
Finally, Range scaling makes all the clusters have different dynamics and
makes improvement for most of the variables.

When analysing the clusters obtained with the algorithm, Pareto shows poor
results and it is not able to identify groups of variables which are coherent
from a chemical point of view. On the other hand, both Vast and Range scaling
provide good results for both clustering and reconstructing the solution,
while Range scaling gives better results overall.

This new iterative algorithm for higher order dynamic mode decomposition has
been proved to be promising to study high-dimensional data such as the ones
obtained from reacting flow simulation carried out with detailed kinetic
mechanisms. In particular, the proposed method for feature selection and
clustering of variables can be useful for the development of new reduced order
models based on the physics knowledge. This allows new strategies to face
increasingly large databases, reducing their associated computational cost.

###### Acknowledgements.

A.C. and S.L.C. acknowledge the grant PID2020-114173RB-I00 funded by MCIN/AEI/
10.13039/501100011033. S.L.C. and A.C. acknowledge the support of Comunidad de
Madrid through the call Research Grants for Young Investigators from
Universidad Politécnica de Madrid. A.C. also acknowledges the support of
Universidad Politécnica de Madrid, under the programme ‘Programa Propio’.
G.D.A. has received funding from the Fonds National de la Recherche
Scientifique (FRS-FNRS) through a FRIA fellowship. A.P. acknowledges funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme under grant agreement No 714605.

## References

  * Cant [2002] S. Cant, High-performance computing in computational fluid dynamics: progress and challenges, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 360, 1211 (2002).
  * Parente  _et al._ [2011] A. Parente, J. Sutherland, B. B. Dally, L. Tognotti, and P. Smith, Investigation of the mild combustion regime via principal component analysis, Proceedings of the Combustion Institute 33, 3333 (2011).
  * Bellemans  _et al._ [2018] A. Bellemans, G. Aversano, A. Coussement, and A. Parente, Feature extraction and reduced-order modelling of nitrogen plasma models using principal component analysis, Computers & chemical engineering 115, 504 (2018).
  * Coussement  _et al._ [2013] A. Coussement, O. Gicquel, and A. Parente, MG-local-PCA method for reduced order combustion modeling, Proceedings of the Combustion Institute 34, 1117 (2013).
  * Jolliffe [2011] I. Jolliffe, Principal component analysis, in  _International Encyclopedia of Statistical Science_, edited by M. Lovric (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011) pp. 1094–1096.
  * Berkooz  _et al._ [1993] G. Berkooz, P. Holmes, and J. L. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annual review of fluid mechanics 25, 539 (1993).
  * Isaac  _et al._ [2014] B. J. Isaac, A. Coussement, O. Gicquel, P. J. Smith, and A. Parente, Reduced-order PCA models for chemical reacting flows, Combustion and flame 161, 2785 (2014).
  * Malik  _et al._ [2021] M. R. Malik, P. O. Vega, A. Coussement, and A. Parente, Combustion modeling using principal component analysis: A posteriori validation on sandia flames d, e and f, Proceedings of the Combustion Institute 38, 2635 (2021).
  * Bellemans  _et al._ [2017] A. Bellemans, T. Magin, A. Coussement, and A. Parente, Reduced-order kinetic plasma models using principal component analysis: Model formulation and manifold sensitivity, Physical Review Fluids 2, 073201 (2017).
  * D’Alessio  _et al._ [2022] G. D’Alessio, S. Sundaresan, and M. E. Mueller, Automated and efficient local adaptive regression for principal component-based reduced-order modeling of turbulent reacting flows, Proceedings of the Combustion Institute (2022).
  * D’Alessio  _et al._ [2020a] G. D’Alessio, A. Parente, A. Stagni, and A. Cuoci, Adaptive chemistry via pre-partitioning of composition space and mechanism reduction, Combustion and Flame 211, 68 (2020a).
  * D’Alessio  _et al._ [2020b] G. D’Alessio, A. Cuoci, G. Aversano, M. Bracconi, A. Stagni, and A. Parente, Impact of the partitioning method on multidimensional adaptive-chemistry simulations, Energies 13, 2567 (2020b).
  * D’Alessio  _et al._ [2020c] G. D’Alessio, A. Attili, A. Cuoci, H. Pitsch, and A. Parente, Analysis of turbulent reacting jets via principal component analysis, in  _Data Analysis for Direct Numerical Simulations of Turbulent Combustion_ (Springer, 2020) pp. 233–251.
  * D’Alessio  _et al._ [2020d] G. D’Alessio, A. Attili, A. Cuoci, H. Pitsch, and A. Parente, Unsupervised data analysis of direct numerical simulation of a turbulent flame via local principal component analysis and procustes analysis, in  _International Workshop on Soft Computing Models in Industrial and Environmental Applications_ (Springer, 2020) pp. 460–469.
  * Tipler  _et al._ [2022] S. Tipler, G. D’Alessio, Q. Van Haute, A. Parente, F. Contino, and A. Coussement, Predicting octane numbers relying on principal component analysis and artificial neural network, Computers & Chemical Engineering 161, 107784 (2022).
  * D’Alessio  _et al._ [2021] G. D’Alessio, A. Cuoci, and A. Parente, Feature extraction and artificial neural networks for the on-the-fly classification of high-dimensional thermochemical spaces in adaptive-chemistry simulations, Data-Centric Engineering 2 (2021).
  * Schmid [2010] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of fluid mechanics 656, 5 (2010).
  * Chatterjee  _et al._ [2013] S. Chatterjee, A. Mukhopadhyay, and S. Sen, Stability study of laminar flame using proper orthogonal decompostion and dynamic mode decomposition, in  _n3l-Int’l Summer School and Workshop on Non-Normal and Nonlinear Effects In Aero-and Thermoacoustics_ (2013) p. 13.
  * Quinlan and Zinn [2014] J. M. Quinlan and B. T. Zinn, Transverse combustion instabilities: Modern experimental techniques and analysis, in  _50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference_ (2014) https://arc.aiaa.org/doi/pdf/10.2514/6.2014-3682 .
  * Huang  _et al._ [2016] C. Huang, W. E. Anderson, M. E. Harvazinski, and V. Sankaran, Analysis of self-excited combustion instabilities using decomposition techniques, AIAA Journal 54, 2791 (2016).
  * Motheau  _et al._ [2014] E. Motheau, F. Nicoud, and T. Poinsot, Mixed acoustic-entropy combustion instabilities in gas turbines, Journal of Fluid Mechanics 749, 542 (2014).
  * Le Clainche and Vega [2018] S. Le Clainche and J. M. Vega, Analyzing nonlinear dynamics via data-driven dynamic mode decomposition-like methods, Complexity 2018 (2018).
  * Le Clainche and Vega [2017] S. Le Clainche and J. M. Vega, Higher order dynamic mode decomposition, SIAM Journal on Applied Dynamical Systems 16, 882 (2017).
  * Corrochano  _et al._ [2021] A. Corrochano, D. Xavier, P. Schlatter, R. Vinuesa, and S. Le Clainche, Flow structures on a planar food and drug administration (fda) nozzle at low and intermediate reynolds number, Fluids 6, 4 (2021).
  * Lazpita  _et al._ [2022] E. Lazpita, A. Martínez-Sánchez, A. Corrochano, S. Hoyas, S. Le Clainche, and R. Vinuesa, On the generation and destruction mechanisms of arch vortices in urban fluid flows, Physics of Fluids 34, 051702 (2022), https://doi.org/10.1063/5.0088305 .
  * Le Clainche  _et al._ [2019a] S. Le Clainche, R. Moreno-Ramos, P. Taylor, and J. M. Vega, New robust method to study flight flutter testing, Journal of Aircraft 56, 336 (2019a).
  * Le Clainche  _et al._ [2020] S. Le Clainche, D. Izbassarov, M. Rosti, L. Brandt, and O. Tammisola, Coherent structures in the turbulent channel flow of an elastoviscoplastic fluid, Journal of Fluid Mechanics 888 (2020).
  * Corrochano  _et al._ [2022] A. Corrochano, G. D’Alessio, A. Parente, and S. L. Clainche, Higher order dynamic mode decomposition to model reacting flows, (2022), arXiv:Change [physics.flu-dyn] .
  * Le Clainche [2019] S. Le Clainche, Prediction of the optimal vortex in synthetic jets, Energies 12, 10.3390/en12091635 (2019).
  * Mendez  _et al._ [2021] C. Mendez, S. Le Clainche, R. Moreno, and J. Vega, A new automatic, very efficient method for the analysis of flight flutter testing data, Aerospace Science and Technology (2021).
  * Le Clainche  _et al._ [2019b] S. Le Clainche, Z.-H. Han, and E. Ferrer, An alternative method to study cross-flow instabilities based on high order dynamic mode decomposition, Physics of Fluids 31, 094101 (2019b), https://doi.org/10.1063/1.5110697 .
  * Vega and Le Clainche [2020] J. M. Vega and S. Le Clainche,  _Higher order dynamic mode decomposition and its applications_ (Academic Press, 2020).
  * Tucker [1966] L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika 31, 279 (1966).
  * Kou  _et al._ [2018] J. Kou, S. Le Clainche, and W. Zhang, A reduced-order model for compressible flows with buffeting condition using higher order dynamic mode decomposition with a mode selection criterion, Physics of Fluids 30, 016103 (2018).
  * Parente and Sutherland [2013] A. Parente and J. C. Sutherland, Principal component analysis of turbulent combustion data: Data pre-processing and manifold sensitivity, Combustion and flame 160, 340 (2013).
  * Ranzi  _et al._ [2012] E. Ranzi, A. Frassoldati, R. Grana, A. Cuoci, T. Faravelli, A. P. Kelley, and C. K. Law, Hierarchical and comparative kinetic modeling of laminar flame speeds of hydrocarbon and oxygenated fuels, Progress in Energy and Combustion Science 38, 468–501 (2012).
  * Smith [1999] G. P. Smith, Gri-mech 3.0, http://www. me. berkley. edu/gri_mech/ (1999).
  * Cuoci  _et al._ [2013] A. Cuoci, A. Frassoldati, T. Faravelli, and E. Ranzi, Numerical modeling of laminar flames with detailed kinetics based on the operator-splitting method, Energy & Fuels 27, 7730 (2013).