# Dynamical Majorana Ising spin response in a topological superconductor-
magnet hybrid by microwave irradiation

Yuya Ominato  Kavli Institute for Theoretical Sciences, University of Chinese
Academy of Sciences, Beijing, 100190, China.  Waseda Institute for Advanced
Study, Waseda University, Shinjuku, Tokyo 169-8050, Japan.  Ai Yamakage
Department of Physics, Nagoya University, Nagoya 464-8602, Japan  Mamoru
Matsuo  Kavli Institute for Theoretical Sciences, University of Chinese
Academy of Sciences, Beijing, 100190, China.  CAS Center for Excellence in
Topological Quantum Computation, University of Chinese Academy of Sciences,
Beijing 100190, China  Advanced Science Research Center, Japan Atomic Energy
Agency, Tokai, 319-1195, Japan  RIKEN Center for Emergent Matter Science
(CEMS), Wako, Saitama 351-0198, Japan

###### Abstract

We study a dynamical spin response of surface Majorana modes in a topological
superconductor-magnet hybrid under microwave irradiation. We find a method to
toggle between dissipative and non-dissipative Majorana Ising spin dynamics by
adjusting the external magnetic field angle and the microwave frequency. This
reflects the topological nature of the Majorana modes, enhancing the Gilbert
damping of the magnet, thereby, providing a detection method for the Majorana
Ising spins. Our findings illuminate a magnetic probe for Majorana modes,
paving the path to innovative spin devices.

Introduction.—The quest for Majoranas within matter stands as one of the
principal challenges in the study of condensed matter physics, more so in the
field of quantum many-body systems Wilczek (2009). The self-conjugate nature
of Majoranas leads to peculiar electrical characteristics that have been the
subject of intensive research, both theoretical and experimental Yazdani  _et
al._ (2023). In contrast, the focus of this paper lies on the magnetic
properties of Majoranas, specifically the Majorana Ising spin Chung and Zhang
(2009); Sato and Fujimoto (2009); Nagato  _et al._ (2009); Shindou  _et al._
(2010); Shiozaki and Sato (2014); Xiong  _et al._ (2017). A distinctive
characteristic of Majorana modes, appearing as a surface state in topological
superconductors (TSC), is its exceedingly strong anisotropy, which makes it
behave as an Ising spin. In particular, this paper proposes a method to
explore the dynamical response of the Majorana Ising spin through the exchange
interaction at the magnetic interface, achieved by coupling the TSC to a
ferromagnet with ferromagnetic resonance (FMR) (as shown in Fig.1 (a)).

FMR modulation in a magnetic hybrid system has attracted much attention as a
method to analyze spin excitations in thin-film materials attached to magnetic
materials Qiu  _et al._ (2016); Han  _et al._ (2020). Irradiating a magnetic
material with microwaves induces dynamics of localized spin in magnetic
materials, which can excite spins in adjacent thin-film materials via the
magnetic proximity effect. This setup is called spin pumping, and has been
studied intensively in the field of spintronics as a method of injecting spins
through interfaces Tserkovnyak  _et al._ (2002, 2005). Recent studies have
theoretically proposed that spin excitation can be characterized by FMR in
hybrid systems of superconducting thin films and magnetic materials Inoue  _et
al._ (2017); Kato  _et al._ (2019); Silaev (2020a, b); Ominato  _et al._
(2022a, b). Therefore, it is expected to be possible to analyze the dynamics
of surface Majorana Ising spins using FMR in hybrid systems.

Figure 1:  (a) The TSC-FI hybrid schematic reveals how, under resonance
frequency microwave irradiation, localized spins commence precessional motion,
consequently initiating the dynamical Majorana Ising spin response at the TSC
interface. (b) In the TSC context, the liaison between a spin-up electron and
a spin-down hole with the surrounding sea of spin-triplet Cooper pairs
drastically modulate their properties; notably, a spin-down hole can engage
with a spin-triplet Cooper pair, thereby inheriting a negative charge. (c)
Notably, spin-triplet Cooper pairs amass around holes and scatter around
electrons, thereby eroding the rigid distinction between the two. (d) The
interplay between the Majorana mode and the localized spin manipulates the FMR
spectrum, triggering a frequency shift and linewidth broadening.

In this work, we consider a TSC-ferromagnetic insulator (FI) hybrid system as
shown in Fig. 1 (a). The FMR is induced by microwave irradiation on the FI. At
the interface between the TSC and the FI, the surface Majorana modes interact
with the localized spins in the FI. As a result, the localized spin dynamics
leads to the dynamical Majorana Ising spin response (DMISR), which means the
Majorana Ising spin density is dynamically induced, and it is possible to
toggle between dissipative and non-dissipative Majorana Ising spin dynamics by
adjusting the external magnetic field angle and the microwave frequency.
Furthermore, the modulation of the localized spin dynamics due to the
interface interaction leads to a frequency shift and a linewidth broadening,
which reflect the properties of the Majorana Ising spin dynamics. This work
proposes a setup for detecting Majorana modes and paves the way for the
development of quantum computing and spin devices using Majoranas.

Model.—We introduce a model Hamiltonian $H$ consisting of three terms

$\displaystyle H=H_{\mathrm{M}}+H_{\mathrm{FI}}+H_{\mathrm{ex}}.$ (1)

The first, second, and third terms respectively describe the surface Majorana
modes on the TSC surface, the bulk FI, and the proximity-induced exchange
coupling. Our focus is on energy regions significantly smaller than the bulk
superconducting gap. This focus allows the spin excitation in the TSC to be
well described using the surface Majorana modes. The subsequent paragraphs
provide detailed explanations of each of these three terms.

The first terms $H_{\mathrm{M}}$ describes the surface Majorana modes,

$\displaystyle H_{\mathrm{M}}=\frac{1}{2}\int d{\bf\it
r}{\bf\it\psi}^{\mathrm{T}}({\bf\it r})\quantity(\hbar
v\hat{k}_{y}\sigma^{x}-\hbar v\hat{k}_{x}\sigma^{y}){\bf\it\psi}({\bf\it r}),$
(2)

where ${\bf\it r}=(x,y)$, $\hat{{\bf\it k}}=(-i\partial_{x},-i\partial_{y})$,
$v$ is a constant velocity, and
${\bf\it\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ are the Pauli matrices.
The two component Majorana field operator is given by ${\bf\it\psi}({\bf\it
r})=\quantity(\psi_{\to}({\bf\it r}),\psi_{\leftarrow}({\bf\it
r}))^{\mathrm{T}}$, with the spin quantization axis along the $x$ axis. The
Majorana field operators satisfy the Majorana condition $\psi_{\sigma}({\bf\it
r})=\psi_{\sigma}^{\dagger}({\bf\it r})$ and the anticommutation relation
$\\{\psi_{\sigma}({\bf\it r}),\psi_{\sigma^{\prime}}({\bf\it
r})\\}=\delta_{\sigma\sigma^{\prime}}\delta({\bf\it r}-{\bf\it r}^{\prime})$
where $\sigma,\sigma^{\prime}=\to,\leftarrow$. We can derive $H_{\mathrm{M}}$
by using surface-localized solutions of the BdG equation based on the bulk TSC
Hamiltonian. The details of the derivation of $H_{\mathrm{M}}$ are provided in
the Supplemental Material sm .

A notable feature of the surface Majorana modes is that the spin density is
Ising like, which we call the Majorana Ising spin Chung and Zhang (2009); Sato
and Fujimoto (2009); Nagato  _et al._ (2009); Shindou  _et al._ (2010);
Shiozaki and Sato (2014); Xiong  _et al._ (2017). The feature follows
naturally from the Majorana condition and the anticommutation relation. The
Majorana Ising spin density operator is given by ${\bf\it s}({\bf\it
r}):={\bf\it\psi}^{\mathrm{T}}({\bf\it
r})\quantity({\bf\it\sigma}/2){\bf\it\psi}({\bf\it
r})=(0,0,-i\psi_{\to}({\bf\it r})\psi_{\leftarrow}({\bf\it r}))$ (See the
Supplemental Material for details sm ). The anisotropy of the Majorana Ising
spin is the hallmark of the surface Majorana modes on the TSC surface.

The second term $H_{\mathrm{FI}}$ descries the bulk FI and is given by the
ferromagnetic Heisenberg model,

$\displaystyle H_{\mathrm{FI}}=$ $\displaystyle-\mathcal{J}\sum_{\langle
n,m\rangle}{\bf\it S}_{n}\cdot{\bf\it S}_{m}-\hbar\gamma
h_{\mathrm{dc}}\sum_{n}S_{n}^{Z},$ (3)

where $\mathcal{J}>0$ is the exchange coupling constant, ${\bf\it S}_{n}$ is
the localized spin at site $n$, $\langle n,m\rangle$ means summation for
nearest neighbors, $\gamma$ is the electron gyromagnetic ratio, and
$h_{\mathrm{dc}}$ is the static external magnetic field. We consider the spin
dynamics of the localized spin under microwave irradiation, applying the spin-
wave approximation. This allows the spin excitation to be described by a free
bosonic operator, known as a magnon Holstein and Primakoff (1940).

The third term $H_{\mathrm{ex}}$ represents the proximity exchange coupling at
the interface between the TSC and the FI,

$\displaystyle H_{\mathrm{ex}}=-\int d{\bf\it r}\sum_{n}J({\bf\it r},{\bf\it
r}_{n}){\bf\it s}({\bf\it r})\cdot{\bf\it
S}_{n}=H_{\mathrm{Z}}+H_{\mathrm{T}},$ (4) $\displaystyle
H_{\mathrm{Z}}=-\cos\vartheta\int d{\bf\it r}\sum_{n}J({\bf\it r},{\bf\it
r}_{n})s^{z}({\bf\it r})S^{Z}_{n},$ (5) $\displaystyle
H_{\mathrm{T}}=-\sin\vartheta\int d{\bf\it r}\sum_{n}J({\bf\it r},{\bf\it
r}_{n})s^{z}({\bf\it r})S_{n}^{X},$ (6)

where the angle $\vartheta$ is shown in Fig. 1 (a). $H_{\mathrm{Z}}$ is the
coupling along the precession axis and $H_{\mathrm{T}}$ is the coupling
perpendicular to the precession axis. In our setup, $H_{\mathrm{Z}}$ leads to
gap opening of the energy spectrum of the surface Majorana modes and
$H_{\mathrm{T}}$ gives the DMISR under the microwave irradiation.

Dynamical Majorana Ising spin response.—We consider the microwave irradiation
on the FI. The coupling between the localized spins and the microwave is given
by

$\displaystyle V(t)=-\hbar\gamma
h_{\mathrm{ac}}\sum_{n}\quantity(S^{X}_{n}\cos{\omega t}-S^{Y}_{n}\sin{\omega
t}),$ (7)

where $h_{\mathrm{ac}}$ is the microwave amplitude, and $\omega$ is the
microwave frequency. The microwave irradiation leads to the precessional
motion of the localized spin. When the frequency of the precessional motion
and the microwave coincide, the FMR occurs. The FMR leads to the DMISR via the
exchange interaction. The DMISR is characterized by the dynamic spin
susceptibility of the Majorana modes, $\tilde{\chi}^{zz}({\bf\it q},\omega)$,
defined as

$\displaystyle\tilde{\chi}^{zz}({\bf\it q},\omega):=\int d{{\bf\it
r}}e^{-i{\bf\it q}\cdot{\bf\it r}}\int dte^{i(\omega+i0)t}\chi^{zz}({\bf\it
r},t),$ (8)

where $\chi^{zz}({\bf\it
r},t):=-(L^{2}/i\hbar)\theta(t)\langle\quantity[s^{z}({\bf\it
r},t),s^{z}({\bf\it 0},0)]\rangle$ with the interface area $L^{2}$ and the
spin density operator in the interaction picture, $s^{z}({\bf\it
r},t)=e^{i(H_{\mathrm{M}}+H_{\mathrm{Z}})t/\hbar}s^{z}({\bf\it
r})e^{-i(H_{\mathrm{M}}+H_{\mathrm{Z}})t/\hbar}$. For the exchange coupling,
we consider configuration average and assume $\langle\sum_{n}J({\bf\it
r},{\bf\it r}_{n})\rangle_{\mathrm{ave}}=J_{1}$ , which means that
$H_{\mathrm{Z}}$ is treated as a uniform Zeeman like interaction and the
interface is specular 111In a realistic sample, correction terms would be
added due to interface roughness Ominato  _et al._ (2022a, b). Even with such
correction terms, the characteristic angular dependence is expected to be
preserved, since the Majorana Ising spin property of having spin density only
in the perpendicular direction is maintained.. Using eigenstates of Eq. (2)
and after a straightforward calculation, the uniform spin susceptibility is
given by

$\displaystyle\tilde{\chi}^{zz}(\mathbf{0},\omega)$
$\displaystyle=-\sum_{{\bf\it k},\lambda}\absolutevalue{\langle{\bf\it
k},\lambda|\sigma^{z}|{\bf\it k},-\lambda\rangle}^{2}\frac{f(E_{{\bf\it
k},\lambda})-f(E_{{\bf\it k},-\lambda})}{2E_{{\bf\it
k},\lambda}+\hbar\omega+i0},$ $\displaystyle\to-\int
dED(E)\frac{E^{2}-M^{2}}{2E^{2}}\frac{f(E)-f(-E)}{2E+\hbar\omega+i0},$ (9)

where $|{\bf\it k},\lambda\rangle$ is an eigenstate of $H_{\mathrm{M}}$ with
eigenenergy $E_{{\bf\it k},\lambda}=\lambda\sqrt{(\hbar vk)^{2}+M^{2}}$,
($\lambda=\pm$). $M=J_{1}S\cos\vartheta$ is the Majorana gap,
$f(E)=1/(e^{E/k_{\mathrm{B}}T}+1)$ is the Fermi distribution function, and
$D(E)$ is the density of states given by

$\displaystyle D(E)=\frac{L^{2}}{2\pi(\hbar
v)^{2}}\absolutevalue{E}\theta\quantity(\absolutevalue{E}-\absolutevalue{M}),$
(10)

with the Heaviside step function $\theta(x)$. It is important to note that the
behavior of the uniform spin susceptibility is determined by the interband
contribution, which is proportional to the Fermi distribution function, i.e.,
the contribution of the occupied states. This mechanism is similar to the Van
Vleck paramagnetism Neil W. and N. David (1976). The contribution of the
occupied states often plays a crucial role in topological responses Murakami
(2006).

Replacing the localized spin operators with their statistical average values,
we find the induced Majorana Ising spin density, to the first order of
$J_{1}S$, is given by

$\displaystyle\int d{\bf\it r}\expectationvalue{s^{z}({\bf\it r},t)}$
$\displaystyle=\tilde{\chi}^{zz}_{0}({\bf\it 0},0)J_{1}S\cos\vartheta$
$\displaystyle+\mathrm{Re}\quantity[\tilde{\chi}^{zz}_{0}({\bf\it
0},\omega)]\frac{h_{\mathrm{ac}}}{\alpha
h_{\mathrm{dc}}}J_{1}S\sin\vartheta\sin\omega t,$ (11)

where $\tilde{\chi}^{zz}_{0}({\bf\it 0},0)$ is the spin susceptibility for
$M=0$. The first term originates from $H_{\mathrm{Z}}$ and gives a static spin
density, while the second term originates from $H_{\mathrm{T}}$ and gives a
dynamic spin density. Figure 2 shows the induced Ising spin density as a
function of time at several angles. As shown in Eq. (11), the Ising spin
density consists of the static and dynamic components. The dynamic component
is induced by the precessional motion of the localized spin, which means one
can induce the DMISR using the dynamics of the localized spin.

The inset in Fig. 2 shows $\mathrm{Im}\tilde{\chi}^{zz}({\bf\it 0},\omega)$ as
a function of $\vartheta$ at a fixed frequency. When the frequency
$\hbar\omega$ is smaller than the Majorana gap,
$\mathrm{Im}\tilde{\chi}^{zz}({\bf\it 0},\omega)$ is zero. Once the frequency
overcomes the Majorana gap, $\mathrm{Im}\tilde{\chi}^{zz}({\bf\it 0},\omega)$
becomes finite. The implications of these behaviors are that if the magnon
energy is smaller than the Majorana gap, there is no energy dissipation due to
the DMISR. However, once the magnon energy exceeds the Majorana gap, finite
energy dissipation associated with the DMISR occurs at the surface of the TSC.
Therefore, one can toggle between dissipative and non-dissipative Majorana
Ising spin dynamics by adjusting the precession axis angle and the microwave
frequency.

Figure 2: The induced Ising spin density, with a unit
$\tilde{\chi}^{zz}_{0}({\bf\it 0},0)J_{1}S$, is presented as a function of
$\omega t$ and $\vartheta$. The frequency and temperature are set to
$\hbar\omega/J_{1}S=1.5$ and $k_{\mathrm{B}}T/J_{1}S=0.1$, respectively. The
coefficient, $h_{\mathrm{ac}}/\alpha h_{\mathrm{dc}}$, is set to $0.3$. The
static Majorana Ising spin density arises from $H_{\mathrm{Z}}$. When the
precession axis deviates from the direction perpendicular to the interface,
the precessional motion of the localized spins results in the dynamical
Majorana Ising spin response (DMISR). Energy dissipation due to the DMISR is
zero for small angles $\vartheta$ as the Majorana gap exceeds the magnon
energy. However, once the magnon energy overcomes the Majorana gap, the energy
dissipation becomes finite. Therefore, one can toggle between dissipative and
non-dissipative DMISR by adjusting $\vartheta$.

FMR modulation.—The retarded component of the magnon Green’s function is given
by $G^{R}({\bf\it
r}_{n},t)=-(i/\hbar)\theta(t)\langle[S^{+}_{n}(t),S^{-}_{0}(0)]\rangle$ with
the interaction picture
$S^{\pm}_{n}(t)=e^{iH_{\mathrm{FI}}t/\hbar}S^{\pm}_{n}e^{-iH_{\mathrm{FI}}t/\hbar}$.
The FMR signal is characterized by the spectral function defined as

$\displaystyle A({\bf\it
q},\omega):=-\frac{1}{\pi}\mathrm{Im}\quantity[\sum_{n}e^{-i{\bf\it
q}\cdot{\bf\it r}_{n}}\int dte^{i(\omega+i0)t}G^{R}({\bf\it r}_{n},t)].$ (12)

For uniform external force, the spectral function is given by

$\displaystyle
A(\mathbf{0},\omega)=\frac{2S}{\hbar}\frac{1}{\pi}\frac{(\alpha+\delta\alpha)\omega}{\quantity[\omega-\gamma(h_{\mathrm{dc}}+\delta
h)]^{2}+\quantity[(\alpha+\delta\alpha)\omega]^{2}}.$ (13)

The peak position and width of the FMR signal is given by
$h_{\mathrm{dc}}+\delta h$ and $\alpha+\delta\alpha$, respectively.
$h_{\mathrm{dc}}$ and $\alpha$ correspond to the peak position and the
linewidth of the FMR signal of the FI alone. $\delta h$ and $\delta\alpha$ are
the FMR modulations due to the exchange interaction $H_{\mathrm{T}}$. We treat
$H_{\mathrm{M}}+H_{\mathrm{FI}}+H_{\mathrm{Z}}$ as an unperturbed Hamiltonian
and $H_{\mathrm{T}}$ as a perturbation. In this work, we assume the specular
interface, where the coupling $J({\bf\it r},{\bf\it r}_{n})$ is approximated
as $\Big{\langle}\sum_{n,n^{\prime}}{J({\bf\it r},{\bf\it r}_{n})J({\bf\it
r}^{\prime},{\bf\it r}_{n^{\prime}})}\Big{\rangle}_{\mathrm{ave}}=J_{1}^{2}$.
The dynamics of the localized spins in the FI is modulated due to the
interaction between the localized spins and the Majorana Ising spins. In our
setup, the peak position and the linewidth of the FMR signal are modulated and
the FMR modulation is given by

$\displaystyle\delta
h=\sin^{2}\vartheta\frac{SJ^{2}_{1}}{2N\gamma\hbar}\mathrm{Re}\tilde{\chi}^{zz}({\bf\it
0},\omega),$ (14)
$\displaystyle\delta\alpha=\sin^{2}\vartheta\frac{SJ^{2}_{1}}{2N\hbar\omega}\mathrm{Im}\tilde{\chi}^{zz}({\bf\it
0},\omega),$ (15)

where $N$ is the total number of sites in the FI. These formulas were derived
in the study of the FMR in magnetic multilayer systems including
superconductors. One can extract the spin property of the Majorana mode from
the data on $\delta h$ and $\delta\alpha$. Because of the Ising spin
anisotropy, the FMR modulation exhibits strong anisotropy, where the FMR
modulation is proportional to $\sin^{2}\vartheta$.

Figure 3 shows the FMR modulations (a) $\delta\alpha$ and (b) $\delta h$. The
FMR modulation at a fixed frequency increases with angle $\vartheta$ and
reaches a maximum at $\pi/2$, as can be read from Eqs. (14) and (15). When the
angle $\vartheta$ is fixed and the frequency $\omega$ is increased,
$\delta\alpha$ becomes finite above a certain frequency at which the energy of
the magnon coincides with the Majorana gap. When $\vartheta<\pi/2$ and
$\hbar\omega\approx 2M$, $\delta\alpha$ linearly increases as a function of
$\omega$ just above the Majorana gap. The localized spin damping is enhanced
when the magnon energy exceeds the Majorana gap. At $\vartheta=\pi/2$ and
$\omega\approx 0$, the Majorana gap vanishes and $\delta\alpha$ is
proportional to $\omega/T$. In the high frequency region
$\hbar\omega/J_{1}S\gg 1$, $\delta\alpha$ converges to its upper threshold.

The frequency shift $\delta h$ is almost independent of $\omega$ and has a
finite value even in the Majorana gap. This behavior is analogous to the
interband contribution to the spin susceptibility in strongly spin-orbit
coupled band insulators, and is due to the fact that the effective Hamiltonian
of the Majorana modes includes spin operators. It is important to emphasize
that although the Majorana modes have spin degrees of freedom, only the $z$
component of the spin density operator is well defined. This is a hallmark of
Majorana modes, which differs significantly from electrons in ordinary solids.
Note that $\delta h$ is proportional to the energy cutoff, which is introduced
to converge energy integral for $\mathrm{Re}\tilde{\chi}^{zz}({\bf\it
0},\omega)$. The energy cutoff corresponds to the bulk superconducting gap,
which is estimated as $\Delta\sim 0.1[\mathrm{meV}]$ $(\sim 1[\mathrm{K}])$.
Therefore, our results are applicable in the frequency region below
$\hbar\omega\sim 0.1[\mathrm{meV}]$ $(\sim 30[\mathrm{GHz}])$. In addition, we
assume that Majorana gap is estimated to be $J_{1}S\sim 0.01[\mathrm{meV}]$
$(\sim 0.1[\mathrm{K}])$.

Figure 3:  The temperature is set to $k_{\mathrm{B}}T/J_{1}S=0.1$. (a) The
damping modulation $\delta\alpha$ only becomes finite when the magnon energy
exceeds the Majorana gap; otherwise, it vanishes. This behavior corresponds to
the energy dissipation of the Majorana Ising spin. (b) The peak shift is
finite, except for $\vartheta=0$, and is almost independent of $\omega$. This
behavior resembles the spin response observed in strongly spin-orbit coupled
band insulators, where the interband contribution to spin susceptibility
results in a finite spin response, even within the energy gap.

Discussion.—Comparing the present results with spin pumping (SP) in a
conventional metal-ferromagnet hybrid, the qualitative behaviors are quite
different. In conventional metals, spin accumulation occurs due to FMR. In
contrast, in the present system, no corresponding spin accumulation occurs due
to the Ising anisotropy. Also, in the present calculations, the proximity-
induced exchange coupling is assumed to be an isotropic Heisenberg-like
coupling. However, in general, the interface interaction can also be
anisotropic. Even in such a case, it is no qualitative change in the case of
ordinary metals, although a correction term due to anisotropy is added Matsuo
_et al._ (2018). Therefore, the Ising anisotropy discussed in the present work
is a property unique to the Majorana modes and can characterize the Majorana
excitations.

Let us comment on the universal nature of the toggling between non-dissipative
and dissipative dynamical spin responses observed in our study. Indeed, such
toggling becomes universally feasible when the microwave frequency and the
energy gap are comparable, and when the Hamiltonian and spin operators are
non-commutative, indicating that spin is not a conserved quantity. The non-
commutativity can be attributed to the presence of spin-orbit couplings Thakur
_et al._ (2018); Zhou and Chang (2018); Ominato and Nomura (2018), and spin-
triplet pair correlations Sigrist and Ueda (1991).

Microwave irradiation leads to heating within the FI, so that thermally
excited magnons due to the heating could influence the DMISR. Phenomena
resulting from the heating, which can affect interface spin dynamics, include
the spin Seebeck effect (SSE) Uchida  _et al._ (2008), where a spin current is
generated at the interface due to a temperature difference. In hybrid systems
of normal metal and FI, methods to separate the inverse spin Hall voltage due
to SP from other signals caused by heating have been well studied Iguchi and
Saitoh (2017). Especially, it has been theoretically proposed that SP and SSE
signals can be separated using a spin current noise measurement Matsuo  _et
al._ (2018). Moreover, SP coherently excites specific modes, which
qualitatively differs from SSE induced by thermally excited magnons Kato  _et
al._ (2019). Therefore, even if heating occurs in the FI in our setup, the
properties of Majorana Ising spins are expected to be captured. Details of the
heating effect on the DMISR will be examined in the near future.

We also mention the experimental feasibility of our theoretical proposals. As
we have already explained, the FMR modulation is a very sensitive spin probe.
Indeed, the FMR modulation by surface states of 3D topological insulators Tang
_et al._ (2018) and graphene Patra  _et al._ (2012); Tang  _et al._ (2013);
Mendes  _et al._ (2015); Indolese  _et al._ (2018); Mendes  _et al._ (2019)
has been reported experimentally. Therefore, we expect that the enhanced
Gilbert damping due to Majorana Ising spin can be observable in our setup when
the thickness of the ferromagnetic insulator is sufficiently thin.

Finally, it is pertinent to mention the potential candidate materials where
surface Majorana Ising spins could be detectable. Notably, $\mathrm{UTe}_{2}$
Ran  _et al._ (2019), $\mathrm{Cu}_{x}\mathrm{Bi}_{2}\mathrm{Se}_{3}$ Matano
_et al._ (2016); Yonezawa  _et al._ (2017), SrxBi2Se3 and NbxBi2Se3 Sharma
_et al._ (2022) are reported to be in a $p$-wave superconducting state and
theoretically can host surface Majorana Ising spins. Recent NMR measurements
indicate that $\mathrm{UTe}_{2}$ could be a bulk $p$-wave superconductor in
the Balian-Werthamer state Matsumura  _et al._ (2023), which hosts the surface
Majorana Ising spins with the perpendicular Ising anisotropy, as considered in
this work. $A_{x}\mathrm{Bi}_{2}\mathrm{Se}_{3}$ ($A=$ Cu, Sr, Nb) is
considered to possess in-plane Ising anisotropy Xiong  _et al._ (2017),
differing from the perpendicular Ising anisotropy explored in this work.
Therefore, we expect that it exhibits anisotropy different from that
demonstrated in this work.

Conclusion.—We present herein a study of the spin dynamics in a topological
superconductor (TSC)-magnet hybrid. Ferromagnetic resonance under microwave
irradiation leads to the dynamically induced Majorana Ising spin density on
the TSC surface. One can toggle between dissipative and non-dissipative
Majorana Ising spin dynamics by adjusting the external magnetic field angle
and the microwave frequency. Therefore, our setup provides a platform to
detect and control Majorana excitations. We expect that our results provide
insights toward the development of future quantum computing and spintronics
devices using Majorana excitations.

Acknowledgments.— The authors are grateful to R. Shindou for valuable
discussions. This work is partially supported by the Priority Program of
Chinese Academy of Sciences, Grant No. XDB28000000. We acknowledge JSPS
KAKENHI for Grants (Nos. JP20K03835, JP21H01800, JP21H04565, and JP23H01839).

## Supplemental Material

###  .1 Surface Majorana modes

In this section, we describe the procedure for deriving the effective
Hamiltonian of the surface Majorana modes. We start with the bulk Hamiltonian
of a three-dimensional topological superconductor. Based on the bulk
Hamiltonian, we solve the BdG equation to demonstrate the existence of a
surface-localized solution. Using this solution, we expand the field operator
and show that it satisfies the Majorana condition when the bulk excitations
are neglected. As a result, on energy scales much smaller than the bulk
superconducting gap, the low-energy excitations are described by surface-
localized Majorana modes. The above procedure is explained in more detail in
the following. Note that we use ${\bf\it r}$ for three-dimensional coordinates
and ${\bf\it r}_{\parallel}$ for two-dimensional ones in the Supplemental
Material.

We start with the mean-field Hamiltonian given by

$\displaystyle H_{\mathrm{SC}}=\frac{1}{2}\int d{\bf\it
r}{\bf\it\Psi}_{\mathrm{BdG}}^{\dagger}({\bf\it
r})H_{\mathrm{BdG}}{\bf\it\Psi}_{\mathrm{BdG}}({\bf\it r}),$ (16)

with ${\bf\it r}=(x,y,z)$. We consider the Balian-Werthamer (BW) state, in
which the pair potential is given by $\Delta_{\hat{{\bf\it
k}}}=\frac{\Delta}{k_{\mathrm{F}}}\quantity(\hat{{\bf\it
k}}\cdot{\bf\it\sigma})i\sigma^{y}$ with the bulk superconducting gap
$\Delta$. Here, we do not discuss the microscopic origin of the pair
correlation leading to the BW state. As a result, the BdG Hamiltonian
$H_{\mathrm{BdG}}$ is given by

with $\hat{k}_{\pm}=\hat{k}_{y}\pm i\hat{k}_{z}$, $\hat{{\bf\it k}}=-i\nabla$,
and $\varepsilon_{\hat{{\bf\it k}}}=\frac{\hbar^{2}\hat{{\bf\it k}}^{2}}{2m}$.
The four component Nambu spinor ${\bf\it\Psi}_{\mathrm{BdG}}({\bf\it r})$ is
given by

$\displaystyle{\bf\it\Psi}_{\mathrm{BdG}}({\bf\it
r}):=\matrixquantity(\Psi_{\to}({\bf\it r})\\\ \Psi_{\leftarrow}({\bf\it
r})\\\ \Psi_{\to}^{\dagger}({\bf\it r})\\\ \Psi_{\leftarrow}^{\dagger}({\bf\it
r})),$ (17)

with the spin quantization axis along the $x$ axis. The matrices of the spin
operators are represented as

The fermion field operators satisfy the anticommutation relations

$\displaystyle\\{\Psi_{\sigma}({\bf\it r}),\Psi_{\sigma^{\prime}}({\bf\it
r}^{\prime})\\}=0,$ (18) $\displaystyle\\{\Psi_{\sigma}({\bf\it
r}),\Psi_{\sigma^{\prime}}^{\dagger}({\bf\it
r}^{\prime})\\}=\delta_{\sigma\sigma^{\prime}}\delta({\bf\it r}-{\bf\it
r}^{\prime}),$ (19)

with the spin indices $\sigma,\sigma^{\prime}=\to,\leftarrow$.

To diagonalize the BdG Hamiltonian, we solve the BdG equation given by

$\displaystyle H_{\mathrm{BdG}}{\bf\it\Phi}({\bf\it r})=E{\bf\it\Phi}({\bf\it
r}).$ (20)

We assume that a solution is written as

$\displaystyle{\bf\it\Phi}({\bf\it r})=e^{i{\bf\it k}_{\parallel}\cdot{\bf\it
r}_{\parallel}}f(z)\matrixquantity(u_{\to}\\\ u_{\leftarrow}\\\ v_{\to}\\\
v_{\leftarrow}),$ (21)

with ${\bf\it k}_{\parallel}=(k_{x},k_{y})$ and ${\bf\it
r}_{\parallel}=(x,y)$. If we set the four components vector to satisfy the
following equation (Majorana condition)

$\displaystyle\matrixquantity(0&0&1&0\\\ 0&0&0&1\\\ 1&0&0&0\\\
0&1&0&0)\matrixquantity(u_{\to}\\\ u_{\leftarrow}\\\ v_{\to}\\\
v_{\leftarrow})=\pm\matrixquantity(u_{\to}\\\ u_{\leftarrow}\\\ v_{\to}\\\
v_{\leftarrow}),$ (22)

we can obtain a surface-localized solution. If we take a positive (negative)
sign, we obtain a solution localized on the top surface (bottom surface). As
we will consider solutions localized on the bottom surface below, we take a
negative sign. Finally, we obtain the normalized eigenvectors of the BdG
equation given by

$\displaystyle{\bf\it\Phi}_{\lambda,{\bf\it k}_{\parallel}}({\bf\it
r})=\frac{e^{i{\bf\it k}_{\parallel}\cdot{\bf\it
r}_{\parallel}}}{\sqrt{L^{2}}}f_{{\bf\it k}_{\parallel}}(z){\bf\it
u}_{\lambda,{\bf\it k}_{\parallel}},$ (23)

with

$\displaystyle f_{{\bf\it k}_{\parallel}}(z)=N_{{\bf\it
k}_{\parallel}}\sin(k_{\perp}z)e^{-\kappa z},$ (24) $\displaystyle N_{{\bf\it
k}_{\parallel}}=\sqrt{\frac{4\kappa(k_{\perp}^{2}+\kappa^{2})}{k_{\perp}^{2}}},$
(25) $\displaystyle\kappa=\frac{m\Delta}{\hbar^{2}k_{\mathrm{F}}},$ (26)
$\displaystyle
k_{\perp}=\sqrt{k_{\mathrm{F}}^{2}-{k_{\parallel}^{2}-\kappa^{2}}},$ (27)

and

$\displaystyle{\bf\it u}_{+,{\bf\it
k}_{\parallel}}=\matrixquantity(u_{+,\to{\bf\it k}_{\parallel}}\\\
u_{+,\leftarrow{\bf\it k}_{\parallel}}\\\ v_{+,\to{\bf\it k}_{\parallel}}\\\
v_{+,\leftarrow{\bf\it
k}_{\parallel}})=\frac{1}{\sqrt{2}}\matrixquantity(\sin\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}\\\ -\cos\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}\\\ -\sin\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}\\\ \cos\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}),$ (28) $\displaystyle{\bf\it u}_{-,{\bf\it
k}_{\parallel}}=\matrixquantity(u_{-,\to{\bf\it k}_{\parallel}}\\\
u_{-,\leftarrow{\bf\it k}_{\parallel}}\\\ v_{-,\to{\bf\it k}_{\parallel}}\\\
v_{-,\leftarrow{\bf\it
k}_{\parallel}})=\frac{1}{\sqrt{2}}\matrixquantity(-\cos\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}\\\ -\sin\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}\\\ \cos\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}\\\ \sin\frac{\phi_{{\bf\it
k}_{\parallel}}+\pi/2}{2}).$ (29)

The eigenenergy is given by $E_{\lambda,{\bf\it k}_{\parallel}}=\lambda\Delta
k_{\parallel}/k_{\mathrm{F}}$. We can show that the eigenvectors satisfy

$\displaystyle{\bf\it u}_{-,-{\bf\it k}_{\parallel}}={\bf\it u}_{+,{\bf\it
k}_{\parallel}}.$ (30)

Consequently, the field operator is expanded as

$\displaystyle{\bf\it\Psi}_{\mathrm{BdG}}({\bf\it r})$
$\displaystyle=\sum_{{\bf\it k}_{\parallel}}\quantity(\gamma_{{\bf\it
k}_{\parallel}}\frac{e^{i{\bf\it k}_{\parallel}\cdot{\bf\it
r}_{\parallel}}}{\sqrt{L^{2}}}+\gamma_{{\bf\it
k}_{\parallel}}^{\dagger}\frac{e^{-i{\bf\it k}_{\parallel}\cdot{\bf\it
r}_{\parallel}}}{\sqrt{L^{2}}})$ $\displaystyle\hskip 28.45274pt\times
f_{{\bf\it k}_{\parallel}}(z){\bf\it u}_{+,{\bf\it
k}_{\parallel}}+(\mathrm{bulk}{~{}}\mathrm{modes}),$ (31)

where $\gamma_{{\bf\it k}_{\parallel}}$ ($\gamma_{{\bf\it
k}_{\parallel}}^{\dagger}$) is the quasiparticle creation (annihilation)
operator with the eigenenergy $E_{+,{\bf\it k}_{\parallel}}$. Substituting the
above expression into Eq. (16) with omission of bulk modes and performing the
integration in the $z$-direction, we obtain the effective Hamiltonian for the
surface states

$\displaystyle H_{\mathrm{M}}=\frac{1}{2}\int d{\bf\it
r}_{\parallel}{\bf\it\psi}^{\mathrm{T}}({\bf\it r}_{\parallel})\quantity(\hbar
v\hat{k}_{y}\sigma^{x}-\hbar v\hat{k}_{x}\sigma^{y}){\bf\it\psi}({\bf\it
r}_{\parallel}),$ (32)

where $v=\Delta/\hbar k_{\mathrm{F}}$ and we introduced the two component
Majorana field operator

$\displaystyle{\bf\it\psi}({\bf\it
r}_{\parallel})=\matrixquantity(\psi_{\to}({\bf\it r}_{\parallel})\\\
\psi_{\leftarrow}({\bf\it r}_{\parallel})),$ (33)

satisfying the Majorana condition

$\displaystyle\psi_{\sigma}({\bf\it
r}_{\parallel})=\psi_{\sigma}^{\dagger}({\bf\it r}_{\parallel}),$ (34)

and the anticommutation relation

$\displaystyle\quantity{\psi_{\sigma}({\bf\it
r}_{\parallel}),\psi_{\sigma^{\prime}}({\bf\it
r}_{\parallel}^{\prime})}=\delta_{\sigma\sigma^{\prime}}\delta({\bf\it
r}_{\parallel}-{\bf\it r}_{\parallel}^{\prime}).$ (35)

The spin density operator of the Majorana mode is given by

$\displaystyle{\bf\it s}({\bf\it
r}_{\parallel})={\bf\it\psi}^{\dagger}({\bf\it
r}_{\parallel})\frac{{\bf\it\sigma}}{2}{\bf\it\psi}({\bf\it r}_{\parallel}).$
(36)

The $x$ component is given by

(37)

In a similar manner, the $y$ and $z$ components are given by

(38)

and

(39)

respectively. As a result, the spin density operator is given by

$\displaystyle{\bf\it s}({\bf\it
r}_{\parallel})=\quantity(0,0,-i\psi_{\to}({\bf\it
r}_{\parallel})\psi_{\leftarrow}({\bf\it r}_{\parallel})).$ (40)

One can see that the spin density of the Majorana mode is Ising like.

###  .2 Majorana Ising spin dynamics

In this section, we calculate the Ising spin density induced on the TSC
surface by the proximity coupling $H_{\mathrm{ex}}$. $H_{\mathrm{ex}}$
consists of two terms, $H_{\mathrm{Z}}$ and $H_{\mathrm{T}}$. $H_{\mathrm{Z}}$
leads to the static spin density and $H_{\mathrm{T}}$ leads to the dynamic
spin density. First, we calculate the static spin density. Next, we calculate
the dynamic spin density.

The total spin density operator is given by

$\displaystyle s^{z}_{\mathrm{tot}}$ $\displaystyle=\int d{\bf\it
r}_{\parallel}s^{z}({\bf\it r}_{\parallel}).$ (41)

The statistical average of the static spin density is calculated as

$\displaystyle\expectationvalue{s^{z}_{\mathrm{tot}}}$
$\displaystyle=-\sum_{{\bf\it k}_{\parallel}}\frac{M}{2E_{{\bf\it
k}_{\parallel}}}\quantity[f(E_{{\bf\it k}_{\parallel}})-f(-E_{{\bf\it
k}_{\parallel}})]$ $\displaystyle\to-\quantity(\frac{L}{2\pi\hbar
v})^{2}\int_{M}^{\Delta}EdE\int_{0}^{2\pi}d\phi\frac{M}{2E}\quantity[f(E)-f(-E)]$
$\displaystyle=-\int_{0}^{\Delta}dED(E)\frac{f(E)-f(-E)}{2E}M.$ (42)

At the zero temperature limit $T\to 0$, the static spin density is given by

$\displaystyle\expectationvalue{s^{z}_{\mathrm{tot}}}$
$\displaystyle=\frac{1}{2}\frac{L^{2}}{2\pi(\hbar
v)^{2}}\quantity(\Delta-M)M\approx\tilde{\chi}^{zz}_{0}({\bf\it 0},0)M,$ (43)

where $\tilde{\chi}^{zz}_{0}({\bf\it 0},0)=D(\Delta)/2$ and we used $\Delta\gg
M$.

The dynamic spin density is given by the perturbative force

$\displaystyle H_{\mathrm{T}}(t)=\int d{\bf\it r}_{\parallel}s^{z}({\bf\it
r}_{\parallel})F({\bf\it r}_{\parallel},t),$ (44)

where $F({\bf\it r}_{\parallel},t)$ is given by

$\displaystyle F({\bf\it r}_{\parallel},t)$
$\displaystyle=-\sin\vartheta\sum_{n}J({\bf\it r}_{\parallel},{\bf\it
r}_{n})\expectationvalue{S_{n}^{X}(t)}$ $\displaystyle\approx-\sin\vartheta
J_{1}S\frac{\gamma h_{\mathrm{ac}}}{\sqrt{(\omega-\gamma
h_{\mathrm{dc}})^{2}+\alpha^{2}\omega^{2}}}\cos\omega t$
$\displaystyle=:F\cos\omega t.$ (45)

The time dependent statistical average of the Ising spin density, to the first
order of $J_{1}S$, is given by

$\displaystyle\int d{\bf\it r}_{\parallel}\expectationvalue{s^{z}({\bf\it
r}_{\parallel},t)}$ $\displaystyle=\int d{\bf\it r}_{\parallel}\int d{\bf\it
r}_{\parallel}^{\prime}\int dt^{\prime}\chi^{zz}({\bf\it
r}_{\parallel}-{\bf\it r}_{\parallel}^{\prime},t^{\prime})F({\bf\it
r}_{\parallel}^{\prime},t-t^{\prime})$
$\displaystyle=\mathrm{Re}\quantity[\tilde{\chi}^{zz}({\bf\it
0},\omega)Fe^{-i\omega t}]$
$\displaystyle\approx\mathrm{Re}\quantity[\tilde{\chi}^{zz}_{0}({\bf\it
0},\omega)]F\cos\omega t,$ (46)

where we used $\mathrm{Re}\tilde{\chi}^{zz}_{0}({\bf\it
0},\omega)\gg\mathrm{Im}\tilde{\chi}^{zz}_{0}({\bf\it 0},\omega)$. The real
part of $\tilde{\chi}^{zz}({\bf\it 0},\omega)$ is given by

$\displaystyle\mathrm{Re}\tilde{\chi}^{zz}({\bf\it 0},\omega)$
$\displaystyle=-\mathcal{P}\int
dED(E)\frac{E^{2}-M^{2}}{2E^{2}}\frac{f(E)-f(-E)}{2E+\hbar\omega},$ (47)

where $\mathcal{P}$ means the principal value. When the integrand is expanded
with respect to $\omega$, the lowest order correction term becomes quadratic
in $\omega$. In the frequency range considered in this work, this correction
term is significantly smaller compared to the static spin susceptibility
$\mathrm{Re}\tilde{\chi}^{zz}({\bf\it 0},0)$. Therefore, the spin
susceptibility exhibits almost no frequency dependence and remains constant as
a function of $\omega$. The imaginary part of $\tilde{\chi}^{zz}({\bf\it
0},\omega)$ is given by

$\displaystyle\mathrm{Im}\tilde{\chi}^{zz}({\bf\it 0},\omega)$
$\displaystyle=\pi
D(\hbar\omega/2)\frac{(\hbar\omega/2)^{2}-M^{2}}{2(\hbar\omega/2)^{2}}\quantity[f(-\hbar\omega/2)-f(\hbar\omega/2)].$
(48)

###  .3 FMR modulation due to the proximity exchange coupling

In this section, we provide a brief explanation for the derivation of the FMR
modulations $\delta h$ and $\delta\alpha$. The FMR modulations can be
determined from the retarded component of the magnon Green’s function, which
is given by

$\displaystyle\tilde{G}^{R}({\bf\it
k},\omega)=\frac{2S/\hbar}{\omega-\omega_{{\bf\it
k}}+i\alpha\omega-(2S/\hbar)\Sigma^{R}({\bf\it k},\omega)},$ (49)

where we introduce the Gilbert damping constant $\alpha$ phenomenologically.
In the second-order perturbation calculation with respect to $H_{\mathrm{T}}$,
the self-energy is given by

$\displaystyle\Sigma^{R}({\bf\it
k},\omega)=-\quantity(\frac{\sin\vartheta}{2})^{2}\sum_{{\bf\it
q}_{\parallel}}|\tilde{J}({\bf\it q}_{\parallel},{\bf\it
k})|^{2}\tilde{\chi}^{zz}({\bf\it q}_{\parallel},\omega),$ (50)

where $\tilde{J}({\bf\it q}_{\parallel},{\bf\it 0})$ is given by

$\displaystyle\tilde{J}({\bf\it q}_{\parallel},{\bf\it
k})=\frac{1}{L^{2}\sqrt{N}}\int d{\bf\it r}_{\parallel}\sum_{n}J({\bf\it
r}_{\parallel},{\bf\it r}_{n})e^{i({\bf\it q}_{\parallel}\cdot{\bf\it
r}_{\parallel}+{\bf\it k}\cdot{\bf\it r}_{n})}$ (51)

The pole of $\tilde{G}^{R}({\bf\it k},\omega)$ signifies the FMR modulations,
including both the frequency shift and the enhanced Gilbert damping. These are
given by

$\displaystyle\delta h=\frac{2S}{\gamma\hbar}\mathrm{Re}\Sigma^{R}({\bf\it
0},\omega),\hskip
5.69054pt\delta\alpha=-\frac{2S}{\hbar\omega}\mathrm{Im}\Sigma^{R}({\bf\it
0},\omega).$ (52)

From the above equations and Eq. (50), it is apparent that FMR modulations
provide information regarding both the properties of the interface coupling
and the dynamic spin susceptibility of the Majorana modes.

The form of matrix element $\tilde{J}({\bf\it q}_{\parallel},{\bf\it 0})$
depends on the details of the interface. In this work, we assume the specular
interface. $|\tilde{J}({\bf\it q}_{\parallel},{\bf\it 0})|^{2}$ is given by

$\displaystyle|\tilde{J}({\bf\it q}_{\parallel},{\bf\it
0})|^{2}=\frac{J_{1}^{2}}{N}\delta_{{\bf\it q}_{\parallel},{\bf\it 0}}.$ (53)

Using Eq. (53), the self-energy for the uniform magnon mode is given by

$\displaystyle\Sigma^{R}({\bf\it
0},\omega)=-\quantity(\frac{\sin\vartheta}{2})^{2}\frac{J_{1}^{2}}{N}\tilde{\chi}^{zz}({\bf\it
0},\omega).$ (54)

## References

  * Wilczek (2009) F. Wilczek, Nat. Phys. 5, 614 (2009).
  * Yazdani  _et al._ (2023) A. Yazdani, F. von Oppen, B. I. Halperin, and A. Yacoby, Science 380, eade0850 (2023).
  * Chung and Zhang (2009) S. B. Chung and S.-C. Zhang, Phys. Rev. Lett. 103, 235301 (2009).
  * Sato and Fujimoto (2009) M. Sato and S. Fujimoto, Phys. Rev. B 79, 094504 (2009).
  * Nagato  _et al._ (2009) Y. Nagato, S. Higashitani, and K. Nagai, J. Phys. Soc. Jpn. 78, 123603 (2009).
  * Shindou  _et al._ (2010) R. Shindou, A. Furusaki, and N. Nagaosa, Phys. Rev. B 82, 180505 (2010).
  * Shiozaki and Sato (2014) K. Shiozaki and M. Sato, Phys. Rev. B 90, 165114 (2014).
  * Xiong  _et al._ (2017) Y. Xiong, A. Yamakage, S. Kobayashi, M. Sato, and Y. Tanaka, Crystals 7, 58 (2017).
  * Qiu  _et al._ (2016) Z. Qiu, J. Li, D. Hou, E. Arenholz, A. T. N’Diaye, A. Tan, K.-i. Uchida, K. Sato, S. Okamoto, Y. Tserkovnyak, Z. Q. Qiu, and E. Saitoh, Nat. Commun. 7, 12670 (2016).
  * Han  _et al._ (2020) W. Han, S. Maekawa, and X.-C. Xie, Nat. Mater. 19, 139 (2020).
  * Tserkovnyak  _et al._ (2002) Y. Tserkovnyak, A. Brataas, and G. E. Bauer, Phys. Rev. Lett. 88, 117601 (2002).
  * Tserkovnyak  _et al._ (2005) Y. Tserkovnyak, A. Brataas, G. E. Bauer, and B. I. Halperin, Reviews of Modern Physics 77, 1375 (2005).
  * Inoue  _et al._ (2017) M. Inoue, M. Ichioka, and H. Adachi, Phys. Rev. B 96, 024414 (2017).
  * Kato  _et al._ (2019) T. Kato, Y. Ohnuma, M. Matsuo, J. Rech, T. Jonckheere, and T. Martin, Phys. Rev. B 99, 144411 (2019).
  * Silaev (2020a) M. A. Silaev, Phys. Rev. B 102, 144521 (2020a).
  * Silaev (2020b) M. A. Silaev, Phys. Rev. B 102, 180502 (2020b).
  * Ominato  _et al._ (2022a) Y. Ominato, A. Yamakage, T. Kato, and M. Matsuo, Phys. Rev. B 105, 205406 (2022a).
  * Ominato  _et al._ (2022b) Y. Ominato, A. Yamakage, and M. Matsuo, Phys. Rev. B 106, L161406 (2022b).
  * (19) See Supplemental Material for the derivation of the effective Hamiltonian of surface Majorana modes, the derivation of the expression for Majorana Ising spin density, the calculation of dynamic spin density and dynamic spin susceptibility, and the formulation of ferromagnetic resonance modulation due to proximity exchange coupling.
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Physical Review 58, 1098 (1940).
  * Note (1) In a realistic sample, correction terms would be added due to interface roughness Ominato  _et al._ (2022a, b). Even with such correction terms, the characteristic angular dependence is expected to be preserved, since the Majorana Ising spin property of having spin density only in the perpendicular direction is maintained.
  * Neil W. and N. David (1976) A. Neil W. and M. N. David,  _Solid State Physics_ (Saunders College Publishing, 1976).
  * Murakami (2006) S. Murakami, Phys. Rev. Lett. 97, 236805 (2006).
  * Matsuo  _et al._ (2018) M. Matsuo, Y. Ohnuma, T. Kato, and S. Maekawa, Phys. Rev. Lett. 120, 037201 (2018).
  * Thakur  _et al._ (2018) A. Thakur, K. Sadhukhan, and A. Agarwal, Phys. Rev. B 97, 035403 (2018).
  * Zhou and Chang (2018) J. Zhou and H.-R. Chang, Phys. Rev. B 97, 075202 (2018).
  * Ominato and Nomura (2018) Y. Ominato and K. Nomura, Phys. Rev. B Condens. Matter 97, 245207 (2018).
  * Sigrist and Ueda (1991) M. Sigrist and K. Ueda, Reviews of Modern physics 63, 239 (1991).
  * Uchida  _et al._ (2008) K. Uchida, S. Takahashi, K. Harii, J. Ieda, W. Koshibae, K. Ando, S. Maekawa, and E. Saitoh, Nature 455, 778 (2008).
  * Iguchi and Saitoh (2017) R. Iguchi and E. Saitoh, journal of the physical society of japan 86, 011003 (2017).
  * Tang  _et al._ (2018) C. Tang, Q. Song, C.-Z. Chang, Y. Xu, Y. Ohnuma, M. Matsuo, Y. Liu, W. Yuan, Y. Yao, J. S. Moodera, S. Maekawa, W. Han, and J. Shi, Sci. Adv. 4, eaas8660 (2018).
  * Patra  _et al._ (2012) A. K. Patra, S. Singh, B. Barin, Y. Lee, J.-H. Ahn, E. del Barco, E. R. Mucciolo, and B. Özyilmaz, Appl. Phys. Lett. 101, 162407 (2012).
  * Tang  _et al._ (2013) Z. Tang, E. Shikoh, H. Ago, K. Kawahara, Y. Ando, T. Shinjo, and M. Shiraishi, Phys. Rev. B Condens. Matter 87, 140401 (2013).
  * Mendes  _et al._ (2015) J. B. S. Mendes, O. Alves Santos, L. M. Meireles, R. G. Lacerda, L. H. Vilela-Leão, F. L. A. Machado, R. L. Rodríguez-Suárez, A. Azevedo, and S. M. Rezende, Phys. Rev. Lett. 115, 226601 (2015).
  * Indolese  _et al._ (2018) D. Indolese, S. Zihlmann, P. Makk, C. Jünger, K. Thodkar, and C. Schönenberger, Phys. Rev. Appl. 10, 044053 (2018).
  * Mendes  _et al._ (2019) J. B. S. Mendes, O. Alves Santos, T. Chagas, R. Magalhães Paniago, T. J. A. Mori, J. Holanda, L. M. Meireles, R. G. Lacerda, A. Azevedo, and S. M. Rezende, Phys. Rev. B 99, 214446 (2019).
  * Ran  _et al._ (2019) S. Ran, C. Eckberg, Q.-P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I.-L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Science 365, 684 (2019).
  * Matano  _et al._ (2016) K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-Q. Zheng, Nat. Phys. 12, 852 (2016).
  * Yonezawa  _et al._ (2017) S. Yonezawa, K. Tajiri, S. Nakata, Y. Nagai, Z. Wang, K. Segawa, Y. Ando, and Y. Maeno, Nat. Phys. 13, 123 (2017).
  * Sharma  _et al._ (2022) M. M. Sharma, P. Sharma, N. K. Karn, and V. P. S. Awana, Supercond. Sci. Technol. 35, 083003 (2022).
  * Matsumura  _et al._ (2023) H. Matsumura, H. Fujibayashi, K. Kinjo, S. Kitagawa, K. Ishida, Y. Tokunaga, H. Sakai, S. Kambe, A. Nakamura, Y. Shimizu, Y. Homma, D. Li, F. Honda, and D. Aoki, J. Phys. Soc. Jpn. 92, 063701 (2023).


# On the Relation Between Autoencoders and Non-negative Matrix Factorization,
and Their Application for Mutational Signature Extraction

Ida Egendal  Center for Clinical Data Science, Aalborg University and Aalborg
University Hospital, Aalborg, Denmark  Clinical Cancer Research Center,
Aalborg University Hospital, Aalborg, Denmark  Rasmus Froberg Brøndum  Center
for Clinical Data Science, Aalborg University and Aalborg University Hospital,
Aalborg, Denmark  Clinical Cancer Research Center, Aalborg University
Hospital, Aalborg, Denmark  Marta Pelizzola  Department of Mathematics, Aarhus
University, Aarhus, Denmark  Asger Hobolth  Department of Mathematics, Aarhus
University, Aarhus, Denmark  Martin Bøgsted  Center for Clinical Data Science,
Aalborg University and Aalborg University Hospital, Aalborg, Denmark  Clinical
Cancer Research Center, Aalborg University Hospital, Aalborg, Denmark

###### Abstract

The aim of this study is to provide a foundation to understand the
relationship between non-negative matrix factorization (NMF) and non-negative
autoencoders enabling proper interpretation and understanding of autoencoder-
based alternatives to NMF. Since its introduction, NMF has been a popular tool
for extracting interpretable, low-dimensional representations of high-
dimensional data. However, recently, several studies have proposed to replace
NMF with autoencoders. This increasing popularity of autoencoders warrants an
investigation on whether this replacement is in general valid and reasonable.
Moreover, the exact relationship between non-negative autoencoders and NMF has
not been thoroughly explored. Thus, a main aim of this study is to investigate
in detail the relationship between non-negative autoencoders and NMF. We find
that the connection between the two models can be established through convex
NMF, which is a restricted case of NMF. In particular, convex NMF is a special
case of an autoencoder. The performance of NMF and autoencoders is compared
within the context of extraction of mutational signatures from cancer genomics
data. We find that the reconstructions based on NMF are more accurate compared
to autoencoders, while the signatures extracted using both methods show
comparable consistencies and values when externally validated. These findings
suggest that the non-negative autoencoders investigated in this article do not
provide an improvement of NMF in the field of mutational signature extraction.

##  1 Introduction

Non-negative matrix factorization (NMF) is a popular tool for unsupervised
learning [[15]]. NMF factorizes a non-negative data matrix into a product of
two non-negative matrices of lower dimension: a basis matrix consisting of
basis vectors and a weight matrix consisting of the basis vector’s weights for
each observation in the data matrix. NMF has gained a strong footing in
different scientific fields due to its high interpretability [[1, 6, 28]].
Specifically, NMF has proven to be a useful tool to derive mutational
signatures from cancer genomics data.

In mutational signature analysis, it is typically assumed that all mutations
in a cancer genome are caused by mutagenic processes that leave a
characteristic pattern of mutations in the genome. These patterns are denoted
mutational signatures. Several signatures have been identified and linked to
different mutagenic processes such as ultraviolet light exposure and tobacco
smoking [[17]]. [1] proposed using NMF on mutational count data from cancer
genomes to decipher the mutational signatures of the processes the patients
have been exposed to throughout the development of the disease. NMF has since
then been the dominating model for mutational signature extraction [[2, 3,
8]]. When extracting mutational signatures with NMF, the data matrix
consisting of a number of patients’ mutational profiles is decomposed into a
matrix representing the signatures of mutagenic processes (basis vectors) and
an exposure matrix dictating the number of mutations that can be attributed to
each specific process in the mutational profiles of each patient (weight
matrix).

Recently, several studies have proposed substituting NMF with non-negative
autoencoders which are increasingly popular for dimensionality reduction [[7,
28, 10, 22, 16]]. This is also the case in mutational signature extraction
[[20, 18]]. [20] suggested using a sparse autoencoder to identify mutational
signatures from cancer genomics data and generated estimates that were not
only in concordance with existing literature but also correlated with observed
exogenous exposures in a meaningful way. [18] suggested a hybrid architecture
with a deep encoding and shallow decoding to relax the assumption of linearity
NMF imposes on mutational signature extraction. However, they did not compare
the results and performance to NMF. This trend of using non-negative
autoencoders as an alternative to NMF with promising results prompts the
question: What is the mathematical relation between autoencoders and NMF? And
how do they compare in applications?

The aim of this study is to compare the performance of shallow, non-negative
autoencoders to NMF in the field of mutational signatures. In particular, we
compare the in- and out-of-sample reconstruction error, as well as the
consistency of the extracted signatures from the tumor-normal whole genome
sequences of 713 ovary, 311 prostate, and 523 uterus tumors from the Genomics
England 100,000 Genomes Project (GEL) [[25, 26]]. Moreover, we show
theoretically that shallow, non-negative autoencoders and NMF is a special
case of NMF where the basis vectors are restricted to be convex combinations
of columns in the data matrix [[5]]. Based on this, we deduce how it impacts
interpretation of the estimates generated by non-negative autoencoders in
general and especially within the field of mutational signatures.

In Section 2, we characterize the mathematical relationship between
autoencoders and NMF and introduce the framework for comparing the two models.
In Section 3, we compare the performance of NMF and autoencoders and
demonstrate the mathematical equivalence between convex NMF and non-negative,
shallow autoencoders. Lastly, we discuss and conclude on the results in
Sections 4 and 5.

##  2 Methods

Consider a non-negative data matrix $\bm{V}\in\mathbb{R}_{+}^{M\times N}$. In
this study the aim is to decompose $\bm{V}$ into a basis matrix
$\bm{H}\in\mathbb{R}_{+}^{M\times K}$ and a weight matrix
$\bm{W}\in\mathbb{R}_{+}^{K\times N}$, where $M$ denotes the number of
features, $N$ denotes the number of observations, and $K$ denotes the number
of basis vectors in the latent representation of $\bm{V}$. A schematic
overview of all considered decompositions in this section can be seen in
Figure 1.

###  2.1 Non-negative matrix factorization

NMF decomposes a matrix with non-negative entries into a matrix product of two
factor matrices with non-negative entries, one containing a set of basis
vectors and one containing a set of weights. The shared dimension, $K$, of the
factor matrices, is typically chosen to be much smaller than the dimensions of
the input matrix, making NMF a dimensionality reduction technique.

Standard $K$-dimensional NMF was introduced by [15] and aims to make a
reconstruction, $\bm{\hat{V}}$, of the original data matrix by a product of
two non-negative matrices:

$\bm{\hat{V}}=\bm{HW},$ (1)

where each column in $\bm{H}$ represents a basis vector and each column in
$\bm{W}$ represents each sample’s weights when being reconstructed as a linear
mixture of the basis vectors, i.e.,

$\bm{\hat{v}}_{n}=\sum_{k=1}^{K}\bm{h}_{k}w_{k,n},\quad n=1,\dots,N,$ (2)

where $\bm{\hat{v}}_{n}$ is the $n$’th column of the reconstructed data
matrix, $\hat{\bm{V}}$, $w_{n,k}$ represents the $(k,n)$’th entry of the
weight matrix $\bm{W}$, and $\bm{h}_{k}$ represents the $k$’th column of the
basis matrix $\bm{H}$ [[15]].

###  2.2 Convex NMF

Convex NMF, as introduced by [5], is a special case of NMF, where the basis
vectors are constrained to be spanned by the columns of the data matrix
$\bm{V}$, thus the data matrix is approximated by:

$\bm{\hat{V}}=\bm{V}\bm{W}_{1}\bm{W}_{2},$ (3)

with $\bm{W}_{1},\bm{W}_{2}^{T}\in\mathbb{R}_{+}^{N\times K}$. Defining
$\bm{H}:=\bm{VW}_{1}$ and $\bm{W}:=\bm{W}_{2}$ gives the model formulation of
NMF as defined in Equation (1). [5] focus mainly on the general case where the
data matrix, $\bm{V}$, can assume values within all real numbers, but for this
paper we consider a version of convex NMF where $\bm{V}$ is constrained to be
non-negative.

###  2.3 Autoencoders

Autoencoders consist of an encoder applied to the input to create a latent
representation and a decoder that maps the latent representation to a
reconstruction of the input [[12]]. Choosing the dimension of the latent
representation to be lower than the dimension of the input makes the
autoencoder a dimensionality reduction technique. A single hidden layer and
fully connected autoencoder’s reconstruction, $\hat{\bm{V}}$, of a data
matrix, $\bm{V}$, is mathematically defined as:

$\hat{\bm{V}}=\phi_{\text{dec}}(\phi_{\text{enc}}(\bm{V}\bm{W}_{\text{enc}}+\bm{b}_{\text{enc}})\bm{W}_{\text{dec}}+\bm{b}_{\text{dec}}),$
(4)

where $\bm{W}_{\text{enc}},\bm{W}_{\text{dec}}^{T}\in\mathbb{R}^{N\times K}$
are the encoding and decoding matrices, $\bm{b}_{\text{enc}}\in\mathbb{R}^{K}$
and $\bm{b}_{\text{dec}}\in\mathbb{R}^{N}$ are the bias terms of the encoding
and the decoding layers and $\phi_{\text{enc}}:\mathbb{R}^{M\times
K}\mapsto\mathbb{R}^{M\times K}$ and $\phi_{\text{dec}}:\mathbb{R}^{M\times
N}\mapsto\mathbb{R}^{M\times N}$ are the activation functions which provide
entry-wise modifications of the affected nodes.

###  2.4 Mathematical Equivalence and Interpretation

Figure 1: Schematic representation of the composition of basis vectors and
weights in NMF (top) and C-NMF and AE-NMF (bottom).

Setting $b_{\text{enc}}=\bm{0}_{K}$, $b_{\text{dec}}=\bm{0}_{N}$, and
$\phi_{\text{enc}},\phi_{\text{dec}}:x\mapsto x$ while constraining the
weights to be non-negative in Equation (4) yields exactly the convex NMF
formulation from Equation (3), where $\bm{W}_{1}$ corresponds to the encoding
matrix and $\bm{W}_{2}$ corresponds to the decoding matrix. Thus, non-negative
autoencoders can be constructed to be mathematically equivalent to convex NMF.
The choice between convex NMF and the autoencoder defined above therefore
reduces to a choice of how to optimize the problem, either by the
multiplicative updating steps derived by [5] or by the gradient descent-based
additive updates of the autoencoder. The core architecture of this autoencoder
is analogous to that of [20], but differ by fixing $b_{\text{enc}}$ and
$b_{\text{dec}}$ to zero instead of estimating them through training and
choosing the identity function as activation functions instead of
$\phi_{\text{enc}}:A\mapsto ReLU(x)=\text{max}(0,A)$ and
$\phi_{\text{dec}}:A\mapsto
Softmax(A)=\\{\exp(a_{m,n})/\sum_{m=1}^{M}\exp(a_{m,n})\\}_{m=1,\dots,M,n=1,\dots,N}$,
for a matrix $\bm{A}\in\mathbb{R}^{M\times N}$.

We will use ’C-NMF’ to refer to convex NMF optimized with the multiplicative
updating steps derived by [5] and use ’AE-NMF’ for the class of autoencoders
constructed equivalently to convex NMF.

The encoded data matrix in AE-NMF and C-NMF, $\bm{VW}_{1}$, which is a convex
combination of the columns of the data matrix, is interpreted as the basis
matrix, $\bm{H}$, in conventional NMF and hereby the interpretation of
$\bm{W}_{2}$ corresponds to that of $\bm{W}$ in conventional NMF.
Specifically, within mutational signature analysis this means that the
signature matrix is a convex combination of the patients’ mutational profiles
which aligns with how mutational signatures are understood. If the data matrix
is transposed to follow conventional orientation in neural networks, i.e.,
with patients as rows and mutation types as columns the weights or exposures
are modelled as convex combinations of the mutations which lacks an equally
straightforward interpretation and direct equivalence with convex NMF. The
equivalence between C-NMF and AE-NMF is the necessary and previously missing
link that enables one to interpret the parameters in AE-NMF similarly to NMF,
therefore this orientation is crucial for proper comparison.

Though the interpretation of AE-NMF, C-NMF, and NMF is similar, there are
still considerable differences between C-NMF, AE-NMF, and standard NMF. In
particular, NMF estimates $N\cdot K+K\cdot M$ parameters whereas AE-NMF and
C-NMF estimates $2\cdot(K\cdot N)$ parameters. Thus, AE-NMF and C-NMF will
estimate a larger number of parameters in the factor matrices than NMF when
the number of observations $N$ surpasses the number of features $M$, which is
often the case within mutational signatures.

###  2.5 Modeling and performance framework

All estimation in this paper, for NMF, C-NMF, and AE-NMF, is done by
minimizing the average Frobenius distance between the original input and the
reconstructed input:

$\displaystyle
L_{F}(\bm{V},\bm{\hat{V}})=\frac{||\bm{V}-\bm{\hat{V}}||_{F}}{M\cdot
N}=\frac{1}{M\cdot
N}\sqrt{\sum_{i=1}^{M}\sum_{j=1}^{N}|v_{i,j}-\hat{v}_{i,j}|^{2}}.$ (5)

NMF is estimated using the multiplicative and iterative updating steps defined
by [14]:

$\displaystyle\bm{H}\leftarrow\bm{H}\frac{\bm{VW}^{T}}{\bm{HWW}^{T}},\quad\quad\bm{W}\leftarrow\bm{W}\frac{\bm{H}^{T}\bm{V}}{\bm{H}^{T}\bm{H}\bm{W}}.$
(6)

C-NMF is estimated using the multiplicative updating scheme defined by [5].
This updating scheme reduces to

$\displaystyle\begin{split}\bm{W}_{1}\leftarrow\bm{W}_{1}\sqrt{\frac{(\bm{V}^{T}\bm{V})\bm{W}_{2}^{T}}{(\bm{V}^{T}\bm{V})\bm{W}_{1}\bm{W}_{2}\bm{W}_{2}^{T}}},\quad\quad\bm{W}_{2}^{T}\leftarrow\bm{W}_{2}^{T}\sqrt{\frac{(\bm{V}^{T}\bm{V})\bm{W}_{1}}{\bm{W}_{2}^{T}\bm{W}_{1}^{T}(\bm{V}^{T}\bm{V})\bm{W}_{1}}},\end{split}$
(7)

in the case where $\bm{V}$ is constrained to be in the non-negative domain.

AE-NMF is estimated by an Adam optimizer, where the $t$’th iteration of a
parameter, $w^{(t)}$, is updated additively and iteratively using the
following method:

$\displaystyle\begin{split}w^{(t+1)}\leftarrow
w^{(t)}-\eta\frac{\frac{m^{(t)}}{1-\beta_{1}^{t}}}{\sqrt{\frac{v^{(t)}}{1-\beta_{2}^{t}}}+\varepsilon},\quad\text{where}\quad&m^{(t)}\leftarrow\beta_{1}m^{(t-1)}+(1-\beta_{1})\nabla_{w}L_{t}\\\
\text{and}\quad&v^{(t)}\leftarrow\beta_{2}v^{(t-1)}+(1-\beta_{2})(\nabla_{w}L_{t})^{2}.\end{split}$
(8)

Here $L_{t}$ is the loss function at the $t$’th iteration and
$\varepsilon=10^{-8}$ is a small scalar to prevent division by zero. The
factors $\beta_{1}$ and $\beta_{2}$ are exponential decay rates initialized to
the default values $\beta_{1}=0.9$ and $\beta_{2}=0.999$. The moment vectors
of parameter $w$ at time $t$, $m^{(t)}$ and $v^{(t)}$, are initialized as
$m^{(0)}\leftarrow 0$ and $v{(0)}\leftarrow 0$ [[11]]. Training of AE-NMF is
done using the PyTorch module, version 1.13.1 [[19]], in Python version 3.11.3
[[27]].

####  2.5.1 Average Cosine Similarity

The similarity between two sets of basis vectors is evaluated using the
average cosine similarity (ACS). The cosine similarity for two equal length
vectors, $\tilde{\bm{h}},\hat{\bm{h}}\in\mathbb{R}^{M}$, is defined as

$S_{C}(\tilde{\bm{h}},\hat{\bm{h}})=\frac{\tilde{\bm{h}}\cdot\hat{\bm{h}}}{||\tilde{\bm{h}}||||\bm{\hat{h}}||}.$
(9)

The ACS for two matrices of matched basis vectors
$\tilde{\bm{H}},\hat{\bm{H}}\in\mathbb{R}_{+}^{M\times K}$ with $K$ basis
vectors is defined as

$ACS(\tilde{\bm{H}},\bm{\hat{H}})=\frac{1}{K}\sum_{k=1}^{K}S_{c}(\bm{\tilde{h}}_{k},\bm{\hat{h}}_{k}),$
(10)

where $\tilde{\bm{h}}_{k}$ denotes the $k$’th column in $\tilde{\bm{H}}$, and
$\hat{\bm{h}}_{k}$ is the corresponding basis vector in $\hat{\bm{H}}$.

####  2.5.2 Signature Matching

Given two matrices of basis vectors $\bm{H}\in\mathbb{R}^{M\times K}$ and
$\tilde{\bm{H}}\in\mathbb{R}^{M\times\tilde{K}}$ where $K\leq\tilde{K}$, the
task is to find $K$ pairs $(i,j)$, where $i=1,\dots,K$ and
$j\in\\{1,\dots,\tilde{K}\\}$, such that the sum of cosine similarities (and
thus the ACS) over all pairs
$\sum_{(i,j)}S_{c}(\bm{h}_{i},\bm{\tilde{h}}_{j})$ is maximized, all vectors
in $\bm{H}$ are matched exactly once, and all vectors in $\tilde{\bm{H}}$ are
matched at most once. This combinatorial problem is a linear assignment
problem, where the cost matrix is all cosine distances between two basis
vectors
$\left\\{1-S_{c}(\bm{h}_{i},\tilde{\bm{h}}_{j})\right\\}_{i=1,\dots,K,\newline
j=1,\dots,\tilde{K}}$. This problem is solved using the Hungarian algorithm in
both the balanced ($K=\tilde{K}$) and unbalanced ($K<\tilde{K}$) case [[13]].

##  3 Results

In this section, we compare NMF, C-NMF, and AE-NMF in a simple simulated
example (Section 3.1) and compare NMF and AE-NMF on the ovary, prostate, and
uterus genomes from the Genomics England 100,000 Genomes Project [[25, 26]] by
the ability to reconstruct the input data accurately and to generate stable
and sensible mutational signatures (Section 3.2).

All training was performed with a relative tolerance of $10^{-10}$ as
convergence criteria. Estimation with AE-NMF was performed with an Adam
optimizer with a learning rate of $10^{-4}$. Non-negativity in AE-NMF was
enforced taking the absolute value of the encoding and decoding weight
matrices in the forward pass. The choice of method to enforce non-negativity
in AE-NMF is elaborated in Supplementary Material Section 1.4.

###  3.1 Simulated example

Consider an example with two basis vectors (signatures) consisting of 6
features (mutation types):

$\displaystyle\begin{split}\bm{h}_{1}&=(2,2,1,1,0,0)^{T}/6\\\
\bm{h}_{2}&=(0,0,0,1,1,1)^{T}/3.\end{split}$ (11)

From these basis vectors we simulated 30 samples (patients),
$\\{\bm{v}_{1},\dots,\bm{v}_{30}\\}$, where each sample was simulated with one
of three distinct compositions of $\bm{h}_{1}$ and $\bm{h}_{2}$ and Poisson
distributed noise:

$\displaystyle\bm{v}_{j}=(v_{j,1},v_{j,2},v_{j,3},v_{j,4},v_{j,5},v_{j,6})^{T},$
(12)

for $j=1,\dots,30$ and

$\displaystyle
v_{i,j}\overset{\text{i.i.d}}{\sim}\begin{cases}\text{Po}(180\cdot\bm{h}_{i,1}+20\cdot\bm{h}_{i,2})&\text{for
}j=1,\dots,10\\\
\text{Po}(100\cdot\bm{h}_{i,1}+100\cdot\bm{h}_{i,2})&\text{for
}j=11,\dots,20\\\
\text{Po}(20\cdot\bm{h}_{i,1}+180\cdot\bm{h}_{i,2})&\text{for
}j=21,\dots,30,\\\ \end{cases}$ (13)

for $i=1,\dots,6$. For the full data matrix, $\bm{V}=\left[v_{i,j}\right]$,
the two basis vectors and the corresponding weights (exposures) were estimated
using NMF, C-NMF, and AE-NMF, with the Frobenius norm as the loss function. To
encourage convergence towards the same minimum, C-NMF and AE-NMF were
initialized with the same matrices, $\bm{W}_{1}^{(0)}$ and $\bm{W}_{2}^{(0)}$,
respectively to a $30\times 2$ and a $2\times 30$ matrix with values sampled
uniformly in $[0,1)$. This was not possible for NMF since its factor matrices
have different dimensions. For NMF the matrix of basis vectors was initialized
uniformly in $[0,1)$, and the weight matrix was initialized as
$\bm{W}_{2}^{(0)}$ from C-NMF and AE-NMF.

Figure 2: Estimated and true weights and basis vectors of the simulated data.
Note the light green lines representing C-NMF coincide with the the dashed,
coral lines of AE-NMF since the estimates found by both methods are almost
identical.

Figure 2 shows the estimated basis vectors and weights for each method. In
this example, identical factor matrices are recovered by the C-NMF and AE-NMF
updates, illustrated by the coinciding dashed green and coral lines, while NMF
yields different results, illustrated by the blue lines. Furthermore, the NMF
solution reconstructs the input with higher accuracy, with a reconstruction
error of $12.95$ compared to the higher values of $13.64$ for AE-NMF and
$13.00$ for C-NMF.

###  3.2 Cancer Data Analysis

To emulate the terminology used in the cancer data analysis, the matrix of
basis vectors will be denoted as the signatures, and the weight matrix will be
denoted the exposures.

The number of signatures in each diagnosis was determined with the method
described in the Supplementary Material Section 1.3 with ten training/test
splits for each diagnosis. The test errors as a function of
$K\in\\{2,\dots,12\\}$ are depicted in Figure S1. Algorithm S2 in the
Supplementary Material yielded three signatures for AE-NMF and four signatures
for C-NMF and NMF in the ovary cohort; four, four and six signatures for AE-
NMF, C-NMF and NMF, respectively in the prostate cohort and four, six, and 11
signatures for AE-NMF, C-NMF, and NMF respectively in the uterus cohort. Thus,
the number of signatures used in the cancer data analyses is chosen, using the
weighted average in Equation (S4) in the Supplementary Material, as four for
the ovary cohort, five for the prostate cohort and eight for the uterus
cohort.

As shown in Section 2.4 C-NMF and AE-NMF are mathematically equivalent.
Furthermore, since the resulting factor matrices of C-NMF and AE-NMF were
identical in the simulated example and performed similarly in the cancer data
analysis with respect to reconstruction error (Figure S1 and Figure S3) and
consistency (Figure S4), we consider C-NMF and AE-NMF as practically
equivalent. Thus, the following analyses will be performed comparing only NMF
and AE-NMF.

####  3.2.1 Extraction performance

Figure 3: Second principal component plotted agains the first principal
component of the de novo extracted signatures from the 30 train/test splits
for AE-NMF and NMF (columns) and each diagnosis (rows). Points are colored by
the PAM clustering assignment, and the cluster mediod is highlighted with a
black outline. NMF: non-negative matrix factorization; AE: autoencoder; PAM:
partition around mediods.

For each cohort, we divided the patients’ mutational profile data into 30
80/20 train/test set splits. De novo extractions by NMF and AE-NMF were
performed on the training set yielding a set of training errors and refits
were performed on the test sets yielding a set of test errors. Plots of the
first and second principal component of all extracted signatures from the 30
training sets are depicted in Figure 3.

The procedure of calculating the test errors is detailed in the Supplementary
Material Section 1.2. Boxplots of the training and test errors when
reconstructing the input matrix for each method and diagnosis are shown in
Figure 4. Table 1 reports the average training and test error across the 30
splits for each method and diagnosis.

Figure 4: Boxplots of the training and test errors of 30 train/test splits of
the ovary, prostate and uterus cohorts. NMF and AE-NMF errors resulting from
the same splits are connected by a black line. The boxes are colored
corresponding to the method used, and a green diamond depicts the average
error. The y-axis is on log 10 scale.

Considering the reconstruction errors in Figure 4 and Table 1, it is apparent
that NMF consistently performs better than AE-NMF in terms of reconstructing
the input data. This is the case both on average and in the majority of splits
as shown by Figure 4. The ratios in Table 1 reveal that the difference is more
expressed in the training splits than the test splits.

|  |  | Average Error |
---|---|---|---|---
Cohort | Split | $\bm{n}$ | NMF | AE-NMF | Ratio
Ovary | Train | $418$ | $\bm{5.92\cdot 10^{3}}$ | $1.29\cdot 10^{4}$ | $2.17$
Test | $105$ | $\bm{1.51\cdot 10^{6}}$ | $1.72\cdot 10^{6}$ | $1.63$
Prostate | Train | $248$ | $\bm{1.06\cdot 10^{2}}$ | $2.39\cdot 10^{2}$ | $2.24$
Test | $63$ | $\bm{2.81\cdot 10^{3}}$ | $3.16\cdot 10^{3}$ | $1.23$
Uterus | Train | $570$ | $\bm{5.33\cdot 10^{5}}$ | $1.39\cdot 10^{6}$ | $2.61$
Test | $143$ | $\bm{1.11\cdot 10^{6}}$ | $2.57\cdot 10^{6}$ | $2.45$
Table 1: Average training and test error between the input and reconstructed
data for each method across the 30 splits of the ovary, prostate, and uterus
cohort and the average ratio between the errors. The lowest reconstruction
error in each cohort is highlighted in bold.

For each method and diagnosis the analyses yielded 30 signature sets, one from
each training set. The consistency of the estimated signatures extracted
within each method is investigated by calculating the ACS between each
pairwise combination of signatures extracted using a given method across the
30 splits of the data matrix. This yields a total of $\binom{30}{2}=435$
comparisons for each cohort and method, and resulted in an average ACS
consistency of $0.91$ for NMF and $0.85$ for AE-NMF in the ovary cohort,
$0.86$ for NMF and $0.88$ for AE-NMF in the prostate cohort, and $0.94$ for
NMF and $0.91$ for AE-NMF in the uterus cohort. To asses whether the
difference in mean consistency between NMF and AE-NMF were significant, two-
sided t-tests for equal means were performed for each cohort. These resulted
in a $p$-value of $3.7\cdot 10^{-18}$ for the ovary cohort, $7.0\cdot
10^{-10}$ for the prostate cohort, and $1.1^{\cdot}10^{-11}$ for the uterus
cohort, thus the mean consistency for NMF and AE-NMF can not be assumed to be
equal for any cohort. The consistencies are depicted by boxplots in Figure 5,
with green diamonds marking the average consistencies listed above. This
reveals generally comparable consistencies of NMF and AE-NMF signatures.

Figure 5: Boxplots of the ACS between each combination of signatures,
extracted using NMF or AE-NMF and across the 30 splits of the data matrix for
the ovary, prostate, and uterus cohort. The average is marked by a green
diamond.

####  3.2.2 COSMIC Validation

To compare the estimated signatures to those of the leading library of
mutational signatures, COSMIC v. 3.4 [[24]], the signatures were clustered to
form sets of consensus signatures using the partitioning around medoids (PAM)
algorithm [[9]] with the number of clusters equal to the number of signatures
used in the initial extractions. The clustering is depicted in Figure 3 where
the first and second principal components of all de novo signatures are
depicted and points are colored by their assigned PAM clustering. The
consensus signatures were subsequently matched to the COSMIC signatures. The
matched COSMIC signatures along with their cosine similarity can be seen in
Table 2.

Cohort | Signature | NMF | AE-NMF
---|---|---|---
Ovary | S1 | SBS10a | 0.93 | SBS10a | 0.93
S2 | SBS3 | 0.95 | SBS40a | 0.80
S3 | SBS10c | 0.70 | SBS10c | 0.74
S4 | SBS44 | 0.84 | SBS44 | 0.82
ACS |  | 0.80 |  | 0.82
Prostate | S1 | SBS17b | 0.84 | SBS17b | 0.77
S2 | SBS33 | 0.95 | SBS33 | 0.93
S3 | SBS1 | 0.94 | SBS6 | 0.90
S4 | SBS44 | 0.86 | SBS44 | 0.85
S5 | SBS40a | 0.90 | SBS40a | 0.87
ACS |  | 0.90 |  | 0.86
Uterus | S1 | SBS10d | 0.96 | SBS10d | 0.95
S2 | SBS10c | 0.84 | SBS10a | 0.95
S3 | SBS10a | 0.97 | SBS10c | 0.75
S4 | SBS6 | 0.88 | SBS6 | 0.84
S5 | SBS20 | 0.93 | SBS20 | 0.90
S6 | SBS10b | 0.78 | SBS2 | 0.71
S7 | SBS28 | 0.96 | SBS28 | 0.70
S8 | SBS26 | 0.83 | SBS26 | 0.82
ACS |  | 0.89 |  | 0.83
Table 2: The COSMIC signatures matched to the consensus signatures and the
corresponding cosine similarity. The average cosine similarity (ACS) between a
consensus signature set and the matched COSMIC signatures is reported in the
last line for each diagnosis. Signatures that have previously been observed in
patients with the same diagnosis are highlighted with bold text.

The ACS between the consensus signatures and their matched COSMIC signatures
is 0.80 for NMF and 0.82 for AE-NMF in the ovary cohort, 0.90 for NMF and 0.86
for AE-NMF in the prostate cohort and 0.90 for NMF and 0.83 for AE-NMF in the
uterus cohort. This reveals a slightly higher average similarity to COSMIC for
NMF than AE-NMF. Additionally, the majority of NMF signatures also have better
COSMIC matches than the corresponding AE-NMF signature. The matched COSMIC
signatures for all splits before clustering for NMF and AE-NMF can be seen in
Supplementary Figure S5-S7. Of all matched consensus signatures, the following
proportion has been previously observed in the corresponding diagnosis for
each cohort: 1/4 of NMF and 0/4 of AE-NMF signatures in the ovary cohort, 2/5
of both NMF and AE-NMF signatures in the prostate cohort, and 6/8 of both NMF
and AE-NMF signatures in the uterus cohort.

Overall both NMF and AE-NMF show a high degree of conformity with the COSMIC
SBS signatures, and the extracted signatures are relevant to the diagnoses in
which they have been identified. In signature cosistency, conformity with
COSMIC and choosing relevant signatures, NMF and AE-NMF perform similarly,
perhaps with a slight advantage to NMF.

##  4 Discussion

In this study, we compare NMF and AE-NMF by their ability to extract valid and
consistent basis vectors and creating accurate reconstructions of the input
data. We assert that such comparisons are theoretically meaningful since we
demonstrate that AE-NMF and C-NMF are mathematically equivalent.

The study focuses on extracting mutational signatures in the ovary, prostate,
and uterus cancer genomes of Genomics England’s 100,000 Genomes cohort. NMF
consistently outperformed AE-NMF in terms of reconstruction error; the
differences being more expressed in the training splits than in the test
splits. While AE-NMF constrains parameters to convex combinations of patients’
profiles, NMF can freely assume non-negative values, giving it an advantage in
training set reconstruction. When reconstructing the test set the signature
matrix is fixed, and the task is thus identical for NMF and AE-NMF. One could
expect the constrained nature of AE-NMF to regularize the signatures such that
they reconstruct the test splits better compared to NMF, but as this was not
the case in this study, it suggests that AE-NMF signatures may be generally
less informative than the corresponding NMF signatures. The COSMIC validation
revealed that both models recovered relevant signatures with high cosine
similarity, with a slight advantage to NMF. This advantage is likely driven by
the fact that the majority of COSMIC signatures are extracted using NMF-based
methods [[2]].

The mathematical equivalence between convex NMF and non-negative autoencoders
enables the interpretation of parameters from AE-NMF to be similar to that of
NMF. Thus, the mathematical equivalence between convex NMF and non-negative
autoencoders identified in this study is a necessary link in properly
comparing AE-NMF and NMF. A link that has been missing in previous attempts to
replace NMF with autoencoders while using the same interpretation. [23] also
came to the conclusion that an autoencoder with the architecture of AE-NMF
yields a hidden layer consisting of convex combinations of the data points but
did not make the connection to convex NMF.

In practice we observed that C-NMF and AE-NMF will yield identical solutions
in a sufficiently simple setup. When increasing the complexity of the problem
both methods perform similarly in terms of reconstruction error and stability
of the basis vectors, albeit finding different solutions within this minimum.

The architecture of AE-NMF is atypical for autoencoders by its shallow and
linear nature and by transposing the input data matrix. By orienting the input
matrix, $\bm{V}$, with features ($M$) as rows and observations ($N$) as
columns the architecture will stray from how data is conventionally passed
through neural networks, but this orientation is not uncommon in the
literature of non-negative autoencoders [[10, 20, 23]]. One direct consequence
of this orientation is that it will not be possible to fit a separate sample
by passing it through the trained network, a compelling attribute of neural
networks, as the trained parameters will be observation-specific. Furthermore,
the shallowness of AE-NMF favours the capture of linear relationships over
more complex patterns. It is these strict architectural choices that yields
the equivalence to convex NMF, ensuring a precise understanding of the
interpretation of the parameters. Other efforts, such as the MUSE-XAE
autoencoder for mutational signature extraction [[18]], have utilized the
benefits of classic autoencoders with deeper extraction and conventional
orientation of the input data matrix but at the cost of jeopardizing the link
to NMF and thus the exact interpretation of the parameters, since this
autoencoder is not equivalent to C-NMF.

All methods were optimized by minimizing the Frobenius loss function since
this is the only loss function for conventional convex NMF by [5] that
presently has updating steps. This loss is only asymptotically efficient given
Gaussian distributed input data, which is a questionable assumption since the
data considered in this study consists of count data. Choosing a loss function
adapted to Poisson may therefore be more adequate for this problem [[2, 4]].
An even more sophisticated error model would be the Negative Binomial
distribution [[21]]. The Negative Binomial distribution can model
overdispersion in the mutational counts, which is often present on the
patient-specific level. In future work, it would be interesting to see how
training with a more appropriate loss function affects how NMF and AE-NMF
compare in reconstruction accuracy and stability and whether the results from
this study generalize to the Kullback-Leibler divergence.

In this study a relative convergence criteria of $10^{-10}$ was used for all
analyses. Using such a low convergence criteria promoted that the three
methods, based on two vastly different updating schemes, converged towards
similar minima and, thus, established a common ground for comparison. A lower
relative tolerance in the real data analyses was computationally infeasible.
In particular, the high tolerance made especially C-NMF extremely time
consuming in cases with many patients and/or signatures. The bootstrap method
for determining the number of signatures were limited to ten splits for the
same reason. If the aim is to just extract signatures using either method, we
do not necessarily recommend training with such a low tolerance.

##  5 Conclusion

This study compares NMF with non-negative autoencoders, facilitated by the
mathematical equivalence between non-negative autoencoders and convex NMF.
This bridges a crucial gap in the comparison of NMF and non-negative
autoencoders by offering insights into parameter interpretation that were
previously lacking. The choice between convex NMF and its autoencoder
equivalent is a question of choosing between a multiplicative or gradient
descent-based optimizing algorithm to solve the same optimization problem.
Thus, the non-negative autoencoder described in this study can be used as a
faster alternative to convex NMF. Non-negative autoencoders exhibit higher
reconstruction errors and similar consistencies compared to NMF, therefore the
non-negative autoencoder investigated in this study is not a suitable
alternative to NMF in mutational signature extraction. Thus, this study
underscores the significance of methodological considerations when replacing
NMF with non-negative autoencoders. On the other hand autoencoders hold
promise for modeling non-linearity, a capability absent in NMF, but such
advancements are made at the cost of exact parameter interpretation.

## Acknowledgements

This work was supported by the Danish Data Science Academy, which is funded by
the Novo Nordisk Foundation (NNF21SA0069429) and VILLUM FONDEN (40516), the
Novo Nordisk Foundation (NNF21OC0069105), and the Research Hive ”REPAIR”,
which is funded by Aalborg University Hospital.

## Code- and Data avaliability

Code used for this study can be found at Github and the Genomics England WGS
data used in this study was provided by [26] on Zendo.

## Supplementary Materials

### Supplementary Methods

The supplementary methods contains a section on how test errors are calculated
for all methods, the algorithms for choosing the optimal number of basis
vectors, and the theoretical overview of methods for enforcing non-negativity
in the autoencoder.

### Supplementary Results

The supplementary results contains the results used to choose the optimal
number of basis vectors, the results used to determine the method to constrain
non-negativity, and additional figures and tables.

## References

  * [1] Ludmil Alexandrov et al.  “Deciphering Signatures of Mutational Processes Operative in Human Cancer”  In _Cell Reports_ 3, 2013, pp. 246–259  DOI: 10.1016/j.celrep.2012.12.008
  * [2] Ludmil Alexandrov et al.  “The repertoire of mutational signatures in human cancer”  In _Nature_ 578, 2020, pp. 94–101  DOI: 10.1038/s41586-020-1943-3
  * [3] Francis Blokzijl, Roel Janssen, Ruben Boxtel and Edwin Cuppen  “MutationalPatterns: Comprehensive genome-wide analysis of mutational processes”  In _Genome Medicine_ 10, 2018  DOI: 10.1186/s13073-018-0539-0
  * [4] Andrea Degasperi et al.  “Substitution mutational signatures in whole-genome–sequenced cancers in the UK population”  In _Science_ 376, 2022  DOI: 10.1126/science.abl9283
  * [5] Chris H.Q. Ding, Tao Li and Michael I. Jordan  “Convex and Semi-Nonnegative Matrix Factorizations”  In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 32, 2010, pp. 45–55  DOI: 10.1109/TPAMI.2008.277
  * [6] Hao Fang, Aihua Li, Huoxi Xu and Tao Wang  “Sparsity-Constrained Deep Nonnegative Matrix Factorization for Hyperspectral Unmixing”  In _IEEE Geoscience and Remote Sensing Letters_ 15, 2018, pp. 1105–1109  DOI: 10.1109/LGRS.2018.2823425
  * [7] Ehsan Hosseini-Asl, Jacek Zurada and Olfa Nasraoui  “Deep Learning of Part-Based Representation of Data Using Sparse Autoencoders With Nonnegativity Constraints”  In _IEEE Transactions on Neural Networks and Learning Systems_ 27, 2016, pp. 2486–2498  DOI: 10.1109/TNNLS.2015.2479223
  * [8] S.M. Ashiqul Islam et al.  “Uncovering novel mutational signatures by de novo extraction with SigProfilerExtractor”  In _Cell Genomics_ 2, 2022, pp. 100179  DOI: https://doi.org/10.1016/j.xgen.2022.100179
  * [9] Leonard Kaufman and Peter J. Rousseeuw  “Partitioning Around Medoids (Program PAM)”  In _Finding Groups in Data_ John Wiley & Sons, Ltd, 1990, pp. 68–125  DOI: https://doi.org/10.1002/9780470316801.ch2
  * [10] Alaa El Khatib, Shimeng Huang, Ali Ghodsi and Fakhri Karray  “Nonnegative Matrix Factorization Using Autoencoders And Exponentiated Gradient Descent”  In _Int. Joint Conf. on Neural Networks (IJCNN)_ , 2018, pp. 1–8  DOI: 10.1109/IJCNN.2018.8489242
  * [11] Diederik P. Kingma and Jimmy Ba  “Adam: A Method for Stochastic Optimization”  In _arXiv preprint arXiv:1412.6980_ , 2017
  * [12] Mark A. Kramer  “Nonlinear principal component analysis using autoassociative neural networks”  In _AIChE Journal_ 37, 1991, pp. 233–243  DOI: 10.1002/aic.690370209
  * [13] Harold Kuhn  “The Hungarian Method for The Assignment Problem”  In _Naval Research Logistics Quarterly_ 2, 1955, pp. 83 –97  DOI: 10.1002/nav.3800020109
  * [14] Daniel Lee and H. Sebastian Seung  “Algorithms for Non-negative Matrix Factorization”  In _Advances in Neural Information Processing Systems_ MIT Press, 2000, pp. 556–562
  * [15] Daniel Lee and Hyunjune Seung  “Learning the Parts of Objects by Non-Negative Matrix Factorization”  In _Nature_ 401, 1999, pp. 788–791  DOI: 10.1038/44565
  * [16] Andre Lemme, René Felix Reinhart and Jochen Jakob Steil  “Online learning and generalization of parts-based image representations by non-negative sparse autoencoders”  In _Neural Networks_ 33, 2012, pp. 194–203  DOI: 10.1016/j.neunet.2012.05.003
  * [17] Serena Nik-Zainal et al.  “The genome as a record of environmental exposure”  In _Mutagenesis_ 30, 2015  DOI: 10.1093/mutage/gev073
  * [18] Corrado Pancotti et al.  “MUSE-XAE: MUtational Signature Extraction with eXplainable AutoEncoder enhances tumour type classification”  In _bioRxiv preprint 2023.10.23.562664_ , 2023  DOI: 10.1101/2023.10.23.562664
  * [19] Adam Paszke et al.  “PyTorch: An Imperative Style, High-Performance Deep Learning Library”  In _Advances in Neural Information Processing Systems_ , 2019
  * [20] Guangsheng Pei et al.  “Decoding whole-genome mutational signatures in 37 human pan-cancers by denoising sparse autoencoder neural network”  In _Oncogene_ 39, 2020, pp. 5031–5041  DOI: 10.1038/s41388-020-1343-z
  * [21] Marta Pelizzola, Ragnhild Laursen and Asger Hobolth  “Model selection and robust inference of mutational signatures using Negative Binomial non-negative matrix factorization”  In _BMC Bioinformatics_ 24, 2023  DOI: 10.1186/s12859-023-05304-1
  * [22] Paris Smaragdis and Shrikant Venkataramani  “A neural network alternative to non-negative audio models”  In _IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP)_ , 2017, pp. 86–90  DOI: 10.1109/ICASSP.2017.7952123
  * [23] Steven Squires, Adam Prügel Bennett and Mahesan Niranjan  “A Variational Autoencoder for Probabilistic Non-Negative Matrix Factorisation”  In _arXiv preprint arXiv:1906.05912_ , 2019
  * [24] John G Tate et al.  “COSMIC: the Catalogue Of Somatic Mutations In Cancer”  In _Nucleic Acids Research_ 47.D1, 2018, pp. D941–D947  DOI: 10.1093/nar/gky1015
  * [25] Clare Turnbull  “Introducing Whole Genome Sequencing into routine cancer care: The Genomics England 100,000 Genomes project”  In _Annals of Oncology_ 29, 2018, pp. 784–787  DOI: 10.1093/annonc/mdy054
  * [26] Ernest Turro et al.  “Whole-genome sequencing of patients with rare diseases in a national health system”  In _Nature_ 583, 2020, pp. 96–102  DOI: 10.1038/s41586-020-2434-2
  * [27] Guido Van Rossum and Fred L. Drake  “Python 3 Reference Manual”  Scotts Valley, CA: CreateSpace, 2009
  * [28] Yigitcan Özer, Jonathan Hansen, Tim Zunner and Meinard Müller  “Investigating Nonnegative Autoencoders for Efficient Audio Decomposition”  In _30th European Signal Processing Conf. (EUSIPCO)_ , 2022, pp. 254–258  DOI: 10.23919/EUSIPCO55093.2022.9909787