# Terahertz Slonczewski propagating spin waves  and large output voltage in
antiferromagnetic spin-Hall nano-oscillators

M. Hamdi<EMAIL_ADDRESS>
<EMAIL_ADDRESS>École Polytechnique Fédérale de Lausanne (EPFL),
Institute of Materials, Laboratory of Nanoscale Magnetic Materials and
Magnonics, CH-1015 Lausanne, Switzerland  D. Grundler<EMAIL_ADDRESS>École Polytechnique Fédérale de Lausanne (EPFL), Institute of Materials,
Laboratory of Nanoscale Magnetic Materials and Magnonics, CH-1015 Lausanne,
Switzerland  École Polytechnique Fédérale de Lausanne (EPFL), Institute of
Electrical and Micro engineering, CH-1015 Lausanne, Switzerland

###### Abstract

We study theoretically antiferromagnet (AFM) based spin-Hall nano-oscillators
(SHNOs) consisting of a nano-constriction (NC) in a thin-film uniaxial AFM. By
solving the derived SW equation we evidence radially propagating spin waves
(SWs) at THz frequencies similar to the Slonczewski SWs known at GHz
frequencies for a ferromagnet-based SHNO. We predict a minimum threshold
current for a specific NC radius accessible by the state-of-the-art
nanotechnology. The exchange interaction enhanced spin pumping for AFMs leads
to a strong thickness dependent threshold frequency. We show that the uniaxial
AFMs generate ac electrical fields via spin pumping that are three orders of
magnitude larger than reported for biaxial AFMs. Our work enhances the
fundamental understanding of current-driven SWs in AFM-SHNOs and enables
optimization of practical devices in terms of material choice, device
geometry, and frequency tunability. The propagating SWs offer remote THz
signal generation and an efficient means for synchronization of SHNOs when
aiming at high power.

Introduction.— Terahertz (THz) radiation attracts tremendous attention due to
its several promising applications including ultrafast communication [1, 2],
biosensing and imaging [3, 4]. However, the widespread use of THz
instrumentation is constrained, in particular, because of the lack of
efficient sources and detectors of THz radiation for frequencies in the range
of 0.1-10 THz (THz gap) [3, 4]. Therefore, there is an urgent need to develop
small-sized solid-state THz sources and detectors. Spin wave (SW) excitations
in antiferromagnetic materials (AFMs) offer frequencies which cover perfectly
the THz gap. These materials can be utilized to close the THz gap in terms of
antiferromagnetic spintronic and magnonic devices [5, 6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30].
One of the most promising devices is the AFM-based spin-Hall nano-oscillator
(AFM-SHNO) consisting of an AFM and a heavy metal (HM) layer [14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. In an AFM-SHNO, a spin
current generated by the spin-Hall effect (SHE) in the HM enters the AFM layer
and exerts the so-called spin-orbit torque (SOT) [31, 32, 33, 34]. It has been
shown theoretically that when the spin current exceeds a certain threshold,
damping can  be compensated, and the magnetization starts self-sustained auto-
oscillations. They can generate THz radiation via the combination of spin
pumping [5, 17, 35, 36] and inverse spin-Hall effect [14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. The AFM-SHNOs have potential not
only as THz sources and detectors [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27, 28], but also as ultra-fast artificial neurons which generate
picosecond-duration spikes [29, 30] and enable ultra-fast neuromorphic
applications [37, 38, 39, 40]. Howver, there has not yet been any report on
the experimental realization of an AFM-SHNOs.
All the theoretical studies mentioned above were carried out by considering an
AFM in the macrospin approximation. They described uniform antiferromagnetic
resonance (AFMR) excitation by a spin current in planar bilayer HM/AFM
structures [14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30]. Still, there have been few micromagnetic studies of SHNOs addressing the
planar HM/AFM geometry and AFMR excitation [22, 23]. These studies do not
describe however a realistic AFM-SHNO which incorporates a nano-constriction
(NC) [41, 42, 43, 44, 45].  NCs are essential for ferromagnet-based SHNOs (FM-
SHNOs) which excite propagating SWs [33, 46, 47]. These SWs enable the
efficient synchronization of FM-SHNOs [48, 49, 50, 51, 52, 53] granting
sufficient signal quality for applications.
In this Letter, we study theoretically and by means of micromagnetic
simulations an NC-based AFM-SHNO and provide a fundamental understanding of
its spin dynamics. We first derive the SW equation for an insulating AFM layer
subjected to a spin current in a circular NC. We show that the solution to
this equation is the antiferromagnetic counterpart of Slonczewski SWs known
for ferromagnetic SHNOs. The Slonczewski spin waves are radially propagating
waves generated by spin torque in a circular NC in a FM layer [33, 34]. We
obtain the threshold frequency and current for such THz antiferromagnetic
Slonczewski SWs (THz-ASSWs) as a function of the material parameters and the
NC radius $R$. Our analytical results show good agreement with the
micromagnetic simulations. A unique feature of these THz-ASSWs is that the
threshold current exhibits a minimum with respect to the NC radius which
allows one to reduce the power consumption to a minimum. Interestingly, the
minimum diameter $2R$ is in the range of few tens to few hundreds of
nanometers for realistic material parameters and is realizable by state-of-
the-art fabrication technology.  Moreover, we show that when used in NC-SHNOs,
the uniaxial AFMs generate THz outputs, which are three orders of magnitude
larger compared to biaxial AFM based uniformly excited SHNOs. This is a
significant improvement on both the power consumption and output signal of the
AFM-SHNOs and draws the attention to the previously abandoned uniaxial AFMs
which were considered not suitable for THz SHNOs [14, 15].
Theory.— We consider an AFM-SHNO as a plane AFM film that is subjected to a
spin current injected locally in a circular region of radius $R$, as depicted
in Fig. 1 (a). Such spin current is generated by restricting the electric
current in the adjacent HM into a nano-constriction [41, 42]. We assume an AFM
with a uniaxial anisotropy in the plane along the $x$-axis and a thickness,
$d_{AF}$, smaller than the exchange length,
$\lambda_{ex}=\sqrt{a^{2}\omega_{E}/4\omega_{A}}$. Thereby we avoid standing
waves along the thickness. Here, $a$ is the magnetic lattice constant of the
material and $\omega_{E}=\gamma H_{E}$ ($\omega_{A}=\gamma H_{A}$) is the
exchange (anisotropy) field in the units of frequency with $H_{E}$ ($H_{A}$)
being the exchange (anisotropy) field and $\gamma$ is the electron
gyromagnetic ratio. The Néel order vector and average magnetization are
defined as $\mathbf{n}=(\mathbf{M}_{1}-\mathbf{M}_{2})/2M_{s}$ and
$\mathbf{m}=(\mathbf{M}_{1}+\mathbf{M}_{2})/2M_{s}$, respectively, where
$\mathbf{n}.\mathbf{m}=0$ and $|\mathbf{n}|^{2}+|\mathbf{m}|^{2}=1$.
$\mathbf{M}_{i}$ ($i=1,2$) is the sublattice magnetization with saturation
value of $M_{s}=|\mathbf{M}_{i}|$. We consider the polarization direction of
the spin current, $\mathbf{\hat{p}}$, to be parallel to $\mathbf{n}$ where the
SOT has the highest efficiency [23]. There is no externally applied magnetic
field. The total magnetic free energy of the AFM is $F=\int fdV$, where the
energy density, $f$, is given by [54, 5, 55, 17, 56]

$\displaystyle f=\dfrac{M_{s}}{\gamma}$
$\displaystyle\left\\{\vphantom{\int}\right.\dfrac{1}{2}\omega_{E}\left(\mathbf{m}^{2}-\mathbf{n}^{2}\right)-\dfrac{1}{2}\omega_{A}\left(m_{x}^{2}+n_{x}^{2}\right)$
(1)
$\displaystyle-\dfrac{1}{2}\omega_{A}\lambda_{ex}^{2}\left[\left(\mathbf{\nabla}\mathbf{m}\right)^{2}-\left(\mathbf{\nabla}\mathbf{n}\right)^{2}\right]\left.\vphantom{\int}\right\\}.$

The dynamics of $\mathbf{m}$ and $\mathbf{n}$ is governed by the Landau-
Lifshitz-Gilbert-Slonczewski (LLGS) equations with the spin-torque term as

$\displaystyle\mathbf{\dot{m}}$ $\displaystyle=$
$\displaystyle\left[\vphantom{\int}\right.\bm{\omega}_{m}\times\mathbf{m}+\bm{\omega}_{n}\times\mathbf{n}\left.\vphantom{\int}\right]+\bm{\tau}_{m}^{\text{D}}+\bm{\tau}_{m}^{\text{SOT}},$
(2a) $\displaystyle\mathbf{\dot{n}}$ $\displaystyle=$
$\displaystyle\left[\vphantom{\int}\right.\bm{\omega}_{m}\times\mathbf{n}+\bm{\omega}_{n}\times\mathbf{m}\left.\vphantom{\int}\right]+\bm{\tau}_{n}^{\text{D}}+\bm{\tau}_{n}^{\text{SOT}},$
(2b)

where $\bm{\omega}_{m}=-\frac{\gamma}{M_{s}}\frac{\delta F}{\delta\mathbf{m}}$
and $\bm{\omega}_{n}=-\frac{\gamma}{M_{s}}\frac{\delta F}{\delta\mathbf{n}}$
are the effective fields in frequency units [54, 5, 17, 56]. The damping
torques are given by
$\bm{\tau}_{m}^{\text{D}}=\alpha\left(\vphantom{\int}\right.\mathbf{m}\times\mathbf{\dot{m}}+\mathbf{n}\times\mathbf{\dot{n}}\left.\vphantom{\int}\right)$
and
$\bm{\tau}_{n}^{\text{D}}=\alpha\left(\vphantom{\int}\right.\mathbf{m}\times\mathbf{\dot{n}}+\mathbf{n}\times\mathbf{\dot{m}}\left.\vphantom{\int}\right)$
where $\alpha=\alpha_{0}+\alpha_{\text{SP}}$ is the effective damping
consisting of the Gilbert damping, $\alpha_{0}$, and the enhanced damping due
to spin pumping [17, 5], $\alpha_{\text{SP}}=\frac{\gamma\hbar g_{r}}{4\pi
M_{s}d_{AF}}$. Here, $\hbar$ and $g_{r}$ are the reduced Planck constant and
interface spin mixing conductance, respectively. The spin-torque terms are
given as
$\bm{\tau}_{m}^{\text{SOT}}=\omega_{s}\left[\vphantom{\int}\right.\mathbf{m}\times\left(\mathbf{m}\times\mathbf{\hat{p}}\right)+\mathbf{n}\times\left(\mathbf{n}\times\mathbf{\hat{p}}\right)\left.\vphantom{\int}\right]H\left(R-r\right)$
and
$\bm{\tau}_{n}^{\text{SOT}}=\omega_{s}\left[\vphantom{\int}\right.\mathbf{m}\times\left(\mathbf{n}\times\mathbf{\hat{p}}\right)+\mathbf{n}\times\left(\mathbf{m}\times\mathbf{\hat{p}}\right)\left.\vphantom{\int}\right]H\left(R-r\right)$
where $\omega_{s}=\frac{2\hbar\gamma\theta_{SH}J_{c}\eta}{eM_{s}d_{AF}}$ with
$\eta=\frac{\lambda_{N}G_{r}\tanh\left(d_{N}/2\lambda_{N}\right)}{\sigma_{N}+2\lambda_{N}G_{r}\coth\left(d_{N}/\lambda_{N}\right)}$is
the SOT strength in frequency units, $H(R-r)$ is the Heaviside step function
and $r$ is the radial distance from the center of the NC. $\theta_{SH}$,
$\lambda_{N}$, $\sigma_{N}$, $d_{N}$ and $J_{c}$ are spin-Hall angle, spin
diffusion length, electrical conductivity, thickness and applied electric
current of the HM, respectively, and $G_{r}=e^{2}g_{r}/2\pi\hbar$ where $-e$
is the electron charge. We consider the ansatz
$\mathbf{n}=n_{0}\mathbf{\hat{x}}+\mathbf{n}_{\bot}\left(x,y\right)e^{i\omega
t}$ with $|\mathbf{n}_{\bot}\left(x,y\right)|\ll 1$ to derive the wave
equation for small amplitude SWs excited by SOT. In the exchange limit
($\omega_{A}/\omega_{E}\ll 1$) which is the case for most antiferromagnets and
the absence of an external field,
$\mathbf{m}\simeq(\omega_{A}/\omega_{E})\mathbf{n}_{\bot}\ll\mathbf{n}_{\bot}$
and $\left|n_{0}\right|\simeq 1$. By keeping the linear terms in $\mathbf{m}$
and $\mathbf{n}_{\bot}$ and eliminating $\mathbf{m}$ between the Eqs. 2a and
2b [57], we obtain two coupled equations for two components of
$\mathbf{n}_{\bot}$ as

$\displaystyle\left(\dfrac{\mp\omega^{2}}{2\omega_{E}+\omega_{A}}\right)n_{z(y)}=\mp\omega_{A}n_{z(y)}\mp
i\alpha\omega n_{z(y)}$
$\displaystyle+\omega_{s}n_{y(z)}\pm\omega_{A}\lambda_{ex}^{2}\nabla^{2}n_{z(y)}.$
(3)

By defining the complex SW amplitude as $u=n_{y}+in_{z}$, and the normalized
distance $\tilde{r}=r/R$ and by adding the dispersive dissipation term ,
$i\alpha\omega_{A}\lambda_{ex}^{2}\nabla^{2}u$,  which describes the increase
of the relaxation rate with the wave vector $k$ of the excited SW [58, 59], we
obtain the SW equation as [57]

$\displaystyle\dfrac{1}{\tilde{r}}\dfrac{\partial}{\partial\tilde{r}}\left(\tilde{r}\dfrac{\partial
u}{\partial\tilde{r}}\right)+Au+iBu+\alpha Bu-i\alpha Cu=0,$ (4a)
$\displaystyle
A=\left(\dfrac{R^{2}}{\omega_{A}\lambda_{ex}^{2}}\right)\left(\dfrac{\omega^{2}}{2\omega_{E}+\omega_{A}}-\omega_{A}\right),$
(4b) $\displaystyle
B=\left(\dfrac{R^{2}}{\omega_{A}\lambda_{ex}^{2}}\right)\omega_{s}H\left(1-\tilde{r}\right),$
(4c) $\displaystyle
C=\left(\dfrac{R^{2}}{\omega_{A}\lambda_{ex}^{2}}\right)\left[\omega+\dfrac{\omega^{2}}{2\omega_{E}+\omega_{A}}-\omega_{A}\right].$
(4d)

One can summarize Eqs. 4a-4d as
$\frac{1}{\tilde{r}}\frac{\partial}{\partial\tilde{r}}\left(\tilde{r}\frac{\partial
u}{\partial\tilde{r}}\right)+k_{i}^{2}u=0$, which is the Bessel equation with
complex wave vector $k_{1}=\sqrt{A+iB+\alpha\left(B-iC\right)}$ for
$0\leq\tilde{r}\leq 1$ and $k_{2}=\sqrt{A-i\alpha C}$ for $\tilde{r}>1$.
Following Slonczewski [60, 46, 61], we take a piecewise Bessel solution of the
form

$\displaystyle
u\left(\tilde{r}\right)=\begin{cases}J_{0}\left(k_{1}\tilde{r}\right);\quad
0\leq r\leq 1,\\\ cH_{0}^{\left(1\right)}\left(k_{2}\tilde{r}\right);\quad
r>1,\text{ with
}c=\dfrac{J_{0}\left(k_{1}\tilde{r}\right)}{H_{0}^{\left(1\right)}\left(k_{2}\tilde{r}\right)},\end{cases}$
(5)

where $J_{0}(x)$ and $H_{0}^{\left(1\right)}(x)$ are the Bessel and Hankel
functions of the first kind, respectively. The continuity of the derivative of
$u$ at the NC boundary $\tilde{r}=1$, leads to a complex transcendental
equation whose real and imaginary parts determine two real eigenvalues
$\omega$ and $\omega_{s}$ at the threshold [57, 60, 61]. By solving this
equation we have derived the frequency of the excited mode and required
current density at the threshold as

$\displaystyle\omega_{\text{th}}=\dfrac{\alpha
A_{12}\left(2\omega_{E}+\omega_{A}\right)}{2\left(1-\alpha
A_{12}\right)}+\sqrt{\left[\dfrac{\alpha
A_{12}\left(2\omega_{E}+\omega_{A}\right)}{2\left(1-\alpha
A_{12}\right)}\right]^{2}+\left(2\omega_{E}+\omega_{A}\right)\left[\omega_{A}+\left(\dfrac{A_{0}+\alpha
A_{11}}{1-\alpha
A_{12}}\right)\left(\dfrac{\omega_{A}\lambda_{ex}^{2}}{R^{2}}\right)\right]},$
(6) $\displaystyle
J_{c,\text{th}}=\left(\dfrac{eM_{s}d_{AF}}{2\gamma\hbar\theta_{SH}\eta}\right)\left\\{\left(\dfrac{\alpha
B_{12}}{1-\alpha
A_{12}}\right)\omega_{\text{th}}+\left(\dfrac{\omega_{A}\lambda_{ex}^{2}}{R^{2}}\right)\left[B_{0}+\alpha
B_{11}+\alpha B_{12}\left(\dfrac{A_{0}+\alpha A_{11}}{1-\alpha
A_{12}}\right)\right]\right\\},$ (7)

respectively, where $A_{0}=1.428$, $A_{11}=-3.087$, $A_{12}=0.815$,
$B_{0}=1.863$, $B_{11}=-1.581$ and $B_{12}=1.004$.
Micromagnetic simulations.— We use the Object Oriented MicroMagnetic Framework
(OOMMF) [62] which solves the LLGS equation based on the finite difference
method to validate the analytical expressions derived from the continuum model
(Eqs. 6 and 7). These simulations are performed on a cuboid-like mesa with
dimensions of 2 $\mu$m×2 $\mu$m×1 nm, divided into 2000×2000×1 cells. A
circular NC with a radius of $R$ is defined at the center of the device for
spin current driven excitations (Fig. 1(a)). In order to compare the
simulations with our theoretical results, we performed the simulations for
different NC radii of $R=10$, 25, 75 and 125 nm. In all the simulations, we
used $M_{s}=350$ kA/m, $K_{x}=20$ kJ/m3, $\alpha=0.001$ as the sublattice
saturation magnetization, magnetic anisotropy and Gilbert damping,
respectively. In order to achieve antiferromagnetic order, we used
$A_{ex}=-10$ pJ/m as the exchange stiffness in OOMMF simulations. This
corresponds to a two-sublattice G-type antiferromagnet with homogeneous
(nonhomogeneous) exchange constant of $A_{h}=-4A_{ex}/a^{2}$
($A_{nh}=-A_{ex}$) with the magnetic lattice constant $a=1$ nm (See
Supplemental Material [57]). These parameters correspond to
$\mu_{0}H_{E}=228.57$ T and $\mu_{0}H_{A}=0.114$ T  which are in the range of
material parameters for the existing AFMs. For HM material parameters we
considered $d_{HM}=10$ nm, $\theta_{SH}=0.15$, $\lambda_{HM}=7$ nm and
$\sigma_{N}=2.083\times 10^{6}$ $\Omega^{-1}$m-1. Considering $g_{r}=6.9\times
10^{18}$ m-2 [5] results in $\eta=0.18$ for the interface. Absorbing
boundaries were used to prevent back-reflection and interference effects [63].
The time dependent Néel order vector of the AFM layer under the d.c. current
was sampled every 0.1 ps for a total number of 1000 samplings for each
simulation.
Results and discussion.— By taking the fast-Fourier-transform (FFT) of
$\mathbf{n}$ we obtain the power spectral density (PSD) of the magnetization
dynamics.

Figure 1: (a) Schematics of the studied system. The spin current with
polarization $\mathbf{p}$ originates from the SHE of the electric current
constrained in a circular NC of radius $R$ in the adjacent HM layer. The spin
polarization of the spin current is assumed to be along the Néel order vector.
(b) Excitation frequency versus applied current density for the simulated
SHNOs with different NC radii $R$.  The black upward (downward) arrow
indicates the mode for $R=25$ nm at threshold $J_{c,\text{th}}=2.12\times
10^{8}$ A/cm2 with $f_{\text{th}}=0.3$ THz ($J_{c,\text{h}}=11.3\times 10^{8}$
A/cm2 with $f_{\text{h}}=1.2$ THz) which is selected to calculate output
electric fields due to spin pumping in Fig. 4. The slope $\beta$ as a function
of $R^{2}$ is shown in the inset. Threshold current density (c) and frequency
(d) obtained from simulations (filled red circles) and Eqs. 7 and 6 (blue
line), respectively. Gray dashed lines indicate the plane film limit. Inset in
(d) shows a snapshot of current induced SW at threshold for SHNO with $R=50$
nm (marked with red circle) and simulation parameters given in the text.

Varying the current density and extracting the peak position of the PSD, we
obtained the frequency versus current density curve of the SHNOs plotted in
Fig. 1(b).  The frequency variation for large currents is attributed to the
excessive SOT beyond threshold according to
$\omega=\omega_{\text{th}}+\beta(J_{c}-J_{\text{th}})$. The inset in Fig. 1(b)
shows the SW frequency shift coefficient $\beta=d\omega/dJ_{c}$ as a function
of the NC area. $\beta$ as determined from micromagnetic simulations increases
with $R$ and seems to asymptotically approach the macrospin limit. The
observed tunability of $\beta$ goes beyond macrospin studies where $\beta$ is
a constant and can be used to engineer NC performance.
Now we focus on the dynamics at the threshold. At the threshold, the PSD peak
jumps from zero to the threshold frequency. Here, the dynamics of the SHNO is
in the linear regime. The threshold current density and corresponding
frequency decrease with the NC radius are plotted in Fig. 1(c) and (d),
respectively. The values $J_{c,\text{th}}$ obtained from micromagnetic
simulations (red symbols) are plotted with an offset of $\Delta J=-0.4\times
10^{8}$ A/cm2. We attribute this systematic error to the Taylor expansion of
the complex transcendental equation that we employed to solve the eigenvalue
problem [57]. The fixed offset does not affect the main conclusion of the
comparison, i.e., the significant dependency of $J_{c,\text{th}}$ on $R$
making a geometrical optimization of an NC device necessary. The inset in Fig.
1(c) shows a zoom in of the region indicated by dashed rectangle. By
considering the constant offset, the remaining discrepancies amount to only
-1$\%$ and 6$\%$, which can be due to the coarse sweeping of the applied
current in the simulations.
For large $R\gg\lambda_{ex}$ the threshold current density $J_{c,\text{th}}$
(frequency $\omega_{\rm th}$) asymptotically approaches the AFMR limit in a
thin film,
$J_{c,\text{th}}=\left(\frac{eM_{s}d_{AF}}{2\gamma\hbar\theta_{SH}\eta}\right)\alpha\omega_{0}$
($\omega_{0}=\sqrt{\left(2\omega_{E}+\omega_{A}\right)\omega_{A}}$), indicated
by the gray dashed line in Fig. 1(c) (1(d)). The inset in Fig. 1(d) shows a
snapshot of the spin-current-generated spin wave simulated for $R=50~{}$nm.
For small $R$ (below almost $\lambda_{ex}=22.3$ nm), the exchange term
dominates and the threshold current density and frequency are large. One can
understand the increase in threshold current density with decreasing $R$ by
considering the energy flow away from the NC by outgoing SWs. The wave vector,
$k$, of the excited spin wave at the threshold is comparable to $2\pi/R$.
Therefore, by decreasing $R$ one increases $k$ and the group velocity of the
excited wave (following the dispersion in Eq. 6) and hence, the propagating
energy away from the NC. To sustain this energy flow, a higher threshold
current is necessary.

Figure 2: (a) Threshold current, $I_{\text{th}}$, as a function of $R$ for the
parameters used in Fig. 1 (blue line). The orange line considers Cr2O3
material parameters with the same HM and interface parameters as Fig. 1. NC
radius (b), threshold current (c) and frequency (d) at the minimum point of
$I_{\text{th}}(R)$ as a function of exchange and anisotropy fields. For all
graphs we considered $d_{AF}=20$ nm as well as the corresponding
$\alpha_{\text{SP}}$.  Red asterisks indicate the material parameters used in
Fig. 1. White circles highlight further AFMs.

For practical applications the power consumption of the SHNO is of crucial
importance. Approximating the current passing part of the NC in HM as a
resistance with the length $L=2R$ and cross-sectional area $S=2Rd_{N}$ in the
$xz-$plane one obtains its resistance as
$\sim(L/\sigma_{N}S)=1/\sigma_{N}d_{N}$ which is independent of $R$. By
multiplying Eq. 7 by $S$ we obtain the threshold current as

$I_{\text{th}}(R)=\alpha\mathcal{A}\omega_{\text{th}}R+(\mathcal{B}\lambda_{ex}^{2}/R).$
(8)

The first term in this equation represents the current needed to overcome the
damping which grows linearly with $R$. The second term reflects the
propagation energy which scales with $R^{-1}$. The combination of these two
terms leads to a minimum threshold current, $I_{\text{th}}^{\min}$, with
respect to $R$. In Fig. 2(a) we depict $I_{\text{th}}$ as a function of $R$
for the parameters used in Fig. 1 (blue line). For comparison, we display
$I_{\text{th}}(R)$ for the uniaxial antiferromagnet Cr2O3 (orange line) with
material parameters of $M_{s}=45.5$ kA/m, $H_{E}=490$ T, $H_{A}=0.035$ T and
$a=1$ nm [64, 65]. Since $\omega_{\text{th}}$ is a slowly varying function of
$R$ for $R\gtrsim 2\lambda_{ex}$ (Fig. 1(d)) away from the exchange dominated
regime, we neglect its dependency on $R$, for large $R$ and obtain
$R_{\text{min}}\simeq\sqrt{\frac{\mathcal{B}\lambda_{ex}^{2}}{\alpha\mathcal{A}\omega_{\text{th}}^{\text{min}}}}$
where $R_{\text{min}}$ and $\omega_{\text{th}}^{\text{min}}$ are the NC radius
and threshold frequency at the threshold current minimum, respectively. This
expression slightly underestimates $R_{\text{min}}$. In Fig. 2(b), (c) and (d)
we depict $R_{\text{min}}$, $I_{\text{th}}^{\min}$ and
$\omega_{\text{th}}^{\text{min}}$, respectively, as a function of $H_{E}$ and
$H_{A}$ (for a set of parameters that hold the exchange limit) which are
numerically obtained by minimizing $I_{\text{th}}(R)$. Figure 2 (b) shows that
$R_{\text{min}}$ increases (decreases) by increasing $H_{E}$ ($H_{A}$).
However, $I_{\text{th}}^{\min}$ (Fig. 2 (c)) and
$\omega_{\text{th}}^{\text{min}}$ (Fig. 2 (d)) increase by increasing both
$H_{E}$ and $H_{A}$.  The parameters used for simulations is indicated with
the red asterisk. White circles indicate some of the most interesting
materials for AFM spintronics [64, 65, 66, 67, 68, 35, 36, 69, 70, 71, 72, 73,
74] which all of them fall in the exchange limit and our theory holds.
Due to the large exchange field in AFMs, the effect of damping on
$I_{\text{th}}(R)$ is enhanced (Eq. 7 and 8). Therefore, we evaluate the
effect of AFM thickness which drastically influences $\alpha_{\text{SP}}$. The
intrinsic Gilbert damping, $\alpha_{0}$, is in the order of $10^{-6}-10^{-4}$
for insulating AFMs [69, 70, 71, 72, 73] which is negligibly small compared to
$\alpha_{\text{SP}}$ for thin films.

Figure 3: Spin pumping induced damping, $\alpha_{\text{SP}}$ (a) threshold
current $I_{\text{th}}^{\min}$ (b), NC size $R_{\text{min}}$ (c) and threshold
frequency, $\omega_{\text{th}}^{\text{min}}$ (d) at the threshold current
minimum as a function of AFM layer thickness, $d_{AF}$.

Figure 3(a) shows $\alpha_{\text{SP}}$ for the AFM simulated in Fig. 1 and
Cr2O3. Since $M_{s}$ for Cr2O3 is small, $\alpha_{\text{SP}}$ depends
prominently on $d_{AF}$. Despite the large $\alpha_{\text{SP}}$ for Cr2O3, its
$I_{\text{th}}^{\min}$ (orange curve) is considerably smaller compared to the
AFM with large $M_{s}$ (blue curve) as is shown in Fig. 3(b). Though
$\alpha_{\text{SP}}$ decreases by increasing $d_{AF}$, $I_{\text{th}}^{\min}$
increases since SOT is an interface effect and scales with $d_{AF}^{-1}$.
$R_{\text{min}}$ and $\omega_{\text{th}}^{\text{min}}$ as a function of
$d_{AF}$ are shown in Fig. 3(c) and (d), respectively. $R_{\text{min}}$
increases as a function of $d_{AF}$ (Fig. 3(c)) since
$R_{\text{min}}\propto\sqrt{\alpha_{\text{SP}}^{-1}}\propto\sqrt{d_{AF}}$ and
$\mathcal{A}$ and $\mathcal{B}$ are slowly varying functions of $\alpha$.
Despite having larger $H_{E}$ and smaller $H_{A}$, Cr2O3 has smaller
$R_{\text{min}}$ due to larger $\alpha_{\text{SP}}$ and
$\omega_{\text{th}}^{\text{min}}$. Though the AFMR frequency of Cr2O3 (164
GHz) is lower than the AFMR of Fig. 1 (202.5 GHz), the threshold frequency at
$I_{\text{th}}^{\min}$ is higher for Cr2O3 (Fig. 3(d)). The higher frequency
is a consequence of low $M_{s}$ which results in a significant change of
$\alpha_{\text{SP}}$ with $d_{AF}$.
In the following we make use of the predicted spin dynamics and evaluate
electrical signals that arise due to spin pumping into an adjacent HM layer.
The dynamic part of the Néel vector, $\mathbf{n}(r,t)$, gives rise to spin
pumping which scales with
$\mathbf{j}_{s}=\left(\mathbf{n}\times\mathbf{\dot{n}}\right)\times\hat{\mathbf{z}}$
[5, 14, 15]. By adding a HM layer either above or below the AFM, an ac
electrical field is generated according to
$\mathbf{E}=\mathcal{E}\mathbf{j}_{s}=\mathcal{E}\left(\mathbf{n}\times\dot{\mathbf{n}}\right)\times\mathbf{\hat{z}}$
with $\mathcal{E}=\hbar\eta\theta_{SH}/ed_{N}$ [14]. Taking the complex SW
amplitude, $u$, we obtain the ac signal [57]

$\displaystyle E_{\rm{ac}}(r,t)$ $\displaystyle=$
$\displaystyle\mathcal{E}n_{x}\dot{n}_{z}=\mathcal{E}\omega\sqrt{1-|u(r)|^{2}}\,\text{Re}\left\\{u(r)e^{i\omega
t}\right\\}.$ (9)

We numerically calculated
$\mathbf{j}_{s}=\left(\mathbf{n}\times\mathbf{\dot{n}}\right)\times\hat{\mathbf{z}}$
for each cell at each moment in time from micromagnetic simulations. We took
the FFT of this quantity, multiplied by $\mathcal{E}=0.7\times 10^{-7}$ V.s/cm
and color-coded the absolute values as shown for two different currents
$J_{c}$ applied to an NC with $R=25$ nm in Fig. 4(a) and (b). Thereby, we
evaluate the spin pumping at $f_{\text{th}}=0.3$ THz and $f_{\text{h}}=1.2$
THz, respectively. We display its spatially inhomogeneous characteristics in
units of V/cm. The panels of Fig. 4(a) and (b) hence show ac electric fields
in the typical HM assumed above. They allow for a straightforward comparison
with the earlier reports on spin pumping from biaxial AFMs. For the
$f_{\text{th}}$ mode, the maximum spin pumping signal is obtained in the NC,
while for the $f_{\text{h}}$ mode, $E_{\text{ac}}$ vanishes under the NC. This
result is explained by Fig. 4(c) illustrating the spin-precessional motion at
$f_{\text{th}}$ (blue) and $f_{\text{h}}$ (red) in the center of the NC. For
$f_{\text{th}}$, the SOT excites a conically precessing Néel vector (See
Supplemental video 1 [57]) with a well-defined almost circular polarization
and a cone angle smaller than $\pi/2$. At large $J_{c}$ (SOT), the cone angle
reaches $\pi/2$ and the Néel vector precesses entirely in the $yz$-plane (See
Supplemental video 2 [57]). In this case, the spin-precessional amplitude
$|u(r)|$ is maximum ($|u(r)|=1$). In Fig. 4(d) we depict the quantity
$j_{s,ac}=|u(r)|\sqrt{1-|u(r)|^{2}}$ as a function of $|u(r)|$. The striking
feature of $j_{s,ac}$ is the maximum ac spin pumping at
$|u|_{\text{max}}=\sqrt{2}/2$. Consistent with the macrospin approximation
applied in Ref. 15, we observe zero spin pumping inside the NC when a large
$J_{c}$ gives rise to $|u(r)|=1$. But, the zero-spin pumping is no longer true
for the discovered Slonzweski waves propagating into the AFM. Outside the NC a
region is found where due to damping the amplitude
$|u|_{\text{max}}=\sqrt{2}/2$ is realized and maximum ac spin pumping occurs
remotely from the NC. The predicted excitation of propagating THz-ASSWs makes
uniaxial AFMs a promising material for SHNOs as they exhibit a lower threshold
current and enable 3 orders of magnitude larger values $|E_{ac}|$ than the
earlier considered biaxial AFMs [14, 15]. In Eq. 9, we find $|E_{ac}|$ to be
proportional to the operational frequency $\omega$, whereas for biaxial AFMs
one finds
$|E_{ac}|\propto\frac{\omega_{e}\omega_{E}}{4\sqrt{\alpha^{2}\omega_{E}^{2}+\omega^{2}}}$.
The latter value is orders of magnitude smaller due to small biaxial
anisotropy $\omega_{e}$.

Figure 4: Spin pumping induced ac output electric field for modes (a)
$f_{\text{th}}$ at threshold and (b) $f_{\text{h}}$ above threshold. Red
circule marks the NC. (c) Time dependent dynamic Néel vector for
$f_{\text{th}}$ (blue) and $f_{\text{h}}$ (red). The arrows indicate the
direction of time. (d) The dimensionless parameter
$j_{s,ac}=|u|\sqrt{1-|u|^{2}}$ as a function of the spin wave amplitude $|u|$.

Conclusion.— We have derived and solved the spin wave equation for a realistic
SHNO considering AFMs with uni-axial anisotropy. We find propagating THz
Slonzewski waves which allow one to generate large ac electrical fields via
spin pumping into an adjacent heavy metal film. At the same time,
synchronization between neighboring SHNOs via propagating SWs might arise. The
predicted waves hence open unprecedented pathways for dc current driven THz
signal generation by means of magnons in thin-film AFMs.

Acknowledgment.— We thank SNSF for financial support via grant 177550.

## References

  * Nagatsuma  _et al._ [2016] T. Nagatsuma, G. Ducournau, and C. C. Renaud, Nat. Photonics 10, 371 (2016).
  * Ma  _et al._ [2019] Z. T. Ma, Z. X. Geng, Z. Y. Fan, J. Liu, and H. D. Chen, Research 2019, 1 (2019).
  * Sirtori [2002] C. Sirtori, Nature 417, 132 (2002).
  * Osborne [2008] I. S. Osborne, Science (80-. ). 320, 1262b (2008).
  * Cheng  _et al._ [2014] R. Cheng, J. Xiao, Q. Niu, and A. Brataas, Phys. Rev. Lett. 113, 057601 (2014).
  * Jungwirth  _et al._ [2016] T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Nat. Nanotechnol. 11, 231 (2016).
  * Baltz  _et al._ [2018] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Rev. Mod. Phys. 90, 015005 (2018).
  * Jungwirth  _et al._ [2018] T. Jungwirth, J. Sinova, A. Manchon, X. Marti, J. Wunderlich, and C. Felser, Nat. Phys. 2018 143 14, 200 (2018).
  * Duine  _et al._ [2018] R. A. Duine, K. J. Lee, S. S. Parkin, and M. D. Stiles, Nat. Phys. 2018 143 14, 217 (2018).
  * Gomonay  _et al._ [2018a] O. Gomonay, V. Baltz, A. Brataas, and Y. Tserkovnyak, Nat. Phys. 2018 143 14, 213 (2018a).
  * Železný  _et al._ [2018] J. Železný, P. Wadley, K. Olejník, A. Hoffmann, and H. Ohno, Nat. Phys. 2018 143 14, 220 (2018).
  * Němec  _et al._ [2018] P. Němec, M. Fiebig, T. Kampfrath, and A. V. Kimel, Nat. Phys. 2018 143 14, 229 (2018).
  * Šmejkal  _et al._ [2018] L. Šmejkal, Y. Mokrousov, B. Yan, and A. H. MacDonald, Nat. Phys. 2018 143 14, 242 (2018).
  * Cheng  _et al._ [2016] R. Cheng, D. Xiao, and A. Brataas, Phys. Rev. Lett. 116, 207603 (2016).
  * Khymyn  _et al._ [2017a] R. Khymyn, I. Lisenkov, V. Tiberkevich, B. A. Ivanov, and A. Slavin, Sci. Rep. 7, 43705 (2017a).
  * Sulymenko  _et al._ [2017] O. R. Sulymenko, O. V. Prokopenko, V. S. Tiberkevich, A. N. Slavin, B. A. Ivanov, and R. S. Khymyn, Phys. Rev. Appl. 8, 064007 (2017).
  * Johansen and Brataas [2017] Ø. Johansen and A. Brataas, Phys. Rev. B 95, 220408 (2017).
  * Chȩciński  _et al._ [2017] J. Chȩciński, M. Frankowski, and T. Stobiecki, Phys. Rev. B 96, 174438 (2017).
  * Khymyn  _et al._ [2017b] R. Khymyn, V. Tiberkevich, and A. Slavin, AIP Adv. 7, 055931 (2017b).
  * Johansen  _et al._ [2018] Ø. Johansen, H. Skarsvåg, and A. Brataas, Phys. Rev. B 97, 054423 (2018).
  * Gomonay  _et al._ [2018b] O. Gomonay, T. Jungwirth, and J. Sinova, Phys. Rev. B 98, 104430 (2018b).
  * Puliafito  _et al._ [2019] V. Puliafito, R. Khymyn, M. Carpentieri, B. Azzerboni, V. Tiberkevich, A. Slavin, and G. Finocchio, Phys. Rev. B 99, 024405 (2019).
  * Lee  _et al._ [2019] D.-K. Lee, B.-G. Park, and K.-J. Lee, Phys. Rev. Appl. 11, 054048 (2019).
  * Lisenkov  _et al._ [2019] I. Lisenkov, R. Khymyn, J. Åkerman, N. X. Sun, and B. A. Ivanov, Phys. Rev. B 100, 100409 (2019).
  * Troncoso  _et al._ [2019] R. E. Troncoso, K. Rode, P. Stamenov, J. M. D. Coey, and A. Brataas, Phys. Rev. B 99, 054433 (2019).
  * Safin  _et al._ [2020] A. Safin, V. Puliafito, M. Carpentieri, G. Finocchio, S. Nikitov, P. Stremoukhov, A. Kirilyuk, V. Tyberkevych, and A. Slavin, Appl. Phys. Lett. 117, 222411 (2020).
  * Parthasarathy  _et al._ [2021] A. Parthasarathy, E. Cogulu, A. D. Kent, and S. Rakheja, Phys. Rev. B 103, 024450 (2021).
  * Shtanko and Prokopenko [2020] O. Shtanko and O. Prokopenko, 2020 IEEE 10th Int. Conf. Nanomater. Appl. Prop. , 01NMM04 (2020).
  * Sulymenko  _et al._ [2018] O. Sulymenko, O. Prokopenko, I. Lisenkov, J. Åkerman, V. Tyberkevych, A. N. Slavin, and R. Khymyn, J. Appl. Phys. 124, 152115 (2018).
  * Khymyn  _et al._ [2018] R. Khymyn, I. Lisenkov, J. Voorheis, O. Sulymenko, O. Prokopenko, V. Tiberkevich, J. Akerman, and A. Slavin, Sci. Rep. 8, 15727 (2018).
  * Stiles and Zangwill [2002] M. D. Stiles and A. Zangwill, Phys. Rev. B 66, 014407 (2002).
  * Ralph and Stiles [2008] D. Ralph and M. Stiles, J. Magn. Magn. Mater. 320, 1190 (2008).
  * Slonczewski [1996] J. Slonczewski, J. Magn. Magn. Mater. 159, L1 (1996).
  * Berger [1996] L. Berger, Phys. Rev. B 54, 9353 (1996).
  * Li  _et al._ [2020a] J. Li, C. B. Wilson, R. Cheng, M. Lohmann, M. Kavand, W. Yuan, M. Aldosary, N. Agladze, P. Wei, M. S. Sherwin, and J. Shi, Nature 578, 70 (2020a).
  * Vaidya  _et al._ [2020] P. Vaidya, S. A. Morley, J. Van Tol, Y. Liu, R. Cheng, A. Brataas, D. Lederman, and E. Del Barco, Science (80-. ). 368, 160 (2020).
  * Daniels  _et al._ [2018] M. W. Daniels, R. Cheng, W. Yu, J. Xiao, and D. Xiao, Phys. Rev. B 98, 134450 (2018).
  * Sulymenko and Prokopenko [2019a] O. R. Sulymenko and O. V. Prokopenko, Proc. Int. Conf. Adv. Optoelectron. Lasers, CAOL 2019-Septe, 533 (2019a).
  * Sulymenko and Prokopenko [2019b] O. R. Sulymenko and O. V. Prokopenko, 2019 IEEE 39th Int. Conf. Electron. Nanotechnology, ELNANO 2019 - Proc. , 132 (2019b).
  * Grollier  _et al._ [2020] J. Grollier, D. Querlioz, K. Y. Camsari, K. Everschor-Sitte, S. Fukami, and M. D. Stiles, Neuromorphic spintronics (2020).
  * Chen  _et al._ [2016] T. Chen, R. K. Dumas, A. Eklund, P. K. Muduli, A. Houshang, A. A. Awad, P. Dürrenfeld, B. G. Malm, A. Rusu, and J. Akerman, Proc. IEEE 104, 1919 (2016).
  * Demidov  _et al._ [2017] V. Demidov, S. Urazhdin, G. de Loubens, O. Klein, V. Cros, A. Anane, and S. Demokritov, Phys. Rep. 673, 1 (2017).
  * Hoefer  _et al._ [2010] M. A. Hoefer, T. J. Silva, and M. W. Keller, Phys. Rev. B 82, 054432 (2010).
  * Mohseni  _et al._ [2013] S. M. Mohseni, S. R. Sani, J. Persson, T. N. Anh Nguyen, S. Chung, Y. Pogoryelov, P. K. Muduli, E. Iacocca, A. Eklund, R. K. Dumas, S. Bonetti, A. Deac, M. A. Hoefer, and J. Åkerman, Science (80-. ). 339, 1295 (2013).
  * Mohseni  _et al._ [2018a] M. Mohseni, M. Hamdi, H. F. Yazdi, S. A. H. Banuazizi, S. Chung, S. R. Sani, J. Åkerman, and M. Mohseni, Phys. Rev. B 97, 184402 (2018a).
  * Hoefer  _et al._ [2005] M. A. Hoefer, M. J. Ablowitz, B. Ilan, M. R. Pufall, and T. J. Silva, Phys. Rev. Lett. 95, 267206 (2005).
  * Mohseni  _et al._ [2018b] S. M. Mohseni, H. F. Yazdi, M. Hamdi, T. Brächer, and S. M. Mohseni, J. Magn. Magn. Mater. 450, 40 (2018b).
  * Kendziorczyk  _et al._ [2014] T. Kendziorczyk, S. O. Demokritov, and T. Kuhn, Phys. Rev. B 90, 054414 (2014).
  * Houshang  _et al._ [2016] A. Houshang, E. Iacocca, P. Dürrenfeld, S. R. Sani, J. Åkerman, and R. K. Dumas, Nat. Nanotechnol. 11, 280 (2016).
  * Kendziorczyk and Kuhn [2016] T. Kendziorczyk and T. Kuhn, Phys. Rev. B 93, 134413 (2016).
  * Awad  _et al._ [2017] A. A. Awad, P. Dürrenfeld, A. Houshang, M. Dvornik, E. Iacocca, R. K. Dumas, and J. Åkerman, Nat. Phys. 13, 292 (2017).
  * Zahedinejad  _et al._ [2020] M. Zahedinejad, A. A. Awad, S. Muralidhar, R. Khymyn, H. Fulara, H. Mazraati, M. Dvornik, and J. Åkerman, Nat. Nanotechnol. 15, 47 (2020).
  * Zahedinejad  _et al._ [2022] M. Zahedinejad, H. Fulara, R. Khymyn, A. Houshang, M. Dvornik, S. Fukami, S. Kanai, H. Ohno, and J. Åkerman, Nat. Mater. 21, 81 (2022).
  * Hals  _et al._ [2011] K. M. D. Hals, Y. Tserkovnyak, and A. Brataas, Phys. Rev. Lett. 106, 107206 (2011).
  * Galkina and Ivanov [2018] E. G. Galkina and B. A. Ivanov, Low Temp. Phys. 44, 618 (2018).
  * Wang  _et al._ [2021] X.-g. Wang, Y.-Z. Nie, L. Chotorlishvili, Q.-l. Xia, J. Berakdar, and G.-h. Guo, Phys. Rev. B 103, 064404 (2021).
  * [57] See Supplemental Material at (link) for derivations of the SW equation as well as the solution to SW equation and derivation of the current density and frequency at the threshold in section S-I. The derivations of continuum model free energy from micromagnetic energy terms used in OOMMF is presented in section S-II. We employ a similar approach to Ref. 75 in section S-II. However, by correctly considering the exchange term for a pair of spins and appropriate choice of the coordinate system we do resolve the issue of manually plugging a factor of 1/2 in Ref. 75. See section S-III for a comparative study of SW dispersion of an AFM obtained by OOMMF and a two sublattice micromagnetic model (Boris Computational Spintronics [76]) as a validation tool for our micromagnetic simulations. See supplemental video 1 and 2 for the visualization of the Slonczewski spin wave obtained from micromagnetic simulations at the threshold and above threshold, respectively. Derivation of Eq. 9 is presented in section S-IV. The scientific colour maps developed by Crameri et. al. [77] is used in this study to prevent visual distortion of the data and exclusion of readers with colour-vision deficiencies [78].
  * Slavin and Kabos [2005] A. N. Slavin and P. Kabos, IEEE Trans. Magn. 41, 1264 (2005).
  * Slavin and Tiberkevich [2008] A. Slavin and V. Tiberkevich, IEEE Trans. Magn. 44, 1916 (2008).
  * Slonczewski [1999] J. C. Slonczewski, J. Magn. Magn. Mater. 195, L261 (1999).
  * Chung  _et al._ [2016] S. Chung, A. Eklund, E. Iacocca, S. M. Mohseni, S. R. Sani, L. Bookman, M. A. Hoefer, R. K. Dumas, and J. Åkerman, Nat. Commun. 7, 1 (2016).
  * Donahue and Porter [1999] M. J. Donahue and D. G. Porter 10.6028/NIST.IR.6376 (1999).
  * Venkat  _et al._ [2018] G. Venkat, H. Fangohr, and A. Prabhakar, J. Magn. Magn. Mater. 450, 34 (2018).
  * Foner [1963] S. Foner, Phys. Rev. 130, 183 (1963).
  * Artman  _et al._ [1965] J. O. Artman, J. C. Murphy, and S. Foner, Phys. Rev. 138, A912 (1965).
  * Saidl  _et al._ [2017] V. Saidl, P. Němec, P. Wadley, V. Hills, R. P. Campion, V. Novák, K. W. Edmonds, F. Maccherozzi, S. S. Dhesi, B. L. Gallagher, F. Trojánek, J. Kuneš, J. Železný, P. Malý, and T. Jungwirth, Nat. Photonics 11, 91 (2017).
  * Lebrun  _et al._ [2018] R. Lebrun, A. Ross, S. A. Bender, A. Qaiumzadeh, L. Baldrati, J. Cramer, A. Brataas, R. A. Duine, and M. Kläui, Nature 561, 222 (2018).
  * Wang  _et al._ [2020] M. Wang, C. Andrews, S. Reimers, O. J. Amin, P. Wadley, R. P. Campion, S. F. Poole, J. Felton, K. W. Edmonds, B. L. Gallagher, A. W. Rushforth, O. Makarovsky, K. Gas, M. Sawicki, D. Kriegner, J. Zubáč, K. Olejník, V. Novák, T. Jungwirth, M. Shahrokhvand, U. Zeitler, S. S. Dhesi, and F. Maccherozzi, Phys. Rev. B 101, 094429 (2020).
  * Fink [1964] H. J. Fink, Phys. Rev. 133, A1322 (1964).
  * Searle and Wang [1968] C. W. Searle and S. T. Wang, J. Appl. Phys. 39, 1025 (1968).
  * Moriyama  _et al._ [2019] T. Moriyama, K. Hayashi, K. Yamada, M. Shima, Y. Ohya, and T. Ono, Phys. Rev. Mater. 3, 051402 (2019).
  * Lebrun  _et al._ [2020] R. Lebrun, A. Ross, O. Gomonay, V. Baltz, U. Ebels, A. L. Barra, A. Qaiumzadeh, A. Brataas, J. Sinova, and M. Kläui, Nat. Commun. 11, 1 (2020).
  * Hamdi  _et al._ [2022] M. Hamdi, F. Posva, and D. Grundler 10.48550/arxiv.2212.11887 (2022).
  * Zhou  _et al._ [2022] Y. Zhou, L. Liao, T. Guo, H. Bai, M. Zhao, C. Wan, L. Huang, L. Han, L. Qiao, Y. You, C. Chen, R. Chen, Z. Zhou, X. Han, F. Pan, and C. Song, Nat. Commun. 13, 3723 (2022).
  * Li  _et al._ [2020b] X. Li, L. Shen, Y. Bai, J. Wang, X. Zhang, J. Xia, M. Ezawa, O. A. Tretiakov, X. Xu, M. Mruczkiewicz, M. Krawczyk, Y. Xu, R. F. Evans, R. W. Chantrell, and Y. Zhou, npj Comput. Mater. 6, 1 (2020b).
  * Lepadatu [2020] S. Lepadatu, J. Appl. Phys. 128, 243902 (2020).
  * [77] F. Crameri, Zenodo 10.5281/ZENODO.5501399.
  * Crameri  _et al._ [2020] F. Crameri, G. E. Shephard, and P. J. Heron, Nat. Commun. 11, 1 (2020).


# Unsupervised Constrained Community Detection via
Self-Expressive Graph Neural Network

Sambaran Bandyopadhyay Sambaran was affiliated to Indian Institute of Science
and IBM Research when the work was conducted. Amazon Science
Bangalore, India  Vishal Peter Both the authors contributed equally. Indian
Institute of Science
Bangalore, India

###### Abstract

Graph neural networks (GNNs) are able to achieve promising performance on
multiple graph downstream tasks such as node classification and link
prediction. Comparatively lesser work has been done to design GNNs which can
operate directly for community detection on graphs. Traditionally, GNNs are
trained on a semi-supervised or self-supervised loss function and then
clustering algorithms are applied to detect communities. However, such
decoupled approaches are inherently sub-optimal. Designing an unsupervised
loss function to train a GNN and extract communities in an integrated manner
is a fundamental challenge. To tackle this problem, we combine the principle
of self-expressiveness with the framework of self-supervised graph neural
network for unsupervised community detection for the first time in literature.
Our solution is trained in an end-to-end fashion and achieves state-of-the-art
community detection performance on multiple publicly available datasets.

##  1 Introduction

Graphs or networks are ubiquitous in our daily life. Graph representation
learning [Perozzi et al., 2014, Hamilton et al., 2017b] is the task of mapping
different components of a graph (such as nodes, edges or the entire graph) to
a vector space to facilitate downstream graph mining tasks. Among various
types of graph representation techniques, graph neural networks (GNNs) [Wu et
al., 2020] have received significant attention as they are able to apply
neural networks directly on the graph structure. Most of the GNNs can be
represented in the form of a message passing network, where each node updates
its vector representation by aggregating messages from neighboring nodes with
its own [Gilmer et al., 2017, Hamilton et al., 2017a]. GNNs are traditionally
trained in a semi-supervised way [Kipf and Welling, 2017] on a node
classification loss when a subset of node labels are available. More recently,
unsupervised and self-supervised graph neural networks have been proposed
where a reconstruction loss [Kipf and Welling, 2016, Bandyopadhyay et al.,
2020] or noise contrastive loss [Veličković et al., 2019, Zhu et al., 2020] is
used to train the networks.

Community detection is one of the most important tasks for network analysis
and has been studied for decades in classical network analysis community
[Fortunato and Hric, 2016, Xie et al., 2013]. However, compared to other tasks
such as node classification [Kipf and Welling, 2017, Veličković et al., 2019]
and link prediction [Kipf and Welling, 2016, Zhang and Chen, 2018], community
detection has not been explored much in the framework of graph neural
networks. Being inherently unsupervised in nature, it is challenging to train
GNNs for community detection directly. Traditionally, methods have been
proposed where a graph representation learning algorithm is trained on a
generic unsupervised loss and then a clustering algorithm is applied as a
post-processing step to discover communities [Perozzi et al., 2014,
Bandyopadhyay et al., 2019]. Such approaches are sub-optimal in nature as the
node representation learning module and the clustering algorithm work
independently. More recently, there have been efforts to train graph neural
networks directly for community detection in graphs [Bo et al., 2020, Zhang et
al., 2020] (Section 2).

In contrast to a fully unsupervised approach, a graph neural architecture is
proposed in [Chen et al., 2019] for a supervised version of community
detection . In classical machine learning, constraint clustering has been
shown to be very efficient where must-link or no-link constraints are given as
input [Wagstaff et al., 2001]. But, obtaining direct ground truth community
labels or such pair-wise constraints is expensive for real-world networks. In
this paper, we aim to derive such constraints in an unsupervised way, by using
the principle of self-expressiveness of data [Ji et al., 2014]. This allows to
express each data point by a linear combination of other data points which
potentially lie in the same subspace. The principle of self-expressiveness has
been successfully applied in computer vision and image processing for object
detection and segregation [Zhang et al., 2019a, Li and Vidal, 2015]. However,
inherent computational demand to build pair-wise similarity matrix and
subsequent use of spectral clustering makes it infeasible to directly apply
the principle of self-expressiveness and subspace clustering in domains like
graphs where number of nodes can be very large [Elhamifar and Vidal, 2009, Ji
et al., 2014, 2017]. We have taken a different approach in this paper to
address these computational challenges. Our solution uses a self-supervised
GNN and generate node communities from the embeddings obtained. To guide the
generated communities, we use the principle of self-expressiveness on randomly
sampled batches of nodes to generate a set of soft must-link and no-link
constraints.

Following are the novel contributions made in this paper:

  * •

We propose a novel community detection algorithm, referred as SEComm (Self-
Expressive Community detection in graph). To the best of our knowledge, we are
the first in literature to combine the principle of self-expressiveness with
the framework of self-supervised graph neural network for unsupervised
community detection. Our solution is able to use both link structure and the
node attributes of a graph to detect node communities.

  * •

To address the computational issues, our solution uses the principle of self-
expressiveness to generate a set of soft must-link or no-link constraints on a
subset of nodes divided into batches. In contrast to existing literature on
self-expressiveness (which typically applies spectral clustering as a post-
processing step), our solution is trained in an end-to-end fashion.

  * •

To show the merit of the proposed algorithm, we conduct experiments with
multiple publicly available graph datasets and compare the results with a
diverse set of algorithms. SEComm is able to improve the state-of-the-art
performance of unsupervised community detection with a significant margin in
almost all the real-world datasets we used. Model ablation study and
sensitivity analysis give further insights of the algorithm. Source code of
SEComm is available at https://github.com/viz27/SEComm.

##  2 Related Work

We have presented the related work into three categories.

Graph Neural Networks: Graph neural networks have gained significant attention
in last few years with their success on a diverse set of applications [Wu et
al., 2020, Desai et al., 2021]. Typically, GNNs are trained on node-
classification, link prediction and graph reconstruction losses [Kipf and
Welling, 2017, Hamilton et al., 2017a]. Recently, self-supervised learning has
been able to achieve performance close to supervised learning for multiple
downstream tasks [Belghazi et al., 2018, Hjelm et al., 2019]. Extending the
concept of information maximization, DGI [Veličković et al., 2019] and GRACE
[Zhu et al., 2020] have been proposed where information between different
graph entities (graph-level to node-level, corrupted versions of a graph etc.)
are maximized. However, none of the GNNs above handles community detection in
their respective objectives.

Principle of Self-Expressiveness: The concept of self expressiveness was
proposed to cluster data drawn from multiple low dimensional linear or affine
subspaces embedded in a high dimensional space [Elhamifar and Vidal, 2009].
Given enough samples, each data point in a union of subspaces can always be
written as a linear or affine combination of all other points [Elhamifar and
Vidal, 2009, Ji et al., 2014]. Subspace clustering exploits this to build a
similarity matrix, from which the segmentation of the data can be easily
obtained using spectral clustering [Lu et al., 2012, Ji et al., 2014].
Recently, a deep learning based subspace clustering method has been proposed
where an encoder is used to map data to some embedding space before building
the pair-wise similarity matrix and applying spectral clustering [Ji et al.,
2017]. However, inherent computational demand to build pair-wise similarity
matrix and subsequent use of spectral clustering makes it infeasible to
directly apply the principle of self-expressiveness in domains like graphs
where the number of nodes can be very large.

Community Detection with GNNs: As explained in Section 1, a disjoint approach
of applying clustering on node embeddings obtained by some representation
algorithm is inherently sub-optimal in nature [Bandyopadhyay et al., 2020]. In
[Zhang et al., 2019b], the authors have used an adaptive graph convolution
method that performs high-order graph convolution to obtain smooth node
embeddings that capture global cluster structure. The node embeddings obtained
are subsequently used to detect communities using spectral clustering. In [Sun
et al., 2019], a probabilistic generative model is proposed to learn community
membership and node representation collaboratively. More recently, researchers
have tried to propose GNN algorithms that can operate directly for community
detection in a graph [Bo et al., 2020]. In [Zhang et al., 2020], authors
propose to derive node community membership in the hidden layer of an encoder
and introduced a community-centric dual decoder to reconstruct network
structures and node attributes in an unsupervised fashion. Our work is towards
this direction of obtaining node communities directly in the framework of
graph neural networks.

##  3 Problem Formulation

Let us denote an input graph by $G=(V,E,X)$, where $V=\\{1,2,\cdots,N\\}$ is
the set of nodes and $E\subseteq V\times V$ is the set of edges. We assume
that each node has some attribute values present in a vector
$x_{i}\in\mathbb{R}^{F}$ and
$X=[x_{1}\;x_{2}\;\cdots\;x_{N}]^{T}\in\mathbb{R}^{N\times F}$ is the node
attribute matrix of the graph. The goal of our work is to learn a function
$f:V\mapsto[K]$, where $[K]=\\{1,2,\cdots,K\\}$ is the set of community (or
cluster) indices, to map each node to a community by exploiting the link
structure and node attributes of the graph. We want to achieve this without
having any ground truth community information of a node. Intuitively, nodes
which are closely connected in the graph or have similar attributes, should be
members of the same community. Important notations in the paper are summarized
in Table 1.

Notations | Explanations
---|---
$G=(V,E,X)$ | Input graph
$i,j\in\\{1,2,\cdots,N\\}=V$ | Indices over nodes
$X\in\mathbb{R}^{N\times F}$ | Node feature matrix
$Z\in\mathbb{R}^{N\times F^{\prime}}$ | Node embedding matrix
$M,\;P$ | Batch size and number of batches sampled
$S_{ij}\in[0,1]$ | Similarity between nodes $i$ and $j$
$\mathfrak{S}\subseteq V\times V$ | node pairs whose similarities are computed
$\mathfrak{S}_{ext}\subseteq\mathfrak{S}$ | After filtering out with $\theta_{low}$ and $\theta_{high}$
$C_{i}\in\mathbb{R}^{K}$ | Community membership vector for $i$th node
$C\in\mathbb{R}^{N\times K}$ | Output node community membership matrix
$W_{SS}$, $W_{MLP}$ | Parameters of SS-GNN and MLP modules

Table 1: Different notations used in the paper Figure 1: Training phases of
SEComm

##  4 Our Solution: SEComm

There are multiple steps in our proposed solution SEComm as shown in Figure 1.
We discuss each of them.

###  4.1 Self-supervised Node Embedding

The first step of SEComm is to learn node representation in an unsupervised
way. Self-supervised learning [Hjelm et al., 2019] has been used recently for
obtaining both node embeddings [Veličković et al., 2019, Zhu et al., 2020] and
graph-level embeddings [Sun et al., 2020]. Potentially, any self-supervised
differentiable approach to obtain node representation can be integrated with
our solution. In the following, we have adopted the principle of mutual
information maximization between two corrupted versions of the given graph,
motivated from [Veličković et al., 2019, Zhu et al., 2020], which is then used
to formulate the final objective of SEComm in Section 4.3.

Given the input graph $G=(V,E,X)$, two graph views $G_{1}$ and $G_{2}$ are
generated from it by employing a corruption function. The corruption function
randomly removes a small portion of edges from the input graph and also
randomly masks a fraction of dimensions with zeros in node features. The
vertex sets of $G_{1}$ and $G_{2}$ remain the same. These views are used for
contrastive learning at both graph topology and node feature levels. We use a
GCN encoder to generate node embeddings for both $G_{1}$ and $G_{2}$. For a
graph $G$, GCN encoder derives node representations as follows:

$\displaystyle
Z=f(X,A)=\text{ReLU}(\hat{A}\;\text{ReLU}(\hat{A}XW^{(0)})\;W^{(1)})$ (1)

where each row of $Z\in\mathbb{R}^{|V|\times F^{\prime}}$ contains the
corresponding node representation. $A$ is the adjacency matrix of the graph
$G$. We compute $\tilde{A}=A+I$, where $I\in\mathbb{R}^{|V|\times|V|}$ is the
identity matrix and the degree diagonal matrix $\tilde{D_{ii}}$ with
$\tilde{D_{ii}}=\sum\limits_{j\in V}\tilde{A}_{ij}$, $\forall i\in V$. We set
$\hat{A}=\tilde{D}^{-\frac{1}{2}}\tilde{A}\tilde{D}^{-\frac{1}{2}}$. $W^{(0)}$
and $W^{(1)}$ are the trainable parameter matrices of GCN. Let us use $Z_{1}$
and $Z_{2}$ to denote the node embedding matrices for the two views $G_{1}$
and $G_{2}$ obtained from the GCN encoder (parameters shared).

Next, the following noise contrastive objective (via a discriminator) is used.
For any node $i\in V$, let us denote the corresponding nodes in $G_{1}$ and
$G_{2}$ as $G_{1}(i)$ and $G_{2}(i)$ respectively. For each $i\in V$, the pair
$(G_{1}(i),G_{2}(i))$ is considered as a positive example. Negative examples
are sampled from both the views for each node $i\in V$. More formally, we
randomly select a set of nodes $V_{-i}=\\{j\in V\;|\;j\neq i\\}$ such that
$|V_{-i}|=N_{-}$ (number of negative samples), $\forall i\in V$. Both
$(G_{1}(i),G_{1}(j))$ and $(G_{1}(i),G_{2}(j))$ are considered as negative
examples. Then, the following objective function is minimized:

$\displaystyle\underset{W_{SS}}{\text{min}}\;\mathcal{L}_{SS}\;=\;$
$\displaystyle\sum\limits_{i\in V}\Bigg{[}-\frac{cos(Z_{1i},Z_{2i})}{\tau}$
(2) $\displaystyle+\log\Big{(}\sum\limits_{j\in
V_{-i}}e^{\frac{cos(Z_{1i},Z_{1j})}{\tau}}+e^{\frac{cos(Z_{1i},Z_{2j})}{\tau}}\Big{)}\Bigg{]}$

where $Z_{1i}$ and $Z_{2i}$ denotes the $i$th row of $Z_{1}$ and $Z_{2}$
respectively, $cos()$ is the cosine similarity between the two embeddings and
$\tau$ is a temperature parameter. This essentially maximizes the agreement
between the embeddings of $i$th node in two views. Both link structure and
node attributes of the graph are considered because of the use of GCN encoder
in Equation 1. Node embedding matrix $Z$ of the original input graph $G$ can
be obtained by Equation 1, where the purpose of $G_{1}$ and $G_{2}$ were just
to enable the training of the self-supervised loss in Equation 2.

###  4.2 Learning Node Similarities through Self-Expressiveness

Vector node representations obtained in Section 4.1 are quite generic. There
is no guarantee that similarities between the nodes captured through such
embeddings are suitable to discover communities in the graph. In contrast to
node classification or other supervised tasks, lack of any ground truth in
training for community detection makes the problem highly non-trivial. To
tackle this, we use the principle of self-expressiveness [Elhamifar and Vidal,
2009, Ji et al., 2017] which aims to approximate a data point by a sparse
linear sum of a subset of other points, which stays in the same subspace. This
is more prevalent when the graph has large number of nodes and embedding
dimension is also reasonably high [Elhamifar and Vidal, 2009]. Based on the
contribution of a point to reconstruct some other point, it is possible to
learn a pair-wise similarity matrix using this principle. Such a pair-wise
similarity can guide the generation of communities from the node embeddings
obtained by the self-supervised layer. However, computation of pair-wise
similarity matrix for a graph can be too expensive as it needs $O(N^{2})$
storage and processing. Hence, we propose a batch-wise learning procedure, as
discussed later.

Given the node embedding matrix $Z\in\mathbb{R}^{N\times F^{\prime}}$ obtained
from Section 4.1, we want to derive a node similarity matrix
$S\in\mathbb{R}^{N\times N}$ using the principle of self-expressiveness. For
each node $i\in V$, we first try to express $Z_{i}$ ($i$th row of $Z$) by a
linear sum of few other node embeddings $Z_{j}$, $j\neq i$. So,
$Z_{i}=\sum\limits_{j\in V}q_{ij}Z_{j}$, where $q_{ij}$ is the $(i,j)$th
element of a coefficient matrix $Q\in\mathbb{R}^{N\times N}$ and we enforce
$Q_{ii}=0$ to avoid the trivial solution of Q being assigned to a identity
matrix. We need to learn this coefficient matrix $Q$ which will be used to
generate similarity matrix $S$. It can be shown [Ji et al., 2014] under the
assumption of subspace independence that, by minimizing certain norms of $Q$,
it is possible to have a block-diagonal structure (up to a permutation) of
$Q$. In that case, each block in $Q$ would contain nodes which belong to the
same subspace. This can be posed as the following optimization problem.

$\displaystyle\underset{Q}{\text{min}}\;||Q||_{p}\;\;\text{such
that,}\;Z=QZ;\;\text{diag}(Q)=0$ (3)

where $||Q||_{p}$ is $p$th matrix norm of $Q$ and diag$(Q)$ denotes the
diagonal entries of $Q$. Based on the choice in some existing literature [Lu
et al., 2012, Ji et al., 2014], we use square Frobenius matrix norm for our
implementation. However, exact reconstruction of $Z$ using the this principle
may not be possible. So, we relax the hard constraint $Z=QZ$ with square
Frobenius norm of $(Z-QZ)$ (soft constraint). This gives us the following
objective function.

$\displaystyle\underset{Q}{\text{min}}\;\mathcal{L}_{SE}\;=\;||Z-QZ||_{F}^{2}+\lambda_{1}||Q||_{F}^{2}\;\;\text{such
that,}\;\text{diag}(Q)=0$ (4)

where $\lambda_{1}$ is a weight parameter of this optimization.

In principle, while a pairwise similarity matrix $S$ can be constructed
trivially as $Q$+$Q^{T}$, many heuristics have been proposed to improve the
clustering performance of $S$ (when using methods such as spectral clustering
directly on $S$). We follow the heuristics proposed in [Ji et al., 2014] to
construct the node similarity matrix $S$ as:

  1. 1.

$Q^{\prime}=\frac{1}{2}(Q+Q^{T})$

  2. 2.

Compute the $r$ rank SVD of $Q^{\prime}$, ie. $Q^{\prime}=U\Sigma V^{T}$,
where $r=dK+1$, $K$ is the number of communities and $d$ is the maximal
intrinsic dimension of subspaces which is set to 4 in all our experiments.

  3. 3.

Compute $L=U\Sigma^{\frac{1}{2}}$ and normalize each row of $L$ to have unit
norm.

  4. 4.

Set negative values in $L$ to zero to obtain $L^{\prime}$.

  5. 5.

Construct similarity matrix $S$ as $S=(L^{\prime}+L^{\prime
T})/||L||_{\infty}$, so that $s_{ij}\in[0,1]$.

As mentioned before, an inherent difficulty to compute the pair-wise
similarity matrix is the computation and storage of $N\times N$ dimensional
matrix $S$. So, instead of computing this matrix for all pairs of nodes, we
use batch-wise learning. We sample batches of randomly selected nodes with
batch size $M$, where $M\leq N$. We train the loss in Equation 4 for each
batch. The required computation in each batch is $O(M^{3})$ (for solving
Equation 4 and the subsequent use of SVD decomposition) which is much lesser
than $O(N^{3})$ for a significantly smaller $M$. However, the problem with
this approach is that one would not get complete similarity matrix for the
graph. It only computes $s_{ij}$ if nodes $i$ and $j$ belong to a same batch.
Let us denote $\mathfrak{S}$ to be the set of node pairs for which the
similarity is computed in the batch-wise learning. Clearly,
$\mathfrak{S}\subsetneq V\times V$ and $|\mathfrak{S}|<<N^{2}$ (when $M<N)$.
This makes it difficult to use with spectral clustering, as most of the
subspace clustering algorithms do [Ji et al., 2017]. But as explained next,
our overall solution does not need all the node-pair similarities. Rather, it
filters the existing similarities computed with the batch-wise solution using
a simple trick explained next.

###  4.3 Constrained Node Community Detection

Instead of applying expensive spectral clustering on the complete matrix $S$
as a post processing step to find node clusters, we use a neural network based
solution which is significantly more scalable. We use a fully-connected
multilayer perceptron (MLP) with the set of trainable parameters $W_{MLP}$ to
map each node embedding to its corresponding soft community membership vector
as follows.

$C_{i}=\text{Softmax}(\text{MLP}(Z_{i}))\in\mathbb{R}^{K}$ (5)

where the MLP maps each $Z_{i}\in\mathbb{R}^{F^{\prime}}$ to a $K$ dimensional
vector, $K$ is the number of communities. We assume to know $K$ beforehand.
The softmax layer converts the $K$ dimensional vector to a probability
distribution such that $c_{ik}$ ($k$th element of $C_{i}$) denotes the
probability that $i$th node belongs to $k$th community, $\forall k=[K]$.
Equation 5 ensures that nodes having similar embeddings will be mapped to
similar positions in the $(K-1)$ dimensional probability simplex. However,
relying completely on embeddings to detect communities is not desirable since
the embeddings are generated with generic objectives. So, they may not be
optimal to generate node communities. Hence, we use the information learned in
node-pair similarities in Section 4.2 to guide both the detection of node
communities by training the parameters of MLP in Equation 5 and updating node
embeddings.

Let us form the community membership matrix $C=[C_{1},\cdots,C_{N}]^{T}$
$\in\mathbb{R}^{N\times K}$. If the complete node similarity matrix $S$ is
available, one may try to minimize the following objective.

$\displaystyle\underset{W_{SS},\;W_{MLP}}{\text{min}}||CC^{T}-S||_{F}^{2}$ (6)

There are multiple drawbacks present in the objective function above. First,
it needs us to compute the complete node similarity matrix $S$ in Section 4.2
which prevents the batch-learning mechanism explained before. Next, the
computation involved is $O(N^{2})$ in Equation 6. Further, there is another
issue if one wants to use all pair-wise node similarities in $S$ to guide the
community detection. Due to noise in the dataset, many of the pair-wise
similarities may not reflect the actual similarities between the nodes. The
similarity values which are around 0.5 neither express a strong similarity nor
a strong dissimilarity between a node pair. So they are less informative
compared to the similarity values which are close to 0 or 1. But they can
still influence the parameters of the neural network because of Equation 4.

Hence, instead of considering all the pair-wise similarity values, we only
consider the ones in $\mathfrak{S}$ computed over the batches as discussed in
Section 4.2. Further, we have observed experimentally in Section 5.6 that for
a larger dataset, it is okay even if some nodes are not part of any of the
batches selected randomly. Thus, the number of batches can be significantly
smaller than $\frac{N}{M}$ for a larger dataset. We also introduce two
thresholds $\theta_{low}$ and $\theta_{high}$ to use only those node-pair
similarities which are extreme in their values, thus more informative in
nature. We set $0<\theta_{low}\leq\theta_{high}<1$. We also set
$\theta_{high}=1-\theta_{low}$, as this choice works well in the experiments
and also reduces the number of hyperparameters. Let us introduce the set
$\mathfrak{S}_{ext}\subseteq\mathfrak{S}$ as follows.

$\mathfrak{S}_{ext}=\Big{\\{}(i,j)\in\mathfrak{S}\;|\;S_{ij}\leq\theta_{low}\;\text{or}\;S_{ij}\geq\theta_{high}\Big{\\}}$
(7)

Here, a node pair $(i,j)$ in $\mathfrak{S}_{ext}$ should be roughly
constrained to be in the same cluster when $S_{ij}$ value is very high or in
different clusters when $S_{ij}$ is very low. Thus, we derive a set of soft
version of must-link and no-link constraints in an unsupervised way to guide
the formation of communities. With these, we formulate the following
optimization to detect the communities:

$\displaystyle\underset{W_{SS},\;W_{MLP}}{\text{min}}\sum\limits_{(i,j)\in\mathfrak{S}_{ext}}\Big{(}C_{i}^{T}C_{j}-S_{ij}\Big{)}^{2}$
(8)

By considering only the node pairs in $\mathfrak{S}_{ext}$, we are able to
ignore the pairs which are neither too similar nor too dissimilar, to
contribute to the learning of community memberships. As $C_{i}$ is a
probability distribution over all the $K$ communities for a node $i$, we want
to avoid trivial community formations where each node is assigned to all the
communities with roughly uniform probabilities, or all the nodes are assigned
to a single community [Bianchi et al., 2020]. So, we update the main objective
in Equation 8 as:

$\displaystyle\underset{W_{SS},\;W_{MLP}}{\text{min}}\;\mathcal{L}_{Com}=$
$\displaystyle\sum\limits_{(i,j)\in\mathfrak{S}_{ext}}\Big{(}C_{i}^{T}C_{j}-S_{ij}\Big{)}^{2}$
(9)
$\displaystyle+\lambda_{2}\biggl{\lvert}\biggl{\lvert}\frac{C^{T}C}{||C^{T}C||_{F}}-\frac{I_{K}}{\sqrt{K}}\biggl{\lvert}\biggl{\lvert}_{F}^{2}$

The second component in the equation above ensures that communities are close
to orthogonal and they are balanced in size. Please note that due to the use
of neural network to generate community membership for each node in Equation
5, the optimization in Equation 9 is not a discrete optimization. Rather, we
solve it with respect to the parameters $W_{SS}$ of the self-supervised layer
(Eq. 2) and $W_{MLP}$ of the MLP (Eq. 5). The total loss to train the node
embeddings and community detection can be written as a weighted sum of self-
supervised loss and community detection loss, which is shown below:

$\displaystyle\underset{W_{SS},\;W_{MLP}}{\text{min}}\;\mathcal{L}_{total}=\alpha\mathcal{L}_{SS}+\mathcal{L}_{Com}$
(10)

where $\alpha$ is a weight factor of the optimization. The node-pair
similarity values are obtained by solving the batch-learning technique in
Section 4.2. The overall algorithm SEComm proceeds in an iterative way by
solving the self-expressive layer for each batch and then updating the
parameters of the neural network by minimizing Equation 10. The pseudo code of
SEComm is presented in 1.

Algorithm 1 SEComm \- Self-Expressive Community Detection

1:Input: The graph $G=(V,E,X)$, $|V|=N$, $K$: Number of communities in the
graph, $M$: Batch size for the self-expressive layer, $P$: Number of batches
used for training self expressive layer, Thresholds $\theta_{low}$ and
$\theta_{high}$.

2:Output: Community membership vector $C_{i}\in\mathbb{R}^{K}$ for each node
$i\in V$.

3:Initialize the parameters of the self-supervised GNN and clustering MLP (in
Eq. 5).

4:pre-training step: Obtain node embeddings $Z\in\mathbb{R}^{N\times
F^{\prime}}$ by training the self-supervised GNN.

5:Initialize $\mathfrak{S}$ as empty.

6:for $batch\in\\{1,2,\cdots,P\\}$ do

7: Sample a batch of $M$ nodes from $V$

8: Learn the pair-wise node similarity matrix $S_{m}$ for the selected nodes
by optimizing Eq. 4.

9: Add all the node-pairs from the batch to $\mathfrak{S}$.

10:end for

11:Construct $\mathfrak{S}_{ext}$ according to Eq. 7

12:for $iter\in\\{1,2,\cdots,T\\}$ do

13: Generate node embedding matrix $Z\in\mathbb{R}^{N\times F^{\prime}}$ using
the self-supervised GNN.

14: Generate cluster membership vector $C_{i}\in\mathbb{R}^{K}$ for each node
$i\in V$.

15: Update the parameters of the GNN and clustering MLP by optimizing Eq. 10

16:end for

###  4.4 Training and Analysis of SEComm

We use ADAM with default parameterization to solve the optimization
formulations in Equations 4 and 10. For the self-expressive loss in Equation
4, we train until the loss saturates. For the total loss in Equation 10, we
particularly focus on the saturation of the regularization
$\biggl{\lvert}\biggl{\lvert}\frac{C^{T}C}{||C^{T}C||_{F}}-\frac{I_{K}}{\sqrt{K}}\biggl{\lvert}\biggl{\lvert}_{F}^{2}$.
Experimentally, using the convergence of this component explicitly as a
stopping criteria for SEComm gives slightly better result for all the
datasets, than checking the total convergence. But as explained in Section
5.3, different components of the loss function have similar contributions.
Hence they saturate almost in the same time for most of the cases. This can
also be observed in Section 5.5.

Time Complexity: Time complexity of the self-supervised GNN in Section 4.1 is
$O(|E|+NFF^{\prime}N_{-})$, where $N_{-}$ is the number of negative samples
used. The self-expressive layer takes another $O(PM^{3})$ time, where $P$ is
the number of batches sampled, and $M$ is the size of each batch. Finally,
generating community membership takes $O(NK^{2})$ time because of solving the
loss in Equation 9. Thus, the overall run time of each iteration of SEComm is
linearly dependent on the number of nodes and number of edges in the graph.

##  5 Experimental Evaluation

This section presents the details of the experiments that we conducted and the
analysis of the results.

Dataset | #Nodes | #Edges | #Features | #Labels
---|---|---|---|---
Cora | 2,708 | 5,429 | 1,433 | 7
Citeseer | 3,327 | 4,732 | 3,703 | 6
Pubmed | 19,717 | 44,338 | 500 | 3
Wiki | 2,405 | 17,981 | 4,973 | 17
Physics | 34,493 | 247,962 | 8,415 | 5

Table 2: Summary of the datasets used

###  5.1 Datasets Used

To show the merit of SEComm, we conduct experiments on 5 publicly available
graph datasets [Kipf and Welling, 2017, Zhang et al., 2019b]. Different
statistics of the datasets are summarized in Table 2. Cora, Citeseer and
Pubmed are citation datasets where nodes correspond to papers and are
connected by an edge if one cites the other. Wiki is a collection of webpages
where nodes are webpages and are connected if one links to other. Physics is a
co-authorship network where nodes are authors, that are connected by an edge
if they have co-authored a paper [Shchur et al., 2018]. Each of these datasets
have attribute vector associated with each node. They also have ground truth
community membership of each node, which we use to evaluate the performance of
our proposed and baseline algorithms.

Methods | Input | Cora | Citeseer | Pubmed | Wiki
---|---|---|---|---|---
|  | Acc% | NMI% | F1% | Acc% | NMI% | F1% | Acc% | NMI% | F1% | Acc% | NMI% | F1%
k-means | Feature | 34.65 | 16.73 | 25.42 | 38.49 | 17.02 | 30.47 | 57.32 | 29.12 | 57.35 | 33.37 | 30.20 | 24.51
Spectral-f | Feature | 36.26 | 15.09 | 25.64 | 46.23 | 21.19 | 33.70 | 59.91 | 32.55 | 58.61 | 41.28 | 43.99 | 25.20
Spectral-g | Graph | 34.19 | 19.49 | 30.17 | 25.91 | 11.84 | 29.48 | 39.74 | 3.46 | 51.97 | 23.58 | 19.28 | 17.21
DeepWalk | Graph | 46.74 | 31.75 | 38.06 | 36.15 | 9.66 | 26.70 | 61.86 | 16.71 | 47.06 | 38.46 | 32.38 | 25.74
DNGR | Graph | 49.24 | 37.29 | 37.29 | 32.59 | 18.02 | 44.19 | 45.35 | 15.38 | 17.90 | 37.58 | 35.85 | 25.38
GAE | Both | 53.25 | 40.69 | 41.97 | 41.26 | 18.34 | 29.13 | 64.08 | 22.97 | 49.26 | 17.33 | 11.93 | 15.35
VGAE | Both | 55.95 | 38.45 | 41.50 | 44.38 | 22.71 | 31.88 | 65.48 | 25.09 | 50.95 | 28.67 | 30.28 | 20.49
MGAE | Both | 63.43 | 45.57 | 38.01 | 63.56 | 39.75 | 39.49 | 43.88 | 8.16 | 41.98 | 50.14 | 47.97 | 39.20
ARGE | Both | 64.00 | 44.90 | 61.90 | 57.30 | 35.00 | 54.60 | 59.12 | 23.17 | 58.41 | 41.40 | 39.50 | 38.27
ARVGE | Both | 63.80 | 45.00 | 62.70 | 54.40 | 26.10 | 52.90 | 58.22 | 20.62 | 23.04 | 41.55 | 40.01 | 37.80
AGC | Both | 68.92 | 53.68 | 65.61 | 67.00 | 41.13 | 62.48 | 69.78 | 31.59 | 68.72 | 47.65 | 45.28 | 40.36
GUCD | Both | 50.5 | 32.3 | NA | 54.47 | 27.43 | NA | 63.13 | 26.98 | NA | NA | NA | NA
SEComm | Both | 75.92 | 56.04 | 73.94 | 69.82 | 42.53 | 60.25 | 74.49 | 36.50 | 73.50 | 53.10 | 51.38 | 44.48
Rank | (SEComm) | 1 | 1 | 1 | 1 | 1 | 2 | 1 | 1 | 1 | 1 | 1 | 1

Table 3: Performance of Community Detection by SEComm and other baseline
algorithms

Methods | Acc% | NMI% | F1% | Runtime (sec.)
---|---|---|---|---
k-means | 44.20 | 44.46 | 37.63 | 28
ARGE | 60.67 | 51.77 | 62.32 | 2112
ARVGE | 61.28 | 53.49 | 65.47 | 2221
AGC | 75.36 | 58.19 | 60.72 | 6931
SEComm | 77.93 | 56.08 | 76.42 | 624

Table 4: Community Detection on Physics Dataset

(a) Cora

(b) Wiki

Figure 2: Loss vs Accuracy progression during training

###  5.2 Baseline Algorithms

We use a diverse set of baselines to compare the performance of SEComm. We
divide them into the following categories.

Using only Node Features: As each node is associated with some attribute
vectors, we use k-means and spectral clustering (Spectral-f) algorithms
directly on the node attributes to cluster nodes into different communities.
Naturally, these approaches ignore the graph structure completely.

Using only Graph Structure: We also use spectral clustering (Spectral-g) on
the graph structure. Here we consider adjacency matrix of a graph as the
similarity matrix between the nodes. We also use popular unsupervised node
embedding techniques DeepWalk [Perozzi et al., 2014], which is a random walk
based technique and DNGR [Cao et al., 2016], which is an auto-encoder based
technique.

Using both Node Features and Graph Structure: We use a set of unsupervised
graph neural network based techniques. GNN based approaches are naturally able
to use both link structure and node attribute of the graph. They are: graph
autoencoder (GAE) and graph variational autoencoder (VGAE) [Kipf and Welling,
2016], marginalized graph autoencoder (MGAE) [Wang et al., 2017],
adversarially regularized graph autoencoder (ARGE) and variational graph
autoencoder (ARVGE) [Pan et al., 2018]. These methods typically learn the node
embeddings and use clustering on the embeddings as a post processing step.
Finally, we use two recently proposed community detection methods - AGC [Zhang
et al., 2019b], which uses high-order graph convolution to get node embeddings
and detect communities via spectral clustering on the embeddings and GUCD
[Zhang et al., 2020], which uses an auto-encoder based framework to obtain
direct community assignments for every node.

###  5.3 Experimental Setup

Our proposed algorithm SEComm generates community membership of each node in a
graph in the framework of graph neural networks. As each node has a single
ground truth community membership in all the datasets that we use, we consider
the index of the maximum value of $C_{i}\in\mathbb{R}^{K}$ (from Equation 5)
as the community index of the node generated by SEComm.

There are multiple hyperparameters present in SEcomm. For weight factors in
optimization such as $\lambda_{1}$, $\lambda_{2}$ and $\alpha$, we check the
contribution of different components in a loss function at the beginning of
the algorithm, and set these parameters to values such that effective
contributions of those components become roughly the same. This ensures that
the optimization pays similar importance to different components of SEComm.
For the temperature parameter $\tau$, we use the same values used in the
literature [Zhu et al., 2020]. For threshold parameters $\theta_{low}$
($0<\theta_{low}\leq 0.5$), we set it to 0.5 for relatively smaller datasets
as we do not want to discard any information for them. For Pubmed, we set it
to 0.05 as considering more node-pair similarity values adds noise and also
increases runtime of SEComm. However on Physics, the training convergence is
not smooth when we set $\theta_{low}$ to a very small number. So, we set it to
$0.3$ on this dataset. As mentioned in Section 4.3, we set
$\theta_{high}=1-\theta_{low}$. We have also conducted sensitivity analysis of
SEComm with respect to some of these hyperparameters in Section 5.6.

###  5.4 Results of Community Detection

Tables 3 and 4 show the performance of community detection by different
baseline algorithms and SEComm. We use three popularly used metrics to
evaluate the quality of community detection. They are clustering accuracy
(Acc), normalized mutual information (NMI), and macro F1-score (F1) [Aggarwal
and Reddy, 2014, Zhang et al., 2019b]. We use ground truth community
information of the nodes only to calculate these quality metrics.

While reporting the performance of baseline algorithms for the first four
datasets in Table 3, we have collected the best results from the available
literature [Zhang et al., 2019b, 2020] which adopted the same experimental
setup. We mark some entry as ‘NA’ if the result of that algorithm for a
dataset is not publicly available. For Physics dataset, the baseline results
are not available in the literature. So, we have run and reported results only
for better-performing and diverse subset of baselines in Table 4 with adequate
hyperparameter tuning. Additionally, we have also reported the runtime on
Physics dataset for these algorithms to give more insight about scalability.

We run SEComm 10 times on each dataset and report the average performance.
Tables 3 and 4 show that SEComm is able to achieve state-of-the-art (SOTA)
performance for all the datasets, and for all the metrics, except on
Citeseer-F1% and Physics-NMI% scores, where SEComm is next to AGC. In terms of
performance improvement by clustering accuracy, SEComm is able to improve SOTA
by 10.1% on Cora, 4.2% on Citeseer, 6.7% on Pubmed, 5.9% on Wiki and 3.4% on
Physics. We also check the standard deviation of the performance of SEComm
over 10 runs in each dataset. Standard deviation lies in the range of 0.5% -
1% on all the datasets, which shows the robustness of SEComm. Among the
baselines, AGC mostly performs better than others. But, AGC is computationally
much expensive because of their adaptive strategy for determining k in k-order
graph convolution. As shown in Table 4, runtime of AGC on Physics dataset is
~16 times more than that of SEComm. As expected, the algorithms which use both
graph structure and node attributes perform better than the ones which use
only one of those. The consistent performance of SEComm on all the datasets
shows the importance of integrating the objective of community detection
directly into the framework of self-supervised graph neural network (the loss
from these components propagates to each other through backpropagation).
Further, use of the principle of self-expressiveness regularizes the
communities formed in SEComm to achieve better performance. Run-time of SEComm
and its various components on all the datasets is shown in Table 5. Usefulness
of the individual components of SEComm are presented in Section 5.7.

Dataset | GNN | Self-Express. | Comm. Module | Total
---|---|---|---|---
Wiki | 5.4s | 34.8s | 44.8s | 85.1s
Cora | 6.6s | 39.2s | 84.2s | 130.1s
Citeseer | 27.8s | 58.4s | 157.8s | 244.1s
Pubmed | 239.6s | 103.6s | 121.9s | 465.1s
Physics | 106.9s | 280.4s | 236.4s | 623.7s

Table 5: Runtime of SEComm on various datasets

###  5.5 Loss and Accuracy of SEComm

Typically for an unsupervised algorithm, the loss that it minimizes and the
metric that is used to evaluate the performance of the algorithm are not
necessarily the same. So, it is important to see if reducing the loss over the
epochs actually increases the performance of the algorithm with respect to the
quality metric. For SEComm, we plot the loss in Equation 10 and the clustering
accuracy that it achieves over different epochs of the algorithm for the
datasets Cora and Wiki in Figure 2. One can see that with the decreasing loss,
overall clustering accuracy improves with some minor fluctuations. Thus, the
unsupervised loss that SEComm minimizes essentially helps to improve the
performance of clustering. This is also another reason of consistent
performance (improved metric scores with less standard deviation) of SEComm on
multiple datasets.

(a) Cora

(b) Wiki

Figure 3: Performance of SEComm with varying embedding dimension

(a) Cora

(b) Pubmed

Figure 4: Performance of SEComm with varying threshold $\theta_{low}$

(a) Pubmed

(b) Pubmed

Figure 5: Performance of SEComm with varying number of batches sampled and
batch size respectively

Methods | Cora | Citeseer | Pubmed | Wiki
---|---|---|---|---
| Acc% | NMI% | F1% | Acc% | NMI% | F1% | Acc% | NMI% | F1% | Acc% | NMI% | F1%
SEComm-GNN | 72.85 | 54.40 | 68.16 | 64.95 | 39.01 | 60.51 | 49.63 | 19.59 | 40.04 | 32.14 | 34.88 | 28.34
SEComm-Spectral | 74.26 | 55.38 | 71.75 | 57.58 | 36.63 | 54.99 | - | - | - | 49.81 | 51.07 | 43.95
SEComm-Embeddings | 73.85 | 57.57 | 65.07 | 60.29 | 35.46 | 54.46 | 68.30 | 34.47 | 67.97 | 35.67 | 36.12 | 29.49
SEComm | 75.92 | 56.04 | 73.94 | 69.82 | 42.53 | 60.25 | 74.49 | 36.50 | 73.50 | 53.10 | 51.38 | 44.48

Table 6: Model Ablation Study of SEComm

###  5.6 Sensitivity to Hyperparameters

In this section, we show the sensitivity of SEComm to different
hyperparameters. We keep all other hyperparameters fixed while changing only
the hyperparameter of interest. We vary the embedding dimension $F^{\prime}$
from 64 to 512 for Cora and Wiki and show the performance of community
detection in Figure 3. We can see some fluctuation in the performance for
Cora. As other hyperparameters are tuned keeping $F^{\prime}=256$ for Cora,
there are some sudden lows around it. For Wiki, the performance was low at
$F^{\prime}=64$ as that is not sufficient enough to preserve all the
information about the graph in the embeddings. It increases at
$F^{\prime}=128$. There is a gradual decrease of performance beyond that as
the embeddings start holding noisy information when dimension increases more.

We check the performance on Cora, Wiki and Pubmed in Figure 4 with varying
$\theta_{low}$ (we set $\theta_{high}=1-\theta_{low}$). As we decrease
$\theta_{low}$, we are filtering out more pair-wise similarities, especially
the ones which lies in the mid zone of the range $[0,1]$. Filtering out such
similarity values might lead to less amount of data to regularize the
communities learned by SEComm in Equation 9 for a smaller dataset. Thus, the
performance on Cora is affected when $\theta_{low}$ is very low in Figure 4.
But on a larger dataset, filtering out such less informative similarity values
(as explained in Section 4.2) can lead to removing noise and help improving
the performance. Thus in Pubmed, better performance is observed around
$\theta_{low}=0.05$ (which implies $\theta_{high}=0.95$). Below which the
amount data becomes too less to train the algorithm properly, and above which
it was adding noise.

As discussed in Section 4.2, we do not need to train the self-expressive layer
(in Eq. 4) on the complete dataset. In Figure 5(a), we vary the number of
batches sampled, where each batch contains 2000 nodes. We can see that the
performance improved initially and then saturates when number of batches is 5
or more. Thus, optimal performance on Pubmed for community detection can be
achieved by using only ~50% (or more) of data points to train the self-
expressive layer. In Figure 5(b), we change the batch size, keeping number of
batches as 6. More is the batch size, more computation resource and time
needed. We observe that SEComm is able to achieve reasonably good performance
when the batch size is 2000 or more.

###  5.7 Model Ablation Study

In this section, we show the usefulness of different components of SEComm. In
particular, we check the community detection performance in the following
scenarios.

SEComm-GNN We run k-means on the node embeddings produced by the self-
supervised GNN used in SEComm, without running the other modules of SEComm.

SEComm-Spectral We run spectral clustering on the complete similarity matrix
$S$ to find node clusters. However, $S$ can be computed only for smaller
graphs and hence this experiment cannot be performed on larger datasets like
Pubmed and Physics.

SEComm-Embeddings We run k-means on the node embeddings generated after the
complete training of SEComm (including the self-expressive and clustering
modules)

We compare the results of the above with the community detection output of the
complete model of SEComm in Table 6. Again we use three metrics clustering
accuracy, NMI and F1 score to evaluate the quality of community detection.
Interestingly, we do not see any clear winners between the three model
variants. But it is clear from the reported performance numbers that the
complete model of SEComm outperforms its variants (except Citeseer-F1%).

##  6 Discussion and Conclusion

In this work, we have proposed a novel graph neural network that can directly
be used for node community detection in a graph. We use the principle of self-
expressiveness to derive a set of soft node-pair constraints to regularize the
formation of the communities. To the best of our understanding, this is the
first work to integrate a self-expressive layer into a self-supervised GNN.
Our approach is highly scalable, without compromising the performance of
community detection. SEComm is able to achieve state-of-the-art performance on
all the datasets that we used for community detection.

Due to the use of graph neural network to directly generate community
memberships of nodes, SEComm can work in an inductive setup. It would be
interesting to analyze the performance of SEComm on newly added nodes or even
new graphs without retraining.

###### Acknowledgements.

We want to thank Prof. M. Narasimha Murty from CSA, IISc for his feedback on
this work.

## References

  * Aggarwal and Reddy [2014] Charu C Aggarwal and Chandan K Reddy.  Data clustering.  _Algorithms and applications. Chapman &Hall/CRC Data mining and Knowledge Discovery series, Londra_, 2014.
  * Bandyopadhyay et al. [2019] Sambaran Bandyopadhyay, N Lokesh, and M Narasimha Murty.  Outlier aware network embedding for attributed networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 12–19, 2019.
  * Bandyopadhyay et al. [2020] Sambaran Bandyopadhyay, N Lokesh, Saley Vishal Vivek, and MN Murty.  Outlier resistant unsupervised deep architectures for attributed network embedding.  In _Proceedings of the 13th International Conference on Web Search and Data Mining_ , pages 25–33, 2020.
  * Belghazi et al. [2018] Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeshwar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm.  Mutual information neural estimation.  In _International Conference on Machine Learning_ , pages 531–540, 2018.
  * Bianchi et al. [2020] Filippo Maria Bianchi, Daniele Grattarola, and Cesare Alippi.  Spectral clustering with graph neural networks for graph pooling.  In _International Conference on Machine Learning_ , pages 874–883. PMLR, 2020.
  * Bo et al. [2020] Deyu Bo, Xiao Wang, Chuan Shi, Meiqi Zhu, Emiao Lu, and Peng Cui.  Structural deep clustering network.  In _Proceedings of The Web Conference 2020_ , pages 1400–1410, 2020\.
  * Cao et al. [2016] Shaosheng Cao, Wei Lu, and Qiongkai Xu.  Deep neural networks for learning graph representations.  In _AAAI_ , volume 16, pages 1145–1152, 2016.
  * Chen et al. [2019] Zhengdao Chen, Lisha Li, and Joan Bruna.  Supervised community detection with line graph neural networks.  In _International Conference on Learning Representations_ , 2019.
  * Desai et al. [2021] Utkarsh Desai, Sambaran Bandyopadhyay, and Srikanth Tamilselvam.  Graph neural network to dilute outliers for refactoring monolith application.  In _Proceedings of 35th AAAI Conference on Artificial Intelligence (AAAI’21)_ , 2021.
  * Elhamifar and Vidal [2009] Ehsan Elhamifar and René Vidal.  Sparse subspace clustering.  In _2009 IEEE Conference on Computer Vision and Pattern Recognition_ , pages 2790–2797. IEEE, 2009.
  * Fortunato and Hric [2016] Santo Fortunato and Darko Hric.  Community detection in networks: A user guide.  _Physics reports_ , 659:1–44, 2016.
  * Gilmer et al. [2017] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , pages 1263–1272, 2017.
  * Hamilton et al. [2017a] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In _Advances in neural information processing systems_ , pages 1024–1034, 2017a.
  * Hamilton et al. [2017b] William L Hamilton, Rex Ying, and Jure Leskovec.  Representation learning on graphs: Methods and applications.  _arXiv preprint arXiv:1709.05584_ , 2017b.
  * Hjelm et al. [2019] R Devon Hjelm, Alex Fedorov, Samuel Lavoie-Marchildon, Karan Grewal, Phil Bachman, Adam Trischler, and Yoshua Bengio.  Learning deep representations by mutual information estimation and maximization.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=Bklr3j0cKX.
  * Ji et al. [2014] Pan Ji, Mathieu Salzmann, and Hongdong Li.  Efficient dense subspace clustering.  In _IEEE Winter Conference on Applications of Computer Vision_ , pages 461–468. IEEE, 2014.
  * Ji et al. [2017] Pan Ji, Tong Zhang, Hongdong Li, Mathieu Salzmann, and Ian Reid.  Deep subspace clustering networks.  In _Advances in Neural Information Processing Systems_ , pages 24–33, 2017.
  * Kipf and Welling [2016] Thomas N Kipf and Max Welling.  Variational graph auto-encoders.  _arXiv preprint arXiv:1611.07308_ , 2016.
  * Kipf and Welling [2017] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In _International Conference on Learning Representations (ICLR)_ , 2017.
  * Li and Vidal [2015] Chun-Guang Li and Rene Vidal.  Structured sparse subspace clustering: A unified optimization framework.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 277–286, 2015.
  * Lu et al. [2012] Can-Yi Lu, Hai Min, Zhong-Qiu Zhao, Lin Zhu, De-Shuang Huang, and Shuicheng Yan.  Robust and efficient subspace segmentation via least squares regression.  In _European conference on computer vision_ , pages 347–360. Springer, 2012.
  * Pan et al. [2018] Shirui Pan, Ruiqi Hu, Guodong Long, Jing Jiang, Lina Yao, and Chengqi Zhang.  Adversarially regularized graph autoencoder for graph embedding.  _arXiv preprint arXiv:1802.04407_ , 2018.
  * Perozzi et al. [2014] Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: Online learning of social representations.  In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 701–710, 2014.
  * Shchur et al. [2018] Oleksandr Shchur, Maximilian Mumme, Aleksandar Bojchevski, and Stephan Günnemann.  Pitfalls of graph neural network evaluation.  _arXiv preprint arXiv:1811.05868_ , 2018.
  * Sun et al. [2019] Fan-Yun Sun, Meng Qu, Jordan Hoffmann, Chin-Wei Huang, and Jian Tang.  vgraph: A generative model for joint community detection and node representation learning.  In _Advances in Neural Information Processing Systems_ , pages 514–524, 2019.
  * Sun et al. [2020] Fan-Yun Sun, Jordan Hoffman, Vikas Verma, and Jian Tang.  Infograph: Unsupervised and semi-supervised graph-level representation learning via mutual information maximization.  In _International Conference on Learning Representations_ , 2020.  URL https://openreview.net/forum?id=r1lfF2NYvH.
  * Veličković et al. [2019] Petar Veličković, William Fedus, William L Hamilton, Pietro Liò, Yoshua Bengio, and R Devon Hjelm.  Deep graph infomax.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=rklz9iAcKQ.
  * Wagstaff et al. [2001] Kiri Wagstaff, Claire Cardie, Seth Rogers, Stefan Schrödl, et al.  Constrained k-means clustering with background knowledge.  In _Icml_ , volume 1, pages 577–584, 2001.
  * Wang et al. [2017] Chun Wang, Shirui Pan, Guodong Long, Xingquan Zhu, and Jing Jiang.  Mgae: Marginalized graph autoencoder for graph clustering.  In _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management_ , pages 889–898, 2017.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  _IEEE Transactions on Neural Networks and Learning Systems_ , 2020\.
  * Xie et al. [2013] Jierui Xie, Stephen Kelley, and Boleslaw K Szymanski.  Overlapping community detection in networks: The state-of-the-art and comparative study.  _Acm computing surveys (csur)_ , 45(4):1–35, 2013\.
  * Zhang et al. [2020] Binbin Zhang, Zhizhi Yu, and Weixiong Zhang.  Community-centric graph convolutional network for unsupervised community detection.  In _IJCAI_ , 2020.
  * Zhang et al. [2019a] Junjian Zhang, Chun-Guang Li, Chong You, Xianbiao Qi, Honggang Zhang, Jun Guo, and Zhouchen Lin.  Self-supervised convolutional subspace clustering network.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 5473–5482, 2019a.
  * Zhang and Chen [2018] Muhan Zhang and Yixin Chen.  Link prediction based on graph neural networks.  In _Advances in Neural Information Processing Systems_ , pages 5165–5175, 2018.
  * Zhang et al. [2019b] Xiaotong Zhang, Han Liu, Qimai Li, and Xiao-Ming Wu.  Attributed graph clustering via adaptive graph convolution.  In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , pages 4327–4333. AAAI Press, 2019b.
  * Zhu et al. [2020] Yanqiao Zhu, Yichen Xu, Feng Yu, Qiang Liu, Shu Wu, and Liang Wang.  Deep graph contrastive representation learning.  _arXiv preprint arXiv:2006.04131_ , 2020.