# Electrical control of magnons in multiferroic NiI2

Andrey Rybakov  Instituto de Ciencia Molecular (ICMol), Universitat de
València, Paterna 46980, Spain  Adolfo O. Fumega  Department of Applied
Physics, Aalto University, 02150 Espoo, Finland  Dorye L. Esteras  Instituto
de Ciencia Molecular (ICMol), Universitat de València, Paterna 46980, Spain
Jose L. Lado  Department of Applied Physics, Aalto University, 02150 Espoo,
Finland<EMAIL_ADDRESS>José J. Baldoví  Instituto de Ciencia Molecular
(ICMol), Universitat de València, Paterna 46980, Spain<EMAIL_ADDRESS>
###### Abstract

Layered van der Waals two-dimensional (2D) magnets are a cornerstone of
ultrathin spintronic and magnonic devices. The recent discovery of a 2D
multiferroic with strong magnetoelectric coupling in NiI2 offers a promising
platform for the electrical control of spin-wave transport. In this work,
using _ab initio_ calculations, we investigate how the magnonic properties of
monolayer NiI2 can be controlled using an external electric field. We show
that the emergence of a ferroelectric polarization leads to an energy
splitting in the magnon spectrum, thus establishing a way to detect the
electric polarization experimentally. We also show the modulation of the
magnon splitting and the energy position of the singularities in magnon DOS by
an electric field due to the strong magnetoelectric coupling. Our results
highlight the interplay between ferroelectricity and magnons in van der Waals
multiferroics and pave the way to design electrically tunable magnetic devices
at the 2D limit.

##  1 Introduction

Magnetism in two-dimensional (2D) materials provides a unique playground for
new fundamental phenomena and novel functional devices1, 2, 3, 4. Magnetic 2D
materials allow modification of their magnetic ground state and magnetic
excited states by applying mechanical strain 5, 6, creating van der Waals
heterostructures 7, electrostatic doping 8, 9, 10, intercalation 11, 12 and
molecular deposition 13. These systems range from conventional ferromagnets14,
15, quantum spin liquids16 and heavy-fermion Kondo insulators17, 18 dominated
by quantum fluctuations, and recently multiferroic 2D materials19, 20, 21
featuring entangled magnetism and ferroelectricity. Importantly, 2D materials
are prominent as building-blocks for the realization of the spin-wave circuits
at the limit of miniaturization, where spin waves (SW) are utilized as an
information carriers in contrast to conventional electronics 22, 23.

Multiferroic materials offer promising building-blocks to achieve electrical
control over magnon transport. This class of materials can display
simultaneously ferroelectric and magnetic order and they exhibit a strong
magnetoelectric coupling between the two, which makes them highly attractive
from a technological point of view 24, 25, 26, 27, 28, 29, 30. Special
attention has been recently given to van der Waals multiferroics due to the
enlarged possibilities that these compounds offer: a way to engineer novel
heterostructures with emergent functionalities 31, 32, 20. In this realm, NiI2
is the paradigmatic van der Waals multiferroic and is attracting a growing
attention 33, 34, 35, 36, 37, 19, 20. Its multiferroic order is a consequence
of non-collinear magnetic spin-spiral and the strong spin-orbit coupling
provided by the iodine atoms 31. Recently, it has been demonstrated that the
multiferroic order of NiI2 can survive in the few-layer and single-layer
limits 37, 19, 20. Therefore, establishing monolayer NiI2 as the first purely
2D multiferroic and adding a new system to the family of van der Waals
compounds to engineer artificial materials 38. However, a theoretical
investigation of the tunability of the magnon spectrum of multiferroic NiI2
with external electric fields has remained unexplored.

In this work, combining effective models and _ab initio_ calculations, we show
how the magnon spectrum of monolayer NiI2 can be controlled with an external
electric field. We observe that the C6 rotation symmetry breaking associated
with the emergence of ferroelectricity is directly manifested as a splitting
$\Theta$ in the magnon spectrum occurring at the rotational image of the
$\mathbf{q}$-vector of the spin spiral. We show that the size of the magnon
splitting can be controlled with an external electric field. Our results
highlight a way to experimentally estimate the magnitude of the electric
dipole through the splitting in the magnon spectrum and demonstrate that
multiferroic NiI2 is a promising platform for electrically tunable magnon
transport.

Figure 1:  (a) Schematic of the model used to study the magnons in
multiferroic NiI2. The structural unit cell is shown in black with lattice
parameter $a$ and Ni (I) atoms in gray (purple). The ferroelectric
displacements on the atoms are depicted as green arrows, in the direction
where the average electric polarization emerges. The first $J_{1}$ and third
$J_{3}$ neighbor magnetic exchange couplings are shown in red and blue
respectively. The ferroelectric distortion breaks the crystal’s C6 symmetry,
leading to two non-equivalent $J^{\prime}$ and $J^{\prime\prime}$ magnetic
exchanges for both – first and third – groups of neighbors. (b) Schematic of
the emergent ferroelectric polarization perpendicular to the plane formed by
the propagation vector $\mathbf{q}$ and the rotation vector $\mathbf{e}$ of
the spin spiral. The ferroelectric forces associated with the emergent
polarization are directed parallel to it. (c) Phonon band structure for
monolayer NiI2. The optical phonon mode associated with the ferroelectric
displacements $\omega_{FE}$ is highlighted in blue. (d,e) Equivalent spin
spiral configurations with opposite emergent electric polarization used to
calculate the ferroelectric forces. Commensurate $\sqrt{3}a\times\sqrt{21}a$
supercells are used to describe the spin spiral associated with the
$J_{3}/J_{1}$ ratio predicted by DFT calculations. The red arrows correspond
to the spins.

##  2 Results and discussion

###  2.1 Magnons in multiferroic NiI_2

We start by addressing the magnon spectrum in monolayer NiI2. The structural
unit cell of monolayer NiI2 is shown in Fig. 1a. The unpaired $d$ electrons of
Ni atoms lead to a spin state $S=1$ with ferromagnetic ($J_{1}$) and
antiferromagnetic ($J_{3}$) couplings between first and third neighbors
respectively. Due to the competition between the exchange couplings, monolayer
NiI2 develops a spin-spiral order at low temperatures. The spin spiral is
characterized by a propagation vector $\mathbf{q}$, that determines its
periodicity and direction, and rotation vector $\mathbf{e}$, that determines
the rotation plane of the spiral. In the case of a spin spiral dominated by
$J_{3}/J_{1}$ ratio, the propagation vector in the reciprocal space lies along
the $\Gamma$K segment (see Fig. 2a for definition of high symmetry points),
corresponding to a propagation vector in real space along the third neighbors’
direction (Fig. S1). In bulk NiI2, it has been shown that interlayer and
Kitaev interactions can lead to an out-of-plane component of the spin spiral
39, 40. However, in the monolayer limit, the magnetic order is dominated by
the $J_{3}/J_{1}$ ratio 20, 41. Therefore, its magnetic order can be described
with the following Heisenberg Hamiltonian

$\mathcal{H}=\sum_{ij}\mathcal{J}_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+\mathcal{H}^{\text{ani}}$
(1)

where the most important exchange couplings correspond to first $J_{1}$ and
third neighbors $J_{3}$ (see SI for details). The term
$\mathcal{H}^{\text{ani}}$ corresponds to the anisotropic exchange couplings
arising as a result of spin-orbit coupling, which have a minimal effect on the
spin-spiral q-vector.

Due to the strong spin-orbit coupling ($\lambda_{SOC}$) originating from the I
atoms 31, NiI2 develops an electric polarization $\mathbf{P}$ in the direction
perpendicular to the propagation and rotation vectors of the spin spiral 42.
Note that $\mathbf{P}$ corresponds to the average emergent polarization.
Experimentally it has been demonstrated that a modulated polarization with
half of the periodicity of the spin spiral emerges, allowing the atomic scale
visualization of the multiferroic order 20. For simplicity, we will consider
the average emergent polarization, although the model could be extended to
consider the effect of the ferroelectric modulation. The emergence of the
ferroelectric order is accompanied by emanating ferroelectric forces
($\mathbf{F}_{Ni}$ and $\mathbf{F}_{I}$ in Fig. 1b) parallel to the
polarization’s direction causing a ferroelectric displacement
$\boldsymbol{\Delta}_{Ni}$ and $\boldsymbol{\Delta}_{I}$ in each of the atoms
(as shown in the schematic of Fig. 1a). All the parameters entering this
process can be computed _ab initio_.

We start estimating the values of $J_{1}$ and $J_{3}$ in the undistorted
structure using density functional theory (DFT) 43, 44, 45, 46, 47, 48, 49,
50, 51, 52, 53 (see SI for computational details). Obtained values allow us to
determine the $\mathbf{q}$ vector of the spin spiral via minimization of the
energy within the first Brillouin zone

$E(\mathbf{q})\propto\sum_{ij}\mathcal{J}_{ij}\cos{\left(\mathbf{q}\cdot\mathbf{r}_{ij}\right)}$
(2)

where $\boldsymbol{r}_{ij}$ is the vector between two corresponding magnetic
centers. Spin vectors are considered to be in the $xy$ plane and the
rotational axis is assumed along the direction of the $z$ axis, due to the
isotropic nature of the Hamiltonian (eq. (1)). Due to the $C_{6}$ symmetry of
the undistorted structure, $6$ degenerate energy minima are exhibited in the
first Brillouin zone (Fig. 2a). With $J_{1}=-2.8$ meV and $J_{3}=2.32$ meV we
obtained $\mathbf{q}=(0.124,0.124,0)$ (written in the basis of reciprocal
vectors) along the $\Gamma$K path, which is in good agreement with recent
theoretical 41 and experimental 20 studies. The approximate commensurate spin-
spiral associated with the predicted $\mathbf{q}$ is a
$\sqrt{3}a\times\sqrt{21}a$ (Fig. 1d). When the spin-orbit coupling is
included in the DFT calculations, the ferroelectricity emerges indicated by
the ferroelectric forces $\mathbf{F}_{Ni}$ and $\mathbf{F}_{I}$ (Fig. 1b) that
induce the ferroelectric displacements $\boldsymbol{\Delta}_{Ni}$ and
$\boldsymbol{\Delta}_{I}$ in each atom. The ferroelectric forces can be
computed from DFT by considering two equivalent spin-spiral configurations
(Figs. 1d and 1e) with opposite $\mathbf{q}$ vectors31. The difference between
the emergent forces in each configuration allows us to factor out other
contributions to get $\mathbf{F}_{Ni}=(0,-8.8,0)$ meV/Å and
$\mathbf{F}_{I}=(0,4.4,0)$ meV/Å . We can estimate the ferroelectric
displacements in the absence of an electric field from the ferroelectric
forces as

$\boldsymbol{\Delta}_{\alpha}=\frac{\mathbf{F}_{\alpha}}{k_{\alpha}},$ (3)

where $k_{\alpha}$ are the elastic constants for each species
$\alpha=(\mathrm{Ni},\mathrm{I})$. The elastic constants can also be computed
from DFT by calculating the harmonic phonon spectrum (Fig. 1c) and identifying
the frequency of the phonon mode associated with the ferroelectric
displacements $\omega_{FE}=4.4$ THz. Taking into account the atomic masses of
each species $m_{\alpha}$ the elastic constants are given by

$k_{\alpha}=m_{\alpha}\omega_{FE}^{2}.$ (4)

We obtain $k_{Ni}=146$ meV/Å2, $k_{I}=315$ meV/Å2, which using eq. (3) lead to
$\Delta_{Ni}=0.06$ Å and $\Delta_{I}=0.014$ Å in the opposite direction (Fig.
1a). An effective ferroelectric dipole $p=de^{*}=0.003$ eÅ can be estimated by
considering an average atomic displacement from the computed ones $d=0.03$ Å
and an effective charge $e^{*}=\beta e$, with $e$ the charge of the electron
and $\beta=0.1$ a constant that can be extracted from the Born effective
charges.

The introduction of the ferroelectric distortion breaks the C6 symmetry of the
crystal, leading to a C2 symmetry. As a result, the six first ($J_{1}$) and
third ($J_{3}$) nearest neighbor exchange interactions in the spin model (1)
break their degeneracy resulting in the new sets of {$J_{1}^{\prime}$,
$J_{1}^{\prime\prime}$} and {$J_{3}^{\prime}$, $J_{3}^{\prime\prime}$}
magnetic exchanges respectively as shown in Fig. 1a. We can compute for this
distorted structure the new set of magnetic exchange interactions (see SI for
the values) and get the ground state $\mathbf{q}$ vector that minimizes the
energy using eq. (2). Due to the symmetry breaking, there are only two energy
minima degenerated (Fig. 2b) in the distorted structure in contrast to the six
degenerate minima (Fig. 2a) occurring in the undistorted structure.

Figure 2:  (a) Ground state energy of (1) under variation of the spin-spiral
vector $\mathbf{q}$ for undistorted structure with C6 symmetry; (b) Same as
(a), but for distorted structure with C2 symmetry; (c) Magnon dispersion for
the first Brillouin zone of the undistorted structure; (d) Modification of
magnon energies via the introduction of the ferroelectric order; (e) Magnon
energies for the undistorted (blue, solid) and distorted (red, dashed)
structures. Inset shows the splitting; (f) Magnon density of states (DOS) for
the undistorted (blue, solid) and distorted (red, dashed) structures. Arrows
show positions of the peaks for the first three van Hove singularities; In
plots (a-d) first Brillouin zone is indicated as a red hexagon. In plots (c,d)
ordering vectors of spin-spiral $\mathbf{q}$ are plotted as white or black
arrows, positions of saddle points that correspond to the first three van Hove
singularities in DOS of plot (f) are marked with white or black markers.

###  2.2 Magnons in multiferroic NiI_2

Magnon excitations appear as fluctuations of the spin spiral magnetic order of
NiI2. The magnetic ground state for the undistorted and the ferroelectric
structures allows us to define the spin-wave excitations. Due to the
magnetoelectric coupling of this material, some magnon branches will couple to
the ferroelectric order, giving rise to electromagnon excitations54, 21. The
magnon Hamiltonian is obtained via the expansion of the Heisenberg Hamiltonian
(1) with a Holstein-Primakoff transformation 55 within the linear
approximation 56 and consecutive diagonalization of bosonic Hamiltonian 57.
For each spin of the lattice the quantization axis is chosen along the
direction defined by the spiral configuration and the transversal excitations
are defined in the plane perpendicular to the spin direction. We arrive at the
expression for the magnon Hamiltonian (see SI for details), whose
eigenfrequencies take the form58
$\mathcal{H}^{LSWT}=\sum_{k}\hbar\omega(\boldsymbol{k})\biggl{[}\beta_{k}^{\dagger}\beta_{k}+\dfrac{1}{2}\biggr{]}$,
with
$\hbar\omega(\mathbf{k})=2S\cdot\sqrt{A^{2}(\boldsymbol{k})-B^{2}(\boldsymbol{k})}$,
$A(\boldsymbol{k})=\sum_{j}J_{j}\biggl{(}\bigl{[}1+\cos(\mathbf{q}\cdot\boldsymbol{\delta}_{j})\bigr{]}\cdot\cos(\mathbf{k}\cdot\boldsymbol{\delta}_{j})-2\cos(\mathbf{q}\cdot\boldsymbol{\delta}_{j})\biggr{)}$,
$B(\boldsymbol{k})=\sum_{j}J_{j}\bigl{[}1-\cos(\mathbf{q}\cdot\boldsymbol{\delta}_{j})\bigr{]}\cdot\cos(\mathbf{k}\cdot\boldsymbol{\delta}_{j})$
where $\boldsymbol{\delta}_{j}$ are the vectors connecting first and third
nearest neighbors and index $j$ runs over the set of pairs
{$(J_{1}^{\prime},\boldsymbol{\delta}_{1})$,
$(J_{1}^{\prime\prime},\boldsymbol{\delta}_{2})$,
$(J_{1}^{\prime\prime},\boldsymbol{\delta}_{3})$,
$(J_{3}^{\prime},2\boldsymbol{\delta}_{1})$,
$(J_{3}^{\prime\prime},2\boldsymbol{\delta}_{2})$,
$(J_{3}^{\prime\prime},2\boldsymbol{\delta}_{3})$}.

The magnon dispersion of the undistorted structure does not exhibit C6
symmetry as it depends on the choice between three equivalent spiral vectors:
along $\Gamma$K1, $\Gamma$K and $\Gamma$K2. However, Brillouin zone-resolved
magnon energies for three spiral vectors are equivalent upon rotation by
$0^{\circ}$, $-60^{\circ}$, $-120^{\circ}$ correspondingly (see Fig. S2). In
Fig. 2c we present the magnon dispersion for the undistorted structure with
the $\mathbf{q}$-vector along $\Gamma$K1 direction, for direct comparison with
magnon energies of the structure with emergent polarization. The introduction
of the spiral keeps the C6 symmetry of the points with minimal energies, but
the excitations with $\mathbf{k}$ vectors along the paths $\Gamma$K and
$\Gamma$K1 are not equivalent in general (see Fig. 2e). Inclusion of the
emerging polarization introduces non-uniform shift of the magnon energies
across the whole Brillouin zone (see Fig. 2d). However, the major effect is
the opening of the splitting due to the breaking of C6 rotation symmetry. The
appearance of a splitting in one of the magnon branches in both cases can be
accounted as characteristic signature of the emerging ferroelectric order in
NiI_2.

The spectra of magnonic excitations can be characterized by the magnon density
of states (DOS). There are $6$ van Hove singularities 59, that correspond to
the six groups of energy degenerate saddle points in the magnon dispersion. We
plot only three singularities in the low energy region (Fig. 2f) that
correspond to the saddle points depicted in Figs. 2c,d as only those points
can be observed at low temperatures (see SI for the full energy range).
Interestingly, the peak that corresponds to the saddle point is defined by the
intersections of two lines. First, the line perpendicular to the spin-spiral
vector and second, the line connecting two energy minima where the splitting
is opened upon ferroelectric distortion (stars in Fig. 2). This second line
changes as the distortion gets stronger,. The magnon state, that would
correspond to this point in momentum space, can be describe as a set of one-
dimensional spirals along the $x$ axis with the phase shift of magnetic
excitations between them and each fifth spiral being in phase. Such magnetic
excitations can be probed via inelastic tunneling spectroscopy60, 61, or
electrically driven paramagnetic resonance 62, 63, allowing to directly image
the magnon DOS at the atomic scale64, 65.

Figure 3:  (a) Sketch of the possible device for the application of electric
field, gray (purple) circles depict Ni (I) atoms. (b) Ferroelectric
displacements vs applied electric field. Values that were used in simulations
are plotted with circles (squares) for Ni (I), lines are the square root fits
($\Delta\propto\sqrt{E}$). Zero corresponds to the undistorted structure; (c)
Magnetic exchange parameters vs applied electric fields. Values are relative
to the ones of undisturbed structure; (d) Splitting in the magnon energies at
the rotational image of the spin-spiral vector $\mathbf{q}$ along the
$\Gamma$K path vs electric field. Circles are for the calculated values, line
is a square root fit (Gap $\propto\sqrt{E}$); (e) Magnon DOS as for three
values of applied electric field: $0$ mV/Å (red, solid), $50$ mV/Å (orange,
dashed), $100$ mV/Å (cyan, dash-dotted). Arrows indicate positions of the
first three van Hove singularities; (f) Magnon DOS for the full range of
energies and applied values of electric field. White line separates DOS of the
C2 symmetric (undistorted) structure.

###  2.3 Electrical control of magnons

Finally, we analyze the control of spin-waves using an external electric
field. The emergence of an electric dipole in multiferroic NiI2 allows
coupling magnons to an external electric field $\mathbf{E}$. The conceptual
scheme of the device for the measurement at the monolayer limit was
demonstrated recently 66 and depicted in Fig. 3a. This phenomena can be
effectively modeled in a linear response regime, where the energy of the
electric dipole is balanced by the energy of the harmonic ferroelectric
displacements
$\mathbf{p}\cdot\mathbf{E}=\sum_{\alpha}\dfrac{k_{\alpha}\Delta_{\alpha}^{2}}{2}.$.
Note that only the component of the electric field that is parallel to the
electric polarization couples to it in lowest order. Therefore, we will
consider an external field parallel to the polarization. Considering the
relationship between the ferroelectric displacements of the two species from
the force balance as $\Delta_{Ni}=2k_{I}/k_{Ni}\Delta_{I}$, we can estimate
the contribution to the two ferroelectric displacements introduced by the
external electric field ($\Delta_{Ni}(E)$ and $\Delta_{I}(E)$) as
$\Delta_{I}(E)=\sqrt{\dfrac{k_{Ni}}{(2k_{I}+k_{Ni})k_{I}}Ep}.$. The total
atomic displacements are the sum of the displacements caused by the emergent
ferroelectric dipole and the distortion caused by the external electric field.
Figure 3b shows the ferroelectric displacements as a function of an external
electric field. Effectively, the distortions induced by the electric field
modify the set of exchange parameters
{$J_{1}^{\prime}$,$J_{1}^{\prime\prime}$} and
{$J_{3}^{\prime}$,$J_{3}^{\prime\prime}$}, keeping the symmetry of magnetic
exchange interactions of the distorted structure as shown in Fig. 3c. As
$\sqrt{E}\propto\Delta$ and changes of the exchange couplings are proportional
to $\sqrt{E}$ as well, we conclude that applied values of electric field
induce linear response of the exchange couplings $J\propto\Delta$. With the
new set of exchange parameters we apply the same analysis of the magnetic
excitations as in the previous section. In Fig. 3d the evolution of the
splitting in the magnon energies follows a similar trend as the exchange
couplings, as evident from the equation (1) (see SI for the changes of the
magnon dispersion over the first Brillouin zone). In the same way as for the
splitting caused by the emerging ferroelectric order, the dependency of the
splitting height with the electric field can be probed using dynamical
measurements of the magnetic response21, or magnetic resonance67, 68.

Following the fact that the electric field causes distortion that is aligned
with the distortion caused by the symmetry breaking, the response of the first
van Hove singularity to the applied electric field is more robust. In Figs.
3e-f we plot the magnon DOS for the low energy region. One can see that the
lower energy peak approaches the second and third, while the latter once show
weak sensitivity to the electric field. That behavior opens a way for the non-
invasive selective control of magnon van Hove singularities in magnetic
system. This electric control of magnons should be contrasted with the case of
electronic systems, where the control of van Hove singularities usually
requires twisting of two layers 69, chemical doping 70 or gating 71. The
application of electric field in multiferroic material offers a controllable
and reversible knob to engineer the magnonic dispersion.

##  3 Conclusions

The electrical control of magnon transport is one of milestones in ultrathin
spintronic devices. Here, we demonstrated a mechanism to control the magnon
dispersion in the monolayer multiferroic NiI2 by leveraging its intrinsic
magnetoelectric coupling. Our mechanism exploits that ferroelectric order in
NiI2 strongly couples to the magnon spectra, providing a strategy to modify
spin waves transport with external electric fields. Our strategy leverages the
microscopic mechanism leading to multiferroicity in NiI2, namely the
combination of the spin-spiral order and strong spin-orbit coupling, which
leads to displacements on the atomic positions. We first showed that the
atomic displacements associated to the ferroelectricity lead to a modification
in the magnetic exchange interactions, directly impacting the magnon spectra.
Thanks to the emergence of in-plane ferroelectric order, an in-plane bias
allows to directly modify the structural distortion by directly coupling to
the ferroelectric polarization. We showed that such electric control of the
magnon dispersion enables controlling the low energy dispersion of the magnon
modes. In particular, a single-layer of NiI2 features a magnon branch whose
splitting is electrically controlled by the in-plane bias. Ultimately, our
mechanism establishes a strategy to control low energy magnon transport
electrically, opening a pathway towards designing ultrathin electrically
controllable magnonic devices.

This work was supported by the European Union (ERC-2021-StG-101042680
2D-SMARTiES) and the Excellence Unit “María de Maeztu” (CEX2019-000919-M).
J.J.B. and A.R. acknowledge the Generalitat Valenciana (grant CIDEXG/2023/1
and PhD scholarship GRISOLIAP/2021/038, respectively). A.O.F. and J.L.L.
acknowledges financial support from the Academy of Finland Projects Nos.
331342, 358088, and 349696, the Jane and Aatos Erkko Foundation, and the
Finnish Quantum Flagship. A.R. thanks Jaime Ferrer for fruitful discussions.
We acknowledge the computational resources provided by the Aalto Science-IT
project and HPC systems Cobra and Raven at the Max Planck Computing and Data
Facility.

Detailed calculation scheme and relevant simulation parameters are available
in the Supporting Information.

## References

  * Wang et al. 2022 Wang, Q. H. et al. The Magnetic Genome of Two-Dimensional van der Waals Materials. _ACS Nano_ 2022, _16_ , 6960–7079
  * Gibertini et al. 2019 Gibertini, M.; Koperski, M.; Morpurgo, A. F.; Novoselov, K. S. Magnetic 2D materials and heterostructures. _Nature Nanotechnology_ 2019, _14_ , 408–419
  * Burch et al. 2018 Burch, K. S.; Mandrus, D.; Park, J.-G. Magnetism in two-dimensional van der Waals materials. _Nature_ 2018, _563_ , 47–52
  * Blei et al. 2021 Blei, M.; Lado, J. L.; Song, Q.; Dey, D.; Erten, O.; Pardo, V.; Comin, R.; Tongay, S.; Botana, A. S. Synthesis, engineering, and theory of 2D van der Waals magnets. _Applied Physics Reviews_ 2021, _8_
  * Cenker et al. 2022 Cenker, J.; Sivakumar, S.; Xie, K.; Miller, A.; Thijssen, P.; Liu, Z.; Dismukes, A.; Fonseca, J.; Anderson, E.; Zhu, X.; Roy, X.; Xiao, D.; Chu, J.-H.; Cao, T.; Xu, X. Reversible strain-induced magnetic phase transition in a van der Waals magnet. _Nature Nanotechnology_ 2022, _17_ , 256–261
  * Esteras et al. 2022 Esteras, D. L.; Rybakov, A.; Ruiz, A. M.; Baldoví, J. J. Magnon Straintronics in the 2D van der Waals Ferromagnet CrSBr from First-Principles. _Nano Letters_ 2022, _22_ , 8771–8778
  * Zhang et al. 2022 Zhang, Y.; Kim, H.; Zhang, W.; Watanabe, K.; Taniguchi, T.; Gao, Y.; Maruyama, M.; Okada, S.; Shinokita, K.; Matsuda, K. Magnon-Coupled Intralayer Moiré Trion in Monolayer Semiconductor–Antiferromagnet Heterostructures. _Advanced Materials_ 2022, _34_ , 2200301
  * Jiang et al. 2018 Jiang, S.; Li, L.; Wang, Z.; Mak, K. F.; Shan, J. Controlling magnetism in 2D CrI3 by electrostatic doping. _Nature Nanotechnology_ 2018, _13_ , 549–553
  * Huang et al. 2018 Huang, B.; Clark, G.; Klein, D. R.; MacNeill, D.; Navarro-Moratalla, E.; Seyler, K. L.; Wilson, N.; McGuire, M. A.; Cobden, D. H.; Xiao, D.; Yao, W.; Jarillo-Herrero, P.; Xu, X. Electrical control of 2D magnetism in bilayer CrI3. _Nature Nanotechnology_ 2018, _13_ , 544–548
  * Zhang et al. 2020 Zhang, X.-X.; Li, L.; Weber, D.; Goldberger, J.; Mak, K. F.; Shan, J. Gate-tunable spin waves in antiferromagnetic atomic bilayers. _Nature Materials_ 2020, _19_ , 838–842
  * Mishra et al. 2024 Mishra, S.; Park, I. K.; Javaid, S.; Shin, S. H.; Lee, G. Enhancement of interlayer exchange coupling via intercalation in 2D magnetic bilayers: towards high Curie temperature. _Mater. Horiz._ 2024, –
  * Iturriaga et al. 2023 Iturriaga, H.; Martinez, L. M.; Mai, T. T.; Biacchi, A. J.; Augustin, M.; Hight Walker, A. R.; Sanad, M. F.; Sreenivasan, S. T.; Liu, Y.; Santos, E. J., et al. Magnetic properties of intercalated quasi-2D Fe3-xGeTe2 van der Waals magnet. _npj 2D Materials and Applications_ 2023, _7_ , 56
  * Ruiz et al. 2024 Ruiz, A. M.; Rivero-Carracedo, G.; Rybakov, A.; Dey, S.; Baldoví, J. J. Towards molecular controlled magnonics. _Nanoscale Adv._ 2024, _6_ , 3320–3328
  * Huang et al. 2017 Huang, B.; Clark, G.; Navarro-Moratalla, E.; Klein, D. R.; Cheng, R.; Seyler, K. L.; Zhong, D.; Schmidgall, E.; McGuire, M. A.; Cobden, D. H.; Yao, W.; Xiao, D.; Jarillo-Herrero, P.; Xu, X. Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit. _Nature_ 2017, _546_ , 270–273
  * Bedoya-Pinto et al. 2021 Bedoya-Pinto, A.; Ji, J.-R.; Pandeya, A. K.; Gargiani, P.; Valvidares, M.; Sessi, P.; Taylor, J. M.; Radu, F.; Chang, K.; Parkin, S. S. P. Intrinsic 2D-XY ferromagnetism in a van der Waals monolayer. _Science_ 2021, _374_ , 616–620
  * Ruan et al. 2021 Ruan, W. et al. Evidence for quantum spin liquid behaviour in single-layer 1T-TaSe2 from scanning tunnelling microscopy. _Nature Physics_ 2021, _17_ , 1154–1161
  * Vaňo et al. 2021 Vaňo, V.; Amini, M.; Ganguli, S. C.; Chen, G.; Lado, J. L.; Kezilebieke, S.; Liljeroth, P. Artificial heavy fermions in a van der Waals heterostructure. _Nature_ 2021, _599_ , 582–586
  * Posey et al. 2024 Posey, V. A. et al. Two-dimensional heavy fermions in the van der Waals metal CeSiI. _Nature_ 2024, _625_ , 483–488
  * Song et al. 2022 Song, Q.; Occhialini, C. A.; Ergeçen, E.; Ilyas, B.; Amoroso, D.; Barone, P.; Kapeghian, J.; Watanabe, K.; Taniguchi, T.; Botana, A. S.; Picozzi, S.; Gedik, N.; Comin, R. Evidence for a single-layer van der Waals multiferroic. _Nature_ 2022, _602_ , 601–605
  * Amini et al. 2024 Amini, M.; Fumega, A. O.; González‐Herrero, H.; Vaňo, V.; Kezilebieke, S.; Lado, J. L.; Liljeroth, P. Atomic‐Scale Visualization of Multiferroicity in Monolayer NiI2. _Advanced Materials_ 2024, _36_
  * Gao et al. 2024 Gao, F. Y. et al. Giant chiral magnetoelectric oscillations in a van der Waals multiferroic. _Nature_ 2024,
  * Barman et al. 2021 Barman, A. et al. The 2021 Magnonics Roadmap. _Journal of Physics: Condensed Matter_ 2021, _33_ , 413001
  * Flebus et al. 2024 Flebus, B. et al. The 2024 magnonics roadmap. _Journal of Physics: Condensed Matter_ 2024, _36_ , 363501
  * Kimura et al. 2003 Kimura, T.; Goto, T.; Shintani, H.; Ishizaka, K.; Arima, T.; Tokura, Y. Magnetic control of ferroelectric polarization. _Nature_ 2003, _426_ , 55–58
  * Hur et al. 2004 Hur, N.; Park, S.; Sharma, P. A.; Ahn, J. S.; Guha, S.; Cheong, S.-W. Electric polarization reversal and memory in a multiferroic material induced by magnetic fields. _Nature_ 2004, _429_ , 392–395
  * Fiebig 2005 Fiebig, M. Revival of the magnetoelectric effect. _Journal of Physics D: Applied Physics_ 2005, _38_ , R123–R152
  * Gajek et al. 2007 Gajek, M.; Bibes, M.; Fusil, S.; Bouzehouane, K.; Fontcuberta, J.; Barthélémy, A.; Fert, A. Tunnel junctions with multiferroic barriers. _Nature Materials_ 2007, _6_ , 296–302
  * Pantel et al. 2012 Pantel, D.; Goetze, S.; Hesse, D.; Alexe, M. Reversible electrical switching of spin polarization in multiferroic tunnel junctions. _Nature Materials_ 2012, _11_ , 289–293
  * Hu et al. 2015 Hu, J.; Chen, L.; Nan, C. Multiferroic Heterostructures Integrating Ferroelectric and Magnetic Materials. _Advanced Materials_ 2015, _28_ , 15–39
  * Matsubara et al. 2015 Matsubara, M.; Manz, S.; Mochizuki, M.; Kubacka, T.; Iyama, A.; Aliouane, N.; Kimura, T.; Johnson, S. L.; Meier, D.; Fiebig, M. Magnetoelectric domain control in multiferroic TbMnO 3. _Science_ 2015, _348_ , 1112–1115
  * Fumega and Lado 2022 Fumega, A. O.; Lado, J. L. Microscopic origin of multiferroic order in monolayer NiI2. _2D Materials_ 2022, _9_ , 025010
  * Fumega and Lado 2023 Fumega, A. O.; Lado, J. L. Moiré-driven multiferroic order in twisted CrCl3, CrBr3 and CrI3 bilayers. _2D Materials_ 2023, _10_ , 025026
  * Billerey et al. 1977 Billerey, D.; Terrier, C.; Ciret, N.; Kleinclauss, J. Neutron diffraction study and specific heat of antiferromagnetic NiI2. _Physics Letters A_ 1977, _61_ , 138–140
  * Kuindersma et al. 1981 Kuindersma, S.; Sanchez, J.; Haas, C. Magnetic and structural investigations on NiI2 and CoI2. _Physica B+C_ 1981, _111_ , 231–248
  * Kurumaji et al. 2013 Kurumaji, T.; Seki, S.; Ishiwata, S.; Murakawa, H.; Kaneko, Y.; Tokura, Y. Magnetoelectric responses induced by domain rearrangement and spin structural change in triangular-lattice helimagnets NiI2 and CoI2. _Physical Review B_ 2013, _87_
  * Liu et al. 2020 Liu, H.; Wang, X.; Wu, J.; Chen, Y.; Wan, J.; Wen, R.; Yang, J.; Liu, Y.; Song, Z.; Xie, L. Vapor Deposition of Magnetic Van der Waals NiI2 Crystals. _ACS Nano_ 2020, _14_ , 10544–10551
  * Ju et al. 2021 Ju, H.; Lee, Y.; Kim, K.-T.; Choi, I. H.; Roh, C. J.; Son, S.; Park, P.; Kim, J. H.; Jung, T. S.; Kim, J. H.; Kim, K. H.; Park, J.-G.; Lee, J. S. Possible Persistence of Multiferroic Order down to Bilayer Limit of van der Waals Material NiI2. _Nano Letters_ 2021, _21_ , 5126–5132
  * Geim and Grigorieva 2013 Geim, A. K.; Grigorieva, I. V. Van der Waals heterostructures. _Nature_ 2013, _499_ , 419–425
  * Li et al. 2023 Li, X.; Xu, C.; Liu, B.; Li, X.; Bellaiche, L.; Xiang, H. Realistic Spin Model for Multiferroic NiI2. _Physical Review Letters_ 2023, _131_
  * Kapeghian et al. 2024 Kapeghian, J.; Amoroso, D.; Occhialini, C. A.; Martins, L. G. P.; Song, Q.; Smith, J. S.; Sanchez, J. J.; Kong, J.; Comin, R.; Barone, P.; Dupé, B.; Verstraete, M. J.; Botana, A. S. Effects of pressure on the electronic and magnetic properties of bulk NiI2. _Physical Review B_ 2024, _109_
  * Sødequist and Olsen 2023 Sødequist, J.; Olsen, T. Type II multiferroic order in two-dimensional transition metal halides from first principles spin-spiral calculations. _2D Materials_ 2023, _10_ , 035016
  * Mostovoy 2006 Mostovoy, M. Ferroelectricity in Spiral Magnets. _Physical Review Letters_ 2006, _96_
  * Giannozzi et al. 2017 Giannozzi, P. et al. Advanced capabilities for materials modelling with Quantum ESPRESSO. _Journal of Physics: Condensed Matter_ 2017, _29_ , 465901
  * Giannozzi et al. 2009 Giannozzi, P. et al. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. _Journal of Physics: Condensed Matter_ 2009, _21_ , 395502
  * Dudarev et al. 1998 Dudarev, S. L.; Botton, G. A.; Savrasov, S. Y.; Humphreys, C. J.; Sutton, A. P. Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study. _Physical Review B_ 1998, _57_ , 1505–1509
  * Ernzerhof and Scuseria 1999 Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew–Burke–Ernzerhof exchange-correlation functional. _The Journal of Chemical Physics_ 1999, _110_ , 5029–5036
  * Hamann 2013 Hamann, D. R. Optimized norm-conserving Vanderbilt pseudopotentials. _Physical Review B_ 2013, _88_
  * Monkhorst and Pack 1976 Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. _Physical Review B_ 1976, _13_ , 5188–5192
  * Timrov et al. 2018 Timrov, I.; Marzari, N.; Cococcioni, M. Hubbard parameters from density-functional perturbation theory. _Physical Review B_ 2018, _98_
  * Timrov et al. 2022 Timrov, I.; Marzari, N.; Cococcioni, M. HP – A code for the calculation of Hubbard parameters using density-functional perturbation theory. _Computer Physics Communications_ 2022, _279_ , 108455
  * Pizzi et al. 2020 Pizzi, G. et al. Wannier90 as a community code: new features and applications. _Journal of Physics: Condensed Matter_ 2020, _32_ , 165902
  * He et al. 2021 He, X.; Helbig, N.; Verstraete, M. J.; Bousquet, E. TB2J: A python package for computing magnetic interaction parameters. _Computer Physics Communications_ 2021, _264_ , 107938
  * Solovyev 2021 Solovyev, I. V. Exchange interactions and magnetic force theorem. _Physical Review B_ 2021, _103_
  * Kim et al. 2023 Kim, J. H.; Jung, T. S.; Lee, Y.; Kim, C.; Park, J.-G.; Kim, J. H. Terahertz evidence of electromagnon excitations in the multiferroic van der Waals insulator $\mathrm{Ni}{\mathrm{I}}_{2}$. _Phys. Rev. B_ 2023, _108_ , 064414
  * Holstein and Primakoff 1940 Holstein, T.; Primakoff, H. Field Dependence of the Intrinsic Domain Magnetization of a Ferromagnet. _Physical Review_ 1940, _58_ , 1098–1113
  * Toth and Lake 2015 Toth, S.; Lake, B. Linear spin wave theory for single-Q incommensurate magnetic structures. _Journal of Physics: Condensed Matter_ 2015, _27_ , 166002
  * White et al. 1965 White, R. M.; Sparks, M.; Ortenburger, I. Diagonalization of the Antiferromagnetic Magnon-Phonon Interaction. _Physical Review_ 1965, _139_ , A450–A454
  * Cong and Shen 2024 Cong, A.; Shen, K. Soft magnons in van der Waals multiferroic NiI. _Physical Review B_ 2024, _109_
  * Van Hove 1953 Van Hove, L. The Occurrence of Singularities in the Elastic Frequency Distribution of a Crystal. _Physical Review_ 1953, _89_ , 1189–1193
  * Heinrich et al. 2004 Heinrich, A. J.; Gupta, J. A.; Lutz, C. P.; Eigler, D. M. Single-Atom Spin-Flip Spectroscopy. _Science_ 2004, _306_ , 466–469
  * Klein et al. 2018 Klein, D. R.; MacNeill, D.; Lado, J. L.; Soriano, D.; Navarro-Moratalla, E.; Watanabe, K.; Taniguchi, T.; Manni, S.; Canfield, P.; Fernández-Rossier, J.; Jarillo-Herrero, P. Probing magnetism in 2D van der Waals crystalline insulators via electron tunneling. _Science_ 2018, _360_ , 1218–1222
  * Baumann et al. 2015 Baumann, S.; Paul, W.; Choi, T.; Lutz, C. P.; Ardavan, A.; Heinrich, A. J. Electron paramagnetic resonance of individual atoms on a surface. _Science_ 2015, _350_ , 417–420
  * Wang et al. 2024 Wang, H.; Fan, P.; Chen, J.; Jiang, L.; Gao, H.-J.; Lado, J. L.; Yang, K. Realizing topological quantum magnets with atomic spins on surfaces. _arXiv e-prints_ 2024, arXiv:2403.14145
  * Spinelli et al. 2014 Spinelli, A.; Bryant, B.; Delgado, F.; Fernández-Rossier, J.; Otte, A. F. Imaging of spin waves in atomically designed nanomagnets. _Nature Materials_ 2014, _13_ , 782–785
  * Ganguli et al. 2023 Ganguli, S. C.; Aapro, M.; Kezilebieke, S.; Amini, M.; Lado, J. L.; Liljeroth, P. Visualization of Moiré Magnons in Monolayer Ferromagnet. _Nano Letters_ 2023, _23_ , 3412–3417
  * Higashitarumizu et al. 2020 Higashitarumizu, N.; Kawamoto, H.; Lee, C.-J.; Lin, B.-H.; Chu, F.-H.; Yonemori, I.; Nishimura, T.; Wakabayashi, K.; Chang, W.-H.; Nagashio, K. Purely in-plane ferroelectricity in monolayer SnS at room temperature. _Nature Communications_ 2020, _11_
  * Qin et al. 2021 Qin, H.; Holländer, R. B.; Flajšman, L.; Hermann, F.; Dreyer, R.; Woltersdorf, G.; van Dijken, S. Nanoscale magnonic Fabry-Pérot resonator for low-loss spin-wave manipulation. _Nature Communications_ 2021, _12_
  * Colombano et al. 2020 Colombano, M. F.; Arregui, G.; Bonell, F.; Capuj, N. E.; Chavez-Angel, E.; Pitanti, A.; Valenzuela, S. O.; Sotomayor-Torres, C. M.; Navarro-Urrios, D.; Costache, M. V. Ferromagnetic Resonance Assisted Optomechanical Magnetometer. _Phys. Rev. Lett._ 2020, _125_ , 147201
  * Li et al. 2009 Li, G.; Luican, A.; Lopes dos Santos, J. M. B.; Castro Neto, A. H.; Reina, A.; Kong, J.; Andrei, E. Y. Observation of Van Hove singularities in twisted graphene layers. _Nature Physics_ 2009, _6_ , 109–113
  * Markiewicz 1990 Markiewicz, R. S. Phase separation near the Mott transition in La2-xSrxCuO4. _Journal of Physics: Condensed Matter_ 1990, _2_ , 665–676
  * Novoselov et al. 2004 Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. _Science_ 2004, _306_ , 666–669


#  Quid Manumit - Freeing the Qubit for Art
Embedded Quantum Simulators for Musical Instruments

Mark Carney1 0000-0001-9372-9033

###### Abstract

This paper describes how to ‘Free the Qubit’ for art, by creating standalone
quantum musical effects and instruments. Previously released quantum simulator
code for an ARM-based Raspberry Pi Pico embedded microcontroller is utilised
here, and several examples are built demonstrating different methods of
utilising embedded resources: The first is a Quantum MIDI processor that
generates additional notes for accompaniment and unique quantum generated
instruments based on the input notes, decoded and passed through a quantum
circuit in an embedded simulator. The second is a Quantum Distortion module
that changes an instrument’s raw sound according to a quantum circuit, which
is presented in two forms; a self-contained Quantum Stylophone, and an effect
module plugin called ’QubitCrusher’ for the Korg Nu:Tekt NTS-1. This paper
also discusses future work and directions for quantum instruments, and
provides all examples as open source. This is, to the author’s knowledge, the
first example of embedded Quantum Simulators for Instruments of Music (another
QSIM).

##  1 Introduction

Quantum Music is a new and inspiring take on the possible futures of digitized
and synthesized sound. Previous work, much of which can be found in
[Miranda2022].

With the advent of free and pay-as-you-go Quantum Compute-as-a-Service
(QCaaS), such as services offered by IBM and on the AWS Braket platform
respectively, it has never been easier nor more accessible to perform quantum
computing tasks on real quantum hardware.

We have also seen an explosion in new simulators, possibly the most impressive
being a 256-qubit quantum simulator developed by Ebadi _et al._ in
[Ebadi2021].

Whilst we await the promise of the coming power offered by real quantum
hardware, high-quality quantum simulators have been indispensable to both
students and professionals alike in learning the inner workings of Quantum
Information Theory, such as found in now classic textbooks like [Nielsen2012].

Previous work from Quantum Village has already dealt with ways in which
quantum circuits can be sonified directly - producing the ability to ‘hear’
entanglement and superposition by means of blending tones and overtones based
on the statevector of a system. This project, called Soniq, is freely
available FOSS, and can be found in [Carney_Soniq].

Although there have been many investigations and research projects leading to
the creation of many quantum synthesizers [Miranda2022], the aim of this
project was to create self-contained standalone devices that can show how to
‘add quantum’ in ways that become natural parts of the audio processing chains
in use today.

To this end, the investigations for this paper focused on creating small
effects units and devices that interact with common audio communications
protocols, the best example being MIDI. Along the way, a number of methods
were found to embed quantum circuit simulations into existing devices - both
by the physical addition of embedded quantum simulators, but also quantum-
based plugins for widely available synthesizers.

Putz and Svozil in [Reck2022] consider whether quantum art will ever reveal
something that classical (in many senses) art could not? Does quantum
expressibility, with access to exclusively quantum notions such as
entanglement and superposition, enhance the musical experience beyond that
which classical or contemporary music has been able to so far? This paper does
not purport to have an answer, but instead to have been driven to create tools
with which musicians may begin to explore such notions.

##  2 Preliminaries

We introduce the various players in the development of our quantum sonic
modules; the raspberry pi pico, quantum circuits and their simulation, digital
signal processing, and MIDI.

###  2.1 Raspberry Pi Pico

Released in January 2021 [raspberrypi], the Raspberry Pi RP2040, a custom-made
ARM-based microcontroller has fast become a tool of choice for many and varied
applications.

The features specifically useful to our purposes include its wide
availability, low price, powerful ARM Cortex-M0+ dual core processor running
at 133MHz, with 264KB RAM and many peripherals, including 8 programmable PIO
pins, 30 GPIOs and built in integrator and integer divider peripherals.

Many of these features will be used in these projects, especially the ability
to run two separate fast ARM cores separately - one of which will handle
signals acquisition and processing, whilst the other runs our bespoke quantum
simulator code.

###  2.2 The Embedded Quantum Simulator

Figure 1: An assembled and running microcontroller based quantum simulator,
here showing a CHSH game quantum circuit.


# H2-optimal model order reduction over a finite time interval

Kasturi Das  EEE Department, IIT Guwahati  S. Krishnaswamy  EEE Department,
IIT Guwahati  S. Majhi  EEE Department, IIT Guwahati

###### Abstract

For a time-limited version of the H2 norm defined over a fixed time interval,
we obtain a closed form expression of the gradients. After that, we use the
gradients to propose a time-limited model order reduction method. The method
involves obtaining a reduced model which minimizes the time-limited H2 norm,
formulated as a nonlinear optimization problem. The optimization problem is
solved using standard optimization software.

##  1 Introduction

Capturing system dynamics accurately requires large- scale, linear dynamical
models. Simulating or analysing such models and designing controllers require
considerable computational effort. Such issues are resolved by replacing the
large model with a lower order approximation based on various performance
measures. Model order reduction techniques are used in a wide range of areas
including computational aerodynamics, large-scale network systems,
microelectronics, electromagnetic systems, chemical processes etc. [2].
References [3], [1] contain a comprehensive discussion on a large number of
model reduction techniques available in literature.

The H2 norm of the error between the original and the reduced system acts as
an important performance measure for obtaining reduced order models. In H2
optimal model reduction the aim is to find a lower order model that minimizes
this norm. Since finding global minimizers is a difficult task, the existing
methods focus on finding local minimizers. These methods are divided into two
categories: optimization-based methods and tangential interpolation methods.
In optimization based methods, the task of model reduction is formulated as an
optimization problem over various manifolds [22, 24, 18, 19]. The solution
yields optimal reduced models. Tangential interpolation methods use Krylov-
based algorithms and work well for large-scale systems. Examples of such
methods include Iterative Rational Krylov Algorithm (IRKA) [10] and Two-Sided
Iteration Algorithm (TSIA) [23]. Reference [15] deals with optimization based
frequency-limited H2 optimal model reduction.

Availability of simulation data for a finite time interval or the need to
approximate the system behaviour over a finite time interval led to the
development of finite time model reduction methods. These include methods such
as Proper Orthogonal Decomposition (POD) [11], Time-Limited Balanced
Truncation (TL-BT) [8] etc. Error bounds for TL-BT are proposed in [17, 16].
References [12] and [7] deal with implementation of TL-BT for large-scale
continuous and discrete systems respectively. Lyapunov based time-limited H2
optimality conditions are obtained in [9] using a time-limited H2 norm. The
same paper proposes an iterative scheme similar to TSIA [23]. We refer to this
scheme as TL-TSIA. Projection-based algorithms like TL-BT and TL-TSIA fail to
exactly satisfy the time-limited optimality conditions. Reference [20] obtains
interpolation based first-order necessary conditions for time-limited H2
optimality and proposes an optimization algorithm named FHIRKA. This algorithm
produces time-limited H2 optimal models but is valid for SISO systems.

To the best of the authors’ knowledge, there are no optimization based
algorithms in the literature that yield time-limited H2 optimal reduced models
for both SISO and MIMO systems. In this letter, we aim to fill this gap. The
time-limited H2 optimal model reduction problem is formulated as an
optimization problem. We derive closed-form expressions of the gradients of
the objective function. These gradients are used with standard quasi-Newton
solvers to propose a time-limited H2 optimal model reduction method. We
initialize the proposed method with reduced models obtained from time-limited
projection based model reduction techniques. Two numerical examples show how
our proposed method significantly improves the objective function compared to
the projection based methods for reduced orders less than a certain upper
bound. Due to space constraints, we have demonstrated a single example.
However, we observe that the bound is different for different models.

The letter is arranged as follows. In Section 2 we discuss some basic concepts
related to model order reduction over a limited time. In Section 3, we
formulate time-limited H2 optimal model reduction as an optimization problem
and derive the gradients of the objection function in Section 3.1. In Section
3.2 we propose a method for solving the optimization problem by using the
gradients and discussing its computational complexity in Section 3.3. The
proposed method is implemented on two numerical examples in Section 4. We
conclude the paper in Section 5.

#### Notations

Let $\mathbb{R}$ and $\mathbb{C}$ be the set of real and complex numbers
respectively. For a matrix $P\in\mathbb{R}^{n\times n}$, Tr(P) denotes the
trace, $P^{T}$ denotes the transpose, $\left\|P\right\|$ denotes the 2-norm
and $\left\|P\right\|_{F}$ denotes the Frobenius norm of the matrix $P$. Let
us consider a function $f:[0,\infty)\to\mathbb{R}^{p\times m}$ whose Laplace
transform $F(s)\in\mathbb{C}^{p\times m}$ exists. The time-limited H2 norm of
$F$, denoted by $\left\|F\right\|_{\text{H}_{2,\tau}}$, is defined as
$\sqrt{\int_{0}^{\tau}\left\|f(t)\right\|_{F}^{2}dt}$ where
$\left\|f(t)\right\|_{F}^{2}=Tr({f(t)}^{T}f(t))$. The $L_{\infty}^{\tau}$ norm
and $L_{2}^{\tau}$ norm of $f(t)$ are defined as follows,

$\displaystyle\left\|f\right\|_{L_{2}^{\tau}}$
$\displaystyle=\sqrt{\int_{0}^{\tau}\left\|f(t)\right\|^{2}dt}$
$\displaystyle\left\|f\right\|_{L_{\infty}^{\tau}}$
$\displaystyle=\operatorname*{sup}_{t\in[0,\tau]}\left\|f(t)\right\|$

##  2 Preliminaries

A stable and strictly proper Linear Time-Invariant (LTI) system, $\Sigma$ is
given by,

$\dot{x}(t)=Ax(t)+Bu(t),\,x(0)=0,\\\ $ (1a) $y(t)=Cx(t),\,t\geq 0$ (1b)

where $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$ and
$C\in\mathbb{R}^{p\times n}$. Let $H(s)$ be the transfer function and $h(t)$
be the impulse response. We assume that the state dimension $n$ is large and
is much larger than the number of inputs and outputs, i.e. $n\gg m,p$.

Over a limited time interval $\begin{bmatrix}0&\tau\end{bmatrix}$ with
$\tau<\infty$, time-limited gramians are defined in [8] as follows,

$\displaystyle P_{\tau}=\int_{0}^{\tau}e^{At}BB^{T}e^{A^{T}t}dt,\quad
Q_{\tau}=\int_{0}^{\tau}e^{A^{T}t}C^{T}Ce^{At}dt$ (2)

The time-limited gramians are solutions of the following Lyapunov equations

$\displaystyle AP_{\tau}+P_{\tau}A^{T}+BB^{T}-e^{A\tau}BB^{T}e^{{A^{T}}\tau}$
$\displaystyle=0$ (3) $\displaystyle
A^{T}Q_{\tau}+Q_{\tau}A+C^{T}C-e^{A^{T}\tau}C^{T}Ce^{A\tau}$ $\displaystyle=0$
(4)

The following expressions for $\|G\|^{2}_{\text{H}_{2,\tau}}$ are derived in
[9].

$\|G\|^{2}_{\text{H}_{2,\tau}}=\text{Tr }CP_{\tau}C^{T}=\text{Tr
}B^{T}Q_{\tau}B$ (5)

Since the H2,τ norm is defined over a strictly finite time interval, a system
need not be asymptotically stable inorder to have a finite H2,τ norm.

##  3 Optimization based time-limited model reduction

Consider a reduced order system $\hat{\Sigma}$ given by,

$\dot{x_{r}}(t)=A_{r}x_{r}(t)+B_{r}u(t),\,x_{r}(0)=0,\\\ $ (6a)
$y_{r}(t)=C_{r}x_{r}(t),\,t\geq 0$ (6b)

where $A_{r}\in\mathbb{R}^{r\times r}$, $B_{r}\in\mathbb{R}^{r\times m}$ and
$C\in\mathbb{R}^{p\times r}$. Let $H_{r}(s)$ be the transfer function and
$h_{r}(t)$ be the impulse response. It is essential that $r\ll n$ and
$y-y_{r}$ is small for an appropriate time limited norm.

For all admissible inputs $u(t)$ with unity $L_{2}^{\tau}$ norm, the following
relation holds [9].

$\left\|y(t)-y_{r}(t)\right\|_{L_{\infty}^{\tau}}\leq\left\|G-G_{r}\right\|_{\text{H}_{2,\tau}}$
(7)

Thus, minimizing the H2,τ error norm ensures that $y_{r}(t)$ is a good
approximation of $y(t)$ over the time interval
$\begin{bmatrix}0&\tau\end{bmatrix}$.

In this paper, we aim to obtain H2,τ optimal reduced models by solving the
following optimization problem:

$\left\|G-G_{r}\right\|_{\text{H}_{2,\tau}}=\operatorname*{minimize\quad}_{\text{dim}(\hat{G})=r}\left\|G-\hat{G}\right\|_{\text{H}_{2,\tau}}$
(8)

The feasible set for the optimization problem formulated above comprises of
all the reduced order systems of the form (6) with state dimension $r$.
$\left\|G-G_{r}\right\|_{\text{H}_{2,\tau}}$ is defined as
$\left(\bigintss_{0}^{\tau}\left\|Ce^{At}B-C_{r}e^{A_{r}t}B_{r}\right\|_{F}^{2}dt\right)^{\frac{1}{2}}$.
The error system $(G-G_{r})$ can be represented by the following state-space
realization.

$\\{A_{e},B_{e},C_{e}\\}=\left\\{\begin{bmatrix}A&0\\\
0&A_{r}\end{bmatrix},\begin{bmatrix}B\\\
B_{r}\end{bmatrix},\begin{bmatrix}C&-C_{r}\end{bmatrix}\right\\}$ (9)

As a consequence of (5), the square of the H2,τ norm of the above realization
can be expressed as

$\left\|G-G_{r}\right\|_{\text{H}_{2,\tau}}^{2}=\text{Tr}({B_{e}}^{T}Q_{e,\tau}B_{e})=\text{Tr}({C_{e}}P_{e,\tau}{C_{e}}^{T})$
(10)

Here, $P_{e,\tau}$ and $Q_{e,\tau}$ are the time-limited controllability and
observability gramians and they satisfy the following Lyapunov equations,

$A_{e}P_{e,\tau}+P_{e,\tau}{A_{e}}^{T}+B_{e}{B_{e}}^{T}-e^{A_{e}\tau}B_{e}{B_{e}}^{T}e^{{A_{e}}^{T}\tau}=0$
(11a)
${A_{e}}^{T}Q_{e,\tau}+Q_{e,\tau}{A_{e}}+{C_{e}}^{T}C_{e}-e^{{A_{e}}^{T}\tau}{C_{e}}^{T}C_{e}e^{A_{e}\tau}=0$
(11b)

For the realization (9), the corresponding gramians $P_{e,\tau}$ and
$Q_{e,\tau}$ can be partitioned as follows,

$P_{e,\tau}=\begin{pmatrix}P_{\tau}&X_{\tau}\\\
X_{\tau}^{T}&P_{r,\tau}\end{pmatrix}\quad
Q_{e,\tau}=\begin{pmatrix}Q_{\tau}&Y_{\tau}\\\
Y_{\tau}^{T}&Q_{r,\tau}\end{pmatrix}$ (12)

Further the matrix $e^{A_{e}\tau}$ can be partitioned as follows,

$e^{A_{e}\tau}=\begin{bmatrix}e^{A\tau}&0\\\ 0&e^{A_{r}\tau}\end{bmatrix}$
(13)

Substituting the partitions (12) and (13) for $P_{e,\tau}$, $Q_{e,\tau}$ and
$e^{A_{e}\tau}$ in equations (11a) and (11b) we get the following time-limited
Lyapunov and time-limited Sylvester equations,

$AP_{\tau}+P_{\tau}A^{T}+BB^{T}-e^{A\tau}BB^{T}e^{A^{T}\tau}=0$ (14a)
$AX_{\tau}+X_{\tau}{A_{r}}^{T}+B{B_{r}}^{T}-e^{A\tau}B{B_{r}}^{T}e^{{A_{r}}^{T}\tau}=0$
(14b)
$A_{r}P_{r,\tau}+P_{r,\tau}{A_{r}}^{T}+B_{r}{B_{r}}^{T}-e^{{A_{r}}\tau}{B_{r}}{B_{r}}^{T}e^{{A_{r}}^{T}\tau}=0$
(14c) $A^{T}Q_{\tau}+Q_{\tau}A+C^{T}C-e^{A^{T}\tau}C^{T}Ce^{A\tau}=0$ (14d)
$A^{T}Y_{\tau}+Y_{\tau}A_{r}-C^{T}C_{r}+e^{A^{T}\tau}C^{T}C_{r}e^{A_{r}\tau}=0$
(14e)
${A_{r}}^{T}Q_{r,\tau}+Q_{r,\tau}A_{r}+{C_{r}}^{T}{C_{r}}-e^{A_{r}^{T}\tau}{C_{r}}^{T}{C_{r}}e^{{A_{r}}\tau}=0$
(14f)

Here, $P_{\tau}$ and $Q_{\tau}$ are the controllability and observability
gramian respectively for the full order model (1). $P_{r,\tau}$ and
$Q_{r,\tau}$ are the controllability and observability gramian respectively
for the reduced model (6) Additionally, substituting the above partitions we
can simplify (10) as follows

$\displaystyle\left\|G-G_{r}\right\|_{\text{H}_{2,\tau}}^{2}$ (15a)
$\displaystyle=\text{Tr}\left(CP_{\tau}C^{T}-2CX_{\tau}C_{r}^{T}+C_{r}P_{r,\tau}C_{r}^{T}\right)$
$=\text{Tr}\left(B^{T}Q_{\tau}B+2B^{T}Y_{\tau}B_{r}+B_{r}^{T}Q_{r,\tau}B_{r}\right)$
(15b)

###  3.1 Gradients of the Cost Function

For a matrix valued function $f:\mathbb{R}^{m\times n}\to\mathbb{R}$, the
gradient at $M\in\mathbb{R}^{m\times n}$ is another matrix
${\bigtriangledown}f(M)\in\mathbb{R}^{m\times n}$ which is given by Definition
3.1 of [21]. The inner product of two matrices is given by $\langle
A,B\rangle=\text{Tr}(A^{T}B)$. We now proceed to derive gradients of the
objective function (15a or 15b) with respect to the reduced system matrices.
The following lemma from [21] is essential for proving the subsequent theorem.

######  Lemma 3.1.

If $AM+MB+C=0$ and $NA+BN+D=0$ then $\text{Tr}(CN)=\text{Tr}(DM)$.

######  Theorem 3.2.

For the cost function $J=\left\|G-G_{r}\right\|_{H_{2,\tau}}^{2}$, the
gradients with respect to $A_{r}$, $B_{r}$ and $C_{r}$ denoted by
${\nabla}_{A_{r}}J$, ${\nabla}_{B_{r}}J$ and ${\nabla}_{C_{r}}J$ respectively
are

${\nabla}_{A_{r}}J=2(Q_{r,\tau}P_{r}+Y^{T}_{\tau}X+\tau(L(A_{r}\tau,S_{\tau})^{T}))$
(16a) ${\nabla}_{B_{r}}J=2(Q_{r,\tau}B_{r}+Y^{T}_{\tau}B)$ (16b)
${\nabla}_{C_{r}}J=2(C_{r}P_{r,\tau}-CX_{\tau})$ (16c)

where $P_{r,\tau}$, $Q_{r,\tau}$, $X_{\tau}$, and $Y_{\tau}$ are solutions of
(14c), (14f), (14b), and (14e) respectively. $P_{r}$ and $X$ are obtained by
solving the following Lyapunov and Sylvester equations.

$P_{r}A_{r}^{T}+A_{r}P_{r}+B_{r}B_{r}^{T}=0$ (17a)
$X^{T}A^{T}+A_{r}X^{T}+B_{r}B^{T}=0$ (17b)

Here, the function $L(X,Y)$ is the Fréchet derivative of the matrix
exponential of $X$ along the direction $Y$ [4]. $S_{\tau}$ is given by

$S_{\tau}=\left(X^{T}e^{A^{T}\tau}C^{T}C_{r}-P_{r}e^{A_{r}^{T}\tau}C_{r}^{T}C_{r}\right)$
(18)

###### Proof.

Consider the expression (15b) of the cost function. For a perturbation of
$\Delta_{A_{r}}$ in $A_{r}$, the corresponding first-order perturbation in $J$
denoted by $\Delta_{J}^{A_{r}}$ is,

$\displaystyle\Delta_{J}^{A_{r}}$
$\displaystyle=\text{Tr}\left(2B_{r}B^{T}\Delta_{Y_{\tau}}+B_{r}B_{r}^{T}\Delta_{Q_{r,\tau}}\right)$
(19)

$\Delta_{Y_{\tau}}$, $\Delta_{Q_{r,\tau}}$, $\Delta_{e^{A_{r}\tau}}$ are the
perturbations in $Y_{\tau}$, $Q_{r,\tau}$ and $e^{A_{r}\tau}$ respectively due
to the perturbation $\Delta_{A_{r}}$ in $A_{r}$. The relation between the
perturbations $\Delta_{Y_{\tau}}$ and $\Delta_{A_{r}}$ is through equation
(14e).

$A^{T}\Delta_{Y_{\tau}}+\Delta_{Y_{\tau}}A_{r}+Y_{\tau}\Delta_{A_{r}}+e^{A^{T}\tau}C^{T}C_{r}\Delta_{e^{A_{r}\tau}}=0$
(20)

Similarly, the relation between the perturbation $\Delta_{Q_{r,\tau}}$ and
$\Delta_{A_{r}}$ is through equation (14f).

$\displaystyle
A_{r}^{T}\Delta_{Q_{r,\tau}}+\Delta_{Q_{r,\tau}}A_{r}+\Delta_{A_{r}}^{T}Q_{r,\tau}+Q_{r,\tau}\Delta_{A_{r}}-$
(21)
$\displaystyle\Delta_{e^{A_{r}\tau}}^{T}C_{r}^{T}C_{r}e^{A_{r}\tau}-e^{A_{r}^{T}}C_{r}^{T}C_{r}\Delta_{e^{A_{r}\tau}}=0$

Applying Lemma 3.1 for equations (17b) and (20) we get,

$\text{Tr}\left(Y_{\tau}\Delta_{A_{r}}X^{T}+e^{A^{T}\tau}C^{T}C_{r}\Delta_{e^{A_{r}\tau}}X^{T}\right)=\text{Tr}\left(B_{r}B^{T}\Delta_{Y_{\tau}}\right)$
(22)

Similarly considering (17a) and (21) and using Lemma 3.1 we get,

$\displaystyle\text{Tr}(\Delta_{A_{r}}^{T}Q_{r,\tau}P_{r}+Q_{r,\tau}\Delta_{A_{r}}P_{r}-\Delta_{e^{A_{r}\tau}}^{T}C_{r}^{T}C_{r}e^{A_{r}\tau}P_{r}-$
(23) $\displaystyle
e^{A_{r}^{T}\tau}C_{r}^{T}C_{r}\Delta_{e^{A_{r}\tau}}P_{r})=\text{Tr}\left(B_{r}B_{r}^{T}\Delta_{Q_{r,\tau}}\right)$

Using (22) and (23) in the expression (19) we get,

$\begin{split}\Delta_{J}^{A_{r}}&=2\text{Tr}\left[\left(X^{T}Y_{\tau}+P_{r}Q_{r,\tau}\right)\Delta_{A_{r}}\right]\\\
&+2\text{Tr}\left[\left(X^{T}e^{A^{T}\tau}C^{T}C_{r}-P_{r}e^{A_{r}^{T}\tau}C_{r}^{T}C_{r}\right)\Delta_{e^{A_{r}\tau}}\right]\end{split}$
(24)

From (18), we get $2\text{
Tr}\left[\left(X^{T}e^{A^{T}\tau}C^{T}C_{r}-P_{r}e^{A_{r}^{T}\tau}C_{r}^{T}C_{r}\right)\Delta_{e^{A_{r}\tau}}\right]=2\text{
Tr}\left[S_{\tau}\Delta_{e^{A_{r}\tau}}\right]$.
The Fréchet derivative of matrix exponential $f(A)=e^{A}$ along a perturbation
matrix $E$ is defined as,

$L(A_{r},E)=\int_{0}^{1}e^{A_{r}(1-s)}Ee^{A_{r}s}ds$ (25)

Using the Fréchet derivative expression (25) we can express $2\text{
Tr}\left[S_{\tau}\Delta_{e^{A_{r}\tau}}\right]$ as follows,

$\begin{split}2\text{
Tr}\left[S_{\tau}\Delta_{e^{A_{r}\tau}}\right]&=2\tau\text{
Tr}\left[L(A_{r}\tau,S_{\tau})\Delta_{A_{r}}\right]\end{split}$ (26)

Using (26), we can rewrite (24) as follows

$\displaystyle\Delta_{J}^{A_{r}}$ $\displaystyle=2\text{
Tr}\left[\left(X^{T}Y_{\tau}+P_{r}Q_{r,\tau}\right)\Delta_{A_{r}}\right]+2\tau\text{Tr
}\left(L(A_{r}\tau,S_{\tau})\Delta_{A_{r}}\right)$ (27) $\displaystyle=\langle
2({Y_{\tau}}^{T}X+Q_{r,\tau}P_{r}+\tau(L(A_{r}\tau,S_{\tau})^{T})),\Delta_{A_{r}}\rangle$

Using the relation
$\Delta_{J}^{A_{r}}=\langle\Delta_{A_{r}}J,\Delta_{A_{r}}\rangle$ and (27), we
obtain (16a).

To get $\nabla_{B_{r}}J$, we perturb $B_{r}$ in the cost expression (15b). The
resulting first order perturbation is given by,

$\displaystyle\Delta_{J}^{B_{r}}$
$\displaystyle=\text{Tr}\left(2B^{T}Y_{\tau}\Delta_{B_{r}}+\Delta_{B_{r}}^{T}Q_{r,\tau}B_{r}+B_{r}^{T}Q_{r,\tau}\Delta_{B_{r}}\right)$
(28) $\displaystyle=\langle
2\left(Y_{\tau}^{T}B+Q_{r,\tau}B_{r}\right),\Delta_{B_{r}}\rangle$

Utilizing the relation
$\Delta_{J}^{B_{r}}=\langle\nabla_{B_{r}}J,\Delta_{B_{r}}\rangle$, we have
(16b).

Consider the error cost expression (15a). The first order perturbation in $J$
due to perturbation $\Delta_{C_{r}}$ of $C_{r}$ is,

$\displaystyle\Delta_{J}^{C_{r}}$
$\displaystyle=\text{Tr}\left(-2CX_{\tau}\Delta_{C_{r}}^{T}+\Delta_{C_{r}}P_{r,\tau}C_{r}^{T}+C_{r}P_{r,\tau}\Delta_{C_{r}}^{T}\right)$
(29)
$\displaystyle=\text{Tr}\left(2\left(C_{r}P_{r,\tau}-CX_{\tau}\right)^{T}\Delta_{C_{r}}\right)$
$\displaystyle=\langle
2\left(C_{r}P_{r,\tau}-CX_{\tau}\right),\Delta_{C_{r}}\rangle$

From the above expression we get $\nabla_{C_{r}}J$ (16c).

∎

######  Remark 1.

Reference [9] derives Lyapunov based and [20] derives interpolation based
first-order necessary conditions for time-limited H2 optimality. However, both
these derivations require that $A_{r}$ should be diagonalizable. Setting the
gradients that we derived in Theorem 3.2 to zero gives us another set of
Lyapunov based or Wilson’s first-order necessary condition for H2,τ
optimality. We do not assume diagonalizability of $A_{r}$ for deriving the
optimality conditions. It can be easily verified that our optimality
conditions are equivalent to the Lyapunov based optimality conditions of [9]
if $A_{r}$ is diagonalizable. Further, diagonalizability of $A_{r}$ also
ensures that our optimality conditions are equivalent to the interpolation
based H2,τ optimality conditions derived in [20]. We prove this in the
Appendix.

###  3.2 A Numerical method for H2,τ model reduction

The optimization problem (8) considers
$J=\left\|G-G_{r}\right\|_{\text{H}_{2,\tau}}^{2}$ as the objective function
and $\\{A_{r},B_{r},C_{r}\\}$ as the optimization variables. The optimization
problem is nonlinear and non-convex. Hence finding global minimizers is
difficult. However, we can use standard nonlinear optimization techniques with
good initial conditions to obtain local minimizers [14].

In this work, we solve the above optimization problem using standard quasi-
Newton solvers by employing the MATLAB function ’fminunc’. The Hessian matrix
is updated by means of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm .
The gradients required for the BFGS algorithm are calculated using the closed
form expressions derived in Section 3.1 (Equations 16a), (16b) and (16c). Due
to the non-convex nature of the optimization problem, a good starting point is
very necessary to solve the optimization problem. We use TL-BT and TL-TSIA to
reduce the original model and use the reduced model to initialize the
optimization problem. When the convergence criteria becomes less than a preset
error tolerance, the iterations are stopped. We name the proposed time-limited
model reduction method as TL-H2Opt.

######  Remark 2.

From Remark 1, we note the equivalence of the Lyapunov based and interpolation
based frameworks of optimality conditions. Thus, for a fixed finite-time
interval $\begin{bmatrix}0&\tau\end{bmatrix}$ we require a minimum of $r(m+p)$
parameters for representing the H2,τ optimality conditions. This is similar to
the infinite interval case [21]. We have $r^{2}+r(m+p)$ parameters in our
optimization problem, which leads to overparametrization. However, this does
not impede obtaining better H2,τ optimal models due to the nature of the
quasi-Newton solvers used for solving the optimization problem as observed in
[15, 13].

###  3.3 Computational Cost

We now discuss the computational cost of the proposed method, TL-H2Opt.
Computation of $P_{\tau}$ and $Q_{\tau}$ requires solving the Lyapunov
equations (14a) and (14d) respectively. The equations are solved by the MATLAB
function ’lyap’. The underlying algorithm for ’lyap’ has a computational
complexity of $\mathcal{O}(n^{3})$. The computation of these quantities is
costly. However, they are independent of optimization variables
($\\{A_{r},B_{r},C_{r}\\}$) and hence need to be computed only once before the
start of the optimization process. The terms $P_{r,\tau}$ and $Q_{r,\tau}$ are
dependent on the optimization variables and need to be computed at every
iteration of the optimization process. Both these terms have computing cost
$\mathcal{O}(n_{r}^{3})$. Since $n_{r}\ll n$, the reduced order gramians are
not computationally heavy. The exponential term $e^{A\tau}$ is computed with
the MATLAB function ’expm’ and has a high computation complexity of
$\mathcal{O}(n^{3})$. However, this term needs to be computed only once since
it doesn’t involve the optimization variables. The terms $e^{A_{r}\tau}$ and
$L(e^{A_{r}\tau},S_{\tau})$ include optimization variables and need to be
computed at every iteration of the optimization process. For computing these
terms, we use Algorithm 3 of [4] which has a computational cost of
$\mathcal{O}(n_{r}^{3})$.

The terms $Y_{\tau}$, $X_{\tau}$ and $X$ are required for calculating the cost
function as well as the gradients and have to be computed at every iteration
of the optimization problem. They are solutions of the Sylvester equations
(14e), (14b) and (17b) respectively. Computing them with the ’lyap’ function
in MATLAB costs $\mathcal{O}(n^{3})$. This method of computing $Y_{\tau}$ and
$X$ works for medium scale systems (order $<$ 1000) but becomes
computationally expensive for large-scale models (order $>$ 1000). We can
speed up the computations by using Algorithm 3 of [5] to compute the Sylvester
matrices. In this case, the cost of solving the Sylvester equations is much
less than $\mathcal{O}(n^{3})$ if the matrix $A$ is diagonal or has some
sparse structure.

##  4 NUMERICAL EXAMPLES

In this section, we investigate the performance of the proposed algorithm
TL-H2Opt using two numerical examples. The first example is a SISO model of a
beam with order 348. The second example is a MIMO model of the International
Space Station (ISS) with three inputs and three outputs and order 270. The
examples are taken from [6]. The simulations are done in MATLAB version
8.3.0.532(R2014a) on a Intel(R) Core(TM) i5-6500 CPU @ 3.20GHz 3.19 GHz system
with 16 GB RAM. We reduce the models over fixed finite time intervals. Using
TL-TSIA and TL-BT, reduced models are obtained. These reduced models are
further used to initialize TL-H2Opt. The improvement in performance is
assessed using the quantity $\Delta\text{Err}(\%)$ defined as

$\Delta\text{Err}(\%)=\frac{\text{Err\textsubscript{Alg}-Err\textsubscript{Opt}}}{\text{Err\textsubscript{Alg}}}\times
100\%$ (30)

where ErrAlg is the H2,τ approximation error obtained by the algorithm Alg and
ErrOpt is the error obtained by TL-H2Opt with Alg initialization. The
algorithm Alg may refer to either TL-TSIA or TL-BT in our case.

###  4.1 Beam Example

The first example is a clamped beam model of order 348 with one input and one
output. We fix a time interval of $\begin{bmatrix}0&1\end{bmatrix}$. For this
time interval, we use algorithms TL-TSIA and TL-BT to obtain reduced order
approximations with $r$ varying from $r=2$ to $r=21$ increasing the value of
$r$ one at a time. We use these low order models to initialize the TL-H2Opt
algorithm. The first subplot of Figure 1 displays the H${}_{2}(\tau)$ errors
for the reduced models obtained by TL-TSIA and TL-H2Opt with TL-TSIA
initialization while the second displays the improvement in performance of the
optimization based algorithm over TL-TSIA given by (30) with Alg = TL-TSIA.
Similarly, the first subplot of Figure 2 compares the approximation errors due
to TL-BT and TL-H2Opt with TL-BT initialization and second shows the
improvement in performance of TL-BT due to the time limited H2 optimization
algorithm.

The reduced models obtained by TL-TSIA and TL-BT do not satisfy the
H${}_{2}(\tau)$ optimality conditions exactly[ref]. The reduced models
obtained using TL-H2Opt with TL-TSIA and TL-BT initialization improve the
H${}_{2}(\tau)$ approximation errors as evident from Figure 1 and Figure 2.
For $r=4,18,20$, the $H_{2}(\tau)$ approximation error due to TL-TSIA is high
or it doesn’t converge and hence errors corresponding to those orders are not
included in the first subplot of Figure 1. For reduced models of order ($r$)
less than $16$ some of the reduced orders show good improvement in the
H${}_{2}(\tau)$ errors; for instance in case of TL-TSIA $r=5,9$ and $10$ show
an improvement of $79.82\%$, $76.26\%$ and $91.92\%$ respectively. For TL-BT
$r=5,8,9$ and $14$ have their H${}_{2}(\tau)$ errors reduced by $70.90\%$,
$74.74\%$, $80.71\%$ and $70.94\%$ percent. Beyond order $r=16$, the
optimization algorithm doesn’t lead to any significant improvement in the
H${}_{2}(\tau)$ approximation errors for both TL-TSIA and TL-BT
initialization.

Figure 1: Performance of TL-H2Opt with TL-TSIA init for Beam Example Figure 2:
Performance of TL-H2Opt with TL-BT init for Beam Example

###  4.2 ISS Example

The second example is a MIMO model of the International Space Station (ISS)
with three inputs and three outputs. For this example we consider a time
interval of $\begin{bmatrix}0&0.5\end{bmatrix}$.

The H2,τ approximation errors for the reduced models of various orders
obtained with TL-TSIA and TL-H2Opt with TL-TSIA initialization along with
performance improvement $\Delta\text{Err}(\%)$(30) due to the optimization
procedure is shown in Figure 3. Similar comparisons for TL-BT and TL-H2Opt
with TL-BT initialization are shown in Figure 4. We observe that the H2,τ
errors are improved considerably due to the application of TL-H2Opt for
reduced orders less than $r=38$ for both TL-TSIA and TL-BT initialization.
There is no substantial improvement in the performance of TL-TSIA and TL-BT
due to the application of TL-H2Opt for reduced models of order greater than
$r=38$.

The reduced models obtained by applying TL-TSIA for orders $r=4$ and $r=12$
have high H2,τ approximation errors and are not shown in Figure 3.
Initializing the optimization algorithm TL-H2Opt by the reduced model obtained
with TL-TSIA and solving the optimization problem improves the H2,τ errors for
all concerned reduced model orders ($r<38$) including $r=4$ and $r=12$ as
evident from the second subplot of Figure 3. Apart from the case of $r=4$ and
$r=12$ where the performance improvement is nearly $100\%$, the reduced order
models with $r=22,26,28,30,32$ have performance improvement of over $70\%$.
Unlike the TL-TSIA initialization, there is no considerably high H2,τ errors
for any reduced order $r$ due to TL-BT initialization. Due to the application
of the time-limited optimization algorithm with TL-BT reduced models as
initial points, the H2,τ error reduces as evident from Figure 4. The decrease
in H2,τ error for $r=18,20,22,24,26,28,32$ is more than $70\%$.

Figure 3: Performance of TL-H2Opt with TL-TSIA init for ISS Example Figure 4:
Performance of TL-H2Opt with TL-BT init for ISS Example

##  5 Conclusion

In this work, we propose an optimization based method to obtain H2,τ optimal
reduced models. We derive closed form expressions of the gradients of an
objective function defined over a limited time interval. The gradients are
used with standard optimization algorithms for minimizing the objective
function. The model reduction method proposed involves two steps. We first
obtain a reduced model via TL-TSIA or TL-BT and then use the reduced model to
initialize the optimization algorithm. Through a numerical example, we
demonstrate the superiority of the TL-H2Opt algorithm over TL-TSIA and TL-BT
in obtaining better H2,τ optimal reduced models.

## References

  * [1] P. Benner  _et al._ Model Order Reduction: Volume 2: Snapshot-Based Methods and Algorithms.  De Gruyter, 2020.
  * [2] P. Benner  _et al._ Model Order Reduction: Volume 3: Applications.  De Gruyter, 2020.
  * [3] P. Benner  _et al._ System- and Data-Driven Methods and Algorithms.  De Gruyter, 2021.
  * [4] Awad H Al-Mohy and Nicholas J Higham.  Computing the Fréchet derivative of the matrix exponential, with an application to condition number estimation.  SIAM J. Matrix Analy. Appl., 30(4):1639–1657, 2009.
  * [5] Peter Benner, Martin Køhler, and Jens Saak.  Sparse-dense Sylvester equations in ${H}_{2}$-model order reduction.  Technical Report MPIMD/11-11, Max Planck Institute, Madeburg, Germany, December 2011.
  * [6] Younes Chahlaoui and Paul Van Dooren.  Benchmark examples for model reduction of linear time-invariant dynamical systems.  In Dimension reduction of large-scale systems, pages 379–392. Jan 2005.
  * [7] Igor Pontes Duff and Patrick Kürschner.  Numerical computation and new output bounds for time-limited balanced truncation of discrete-time systems.  Lin. Alg. Appl., 623:367–397, 2021.
  * [8] Wodek Gawronski and Jer-Nan Juang.  Model reduction in limited time and frequency intervals.  Int. J. Syst. Sci., 21(2):349–376, 1990.
  * [9] Pawan Goyal and Martin Redmann.  Time-limited ${H}_{2}$-optimal model order reduction.  App. Math. Comp., 355:184–197, 2019.
  * [10] Serkan Gugercin, Athanasios C Antoulas, and Christopher Beattie.  ${H}_{2}$ model reduction for large-scale linear dynamical systems.  SIAM J. Matrix Analy. Appl., 30(2):609–638, 2008.
  * [11] Philip Holmes, John L Lumley, Gahl Berkooz, and Clarence W Rowley.  Turbulence, coherent structures, dynamical systems and symmetry.  Cambridge Uni. Press, 2012.
  * [12] Patrick Kürschner.  Balanced truncation model order reduction in limited time intervals for large systems.  Adv. Comp. Math., 44(6):1821–1844, 2018.
  * [13] Tomas McKelvey and Anders Helmersson.  System identification using an over-parametrized model class-improving the optimization algorithm.  In Proc. $36^{th}$ IEEE Conf. Decision Control, volume 3, pages 2984–2989, 1997.
  * [14] Jorge Nocedal and Stephen Wright.  Numerical optimization.  Springer Science & Business Media, 2006.
  * [15] Daniel Petersson and Johan Löfberg.  Model reduction using a frequency-limited H2-cost.  Syst. Contr. Lett., 67:32–39, 2014.
  * [16] Martin Redmann.  An $L_{T}^{2}$-error bound for time-limited balanced truncation.  Syst. Contr. Lett., 136:104620, 2020.
  * [17] Martin Redmann and Patrick Kürschner.  An output error bound for time-limited balanced truncation.  Syst. & Contr. Lett., 121:1–6, 2018.
  * [18] Hiroyuki Sato and Kazuhiro Sato.  Riemannian trust-region methods for ${H}_{2}$ optimal model reduction.  In  $54^{th}$ IEEE Conf. Decision Control, pages 4648–4655, 2015\.
  * [19] Kazuhiro Sato.  Riemannian optimal model reduction of linear second-order systems.  IEEE Contr. Syst. Lett., 1(1):2–7, 2017.
  * [20] Klajdi Sinani and Serkan Gugercin.  H${}_{2}(t_{f})$ optimality conditions for a finite-time horizon.  Automatica, 110:108604, 2019.
  * [21] Paul Van Dooren, Kyle A Gallivan, and P-A Absil.  H2-optimal model reduction of MIMO systems.  Appl. Math. Lett., 21(12):1267–1273, 2008.
  * [22] DA Wilson.  Optimum solution of model-reduction problem.  In Proc. of the Inst. of Electr. Eng., volume 117, pages 1161–1165, 1970.
  * [23] Yuesheng Xu and Taishan Zeng.  Optimal $H_{2}$ model reduction for large scale MIMO systems via Tangential Interpolation.  Int. J. Numer. Analy. Mod., 8(1), 2011.
  * [24] Wei-Yong Yan and James Lam.  An approximate approach to H2 optimal model reduction.  IEEE Trans. Automat. Contr., 44(7):1341–1358, 1999.

## Appendix

Equivalence of Lyapunov and Tangential Interpolation based H${}_{2}(\tau)$
Optimality conditions

The partial fraction expansion of the transfer function matrix $H_{r}(s)$ of
system (6) is given by

$H_{r}(s)=\sum_{i=1}^{r}\frac{c_{i}b_{i}}{s-\lambda_{i}}$ (31)

where $\lambda_{i}\in\mathbb{C}$, $c_{i}\in\mathbb{C}^{p\times 1}$,
$b_{i}\in\mathbb{C}^{1\times m}$ and $(\lambda_{i},b_{i},c_{i})$ for
$i=1,2,\ldots,r$ is self conjugate. Let, $\lambda_{i}^{*}$, $c_{i}^{T}$ and
$b_{i}^{T}$ be the conjugate transpose of $\lambda_{i}$, $c_{i}$ and $b_{i}$
respectively. Let $v_{i}$ and $w_{i}$ be the right and left eigenvector
respectively of the matrix $A_{r}$ corresponding to the eigenvalue
$\lambda_{i}$. The following relations hold.

$A_{r}v_{i}=\lambda_{i}v_{i},C_{r}v_{i}=c_{i},w_{i}^{T}A_{r}=\lambda_{i}w_{i}^{T},w_{i}^{T}B_{r}=b_{i}$
(32)

For the transfer function $H(s)$ of system (1), the time-limited counterpart
over the time-interval $\begin{bmatrix}0&\tau\end{bmatrix}$ is given by [20]
as,

$H_{\tau}(s)=C(sI_{n}-A)^{-1}(I_{n}-e^{-s\tau}e^{A\tau})B$ (33)

$H_{r,\tau}(s)$ is defined similarly. We denote the identity matrix of size
$n$ and $r$ as $I_{n}$ and $I_{r}$ respectively and define the matrix
$V=\begin{bmatrix}v_{1}&v_{2}&\ldots&v_{r}.\end{bmatrix}$.

######  Theorem 5.1.

Let $H_{r}(s)$ given by (31) have $r$ distinct first order poles.
$H_{\tau}(s)$ and $H_{r,\tau}(s)$ are the time limited transfer functions of
$H(s)$ and $H_{r}(s)$ respectively over the time interval [0 $\tau$]. Then for
$i=1,2,\ldots,r$, $j=1,2,\ldots,r$ and $i\neq j$

$\frac{1}{2}\left(\nabla_{B_{r}}J\right)^{T}v_{i}=\left[H_{r,\tau}^{T}(-\lambda_{i}^{*})-H_{\tau}^{T}(-\lambda_{i}^{*})\right]c_{i}$
(34a)
$\frac{1}{2}w_{i}^{T}\left(\nabla_{C_{r}}J\right)^{T}=b_{i}\left[H_{r,\tau}^{T}(-\lambda_{i}^{*})-H_{\tau}^{T}(-\lambda_{i}^{*})\right]$
(34b)
$\frac{1}{2}w_{i}^{T}\left(\nabla_{A_{r}}J\right)^{T}v_{i}=b_{i}\left.\frac{d}{ds}\left[H_{\tau}^{T}(s)-H_{r,\tau}^{T}(s)\right]\right|_{s=-\lambda_{i}^{*}}c_{i}$
(34c)
$\frac{1}{2}w_{i}^{T}\left(\nabla_{A_{r}}J\right)^{T}v_{j}=\frac{1}{2(\lambda_{i}-\lambda_{j})}\left[b_{i}\left(\nabla_{B_{r}}J\right)^{T}v_{j}-w_{i}^{T}\left(\nabla_{C_{r}}J\right)^{T}c_{j}\right]$
(34d)

###### Proof.

Let us define
$x_{i,\tau}=X_{\tau}w_{i},x_{i}=Xw_{i},p_{i,\tau}=P_{r,\tau}w_{i},p_{i}=P_{r}w_{i},y_{i,\tau}=Y_{\tau}v_{i},q_{i,\tau}=Q_{r,\tau}v_{i}$.
Using (14a), (14b), (14c), (14d), (14e) and (14f) we get

$x_{i,\tau}=-(A+\lambda_{i}^{*}I_{n})^{-1}(I_{n}-e^{\lambda_{i}^{*}\tau}e^{A\tau})Bb_{i}^{T}$
(35a) $x_{i}=-(A+\lambda_{i}^{*}I_{n})^{-1}Bb_{i}^{T}$ (35b)
$p_{i,\tau}=-(A_{r}+\lambda_{i}^{*}I_{r})^{-1}(I_{r}-e^{\lambda_{i}^{*}\tau}e^{A_{r}\tau})Bb_{i}^{T}$
(35c) $p_{i}=-(A_{r}+\lambda_{i}^{*}I_{r})^{-1}B_{r}b_{i}^{T}$ (35d)
$y_{i,\tau}=(A^{T}+\lambda_{i}I_{n})^{-1}(I_{n}-e^{\lambda_{i}\tau}e^{A^{T}\tau})C^{T}c_{i}$
(35e)
$q_{i,\tau}=-(A_{r}^{T}+\lambda_{i}I_{r})^{-1}(I_{r}-e^{\lambda_{i}\tau}e^{A_{r}^{T}\tau})C_{r}^{T}c_{i}$
(35f)

We use (35f) and (35e) and obtain (34a) as follows,

$\displaystyle\frac{1}{2}\left(\nabla_{B_{r}}J\right)^{T}v_{i}$
$\displaystyle=B_{r}^{T}q_{i,\tau}+B^{T}y_{i,\tau}$ (36)
$\displaystyle=-B_{r}^{T}(A_{r}^{T}+\lambda_{i}I_{r})^{-1}(I_{r}-e^{\lambda_{i}\tau}e^{A_{r}^{T}\tau})C_{r}^{T}c_{i}+$
$\displaystyle
B^{T}(A^{T}+\lambda_{i}I_{n})^{-1}(I_{n}-e^{\lambda_{i}\tau}e^{A^{T}\tau})C^{T}c_{i}$
$\displaystyle=[H_{r,\tau}^{T}(-\lambda_{i}^{*})-H_{\tau}^{T}(-\lambda_{i}^{*})]c_{i}$

Similarly, using (35c) and (35a) the second expression becomes

$\displaystyle\frac{1}{2}(\nabla_{C_{r}}J)w_{i}$
$\displaystyle=[H_{r,\tau}(-\lambda_{i}^{*})-H_{\tau}(-\lambda_{i}^{*})]b_{i}^{T}$
(37)

The third expression is derived as follows.

$\displaystyle\frac{1}{2}w_{i}^{T}\left(\nabla_{A_{r}}J\right)^{T}v_{j}$ (38)
$\displaystyle=x_{i}^{T}y_{j,\tau}+p_{i}^{T}q_{j,\tau}+w_{i}^{T}S_{\tau}v_{j}\int_{0}^{1}e^{\lambda_{i}(\tau-\tau
s)}e^{\lambda_{j}\tau s}(\tau ds)$

Using (33), we have
$\frac{d}{ds}H_{\tau}(s)=-C(sI_{n}-A)^{-2}(I_{n}-e^{-s\tau}e^{A\tau})B+\tau
e^{-s\tau}C(sI_{n}-A)^{-1}e^{A\tau}B$ and
$\frac{d}{ds}H_{r,\tau}(s)=-C_{r}(sI_{r}-A_{r})^{-2}(I_{r}-e^{-s\tau}e^{A_{r}\tau})B_{r}+\tau
e^{-s\tau}C_{r}(sI_{r}-A_{r})^{-1}e^{A_{r}\tau}B_{r}$. Substituting $S_{\tau}$
for $i=j$, the expression (38) reduces to

$\displaystyle x_{i}^{T}y_{i,\tau}+p_{i}^{T}q_{i,\tau}+$ (39)
$\displaystyle(x_{i}^{T}e^{A^{T}\tau}C^{T}c_{i}-p_{i}^{T}e^{A_{r}^{T}\tau}C_{r}^{T}c_{i})\int_{0}^{1}e^{\lambda_{i}(\tau-\tau
s)}e^{\lambda_{i}\tau s}\tau ds$

Using (35b) and (35e), the term $x_{i}^{T}y_{i,\tau}$ becomes
$-b_{i}B^{T}(A^{T}+\lambda_{i}I_{n})^{-2}(I_{n}-e^{\lambda_{i}\tau}e^{A^{T}\tau})C^{T}c_{i}$,
using (35c) and (35f) the term $p_{i}^{T}q_{i,\tau}^{T}$ becomes
$b_{i}B_{r}^{T}(A_{r}^{T}+\lambda_{i}I_{r})^{-2}(I_{r}-e^{\lambda_{i}\tau}e^{A_{r}^{T}\tau})C_{r}^{T}c_{i}$
and after substituting $x_{i}$ and $p_{i}$, the third term in the expression
(39) becomes $-\tau
e^{\lambda_{i}\tau}b_{i}B^{T}(A^{T}+\lambda_{i}I_{n})^{-1}e^{A^{T}\tau}C^{T}c_{i}+\tau
e^{\lambda_{i}\tau}b_{i}B_{r}^{T}(A_{r}^{T}+\lambda_{i}I_{r})^{-1}e^{A_{r}^{T}\tau}C_{r}^{T}c_{i}$.
Combining the three terms we get (34c).
For the case $i\neq j$, the expression (38) becomes

$\displaystyle\frac{1}{2}w_{i}^{T}\left(\nabla_{A_{r}}J\right)^{T}v_{j}$ (40)
$\displaystyle=x_{i}^{T}y_{j,\tau}+p_{i}^{T}q_{j,\tau}+w_{i}^{T}S_{\tau}v_{j}\left(\frac{e^{\lambda_{i}\tau}-e^{\lambda_{j}\tau}}{\lambda_{i}-\lambda_{j}}\right)$

Using the identity
$\left(A^{T}+\lambda_{i}I_{n}\right)^{-1}\left(A^{T}+\lambda_{j}I_{n}\right)^{-1}=\frac{1}{\lambda_{i}-\lambda_{j}}\left[\left(A^{T}+\lambda_{j}I_{n}\right)^{-1}-\left(A^{T}+\lambda_{i}I_{n}\right)^{-1}\right]$
and (35b), (35e) , the term $x_{i}^{T}y_{j,\tau}$ becomes

$\displaystyle-\frac{-1}{\lambda_{i}-\lambda_{j}}b_{i}B^{T}\left[\left(A^{T}+\lambda_{j}I_{n}\right)^{-1}-\left(A^{T}+\lambda_{i}I_{n}\right)^{-1}\right]$
(41)
$\displaystyle\left(I_{n}-e^{\lambda_{j}\tau}e^{A^{T}\tau}\right)C^{T}c_{j}$

Similarly, the term $p_{i}^{T}q_{j,\tau}$ becomes

$\displaystyle\frac{-1}{\lambda_{i}-\lambda_{j}}b_{i}B_{r}^{T}\left[\left(A_{r}^{T}+\lambda_{j}I_{r}\right)^{-1}-\left(A_{r}^{T}+\lambda_{i}I_{r}\right)^{-1}\right]$
(42)
$\displaystyle\left(I_{r}-e^{\lambda_{j}\tau}e^{A_{r}^{T}\tau}\right)C_{r}^{T}c_{j}$

The third right hand side term of (40) becomes

$\displaystyle\frac{e^{\lambda_{i}\tau}-e^{\lambda_{j}\tau}}{\lambda_{i}-\lambda_{j}}\left[\left(-b_{i}B^{T}(A^{T}+\lambda_{i}I_{n})^{-1}e^{A^{T}\tau}C^{T}c_{j}\right)\right]+$
(43)
$\displaystyle\frac{e^{\lambda_{i}\tau}-e^{\lambda_{j}\tau}}{\lambda_{i}-\lambda_{j}}\left[\left(b_{i}B_{r}^{T}(A_{r}^{T}+\lambda_{i}I_{r})^{-1}e^{A_{r}^{T}\tau}C_{r}^{T}c_{j}\right)\right]$

Adding the expressions (41), (42) and (43), we obtain (34d). ∎

The above theorem shows that setting diag
$V^{-1}\left(\nabla_{A_{r}}J\right)^{T}V$, $\left(\nabla_{B_{r}}J\right)^{T}$
and $\left(\nabla_{C_{r}}J\right)^{T}$ to 0 gives the interpolation based H2,τ
optimality conditions.This proves that Lyapunov based and interpolation based
optimality conditions are equivalent.