# Laser induced enhanced coupling between photons and squeezed magnons in
antiferromagnets

J. C. G. Henriques1, T. V. C. Antão1 and N. M. R. Peres1,2 1Department and
Centre of Physics, University of Minho, Campus of Gualtar, 4710-057, Braga,
Portugal  2International Iberian Nanotechnology Laboratory (INL), Av. Mestre
Jose Veiga, 4715-330, Braga, Portugal

###### Abstract

In this paper we consider a honeycomb antiferromagnet subject to an external
laser field. Obtaining a time-independent effective Hamiltonian, we find that
the external laser renormalizes the exchange interaction between the in-plane
components of the spin-operators, and induces a synthetic Dzyaloshinskii-Moria
interaction (DMI) between second neighbors. The former allows the control of
the magnon dispersion’s bandwidth and the latter breaks time-reversal symmetry
inducing non-reciprocity in momentum space. The eigen-excitations of the
system correspond to squeezed magnons whose squeezing parameters depend on the
properties of the laser. When studying how these spin excitations couple with
cavity photons, we obtain a coupling strength which can be enhanced by an
order of magnitude via careful tuning of the laser’s intensity, when compared
to the case where the laser is absent. The transmission plots through the
cavity are presented, allowing the mapping of the magnons’ dispersion
relation.

##  I Introduction

The field of magnon spintronics deals with the study and manipulation of
magnetic excitations, also known as magnons, in ordered magnets [1, 2]. These
spin excitations attract considerable interest from the scientific community
due to the possibility of being used as information carriers, since they
present nanometre wavelengths and reduced losses due to Joule heating when
compared with traditional electronic devices [1, 3, 4, 5, 6]. Due to the
possibility of combining magnons with cavity photons [7, 8, 9] and
superconducting qubits [10, 11], the study of spin excitations also found its
way to the field of quantum information [12]. The combined study of magnon
spintronics and its applications in quantum information science gave rise to
the field of quantum magnonics.

Due to their bosonic nature, magnons share various features with other bosons,
such as phonons and photons. In fact, when studying the quantum properties of
magnons, similar ideas and methods to the ones usually found in quantum optics
[13, 14] appear; one such idea is that of squeezed states. This type of
quantum state has been thoroughly studied with photons, where these states are
usually obtained out of equilibrium through four wave mixing or parametric
processes [13, 14, 15]. The signature feature of squeezed states is the
possibility of reducing the uncertainty associated with a given observable by
increasing the uncertainty of another one, in such a way that, when combined,
the two uncertainties still respect Heisenberg’s uncertainty principle. This
type of quantum state of light has bee used, for example, in the detection of
gravitational waves [16].

Contrarily to what is found with photons, where the study and generation of
squeezed states is a mature field [17], the study of magnon squeezing has only
began to gain traction recently. This topic, however, is a rather interesting
one due to the large values of squeezing found in these systems [18], the
increased spin carried by a squeezed magnon [19], enhanced magnon-magnon
coupling [20], magnon entanglement [21], among others [22]. Moreover,
contrarily to photons, the exploration of such phenomena in magnetically
ordered systems allows for its implementation in on-chip nanodevices [23].

Regarding the generation of squeezed magnons, two distinct approaches can be
employed. On the one hand, squeezing can be obtained by driving the magnetic
system out of equilibrium while coupled with an optomechanical cavity [24, 25,
26]. On the other hand, and in stark contrast with photons, anisotropic
ferromagnets and isotropic antiferromagnets host equilibrium magnon squeezing,
that is, the eigen-excitations of such systems naturally present squeezing,
making them robust against environment perturbations [2, 18, 23].

Recent works focus on the manipulation of magnetic materials via the
application of high-frequency laser fields. In Ref. [27], the authors uncover
an ultrafast Floquet magnonic topological phase transition in a laser-driven
skyrmion crystal, and demonstrate how single skyrmions can be set in motion
with a velocity and propagation direction that can be tuned by the laser. In
Ref. [28] a study of Floquet topological magnons in ferromagnets was
performed, where the author explored how the application of an external laser
field may be used to generate a synthetic Dzyaloshinskii-Moria interaction
(DMI) [29] through the appearance of a time-dependent Aharonov-Casher phase
[30]. This laser induced effect lead to the transformation of Dirac magnons
into magnon Chern insulators. Inspired by this work, and motivated by the
growing interest of the scientific community on the magnonic response of
antiferromagnets, in this paper we study the effect of applying an external
laser field to a honeycomb antiferromagnet, and discuss how it can be used to
tune the properties of the squeezed magnons hosted by the system.

The text is organized as follows: In Sec. II we start by introducing the
Hamiltonian of a honeycomb antiferromagnt, and how it is modified by the
presence of the laser. Afterwards the Hamiltonian is diagonalized with a
Bogoliubov transformation. In Sec. III we consider the system to be placed in
an optical cavity, and study how the laser can be used to enhance the magnon-
photon coupling; plots of the transmission through the cavity are also given.
In Sec. IV we discuss the feasibility of experimentally realizing this type of
system, and in Sec. V we give our final remarks.

Figure 1: (a) Schematic representation of the antiferromagnet honeycomb
lattice. The distance between nearest neighbors is $a$. (b) Representation of
the three $\boldsymbol{\delta}$ vectors connecting one site to its three
nearest neighbors. Also depicted are the 6 vectors $\boldsymbol{b}_{i}$ which
enter in the definition of the Dzyaloshinskii-Moria interaction. (c) Schematic
representation of the first Brillouin zone (1BZ), with the Dirac points
$a\boldsymbol{K}_{\pm}=\pm 2\pi\left(2/3\sqrt{3},0\right)$ highlighted.

##  II Model Hamiltonian

In this section, we introduce the Hamiltonian of a honeycomb antiferromagnet.
We focus on this type of lattice since it is a common one in two-dimensional
magnetic materials, and as such is relevant in this field of research.
Afterwards, the effect of applying a circularly polarized laser field to the
system will be explored, with the introduction of a Floquet effective
Hamiltonian. This process is analogous to the one presented for a ferromagnet
in [28]. Finally, the resulting Hamiltonian will be diagonalized, and its
properties discussed.

###  II.1 Antiferromagnetic Synthetic Hamiltonian

Consider the antiferromagnetic honeycomb lattice depicted in Fig. 1.
Accounting only for nearest neighbor interactions, and assuming an identical
easy axis anisotropy for both sublattices ($A$ and $B$), we write

$\displaystyle H$ $\displaystyle=J\sum_{i\in
A}\sum_{j=1}^{3}\boldsymbol{S}_{A}(\boldsymbol{r}_{i})\cdot\boldsymbol{S}_{B}(\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j})$
$\displaystyle+K\sum_{i\in
A}S_{A}^{z}(\boldsymbol{r}_{i})S_{A}^{z}(\boldsymbol{r}_{i})+K\sum_{i\in
B}S_{B}^{z}(\boldsymbol{r}_{i})S_{B}^{z}(\boldsymbol{r}_{i})$ (1)

where $J>0$ is the antiferromagnetic exchange coupling, $K$ quantifies the
easy axis anisotropy,
$\boldsymbol{S}_{A/B}=(S_{A/B}^{x},S_{A/B}^{y},S_{A/B}^{z})$ are the spin
operators of the sublattice $A/B$ and
$\boldsymbol{\delta}_{1}=a\left(0,1\right)$,
$\boldsymbol{\delta}_{2}=a\left(1,-\sqrt{3}\right)/2$,
$\boldsymbol{\delta}_{3}=-a\left(1,\sqrt{3}\right)/2$ are the vectors
connecting each site of sublattice $A$ site to its three nearest neighbors,
belonging to sublattice $B$, with $a$ the nearest neighbor distance. This
Hamiltonian is commonly employed to study a variety of antiferromagnets, and
can be extended to include effects we neglect here, such as an in-plane
anisotropy.

Inspired by [28], we now consider that a circularly polarized laser field is
applied to this system. The electric field component associated with the laser
is given by
$\boldsymbol{E}_{0}=E_{0}\left(\tau\cos\omega_{l}t,\sin\omega_{l}t,0\right)$,
where $\tau=\pm 1$ encodes the two possible circular polarizations and
$\omega_{l}$ stands for the laser’s frequency. Now, following [28], we
consider that the laser couples with the system through the Aharonov-Casher
effect [30]. This effect, which can be viewed as complementary to the
Aharonov-Bohm effect, presents itself by giving rise to a complex phase due to
the interaction of a magnetic momentum with an electric field. For the system
at hand, this translates to the appearance of a time dependent complex
exponential, termed the Aharanov-Casher phase. This time dependent phase is a
generalization of what is found in the static limit [31], and has been
discussed before in several works [32, 33, 34, 35, 36, 37, 38]. The
Hamiltonian of the system in such a situation becomes time dependent, and may
be expressed as

$\displaystyle H(t)$ $\displaystyle=J\sum_{i\in
A}\sum_{j=1}^{3}S_{A}^{z}(\boldsymbol{r}_{i})S_{B}^{z}(\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j})$
$\displaystyle+\frac{J}{2}\sum_{i\in
A}\sum_{j=1}^{3}\left[S_{A}^{+}(\boldsymbol{r}_{i})S_{B}^{-}(\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j})e^{-i\theta_{\boldsymbol{r}_{i},\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j}}}+\textrm{h.c.}\right]$
$\displaystyle+K\sum_{i\in
A}S_{A}^{z}(\boldsymbol{r}_{i})S_{A}^{z}(\boldsymbol{r}_{i})+K\sum_{i\in
B}S_{B}^{z}(\boldsymbol{r}_{i})S_{B}^{z}(\boldsymbol{r}_{i})$ (2)

where we introduced $S_{\alpha}^{\pm}=S_{\alpha}^{x}\pm S_{\alpha}^{y}$ with
$\alpha=A,B$. The time dependent Aharonov-Casher phase is defined as

$\theta_{\boldsymbol{r}_{i},\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j}}=\lambda\sin\left(\omega_{l}t-\tau\phi_{\boldsymbol{r}_{i},\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j}}\right),$
(3)

with $\lambda=\tau g\mu_{B}E_{0}a/\hbar c^{2}$ where $\mu_{B}$ is the Bohr
magneton and $g$ the Landé g-factor;
$\phi_{\boldsymbol{r}_{i},\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j}}$
corresponds to the angle between the vectors $\boldsymbol{r}_{i}$ and
$\boldsymbol{r}_{i}+\boldsymbol{\delta}_{j}$. Note that since $\lambda\propto
E_{0}$ the value of $\lambda$ is proportional to the laser’s intensity. To
avoid working with a time dependent Hamiltonian we shall make use of Floquet
theory [39, 40, 41, 42]. This is a perturbative approach allows us to obtain
an effective time independent Hamiltonian from $H(t)$. This is achieved by
expressing the effective Hamiltonian in a high frequency series expansion,
which up to first order reads

$H^{\textrm{eff}}=H^{(0)}+\sum_{m=1}^{\infty}\frac{\left[H^{(m)},H^{(-m)}\right]}{m\hbar\omega_{l}},$
(4)

where $H^{(m)}$ is the $m$-th Fourier component of $H(t)$. In order for this
effective Hamiltonian to be valid, the energy of the laser should be larger
than the energy scale of the initial system (set by $J$). For lower laser
energies more terms would have to be accounted for when defining
$H^{\textrm{eff}}$. The leading term $H^{(0)}$ corresponds, perhaps quite
unsurprisingly, to an average of the Hamiltonian time dependent Hamiltonian
over a period of the driving laser field. It is the most relevant term for our
purposes, with the first order correction introducing only small effects.

We shall work within the linear spin-wave theory. Thus, we introduce the
linearized Holstein-Primakoff transformations, which for the antiferromagnetic
case read [43]

$\displaystyle S_{A}^{+}(\boldsymbol{r}_{i})$
$\displaystyle=\sqrt{2S}a_{\boldsymbol{r}_{i}},\quad
S_{A}^{z}(\boldsymbol{r}_{i})=S-a_{\boldsymbol{r}_{i}}^{\dagger}a_{\boldsymbol{r}_{i}}$
(5) $\displaystyle S_{B}^{+}(\boldsymbol{r}_{i})$
$\displaystyle=\sqrt{2S}b_{\boldsymbol{r}_{i}}^{\dagger},\quad
S_{B}^{z}(\boldsymbol{r}_{i})=-S+b_{\boldsymbol{r}_{i}}^{\dagger}b_{\boldsymbol{r}_{i}},$
(6)

where $a_{\boldsymbol{r}_{i}}/a_{\boldsymbol{r}_{i}}^{\dagger}$ and
$b_{\boldsymbol{r}_{i}}/b_{\boldsymbol{r}_{i}}^{\dagger}$ are bosonic
operators, which we refer to as the annihilation/creation operators of
sublattice magnons, and we note
$\left(S_{A/B}^{\pm}\right)^{\dagger}=S_{A/B}^{\mp}$. Introducing the Fourier
representation of the Holstein-Primakoff annihilation/creation operators,
$a_{\boldsymbol{r}_{i}}=\frac{1}{\sqrt{N}}\sum_{\boldsymbol{k}}e^{i\boldsymbol{k}\cdot\bm{r}_{i}}a_{\boldsymbol{k}}$,
we obtain the following time-independent effective Hamiltonian in momentum
space

$\displaystyle H^{\textrm{eff}}=$ $\displaystyle
JS\sum_{\boldsymbol{k}}\left(3+\frac{2K}{J}\right)\left(a_{\boldsymbol{k}}^{\dagger}a_{\boldsymbol{k}}+b_{\boldsymbol{k}}^{\dagger}b_{\boldsymbol{k}}\right)$
$\displaystyle+$ $\displaystyle
JS\mathcal{J}_{0}\left(\lambda\right)\sum_{\boldsymbol{k}}\left(a_{\boldsymbol{k}}b_{\boldsymbol{-k}}\phi_{\boldsymbol{k}}+a_{\boldsymbol{k}}^{\dagger}b_{-\boldsymbol{k}}^{\dagger}\phi_{\boldsymbol{k}}^{*}\right)$
$\displaystyle-$ $\displaystyle
JS\sum_{\boldsymbol{k}}D_{\boldsymbol{k}}\left(\omega_{l},\lambda\right)b_{\boldsymbol{k}}^{\dagger}b_{\boldsymbol{k}}+JS\sum_{\boldsymbol{k}}D_{\boldsymbol{k}}\left(\omega_{l},\lambda\right)a_{\boldsymbol{k}}^{\dagger}a_{\boldsymbol{k}}$
(7)

with
$\phi_{\boldsymbol{k}}=\sum_{j=1}^{3}e^{-i\boldsymbol{k}\cdot\boldsymbol{\delta}_{j}}$
a geometric factor stemming from the nearest neighbor interactions. The first
order correction arises in the previous equation in the form of the
coefficients

$D_{\boldsymbol{k}}(\omega_{l},\lambda)=\frac{4JS}{\hbar\omega_{l}}\mathcal{V}_{\textrm{DM}}\left(\lambda\right)\textrm{Im}\sum_{j=1}^{6}(-1)^{j}e^{-i\boldsymbol{k}\cdot\boldsymbol{b}_{j}},$
(8)

where the $\boldsymbol{b}_{j}$ are shown in Fig. 1 and correspond to the
vectors connecting a given site to its six next nearest neighbors;
$\mathcal{V}_{\textrm{DM}}\left(\lambda\right)=\sum_{m\geq
1}^{\infty}\left[\mathcal{J}_{m}^{2}(\lambda)/2m\right]\sin\left(\tau
m2\pi/3\right)$ and $\mathcal{J}_{m}(x)$ is the cylindrical Bessel function of
the first kind of order $m$. The terms proportional to $D_{\boldsymbol{k}}$
correspond to a synthetic Dzyaloshinskii-Moria interaction (DMI), which
appears from laser induced couplings between second neighbors; its magnitude
can be tuned simultaneously through the laser’s intensity and frequency. Thus,
comparing this Hamiltonian with the one usually employed to describe an
antiferromagnet [18], we realize that: (i) the exchange coupling between the
$z$-component of the spins and the easy axis anisotropy are unaffected by the
laser and originate the first line of Eq. (7); (ii) the laser is responsible
for the renormalization of the in-plane exchange coupling between nearest
neighbors, reducing its strength, i.e. $J\rightarrow
J\mathcal{J}_{0}(\lambda)$ (this renormalization would also appear on other
lattice configurations), and finally (iii) a synthetic DMI arises from the
laser-induced interaction between next nearest neighbors. We note in passing
that although the DMI is finite for the honeycomb lattice, that is not always
the case, since, for example, it is absent for square lattices. How these
modifications manifest themselves on the properties of the system is studied
below.

###  II.2 Diagonalization

To diagonalize the effective Hamiltonian $H^{\textrm{eff}}$ a Bogoliubov
transformation shall be used, as often is the case when working with
antiferromagnets. Hence, to achieve this, we introduce two new sets of
operators, $\alpha_{\boldsymbol{k}}/\alpha_{\boldsymbol{k}}^{\dagger}$ and
$\beta_{\boldsymbol{k}}/\beta_{\boldsymbol{k}}^{\dagger}$, which we simply
refer to as magnon annihilation/creation operators. These are defined as:

$\displaystyle\alpha_{\boldsymbol{k}}$
$\displaystyle=u_{\boldsymbol{k}}a_{\boldsymbol{k}}+v_{\boldsymbol{k}}b_{-\boldsymbol{k}}^{\dagger}$
(9) $\displaystyle\beta_{\boldsymbol{k}}$
$\displaystyle=v_{\boldsymbol{k}}a_{\boldsymbol{k}}^{\dagger}+u_{\boldsymbol{k}}b_{-\boldsymbol{k}},$
(10)

with $u_{\boldsymbol{k}}=\cosh\frac{\xi_{\boldsymbol{k}}}{2}$ and
$v_{\boldsymbol{k}}=e^{i\theta_{\boldsymbol{k}}}\sinh\frac{\xi_{\boldsymbol{k}}}{2}$
where $\xi_{\boldsymbol{k}}>0$ and $\theta_{\boldsymbol{k}}\in[0,2\pi)$.
Imposing that $H^{\textrm{eff}}$ is diagonal when expressed in terms of
$\alpha_{\boldsymbol{k}}$ and $\beta_{\boldsymbol{k}}$, that is,
$H^{\textrm{eff}}=\sum_{\boldsymbol{k}}\epsilon_{\alpha,\boldsymbol{k}}\alpha_{\boldsymbol{k}}^{\dagger}\alpha_{\boldsymbol{k}}+\epsilon_{\beta,\boldsymbol{k}}\beta_{\boldsymbol{k}}^{\dagger}\beta_{\boldsymbol{k}}$,
we find the following dispersion relation up to an overall constant factor:

$\epsilon_{\alpha,\boldsymbol{k}}=\epsilon_{\beta,\boldsymbol{k}}=JS\left(3+D_{\boldsymbol{k}}+\frac{2K}{J}\right)\sqrt{1-\tanh^{2}\xi_{\boldsymbol{k}}}$
(11)

for

$\tanh^{2}\xi_{\boldsymbol{k}}=\left(\frac{\mathcal{J}_{0}\left(\lambda\right)}{3+D_{\boldsymbol{k}}+2K/J}\right)^{2}|\phi_{\boldsymbol{k}}|^{2},$
(12)

and

$\theta_{\boldsymbol{k}}=n\pi-\arctan\frac{\textrm{Im{\hbox{\phi_{\boldsymbol{k}}}}}}{\textrm{Re}\phi_{\boldsymbol{k}}},$
(13)

where $n=0,1$ is fixed through the condition
$\textrm{sign}\left[\cos\theta_{\boldsymbol{k}}\right]=\textrm{sign}\left[\mathcal{J}_{0}\left(\lambda\right)\right]$.

####  II.2.1 Energy dispersion

From the diagonalization procedure, we find that the two magnon modes of the
system are degenerate, i.e.
$\epsilon_{\alpha,\boldsymbol{k}}=\epsilon_{\beta,\boldsymbol{k}}$. This is
known to be the case for usual antiferromagntes, and we verify that the
introduction of the laser field does not change this aspect of the problem.
This degeneracy can, however, be broken in several manners, for example by
defining the easy axis anisotropy as being different for the two sublattices,
applying a magnetic field along the $z$-direction, or accounting for magnetic
dipolar interactions [44].

The effect of the applied laser field on the magnon dispersion is twofold. If
the laser energy is tuned such that $\hbar\omega_{l}\gg J$, then the synthetic
DMI essentially vanishes, $D_{\boldsymbol{k}}\approx 0$, and the laser
manifests itself through the term $\mathcal{J}_{0}(\lambda)$ only. This term
is responsible for modifying the bandwidth of the magnon dispersion, which now
oscillates with steadily decaying amplitude as $\lambda$ increases, similarly
to the behavior of a Bessel function. Since $\phi_{\boldsymbol{k}}$ is maximal
at $\boldsymbol{k}=0$, and vanishes at the Dirac points, the effect of
$\mathcal{J}_{0}(\lambda)$ is more pronounced at the center of the Brillouin
zone, and does nothing at its vertices, $\boldsymbol{K}_{\pm}$.

Figure 2: (a) Energy dispersion
$\epsilon_{\boldsymbol{k}}\equiv\epsilon_{\alpha,\boldsymbol{k}}=\epsilon_{\beta,\boldsymbol{k}}$
along the $k_{y}=0$ line for three different values of $\lambda$. (b)
Dependence of the $\boldsymbol{k}=0$ magnon mode energy,
$\epsilon_{\boldsymbol{0}}$, with $\lambda$. In both panels we considered
$J=1$meV, $K=10^{-4}J$ and $\hbar\omega_{l}=10J$.

This behavior is illustrated in Fig. 2a where we depict the magnon dispersion
in momentum space along the line $k_{y}=0$ for three distinct values of
$\lambda$. In the same figure, we also plot the magnon dispersion at
$\boldsymbol{k}=0$ as a function of $\lambda$, perfectly illustrating the
aforementioned oscillating behavior. Furthermore, we also note that for a
judicious choice of $\lambda$, corresponding to the zeros of
$\mathcal{J}_{0}(\lambda)$, one obtains an almost flat band for the magnon
dispersion. This is easily understood from the inspection of Eq. (11), where
one sees that if $\mathcal{J}_{0}(\lambda)=0$, then the square root term
simply equals one, and the momentum dependence comes only from the term
$D_{\boldsymbol{k}}$, which can be made arbitrarily small by increasing the
energy of the laser.

Consider now the case where the laser’s energy is reduced, such that the
synthetic DMI becomes relevant. In this limit, time-reversal symmetry is
broken, and the dispersion relation is no longer even in momentum space,
meaning that
$\epsilon_{\alpha,\boldsymbol{k}}\neq\epsilon_{\alpha,-\boldsymbol{k}}$. Since
$D_{\boldsymbol{k}=0}=0$, the effect of the DMI is mainly noticeable near the
Dirac points, where it induces an energy difference between the magnon
dispersion at $\boldsymbol{K}_{+}$ and $\boldsymbol{K}_{-}$. This effect is
visible in Fig. 2a, where we observe that as
$\boldsymbol{k}\rightarrow\boldsymbol{K}_{+}$ the magnon dispersion increases
when compared with the case where the laser is turned off $(\lambda=0$); the
opposite statement is valid when
$\boldsymbol{k}\rightarrow\boldsymbol{K}_{-}$. The presence of the DMI is not
guaranteed to always yield this effect, since it depends on the value of
$\mathcal{V}_{\textrm{DM}}$ which may be positive, negative, or zero,
depending on the value of $\lambda$.

####  II.2.2 Bogoliubov coefficients

Now that the energy dispersion of the system was studied, let us discuss the
Bogoliubov transformation itself, and how its coefficients depend on the laser
field.

We start by noting that the definitions we gave for $\alpha_{\boldsymbol{k}}$
and $\beta_{\boldsymbol{k}}$ in terms of $a_{\boldsymbol{k}}$ and
$b_{-\boldsymbol{k}}^{\dagger}$ could be alternatively expressed in terms of a
two mode squeeze operator
$S_{2}(\zeta_{\boldsymbol{k}})=\exp\left(\zeta_{\boldsymbol{k}}a_{\boldsymbol{k}}b_{-\boldsymbol{k}}-\zeta_{\boldsymbol{k}}^{*}a_{\boldsymbol{k}}^{\dagger}b_{-\boldsymbol{k}}^{\dagger}\right)$,
where
$\zeta_{\boldsymbol{k}}=\left(\xi_{\boldsymbol{k}}/2\right)e^{-i\theta_{\boldsymbol{k}}}$
[13, 14]. Using this operator, we could have defined
$\alpha_{\boldsymbol{k}}=S_{2}(\zeta_{\boldsymbol{k}})a_{\boldsymbol{k}}S_{2}^{\dagger}(\zeta_{\boldsymbol{k}})$
and
$\beta_{\boldsymbol{k}}=S_{2}(\zeta_{\boldsymbol{k}})b_{-\boldsymbol{k}}S_{2}^{\dagger}(\zeta_{\boldsymbol{k}})$
(the proof of this statement is given in the Supplementary Information).
Hence, based on these new definitions, we see that the magnonic excitations of
the system ($\alpha_{\boldsymbol{k}}$ and $\beta_{\boldsymbol{k}}$), which are
linear combinations of sublattice magnons ($a_{\boldsymbol{k}}$ and
$b_{-\boldsymbol{k}}$), correspond, in fact, to two-mode squeezed magnons,
with a squeezing parameter $\xi_{\boldsymbol{k}}$. Although the identification
of the eigenmodes of the antiferromagnet as being squeezed magnons is
independent of the external laser field, the presence of the laser introduces
new interesting features in the system due to the possibility of tuning both
$\xi_{\boldsymbol{k}}$ and $\theta_{\boldsymbol{k}}$, as we shall see below.

Let us now study how the parameters $\xi_{\boldsymbol{k}}$ and
$\theta_{\boldsymbol{k}}$ depend on the external laser field, and how their
values change with the momentum $\boldsymbol{k}$.

Figure 3: (a) Depiction of the parameters defining the Bogoliubov
transformation, $\xi_{\boldsymbol{k}}$ and $\theta_{\boldsymbol{k}}$, as a
function of $\lambda$ for the $\boldsymbol{k}=0$ mode. The dashed line
represents the Bessel function $\mathcal{J}_{0}(\lambda)$ and serves as a
guide to the eye. (b) Representation of the momentum dependence of
$\xi_{\boldsymbol{k}}$ along the $k_{y}=0$ line for three different values of
$\lambda$. The parameters $J=1$meV, $K=10^{-4}J$ and $\hbar\omega_{l}=10J$
were used.

In Fig. 3a we depict $\xi_{\boldsymbol{k}}$ for the $\boldsymbol{k}=0$ mode,
as a function of $\lambda$. Since we are considering only the isotropic magnon
mode, the DMI vanishes automatically, i.e. $D_{\boldsymbol{k}=0}=0$. For
$\lambda=0$, that is, when the laser field is turned off, we find a large
value for $\xi_{\boldsymbol{k}=0}$; in fact, its value would be even larger if
the easy axis anisotropy had been ignored, since in that case the squeezing
parameter would diverge in the limit $\boldsymbol{k}\rightarrow 0$. As
$\lambda$ increases the squeezing parameter decreases, until it vanishes when
$\lambda$ reaches the first zero of the Bessel function $\mathcal{J}_{0}$.
Afterwards, we find that $\xi_{\boldsymbol{k}=0}$ increases slightly and an
oscillatory behaviour sets in. It attains successive maxima when maximum or
minimum values of $\mathcal{J}_{0}$ are hit, and vanishes at the zeros of the
Bessel function. The maxima of $\xi_{\boldsymbol{k}=0}$ become progressively
smaller as $\lambda$. Performing an identical analysis for $\boldsymbol{k}\neq
0$ one finds a similar pattern, albeit with an overall smaller magnitude. When
$\boldsymbol{k}\approx\boldsymbol{K}_{\pm}$ we find
$\xi_{\boldsymbol{k}}\approx 0$, regardless of $\lambda$, since
$\phi_{\boldsymbol{k}}\approx 0$ near the vertices of the Brillouin zone.
Thus, we see that the magnon squeezing tends to decrease as the magnon
momentum increases, until it vanishes at the Dirac points, as depicted in Fig.
3b. Also depicted in Fig. 3a, is the dependence of $\theta_{\boldsymbol{k}}$
with $\lambda$ for the mode $\boldsymbol{k}=0$. When the laser is turned off,
i.e. $\lambda=0$, one finds $\theta_{\boldsymbol{k=}0}=0$. Interestingly, when
$\lambda$ crosses the first zero of the Bessel function $\mathcal{J}_{0}$, the
phase $\theta_{\boldsymbol{k}=0}$ jumps to $\pi$. This phase jumping proves to
be the crucial ingredient to enhance the magnon-photon coupling, which we
shall discuss in the next section. For modes with $\boldsymbol{k}\neq 0$ other
values are obtained for $\theta_{\boldsymbol{k}}$, in agreement with Eq. (3),
but the $\pi$-phase jumps remain.

##  III Magnon-photon coupling

In this section we consider that the antiferromagnet is placed inside an
optical cavity, and study how the squeezed magnons couple with the cavity
photons. The dependence of the magnon-photon coupling on the external laser
field will be studied, and the transmission spectrum through the cavity
computed.

###  III.1 System Hamiltonian

The first step to study the interaction of the squeezed magnons with the
cavity photons is to define the Hamiltonian of the system. This Hamiltonian is
composed of three distinct contributions

$H=H_{\textrm{cav}}+H_{\textrm{AFM}}+H_{\textrm{int}}$ (14)

where
$H_{\textrm{cav}}=\sum_{\boldsymbol{q}}\hbar\omega_{\boldsymbol{q}}p_{\boldsymbol{q}}^{\dagger}p_{\boldsymbol{q}}$
is the Hamiltonian of the cavity photons, with
$p_{\boldsymbol{q}}/p_{\boldsymbol{q}}^{\dagger}$ the annihilation/creation
operator of a photon with momentum $\boldsymbol{q}$ and energy
$\hbar\omega_{\boldsymbol{q}}=\hbar c|\boldsymbol{q}|$; $H_{\textrm{AFM}}$ is
the effective Hamiltonian we introduced in the previous section to describe
the antiferromagnet under the incidence of the external laser field, and
$H_{\textrm{int}}=g\mu_{B}\sum_{i\in
A}\boldsymbol{b}(\boldsymbol{r}_{i})\cdot\boldsymbol{S}_{A}(\boldsymbol{r}_{i})+g\mu_{B}\sum_{i\in
B}\boldsymbol{b}(\boldsymbol{r}_{i})\cdot\boldsymbol{S}_{B}(\boldsymbol{r}_{i})$
is the interaction Hamiltonian between the magnetic field of the cavity
photons, $\boldsymbol{b}$, and the spins of the antiferromagnet.

We now emphasize that up to this point all the momenta we considered were two-
dimensional, due to the in-plane nature of the magnons we are studying.
However, since the cavity photons carry a three dimensional momentum, we must
differentiate between the two. To that end, henceforth, when considering a
two-dimensional momentum we shall label it with the index $||$, indicating its
in-plane configuration (e.g. $\boldsymbol{q}_{||}$ corresponds to the in-plane
component of the three-dimensional momentum $\boldsymbol{q}$).

Considering the magnetic field of the cavity photons to be circularly
polarized, we write its quantized form as [8]

$\displaystyle\boldsymbol{b}^{\pm}\left(\boldsymbol{r}\right)=\frac{1}{c}\sum_{\boldsymbol{q}}\sqrt{\frac{\hbar\omega_{\boldsymbol{q}}}{4\epsilon_{0}V}}\bigg{[}$
$\displaystyle\left(p_{\boldsymbol{q}}e^{i\boldsymbol{q}\cdot\boldsymbol{r}}+p_{\boldsymbol{q}}^{\dagger}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}\right)\hat{\boldsymbol{x}}$
$\displaystyle\pm i$
$\displaystyle\left(p_{\boldsymbol{q}}e^{i\boldsymbol{q}\cdot\boldsymbol{r}}-p_{\boldsymbol{q}}^{\dagger}e^{-i\boldsymbol{q}\cdot\boldsymbol{r}}\right)\hat{\boldsymbol{y}}\bigg{]},$
(15)

where $\pm$ stands for the two possible circular polarizations, $c$ is the
speed of light, $\epsilon_{0}$ is the vacuum permittivity and $V$ the volume
of the cavity. To express $H_{\textrm{int}}$ in terms of the eigenmodes of the
antiferromagnet, that is, using the operators $\alpha_{\boldsymbol{k}}$ and
$\beta_{\boldsymbol{k}}$, we start by expressing
$\boldsymbol{S}_{A/B}(\boldsymbol{r}_{i})$ in terms of the sublattice magnon
operators, $a_{i}$ and $b_{i}$, through the linearized Holstein-Primakoff
relations given in the previous section. Then, the Fourier components of these
operators are introduced, and the Bogoliubov transformation inverted, in order
to express $a_{\boldsymbol{k}}$ and $b_{\boldsymbol{k}}$ in terms of
$\alpha_{\boldsymbol{k}}$ and $\beta_{\boldsymbol{k}}$. At last, dropping
terms with the product of two annihilation or two creation operators, one
finds

$\displaystyle H_{\textrm{int}}^{+}$
$\displaystyle=\sum_{\boldsymbol{q}}\mathcal{U}_{-\boldsymbol{q}}p_{\boldsymbol{q}}\beta_{-\boldsymbol{q}_{||}}^{\dagger}+\textrm{h.c.}$
(16) $\displaystyle H_{\textrm{int}}^{-}$
$\displaystyle=\sum_{\boldsymbol{q}}\mathcal{U}_{\boldsymbol{q}}p_{\boldsymbol{q}}\alpha_{\boldsymbol{q}_{||}}^{\dagger}+\textrm{h.c.}$
(17)

where, once again, the superscript $\pm$ refers to the two circular
polarization of the cavity photons, $\boldsymbol{q}_{||}$ refers to the in-
plane component of the 3D-momentum $\boldsymbol{q}$, and

$\mathcal{U}_{\boldsymbol{q}}=\frac{g\mu_{B}}{c}\sqrt{\frac{NS\hbar\omega_{\boldsymbol{q}}}{2\epsilon_{0}V}}\left(\cosh\frac{\xi_{\boldsymbol{\boldsymbol{q}}_{||}}}{2}-e^{i\theta_{\boldsymbol{\boldsymbol{q}}_{||}}}\sinh\frac{\xi_{\boldsymbol{\boldsymbol{q}}_{||}}}{2}\right)$
(18)

is the magnon-photon coupling. Notice how according to Eqs. (16) and (17) the
two orthogonal circular polarizations couple selectively with just one of the
magnon modes each, with equal coupling strength. Although in the system we are
considering the two magnon modes are degenerate, the polarization of the
cavity photons could be used to select a specific magnon branch in a system
where said degeneracy is broken. Also, we note that for a fixed energy of the
cavity photons (and thus for a fixed $\boldsymbol{q}$), the value of
$\boldsymbol{q}_{||}$ can be tuned by changing the relative orientation of the
antiferromagnet and the magnetic field of the cavity mode.

###  III.2 Coupling strength

Let us now study in more detail the magnon-photon coupling strength. The
coupling $\mathcal{U}_{\boldsymbol{q}}$ is composed of two distinct
contributions: (i) a numerical pre-factor
$\mathcal{A}_{\boldsymbol{q}}=\left(g\mu_{B}/c\right)\sqrt{NS\hbar\omega_{\boldsymbol{q}}/2\epsilon_{0}V}$
determined by the properties of the system, namely the cavity photon energy,
cavity volume and number of spins, and (ii) an additional multiplicative term
determined by the coefficients of the Bogoliubov transformation,
$f_{\boldsymbol{q}_{||}}=\cosh\frac{\xi_{\boldsymbol{\boldsymbol{q}}_{||}}}{2}-e^{i\theta_{\boldsymbol{\boldsymbol{q}}_{||}}}\sinh\frac{\xi_{\boldsymbol{\boldsymbol{q}}_{||}}}{2}$,
which depends on the external laser field. The latter contribution is the one
we are interested in studying, in particular how its modulus
$|f_{\boldsymbol{q}_{||}}|=\sqrt{\cosh\xi_{\boldsymbol{\boldsymbol{q}}_{||}}-\cos\theta_{\boldsymbol{\boldsymbol{q}}_{||}}\sinh\xi_{\boldsymbol{\boldsymbol{q}}_{||}}}$
is affected by the momentum dependence and the applied laser.

Consider first the case where the laser is turned off, $\lambda=0$. According
to Fig. 3a, we find $\theta_{\boldsymbol{q}_{||}}=0$, which leads to
$|f_{\boldsymbol{q}_{||}}|=e^{-\xi_{\boldsymbol{q}_{||}}/2}$. Following the
analysis of the previous section regarding the squeezing parameter, we know
that it takes its largest value for $\boldsymbol{q}_{||}=0$ and monotonically
decreases as the momentum increases, until it vanishes at the Dirac points
(where the magnons are no longer squeezed). Hence, in the absence of the
external laser, we find that the magnon-photon coupling decreases
exponentially as the magnon momentum approaches the center of the Brillouin
zone. In particular, for the parameters $J=1$meV, $K=10^{-4}J$ and
$\hbar\omega_{l}=10J$, we find
$|f_{\boldsymbol{q}_{||}=0}|/|f_{\boldsymbol{q}_{||}=\boldsymbol{K}_{\pm}}|\approx
0.05$, that is, the magnon-photon coupling is 20 times smaller near the
Brillouin zone center than at the Dirac points.

Figure 4: Plot of
$|f_{\boldsymbol{q}_{||}}|=\sqrt{\cosh\xi_{\boldsymbol{\boldsymbol{q}}_{||}}-\cos\theta_{\boldsymbol{\boldsymbol{q}}_{||}}\sinh\xi_{\boldsymbol{\boldsymbol{q}}_{||}}}$
as a function of $\lambda$ for three different momenta,
$\boldsymbol{q}_{||}=0$, $\boldsymbol{q}_{||}=\boldsymbol{K}_{+}/2$ and
$\boldsymbol{q}_{||}=3\boldsymbol{K}_{+}/4$, with
$\boldsymbol{K}_{+}=\left(4\pi/3\sqrt{3}a,0\right)$. The dashed line shows the
the evolution of $\xi_{\boldsymbol{\boldsymbol{q}}_{||}=0}$ with $\lambda$
(also shown in Fig. 3a). The areas shaded in blue/orange refer to the range of
$\lambda$ where the magnon-photon coupling is below/above its value at the
Dirac points when the laser is absent. The parameters $J=1$meV, $K=10^{-4}J$
and $\hbar\omega_{l}=10J$ were used.

Let us now consider that $\lambda\neq 0$, corresponding to the scenario where
the laser is turned on. Initially we have $\theta_{\boldsymbol{q}_{||}}=0$,
and as $\lambda$ increases the squeezing parameter $\xi_{\boldsymbol{q}_{||}}$
decreases, leading to an increase of the coupling strength
$|f_{\boldsymbol{q}_{||}}|$ when compared to the case where the laser is
absent. When $\lambda$ reaches the first zero of the Bessel function
$\mathcal{J}_{0}$, the squeezing parameter vanishes on the entire Brillouin
zone and
$|f_{\boldsymbol{q}_{||}}|=1\forall\boldsymbol{q}_{||}\in\textrm{1BZ}$. We
note, however, that a similar effect could be achieved without the laser by
increasing the easy axis anisotropy, which for a large enough value would
yield a similar result. The unique feature introduced by the laser arises when
the system is driven in such a way that $\lambda$ is located between the first
two zeroes of the $\mathcal{J}_{0}$. In that situation, and focusing on the
$\boldsymbol{q}_{||}=0$ mode, we find $\theta_{\boldsymbol{q}=0}=\pi$, leading
to a coupling strength which grows exponentially with
$\xi_{\boldsymbol{q}_{||}}$, that is
$|f_{\boldsymbol{q}_{||}}|=e^{\xi_{\boldsymbol{q}_{||}}/2}$. Since in the
considered range of values for $\lambda$ we have $\xi_{\boldsymbol{q}}>0$,
then $|f_{\boldsymbol{q}_{||}=0}|>1$. In fact, one finds that
$|f_{\boldsymbol{q}_{||}=0}|$ reaches a peak of approximately 1.25 when
$\lambda\approx 3.8$ (corresponding to the first minimum of
$\mathcal{J}_{0}$). The peak value $|f_{\boldsymbol{q}_{||}=0}|\approx 1.25$
is obtained in the limit where the exchange coupling $J$ dominates the
anisotropy term $K$. If $K$ becomes comparable with $J$, then the maximum
value for $|f_{\boldsymbol{q}_{||}=0}|$ decreases.

Hence, by driving the system with a laser field, we are able to enhance the
coupling of the $\boldsymbol{q}_{||}=0$ magnons with the cavity photons by a
factor of approximately 25 when compared to the case where no laser is
applied. Not only that, but we are also able to push the coupling strength
beyond the limit of what is found for unsqueezed magnons, whilst retaining
some magnon squeezing ($\xi_{\boldsymbol{q}_{||}=0}\approx 0.2$ for
$\lambda\approx 3.8$); a result only attainable by driving the system with the
laser field. For magnon modes with $\boldsymbol{q}_{||}\neq 0$ the laser
induced enhanced coupling is also present, although the effect becomes
progressively weaker as we approach the Dirac points. These results are
summarized in Fig. 4 where we depict $|f_{\boldsymbol{q}_{||}}|$ as a function
of $\lambda$ for different momenta $\boldsymbol{q}_{||}$. The shaded
blue/orange areas indicate the range of $\lambda$ where the magnon-photon
coupling is smaller/larger than for unsqueezed magnons.

###  III.3 Transmission spectra

Now that the Hamiltonian of the system has been determined, and the details of
the magnon-photon coupling have been discussed, we move on to the computation
of the transmission spectra through the cavity. Following [45], we state that
the transmission through the cavity should be proportional to the spectral
function of the cavity photons, which follows directly from the retarded
Green’s function $G^{R}(t)=-\frac{i}{\hbar}\theta(t)\langle
p_{\boldsymbol{q}}(t)p_{\boldsymbol{q}}^{\dagger}(0)\rangle$. Considering now
that the cavity has only one mode, and that its magnetic field is polarized
with the “$-$” circular polarization, we obtain the following expression for
the retarded Green’s function [45, 46]

$G^{R}(\hbar\omega)=\frac{1}{\hbar\omega-\hbar\omega_{\boldsymbol{q}}+i\Gamma-\frac{|\mathcal{U}_{\boldsymbol{q}}|^{2}}{\hbar\omega-\epsilon_{\boldsymbol{q}_{||}}+i\delta}},$
(19)

where $\Gamma$ and $\delta$ are the cavity and antiferromagnet intrinsic
damping factors, and $\hbar\omega$ is the energy of the incident photons. Once
more we have $\hbar\omega_{\boldsymbol{q}}=\hbar c|\boldsymbol{q}|$ the energy
of the cavity mode and $\epsilon_{\boldsymbol{q}_{||}}$ the energy of the
magnon mode that couples with the cavity photons. Although the circular
polarization we chose only allows coupling to the $\alpha-$magnons, Eq. (17),
a similar result would be obtained for the other circular polarization where
coupling to $\beta-$magnons would appear instead. Since for the system we are
considering the two magnon modes are degenerate, the transmission spectrum is
identical in both cases.

Having determined $G^{R}(\hbar\omega_{\boldsymbol{q}})$ we define the
transmission amplitude as
$t(\hbar\omega)\propto-2\textrm{Im}G^{R}(\hbar\omega)$. A simple way to avoid
the proportionality relation, is to normalize the transmission amplitude by
its value when the antiferromagnet is not present in the cavity
($\mathcal{U}_{\boldsymbol{q}}=0$) and the incident photons are in resonance
with the cavity mode $\hbar\omega=\hbar\omega_{\boldsymbol{q}}$. Doing so, we
find

$\bar{t}(\hbar\omega)=-\textrm{Im}\frac{\Gamma}{\hbar\omega-\hbar\omega_{\boldsymbol{q}}+i\Gamma-\frac{|\mathcal{U}_{\boldsymbol{q}}|^{2}}{\hbar\omega-\epsilon_{\boldsymbol{q}_{||}}+i\delta}},$
(20)

where $\bar{t}(\hbar\omega)$ is the normalized transmission amplitude. The
poles of $\bar{t}(\hbar\omega)$ correspond to the excitations of the system,
which are termed magnon-polaritons as they correspond to the hybridization of
the antiferromagnet magnons with the cavity photons.

Figure 5: Density plots of the normalized transmission $\bar{t}(\hbar\omega)$
as a function of the incident photon’s energy and the in-plane magnon momentum
(along the line with $k_{y}=0$). In each plot a different value for $\lambda$
was considered and the energy of the cavity mode (horizontal dashed line) was
chosen to match the minimum of the magnon dispersion (curved dashed line). The
parameters $J=1$meV, $K=10^{-4}J$ and $\hbar\omega_{l}=50J$ were used. We also
considered $\Gamma=\delta=10^{-3}J$, and $\mathcal{A}_{\boldsymbol{q}}=8\times
10^{-4}\sqrt{\hbar\omega_{\boldsymbol{q}}}$ (in eV).

In Fig. 5 we depict a density plot of $\bar{t}(\hbar\omega)$ as a function of
the incident photon energy $\hbar\omega$ and the magnon momentum
$\boldsymbol{k}_{||}$, along the line $k_{y}=0$ on the first Brillouin zone,
for three different values of $\lambda$; for each case the energy of the
cavity mode was chosen to match the minimum of the magnon dispersion. The
first thing we note regarding this result is that the bandwidth of the magnon
dispersion is modified by changing $\lambda$ in agreement with what was
previously found in Fig. 2 (the effect of the DMI is small since we are
considering $\hbar\omega_{l}=50J$). Focusing on the leftmost panel, we see
that the hybridization between the cavity mode and the zero-momentum magnon is
almost nonexistent; this is simultaneously due to the small energy of the
cavity photon (leading to a small $\mathcal{A}_{\boldsymbol{k}}$), and due to
the fact that $|f_{\boldsymbol{k}_{||}}|$ is small for $\lambda=0.5$.
Regarding the middle panel, we see that, due to the higher value used for
$\lambda$, the bottom of the magnon dispersion moves to higher energies, and,
at the same time, the magnitude of $|f_{\boldsymbol{k}_{||}}|$ increases (see
Fig. 4). Because of this, one finds a much clearer coupling between the
magnons and the cavity photons, which is identified by the anticrossing of the
magnon-polariton branches. At last, on the rightmost panel, where
$\lambda=3.8$, we observe the largest hybridization between magnons and
photons, with a clear separation between the two magnon-polariton modes. This
strong coupling was expected, since $\lambda$ was carefully chosen to maximize
$|f_{\boldsymbol{k}_{||}}|$, in agreement with the discussion following Fig.
4.

##  IV Plausible parameters

Up to this point we have discussed the general physical properties of a laser
driven antiferromagnet without focusing on a particular implementation of the
system. In this section we wish to estimate the values the different
parameters the system should have in order to observe the phenomena discussed
so far. To the best of our knowledge, an experiment focused on the phenomena
we describe in this paper has not been performed so far.

To do so, we start by setting $\lambda=2.4$, corresponding to the minimum
value of $\lambda$ which allows an enhanced magnon-photon coupling above the
unsqueezed limit (see Fig. 4). Recalling $\lambda=g\mu_{B}E_{0}a/\hbar c^{2}$,
and setting $g=2$, we find that in order to have $\lambda=2.4$ we must have
$E_{0}a\approx 1.2\times 10^{6}$ V. For $a=3.3\text{\r{A}}$ [47],
corresponding to a value in the order of magnitude that is usually found for
hexagonal boron nitride, transition metal dichalcogenides, etc., we obtain
$E_{0}\approx 3\times 10^{13}\textrm{V/cm}$. If a larger $a$ is used, a
smaller electric field is required to produce the same $\lambda$. In Ref. [48]
a laser intensity of $10^{23}\textrm{W/cm}^{2}$ has been reported, which in
terms of the electric field’s magnitude corresponds to roughly $E_{0}\approx
9\times 10^{12}\textrm{V/cm}$, within the order of magnitude of the required
value. Hence, current state of the art laser technology is already compatible
with the requirements to significantly increase the magnon-photon coupling. To
estimate the numerical value of $\mathcal{U}_{\boldsymbol{q}}$ we focus on the
magnitude of
$\mathcal{A}_{\boldsymbol{q}}=\left(g\mu_{B}/c\right)\sqrt{NS\hbar\omega_{\boldsymbol{q}}/2\epsilon_{0}V}$.
We now set $S=1/2$ and note that $N/V=\rho/2h$ where $\rho/2$ is the density
of spins and $h$ is the height of the cavity. Considering $h=40\textrm{nm}$
[49] and $\rho/2\sim 10^{14}\textrm{cm}^{-2}$ (in agreement with the measured
values of [47], and compatible with the value of $a$ considered above), we
find $\mathcal{A}_{\boldsymbol{q}}\approx
10^{-4}\sqrt{\hbar\omega_{\boldsymbol{q}}}$ (in eV), which is comparable with
what we used in Fig. 5. For $\hbar\omega\sim 1\textrm{meV}$ we find
$\mathcal{A}_{\boldsymbol{q}}\sim 3\times 10^{-6}$ eV, corresponding to
approximately 0.7 GHz, which is well within the resolution of current
experiments [8, 50, 51].

##  V Final remarks

In this paper we studied a honeycomb antiferromagnet subject to an external
laser field. After a Floquet effective Hamiltonian was introduced, we
diagonalized it using a Bogoliubov transformation, and found that the
eigenexcitations of the system correspond to squeezed magnons. The presence of
the laser allows the control of several properties of these spin excitations,
namely their dispersion relation, and squeezing parameters.

When studying how the magnons couple with photons in a cavity, we found that
the two magnon modes of the system couple selectively with the two orthogonal
circular polarizations. Furthermore, we found that by tuning the intensity of
the applied laser field the strength of the magnon-photon coupling can be
significantly enhanced. In particular, we found that, for the isotropic magnon
mode, the coupling strength with the cavity photons can be enhanced by an
order of magnitude when compared to what is found in the absence of the laser.
When the system is not driven, the largest magnon-photon coupling is achieved
at the Dirac points, where the magnons do not present squeezing. When the
laser is present, however, the maximum coupling is reached for the magnons
near the Brillouin zone center, which preserve some of the original squeezing;
in this case the magnon-photon coupling surpasses the maximum value found for
the undriven system.

To the best of our knowledge, an experiment probing these physical effects has
not yet been realized. However, from our estimates, state of the art equipment
may be able to probe the discussed phenomena. Since our theoretical
predictions were derived for a strictly 2D monolayer, the results may differ
if experiments are realized on thin films, due to the sensitivity of these
materials regarding the number of layers.

Finally, we note that although we focused on the magnon-photon interaction,
other options could be explored. An example of this is the interaction of
magnons with conduction electrons. According to Ref. [18] the transition rate
for the electron-magnon scattering takes a similar form to the function
$f_{\boldsymbol{q}_{||}}$ we identified in the magnon-photon coupling. Due to
this similarity, the ideas discussed here regarding the coupling strength
could, in principle, also be applied in that type of system.

###### Acknowledgements.

J.C.G.H. acknowledges the Center of Physics for a grant funded by the
UIDB/04650/2020 strategic project. T. V. C. A. acknowledges funding through
the project QML-HEP - CERN/FIS-COM/0004/2021 . N. M. R. P. acknowledges
financial support from FCT (grant PTDC/FIS-MAC/2045/2021) and from the
European Commission through the project GrapheneDriven Revolutions in ICT and
Beyond (Ref. No. 881603, CORE 3).

## References

  * Chumak  _et al._ [2015] A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Nature Physics 11, 453 (2015).
  * Yuan  _et al._ [2021] H. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan, arXiv preprint arXiv:2111.14241 (2021).
  * Chumak  _et al._ [2014] A. V. Chumak, A. A. Serga, and B. Hillebrands, Nature communications 5, 1 (2014).
  * Demidov  _et al._ [2011] V. Demidov, M. Kostylev, K. Rott, J. Münchenberger, G. Reiss, and S. Demokritov, Applied Physics Letters 99, 082507 (2011).
  * Khitun  _et al._ [2010] A. Khitun, M. Bao, and K. L. Wang, Journal of Physics D: Applied Physics 43, 264005 (2010).
  * Schneider  _et al._ [2008] T. Schneider, A. A. Serga, B. Leven, B. Hillebrands, R. L. Stamps, and M. P. Kostylev, Applied Physics Letters 92, 022505 (2008).
  * Soykal and Flatté [2010] Ö. O. Soykal and M. Flatté, Physical review letters 104, 077202 (2010).
  * Harder  _et al._ [2016] M. Harder, L. Bai, C. Match, J. Sirker, and C. Hu, Science China Physics, Mechanics & Astronomy 59, 1 (2016).
  * Harder  _et al._ [2021] M. Harder, B. Yao, Y. Gui, and C.-M. Hu, Journal of Applied Physics 129, 201101 (2021).
  * Tabuchi  _et al._ [2015] Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Science 349, 405 (2015).
  * Liu  _et al._ [2019] Z.-X. Liu, H. Xiong, and Y. Wu, Physical Review B 100, 134421 (2019).
  * Bittencourt  _et al._ [2019] V. A. Bittencourt, V. Feulner, and S. V. Kusminskiy, Physical Review A 100, 013810 (2019).
  * Gerry  _et al._ [2005] C. Gerry, P. Knight, and P. L. Knight,  _Introductory quantum optics_ (Cambridge University Press, New York, 2005).
  * Fox [2006] M. Fox,  _Quantum optics: an introduction_ , Vol. 15 (Oxford University Press Inc., New York, 2006).
  * Andersen  _et al._ [2016] U. L. Andersen, T. Gehring, C. Marquardt, and G. Leuchs, Physica Scripta 91, 053001 (2016).
  * Aasi  _et al._ [2013] J. Aasi, J. Abadie, B. Abbott, R. Abbott, T. Abbott, M. Abernathy, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Nature Photonics 7, 613 (2013).
  * Walls [1983] D. F. Walls, nature 306, 141 (1983).
  * Kamra  _et al._ [2019] A. Kamra, E. Thingstad, G. Rastelli, R. A. Duine, A. Brataas, W. Belzig, and A. Sudbø, Physical Review B 100, 174407 (2019).
  * Kamra and Belzig [2016] A. Kamra and W. Belzig, Physical review letters 116, 146601 (2016).
  * Liensberger  _et al._ [2019] L. Liensberger, A. Kamra, H. Maier-Flaig, S. Geprägs, A. Erb, S. T. Goennenwein, R. Gross, W. Belzig, H. Huebl, and M. Weiler, Physical review letters 123, 117204 (2019).
  * Zou  _et al._ [2020] J. Zou, S. K. Kim, and Y. Tserkovnyak, Physical Review B 101, 014416 (2020).
  * Erlandsen  _et al._ [2019] E. Erlandsen, A. Kamra, A. Brataas, and A. Sudbø, Physical Review B 100, 100503 (2019).
  * Kamra  _et al._ [2020] A. Kamra, W. Belzig, and A. Brataas, Applied Physics Letters 117, 090501 (2020).
  * Li  _et al._ [2019] J. Li, S.-Y. Zhu, and G. Agarwal, Physical Review A 99, 021801 (2019).
  * Yang  _et al._ [2021] Z.-B. Yang, H. Jin, J.-W. Jin, J.-Y. Liu, H.-Y. Liu, and R.-C. Yang, Physical Review Research 3, 023126 (2021).
  * Zhang  _et al._ [2021] W. Zhang, D.-Y. Wang, C.-H. Bai, T. Wang, S. Zhang, and H.-F. Wang, Optics Express 29, 11773 (2021).
  * Hirosawa  _et al._ [2022] T. Hirosawa, J. Klinovaja, D. Loss, and S. A. Díaz, Physical Review Letters 128, 037201 (2022).
  * Owerre [2017] S. Owerre, Journal of Physics Communications 1, 021002 (2017).
  * Kim  _et al._ [2016] S. K. Kim, H. Ochoa, R. Zarzuela, and Y. Tserkovnyak, Physical review letters 117, 227201 (2016).
  * Aharonov  _et al._ [1988] Y. Aharonov, P. Pearle, and L. Vaidman, Physical Review A 37, 4052 (1988).
  * Meier and Loss [2003] F. Meier and D. Loss, Physical review letters 90, 167204 (2003).
  * Elyasi  _et al._ [2019] M. Elyasi, K. Sato, and G. E. Bauer, Physical Review B 99, 134402 (2019).
  * Kar and Basu [2018] S. Kar and B. Basu, Physical Review B 98, 245119 (2018).
  * Owerre [2019a] S. Owerre, Annals of Physics 406, 14 (2019a).
  * Owerre  _et al._ [2019] S. Owerre, P. Mellado, and G. Baskaran, EPL (Europhysics Letters) 126, 27002 (2019).
  * Owerre [2019b] S. Owerre, Scientific Reports 9, 1 (2019b).
  * Proskurin  _et al._ [2019] I. Proskurin, R. Macêdo, and R. L. Stamps, New Journal of Physics 21, 095003 (2019).
  * Vinas Boström  _et al._ [2020] E. Vinas Boström, M. Claassen, J. McIver, G. Jotzu, A. Rubio, and M. Sentef, SciPost Physics 9, 061 (2020).
  * Sentef  _et al._ [2015] M. Sentef, M. Claassen, A. Kemper, B. Moritz, T. Oka, J. Freericks, and T. Devereaux, Nature communications 6, 1 (2015).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, J. M. Zeuner, Y. Plotnik, Y. Lumer, D. Podolsky, F. Dreisow, S. Nolte, M. Segev, and A. Szameit, Nature 496, 196 (2013).
  * Wang  _et al._ [2013] Y. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013).
  * Lindner  _et al._ [2011] N. H. Lindner, G. Refael, and V. Galitski, Nature Physics 7, 490 (2011).
  * Pires [2021] A. S. T. Pires,  _Theoretical Tools for Spin Models in Magnetic Systems_ (IOP Publishing, 2021).
  * Kamra  _et al._ [2017] A. Kamra, U. Agrawal, and W. Belzig, Physical Review B 96, 020411 (2017).
  * Harder [2018] M. Harder,  _Cavity spintronics: foundations and applications of spin-photon hybridization_ , Ph.D. thesis, Faculty of Graduate Studies of The University of Manitoba (2018).
  * Mahan [2000] G. D. Mahan,  _Many-Particle Physics_ , 3rd ed. (Springer Science+Business Media, New York, 2000).
  * Gong  _et al._ [2018] S.-J. Gong, C. Gong, Y.-Y. Sun, W.-Y. Tong, C.-G. Duan, J.-H. Chu, and X. Zhang, Proceedings of the National Academy of Sciences 115, 8511 (2018).
  * Yoon  _et al._ [2021] J. W. Yoon, Y. G. Kim, I. W. Choi, J. H. Sung, H. W. Lee, S. K. Lee, and C. H. Nam, Optica 8, 630 (2021).
  * Noginov  _et al._ [2009] M. Noginov, G. Zhu, A. Belgrave, R. Bakker, V. Shalaev, E. Narimanov, S. Stout, E. Herz, T. Suteewong, and U. Wiesner, Nature 460, 1110 (2009).
  * Harder  _et al._ [2018] M. Harder, Y. Yang, B. Yao, C. Yu, J. Rao, Y. Gui, R. Stamps, and C.-M. Hu, Physical review letters 121, 137203 (2018).
  * Yao  _et al._ [2019] B. Yao, T. Yu, X. Zhang, W. Lu, Y. Gui, C.-M. Hu, and Y. M. Blanter, Physical Review B 100, 214426 (2019).


* [17] X. Yuan, S. Zhang, and W. Zhang, _The gross-zagier formula on shimura curves_ , No. 184. Princeton University Press, 2013.


11institutetext: Webroot 111Webroot is an OpenText company.
11email<EMAIL_ADDRESS>
# Classifying Malware Using Function Representations in a Static Call Graph

Thomas Dalton  Mauritius Schmidtler  Alireza Hadj Khodabakhshi

###### Abstract

We propose a deep learning approach for identifying malware families using the
function call graphs of x86 assembly instructions. Though prior work on static
call graph analysis exists, very little involves the application of modern,
principled feature learning techniques to the problem. In this paper, we
introduce a system utilizing an executable’s function call graph where
function representations are obtained by way of a recurrent neural network
(RNN) autoencoder which maps sequences of x86 instructions into dense, latent
vectors. These function embeddings are then modeled as vertices in a graph
with edges indicating call dependencies. Capturing rich, node-level
representations as well as global, topological properties of an executable
file greatly improves malware family detection rates and contributes to a more
principled approach to the problem in a way that deliberately avoids tedious
feature engineering and domain expertise. We test our approach by performing
several experiments on a Microsoft malware classification data set and achieve
excellent separation between malware families with a classification accuracy
of 99.41%.

###### Keywords:

neural networks representation learning malware detection function call graph
reverse engineering

##  1 Introduction

Malware is often classified into families based on certain shared
characteristics between samples. It is often very useful to distinguish
between malware families in order to detect trends in malware infections over
time and to attribute authorship. Traditionally, classifying malware has
required teams of threat researchers to perform advanced reverse engineering
techniques in order to identify various unique characteristics that define a
family. However, cyber threats have exploded in recent years making it
difficult for threat researchers to keep up. Malware in particular continues
to grow in sophistication with new strains released daily. The practice of
malware polymorphism renders traditional automated signature-based approaches
ineffective for identifying novel instances of malware.

Figure 1: Variable-length sequences of x86 instructions found in functions are
embedded into fixed-length vectors using a GRU-based sequence-to-sequence
autoencoder.

In this work, we propose a new approach to malware classification that is
inspired by reverse engineering techniques yet requires no domain-specific
feature engineering and is invariant to polymorphism. Specifically, we devise
a function call-graph framework in which the function representations are
learned. By framing the problem through the lens of representation learning,
we are able to greatly improve automatic classification while also
contributing to human insight which is helpful to determine authorship and
intent. By incorporating rich, node-level representations as well as global,
structural properties of an executable’s call graph, we are able to classify
malware families with very high accuracy. Our approach consists of several
composite models aimed at learning robust function representations that when
employed together form the full classification system.

##  2 Related Work

Call graphs are commonly used by malware analysts and reverse engineers to
manually analyze and inspect executable files [14]. Indeed, many real-world
data is naturally represented using graphs. Graphs have been successfully
utilized in analyzing data from a wide variety of domains including social
network link prediction [11], protein-protein interactions [1], and
communication networks [12]. Due to their expressive ability, there is growing
interest in applying machine learning techniques directly to graph-represented
data to bypass tedious feature engineering. Graph kernels have been proposed
to allow for kernel-based methods (such as support vector machines) to be
applied directly to graph classification problems. Kernels based on the
Weisfeiler-Lehman test of graph isomorphism have grown in popularity in recent
years [21] owing to their relative simplicity and strong discriminative
ability.

In prior works, call graphs have been used to automatically classify malware
but typically these works employ relatively simple graph similarity measures
such as graph edit distance or rely on heavy feature engineering involving
summary statistics to describe functions in the graph [9, 20, 6, 5]. We build
on this call graph approach by incorporating certain representation learning
techniques such as autoencoding and clustering [7] to obtain an improved
function representation. By extending the well-established call graph
strategies with a principled representation learning approach, we forego the
tedious and heuristic feature engineering steps of prior work, giving us much
better graph representations.

Figure 2: $G_{0}$ and $G_{1}$ are considered similar since they both share a
similar graph topology and near-identical vertex label sets. On the other
hand, $G_{2}$ is considered very different from both $G_{0}$ and $G_{1}$ since
it neither shares a graph topology nor a vertex set with the other graphs. The
notion of graph “similarity” is clarified in Section 3.4

##  3 Call Graph Framework

A call graph describes the logical control flow of an executable file where
functions (or subroutines) are expressed as vertices and edges represent a
dependency or call relationship. Call graphs have proven to be extremely
useful to security researchers for the analysis and classification of
malicious software.

Actual | Predicted
---|---
mov edi, edi | mov edi, edi
push ebp | push ebp
mov ebp, esp | mov ebp, esp
sub esp, 48 | sub esp, 24
lea ecx, 0x | lea ecx, 0x
push eax | <unknown>
push ecx | push ecx
call <addr> | call <addr>
pop ecx | pop ecx
… | …
… | …
Figure 3:  An example sequence decoding given a latent embedding. The decoder
is able to re-create the original sequence with high accuracy, indicating that
the latent embedding has captured sufficient information.

The main intuition behind the call graph approach to malware classification is
that files sharing similar call graphs are likely to have been generated from
the same family. By understanding the logical flow of the executable, we can
gain significant insight into the intent of the malware. It is important,
therefore, to represent the graph such that we capture rich vertex-level
information as well as global, topological properties of the graph. The
intuition is illustrated in Figure 2.

###  3.1 Overview

We break the malware classification task down into three primary subtasks
which we summarise here. In order to obtain a good whole-graph representation
of the executable, it is important to first obtain high quality embeddings for
the functions contained within the file. For this, we use a sequence-to-
sequence autoencoder which captures the sequential nature of the x86 code
instructions into a low-dimensional, latent representation of the function.

This function embedding helps to make our model more robust to polymorphic
techniques since a perturbation in the x86 instruction space results in a
proportional perturbation in the embedding space. It can also be useful for
identifying the specific functions that make the file malicious. This function
embedding approach is one of the key differentiators between our approach and
prior call graph approaches to malware classification.

Having obtained function representions, we cluster the embeddings to obtain
discrete labels for the functions and re-label the graph vertices according to
their respective cluster IDs. Finally, the whole-graph representation is
obtained using a graph kernel inspired by the Weisfeiler-Lehman test of graph
isomorphism. The message-passing property of the Weisfeiler-Lehman framework
allow us to efficiently capture the global structure of the graph.

The executable binary files are disassembled into plain text .asm files using
IDA [10], a popular disassembler widely used by security researchers. Due to
the tendency of code sections to contain very long sequences (sometimes
upwards of hundreds of thousands of instructions), we break up the sequences
into functions or subroutines which provide natural delimiters much like
sentences and paragraphs are in a document. These shorter length sequences
enable us to use recurrent neural units such as long short-term memory (LSTM)
or gated recurrent units (GRU) where the training samples are individual
functions with sequences ranging from very short (fewer than 5 instructions)
to as long as a few hundred instructions. By following the call instructions
in the code, we can construct the file’s call graph,
$G=(\mathcal{V},\mathcal{E})$ where vertices $v\in\mathcal{V}$ represent
functions and edges $e=(v,v^{\prime})\in\mathcal{V}\times\mathcal{V}$
represent a call dependency.

Figure 4: The sequence-to-sequence autoencoder architecture. The left side of
the network encodes the sequence into a fixed-length latent representation.
The right side of the network decodes the latent representation back into the
original sequence.

The model considers only portions of the executable containing valid x86 code
instructions. Traditionally, these code sections are identified by their
section header name (e.g. .text, .code, etc.) but malware often obfuscates
intent by using a packer which may result in non-standard section names such
as .brick or iuagwws. In addition to all code found in the standard code
sections, our approach also considers such nonstandard sections containing
valid x86 code to construct the call graph.

We make a distinction between two kinds of vertices: internal functions and
external functions. Internal functions are those that are present in the
executable and subsequently can be disassembled directly. External functions
are those which are imported from external libraries and thus the code is not
readily available for disassembly. Our graph therefore consists of edges
between both internal and external functions. It is worth noting that
executable files also contain sections that do not typically contain code such
as .data or .reloc. While these sections provide additional data that is often
quite useful for malware classification, we ignore any non-code data for the
call graph construction task. It is possible to attribute the graph with the
information contained in such non-code sections but that is beyond the scope
of this work as our principle concern is that of malware classification using
x86 code representations.

A sequence-to-sequence [23] GRU-based [3] autoencoder architecture was chosen
for the task of embedding variable-length sequences of x86 code instructions
into fixed-length, continuous vectors. The function embedding model is
comprised of an encoder and a decoder with the encoder being responsible for
compressing sequences into low-dimensional, latent representations. A decoder
is used to decompress the fixed-length vector back into the original variable-
length sequence. Because the autoencoder model must recreate the original
sequence from its bottleneck representation, the model learns an efficient,
latent representation of the original sequence. After the autoencoder model is
trained, the decoder is discarded and only the encoder portion is used to
encode new sequences. Sequence-to-sequence architectures have been used
successfully in language modeling tasks such as machine translation. Often,
the goal is to translate a sequence of words from one language, such as
English, into another language, such as French.

Figure 5: For each iteration of the Weisfeiler-Lehman algorithm, a new label
is generated for each vertex. The new label is derived by concatenating the
labels of the adjacent vertices.

The input sequence is encoded into a bottleneck layer which captures a latent
representation of the sequence irrespective of language. In machine
translation tasks, the input and output sequences are usually composed of
words drawn from disjoint vocabularies. However, our sequence-to-sequence task
involves reconstructing the original input sequence from the bottleneck
representation, so the same vocabulary is used for both the input and the
output sequences. By reconstructing the original input sequence from the
latent representation, the sequence-to-sequence network becomes an
autoencoder.

###  3.2 Function Representations

Given a GRU with the following definitions,

$\displaystyle z_{t}$
$\displaystyle=\sigma_{g}(W_{z}x_{t}+U_{z}h_{t-1}+b_{z})$ $\displaystyle
r_{t}$ $\displaystyle=\sigma_{g}(W_{r}x_{t}+U_{r}h_{t-1}+b_{r})$
$\displaystyle h_{t}$ $\displaystyle=(1-z_{t})\circ
h_{t-1}+z_{t}\circ\sigma_{h}(W_{h}x_{t}+U_{h}(r_{t}\circ h_{t-1})+b_{h})$

the hidden representation of a sequence with length $T$ is taken to be
$h_{T}$.

During training, when the encoder receives the last token in the sequence
$x_{T}=$ <end>, the decoder is initialized with $h_{T}^{(enc)}$, thus
transferring the compressed sequence information to the decoder. In addition
to the final hidden state of the encoder, the decoder is also supplied with
the original input sequence with a one-step delay. That is, at time step $t$,
the decoder receives the true input $x_{t-1}$ for $t>1$ where $x_{0}=$
<start>. This technique of supplying the original sequence with a delay into
the decoder is known as teacher forcing.

The decoder, therefore, is trained to predict the next token of the sequence
given the hidden state of the encoder and the previous time step. This helps
to greatly speed up the training of the autoencoder. An example decoding is
illustrated in Figure 3. After training is completed, the decoder portion of
the model is thrown away and only the encoder is used to obtain latent
representations for function sequences.

It is common for authors of malware to obfuscate the intent of a file by
adding junk instructions such as no-op instructions. Because the function is
represented in a latent space, it is relatively immune to such common
obfuscation tactics which can often thwart signature-based or count-based
solutions.

In our model, input sequences are tokenized as the full x86 assembly
instruction after replacing relative memory address locations with an <addr>
string. The vocabulary is thus composed of the top 20,000 most common x86
instructions found in the corpus plus the <start> and <end> tokens.

Input: $G$, $G^{\prime}$, $h$

Output: $k^{(h)}(G,G^{\prime})$

for  _$i\leftarrow 1$ to $h$_ do

for  _$v\in\mathcal{V}$_ do

$N_{v}=neighbors(v)$

$S_{v}=sort(N_{v})$

$v,\sigma_{v}=hash(S_{v})$

end for

for  _$v^{\prime}\in\mathcal{V^{\prime}}$_ do

$N_{v^{\prime}}=neighbors(v^{\prime})$

$S_{v^{\prime}}=sort(N_{v^{\prime}})$

$v^{\prime},\sigma_{v^{\prime}}=hash(S_{v^{\prime}})$

end for

end for

$\phi(G)=\left[c(G,\sigma_{0}),\dots,c(G,\sigma_{s})\right]$

$\phi(G^{\prime})=\left[c(G^{\prime},\sigma_{0}),\dots,c(G^{\prime},\sigma_{s})\right]$

return _$\langle\phi(G),\phi(G^{\prime})\rangle$_

Algorithm 1 Weisfeiler-Lehman Subtree Graph Kernel

The sequence-to-sequence autoencoder network is constructed as in Figure 4
with three hidden layers – two in the encoder and one in the decoder. The
first hidden layer in the encoder is a fully connected layer for learning x86
instruction embedding vectors of length 64. The weights in this layer can be
pre-trained using a Word2vec [13], GloVe [16], or some other similar
unsupervised word embedding technique. However in our experiments we observed
no ill effects from initializing the weights randomly and learning the
instruction embedding directly as part of the autoencoder training.

###  3.3 Function Clustering

Having obtained an x86 instruction sequence encoder, we embed all internal
functions found in the function call graph and attribute the vertices with
their respective embeddings. External vertices representing imported functions
or APIs remain non-attributed since we are unable to obtain function
embeddings for these functions. They therefore retain their discrete external
labels. The graphs then are composed of two types of vertices: internal
functions with continuous attributes but arbitrary discrete labels and
external functions with discrete labels but no continuous attributes. In order
to carry out the graph classification task with established Weisfeiler-Lehman
kernels, we must first obtain learned discrete labels for internal functions.
We achieve this by clustering their continuous embeddings and labeling
clusters according to their cluster identifier. We can then label the internal
functions in the graph with their discrete cluster identifier. For the
purposes of graph classification, our labels can be non-descriptive (e.g.
$C_{1},C_{2},...$) but it may be beneficial in future work to consider
assigning descriptive cluster labels to aid humans in manual threat research
analysis.

Even modestly-sized files can contain tens of thousands of individual
functions so it is important to be able to scale our clustering task for many
function samples. Thus, we use the popular mini-batch K-means clustering
algorithm [19] as it is easily able to scale to many millions of samples.
Using the so-called “elbow method”, we found that $k=7000$ was a reasonable
choice for the number of clusters for our data sets. Although several
hierarchical- and density-based clustering algorithms have been found to yield
superior results over K-means, these algorithms are generally unable to scale
to many samples.

###  3.4 Graph Classification

A graph kernel is a function that computes an inner product between graphs and
can be thought of as a way to measure graph similarity. Graph kernels are
widely studied since they allow kernel-based machine learning algorithms such
as SVMs to be applied directly to graph-structured data. Most graph kernels
are based on the Weisfeiler-Lehman test of graph isomorphism. Indeed, we
employ the popular Weisfeiler-Lehman subtree kernel algorithm to obtain the
whole-graph feature vector. Once a feature vector for the graph is obtained,
we can compute a pairwise kernel matrix and train an SVM for malware
classification.

Let the Weisfeiler-Lehman kernel with base kernel $k$ be defined as

$K^{(h)}_{WL}(G,G^{\prime})=\sum_{i=0}^{h}{\alpha_{i}k(G_{i},G^{\prime}_{i})}$

where $\\{G_{0},G_{1},...,G_{h}\\}$ and
$\\{G^{\prime}_{0},G^{\prime}_{1},...,G^{\prime}_{h}\\}$ are sequences of
graphs that the Weisfeiler-Lehman algorithm generates from $G$ and
$G^{\prime}$ respectively after $h$ iterations. For each iteration of the
Weisfeiler-Lehman algorithm, each vertex obtains a new label by concatenating
the labels of the adjacent vertices. This process is illustrated for one
iteration of the algorithm in Figure 5. The Weisfeiler-Lehman kernel, then, is
simply the weighted sum of the base kernel function applied to the graphs
generated by the Weisfeiler-Lehman algorithm. Let the subtree base kernel $k$
be defined as the inner product between $\phi(G)$ and $\phi(G^{\prime})$

$k(G,G^{\prime})=\langle\phi(G),\phi(G^{\prime})\rangle$

where
$\phi(Q)=\left[c(Q,\sigma_{0}),c(Q,\sigma_{1}),...,c(Q,\sigma_{s})\right]$ and
$c(Q,\sigma_{i})$ is the count of vertex label $\sigma_{i}\in\Sigma^{(h)}$
occurring in the graph $Q$. The set of vertex labels obtained after $h$
iterations of the Weisfeiler-Lehman algorithm is denoted as $\Sigma^{(h)}$. To
compute the kernel matrix, we compute the pairwise kernels for all graphs as

$K^{(h)}=\begin{bmatrix}k^{(h)}(G_{1},G_{1})&\dots&k^{(h)}(G_{1},G_{N})\\\
k^{(h)}(G_{2},G_{1})&\dots&k^{(h)}(G_{2},G_{N})\\\ \vdots&\ddots&\vdots\\\
k^{(h)}(G_{N},G_{1})&\dots&k^{(h)}(G_{N},G_{N})\end{bmatrix}$

Shervashidze, et al (2011) have shown that for $N$ graphs with $n$ vertices
and $m$ edges, the Weisfeiler-Lehman subtree kernel of height $h$ can be
computed in $O(Nhm+N^{2}hn)$ time. This kernel matrix can be supplied directly
to a one-versus-all support vector machine with Platt scaling in order to
obtain the class probabilities [17].

##  4 Experiments

###  4.1 Setup

We performed three separate end-to-end experiments of our malware classifier
system which is composed of the individually trained components below. Each
component feeds into the next to form the final multiclass classifier system.

  * •

Sequence-to-sequence autoencoder model

  * •

K-means clustering model

  * •

Weisfeiler-Lehman subtree kernel model

Table 1: Malware descriptions Family Name | Samples | Type
---|---|---
Ramnit | 1263 | Worm
Lollipop | 2306 | Adware
Kelihos_ver3 | 2931 | Backdoor
Vundo | 344 | Trojan
Simda | 33 | Backdoor
Tracur | 663 | TrojanDownloader
Kelihos_ver1 | 382 | Backdoor
Obfuscator.ACY | 1158 | Obfuscated Malware
Gatak | 954 | Backdoor

Ten percent of the original data set was withheld from training altogether,
not seen by any of the individual sub-models and was used only for testing the
composite classifier models. Ten percent of the training set was used for
validation. Deep learning library Keras [4] was used to construct and train
the sequence-to-sequence network. The network was trained on an NVIDIA Tesla
K80 GPU. We used the GraKeL [22] implementation of the Weisfeiler-Lehman
algorithm with three iterations to obtain the Weisfeiler-Lehman graphs and the
kernel matrix. The scikit-learn [15] implementation for support vector
machines based on LIBSVM [2] was used for the graph classification task with
hyperparameters $C$ and $\gamma$ being obtained through grid search.

Table 2: Results summary Family Name | Precision | Recall | F1-Score | TP | FP | FN | TN | Support
---|---|---|---|---|---|---|---|---
Ramnit | 0.970 | 0.995 | 0.982 | 382 | 12 | 2 | 2682 | 384
Lollipop | 0.997 | 1.000 | 0.998 | 708 | 2 | 0 | 2358 | 708
Kelihos_ver3 | 1.000 | 1.000 | 1.000 | 861 | 0 | 0 | 2205 | 861
Vundo | 0.989 | 1.000 | 0.994 | 93 | 1 | 0 | 2973 | 93
Simda | 1.000 | 1.000 | 1.000 | 15 | 0 | 0 | 3051 | 15
Tracur | 0.984 | 1.000 | 0.992 | 180 | 3 | 0 | 2886 | 180
Kelihos_ver1 | 1.000 | 0.967 | 0.983 | 87 | 0 | 3 | 2979 | 90
Obfuscator.ACY | 1.000 | 0.974 | 0.987 | 368 | 0 | 10 | 2698 | 378
Gatak | 1.000 | 0.992 | 0.996 | 354 | 0 | 3 | 2712 | 357

###  4.2 Data Set

We use the Microsoft Malware Classification data set [18] to evaluate our
approach. The data set consists of samples from nine different malware
families. Each sample in the data set is composed of a pre-disassembled .asm
file generated by IDA and a sanitized hexadecimal representation of the
original executable. Our approach makes use of only the .asm and furthermore
only takes advantage of the parsable code sections. Having obtained only the
publicly available training set of 10,867 samples, we were able to extract
function call graphs for 10,152 samples due, in part, to certain samples being
packed or otherwise obfuscated.

###  4.3 Results

We achieved a prediction accuracy of 99.41% in the malware family
classification task. Our approach outperforms other malware classifiers that
involve extensive feature engineering or extract significantly more data from
the executable such as non-code data [8, 9, 20]. Since we only use the code
sections of the executable, we expect that incorporating additional data such
as the .rsrc and .idata sections would help to further improve classification
results. Table 2 summarizes the results across the three experiments.

##  5 Conclusion

In this work we applied several machine learning techniques to the problem of
malware detection and achieved over 99% accuracy in the malware classification
task using a composite model. A sequence-to-sequence autoencoder was used to
obtain dense, latent representations of x86 code which helped our model
account for anti-malware evasion practices. We then clustered the function
representations of the functions and obtained discrete function labels. Using
the discrete labels, we constructed a function call graph where vertices
represent functions and are labeled according to their cluster IDs. The
Weisfeiler-Lehman graph kernel framework was used to obtain the Weisfeiler-
Lehman graphs and to construct a kernel matrix which allowed us to ultimately
perform the graph classification task using a support vector machine.

##  6 Acknowledgments

Special thanks to our colleague, Andrew Sandoval, for his valuable
contributions and feedback throughout this work.

Figure 6: Normalized confusion matrix

## References

  * [1] Airola, A., Pyysalo, S., Björne, J., Pahikkala, T., Ginter, F., Salakoski, T.: A graph kernel for protein-protein interaction extraction. In: Proceedings of the Workshop on Current Trends in Biomedical Natural Language Processing. pp. 1–9. BioNLP ’08, Association for Computational Linguistics, Stroudsburg, PA, USA (2008), http://dl.acm.org/citation.cfm?id=1572306.1572308
  * [2] Chang, C.C., Lin, C.J.: LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology 2, 27:1–27:27 (2011), software available at http://www.csie.ntu.edu.tw/~cjlin/libsvm
  * [3] Cho, K., van Merrienboer, B., Gulcehre, C., Bahdanau, D., Bougares, F., Schwenk, H., Bengio, Y.: Learning Phrase Representations using RNN Encoder-Decoder for Statistical Machine Translation. arXiv e-prints arXiv:1406.1078 (Jun 2014)
  * [4] Chollet, F., et al.: Keras. https://keras.io (2015)
  * [5] Dam, K.H.T., Touili, T.: Malware detection based on graph classification. In: ICISSP (2017)
  * [6] Dullien, T.: Graph-based comparison of executable objects (2005)
  * [7] Goodfellow, I., Bengio, Y., Courville, A.: Deep Learning. The MIT Press (2016)
  * [8] Hassen, M., Carvalho, M.M., Chan, P.K.: Malware classification using static analysis based features. In: 2017 IEEE Symposium Series on Computational Intelligence (SSCI). pp. 1–7 (Nov 2017). https://doi.org/10.1109/SSCI.2017.8285426
  * [9] Hassen, M., Chan, P.: Scalable function call graph-based malware classification. In: Proceedings of the Seventh ACM on Conference on Data and Application Security and Privacy, CODASPY ’17. pp. 239–248 (Mar 2017). https://doi.org/10.1145/3029806.3029824
  * [10] Hex-Rays: The ida pro disassembler and debugger. https://www.hex-rays.com/products/ida/ (Jan 2011)
  * [11] Liben-nowell, D., Kleinberg, J.: The link prediction problem for social networks. Journal of the American Society for Information Science and Technology 58 (Jan 2003). https://doi.org/10.1002/asi.20591
  * [12] Mesbahi, M., Egerstedt, M.: Graph theoretic methods in multiagent networks. Princeton University Press (2010)
  * [13] Mikolov, T., Chen, K., Corrado, G., Dean, J.: Efficient estimation of word representations in vector space. arXiv preprint arXiv:1301.3781 (2013)
  * [14] Murphy, G.C., Notkin, D., Griswold, W.G., Lan, E.S.: An empirical study of static call graph extractors. ACM Trans. Softw. Eng. Methodol. 7(2), 158–191 (Apr 1998). https://doi.org/10.1145/279310.279314, http://doi.acm.org/10.1145/279310.279314
  * [15] Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., Vanderplas, J., Passos, A., Cournapeau, D., Brucher, M., Perrot, M., Duchesnay, E.: Scikit-learn: Machine learning in Python. Journal of Machine Learning Research 12, 2825–2830 (2011)
  * [16] Pennington, J., Socher, R., Manning, C.: Glove: Global vectors for word representation. In: Proceedings of the 2014 conference on empirical methods in natural language processing (EMNLP). pp. 1532–1543 (2014)
  * [17] Platt, J.C.: Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. In: Advances In Large Margin Classifiers. pp. 61–74. MIT Press (1999)
  * [18] Ronen, R., Radu, M., Feuerstein, C., Yom-Tov, E., Ahmadi, M.: Microsoft Malware Classification Challenge. arXiv e-prints arXiv:1802.10135 (Feb 2018)
  * [19] Sculley, D.: Web-scale k-means clustering. In: Proceedings of the 19th International Conference on World Wide Web. pp. 1177–1178. WWW ’10, ACM, New York, NY, USA (2010). https://doi.org/10.1145/1772690.1772862, http://doi.acm.org/10.1145/1772690.1772862
  * [20] Searles, R., Xu, L., Killian, W., Vanderbruggen, T., Forren, T., Howe, J., Pearson, Z., Shannon, C., Simmons, J., Cavazos, J.: Parallelization of machine learning applied to call graphs of binaries for malware detection. In: 25th Euromicro International Conference on Parallel, Distributed and Network-based Processing (PDP). pp. 69–77. IEEE (2017)
  * [21] Shervashidze, N., Schweitzer, P., Leeuwen, E.J.v., Mehlhorn, K., Borgwardt, K.M.: Weisfeiler-lehman graph kernels. Journal of Machine Learning Research 12(Sep), 2539–2561 (2011)
  * [22] Siglidis, G., Nikolentzos, G., Limnios, S., Giatsidis, C., Skianis, K., Vazirgiannis, M.: Grakel: A graph kernel library in python. arXiv preprint arXiv:1806.02193 (2018)
  * [23] Sutskever, I., Vinyals, O., Le, Q.V.: Sequence to sequence learning with neural networks. In: Advances in neural information processing systems. pp. 3104–3112 (2014)