# Tailoring High-Frequency Magnonics in Monolayer Chromium Trihalides

Raí M. Menezes  Department of Physics, University of Antwerp,
Groenenborgerlaan 171, B-2020 Antwerp, Belgium  NANOlab Center of Excellence,
University of Antwerp, Belgium  Departamento de Física, Universidade Federal
de Pernambuco, Cidade Universitária, 50670-901, Recife-PE, Brazil  Denis
Šabani  Department of Physics, University of Antwerp, Groenenborgerlaan 171,
B-2020 Antwerp, Belgium  NANOlab Center of Excellence, University of Antwerp,
Belgium  Cihan Bacaksiz  Department of Physics, University of Antwerp,
Groenenborgerlaan 171, B-2020 Antwerp, Belgium  NANOlab Center of Excellence,
University of Antwerp, Belgium  Clécio C. de Souza Silva  Departamento de
Física, Universidade Federal de Pernambuco, Cidade Universitária, 50670-901,
Recife-PE, Brazil  Milorad V. Milošević<EMAIL_ADDRESS>Department of Physics, University of Antwerp, Groenenborgerlaan 171, B-2020
Antwerp, Belgium  NANOlab Center of Excellence, University of Antwerp, Belgium

###### Abstract

Monolayer chromium-trihalides, the archetypal two-dimensional (2D) magnetic
materials, are readily suggested as a promising platform for high-frequency
magnonics. Here we detail the spin-wave properties of monolayer CrBr3 and
CrI3, using spin-dynamics simulations parametrized from the first principles.
We reveal that spin-wave dispersion can be tuned in a broad range of
frequencies by strain, paving the way towards flexo-magnonic applications. We
further show that ever-present halide vacancies in these monolayers host
sufficiently strong Dzyaloshinskii–Moriya interaction to scatter spin-waves,
which promotes design of spin-wave guides by defect engineering. Finally we
discuss the spectra of spin-waves propagating across a moiré-periodic
modulation of magnetic parameters in a van der Waals heterobilayer, and show
that the nanoscale moiré periodicities in such samples are ideal for
realization of a magnonic crystal in the terahertz frequency range. Recalling
the additional tunability of magnetic 2D materials by electronic gating, our
results situate these systems among the front-runners for prospective high-
frequency magnonic applications.

###### pacs:

Valid PACS appear here

## 1\. Introduction

Two-dimensional (2D) magnetic materials, such as monolayer chromium trihalides
and manganese dichalcogenides, have recently drawn immense attention of both
theoretical and experimental research, due to their fundamental significance
and promising technological applications. The high tunability of the magnetic
interactions in such materials, e.g., by lattice straining Bacaksiz  _et al._
(2021); Vishkayi  _et al._ (2020), electronic gating Huang  _et al._ (2018);
Jiang  _et al._ (2018), and layer stacking Chen  _et al._ (2019); Sivadas  _et
al._ (2018); Shang  _et al._ (2019), among other techniques, directly leads to
manipulation of magnetic textures, such as domain-walls, spin-waves (SWs) and
magnetic skyrmions, thus opening a range of possibilities for novel device
concepts Gong and Zhang (2019). Arguably, the largest impact of these
manipulations can be realized in the field of magnonics, considering that
magnons have been readily detected in atomically-thin chromium tri-iodide
Cenker  _et al._ (2021), with documented spin-wave modes even in terahertz
frequencies Jin  _et al._ (2018). Further advances on that front would make 2D
magnetic materials the prime candidates for ultra-fast information transport
and processing based on magnonics.

In a chromium trihalide monolayer CrX3 (with X = I, Br or Cl), Cr atoms form a
planar honeycomb structure sandwiched between two atomic planes of the halogen
atoms, as illustrated in Fig. 1. The ferromagnetic super-exchange across the
Cr-X-Cr bonds is anisotropic, which together with the weak single-ion
anisotropy of Cr results in a ferromagnetic order with off-plane easy axis
Lado and Fernández-Rossier (2017); Xu  _et al._ (2018). In addition to the
symmetric exchange, the antisymmetric exchange, also known as the
Dzyaloshinskii–Moriya interaction (DMI), can be awaken in such systems with
any structural or electronic asymmetry seen by Cr atoms, such as present in
e.g. Janus monolayers Xu  _et al._ (2020) or in applied electric field Liu
_et al._ (2018); Behera  _et al._ (2019).

In this work we investigate the SW propagation in monolayer CrBr3 and CrI3 as
exemplary representatives of 2D magnetic materials. We obtain all relevant
magnetic parameters from first-principle calculations on the considered
structures Bacaksiz  _et al._ (2021); Sabani  _et al._ (2020), and perform
spin-dynamics simulations of the SW propagation. We proceed to calculate the
SW dispersion relation in the chromium trihalides under strain, revealing the
tunability of magnonic behavior by strain-engineering. We further analyse the
role of ever-present lattice defects (halide vacancies) in considered
monolayers, usually assumed detrimental for the spin-wave propagation. We show
that such lattice defects induce local DMI, which can indeed strongly affect
the SW dynamics, but can be proven useful rather than detrimental – for
example, a tailored pattern of defects can confine the SWs and serve as a SW
guide. Finally, we report the spectra of SWs propagating across periodic
modulation of the magnetic parameters, induced by the moiré pattern of a van
der Waals heterobilayer Xiao  _et al._ (2021); Tong  _et al._ (2018). We show
that such samples, readily made experimentally, can function as a magnonic
crystal Chumak  _et al._ (2017) for the high-frequency SWs, exhibiting band
gaps of propagating SW frequencies. Coupled to the wide range of manipulations
available in 2D materials, our results clearly indicate a number of promising
directions for terahertz magnonics in magnetic monolayers.

Figure 1: (a) Top view of the CrX3 lattice. $\hat{\textbf{n}}_{1}$ and
$\hat{\textbf{n}}_{2}$ represent the main symmetry axes discussed in Sec. 3.
(b) Side view of the CrX3 monolayer.

## 2\. Atomistic spin simulations

In the simulations we consider a spin system arranged in a honeycomb structure
in order to mimic the magnetic moments of Cr atoms in the chromium trihalide
monolayer. To reduce the magnetic behavior of the entire monolayer to just the
simulation within the Cr plane is made possible by ab initio parametrization
of microscopic magnetic interactions for different considered cases, done
according to Refs. Bacaksiz  _et al._ , 2021; Sabani  _et al._ , 2020, where
we performed ab initio four-state energy mapping (4SM) based on density
functional theory (DFT) calculations. For more details on the method we refer
to Ref. Šabani  _et al._ , 2020.

In the spin simulations we consider the quadratic Heisenberg spin Hamiltonian,
given by

$\mathcal{H}=\frac{1}{2}\sum_{i,j}\textbf{S}_{i}\mathcal{J}_{ij}\textbf{S}_{j}+\sum_{i}\textbf{S}_{i}\mathcal{A}_{i}\textbf{S}_{i}.$
(1)

Here $\mathcal{J}_{ij}$ and $\mathcal{A}_{i}$ are the exchange and single ion
anisotropy (SIA) matrices, respectively, and
$\textbf{S}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ is the spin vector at the
$i^{\rm th}$ site. We consider $S=3/2$ for the Cr3+ ions, with $3$ unpaired
valence electrons and quenched orbital moment ($L=0$), which yields a magnetic
moment of $\sim 3\mu_{B}$ per Cr atom, in agreement with the experimental
observationsMcGuire  _et al._ (2015); Kim  _et al._ (2019); Lado and
Fernández-Rossier (2017); Xu  _et al._ (2018); Šabani  _et al._ (2020).

The sum over $i$ in Eq. (1) runs over all Cr sites, while the sums over $i,j$
run over all nearest-neighbor Cr pairs. The exchange matrix can be further
decomposed into a symmetric exchange $J$ and the antisymmetric DMI vector D.
The Hamiltonian then becomes

$\displaystyle\mathcal{H}=$
$\displaystyle\frac{1}{2}\sum_{i,j}\left[J_{\alpha}S_{i}^{\alpha}S_{j}^{\alpha}+J_{\beta}S_{i}^{\beta}S_{j}^{\beta}+J_{\gamma}S_{i}^{\gamma}S_{j}^{\gamma}\right.$
(2)
$\displaystyle+\textbf{D}_{ij}\cdot(\textbf{S}_{i}\times\textbf{S}_{j})\Big{]}$
$\displaystyle+\sum_{i}\left[\mathcal{A}_{\alpha^{\prime}}(S_{i}^{\alpha^{\prime}})^{2}+\mathcal{A}_{\beta^{\prime}}(S_{i}^{\beta^{\prime}})^{2}+\mathcal{A}_{\gamma^{\prime}}(S_{i}^{\gamma^{\prime}})^{2}\right],$

where {$\alpha\beta\gamma$} and
{$\alpha^{\prime}\beta^{\prime}\gamma^{\prime}$} are the local bases of
eigenvectors that diagonalize $J$ and $\mathcal{A}$, respectively.
$S_{i}^{u}=\frac{\textbf{u}\cdot\textbf{S}_{i}}{\textbf{u}\cdot\textbf{u}}$ is
the projection of the $i^{th}$ spin along the vector u. For simplicity, and
since in this work we are interested in high frequency SWs dominated by the
short-range exchange interactions, we neglect the contributions of dipole-
dipole interactions. Notice that the above Hamiltonian is not limited to
isotropic exchange interactions, and can be applied to different atomistic
structures, e.g. including lattice defects where the bases of eigenvectors can
change for specific bonds. For either considered structure, corresponding
exchange matrices for different Cr pairs are obtained from first principles
calculations. Table 1 lists the magnetic parameters obtained for pristine CrI3
and CrBr3 monolayers. Parameters obtained for other specific cases considered
in this paper will be shown where needed.

To simulate the magnonic behavior, we resort to spin dynamics captured by the
Landau-Lifshitz-Gilbert (LLG) equation

$\frac{\partial\hat{\textbf{S}}_{i}}{\partial
t}=-\frac{\gamma}{(1+\alpha^{2})\mu}\left[\hat{\textbf{S}}_{i}\times\textbf{B}_{i}^{\text{eff}}+\alpha\hat{\textbf{S}}_{i}\times(\hat{\textbf{S}}_{i}\times\textbf{B}_{i}^{\text{eff}})\right],$
(3)

where $\gamma$ is the electron gyromagnetic ratio, $\alpha$ is the damping
parameter, $\mu$ is the magnetic moment per Cr atom, and
$\textbf{B}_{i}^{\text{eff}}=-\partial\mathcal{H}/\partial\hat{\textbf{S}}_{i}$
is the effective field. In this work, the LLG spin dynamics simulations are
primarily based on the simulation package Spirit Müller  _et al._ (2019),
adapted to accommodate our (Kitaev) Hamiltonian [Eq. (2)].

| pair | $J^{xx}$ | $J^{yy}$ | $J^{zz}$ | $J^{xy}$ | $J^{xz}$ | $J^{yz}$
---|---|---|---|---|---|---|---
| ($i$-$j$) | (meV) | (meV) | (meV) | (meV) | (meV) | (meV)
CrI3 |  |  |  |  |  |  |
| (1-2) | -5.10 | -3.72 | -4.63 | 0.00 | 0.00 | 0.84
| (2-3) | -4.07 | -4.76 | -4.63 | -0.60 | 0.72 | -0.42
| (2-5) | -4.07 | -4.76 | -4.63 | 0.60 | -0.72 | -0.42
| $\langle J\rangle$ | -4.41 | -4.41 | -4.63 | 0.00 | 0.00 | 0.00
CrBr3 |  |  |  |  |  |  |
| (1-2) | -3.45 | -3.29 | -3.42 | 0.00 | 0.00 | 0.10
| (2-3) | -3.33 | -3.41 | -3.42 | -0.07 | 0.08 | -0.05
| (2-5) | -3.33 | -3.41 | -3.42 | 0.07 | -0.08 | -0.05
| $\langle J\rangle$ | -3.37 | -3.37 | -3.42 | 0.00 | 0.00 | 0.00
| $\Delta$ | $\mathcal{A}^{zz}$ |  |  |  |  |
| (meV) | (meV) |  |  |  |  |
CrI3 | -0.22 | -0.07 |  |  |  |  |
CrBr3 | -0.04 | -0.01 |  |  |  |  |
Table 1: Magnetic parameters for pristine CrI3 and CrBr3 obtained from first-
principle calculations [see Ref. Bacaksiz  _et al._ , 2021 for details].
$J^{xx}$, $J^{yy}$, and $J^{zz}$ are diagonal elements, and $J^{xy}=J^{yx}$,
$J^{xz}=J^{zx}$, $J^{yz}=J^{zy}$ the off-diagonal elements of the exchange
matrix. $\langle J\rangle$ is the average exchange over the three nearest-
neighbour pairs and $\Delta=\langle J^{xx}\rangle-\langle J^{zz}\rangle$ is
the out-of-plane anisotropy. $\mathcal{A}^{zz}$ is SIA parameter, same for
each Cr site. Pairs (i-j) are indicated in Fig. 1 (a).

## 3\. Strain-tunable spin-wave propagation in CrX3 monolayers

The magnetic interactions in 2D materials are very sensitive to any
deformation in the atomic lattice. In particular, it has been shown that
exchange interactions in CrI3 and CrBr3 are significantly affected by either
tensile or compressive strain Bacaksiz  _et al._ (2021). Therefore, in this
section, we discuss how straining the material can affect the propagation of
SWs in magnetic monolayers. For the simulations, we consider both CrBr3 and
CrI3 monolayers in their pristine form, as well as under uniform biaxial
strain. The SW beams are artificially created by a sinusoidal in-plane
oscillating field $\textbf{B}_{\text{input}}=b_{0}\sin(2\pi
f_{\text{in}}t)\hat{k}$ applied in a narrow rectangular region Gruszecki  _et
al._ (2016); Song  _et al._ (2019), where $f_{\text{in}}$ is the input
frequency, $b_{0}$ the field amplitude, and $\hat{k}$ the SW propagation
direction. In the simulations we consider $b_{0}=0.1$ T and damping parameter
$\alpha=0.001$. The SW frequency $f$ and wavelength $\lambda$ are calculated
by fitting a sine function to the magnetization oscillations [see e.g. Fig.2
(a)] as a function of time and space respectively. Fig. 2 (b,c) shows the
dispersion relation, i.e., the relation between the SW frequency and the
wavenumber $k=2\pi/\lambda$ for CrBr3 and CrI3, respectively, obtained for
materials under different strain. The calculated dispersion curves correspond
to the lower-energy magnon modes. Notice that the SW dispersion can be tuned
throughout a wide range of frequencies by straining the honeycomb structure,
which actually demonstrates the ease of tuning magnonics in such 2D materials.
Solid lines in Fig. 2 (b,c) depict the fit of the quadratic expression
$f(k)=Ak^{2}+f_{0}$ to the numerical data. The values of the fitting
parameters $A$ and $f_{0}$ as a function of strain are shown in Fig. 2 (d).
One can clearly notice the distinctly different response of the parameter
$f_{0}$, which corresponds to the zero-momentum SW mode, to strain in CrBr3
and CrI3. As discussed below, this behaviour is directly related to the out-
of-plane exchange anisotropy, whose magnitude for CrI3 increases for either
tensile or compressive strain, while for CrBr3 it changes linearly with strain
[see, e.g., Fig S1 in the supplemental material SM ]. The corresponding energy
gap at $k=0$, given by $\Delta E=hf_{0}$, where $h$ is the Plank constant,
varies in the range of $1.25$ to $1.74$ meV for CrI3 and from $0.18$ to $0.36$
meV for CrBr3, under the $-5$ to $5$ % range of strain. Those values have the
same order of magnitude, but are significantly smaller than those measured for
the case of few layer CrI3 in Ref. Jin  _et al._ , 2018 ($9.4$ meV) and Ref.
Klein  _et al._ , 2018 ($3$ meV). Despite this quantitative discrepancy, we
expect that the here reported strain engineering of spin-wave excitations,
manifesting differently in different Cr-trihalide monolayers, can be validated
by experimental techniques Chen  _et al._ (2017); Qin  _et al._ (2019); Klein
_et al._ (2018); Jin  _et al._ (2018); Chen  _et al._ (2018).

Figure 2: (a) Temporal evolution of the spin components for the case of a $1$
THz spin wave. The $S_{y}$ spin component is vertically shifted for better
visualization. (b,c) Spin-wave dispersion relation obtained in the simulations
for CrBr3 (b) and CrI3 (c), under different strain. Solid lines show the
numerical fit to the quadratic expression $f(k)=Ak^{2}+f_{0}$. (d) The values
of the fitted parameters $A$ and $f_{0}$ obtained as a function of strain.
Solid and dashed lines are the analytical expressions (obtained from Eq. (5))
for $f_{0}$ and $A$, respectively.  Figure 3:  (a) Schematic view of a
wrinkled 2D magnetic material. The induced curvature is nonuniformly straining
the 2D material, while the magnetic anisotropy remains orthogonal to the
monolayer surface. Here, the red shade shows tensile strain and the blue shade
represents the compressive strain. (b) Illustration of the DMI pattern
expected in the case of nonuniform strain in the CrX3 monolayer Sabani  _et
al._ (2020). (c) Snapshot of the simulated magnetic ground state around a
wrinkle in the 2D material. Dashed lines indicate the maximal vertical slope
in the 2D material. (d-e) Snapshots of simulated SW propagation (d) across and
(e) along a wrinkle in CrI3, with $f=0.5$ THz. (f) Snapshot of the SW
simulation corresponding to (e), but for pristine (flat) monolayer.

The SW dispersion can be obtained analytically in the limit of small SW
amplitude by solving the linearized LLG equation Rezende (2020). For that
purpose, let us simplify the Hamiltonian [Eq. (2)] as follows. First notice
that, due to the symmetry of the lattice, the DMI contribution is annulled,
and the only nonzero term of SIA matrix is $\mathcal{A}^{zz}$, Šabani  _et
al._ (2020) such that the SIA contribution to the energy of the $i^{th}$ spin
becomes $\mathcal{A}^{zz}(S_{i}^{z})^{2}$. Assuming that both CrBr3 and CrI3
have strong out-of-plane anisotropy, the exchange energy can be rewritten in
terms of an isotropic exchange $J_{0}=\langle J^{xx}\rangle$, and the
anisotropic term $\Delta=\langle J^{zz}\rangle-\langle J^{xx}\rangle$, where
$\langle...\rangle$ represents the average over the three nearest-neighbour
pairs. The lattice symmetry guarantees that $\langle J^{yy}\rangle=\langle
J^{xx}\rangle=J_{0}$ and $\langle J^{ab}\rangle=0$ for $a\neq b$ [see, e.g.,
Table 1]. The Hamiltonian for the $i^{th}$ spin takes the simple form of the
XXZ model Zhang  _et al._ (2021); Torelli  _et al._ (2019); Lado and
Fernández-Rossier (2017)

$\mathcal{H}_{i}=J_{0}\sum_{j}\textbf{S}_{i}\cdot\textbf{S}_{j}+\Delta\sum_{j}S_{i}^{z}S_{j}^{z}+\mathcal{A}^{zz}(S_{i}^{z})^{2},$
(4)

where the sum in $j$ runs over the three nearest-neighbours of the $i^{th}$
spin. The SW dispersion is then calculated by assuming the linearized solution
$\hat{S}^{z}\approx 1$, $\hat{S}^{x}\approx
A_{0}e^{i(\textbf{k}\cdot\textbf{r}-\omega t)}$ and $\hat{S}^{y}\approx
iA_{0}e^{i(\textbf{k}\cdot\textbf{r}-\omega t)}$, where $A_{0}\ll 1$
represents the SW amplitude; $\omega$ is the SW angular frequency and r is the
position of the considered spin. Substituting that into the LLG equation [Eq.
(3)], with the effective field derived from Eq. (4), we obtain SM

$\omega(k)=\frac{\gamma
S^{2}}{(1+\alpha^{2})\mu}\left[3\Delta+2\mathcal{A}^{zz}+J_{0}g(\textbf{k})\right],$
(5)

with

$\displaystyle g(\textbf{k})=$ $\displaystyle
3-\cos(ka)-2\cos(ka/2),\quad\text{if
}\textbf{k}\parallel\hat{\textbf{n}}_{1},$ (6) $\displaystyle g(\textbf{k})=$
$\displaystyle 2-2\cos(ka\sqrt{3}/2),\quad\text{if
}\textbf{k}\parallel\hat{\textbf{n}}_{2},$

where $a$ is the distance between nearest-neighbour pairs, and
$\hat{\textbf{n}}_{1}$ and $\hat{\textbf{n}}_{2}$ are the main symmetry axis
of the honeycomb lattice, as shown in Fig. 1.

Notice that in the limit of small $ka$ Eq. (6) can be approximated as
$g(\textbf{k})\approx\frac{3a^{2}}{4}k^{2}$ for both
$\textbf{k}\parallel\hat{\textbf{n}}_{1}$ and
$\textbf{k}\parallel\hat{\textbf{n}}_{2}$. The SW frequency $f=\omega/2\pi$
then assumes quadratic dependence on the wavenumber $k$. The solid and dashed
lines in Fig. 2(d) show the analytical solutions (with $\alpha^{2}\ll 1$) for
the zero-momentum SW frequency $f_{0}=\frac{\gamma
S^{2}}{2\pi\mu}(3\Delta+2\mathcal{A}^{zz})$ and the quadratic coefficient
$A=\frac{\gamma S^{2}}{2\pi\mu}\frac{3a^{2}J_{0}}{4}$, respectively, which are
in very good agreement with the numbers obtained in the simulations. The
values of $\Delta$, $\mathcal{A}^{zz}$ and $J_{0}$ as a function of strain are
shown in Fig. S1 of the supplemental material SM . One thus concludes that the
XXZ model [Eq. (4)] is a suitable approximation to describe properties of SWs
in uniform CrI3 and CrBr3 2D magnets.

In the case where the CrX3 lattice experiences non-zero DMI, such as in Janus
structures Xu  _et al._ (2020) or in the presence of out-of-plane applied
electric field Liu  _et al._ (2018); Behera  _et al._ (2019), the linearized
solution for the SW dispersion [Eq. (5)] results in an extra term that has
linear dependence on the wavevector $k$. Cortés-Ortuño and Landeros (2013) The
contribution of DMI to the SW dispersion depends on the angle between the DMI
vector and the axis around which the spins rotate, and the strongest effect is
observed when the magnetization lies in the same plane as the DMI vectors.
Details on the calculation of SW dispersion in presence of uniform DMI are
provided in the supplemental material SM .

## 4\. Flexo-magnonics

The strong response of SW properties to strain in the 2D material suggests the
possibility of flexo-magnonic applications. In fact, 2D materials are mostly
very flexible and can be strain engineered in a multitude of ways. For
example, localized strain can be induced by growing the 2D material on top of
a patterned substrate Wang  _et al._ (2019); Peng  _et al._ (2020), or by
placing it onto elastometric or piezoelectric substrate, whose compression may
lead to wrinkling and buckling of the 2D material Castellanos-Gomez  _et al._
(2013); Sun and Liu (2019); Peng  _et al._ (2020). In addition, bubbles and
tents can be formed in the 2D material by trapping water, gas or solid
nanoparticles at the interface between the magnetic monolayer and the
substrate Khestanova  _et al._ (2016); Dai  _et al._ (2019).

To exemplify the interaction of SWs with such localized heterostrained
structures we simulate the interaction of a propagating SW with a wrinkle in
the CrI3 monolayer. Fig. 3(a) illustrates the considered system. In such a
configuration, the position of a generic Cr atom in the curved monolayer can
be parametrized as $\textbf{r}^{\prime}=\textbf{r}+\textbf{u}(\textbf{r})$,
where $\textbf{r}=(x,y,0)$ is the atom position in the pristine (flat) lattice
and
$\textbf{u}(\textbf{r})=(u_{x}(\textbf{r}),u_{y}(\textbf{r}),h(\textbf{r}))$
the position displacement. Here, $u_{x}$ and $u_{y}$ represent the in-plane
displacement and $h(\textbf{r})$ is a scalar field accounting for out-of-plane
deformations. The red and blue shades in Fig. 3 (a) illustrate the expected
regions of tensile and compressive strain respectively Quereda  _et al._
(2016). In addition to the local strain, the curvature also induces a rotation
in the direction of the magnetic anisotropy, which points orthogonal to the
surface of the magnetic film [Fig. 3(a)]. Moreover, due to the non-uniform
deformation, the corresponding symmetry-breaking gives rise to localized DMI
Sabani  _et al._ (2020). Figure 3(b) illustrates the expected DMI
configuration in the case of a non-uniform uniaxial strain in the wrinkled
CrX3 monolayer. DFT calculations show that the DMI increases linearly with
local uniaxial strain, with relatively small magnitude of up to $D=0.08$ meV
for $6\%$ of local strain Sabani  _et al._ (2020).

Figure 4:  (a) CrX3 lattice in the vicinity of a single halide vacancy.
Induced DMI vectors are shown as red arrows. (b-d) Snapshots of simulated SW
propagation in presence of randomly distributed defects in CrI3, for defect
densities $\eta=0$, $\eta=0.01$, $\eta=0.02$ vacancies per Cr atom,
respectively. The input SW frequency is $0.5$ THz. (e-g) Corresponding output
SW amplitudes for the cases shown in (b-d), respectively, measured at a
distance of $100$ nm from the SW source. Colors in (e-g) correspond to
different time instances during the $2$ ps SW period.  Figure 5: (a) The
normalized transmission of SWs as a function of the density of defects,
$\eta$, measured at a distance of $y=100$ nm from the SW source, with $f=0.5$
THz. (b) SW transmission as a function of the distance to the output, for
$\eta=0.01$ and $f=0.5$ THz. Here, the transmission amplitudes are normalized
with respect to the case of $\eta=0$, measured at the same distance from the
SW source. Figure 6: (a) CrX3 lattice with a line of halides vacancies.
Induced DMI vectors are illustrated as red arrows. (b) Reflection of a SW by a
defect line shown in (a). Snapshot captures a spin-wave of frequency $0.5$ THz
reaching the defect line in CrI3 under an incident angle of $\phi=\pi/4$.
Arrows indicate the SW propagation direction. (c-d) Transmission and
reflection coefficients as a function of the incident angle $\phi$ for CrBr3
(c) and CrI3 (d), obtained for SW frequency $f=0.5$ THz and damping
$\alpha=0.001$. (e) Frequency dependence of the transmission coefficient $T$
for both CrBr3 and CrI3 at fixed incident angle $\phi=\pi/6$. (f) Snapshot of
the simulations for the SW confinement between two realistically broadened
defect lines in the magnetic 2D material.

First-principles calculations of such a curved structure are very challenging,
requiring relaxation of the atomic structure over a large area before
calculating locally the magnetic parameters for a number of distinctly
different Cr-Cr pairs due to heterostraining. Instead, as a first step in
understanding the interaction of SWs in such systems, here we model the
material curvature by rotating the local bases of eigenvectors of the magnetic
Hamiltonian [Eq. (2)] accordingly to the surface normal, and vary the
magnitude of interactions accordingly to the local strain profile of the
curvature. Note that the considered magnetic Hamiltonian is completely
determined by the exchange and anisotropy matrices and the angle between
neighbouring spins, therefore, local deformations can be modeled by
correspondingly modifying the interaction matrices and precise atomic
positions do not play a role. Fig. 3(c) shows the magnetic ground state
obtained in such simulations for a single wrinkle of Gaussian shape
$h(\textbf{r})=h_{0}e^{-\frac{x^{2}}{w^{2}}}$, where $h_{0}=20$ nm and $w=15$
nm define the height and width of the wrinkle, respectively. The considered
strain profile has maximal deformation of $\epsilon=5\%$ on top of the wrinkle
Yang  _et al._ (2021); Quereda  _et al._ (2016); Deng  _et al._ (2018) and
satisfies $\epsilon(x)\propto\partial^{2}h(\textbf{r})/\partial x^{2}$. Notice
that the curvature induces local canting in the ground-state magnetization and
increases spin-stiffness due to emergent DMI. Fig. 3(d-f) shows snapshots of
the simulated propagation of a $0.5$ THz SW (d) across and (e) along the
wrinkle, where the SW can either be reflected by or confined along the flexed
region, respectively, in comparison with the case of a pristine (flat) sample
(f). Recalling the wide range of possibilities for nanoengineered manipulation
of 2D materials, this basic example suffices to promote ideas towards flexo-
magnonic applications.

## 5\. Scattering of spin-waves on lattice defects

Isolated defect |  |  |  |  |
---|---|---|---|---|---
| pair | $D^{x}$ | $D^{y}$ | $D^{z}$ | $|D|$
| ($i$-$j$) | (meV) | (meV) | (meV) | (meV)
CrI3 |  |  |  |  |
| (1-2) | 0.29 | -0.20 | -0.53 | 0.64
| (2-3) | 4.05 | -2.33 | -0.39 | 4.69
| (2-5) | 0.17 | 0.09 | -0.46 | 0.50
| (3-4) | 0.01 | 0.11 | -0.06 | 0.12
| (3-6) | 0.01 | -0.38 | 0.19 | 0.42
CrBr3 |  |  |  |  |
| (1-2) | -0.05 | -0.13 | 0.26 | 0.29
| (2-3) | 2.24 | -1.31 | -0.08 | 2.60
| (2-5) | 0.00 | 0.06 | -0.04 | 0.07
| (3-4) | 0.02 | 0.03 | 0.07 | 0.08
| (3-6) | -0.09 | -0.14 | 0.07 | 0.18
Line of defects |  |  |  |  |
| pair | $D^{x}$ | $D^{y}$ | $D^{z}$ | $|D|$
| ($i$-$j$) | (meV) | (meV) | (meV) | (meV)
CrI3 |  |  |  |  |
| (1-2) | 0.01 | -0.43 | 0.06 | 0.44
| (2-3) | -0.05 | -0.55 | 2.42 | 2.48
| (2-5) | 3.61 | -0.27 | 4.23 | 5.57
| (3-4) | -0.09 | 0.03 | -0.23 | 0.25
| (3-6) | 3.61 | -0.27 | 4.23 | 5.57
CrBr3 |  |  |  |  |
| (1-2) | 0.00 | -0.09 | 0.01 | 0.09
| (2-3) | 0.12 | 0.14 | 0.47 | 0.50
| (2-5) | 2.09 | -0.02 | 1.50 | 2.57
| (3-4) | -0.04 | -0.03 | 0.01 | 0.05
| (3-6) | 2.09 | -0.02 | 1.50 | 2.57
Table 2: Induced DMI parameters at an isolated defect (halide vacancy) and at
a line of such defects in CrI3 and CrBr3, obtained from first-principle
calculations Sabani  _et al._ (2020). The complete exchange matrix for both
cases is provided in the supplemental material SM . Pairs (i-j) are indicated
in Fig. 4 (a) and Fig. 6 (a).

Vacancies are the most commonly observed defect in 2D materials. They are not
only an unwanted product of synthesis or manipulation of the sample, but can
also be deliberately artificially induced by vacancy engineering Jiang  _et
al._ (2019); Schuler  _et al._ (2019). In this section we therefore
investigate the interaction of SWs with such defects (halides vacancies) in
CrX3 lattice. The local breaking of the inversion symmetry seen by Cr atoms in
the vicinity of the lattice defect results in appearance of DMI, as shown in
Fig.4(a). It is well known that SWs are strongly affected by the DMI
interaction and can even be controllably refracted/reflected at an interface
where DMI changesMulkers  _et al._ (2018); Wang  _et al._ (2018). In the case
of lattice defects, ab initio calculations reveal DMI vectors of magnitudes up
to $D=4.69$ meV for CrI3 and $D=2.60$ meV for CrBr3 in the vicinity of halide
vacancies Sabani  _et al._ (2020). One thus expects pronounced effects of such
strong variations of the local magnetic interactions on the SW propagation in
the 2D magnetic materials.

To start with, we consider the scattering of SWs in CrBr3 and CrI3 with
different density of randomly distributed halide vacancies $\eta$ (density
expressed as the number of vacancies per Cr atom). The used magnetic
parameters, as obtained from first-principles calculations, are listed in
Table 2 and in the supplemental material SM . Fig. 4(b-d) shows the snapshots
of simulated SW propagation through the CrI3 monolayer for $\eta=0$, $0.01$
and $0.02$, respectively. The corresponding output SW amplitudes, measured at
a distance of $100$ nm from the SW source [Fig.4 (e-g)], reveal that the
transmitted wave is strongly affected in the case of high density of defects.
In order to quantify and compare the SW transmission for different defect
densities, we integrate the output SW amplitudes along the direction
perpendicular to the SW propagation and normalize it with respect to the case
without defects, i.e., for $\eta=0$. Fig. 5 (a) shows such normalized
transmission for the $0.5$ THz SWs in CrBr3 and CrI3. Notice that SWs in CrI3
are more affected by the halide vacancies when compared to CrBr3. We relate
this property to the large DMI interactions induced in CrI3, which are
approximately twice as large as in the CrBr3 sample, predominantly due to
stronger spin-orbit coupling on iodine compared to bromine.

Scattering on defects, as shown in Fig. 4, especially over large distances,
must clearly be taken into account when designing novel magnonic devices based
on 2D materials. Therefore in Fig. 5(b) we show the SW transmission as a
function of the output distance, normalized with respect to the case of
$\eta=0$, measured at the same distance from the SW source. Notice that the
normalized transmission decreases quasi-linearly as a function of the output
distance, where the transmission in CrI3 decreases at approximately twice
higher rate than in CrBr3, for $\eta=0.01$, thus again reflecting the stronger
response of SWs to defects in CrI3. Finally we note to not have observed any
significant frequency dependence (for the range of $0.3\text{\,}1.5\text{\,}$
THz) in the results presented in this section.

### 5.1.Defect-engineering of magnonic circuitry

As demonstrated above, the induced strong variations in magnetic parameters
around the lattice defects strongly affect the SW propagation. Therefore, a
designed pattern of defects, such as a chain of adjacent defects, may be able
to e.g. confine the SWs and serve as a waveguide. To illustrate this point, we
investigate the interaction of SWs with a row of consecutive halide vacancies
in CrX3 lattices, as a possible way of controlling the SW propagation
direction in magnetic 2D materials.

The emergent DMI in the considered structure, as illustrated in Fig. 6(a), is
similar in magnitude to DMI found in the vicinity of isolated defects. The
exact values are listed in Table 2 for comparison, while complete exchange
matrices are provided in the supplemental material SM . Fig. 6(b) shows a
snapshot of the simulation for a SW directed across the defect line in CrI3,
where we consider a $0.5$ THz SW reaching the defect line under an incident
angle of $\phi=\pi/4$. Notice that only a small fraction of the SW is able to
cross the defect line while the major part is reflected. In Fig. 6(c-d) we
show the transmission coefficients, defined as Macke and Goll (2010); Wang
_et al._ (2013)

$T=\left(\frac{\Delta_{T}}{\Delta_{0}}\right)^{2},R=\left(\frac{\Delta_{R}}{\Delta_{0}}\right)^{2},$
(7)

where $\Delta_{T}$ and $\Delta_{R}$ are the amplitudes of the transmitted and
reflected waves, respectively, and $\Delta_{0}$ is the corresponding SW
amplitude in the absence of defects, calculated at the same distance from the
SW source. The reflection and transmission coefficients are defined to satisfy
the relation $R=1-T$. In Fig. 6(c-d) the transmission and reflection
coefficients are shown as a function of incident angle $\phi$ for CrBr3 [Fig.
6(c)] and CrI3 [Fig. 6(d)]. Notice that SWs experience stronger reflection by
the defect line in CrI3 while most part of the SW is transmitted across the
vacancy defect line in CrBr3. Once again we relate this property to the larger
DMI induced at defects in CrI3. Similar to an interface where DMI changes
Mulkers  _et al._ (2018), the reflection is enhanced with increasing the
incident angle from $\phi=0$ towards $\phi=\pi/2$. However, total reflection
is not found in our simulations, which is due to a finite cross-talk of spins
on two sides of the defect line, where DMI is not present, which is different
in the case of a DMI interface. Fig. 6(e) shows the frequency dependence of
the transmission coefficient $T$ for both CrBr3 and CrI3 for a fixed incident
angle $\phi=\pi/6$, which shows a weak increase of transmission with
increasing SW frequency.

Of course, such a precise pattern of defects is unlikely to be reproduced
experimentally by vacancy engineering in a 2D material Jiang  _et al._ (2019);
Schuler  _et al._ (2019).

Figure 7: Example of a moiré magnonic crystal. (a-b) Illustration of a
CrBr3/MoSe2 vdW heterobilayer, viewed from the side (a) and the top (b). (c)
moiré pattern created between adjacent Br and Se sublayers. (d) Corresponding
DMI profile emerging in the Cr sublayer from the moiré pattern shown in (c).
Here, $\delta_{1}$ and $\delta_{2}$ denote the two main periodicities of the
system. (e) Simulated SW transmission spectra after passing through the moiré
magnonic crystal (with maximal DMI of 3 meV), in comparison with the reference
waveguide (pristine film), where we consider the incident SW propagating along
the $\delta_{1}$ direction. Black arrows indicate the critical frequencies
where SWs suffer destructive interference.

However, a broadened defective region can also be made useful in a similar
sense, even if induced by e.g. low-energy electron-beam Jiang  _et al._ (2019)
or by “scratching” the material with a scanning tip. To illustrate such a
case, in Fig. 6(f) we show the snapshot of the simulation for the SW
propagating in CrI3, across a distribution of defects, randomly placed along 5
nm wide stripes, roughly mimicking surface scratches on the magnetic 2D
material. Notice that such a pattern of defects can also confine the SWs, and
can be engineered to (de)stimulate self-interference patterns, or control
interactions of two incoming spin-waves. The key ingredient for confining and
controlling the SWs remains the local variation of the magnetic parameters,
which can also be induced by electric gating Liu  _et al._ (2018) or strain-
engineering the 2D material, as demonstrated in Sec. 4, albeit on a
significantly larger scale than what is offered by defect engineering.

## 6\. Moiré van der Waals magnonic crystals

Figure 8:  (a-b) Snapshots of SW simulation for the $0.5$ and $0.65$ THz SWs
propagating across the moiré pattern for the same period of time. The incident
SW propagates parallel to the $\delta_{1}$ direction. (c) Snapshots of the SW
simulation for the $0.4$ THz SW propagating across the moiré pattern, after
being emitted in a radially symmetric fashion from a point source. SW
propagation is clearly prohibited along the $\delta_{1}$ directions, but
remains enabled along the $\delta_{2}$ direction.

Magnonic crystals are artificial materials designed in such a way that the
magnetic properties of the media are characterized by a periodic lateral
variationChumak  _et al._ (2017). The SW spectra in such materials exhibit
features such as band gaps of wavelengths that are not allowed to propagate,
which is useful for magnonic filters as well as information transport and
processing based on magnons. Periodic variations of the magnetic parameters in
the 2D magnetic materials can be induced in multiple ways, for example, by
periodic electric gating de Vries  _et al._ (2021), strain engineering Sun and
Liu (2019); Dai  _et al._ (2019); Peng  _et al._ (2020), by growing the 2D
material on top of a patterned substrate Wang  _et al._ (2019), or by inducing
a checkerboard buckling of the monolayer Mao  _et al._ (2020). However, a
magnonic crystal for THz frequencies requires a modulation period of just a
few nanometers Qin  _et al._ (2019) [see SW dispersion relation in Fig. 2],
which is beyond the ability of the fabrication methods mentioned above.
Instead, we propose here to employ a moiré pattern in a van der Waals (vdW)
heterobilayer as the source of the periodic modulation in magnetic parameters
for the design of magnonic crystals in the high-frequency regime. It is
readily established that a moiré pattern can induce periodic variation in the
magnetic parameters of a 2D material Xiao  _et al._ (2021); Tong  _et al._
(2018). Therefore, let us consider the case where the CrX3 monolayer is
stacked on top of another, non-magnetic material, in order to produce a moiré
pattern. For instance, we consider the vdW heterobilayer CrBr3-MoSe2, as
illustrated in Fig. 7(a,b). The superexchange interaction between Cr atoms,
mediated by the Cr-Br-Cr bonds, will be affected by the Br-Se interaction at
the interface between monolayer constituents of the heterobilayer, which in
turn varies with the local stacking throughout the moiré pattern. In other
words, the periodicity of the moiré structure is directly reflected in a
periodic modulation of the resultant magnetic parameters. Fig. 7(c) shows the
moiré pattern created between the Br sublayer (a trigonal lattice with lattice
constant $a_{\text{Br}}=3.698$ Å) and the Se sublayer (a trigonal lattice with
lattice constant $a_{\text{Se}}=3.289$ Å) Wang  _et al._ (2016). Here the
symmetry axes of the two lattices are considered to be aligned, which
corresponds to an angle of $30^{\circ}$ between CrBr3 and MoSe2 supercells
Zollner  _et al._ (2019). The mismatch in the lattice constants creates the
moiré pattern with two characteristic periodicities of $\delta_{1}\approx
5.13$ nm and $\delta_{2}\approx 2.96$ nm, as indicated in Fig. 7(d). Such
short periodicities found in the moiré heterostructures are, therefore,
capable of interacting with high frequency SWs and can be adjusted, for
example, by fine interlayer twisting, or replacing the MoSe2 layer with a 2D
material with a different lattice constant and/or different lattice symmetry.

Therefore, in the present case the interaction between Br and Se atoms breaks
the inversion symmetry seen by the Cr sublayer where magnetism resides, and
significantly affects the superexchange owing to the strong spin-orbit
coupling in MoSe2 Zhu  _et al._ (2011); Reyes-Retana and Cervantes-Sodi
(2016). As a result, one can expect a periodic modulation in the DMI parameter
throughout the magnetic layer, directly mapped on the moiré pattern between Br
and Se sublayers. As a first step in understanding the magnonics in such a
complex system, and considering that DMI has the main impact on SWs, we will
consider a moiré periodic modulation only in the DMI parameter, and assume the
other magnetic parameters in CrX3 monolayer remain unchanged by vdW
heterostructuring. Fig. 7(d) shows the used DMI profile, reflecting the moiré
pattern, according to the expression $D=D_{0}\exp(-\sqrt{r^{2}+d^{2}}/r_{0})$,
where $r$ is the lateral distance between Br and Se atoms, $d=3.5$ Å is the
separation between the layers Zollner  _et al._ (2019) [see Fig. 7 (a)] and
$D_{0}$ and $r_{0}$ are constants that define the DMI magnitude. Here we
consider $D_{0}=3\exp(d/r_{0})$ meV and $r_{0}=0.2$ Å, such that the maximum
magnitude of DMI is $3$ meV. The small value of $r_{0}$ reflects the short
range interaction between neigbouring orbitals, such that the maximal
interaction is given when Se atoms are on top of Br ones. Notice that a SW
propagating in such a system will experience different periodic modulations,
depending on the propagation direction, with respect to the two main
periodicities in the system $\delta_{1}$ and $\delta_{2}$ [cf. Fig. 7(d)].

Fig. 7(e) shows the simulated SW transmission spectrum of the moiré magnonic
crystal described above, in comparison with the reference waveguide (pristine
CrBr3 monolayer), where we considered the incident SW propagating parallel to
the $\delta_{1}$ direction. Notice that the transmission spectrum of the moiré
system significantly differs from the pristine medium, as three distinct
valleys of prohibited magnon frequencies are visible. Such a spectrum is
typical of magnonic crystalsChumak  _et al._ (2017). Fig. 8(a-b) shows the
snapshots of SW simulation for the $0.5$ and $0.65$ THz SWs propagating across
the moiré DMI pattern, where in the latter case the SWs suffer destructive
interference and do not propagate.

Finally, one notices that the hexagonal symmetry in the modulation of magnetic
parameters results in the anisotropic dispersion of the spin waves. Fig. 8(c)
shows the snapshots of SW simulation for the $0.4$ THz SW emitted radially
across the moiré pattern, from a point source. In this case, SW propagation is
effectively blocked along the $\delta_{1}$ directions, as expected from the
transmission spectra in Fig. 7(e), but remains enabled along the $\delta_{2}$
directions – demonstrating the anisotropic SW dispersion characteristic of
moiré magnonic systems.

## 7\. Conclusions

In summary, we presented an overview of magnonic properties in monolayer
chromium trihalides CrBr3 and CrI3, including selected possibilities for
tuning the magnonic behavior in these and similar magnetic 2D materials. We
revealed that the spin-wave (SW) dispersion relation can be tuned in a broad
range of frequencies by straining the 2D magnetic monolayer, and that the SW
dispersion has different response to strain in CrBr3 compared to CrI3, which
is directly related to differences in their out-of-plane anisotropy. We
further demonstrated the possibility of controlling the SW propagation by
strain-engineering the 2D material, paving the way towards flexo-magnonic
applications. As a next step, we discussed the effect of ever-present
structural defects (in this case halide vacancies) on SW propagation in
chromium trihalide monolayers. Our calculations show that such lattice defects
induce sufficiently large DMI to strongly disturb the SW dynamics, especially
in CrI3. However, we also show that a pattern of defects can be designed
favorably for magnonic circuitry, as extended defect lines are able to confine
SWs and may be useful for SW-guides. Finally, we discuss the spectra of SWs
propagating across periodic modulation of the magnetic parameters, induced by
a moiré pattern in a van der Waals heterobilayer of a magnetic and a non-
magnetic 2D material. Such moiré heterostructures exhibit the necessary
nanoscale modulation period to work as a magnonic crystal for the terahertz
SWs characteristic for CrX3 monolayers. Even though we presented only a small
fraction of possible manipulations of 2D materials, their pronounced impact on
local magnetism and spin excitations, facilitated heterostructuring, as well
as the ease of employing additional electronic gating and nanolithography in
any desired geometry, promotes magnetic 2D materials to the prime echelon of
candidates for terahertz magnonics by design.

## Acknowledgements

This work was supported by the Research Foundation - Flanders (FWO-
Vlaanderen), the Special Research Funds of the University of Antwerp (BOF-UA),
and Brazilian agencies FACEPE, CAPES and CNPq.

## References

  * Bacaksiz  _et al._ (2021) C. Bacaksiz, D. Šabani, R. M. Menezes, and M. V. Milošević, “Distinctive magnetic properties of CrI3 and CrBr3 monolayers caused by spin-orbit coupling,” Physical Review B 103, 125418 (2021).
  * Vishkayi  _et al._ (2020) S. I. Vishkayi, Z. Torbatian, A. Qaiumzadeh, and R. Asgari, “Strain and electric-field control of spin-spin interactions in monolayer CrI3,” Physical Review Materials 4, 094004 (2020).
  * Huang  _et al._ (2018) B. Huang, G. Clark, D. R. Klein, D. MacNeill, E. Navarro-Moratalla, K. L. Seyler, N. Wilson, M. A. McGuire, D. H. Cobden, D. Xiao, _et al._ , “Electrical control of 2D magnetism in bilayer CrI3,” Nature Nanotechnology 13, 544–548 (2018).
  * Jiang  _et al._ (2018) S. Jiang, L. Li, Z. Wang, K. F. Mak, and J. Shan, “Controlling magnetism in 2D CrI3 by electrostatic doping,” Nature nanotechnology 13, 549–553 (2018).
  * Chen  _et al._ (2019) W. Chen, Z. Sun, Z. Wang, L. Gu, X. Xu, S. Wu, and C. Gao, “Direct observation of van der waals stacking–dependent interlayer magnetism,” Science 366, 983–987 (2019).
  * Sivadas  _et al._ (2018) N. Sivadas, S. Okamoto, X. Xu, C. J. Fennie, and D. Xiao, “Stacking-dependent magnetism in bilayer cri3,” Nano Letters 18, 7658–7664 (2018).
  * Shang  _et al._ (2019) J. Shang, X. Tang, X. Tan, A. Du, T. Liao, S. C. Smith, Y. Gu, C. Li, and L. Kou, “Stacking-dependent interlayer magnetic coupling in 2d cri3/crgete3 nanostructures for spintronics,” ACS Applied Nano Materials 3, 1282–1288 (2019).
  * Gong and Zhang (2019) C. Gong and X. Zhang, “Two-dimensional magnetic crystals and emergent heterostructure devices,” Science 363 (2019).
  * Cenker  _et al._ (2021) J. Cenker, B. Huang, N. Suri, P. Thijssen, A. Miller, T. Song, T. Taniguchi, K. Watanabe, M. A. McGuire, Di. Xiao, _et al._ , “Direct observation of two-dimensional magnons in atomically thin CrI3,” Nature Physics 17, 20–25 (2021).
  * Jin  _et al._ (2018) W. Jin, H. Ho Kim, Z. Ye, S. Li, P. Rezaie, F. Diaz, S. Siddiq, E. Wauer, B. Yang, C. Li, _et al._ , “Raman fingerprint of two terahertz spin wave branches in a two-dimensional honeycomb ising ferromagnet,” Nature communications 9, 1–7 (2018).
  * Lado and Fernández-Rossier (2017) J. L. Lado and J. Fernández-Rossier, “On the origin of magnetic anisotropy in two dimensional CrI3,” 2D Materials 4, 035002 (2017).
  * Xu  _et al._ (2018) C. Xu, J. Feng, H. Xiang, and L. Bellaiche, “Interplay between kitaev interaction and single ion anisotropy in ferromagnetic CrI3 and CrGeTe3 monolayers,” npj Computational Materials 4, 1–6 (2018).
  * Xu  _et al._ (2020) C. Xu, J. Feng, S. Prokhorenko, Y. Nahas, H. Xiang, and L. Bellaiche, “Topological spin texture in janus monolayers of the chromium trihalides cr (i, x) 3,” Physical Review B 101, 060404 (2020).
  * Liu  _et al._ (2018) J. Liu, M. Shi, J. Lu, and M. P. Anantram, “Analysis of electrical-field-dependent dzyaloshinskii-moriya interaction and magnetocrystalline anisotropy in a two-dimensional ferromagnetic monolayer,” Physical Review B 97, 054416 (2018).
  * Behera  _et al._ (2019) A. K. Behera, S. Chowdhury, and S. R. Das, “Magnetic skyrmions in atomic thin CrI3 monolayer,” Applied Physics Letters 114, 232402 (2019).
  * Sabani  _et al._ (2020) D. Sabani, C. Bacaksiz, and M. V. Milosevic, “Releasing latent chirality in magnetic two-dimensional materials,” arXiv preprint arXiv:2009.03421 (2020).
  * Xiao  _et al._ (2021) F. Xiao, K. Chen, and Q. Tong, “Magnetization textures in twisted bilayer CrX3 (X= Br, I),” Physical Review Research 3, 013027 (2021).
  * Tong  _et al._ (2018) Q. Tong, F. Liu, J. Xiao, and W. Yao, “Skyrmions in the moiré of van der waals 2d magnets,” Nano Letters 18, 7194–7199 (2018).
  * Chumak  _et al._ (2017) A. V. Chumak, A. A. Serga, and B. Hillebrands, “Magnonic crystals for data processing,” Journal of Physics D: Applied Physics 50, 244001 (2017).
  * Šabani  _et al._ (2020) D. Šabani, C. Bacaksiz, and M. V. Milošević, “Ab initio methodology for magnetic exchange parameters: Generic four-state energy mapping onto a heisenberg spin hamiltonian,” Physical Review B 102, 014457 (2020).
  * McGuire  _et al._ (2015) M. A. McGuire, H. Dixit, V. R. Cooper, and B. C. Sales, “Coupling of crystal structure and magnetism in the layered, ferromagnetic insulator cri3,” Chemistry of Materials 27, 612–620 (2015).
  * Kim  _et al._ (2019) M. Kim, P. Kumaravadivel, J. Birkbeck, W. Kuang, S. G. Xu, D. G. Hopkinson, J. Knolle, P. A. McClarty, A. I. Berdyugin, M. B. Shalom, _et al._ , “Micromagnetometry of two-dimensional ferromagnets,” Nature Electronics 2, 457–463 (2019).
  * Müller  _et al._ (2019) G. P. Müller, M. Hoffmann, C. Dißelkamp, D. Schürhoff, S. Mavros, M. Sallermann, N. S. Kiselev, H. Jónsson, and S. Blügel, “Spirit: Multifunctional framework for atomistic spin simulations,” Physical Review B 99, 224414 (2019).
  * Gruszecki  _et al._ (2016) P. Gruszecki, M. Kasprzak, A. E. Serebryannikov, M. Krawczyk, and W. Śmigaj, “Microwave excitation of spin wave beams in thin ferromagnetic films,” Scientific Reports 6, 1–8 (2016).
  * Song  _et al._ (2019) M. Song, K.-W. Moon, C. Hwang, and K.-J. Kim, “Omnidirectional spin-wave array antenna,” Physical Review Applied 11, 024027 (2019).
  * (26) “See supplemental information at [url].” .
  * Klein  _et al._ (2018) D. R. Klein, D. MacNeill, J. L. Lado, D. Soriano, E. Navarro-Moratalla, K. Watanabe, T. Taniguchi, S. Manni, P. Canfield, J. Fernández-Rossier, _et al._ , “Probing magnetism in 2d van der waals crystalline insulators via electron tunneling,” Science 360, 1218–1222 (2018).
  * Chen  _et al._ (2017) Y.-J. Chen, Kh. Zakeri, A. Ernst, H. J. Qin, Y. Meng, and J. Kirschner, “Group velocity engineering of confined ultrafast magnons,” Physical Review Letters 119, 267201 (2017).
  * Qin  _et al._ (2019) H. J. Qin, S. Tsurkan, A. Ernst, and Kh. Zakeri, “Experimental realization of atomic-scale magnonic crystals,” Physical Review Letters 123, 257202 (2019).
  * Chen  _et al._ (2018) L. Chen, J.-H. Chung, B. Gao, T. Chen, M. B. Stone, A. I. Kolesnikov, Q. Huang, and P. Dai, “Topological spin excitations in honeycomb ferromagnet CrI3,” Physical Review X 8, 041028 (2018).
  * Rezende (2020) S. M. Rezende,  _Fundamentals of Magnonics_ , Vol. 969 (Springer, 2020).
  * Zhang  _et al._ (2021) S. Zhang, R. Xu, N. Luo, and X. Zou, “Two-dimensional magnetic materials: structures, properties and external controls,” Nanoscale 13, 1398–1424 (2021).
  * Torelli  _et al._ (2019) D. Torelli, K. S. Thygesen, and T. Olsen, “High throughput computational screening for 2D ferromagnetic materials: the critical role of anisotropy and local correlations,” 2D Materials 6, 045018 (2019).
  * Cortés-Ortuño and Landeros (2013) D. Cortés-Ortuño and P. Landeros, “Influence of the dzyaloshinskii–moriya interaction on the spin-wave spectra of thin films,” Journal of Physics: Condensed Matter 25, 156001 (2013).
  * Wang  _et al._ (2019) K. Wang, A. A Puretzky, Z. Hu, B. R. Srijanto, X. Li, N. Gupta, H. Yu, M. Tian, M. Mahjouri-Samani, X. Gao, _et al._ , “Strain tolerance of two-dimensional crystal growth on curved surfaces,” Science advances 5, eaav4028 (2019).
  * Peng  _et al._ (2020) Z. Peng, X. Chen, Y. Fan, D. J. Srolovitz, and D. Lei, “Strain engineering of 2d semiconductors and graphene: from strain fields to band-structure tuning and photonic applications,” Light: Science & Applications 9, 1–25 (2020).
  * Castellanos-Gomez  _et al._ (2013) A. Castellanos-Gomez, R. Roldán, E. Cappelluti, M. Buscema, F. Guinea, H. S. J. van-der Zant, and G. A. Steele, “Local strain engineering in atomically thin mos2,” Nano letters 13, 5361–5366 (2013).
  * Sun and Liu (2019) Y. Sun and K. Liu, “Strain engineering in functional 2-dimensional materials,” Journal of Applied Physics 125, 082402 (2019).
  * Khestanova  _et al._ (2016) E. Khestanova, F. Guinea, L. Fumagalli, A. K. Geim, and I. V. Grigorieva, “Universal shape and pressure inside bubbles appearing in van der waals heterostructures,” Nature Communications 7, 1–10 (2016).
  * Dai  _et al._ (2019) Z. Dai, L. Liu, and Z. Zhang, “Strain engineering of 2d materials: issues and opportunities at the interface,” Advanced Materials 31, 1805417 (2019).
  * Quereda  _et al._ (2016) J. Quereda, P. San-Jose, V. Parente, L. Vaquero-Garzon, A. J. Molina-Mendoza, N. Agraït, G. Rubio-Bollinger, F. Guinea, R. Roldán, and A. Castellanos-Gomez, “Strong modulation of optical properties in black phosphorus through strain-engineered rippling,” Nano letters 16, 2931–2937 (2016).
  * Yang  _et al._ (2021) S. Yang, Y. Chen, and C. Jiang, “Strain engineering of two-dimensional materials: Methods, properties, and applications,” InfoMat 3, 397–420 (2021).
  * Deng  _et al._ (2018) S. Deng, A. V. Sumant, and V. Berry, “Strain engineering in two-dimensional nanomaterials beyond graphene,” Nano Today 22, 14–35 (2018).
  * Jiang  _et al._ (2019) J. Jiang, T. Xu, J. Lu, L. Sun, and Z. Ni, “Defect engineering in 2D materials: precise manipulation and improved functionalities,” Research 2019 (2019).
  * Schuler  _et al._ (2019) B. Schuler, D. Y. Qiu, S. Refaely-Abramson, C. Kastl, C. T. Chen, S. Barja, R. J. Koch, D. F. Ogletree, S. Aloni, A. M. Schwartzberg, _et al._ , “Large spin-orbit splitting of deep in-gap defect states of engineered sulfur vacancies in monolayer WS2,” Physical Review Letters 123, 076801 (2019).
  * Mulkers  _et al._ (2018) J. Mulkers, B. Van W., and M. V. Milošević, “Tunable snell’s law for spin waves in heterochiral magnetic films,” Physical Review B 97, 104422 (2018).
  * Wang  _et al._ (2018) Z. Wang, B. Zhang, Y. Cao, and P. Yan, “Probing the dzyaloshinskii-moriya interaction via the propagation of spin waves in ferromagnetic thin films,” Physical Review Applied 10, 054018 (2018).
  * Macke and Goll (2010) S. Macke and D. Goll, “Transmission and reflection of spin waves in the presence of néel walls,”  _Journal of Physics: Conference Series_ , 200, 042015 (2010).
  * Wang  _et al._ (2013) X-G Wang, G-H Guo, G-F Zhang, Y-Z Nie, and Q-L Xia, “Spin-wave resonance reflection and spin-wave induced domain wall displacement,” Journal of Applied Physics 113, 213904 (2013).
  * de Vries  _et al._ (2021) F.t K de Vries, E. Portolés, G. Zheng, T. Taniguchi, K. Watanabe, T. Ihn, K. Ensslin, and P. Rickhaus, “Gate-defined josephson junctions in magic-angle twisted bilayer graphene,” Nature Nanotechnology , 1–4 (2021).
  * Mao  _et al._ (2020) J. Mao, S. P. Milovanović, M. Anđelković, X. Lai, Y. Cao, K. Watanabe, T. Taniguchi, L. Covaci, F. M. Peeters, A. K. Geim, _et al._ , “Evidence of flat bands and correlated states in buckled graphene superlattices,” Nature 584, 215–220 (2020).
  * Wang  _et al._ (2016) X. Wang, Q. Weng, Y. Yang, Y. Bando, and D. Golberg, “Hybrid two-dimensional materials in rechargeable battery applications and their microscopic mechanisms,” Chemical Society Reviews 45, 4042–4073 (2016).
  * Zollner  _et al._ (2019) K. Zollner, P. E. F. Junior, and J. Fabian, “Proximity exchange effects in MoSe2 and WSe2 heterostructures with CrI3: Twist angle, layer, and gate dependence,” Physical Review B 100, 085128 (2019).
  * Zhu  _et al._ (2011) Z. Y. Zhu, Y. C. Cheng, and U. Schwingenschlögl, “Giant spin-orbit-induced spin splitting in two-dimensional transition-metal dichalcogenide semiconductors,” Physical Review B 84, 153402 (2011).
  * Reyes-Retana and Cervantes-Sodi (2016) J. A. Reyes-Retana and F. Cervantes-Sodi, “Spin-orbital effects in metal-dichalcogenide semiconducting monolayers,” Scientific Reports 6, 1–10 (2016).


# Predicting Surface Heat Flux on Complex Systems via Conv-LSTM

Yinpeng Wang
School of Electronic and Information
Beihang University
Beijing, 100191
<EMAIL_ADDRESS>
Nianru Wang
School of Electronic and Information
Beihang University
Beijing, 100191
<EMAIL_ADDRESS>
Qiang Ren
School of Electronic and Information
Beihang University
Beijing, 100191
<EMAIL_ADDRESS>

###### Abstract

Existing algorithms with iterations as the principle for 3D inverse heat
conduction problems (IHCPs) are usually time-consuming. With the recent
advancements in deep learning techniques, it is possible to apply the neural
network to compute IHCPs. In this paper, a new framework based on
Convolutional-LSTM is introduced to predict the transient heat flux via
measured temperature. The inverse heat conduction models concerned in this
work have 3D complex structures with non-linear boundary conditions and
thermophysical parameters. In order to reach high precision, a forward solver
based on the finite element method is utilized to generate sufficient data for
training. The fully trained framework can provide accurate predictions
efficiently once the measured temperature and models are acquired. It is
believed that the proposed framework offers a new pattern for real-time heat
flux inversion.

##  1 Introduction

Real-time identifications of the heat flux distribution at a given domain are
essential in many fields like engine design, geological prospecting and so on.
Direct ways of heat flux measurements to the external surface of a hypersonic
aircraft are expensive and difficult to maintain a high accuracy due to the
severe condition [1]. As a result, there is a strong demand to develop a more
precise and more efficient measurement technique for determining the heat flux
distribution mediately. One reasonable approach is to take advantage of the
numerous temperature measurement on the interior surface of the capsule to
invert the heat flux on the exterior surface. Measuring the temperature is
more economical and robust than directly measuring the heat flux. Hence, we
need to solve the inverse heat conduction problems (IHCPs) to acquire the
anticipated heat flux.

An IHCP generally refers to the problems in identifying the unknown boundary
or the thermal properties in a thermal system. As this problem is super
sensitive to the measurement variation, it is mathematically ill-posed [2].
Therefore, many valuable techniques have been proposed to transform it into a
relatively well-posed problem. Traditional algorithms used in the field of
IHCPs include sequential function specification method (SFSM) [3], iterative
techniques [4, 5, 6, 7, 8], Tikhonov regularization [9, 10], singular value
decomposition with model-reduction [11], digital filter method [1, 12, 13],
genetic algorithm [14, 15], Bayesian framework [16], neural networks [17, 18,
19, 20, 21, 22, 23] and so on.

Many of the inverse techniques mentioned above perform well in heat conduct
systems with simple geometries (1D or 2D models) and linear thermal
parameters. However, in many practical cases, the problems are usually more
sophisticated. Some works about those IHCPs with nonlinearity have been
published. A regular way is to divide the nonlinear model into several
linearized sub models. Keanini et al. [24] proposed a method based on the non-
iterative finite element algorithms and successfully reconstruct the time-
varying heat flux. The sequential function specification algorithm was adopted
in the study of Huang et al. [10] predict the spatial and temporal related
heat flux distribution by supposing that changes are linear at every time step
while the sensitivity coefficients were obtained by finite element method
calculations. Lv et al. [25] proposed a multiple model adaptive inverse method
to forecast the heat flux condition of a nonlinear IHCP. In fact, the accuracy
of sub-model algorithms relies on how many sub-models are constructed. It is
hard for these techniques to solve IHCPs with severe nonlinearity, as they
should use a considerable amount of sub-models and may cost unaffordable
computational resources. A novel approximate Bayesian computation algorithm is
proposed by Zeng et al. [16] for both linear and non-linear systems, while
massive samples are needed. Therefore, it is hard for this method to solve
problems within seconds. Duda et al. [26, 27]. efficiently solved the
nonlinear IHCPs by dividing a complex domain into sample parts. However, the
cases are limited to 2D models and the number of sensors needs to be specified
in different cases.

In industrial applications, the nonlinearity of the IHCP are more severe due
to the following reasons: (1) The geometry of the model can be three-
dimensional with intricate shapes and curved surfaces. In most practical
cases, it is not appropriate to simplify the model into 1D or 2D models as it
will lead to unbearable errors. (2) The nonlinearity of the model should be
considered, especially the sophisticated radiation boundary condition and
temperature-dependent materials. (3) The heat flux computation should be
efficient, especially in a long-period inversion process. In such a case, the
iteration-based methods with finite element methods solving forward problem
[5, 7] are improper because of their unbearable time cost. In summary,
handling IHCPs involving all the above aspects are much more meaningful but
more complicated. Therefore, an efficient inverse algorithm is proposed to
solve IHCPs with those characteristics.

The neural network has two advantages compared to the traditional method. One
is that a fully trained framework spends less time to forecast the heat flux.
The other is that the neural network can simulate the nonlinearity of the IHCP
better than the traditional methods. Most of the studies have applied the
artificial neural network (ANN) to solve IHCPs and used a backpropagation
algorithm to optimize the framework. Ghadimi et al. [19] combined the
Sequential Function Specification method and ANNs to handle the heat flux
inversion of a locomotive brake disc. Via two layers of the fully connected
network, he succeeded to forecast the linearized heat flux which is
homogeneously distributed. The temporal-related heat flux is also predicted
via perturbed temperature by Cortés et al. [20]. They also used the fully
connected network to capture the temporal features of the IHCP and simplified
the question by assuming that the heat flux is homogenously located. However,
in real industrial scenarios, it is impractical to utilize fully connected
networks to solve complex IHCPs since the trainable parameters are too large.
Thanks to the development of the research on deep learning, some special
neural networks have been proposed. The long short-term memory networks (LSTM)
proposed by Hochreite et al. [28] is an efficient framework handling these
temporal-related systems. LSTM has been proven to be stable and powerful for
modeling long-range dependencies in various precious studies as a special RNN
structure [29, 30, 31]. One superiority of employing the memory cell and gates
to bridle information flow is that the gradient will be saved in the cell and
be protected from fading too quickly while this is an essential problem for
the vanilla RNN model [29, 30]. Multiple LSTMs can be coded and temporally
linked to propose more complex structures, which have already been applied to
solve many real-time sequence modeling problems [32]. However, such a network
can only accept 1D tensors, so it could only solve problems that contain
information about time. Very recently, a new method is proposed by Shi et al.
[33] for precipitation nowcasting which modifies the common LSTM and puts the
convolution into the layer to change the input tensor and adds the spatial
information. Hence, the convolutional LSTM (ConvLSTM) is introduced as the
base of our framework.

This paper is organized as follows. The second section firstly introduces the
physical model of the IHCP to be investigated. Then the workflow and the
method to generate the data set are illustrated. Furthermore, the architecture
of the framework is also exhibited in detail. The third section displays and
analyzes the outcome of the proposed framework. A brief conclusion is provided
in the last section.

Table 1: Nomenclature Heat conduction |  | $\left[{{K_{1}}^{(e)}}\right]$ | Coefficient matrix
---|---|---|---
$\Omega$ | Solid body domain | ${\left\\{\overline{P}\right\\}}$ | Load matrix
$r$ | Space coordinate | ConvLSTM |
${{S}_{1,2,3}}$ | Boundary surface | ${{\chi}_{1}},\cdots,{{\chi}_{t}}$ | Inputs of the ConvLSTM
${{T}_{0}}$ | Initial temperature | ${{C}_{1}},\cdots,{{C}_{t}}$ | Cell outputs of the ConvLSTM
$t$ | Time | ${{H}_{1}},\cdots,{{H}_{t}}$ | Hidden states of the ConvLSTM
$q$ | Heat flux | ${{i}_{t}},{{f}_{t}},{{o}_{t}}$ | Inside gates of the ConvLSTM
${{T}_{c}}$ | Measured temperature | $W$ | Coefficient of the ConvLSTM
$k(T)$ | Thermal conductivity | Subscript |
${{C}_{p}}(T)$ | Heat capacity | $j$ | Sensor index
$\rho$ | Mass density | $w$ | Temperature on the radiation boundary
$\sigma$ | Stefan Boltzmann constant | $x$ | Input index
$\varepsilon$ | Emissivity | $i,f,o$ | Gates index
FEM |  | $H$ | Hidden states index
$N_{i}$ | Weight of the $i_{th}$ node | $C$ | Cell output index

##  2 Methodology and Formulation

The whole workflow of the inversion process is demonstrated in Fig. 1. Here,
we first generate various geometries in MATLAB and combine the thermal
properties with the geometries to produce the physical models. Then, the FEM
solver utilizes the physical models with the heat flux generated by MATLAB to
calculate the temperature distribution. After that, the generated geometries
and temperatures in each time step are packed together to serve as the input
of the framework. The framework is then optimized in the training process
where the output is supervised by the generated heat flux.

Figure 1: The workflow to train and test the framework.

###  2.1 Nonlinear heat conduction model

The developed method aims to forecast a time-varying heat flux distribution of
a nonlinear 3D heat conduction model with temperature-dependent heat
conduction parameters. Fig. 2 displays the schematic of a typical example of
such a model, where the solid body is a isotropic homogeneous domain denoted
as $\Omega$. $S_{1}$ is a curved surface where the unknown temporal and
spatial varying heat flux $q(r,t)$ is imposed. The bottom surface is
designated as $S_{2}$ and the physical boundary is radiation boundary
conditions. The rest of the surfaces are insulated and symbolized as $S_{3}$.
The temperature sensors are uniformly scattered on the temperature surface
$S_{2}$ which are denoted as $T_{j}\ (1\leq j\leq I)$, where $I$ is the total
number of the sensors.

Figure 2: A typical example of the 3D nonlinear heat conduction systems.

The governing equations, related boundaries and initial conditions of the
system are

$k\nabla^{2}T=\rho C_{p}\frac{\partial T}{\partial t},r\in\Omega,t\geq 0$ (1)

$-k\frac{\partial T(r;t)}{\partial z}=q(r;t),r\in{{S}_{1}},t\geq 0$ (2)

$-k\frac{\partial T(r;t)}{\partial
z}=\sigma\varepsilon\left(T_{w}^{4}(r;t)-T_{0}^{4}\right),r\in{{S}_{2}},t\geq
0$ (3)

$-k\frac{\partial T(r;t)}{\partial n}=0,r\in{{S}_{3}},t\geq 0$ (4)

$T(r;0)={{T}_{0}},r\in\Omega,t=0$ (5)

where $T(r;t)$ represents the temperature and $\rho$ designates the mass
density of the solid body. Temperature-dependent heat capacity and thermal
conductivity are designated as ${{C}_{p}}(T)$ and $k(T)$. And $\varepsilon$ is
the emissivity while $\sigma$ is the Stefan Boltzmann Constant.

The physic boundary settings, the detailed geometry parameters and thermal
properties expressed in Tabel. 1 will be the foundation of the test case. Fig.
3 shows the temperature-dependent thermal conductivity $k(T)$ and heat
capacity ${{C}_{p}}(T)$ of an alloy used on the aircraft. Both the initial
temperature of the solid body and ambient temperature are 273.15K. 900
temperature sensors are uniformly located on the bottom surface $S_{2}$.

Figure 3: The temperature-dependent thermal conductivity and heat capacity.

###  2.2 Data set generation

After constructing the 3D nonlinear model, we evolve COMSOL based on the
finite element method to generate sufficient data for training and testing.
The detailed procedure of this algorithm is as follows: Firstly, the 3D model
is divided into $E$ sub grids, each of which contains $p$ nodes. The
temperature at any given position $(x,y,z)$ can be expressed approximately by
the interpolation of node temperature.

$T=\sum\limits_{i=1}^{p}{{N_{i}}(x,y,z){T_{i}}^{e}=\left[N\right]{{\left\\{T\right\\}}^{\left(e\right)}}}$
(6)

where $N_{i}$ is the weight of the $i_{th}$ node. Applying the functional
theory, the heat conduction equation can be rewritten as

$\left[{{K_{1}}^{(e)}}\right]{\left\\{T\right\\}^{\left(e\right)}}+\left[{{K_{2}}^{(e)}}\right]{\left\\{T\right\\}^{\left(e\right)}}+\left[{{K_{3}}^{(e)}}\right]{\left\\{{\dot{T}}\right\\}^{\left(e\right)}}-{\left\\{P\right\\}^{\left(e\right)}}=0$
(7)

where $\left[{{K_{1}}^{(e)}}\right]$, $\left[{{K_{2}}^{(e)}}\right]$,
$\left[{{K_{3}}^{(e)}}\right]$ and $\left\\{P\right\\}^{\left(e\right)}$
denote the heat conduction, the boundary condition, the unsteady and the load
matrix, respectively. By assembling each sub matrix, it can be obtained that

$\left[K\right]\left\\{{\overline{T}}\right\\}=\left\\{{\overline{P}}\right\\}$
(8)

where the coefficient matrix are

$\left[K\right]=\sum\limits_{e=1}^{E}{\left({\left[{{K_{1}}^{\left(e\right)}}\right]+\left[{{K_{2}}^{\left(e\right)}}\right]}+\left[{{K_{3}}^{\left(e\right)}}\right]\right)}$
(9)

${\left\\{\overline{P}\right\\}}=\sum\limits_{e=1}^{E}{{{\left\\{P\right\\}}^{\left(e\right)}}}$
(10)

After imposing (2)-(5) to (8), the temperature at the $S_{2}$ can be acquired.
With various models generated by MATLAB, it is feasible to produce sufficient
data for training and testing.

###  2.3 Architecture of the framework

The 3D mode in our model has two tough points: (1) there are 900 sampling
points of heat flux to be predicted, so that it is inefficient to use fully
connected-LSTM to predict one by one; (2) there are various shapes of the
solid body, so that the temperature distribution is influenced by both the
heat flux and the body shapes. To address these problems, the ConvLSTM [33],
containing both convolutional structures in input-to-state and state-to-state
transitions, is introduced to the framework. By assembling multiple ConvLSTM
layers with LSTM and fully connected layers, we could build a network model
for 3D IHCPs.

A major drawback of fully connected-LSTM in handling spatiotemporal data is
that it uses only fully connected layers in transitions so that it does not
encode any spatial information. However, since all the inputs
${{\chi}_{1}},\cdots,{{\chi}_{t}}$, cell outputs ${{C}_{1}},\cdots,{{C}_{t}}$,
hidden states ${{H}_{1}},\cdots,{{H}_{t}}$ and gates
${{i}_{t}},{{f}_{t}},{{o}_{t}}$ of the ConvLSTM are 3D tensors , it is
feasible to store different spatial information in different channels. This is
achieved by using a convolution operator in the transitions. The key equations
of ConvLSTM are shown below, where ‘$*$’ denotes the convolution operator and
‘$\circ$’ denotes the Hadamard product:

${{i}_{t}}=\sigma({{W}_{xi}}*{{\chi}_{t}}+{{W}_{hi}}*{{H}_{t-1}}+{{W}_{ei}}\circ{{C}_{t-1}}+{{b}_{i}})$
(11)

${{f}_{t}}=\sigma({{W}_{xf}}*{{x}_{t}}+{{W}_{hi}}*{{H}_{t-1}}+{{W}_{cf}}\circ{{c}_{t-1}}+{{b}_{f}})$
(12)

${{C}_{t}}={{f}_{t}}\circ{{C}_{t-1}}+{{i}_{t}}\circ\tanh({{W}_{xc}}*{{x}_{t}}+{{W}_{hc}}*{{h}_{t-1}}+{{b}_{c}})$
(13)

${{o}_{t}}=\sigma({{W}_{xo}}*{{x}_{t}}+{{W}_{ho}}*{{h}_{t-1}}+{{W}_{co}}\circ{{c}_{t}}+{{b}_{o}})$
(14)

${{H}_{t}}={{o}_{t}}\circ\tanh({{C}_{t}})$ (15)

ConvLSTM can serve as a layer for a more complex framework to capture spatial-
temporal characteristics in practical scenarios. In this work, the proposed
framework shown in Fig. 4 with three layers of ConvLSTM, three layers of FC-
LSTM and one layer of fully connected network is applied to solve the complex
IHCPs. By minimizing the kernel size gradually in different layers, this
structure could capture the different sizes of spatial relationships
efficiently. We set two kernel filters of the first layers as $4\times 4$ to
separately keep the geometry shapes and temperature distributions. Considering
the nonlinearity of the system, we add an extra layer in which the kernel
filter size is $2\times 2$ to better fit the spatial influence of the heat
flux. And the last layer of ConvLSTM is to unify the result of two filters
into one. By minimizing the kernel filters and unifying two channels, the
output of the last ConvLSTM layer is a combination of geometry shape and
temperature. Therefore, the input of LSTM could be imagined as a 3D
temperature which is similar to the heat flux in spatial so that we only need
to capture the temporal information. The first two layers of LSTM are set to
return the sequence in each time step, as this will help find the relationship
between the temperature and heat flux more precisely. At the last layer, the
LSTM is used to set the sensitive coefficient for every time step to make sure
the main heat flux influence the temperature at this time. To avoid missing
the information included in former layers, the number of the neurons is the
same as that of the sampling points of heat flux. After the LSTM, we also add
a fully connected layer to resize the output of the LSTM as it may be unable
to match the supervised data. Since the network has multiple stacked ConvLSTM
layers, it has striking power to give predictions in complex dynamic systems
like IHCP in 3D complex models.

Figure 4: Architecture of the framework in solving inverse heat conduction
problems.

###  2.4 Training process

The training process is performed on an NVIDIA RTX 2080 Ti graphic card on
Dell Precision 7920 Tower. During the process, the loss of the training set
(training loss) and that of the test set (validation loss) decrease gradually.
After 300 iterations, the loss reduced to an acceptable level. The loss curves
are exhibited in Fig. 5.

Figure 5: The loss curve of the training process.

##  3 Results

In this part, some specimens with various heat flux distributions are chosen
from the complex 3D data set to assess the performance of the proposed
framework. Here, the ground truth of heat flux distribution is defined by the
data generated through MATLAB.

###  3.1 Regular 3D models

Before starting to handle 3D complex models, we first build up a simple 3D
model with flat surfaces to evaluate the ability of the proposed framework in
capturing the temporal and spatial features. In such models, as the heat flux
surface and the temperature measured surface are parallel, the spatial
distribution of the temperature is similar to that of the heat flux and this
kind of models is denoted as regular. The main difficulty is capturing the
temporal features. Here we use three heat flux cases to test the network and
the outcome is visualized in Fig. 6. We randomly choose two different time
steps to assess the framework and the result exhibits that the difference in
each model is slightly affected by the time steps which demonstrates the
proposed framework is good at predicting the temporal variation of the heat
flux. To show the precise difference between the ground truth of heat flux and
predicted ones, we add the distribution of relative error in each case which
is defined as

$error=\left|\frac{{{q}_{predicted}}-{{q}_{real}}}{{{q}_{real}}}\right|$ (16)

Figure 6: Samples in the regular 3D test set. (a) The predicted heat flux, (b)
the ground truth of heat flux and (c) the relative errors.

###  3.2 Complex 3D models

After the 3D regular models, we change the training data of the framework and
measure the ability of the ConvLSTM to reconstruct the distribution of heat
flux combining the features of complex 3D models and temperature
distributions. In this kind of models which is symbolized as complex, the heat
flux surface is no longer parallel to the temperature measured surface and
contains curved surface. To demonstrate the performance of the ConvLSTM
framework in a better way, we randomly choose twelve cases to exhibit the
temporal and spatial inversion ability separately. Fig.6 shows three cases
that include different complex models and heat flux distributions. The
geometry shapes of the models varies in many different aspects. The max
height, min height and curvature are different while the relative error is
stable. Therefore, we could figure out that the proposed framework performed
pretty well on complex 3D models and have a good generalization ability in
geometry shapes.

Figure 7: Samples in the complex 3D test set. (a) The geometry of the complex
models, (b) the forecast of the heat flux and the measured temperature (t=70s)
(c) the ground truth of the heat flux and the measured temperature (t=70s) and
(d) the relative errors.

In Fig. 8, there are three complex models and three different heat fluxes for
each model. These heat fluxes are unevenly distributed in space, and we
randomly select three points to illustrate the temporal fitting capability of
the framework. The three points are located on the diagonal of the surface and
have different ranges of the heat flux variation. The result shows that the
range of the heat flux variation and location of the point have few influence
on the prediction and demonstrates that our framework has a good
generalization ability in predicting temporal and spatial related heat flux.

Figure 8: Samples in the complex 3D test set. (a) The geometry of the complex
models, (b), (c) and (d) three different heat flux distributions imposed on
(a). $R_{i}(i=1,2,3)$: ground truth, $P_{i}(i=1,2,3)$: prediction.

###  3.3 Calculation accuracy and speed

We compute the relative error for the test case in Fig. 9. As mentioned
before, the regular refers to the regular models of which the two surfaces are
parallel while the complex denote the complex models where two surfaces are
not parallel. Since the regular models offer more information about the heat
flux, the difficulty of the IHCP for regular ones is less and the loss
distribution is better than the complex ones. These figures illustrate that
for regular 3D models it is easier to predict the heat flux more precisely and
stable.

Figure 9: Statistical distribution for test cases separately produced by two
kinds of generators. The regular one represents the flat heat flux surface
while the complex one represents the complex surface.

Also, since the framework is based on the GPU and CUDA is employed to forecast
the heat flux, the average time cost to predict a case which contains 71 time
steps and represents a 35s period is only 1.27s. In other words, for each time
step, it takes time within 20 ms.

##  4 Conclusion

In this work, we have proposed a new framework to forecast the heat flux of
complex 3D models, which have temperature-dependent thermal parameters and
nonlinear physical boundaries. We introduce a new variation of ConvLSTM which
is designed to solve spatial and temporal-related problems. Compared to the
regular iterative methods, such a framework can work efficiently after being
trained properly. As it usually uses CUDA to calculate in parallel, it could
figure out the heat flux preciously within 20 ms, which is extremely superior
to those traditional algorithms. Also, as this framework combines ConvLSTM
with LSTM, it takes few epochs to reach convergence. The result illustrate
that our framework have a good generalization ability in forecasting temporal
and spatial related heat flux. Hence, the proposed framework can replace the
direct heat flux measurements in real-time heat flux inversion problems with
advanced temperature measurements and have multiple applications in many
industrial scenarios like engine design or furnace operation in the future.

## Acknowledgments

This work is supported in part by the National Nature Science Foundation
(Grant No. 61801009), in part by the National Science and Technology Major
Project (Grant No. 2019-VIII-0009-0170), in part by Defense Industrial
Technology Development Program (Grant No. JCKY2020601B011) and in part by
Stable Operation Project.

## References

  * [1] Hamidreza Najafi and Keith A Woodbury.  Online heat flux estimation using artificial neural network as a digital filter approach.  International Journal of Heat and Mass Transfer, 91:808–817, 2015\.
  * [2] James V Beck, Ben Blackwell, and Charles R St Clair Jr.  Inverse heat conduction: Ill-posed problems.  James Beck, 1985.
  * [3] S Chantasiriwan.  Comparison of three sequential function specification algorithms for the inverse heat conduction problem.  International communications in heat and mass transfer, 26(1):115–124, 1999.
  * [4] OM Alifanov.  Solution of an inverse problem of heat conduction by iteration methods.  Journal of Engineering Physics, 26(4):471–476, 1974.
  * [5] Mattia Bergagio, Haipeng Li, and Henryk Anglart.  An iterative finite-element algorithm for solving two-dimensional nonlinear inverse heat conduction problems.  International Journal of Heat and Mass Transfer, 126:281–292, 2018\.
  * [6] Miao Cui, Kai Yang, Xiao-liang Xu, Sheng-dong Wang, and Xiao-wei Gao.  A modified levenberg–marquardt algorithm for simultaneous estimation of multi-parameters of boundary heat flux by solving transient nonlinear inverse heat conduction problems.  International Journal of Heat and Mass Transfer, 97:908–916, 2016\.
  * [7] Piotr Duda.  A general method for solving transient multidimensional inverse heat transfer problems.  International Journal of Heat and Mass Transfer, 93:665–673, 2016\.
  * [8] JV Beck, B Litkouhi, and CR St Clair Jr.  Efficient sequential solution of the nonlinear inverse heat conduction problem.  Numerical Heat Transfer, Part A Applications, 5(3):275–286, 1982\.
  * [9] Andrey N Tikhonov and Vasiliy Y Arsenin.  Solutions of ill-posed problems.  New York, 1:30, 1977.
  * [10] Shuwen Huang, Bo Tao, Jindang Li, Yajun Fan, and Zhouping Yin.  Estimation of the time and space-dependent heat flux distribution at the tool-chip interface during turning using an inverse method and thin film thermocouples measurement.  The International Journal of Advanced Manufacturing Technology, 99(5):1531–1543, 2018.
  * [11] JR Shenefelt, R Luck, RP Taylor, and JT Berry.  Solution to inverse heat conduction problems employing singular value decomposition and model-reduction.  International Journal of Heat and Mass Transfer, 45(1):67–74, 2002\.
  * [12] Hongchu Chen, Jay I Frankel, and Majid Keyhani.  Nonlinear inverse heat conduction: Digitally filtered space marching with phase-plane and cross-correlation analyses.  Numerical Heat Transfer, Part B: Fundamentals, 72(2):109–129, 2017\.
  * [13] Farshad Kowsary, Morteza Mohammadzaheri, and Saeed Irano.  Training based, moving digital filter method for real time heat flux function estimation.  International communications in heat and mass transfer, 33(10):1291–1298, 2006.
  * [14] Balázs Czél and Gyula Gróf.  Inverse identification of temperature-dependent thermal conductivity via genetic algorithm with cost function-based rearrangement of genes.  International Journal of Heat and Mass Transfer, 55(15-16):4254–4263, 2012.
  * [15] Keith A Woodbury.  Neural networks and genetic algorithms in the solution of inverse problems.  Bulletim of Brazilian Society for Computing and Applied Mathematics, 2000.
  * [16] Yang Zeng, Hu Wang, Shuai Zhang, Yong Cai, and Enying Li.  A novel adaptive approximate bayesian computation method for inverse heat conduction problem.  International Journal of Heat and Mass Transfer, 134:185–197, 2019\.
  * [17] FJCP Soeiro, Patrícia Oliva Soares, HF Campos Velho, and AJ Silva Neto.  Using neural networks to obtain initial estimates for the solution of inverse heat transfer problems.  In Inverse Problems, Design an Optimization Symposium, pages 358–363, 2004.
  * [18] J Krejsa, KA Woodbury, JD Ratliff, and M Raudensky.  Assessment of strategies and potential for neural networks in the inverse heat conduction problem.  Inverse Problems in Engineering, 7(3):197–213, 1999.
  * [19] B Ghadimi, F Kowsary, and M Khorami.  Heat flux on-line estimation in a locomotive brake disc using artificial neural networks.  International Journal of Thermal Sciences, 90:203–213, 2015.
  * [20] Obed Cortés, Gustavo Urquiza, JA Hernandez, and Marco A Cruz.  Artificial neural networks for inverse heat transfer problems.  In Electronics, Robotics and Automotive Mechanics Conference (CERMA 2007), pages 198–201. IEEE, 2007.
  * [21] Ali Mirsepahi, Lei Chen, and Brian O’Neill.  A comparative artificial intelligence approach to inverse heat transfer modeling of an irradiative dryer.  International communications in heat and mass transfer, 41:19–27, 2013.
  * [22] Balázs Czél, Keith A Woodbury, and Gyula Gróf.  Simultaneous estimation of temperature-dependent volumetric heat capacity and thermal conductivity functions via neural networks.  International Journal of Heat and Mass Transfer, 68:1–13, 2014\.
  * [23] Balázs Czél, Keith A Woodbury, and Gyula Gróf.  Inverse identification of temperature-dependent volumetric heat capacity by neural networks.  International Journal of Thermophysics, 34(2):284–305, 2013.
  * [24] Russell G Keanini, Xianwu Ling, and Harischandra P Cherukuri.  A modified sequential function specification finite element-based method for parabolic inverse heat conduction problems.  Computational Mechanics, 36(2):117–128, 2005.
  * [25] Cai Lv, Guangjun Wang, Hong Chen, and Shibin Wan.  Estimation of the moving heat source intensity using the multiple model adaptive inverse method.  International Journal of Thermal Sciences, 138:576–585, 2019.
  * [26] Piotr Duda and Mariusz Konieczny.  A new algorithm for solving an inverse transient heat conduction problem by dividing a complex domain into parts.  International Journal of Heat and Mass Transfer, 128:865–874, 2019\.
  * [27] Piotr Duda and Mariusz Konieczny.  Solution of an inverse transient heat conduction problem in a part of a complex domain.  International Journal of Heat and Mass Transfer, 127:821–832, 2018\.
  * [28] Sepp Hochreiter and Jürgen Schmidhuber.  Long short-term memory.  Neural computation, 9(8):1735–1780, 1997.
  * [29] Alex Graves.  Generating sequences with recurrent neural networks.  arXiv preprint arXiv:1308.0850, 2013.
  * [30] Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio.  On the difficulty of training recurrent neural networks.  In International conference on machine learning, pages 1310–1318. PMLR, 2013.
  * [31] Ilya Sutskever, Oriol Vinyals, and Quoc V Le.  Sequence to sequence learning with neural networks.  arXiv preprint arXiv:1409.3215, 2014.
  * [32] Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhudinov, Rich Zemel, and Yoshua Bengio.  Show, attend and tell: Neural image caption generation with visual attention.  In International conference on machine learning, pages 2048–2057. PMLR, 2015.
  * [33] Xingjian Shi, Zhourong Chen, Hao Wang, Dit-Yan Yeung, Wai-Kin Wong, and Wang-chun Woo.  Convolutional lstm network: A machine learning approach for precipitation nowcasting.  arXiv preprint arXiv:1506.04214, 2015.