# A Reduced Model for a Phoretic Swimmer

Alexander Farutin<EMAIL_ADDRESS>Univ. Grenoble Alpes,
CNRS, LIPhy, F-38000 Grenoble, France  Suhail M. Rizvi  Univ. Grenoble Alpes,
CNRS, LIPhy, F-38000 Grenoble, France  Current address: Department of
Biomedical Engineering, Indian Institute of Technology Hyderabad, Sangareddy,
Telangana 502285, India  Wei-Fan Hu  Department of Mathematics, National
Central University, 300 Zhongda Road, Taoyuan 320, Taiwan  Te-Sheng Lin
Department of Applied Mathematics, National Chiao Tung University, 1001 Ta
Hsueh Road, Hsinchu 300, Taiwan  Salima Rafai  Univ. Grenoble Alpes, CNRS,
LIPhy, F-38000 Grenoble, France  Chaouqi Misbah  Univ. Grenoble Alpes, CNRS,
LIPhy, F-38000 Grenoble, France

###### Abstract

We consider a 2D model of an autophoretic particle in which the particle has a
circular shape and emits/absorbs a solute that diffuses and is advected by the
suspending fluid. Beyond a certain emission/absorption rate (characterized by
a dimensionless Péclet number, $Pe$) the particle is known to undergo a
bifurcation from a non motile to a motile state, with different trajectories,
going from a straight to circular and to a chaotic motion by progressively
increasing $Pe$. From the full model involving solute diffusion and advection,
we derive a reduced closed model which involves only two time-dependent
amplitudes $C_{1}(t)$ and $C_{2}(t)$ corresponding to the first two Fourier
modes of the solute concentration field. This model consists of two coupled
nonlinear ordinary differential equations for $C_{1}$ and $C_{2}$ and presents
several great advantages:(i) the straight and circular motions can be handled
fully analytically, (ii) complex motions such as chaos can be analyzed
numerically very efficiently in comparison to the numerically expensive full
model involving partial differential equations, (iii) the reduced model has a
universal form dictated only by symmetries, (iv) the model can be extended to
higher Fourier modes. The derivation method is exemplified for a 2D model, for
simplicity, but can easily be extended to 3D, not only for the presently
selected phoretic model, but also for any model in which chemical activity
triggers locomotion. A typical example can be found, for example, in the field
of cell motility involving acto-myosin kinetics. This strategy offers an
interesting way to cope with swimmers on the basis of ordinary differential
equations, allowing for analytical tractability and efficient numerical
treatment.

##  I Introduction

Phoretic particles (rigid particles or drops) have been attracting recently an
increased interest theoretically, numerically and experimentally Izri  _et
al._ (2014); Michelin, Lauga, and Bartolo (2013); Michelin and Lauga (2014);
Schmitt and Stark (2013); Jin, Krüger, and Maass (2017); Morozov and Michelin
(2019a, b); Hu  _et al._ (2019); Morozov (2020); Chen  _et al._ (2021);
Hokmabad  _et al._ (2021); Izzet  _et al._ (2020). In its simplest version,
the model consists of a particle that emits or absorbs (with an
emission/absorption rate $A$) a solute which diffuses, with bulk diffusion
constant $D$, and is advected by the suspending fluid. The interaction between
the solute and the particle can be shown Michelin and Lauga (2014) to result
in a tangential flow along the particle, in the form of (in the frame moving
with particle) ${\bf u}\sim M\nabla_{s}c$, where $\bf u$ is the velocity, $c$
the solute concentration, $\nabla_{s}$ is the gradient along the particle
surface, and $M$ is a mobility factor involving fluid viscosity and the
interaction potential between the solute and the particle. This problem can be
characterized by a Péclet number $Pe=AMa/D^{2}$ where $a$ is the particle
radius. The solute $c({\bf r},t)$ (${\bf r}$ is a position in space) and
velocity ${\bf u}({\bf r},t)$ obey the advection-diffusion and the Stokes
equations. It has been shown analytically, from a linear stability analysis,
both in 2D Hu  _et al._ (2019) and 3D Michelin, Lauga, and Bartolo (2013) that
if $Pe<Pe_{1}$ (where $Pe_{1}$ is a critical number), the particle undergoes a
bifurcation from a non-motile state (swimming velocity $u_{0}=0$) into a
motile state ($u_{0}\neq 0$). A linear stability only informs us on the
instability onset from one state to another, but is not sufficient to describe
how the velocity behaves with $Pe$, where a nonlinear analysis is needed.
Numerical simulations showedMichelin, Lauga, and Bartolo (2013); Hu  _et al._
(2019) that the behavior of $u_{0}$ is well represented by ${u}_{0}\sim\pm(Pe-
Pe_{1})^{1/2}$ in the vicinity of the bifurcation point; the solution
$u_{0}=0$ always exists, and is stable for $Pe<Pe_{1}$ and becomes unstable
for $Pe>Pe_{1}$, in favor of two stable branches of solutions
${u}_{0}\sim\pm(Pe-Pe_{1})^{1/2}$. This is a classical pitchfork bifurcation.
Numerical simulations in 2D Hu  _et al._ (2019) also showed that by increasing
$Pe$ further the straight moving solution becomes unstable in favor of various
states (meandering, circular, and chaotic solutions). Simulations in 3D under
axisymmetric constraint (imposing to the particle to move along a line) also
reported Morozov and Michelin (2019a)on chaotic solutions (the particle goes
back and forth in a chaotic manner). Relaxing the axisymmetric constraint in
3D also revealed, among other motion, meandering, and chaotic motion in the
form of persistent random walk Hu  _et al._ (2021). The transition from a non
motile to a motile state, as well as irregular motion (apparently chaotic) has
also been reported experimentallyHokmabad  _et al._ (2021); Jin, Krüger, and
Maass (2017); Izri  _et al._ (2014); Izzet  _et al._ (2020), and in some
numerical simulations Chen  _et al._ (2021). The transition from a non motile
to a motile state also takes place in the presence of a Marangoni stress Izri
_et al._ (2014); Morozov and Michelin (2019a). In a completely different
context, that of cell motility driven by acto-myosin rich dynamics (transition
from a motile to a non motile state, Hopf bifurcation…) have also been
reported when the activity of cell (represented by myosin contractility)
exceeds a certain critical value Hawkins  _et al._ (2011); Callan-Jones  _et
al._ (2016); Farutin  _et al._ (2019).

The above non trivial behaviors have been obtained in most cases by numerical
simulations of model equations involving diffusion, and advection with a
moving boundary. It is thus highly desirable to see wether or not these
features can be translated into a more universal language. One way is to seek
for a reductive perturbative scheme that allows to extract from the full model
prototypical simplified equations which are universal (in the sense that their
form does not depend on the considered swimmer model). This is the main goal
of this paper. Our main objective here is, by starting from the above
described phoretic model, to extract, by following a nonlinear perturbative
scheme, two coupled (weakly nonlinear) equations for two complex amplitudes of
the first two Fourier modes, which depend only on time, while the spatial
dependence can be handled analytically. These equations will enable us to
find, besides the first (primary) bifurcation (transition from a motile to a
nonmotile state) a secondary bifurcation in the form of a circular trajectory.
We will also exemplify its power in capturing irregular motion. We will
exemplify the method for a 2D model (for simplicity), but the technique should
work perfectly well in 3D, not only for the phoretic model, but also for any
model where locomotion is assisted by chemical activity.

##  II Problem formulation

Our starting point is a simple phoretic model in 2D. We first recall the
problem formulationMichelin, Lauga, and Bartolo (2013); Michelin and Lauga
(2014). We consider a concentration field $c(\bm{r},t)$ subject to an
advection-diffusion equation

$\partial_{t}c(\bm{r},t)=D\nabla^{2}c(\bm{r},t)-\bm{\nabla}\bm{\cdot}[c(\bm{r},t)\bm{u}(\bm{r},t)],$
(1)

where the advection velocity field $\bm{u}(\bm{r},t)$ in the co-moving frame
is obtained by solving the Stokes equations in the fluid outside the particle.
Here $t$ is the time, $\bm{r}$ is the 2D position vector, measured relative to
the particle center, and $D$ is the diffusion coefficient. We use polar
coordinates $(r,\phi)$, centered at the particle position in the following
derivation. The boundary conditions for the concentration field read

$\partial_{r}c(\bm{r},t)|_{r=a}=-A/D,\,\,\,c(\bm{r},t)|_{r=R}=0,$ (2)

where $a$ is the particle size and $R$ is the system size. The flow velocity
satisfies the following boundary conditions:

$\bm{u}(\bm{r},t)|_{r=\infty}=-\bm{u}_{0}(t),\,\,\,\bm{u}(\bm{r},t)|_{r=a}={M\over
a}\bm{\nabla}^{s}c(\bm{r},t)|_{r=a},$ (3)

where $\bm{u}_{0}(t)$ is the swimming speed of the particle, $\bm{\nabla}^{s}$
is the surface gradient operator, and $M$ is the particle mobility. The system
is closed by requiring the net force acting on the particle to be zero, which
allows one to compute the swimming velocity $u_{0}(t)$. It is conventional to
express the fluid velocity in terms of the stream function $\psi$

$u_{r}(\bm{r},t)=\frac{\partial_{\phi}\psi(\bm{r},t)}{r},\,\,\,u_{\phi}(\bm{r},t)=-\partial_{r}\psi(\bm{r},t).$
(4)

The stream function is related to the concentration distribution along the
particle boundary by the boundary condition (3), as explained below.

Hereafter, we use a non-dimensional form of the problem, in which the size of
the particle $a$, the diffusion coefficient $D$, and the release rate $A$ are
all set to 1. The remaining two non-dimensional numbers are $Pe=AMa/D^{2}$ and
system size $R/a$. We will keep notation $R$ for this dimensionless quantity.

##  III Expansion in Fourier harmonics

Following the previous works Michelin, Lauga, and Bartolo (2013); Hu  _et al._
(2019), the concentration field and the stream function are expanded in
Fourier harmonics of the polar angle $\phi$.

$c(\bm{r},t)=\sum\limits_{l=-\infty}^{\infty}c_{l}(r,t)e^{il\phi},\,\,\,\psi(\bm{r},t)=\sum\limits_{l=-\infty}^{\infty}\psi_{l}(r,t)e^{il\phi},$
(5)

where we have $c_{l}(r)=c_{-l}(r)^{*}$ and $\psi_{l}(r)=\psi_{-l}(r)^{*}$ as
requirement for $c$ and $\psi$ to be real. The Stokes equations and the
boundary condition (3) allow one to express the amplitudes $\psi_{l}$ as

$\psi_{l}(r)=\frac{ilPe(1-r^{2})}{2r^{|l|}}c_{l}(1).$ (6)

Exploiting the force-free condition on the particle, the swimming velocity is
obtained as Hu  _et al._ (2019) $v_{x}=-Pe\Re c_{1}(1)$, $v_{y}=Pe\Im
c_{1}(1)$.

The advection-diffusion equation (1) is expanded in Fourier harmonics as

$\partial_{t}c_{l}(r)=\hat{D}_{l}c_{l}(r)-[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l},$
(7)

where $\hat{D}_{l}$ is the diffusion operator applied to the $l$-th Fourier
harmonic:

$\hat{D}_{l}=\frac{1}{r}\partial_{r}(r\partial_{r})-\frac{l^{2}}{r^{2}}$ (8)

and $[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}$ is the l-th Fourier
harmonic of the advection term.

##  IV Stationary state and linear approximation

The linear stability of the isotropic solution was already analyzed in
previous works in 3D Michelin, Lauga, and Bartolo (2013) and in 2D Hu  _et
al._ (2019). Nevertheless, we present it briefly as the first step of the
solution of the non-linear problem. The isotropic solution corresponds to

$c_{0}(r)=-\ln\frac{r}{R}$ (9)

Substituting eq. (5) with $c_{0}(r)$ given by eq. (9) into eq. (7) and
neglecting terms non-linear in $c_{l}$, yields the following evolution
equation for $c_{l}$

$\partial_{t}c_{l}(r)=\hat{L}_{l}(Pe)c_{l}(r)=\hat{D}_{l}c_{l}(r)+Peu_{l}(r)c_{l}(1),$
(10)

where $\hat{L}_{l}$ is the linear stability operator for the $l$-th harmonic.
Here $Peu_{l}e^{il\phi}$ is the advection of the unperturbed solution (9) with
velocity field driven by the surface concentration $c_{l}(1)e^{il\phi}$:

$u_{l}(r)=\frac{l^{2}(r^{2}-1)}{2r^{l+2}}.$ (11)

The linear stability of the isotropic solution is governed by the eigenvalues
$\lambda_{l,k}$ of the operators $\hat{L}_{l}$, defined as

$\hat{L}_{l}f_{l,k}(r)=\lambda_{l,k}f_{l,k}(r),$ (12)

where $f_{l,k}$ are the corresponding proper functions. The subscript $l$ is
associated with the $\phi$ variable and $k$ with the $r$ one. We index the
eigenvalues for given $l$ in descending order with respect to their real part,
starting with $\lambda_{l,0}\equiv\lambda_{l}$, which corresponds to the most
unstable mode. It is essential for our analysis to assume the eigenvalue
spectrum to be discrete, which is the case for a finite domain. In general,
neither the eigenvalues nor the proper functions of $\hat{L}_{l}$ have an
elementary expression. However, it is known that for each $l\geq 1$, there is
exactly one critical Peclet number $Pe_{l}$ such that the operator
$\hat{L}_{l}(Pe_{l})$ has an eigenvalue equal to zeroMichelin and Lauga
(2014); Hu  _et al._ (2019). Both $Pe_{l}$ (i.e. $Pe_{1}$ and $Pe_{2}$) as a
function of $R$ and the proper function of $\hat{L}_{l}(Pe_{l})$ corresponding
to the zero eigenvalue have an explicit expression in elementary functions,
given below.

For low enough $Pe$, the eigenvalues of $\hat{L}_{l}$ are close to those of
the diffusion operator $\hat{D}_{l}$ and are thus all negative. The isotropic
solution is stable in this case (growth rate of all modes is negative for all
$l$). As $Pe$ increases, one of the proper values for given $l$
($\lambda_{l,0}\equiv\lambda_{l}$) becomes equal to zero at $Pe=Pe_{l}$ and
positive for $Pe>Pe_{l}$. The corresponding proper function $f_{l,0}(r)\equiv
f_{l}(r)$ defines the perturbation mode that becomes unstable at $Pe=Pe_{l}$.
The linear stability analysis allows us to determine the angular dependence
(given by $l$ corresponding to the lowest $Pe_{l}$) and the radial dependence
(given by $f_{l}(r)$) of the solution in the anisotropic phase. Further
analysis is performed by higher-order perturbation expansion in order to
compute the swimming speed and the angular velocity of the particle in the
anisotropic phase, as explained below.

Here we exclude the possibility of transition to an anisotropic concentration
distribution by a Hopf bifurcation. Since the operators $\hat{L}_{l}$ are not
self-adjoint, their proper values can be complex. It is therefore possible for
the isotropic solution to become unstable due to one of the complex
eigenvalues seeing its real part to become positive. We exclude this
possibility based on the results of full numerical simulations, shown below.

##  V Weakly non-linear expansion

This section presents the derivation of the weakly non-linear equations
governing the particle dynamics. We consider the situation where both the
first and the second harmonic are close to the instability, or, equivalently,
$Pe$ is simultaneously close to $Pe_{1}$ and $Pe_{2}$. This allows one to keep
only thee first two leading harmonics. But It will appear that extension to
other higher harmonics is straightforward, albeit we restrict ourselves here
to keeping only the two first harmonics. This is made possible by choosing an
appropriate size of the system $R$. We find that $Pe_{1}$ becomes equal to
$Pe_{2}$ at $R=R_{c}$, where $R_{c}=3.17493$. We thus have to assume that
$R=R_{c}+O(\epsilon)$, but also $R>R_{c}$, the second condition guaranteeing
$Pe_{1}<Pe_{2}$. Here $\epsilon$ is the small parameter with respect to which
the perturbation expansion is made. With these assumptions, $c=c_{0}(r)$ is
the only possible solution for $Pe<Pe_{1}$. This solution becomes unstable at
$Pe=Pe_{1}$ and a straight motion is expected to emerge at this point. At some
value of $Pe$ between $Pe_{1}$ and $Pe_{2}$, the straight motion is expected
to become unstable in favor of a circular motionHu  _et al._ (2019). We have
been able to compute the $O(\epsilon^{2})$ terms analytically for an arbitrary
value of $R$ and also all $O(\epsilon^{3})$ terms for a given numerical value
of $R$. This is sufficient to reproduce quantitatively the straight and
circular motions, as shown below. Our reduced model accounts also for
irregular motions, reported inHu  _et al._ (2019).

###  V.1 General strategy

The goal of this derivation is to reduce the full dynamics of the autophoretic
particle to a simplified system of ordinary differential equations for two
complex variables, $C_{1}$ and $C_{2}$, where $C_{1}$ is the amplitude of the
first harmonic and $C_{2}$ is the amplitude of the second harmonic in the
concentration field. The main challenge is that in general, each Fourier
harmonic $c_{l}$ of the concentration field is a function of the distance from
the center of the particle $r$. There is thus no straightforward way to
represent the whole function $c_{l}$ by a single scalar variable. We overcome
this problem by decomposing the functions $c_{1}(r)$ and $c_{2}(r)$ as

$c_{l}(r,t)=C_{l}(t)f_{l}(r)+\delta c_{l}(r,t)\,\,\,l\in\\{1,2\\},$ (13)

where $C_{l}(t)$ is the complex amplitude, $f_{l}(r)$ is the proper function
such that $\hat{L}_{l}(Pe_{l})f_{l}(r)=0$, and $\delta c_{l}(r,t)$ is a
projection of the function $c_{l}(r,t)$ on the space of all other proper
functions of the operator $\hat{L}_{l}(Pe_{l})$:

$\delta c_{l}(r,t)=\sum\limits_{k>0}C_{l,k}(t)f_{l,k}(r).$ (14)

The decomposition (13) is made using the proper functions of
$\hat{L}_{l}(Pe_{l})$ even for $Pe\neq Pe_{l}$. This is to capitalize on the
explicit expression for $Pe_{l}$ and $f_{l}$ for $\hat{L}_{l}(Pe_{l})$. The
difference between the linear stability operator $\hat{L}_{l}(Pe)$ and
$\hat{L}_{l}(Pe_{l})$ is of order $O(\epsilon)$ and is duly accounted for as a
higher-order correction during the derivation.

The eigenvalues corresponding to the functions $f_{l,k}$ with $k>0$ all have
negative real parts. Therefore, the relaxation time scale of the amplitudes
$C_{l,k}(t)$ is defined by the absolute values of $\lambda_{l,k}$. Since the
set of $\lambda_{l,k}$ is discrete and we have assumed $\lambda_{l,k}$ to be
ordered decreasingly, the longest relaxation time scale is defined by
$|\lambda_{l,1}|$. The relaxation time of the $f_{l}(r)$ mode scales as
$1/|Pe-Pe_{l}|$ for $Pe$ close to $Pe_{l}$, which determines the main time
scale of the dynamics close to the critical point. We thus have been able to
reexpress the functions $\delta c_{l}(r,t)$ as a non-linear function of
$C_{1}$ and $C_{2}$, with $r$-dependent coefficients. This is done by
adiabatic elimination, as explained below. Since the eigenvalues and the
proper functions of the operator $\hat{L}_{l}(Pe_{l})$ do not have an
elementary expression, the following procedure relies on the adiabatic
elimination of the amplitude $\delta c_{l}(r,t)$ as a function, instead of
eliminating each amplitude $C_{l,k}$ separately.

A similar strategy is employed to compute the $r$ dependence of the amplitude
of the other Fourier harmonics $c_{l}(r,t)$ for $l\not\in\\{1,2\\}$ as a
function of $C_{1}$ and $C_{2}$, as also explained below. This effectively
shows that the whole concentration field $c(r,\phi,t)$ can be reexpressed as a
perturbation expansion in powers of $C_{1}$ and $C_{2}$ after an initial
transient relaxation. We present below a procedure to compute the coefficients
of this expansion as explicit functions of $r$ and $\phi$.

###  V.2 Perturbation expansion

The following derivation is based on the assumption of $C_{1}$ and $C_{2}$
being small. It is usual for a pitchfork bifurcation to show the perturbation
amplitude scaling as $(Pe-Pe_{1})^{1/2}$. Here, however, we have two coupled
harmonics that are close to instability, which results in scaling
$C_{1}=O(\epsilon)$ and $C_{2}=O(\epsilon)$. The consistency of the evolution
equations requires us to admit the following assumptions: $|Pe-
Pe_{1}|=O(\epsilon)$, $|Pe-Pe_{2}|=O(\epsilon)$,
$\partial_{t}C_{1}=O(\epsilon^{2})$, $\partial_{t}C_{2}=O(\epsilon^{2})$,
$\delta c_{1}(r,t)=O(\epsilon^{2})$, $\delta c_{2}(r,t)=O(\epsilon^{2})$,
$\partial_{t}\delta c_{1}(r,t)=O(\epsilon^{3})$, $\partial_{t}\delta
c_{2}(r,t)=O(\epsilon^{3})$, $\delta c_{0}(r,t)=O(\epsilon^{2})$,
$\partial_{t}\delta c_{0}(r,t)=O(\epsilon^{3})$. For $l>2$, we have $\delta
c_{l}(r,t)=O(\epsilon^{\lceil l/2\rceil})$ and $\partial_{t}\delta
c_{l}(r,t)=O(\epsilon^{\lceil l/2\rceil+1})$. The smallness of $\partial_{t}$
is due to the critical slowing down at the bifurcation point.

###  V.3 Adiabatic elimination

The problem requires us to reduce the partial derivative equations for
functions $c_{l}(r)$ to ordinary differential equations of two scalar
amplitudes $C_{1}$ and $C_{2}$. This is performed by applying adiabatic
elimination, as explained in this Section.

First, we introduce the right proper function $g_{l}$ of the operator
$\hat{L}_{l}(Pe_{l})$ such that the corresponding eigenvalue is zero:

$\hat{L}_{l}^{+}(Pe_{l})g_{l}(r)=0,$ (15)

where $\hat{L}_{l}(Pe_{l})^{+}$ is the adjoint operator. The adjoint operator
is defined with respect to the inner product

$\langle f,g\rangle=\int\limits_{1}^{R}f(r)g(r)^{*}rdr,$ (16)

which is chosen to maintain the self-adjoint property of the diffusion
operators $\hat{D}_{l}$ subject to the boundary conditions of the functions
$c_{l}(r)$. The projection condition then reduces to

$\langle g_{l},\delta c_{l}\rangle=0.$ (17)

The concentration evolution is split as

$\partial_{t}c_{l}(r,t)=\partial_{t}C_{l}f_{l}(r)+\partial_{t}\delta
c_{l}(r,t)\textrm{ for }l\in\\{1,2\\}$ (18)

by the definition (13). We use $g_{l}$ to isolate the $C_{l}$ expression from
$\partial_{t}c_{l}(r,t)$:

$\partial_{t}C_{l}(t)=\frac{\langle
g_{l},\partial_{t}c_{l}(r,t)\rangle}{\langle g_{l},f_{l}\rangle},$ (19)

where $\partial_{t}c(r,t)$ is computed according to eq. (7).

Equation (7) for $l\in\\{1,2\\}$ can be rewritten as

$\begin{split}\partial_{t}c_{l}(r,t)=&\hat{L}_{l}(Pe_{l})c_{l}(r,t)+(Pe-
Pe_{l})u_{l}(r)c_{l}(1)-[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}-Peu_{l}(r)c_{l}(1)\\\
=&\hat{L}_{l}(Pe_{l})\underbrace{\delta
c_{l}(r,t)}_{O(\epsilon^{2})}+\underbrace{(Pe-
Pe_{l})}_{O(\epsilon)}u_{l}(r)\underbrace{c_{l}(1)}_{O(\epsilon)}\underbrace{-[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}-Peu_{l}(r)c_{l}(1)}_{O(\epsilon^{2})},\end{split}$
(20)

where the last term is $O(\epsilon^{2})$ because $Peu_{l}(r)c_{l}(1)$
represents the $O(\epsilon)$ part of
$[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}$, as defined in the
linearization procedure. Combining eqs. (18) and (20), we obtain

$f_{l}(r)\partial_{t}C_{l}(t)=\hat{L}_{l}(Pe_{l})\delta
c_{l}(r,t)+q_{l}(r,t)\textrm{ for }l\in\\{1,2\\},$ (21)

where

$q_{l}(r,t)=(Pe-
Pe_{l})u_{l}(r)c_{l}(1)-[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}-Peu_{l}(r)c_{l}(1)-\partial_{t}\delta
c_{l}(r,t)\textrm{ for }l\in\\{1,2\\}.$ (22)

Equation (21) is the main equation of the derivation, which we have solved
together with the orthogonality condition (17) to express both
$\partial_{t}C_{l}$ and $\delta c(r,t)$ as a function of $q_{l}$. The
definition of $q_{l}$ in (22) includes the time derivative of $\delta
c_{l}(r)$. This is not a problem for the subsequent derivation because
$\partial_{t}\delta c_{l}(r)=O(\epsilon^{3})$ and therefore this contribution
can be neglected to compute the $O(\epsilon^{2})$ terms of $\delta
c_{l}(r,t)$. The present work only uses the $O(\epsilon^{2})$ terms of $\delta
c_{l}(r,t)$ but the proposed method is general enough to compute the higher-
order terms as well. Indeed, once the $O(\epsilon^{2})$ terms of $\delta
c_{l}(r,t)$ are obtained as a function of $C_{1}$ and $C_{2}$ with
coefficients that depend on $r$, the dependence of $\delta c_{l}(r,t)$ on time
is contained only in $C_{1}(t)$ and $C_{2}(t)$. It is thus possible to express
the $O(\epsilon^{3})$ terms in $\partial_{t}\delta c_{l}(r)$ as a function of
$C_{1}$, $C_{2}$, and their time derivatives. Those time derivatives can be
further reexpressed through $C_{1}$ and $C_{2}$ by using the final evolution
equations.

The remaining step is to solve the equation (21). The expression of
$\partial_{t}C_{l}$ is given by eq. (19). The main challenge is to extract the
quasi-static amplitudes $\delta c_{l}(r,t)$. Since we have no explicit proper
functions for $\hat{L}_{l}(Pe_{l})$, finding its inverse is far from trivial.
Here the situation is even more complicated because the operator
$\hat{L}_{l}(Pe_{l})$ is not invertible in the first place.

Luckily, here the operator $\hat{L}_{l}(Pe_{l})$ can be written as

$\hat{L}_{l}(Pe_{l})=\hat{D}_{l}+\bm{u}_{l}\otimes\bm{v}_{l},$ (23)

where $\bm{v}_{l}$ is a distribution defined by the relation
$\langle\bm{f},\bm{v}_{l}\rangle=Pe_{l}f(1)$ (aka Dirac delta function). It is
classically known that it is possible to invert explicitly the operators of
form (23), provided the inverse operator for $D_{l}$ is known:

$(\hat{D}_{l}+\bm{u}_{l}\otimes\bm{v}_{l})^{-1}=\hat{D}_{l}^{-1}-\frac{(\hat{D}_{l}^{-1}\bm{u}_{l})\otimes(\bm{v}_{l}\hat{D}_{l}^{-1})}{1+\langle
v_{l},\hat{D}_{l}^{-1}\bm{u}_{l}\rangle}$ (24)

This is the case for our study as the inverse of the diffusion operator can be
represented in an integral form with a simple kernel, as shown in Appendix A.
The remaining difficulty is the lack of a well-defined inverse for the
operator $\hat{L}_{l}(Pe_{l})$, which we overcome, as explained in Appendix B.
The resulting solution of eq. (21) reads:

$\delta c_{l}=-\hat{D}_{l}^{-1}q_{l}+\frac{\langle
v_{l},\hat{D}_{l}^{-1}q_{l}\rangle}{\langle
v_{l},\hat{D}_{l}^{-1}f_{l}\rangle}\hat{D}_{l}^{-1}f_{l}+p_{l}f_{l},$ (25)

where the coefficient $p_{l}$ is chosen to satisfy eq. (17):

$p_{l}=\frac{\langle v_{l},D^{-2}q_{l}\rangle\langle
v_{l},D_{l}^{-1}f_{l}\rangle-\langle v_{l},D_{l}^{-1}q_{l}\rangle\langle
v_{l},D_{l}^{-2}f_{l}\rangle}{\langle v_{l},D_{l}^{-1}f_{l}\rangle^{2}}$ (26)

Note that we have used here that $f_{l}=D_{l}^{-1}u_{l}$ and
$g_{l}=D^{-1}v_{l}$, as shown in the Appendix B.

The functions $\delta c_{l}(r,t)$ for $l\in\\{1,2\\}$ satisfy equation (20),
which can be solved according to eqs. (25) and (26). The remaining functions
$\delta c_{l}(r,t)$ for $l\not\in\\{1,2\\}$ satisfy similar equations

$0=\hat{L}_{l}(Pe)\delta c_{l}(r,t)+q_{l}(r,t)\textrm{ for
}l\not\in\\{1,2\\},$ (27)

where the function $q_{l}$ is given by

$q_{l}(r,t)=-[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}-Peu_{l}(r)c_{l}(1)-\partial_{t}\delta
c_{l}(r,t)\textrm{ for }l\not\in\\{1,2\\}.$ (28)

Since the kernel $L_{l}(Pe_{l})$ is not singular for $l\not\in\\{1,2\\}$ and
$Pe$ close to $Pe_{1}$ and $Pe_{2}$, eq. (24) is used to compute the functions
$\delta c_{l}(r,t)$ as the solutions of eq. (27) for all $l\not\in\\{1,2\\}$.

This completes the proof that the full partial differential equations
governing the evolution of the concentration field $c(r,t)$ can be reduced to
a system of two differential equations for $C_{1}$ and $C_{2}$. The practical
implementation of the derivation procedure and the explicit expressions for
some intermediate results and the final equations are given in the next
Section.

##  VI The derivation procedure

We use an iterative procedure obtaining the expansions of $\partial_{t}C_{1}$,
$\partial_{t}C_{2}$, and $\delta c_{l}$ one order of $\epsilon$ at a time. The
first step is to neglect the $O(\epsilon^{2})$ terms in $c(r,\phi)$, which
allows us to compute the $O(\epsilon^{2})$ terms in $\partial_{t}C_{1}$,
$\partial_{t}C_{2}$, and $\delta c_{l}$. We then substitute the obtained
expressions of $\delta c_{l}$ into $c(r,\phi)$, keeping the terms up to
$O(\epsilon^{2})$ order this time, which is sufficient to obtain the order of
$O(\epsilon^{3})$ terms in $\partial_{t}C_{1}$, $\partial_{t}C_{2}$, and
$\delta c_{l}$. This procedure is continued ad infinitum but we stop at
computing the $O(\epsilon^{3})$ terms of $\partial_{t}C_{1}$ and
$\partial_{t}C_{2}$, which is sufficient for our purposes.

According to the eq. (13) and the order of magnitude analysis (see section
V.2), the concentration field can be written as

$\begin{split}c(r,\phi,t)=&c_{0}(r)+\left[\left(C_{1}e^{i\phi}+C_{1}^{*}e^{-i\phi}\right)f_{1}(r)+\left(C_{2}e^{2i\phi}+C_{2}^{*}e^{-2i\phi}\right)f_{2}(r)\right]+\\\
&\left[\delta c_{0}(r)+\sum\limits_{l=1}^{4}\left(\delta
c_{l}(r)e^{il\phi}+\delta
c_{l}(r)^{*}e^{-il\phi}\right)\right]+O(\epsilon^{3})\end{split}$ (29)

where $C_{1}$ and $C_{2}$ are time-dependent complex amplitudes of the two
modes that are close to instability and $\delta c_{k}(r)$ are functions of $r$
with coefficients that depend on $C_{1}$, $C_{2}$, $R$, and $Pe$.

Each iteration consists in the following steps:

  1. 1.

The concentration field (29) defines the stream function according to eq. (6),
from which the fluid velocity field is calculated.

  2. 2.

Knowing the velocity field and the concentration, the advection term
$\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))$ is computed.

  3. 3.

The advection term is then decomposed into Fourier harmonics which yields the
amplitudes $[\bm{\nabla}\bm{\cdot}(c(\bm{r})\bm{u}(\bm{r}))]_{l}$.

  4. 4.

These amplitudes are then used to compute the $q_{l}$ terms according to eqs.
(22) and (28).

  5. 5.

$\partial_{t}C_{1}$ and $\partial_{t}C_{2}$ are computed at a given order
according to eqs. (19).

  6. 6.

The quasi-static values of $\delta c_{l}$ are computed according to eqs. (25)
and (26) for $l\in\\{1,2\\}$ and according to eq. (24) for $l\in\\{1,2\\}$.

The last step can be omitted for the final iteration.

##  VII Evolution equations for $C_{1}(t)$ and $C_{2}(t)$ to leading order

We first start with the derivation to order $\epsilon^{2}$. The starting point
is to use (21). Multiplying this equation by $g_{l}$ and integrating both
sides according to the scalar product (16) we obtain

$\partial_{t}C_{l}(t)={\langle g_{l},q_{l}\rangle\over\langle
g_{l},f_{l}\rangle},$ (30)

where we have used $\langle g_{l},\hat{L}_{l}\delta
c_{l}\rangle=\langle\hat{L}_{l}^{+}g_{l}(r),\delta c_{l}\rangle=0$, by virtue
of (15). Since $\partial_{t}\delta c_{l}(r,t)=O(\epsilon^{3})$, this term in
$q_{l}$ (Eq. (22)) does not enter to order $\epsilon^{2}$. The next step is to
insert (29) into $q_{l}$ (Eq. (22)) and report the resulting expression into
(30). Note that $\delta c_{l}$ in (29) does not enter to this order neither
since it produces $O(\epsilon^{3})$ contribution to (30) (see also (21) for
orders in $\epsilon$). We are thus left, on the right hand side of (30), with
scalar products involving $f_{l}$ and $g_{l}$ (and $u_{l}$, which is a known
function, see (11) ) with prefactors containing linear and quadratic terms of
$C_{l}(t)^{\prime}s$. The inner product calculation requires, in principle,
the knowledge of $f_{l}$ and $g_{l}$.

The linear stability of the solution (9) is governed by the eigenvalues of the
linear operators $\hat{L}_{l}$ in eq. (10). The critical Peclet numbers are

$Pe_{1}=-\frac{2\left(R^{2}+1\right)}{R^{2}-\left(R^{2}+1\right)\ln{\left(R\right)}-1}$
(31)
$Pe_{2}=-\frac{R^{4}+1}{-\frac{R^{4}}{4}+R^{2}-\ln{\left(R\right)}-\frac{3}{4}}$
(32)

and the corresponding proper functions are

$f_{1}(r)=\frac{R^{2}-r^{2}}{2r\left(R^{2}+1\right)}+\frac{\left(r^{2}+1\right)\ln{\left(\frac{r}{R}\right)}}{4r}$
(33)
$f_{2}(r)=\frac{\left(-R^{2}+r^{2}\right)\left(2R^{2}-2\left(\frac{R^{2}}{r^{2}}+1\right)\ln{\left(R\right)}-1-\frac{R^{2}+2}{r^{2}}\right)}{4R^{4}+4}+\frac{\ln{\left(\frac{r}{R}\right)}}{2r^{2}}.$
(34)

Here we have computed $f_{l}$ as $D_{l}^{-1}u_{l}$ and $Pe_{l}$ is obtained
from condition $\langle v_{l},D^{-1}u_{l}\rangle=1$, as shown in Appendix B.
It turns out we do not need the explicit expression of $g_{l}$. Indeed, noting
that because $g_{l}=D^{-1}v_{l}$ (see Appendix B), any inner product in the
form $\langle g_{l},H(r)\rangle$ can be written as

$\langle g_{l},H(r)\rangle=\langle\hat{D}_{l}^{-1}v_{l},H(r)\rangle=\langle
v_{l},\hat{D}_{l}^{-1}H(r)\rangle=Pe_{l}[\hat{D}_{l}^{-1}H(r)]_{r=1}$ (35)

Recall that $\hat{D}_{l}$ is self adjoint and so is its inverse. The
calculation of $\hat{D}_{l}^{-1}$ is performed in Appendix A. We are now in a
position to calculate the inner product in (30). Consider the case $l=1$, and
collect the linear term in $c_{1}$ in (22), which is given by $(Pe-
Pe_{1})u_{1}c_{1}(1)=(Pe-Pe_{1})C_{1}(t)u_{1}(r)f_{1}(1)$ (where $u_{1}(r)$ is
given by given by (11)). Using (30), we obtain from the right hand side

$(Pe-Pe_{1})C_{1}(t)f_{1}(1){\langle g_{1}(r),u_{1}(r)\rangle\over\langle
g_{1}(r),f_{1}(r)\rangle}=(Pe-
Pe_{1})C_{1}(t)f_{1}(1){[\hat{D}_{1}^{-1}u(r)]_{r=1}\over[\hat{D}_{1}^{-1}f_{1}(r)]_{r=1}}$
(36)

where we have used (35). In Appendix A wee show how to calculate $D_{1}^{-1}$,
and the above expression can easily be evaluated as a function of $R$.
Reporting this into (30) yields the $C_{1}(t)$ equation to linear order. The
next term from $q_{l}$ is quadratic in $c_{l}$ and comes from the combination
of second (advection term) and third terms in (22). It reads, after using the
expressions of $u_{r}$ and $u_{\theta}$ (see (4), where $\psi_{l}$ is given by
(6))

$\sum_{m\neq 1}{mc_{m}(1,t)Pe\over 2r^{|m|+1}}\left[(r^{2}-1)m\frac{\partial
c_{1-m}(r,t)}{\partial r}+(m-1)(2r^{2}+(1-r^{2})|m|)c_{1-m}(r,t)\right]$ (37)

By using $c_{l}(r,t)=C_{l}(t)f_{l}(r)$, and retaining only first and second
harmonics, it is easy to see that the result reads as $C_{1}^{*}C_{2}h(r)$,
where $h(r)$ is a function of $r$ only (it is a combination of $f_{1}$ and
$f_{2}$ and their derivatives with respect to $r$). Once this expression is
injected into (30) on the right hand side we obtain

$C_{1}(t)^{*}C_{2}(t){\langle g_{1}(r),h(r)\rangle\over\langle
g_{1}(r),f_{1}(r)\rangle}=C_{1}(t)^{*}C_{2}(t){[\hat{D}_{1}^{-1}h(r)]_{r=1}\over[\hat{D}_{1}^{-1}f_{1}(r)]_{r=1}}$
(38)

The last term can easily be obtained as algebraic rational functions of $R$
and $\ln R$ (see Appendix B). The same reasoning can be made for the equation
of $C_{2}$. The nonlinear term is found to be proportional to $C_{1}(t)^{2}$.
Collecting linear and nonlinear terms in $C_{1}$ and $C_{2}$, the resulting
system of equations to the second order is found to be given by

$\displaystyle\dot{C}_{1}$
$\displaystyle=\sigma_{1}C_{1}+\alpha_{1}C_{2}C_{1}^{*}+O(\epsilon^{3}),$
(39a) $\displaystyle\dot{C}_{2}$
$\displaystyle=\sigma_{2}C_{2}+\alpha_{2}C_{1}^{2}+O(\epsilon^{3}),$ (39b)

where different coefficients are functions of $R$ and are listed in Appendix
C.

##  VIII Equations for $C_{1}$ and $C_{2}$ to next order

The next order terms turn out to be essential for nonlinear saturation. We
need thus to extend the derivation to $O(\epsilon^{3})$. We use (29) and
insert it into (30). Taking into account $\delta c_{l}$ in (29) will lead to
higher order terms. $\delta c_{l}$ is given by (25) and (26) for
$l\in\\{1,2\\}$ and by (24) for $l\not\in\\{1,2\\}$. Inserting these solutions
into (29) allows us, by using (30), to obtain the desired terms. The cubic
terms are of the form $\lvert C_{2}\rvert^{2}C_{1}$ for the equation of
$C_{1}$ and $\lvert C_{1}\rvert^{2}C_{2}$ for the equation of $C_{2}$. There
is also a contribution of $\delta c_{l}$ in the form of $\partial_{t}\delta
c_{l}$ in definition of $q_{l}$. Since the time derivative is small (critical
slowing down), the quadratic contribution arising from $\delta c_{l}$ are
sufficient. For example the equation for $C_{1}$ yields terms in the form
$\partial_{t}(C_{1}^{*}C_{2})=C_{1}^{*}\partial_{t}C_{2}+C_{2}\partial_{t}C_{1}^{*}$.
Using (39) we can express these terms as quadratic and cubic terms. The final
set of equation takes the form

$\displaystyle\dot{C}_{1}$
$\displaystyle=\sigma_{1}C_{1}+\alpha_{1}C_{2}C_{1}^{*}-\beta_{1}|C_{2}|^{2}C_{1}-\gamma_{1}|C_{1}|^{2}C_{1},$
(40a) $\displaystyle\dot{C}_{2}$
$\displaystyle=\sigma_{2}C_{2}+\alpha_{2}C_{1}^{2}-\beta_{2}|C_{1}|^{2}C_{2}-\gamma_{2}|C_{2}|^{2}C_{2},$
(40b)

The coefficients of the cubic terms have expressions in $R$ which become too
involved. For this reason, we have resorted to computing them for a given
numerical value of $R$, becoming just real numbers. In Appendix C we plot
these coefficients as functions of $R$.

##  IX Comparison between the full model and the reduced one

Here we would like to confront the reduced model (40) with the full model (as
described in Section II). The full model is solved using finite difference
discretization. The reduced model (40) is a set of ordinary differential
equations, the numerical solution of which is straightforward. The results of
the comparison are shown in Fig. 1 for the steady-state values of the swimming
speed and the angular velocity of the particle (a nonzero value of the angular
velocity corresponds to a circular trajectory). As can be seen, the reduced
$C_{1}$ and $C_{2}$ model is in quantitative agreement with the full numerical
simulation.

The solution of the $C_{1}$ and $C_{2}$ model shows that as the Péclet number
is increased, the particle shows non-motile solution as the stable fixed point
for $Pe<Pe_{1}$, straight motion (velocity is constant in time and finite, but
angular velocity is equal to zero) for $Pe^{*}>Pe>Pe_{1}$, and a circular
trajectory (both the swimming speed and the angular velocity remain fixed in
time) for $Pe^{**}>Pe>Pe^{*}$. No stable fixed point solution could be
obtained for $Pe>Pe^{**}$. Here we have $Pe^{*}\approx 5.77$ and
$Pe^{**}\approx 5.9$. By increasing further $Pe$ we find that the solution
becomes irregular, pointing to the occurrence of chaos (Fig.2 . This
succession of transitions from one trajectory type to another reproduces the
prediction of the phenomenological theory. The circular trajectory and chaotic
one have been also obtained by numerical simulation in Hu  _et al._ (2019)
using the full model. In a recent work Misbah  _et al._ (2021) by taking the
set of equations (40) as a phenomenological model (extracted on the basis of
symmetries, without reference to any given basic model) it has been shown that
it reproduces a variety of solutions going from straight, circular to chaotic
trajectoriees. The circular trajectory could even be obtained analytically
from (40).

Figure 1:  Comparison between the analytical prediction and the direct
numerical simulation of the particle velocity and angular velocity in
stationary, straight, and circular phases. The system size is set to $R=3.25$.
The solid curves are obtained by a direct numerical solution of the $C_{1}$,
$C_{2}$ equations. The symbols are the results of the full numerical
simulations. In both cases the steady-state values are obtained by running the
simulations for a long time until (on the order of $10^{5}a^{2}/D$). The
horizontal axis is cut at $Pe=5.9$ above which the absolute value of the
velocity does not seem to reach a steady state.  Figure 2:  An irregular
trajectory obtained from solution of (40) for $Pe=5.95$.

##  X Conclusion

We have extracted from a phoretic model in 2D, including advection and
diffusion, a reduced model in terms of ordinary differential equations for the
first two Fourier modes. This model captures, in a much less numerically
expensive manner, essential features encountered in the full model. The
derivation has been performed for the two first modes to cubic order.
Extension to a higher number of modes can also be performed, as well as going
to higher orders than cubic terms. This will be essential if one wishes to
have a wider applicability of the method. The derivation can be extended to 3D
without any additional conceptual complications. The method has been
exemplified for a specific phoretic model. However, the same technique can be
used for any swimmer powered by one or many chemical fields. Examples of great
topicality are found in models of cell motility, involving acto-myosin
kinetics inside the cell, outside and on the cell membrane, cortex flow as
well as flow in the suspending and internal fluids Hawkins  _et al._ (2011);
Callan-Jones  _et al._ (2016); Farutin  _et al._ (2019). The present theory
constitutes a precious alternative to full numerical simulations. In addition,
as the reduced model equations have a form which depends only on symmetries,
different active entities can be put together in a universal framework, where
different models can be described by the same reduced evolution equations
(40), and their differences will be encoded in the values of the numerical
coefficients only.

##  XI Acknowledgments

We thank CNES (Centre National d’Etudes Spatiales) for financial support and
for having access to microgravity data, and the French-German university
program ”Living Fluids” (grant CFDA-Q1-14) for financial support.

##  Appendix A Inverting the diffusion operators

A solution to the equation

$\frac{1}{r}\frac{d}{dr}\left(r\frac{dF}{dr}\right)-l^{2}F=G(r)$ (41)

can be written as

$F(r)=\begin{cases}\frac{1}{2l}\left[r^{-l}\int^{r}G(\rho)\rho^{1+l}d\rho-r^{l}\int^{r}G(\rho)\rho^{1-l}d\rho\right]+A_{1}r^{l}+A_{2}r^{-l}&\textrm{
for }l>0\\\ \ln r\int^{r}G(\rho)\rho d\rho-\int^{r}\rho\ln\rho
G(\rho)d\rho+A_{1}+A_{2}\ln r&\textrm{ for }l=0\end{cases},$ (42)

where the constants $A_{1}$ and $A_{2}$ are calculated by imposing the
boundary conditions

$\partial_{r}F(1)=0\,\,\,F(R)=0.$ (43)

The problem considered here requires solving eq. (41) with functions $G(r)$
written as a combination of integer powers of $r$ and $\ln r$. This allows us
to calculate the integrals in (42) analytically.

##  Appendix B Inverting the advection-diffusion operator

Suppose we have a system of equations

$\begin{split}&\mu\bm{f}=(M+\bm{u}\otimes\bm{v}).\bm{x}+\bm{y}\\\
&\langle\bm{g},\bm{x}\rangle=0,\end{split}$ (44)

where $\mu$ is a number, $M$ is a matrix, $\bm{u},$ $\bm{v}$, $\bm{x}$,
$\bm{y}$, $\bm{f}$, $\bm{g}$ are vectors, and

$\begin{split}&(M+\bm{u}\otimes\bm{v}).\bm{f}=0,\\\
&\bm{g}.(M+\bm{u}\otimes\bm{v})=0.\end{split}$ (45)

The first equation (44) corresponds to eq.(21) where $\mu$ represents
$\partial_{t}C_{l}$, $M+\bm{u}\otimes\bm{v}$ is the linear operator
$\hat{L}(Pe_{l})$ (see eq; 23), ${\bm{x}}$ represents $\delta c_{l}$ and
${\bm{y}}$ represents $q_{l}$. The second equation represents the
orthogonality condition (17). Equation (45) represents the eigenvalue problem
$\hat{L}(Pe_{l})f_{l}(r)=0$ and its adjoint (15).

The goal of this Section is to find the solution of the system (44)
representing $\bm{x}$ and $\mu$ as a function of $M^{-1}$, $\bm{u}$, $\bm{v}$,
and $\bm{y}$. First, we note the following relations:

$\bm{f}\propto
M^{-1}.\bm{u},\,\,\,\bm{g}\propto\bm{v}.M^{-1},\,\,\,\langle\bm{v},M^{-1}.\bm{u}\rangle=-1,$
(46)

which we obtain from eqs. (45). Multiplying the first equation in (44) by
$M^{-1}$, we get

$\mu M^{-2}.\bm{u}=[I+(M^{-1}.\bm{u})\otimes\bm{v}].\bm{x}+M^{-1}.\bm{y},\\\ $
(47)

whence

$\bm{x}=\mu
M^{-2}.\bm{u}-(M^{-1}.\bm{u})\langle\bm{v},\bm{x}\rangle-M^{-1}.\bm{y},\\\ $
(48)

Equation (48) implies the following ansatz for $\bm{x}$:

$\bm{x}=-M^{-1}.\bm{y}+pM^{-1}.\bm{u}+qM^{-2}.\bm{u},\\\ $ (49)

where $p$ and $q$ are two numbers to be determined. It is then straightforward
to substitute eq. (49) into eq. (48) to get the values of $\mu$ and $q$ as:

$q=\mu=\frac{\langle\bm{v},M^{-1}.\bm{y}\rangle}{\langle\bm{v},M^{-2}.\bm{u}\rangle}$
(50)

The value of $p$ is determined from the second equation of (44). It is
convenient that the solution (50) use the vector $\bm{v}$ only as part of
$\langle\bm{v},\bm{\cdot}\rangle$ which allows us to use its definition
$\langle\bm{f},\bm{v}\rangle=Pef(1)$.

##  Appendix C  Coefficients of the reduced model

Here we list the coefficients of the linear and quadratic terms and plot those
of the cubic term.

$\sigma_{1}=\frac{16(Pe-
Pe_{1})\left(R^{2}\ln{\left(R\right)}-R^{2}+\ln{\left(R\right)}+1\right)^{2}}{3R^{6}+8R^{4}\ln{\left(R\right)}^{2}-40R^{4}\ln{\left(R\right)}+19R^{4}+8R^{2}\ln{\left(R\right)}^{2}+4R^{2}\ln{\left(R\right)}-35R^{2}+12\ln{\left(R\right)}+13}.$
(51)
$\alpha_{1}=-\frac{\begin{matrix}2\left(R^{2}+1\right)\left(20R^{10}\ln{\left(R\right)}^{2}-40R^{10}\ln{\left(R\right)}+21R^{10}+8R^{8}\ln{\left(R\right)}^{2}+48R^{8}\ln{\left(R\right)}-63R^{8}\right.&\\\
\left.+16R^{6}\ln{\left(R\right)}^{3}-48R^{6}\ln{\left(R\right)}^{2}+56R^{6}\ln{\left(R\right)}+46R^{6}+32R^{4}\ln{\left(R\right)}^{3}-56R^{4}\ln{\left(R\right)}^{2}\right.&\\\
\left.-8R^{4}\ln{\left(R\right)}+6R^{4}+16R^{2}\ln{\left(R\right)}^{3}-20R^{2}\ln{\left(R\right)}^{2}-48R^{2}\ln{\left(R\right)}-3R^{2}-8\ln{\left(R\right)}-7\right)\end{matrix}}{\begin{matrix}R^{2}\left(R^{4}+1\right)\left(R^{2}\ln{\left(R\right)}-R^{2}+\ln{\left(R\right)}+1\right)\left(3R^{6}+8R^{4}\ln{\left(R\right)}^{2}-40R^{4}\right.&\\\
\left.\ln{\left(R\right)}+19R^{4}+8R^{2}\ln{\left(R\right)}^{2}+4R^{2}\ln{\left(R\right)}-35R^{2}+12\ln{\left(R\right)}+13\right)\end{matrix}}$
(52)

$\sigma_{2}=\frac{3(Pe-
Pe_{2})\left(R^{4}-4R^{2}+4\ln{\left(R\right)}+3\right)^{2}}{12R^{8}\ln{\left(R\right)}-13R^{8}+16R^{6}\ln{\left(R\right)}+8R^{6}+12R^{4}\ln{\left(R\right)}+6R^{4}-8R^{2}+8\ln{\left(R\right)}+7}$
(53)

$\alpha_{2}=\frac{\begin{matrix}3\left(R^{4}+1\right)^{2}\left(R^{2}\ln{\left(R\right)}-R^{2}+\ln{\left(R\right)}+1\right)\left(6R^{6}\ln{\left(R\right)}-7R^{6}+6R^{4}\ln{\left(R\right)}\right.&\\\
\left.+9R^{4}+6R^{2}\ln{\left(R\right)}-9R^{2}+6\ln{\left(R\right)}+7\right)\end{matrix}}{\begin{matrix}\left(R^{2}+1\right)^{2}\left(R^{4}-4R^{2}+4\ln{\left(R\right)}+3\right)\left(12R^{8}\ln{\left(R\right)}\right.&\\\
\left.-13R^{8}+16R^{6}\ln{\left(R\right)}+8R^{6}+12R^{4}\ln{\left(R\right)}+6R^{4}-8R^{2}+8\ln{\left(R\right)}+7\right)\end{matrix}}$
(54)

##  Appendix D Explicit expression

Here we list the explicit epxressions for $\partial_{t}C_{1}$ and
$\partial_{t}C_{2}$ for $R=3.25$. We first have simplified the expressions
assuming $Pe$ to be close to the critical Péclet number $Pe_{c}=5.9561$, which
is the critical Péclet number for $R=R_{c}=3.17493$, such that
$Pe_{1}(R_{c})=Pe_{2}(R_{c})=Pe_{c}$. We then define $\Delta P=Pe-
Pe_{c}=O(\epsilon)$ and truncate all expressions keeping only terms of order
$O(\epsilon^{3})$ or higher.

$\displaystyle\partial_{t}C_{1}$ $\displaystyle=C_{1}\left(0.0484+0.1799\Delta
Pe-0.00244\Delta Pe^{2}\right)-C^{*}_{1}C_{2}\left(1.227+0.2467\Delta
Pe\right)-$ (55a) $\displaystyle
0.2016C_{1}|C_{1}|^{2}-3.077C_{1}|C_{2}|^{2}+O(\epsilon^{4}),$
$\displaystyle\partial_{t}C_{2}$ $\displaystyle=C_{2}\left(0.0411+0.3701\Delta
Pe-0.00123\Delta Pe^{2}\right)+C_{1}^{2}\left(0.5921+0.0905\Delta Pe\right)-$
(55b) $\displaystyle 3.1214|C_{1}|^{2}C_{2}-1.7893|C_{2}|^{2}C_{2}.$

## References

  * Callan-Jones  _et al._ (2016) Callan-Jones, A. C., Ruprecht, V., Wieser, S., Heisenberg, C. P., and Voituriez, R., Phys. Rev. Lett. 116, 028102 (2016).
  * Chen  _et al._ (2021) Chen, Y., Chong, K. L., Liu, L., Verzicco, R., and Lohse, D., Journal of Fluid Mechanics 919, A10 (2021).
  * Farutin  _et al._ (2019) Farutin, A., Etienne, J., Misbah, C., and Récho, P., Phys. Rev. Lett. 123, 118101 (2019).
  * Hawkins  _et al._ (2011) Hawkins, R. J., Poincloux, R., Bénichou, O., Piel, M., Chavrier, P., and Voituriez, R., Biophys. J. 101, 1041 (2011).
  * Hokmabad  _et al._ (2021) Hokmabad, B. V., Dey, R., Jalaal, M., Mohanty, D., Almukambetova, M., Baldwin, K. A., Lohse, D., and Maass, C. C., Phys. Rev. X 11, 011043 (2021).
  * Hu  _et al._ (2021) Hu, W.-F., Lin, T. S., Rafai, S., and Chaouqi, M., arXiv preprint (2021).
  * Hu  _et al._ (2019) Hu, W. F., Lin, T. S., Rafai, S., and Misbah, C., Phys. Rev. Lett. 123, 238004 (2019).
  * Izri  _et al._ (2014) Izri, Z., Van Der Linden, M. N., Michelin, S., and Dauchot, O., Phys. Rev. Lett. 113, 248302 (2014).
  * Izzet  _et al._ (2020) Izzet, A., Moerman, P. G., Gross, P., Groenewold, J., Hollingsworth, A. D., Bibette, J., and Brujic, J., Phys. Rev. X 10, 021035 (2020).
  * Jin, Krüger, and Maass (2017) Jin, C., Krüger, C., and Maass, C. C., Proceedings of the National Academy of Sciences 114, 5089 (2017), https://www.pnas.org/content/114/20/5089.full.pdf .
  * Michelin and Lauga (2014) Michelin, S. and Lauga, E., Journal of Fluid Mechanics 747, 572 (2014).
  * Michelin, Lauga, and Bartolo (2013) Michelin, S., Lauga, E., and Bartolo, D., Phys. Fluids 25, 061701 (2013).
  * Misbah  _et al._ (2021) Misbah, C., Rizvi, M. S., Hu, W.-F., Lin, T. S., Rafai, S., and Farutin, A., arXiv preprint (2021).
  * Morozov (2020) Morozov, M., Soft Matter 16, 5624 (2020).
  * Morozov and Michelin (2019a) Morozov, M. and Michelin, S., J. Chem. Phys. 150, 044110 (2019a).
  * Morozov and Michelin (2019b) Morozov, M. and Michelin, S., Journal of Fluid Mechanics 860, 711 (2019b).
  * Schmitt and Stark (2013) Schmitt, M. and Stark, H., Euro. Phys. Lett. 101, 44008 (2013).


11institutetext: National University of Defense Technology
11email<EMAIL_ADDRESS>
# CoG: a Two-View Co-training Framework for Defending Adversarial Attacks on
Graph

Xugang Wu 11 Huijun Wu 11 Xu Zhou 11 Kai Lu 11

###### Abstract

Graph neural networks (GNNs) exhibit remarkable performance in graph data
analysis. However, the robustness of GNN models remains a challenge. As a
result, they are not reliable enough to be deployed in critical applications.
Recent studies demonstrate that GNNs could be easily fooled with adversarial
perturbations, especially structural perturbations. Such vulnerability is
attributed to the model’s excessive dependence on the structure information to
make predictions. To achieve better robustness, it is desirable to build the
prediction of GNNs with more comprehensive features. Graph data, in most
cases, has two views of information, namely structure information and feature
information. In this paper, we propose CoG, a simple yet effective co-training
framework to combine these two views for the purpose of robustness. CoG trains
sub-models from the feature view and the structure view independently and
allows them to distill knowledge from each other by adding their most
confident unlabeled data into the training set. The orthogonality of these two
views diversifies the sub-models, thus enhancing the robustness of their
ensemble. We evaluate our framework on three popular datasets, and results
show that CoG significantly improves the robustness of graph models against
adversarial attacks without sacrificing their performance on clean data. We
also show that CoG still achieves good robustness when both node features and
graph structures are perturbed.

##  1 Introduction

Graph Neural Networks (GNNs) have achieved remarkable success in the learning
representations for graph-structured data across many domains, such as
citation networks, biological networks, and social networks. Among these
models, Graph Convolutional Network (GCN) and its variants [17, 30, 12, 10]
have gained much attention because of their performance and efficiency.
However, recent studies demonstrate that these message-passing-based models
suffer from adversarial perturbations. The accuracy of the learned model drops
drastically when the attacker conducts unnoticeable modifications on the graph
data [6, 37, 2, 38]. Further results [33, 38] illustrate that structural
perturbations lead to more effective attacks compared to feature
perturbations.

Figure 1: Performance of ensemble of GCN and GAT under Metattack.

Some defense techniques are proposed to enhance the robustness of GNNs against
structural perturbations. The general principle is to eliminate the negative
impacts brought by perturbed edges based on certain prior knowledge. For
example, some defense methods [36, 4, 20, 33] assume that all the connected
nodes should be similar in the feature space. Based on it, they remove or pay
less attention to edges that connect dissimilar nodes. Other methods [15, 9]
observe that structural perturbations tend to affect the high-rank portion of
the graph. Therefore, they replace the graph topology with its low-rank
approximation to purify the graph. However, these heuristic knowledge may not
always hold in all scenarios. For example, in heterophilic graphs, removing
edges between dissimilar nodes could impair the model’s performance. Another
drawback of heuristic defenses is their weakness against adaptive attacks.
Among these methods, few take adaptive attacks into consideration. In this
paper, instead of introducing heuristic knowledge, we propose to defend
against poisoning attack by learning an ensemble model that fully exploits the
feature and structure information.

Ensemble method is a prevalent way to defend against adversarial attacks in
the field of computer vision. Recent works [23, 16, 35] point out that output
diversity is the key to the success of these ensemble-based defense methods.
However, applying ensemble techniques to GNN straightforwardly may not improve
the robustness of the ensemble model. Fig. 1 shows how the ensemble of GCN and
GAT performs under Metattack. Metattack can transfer between GCN and GAT well,
leading to poor robustness of their ensemble.

In this paper, We argue that introducing diversified models that share
different views of graph data will help to improve the overall robustness of
the ensemble. Compered to image data, graph data have two orthogonal views:
the node feature view and the graph structure view. The orthogonality of these
two views makes it difficult for adversarial attacks to transfer between them.
To fully exploit this property, our work uses a two-view co-training framework
(CoG) to learn an ensemble of sub-models from both feature view and structure
view. Co-training [1, 32] is a simple yet effective technique for learning an
ensemble of sub-models in two different views under semi-supervised setting.
It trains a separate classifier for each view and adds the most confident
predictions of each sub-model to the training data set. During this process,
sub-models distill knowledge from each other and get strengthened.
Nevertheless, applying vanilla co-training framework to graph data meets two
challenges: (1) Co-training uses the softmax outputs as the indicator of sub-
models’ confidence. However, this can be inaccurate since neural networks are
often miscalibrated, especially when sub-models are heterogenous. (2) Co-
training selects unlabeled data simply based on their confidence. If dominant
classes exist, the co-training process will amplify the imbalance of classes
and force the sub-models to overfit to the dominant class. To address these
issues, we use temperature scaling to calibrate the outputs, and enforce the
constant of class distribution when adding predictions during the co-training
process.

Experiment results show that the diverse outputs between feature and structure
view prevent adversarial examples from transferring between them, thus making
their ensembles more robust against adversarial attacks, even in the scenario
of adaptive attacks. To conclude, our contributions are summarized as follows:

  * •

To employ the orthogonality between feature and structure views, we propose a
two-view co-training framework (CoG) to combine the feature information and
the structure information during the training process.

  * •

We show that model calibration and class balancing play significant roles in
enabling effective co-training.

  * •

Experiments show that CoG can achieve better robustness against poisoning
attacks compared to other state-of-the-art methods, without sacrificing the
performance on clean data.

  * •

We conduct adaptive attacks and results indicate that our method can still
work well under the adaptive setting.

##  2 Related work

###  2.1 Attack and defense on graph data

Despite the great success of GNNs, recent works indicate that these graph-
based models are vulnerable to unnoticeable modifications [37, 38, 6]. Nettack
[37] conducts its attack on a surrogate model and ensures their modifications
to be unnoticeable via considering the degree distribution and feature co-
occurrence. RL-S2V [6] applies reinforcement learning to generate adversarial
examples. Metattack [38] addresses the global attack problem. It uses the
meta-gradient to generate a perturbed graph that leads to an overall decrease
in model’s performance. All these attack methods suggest that structural
attacks are more effective than feature attacks when applied to graph data.

Several techniques have been proposed to defend these topological attacks [33,
36, 15, 14, 4, 20]. GCN-Jaccard [33] assumes that connected nodes should have
similar features so that edges between dissimilar nodes are more likely to be
the perturbed edges. Based on this assumption, they calculate the Jaccard
similarity scores between all connected nodes and remove edges that connect
nodes with scores below the threshold. Similarly, GNNGuard [36] employs the
attention mechanism to assign higher weights to the edges between similar
nodes and less to the edges between unrelated nodes. Other works stem from the
observation that structural attacks tend to affect the high-rank portion of a
graph. GCN-SVD [9] is proposed to purify the perturbed graph via replacing it
with its low-rank approximation, while Pro-GNN [15] introduces a
regularization term to generate a low-rank and sparse graph during the
training process. However, these methods rely on the validity of their
heuristic knowledge. SimP-GCN [14] attempts to integrate the structure
information and node features by combining the kNN graph and the original
graph. Nevertheless, the scoring function that balances these two graphs
merely depends on the hidden representation so that the graph used for
learning could be unstable, thus leading to high variance for the results.

The most related defense method is UM-GNN [26]. UM-GNN is proposed to learn a
feature-based model via distilling knowledge from GNN models using an
uncertainty matching strategy. However, its knowledge distillation is one-
directional: only the feature-based model can distill knowledge from GNN
models. Consequently, GNN models cannot enhance itself using the information
from the feature-based model. As the perturbation rate grows, the knowledge
transferred from GNN models gets less effective, which impairs the performance
of UM-GNN.

###  2.2 Ensemble training for enhanced robustness

Although ensemble training was initially proposed to improve models’
performance [13, 3, 8, 18], a recent line of works [23, 16, 35] show that it
can be applied to enhance adversarial robustness. The intuition behind these
methods is that a small overlap between adversarial subspace (Adv-SS) of
different sub-models can prevent adversarial attacks from transferring between
sub-models. Pang et al. [23] employ an adaptive diversity-promoting
regularizer to encourage diversity among non-maximal predictions. Kariyappa et
al. [16] propose diversity training to reduce the correlation of loss
functions between sub-models. Yang et al. [35] distill the non-robust features
from each sub-model and teach other sub-models to be robust against these non-
robust features. Although these defense approaches have been well studied in
image recognition tasks, their application in graph-based tasks remains to be
explored.

##  3 Preliminary study

In this section, we present why two-view co-training training is a desirable
defense technique for graph data. As is discussed in the previous section,
output diversity plays a crucial role in the robustness of the ensemble model.
Recent work [29] shows that adversarial examples are more likely to transfer
between models that have a large adversarial subspace (Adv-SS) overlap. Sub-
models with diverse outputs have a smaller Adv-SS overlap, making it harder
for attackers to craft adversarial examples that fool both two models [23, 16,
35].

To minimize the Adv-SS overlap, using different aspects of the input data is a
promising approach. Graph data contains two orthogonal views: a node feature
view and a graph structure view. Models in these two views naturally share
little adversarial subspace. An attack on node features would barely affect a
structure-based model and vice versa, which implies the potential for
improving robustness via an ensemble of sub-models from these two views.
Furthermore, from attackers’ perspective, conducting attacks on the feature
view is more difficult due to the following reasons: (1) in most graph data,
node features are sparse and high-dimensional, which makes the perturbations
on node features detectable; (2) in most cases, node features are discrete and
practically interpretable, which causes restrictions on feature modifications.
Suppose an attacker attempts to attack a social network via modifying a
person’s birthday, his/her age should be simultaneously changed to ensure the
data validity. Such restrictions are difficult to model in a differentiable
manner, which constitutes a challenge for gradient-based attacks.

Although existing message-passing-based GNN models are designed to integrate
the feature and structure information of graph data, they do not pay enough
attention to the feature information. In an empirical study, Jin et al. [15]
point out that GCN prefers to preserve structure information rather than
feature information during the message passing process. We call it a
structure-dominant model. For these models, the node’s representation is
highly related to its neighboring structure rather than its feature. This
behavior is consistent with the empirical conclusions in previous work [33]:
(1) Topological attacks are more effective than feature attacks when attacking
models like GCN; (2) Attackers tend to connect nodes with dissimilar node
features to achieve a successful attack.

Therefore, to fully exploit the feature information and achieve better
robustness, we propose a two-view co-training framework, named Co-training on
Graph (CoG), to learn a robust ensemble of the feature-dominant models and the
structure-dominant models.

Figure 2: The overall framework of CoG. Graph data is seperated into a feature
view and a structure view. During the co-training process, sub-models from
each view label their most confident nodes and add them into the training set.

##  4 Proposed approach

###  4.1 Problem formulation

In this paper, we focus on the semi-supervised node classification problem.
Let $\mathcal{G}=(\mathcal{V},\mathcal{E},X)$ be a graph with n nodes, where
$\mathcal{V}$ is the set of nodes $\left\\{v_{1},\cdots,v_{n}\right\\}$ with
$|\mathcal{V}|=n,\mathcal{E}$ is the set of edges, and
$X=\left[x_{1},x_{2},\cdots,x_{n}\right]^{T}\in\mathbb{R}^{n\times m}$ is a
feature matrix. Then we can separate $\mathcal{G}$ into two views. In the node
feature view, for each node $v\in\mathcal{V}$, its feature
$x_{v}\in\mathbb{R}^{m}$ is a $m$-dimensional row vector. In the graph
structure view, the adjacency matrix $A\in\mathbb{R}^{n\times n}$ can be
formulated by setting $A_{ij}=1$ if $\left(v_{i},v_{j}\right)\in\mathcal{E}$
and $A_{ij}=0$, otherwise. In the semi-supervised node classification problem,
nodes are separated into two sets $\mathcal{V}=\mathcal{S}\cup\mathcal{U}$,
where nodes in $\mathcal{S}$ are labeled and nodes in $\mathcal{U}$ are not.
Our goal is to learn an high-performance ensemble of a feature-dominant model
$f_{feat}$ and a structure-dominant model $f_{struct}$, using both the labeled
and unlabeled data from $\mathcal{V}$.

###  4.2 The overall framework

Graph data can be separated into two views, i.e., the feature view and the
structure view. These two views provide different information about the nodes.
The feature view describes nodes’ intrinsic properties, while the structure
view tells the relationship between each node pair. Furthermore, they meet the
requirements[18] of a co-training framework well. (1) Models trained in any of
these two views can classify the nodes effectively. (2) Information provided
by these two views is independent so that knowledge distilled from one model
can promote the performance of the other model. Therefore, we apply the two-
view co-training framework on graph data (CoG) and learn an ensemble of sub-
models in these two views.

The overall framework of CoG is shown in Fig. 2. For each node $v$, we
separate its information into a feature view $X$ and a structure view $A$.
After that, two classifiers are trained separately for each view. Feature-
dominant models primarily use node features as its input to classify the
nodes, while structure-dominant models usually use the graph structure as its
clues for node labels. In each iteration, the most confident unlabeled nodes
from each view are added to the training set. After a preset number of
iterations, we obtain two well-trained models, $f_{feat}$ and $f_{struct}$,
and generate an ensemble of them by averaging their predictions. In the
remaining part of this section, we introduce more details about the framework,
including (1) the structure-dominant models and feature-dominant models we use
to learn representation for nodes; (2) the techniques we use to improve the
performance of the co-training framework, including model calibration and
class balancing.

###  4.3 Structure-dominant models

####  4.3.1 GCN

Graph convolutional network (GCN) [17] achieves remarkable success on learning
representation for graph data. Its message-passing mechanism learns the
representation of each node via propagating information and aggregating with a
permutation invariant function. Given a graph $\mathcal{G}=(X,A)$, the updates
of node embeddings can be derived by the following formulation:

$X^{(l+1)}=\sigma\left(\hat{A}X^{(l)}W^{(l)}\right),$ (1)

where $\hat{A}$ is the normalized adjacency matrix with self loops
$\hat{A}=\tilde{D}^{-\frac{1}{2}}(A+I)\tilde{D}^{-\frac{1}{2}}$, $I$ is the
identity matrix, $\tilde{D}_{ii}=\sum_{j}\tilde{A}_{ij}$, $W^{(l)}$ is the
feature transformation matrix of the $l$th layer, and $\sigma(\cdot)$ is the
activation function, e.g., ReLU(x). Although GCN takes both features and graph
structure as its input, the representation it learns for each node largely
depends on the local graph structure of the node.

####  4.3.2 Spectral-based method (S-MLP)

Spectral-based method is another effective method to distill structure
information for each node. Utilizing the spectral decomposition of the graph
laplacian to represent both local and global structure has been thoroughly
studied in the graph spectral theory [5, 31]. Given the graph topology $A$, we
first compute its eigenvalues and eigenvectors of laplacian by solving:

$D^{-1}L\mathbf{y}=\lambda\mathbf{y}$ (2)

where $D^{-1}L=I-D^{-1}A$ is the normalized laplacian matrix. Let
$\lambda_{0},\lambda_{1},\cdots,\lambda_{k}$ be the k smallest eigenvalues of
$D^{-1}L$ and $\mathbf{y}_{0},\mathbf{y}_{1},\cdots,\mathbf{y}_{k}$ be the
respective eigenvectors, each node $u$ can be embedded into a k-dimension
space as $a_{u}=(\mathbf{y}_{1u},\mathbf{y}_{2u},\cdots,\mathbf{y}_{ku})$.

To enrich the information of $a_{u}$, we also compute the k-dimension
laplacian eigenmaps of $A^{2}$. By concatenating them together, we get the
enhanced $a_{u}\in\mathbb{R}^{2\times k}$. After that, $a_{u}$ is sent as
inputs to an MLP model for classification.

###  4.4 Feature-dominant models

####  4.4.1 MLP (F-MLP)

Multi-layer perception (MLP) model is the simplest while effective method when
we consider the node features alone. The layerwise forwarding process of MLP
can be formulated as,

$x_{v}^{(l+1)}=\sigma(x_{v}^{l}\Theta^{(l)}+b_{l})$ (3)

where $x_{v}\in\mathbb{R}^{m}$ is a $m$-dimensional row vector denotes the
features of node $v$, $\Theta^{(l)}$ and $b_{l}$ are the layer-wise
parameters. $\sigma(\cdot)$ denotes the activation function.

####  4.4.2 kNN-based GCN (kNN-GCN)

In kNN-based model, we first construct a feature graph from feature matrix $X$
by conducting k-nearest-neighbor algorithm based on the cosine similarity. For
each node pair $(v_{i},v_{j})$, we calculate their feature similarity as:

$\mathrm{s}_{ij}=\frac{\mathrm{x}_{i}^{\top}\mathrm{x}_{j}}{\left\|\mathrm{x}_{i}\right\|\left\|\mathrm{x}_{j}\right\|}.$
(4)

After that, the kNN graph $A_{k}=kNN(X)$ can be constructed by connecting the
top-k similar node pairs. Finally, the whole graph data
${\mathcal{G}=(X,A_{k})}$ is sent to a GCN model to classify the nodes. The
structure information here is generated from the node features, so the
original graph structure is not used in this model.

###  4.5 Model calibration

Models’ confidence is the most important indicator during the co-training
process. On one hand, the co-training framework picks up the most confident
predictions in each sub-model. When sub-models attempt to add the same node to
the training set, the co-training framework decides which pseudo label to use
based on their confidence. On the other hand, in the inference stage, CoG
average each sub-model’s confidence to obtain the predictions of their
ensemble.

Generally, the softmax output of each sub-model is used to measure the
confidence. However, [11, 28] show that modern neural networks, including GNN,
can be miscalibrated. Since the sub-models are heterogeneous in CoG, this
miscalibration can impair the performance of the co-training method. To
mitigate the impacts of model’s miscalibration, we use the temperature scaling
method [11] to calibrate the output of each sub-model. Given the logit vector
$z$, the calibrated prediction is obtained as:

$q=\operatorname{softmax}(z/T),$ (5)

where the temperature $T$ is learned by optimization with respect to negative
log likelihood on the validation set.

###  4.6 Class balancing

In CoG, the most confident unlabelled nodes from each view are added to the
training set during the co-training process. The confidence of each node is
measured by the softmax outputs of the sub-models. However, this strategy can
lead to class imbalance if sub-models perform better on one particular class.
What’s worse, an imbalanced class distribution will cause the overfitting
problem as the co-traning process moves on and finally impair the sub-models’
performance.

To keep a balanced training dataset, the added data should follow the class
distribution of the initial training dataset. Formally, suppose we have $N$
inputs with an initial distribution $(N_{1},N_{2},\cdots,N_{C})$, the number
of added proposals for each class $C$ is:

$N^{add}_{c}=\frac{N_{c}}{N}\cdot N^{add}$ (6)

where $N^{add}$ denotes the number of data we add to the training data in each
iteration.

The co-training process will stop when reaching the preset number of
iterations or there is no more data in the test set. In the inference phase,
we average sub-models’ output to get the predictions of their ensemble.
Compared to other defense methods, the advantages of CoG can be summarized as
follows:

  * •

Rather than rely on heuristic knowledge to purify the perturbed graph, CoG
enhances the robustness of GNNs via knowledge distillation. Therefore, it
could be applied to more scenarios given the fact that it does not make any
assumptions about the graph.

  * •

Compared to existing ensemble approaches, where the knowledge of the graph
only flows from GNN to MLP in an one-direction manner[26] or model diversity
is not considered [27], CoG enables a bi-directional knowledge distillation
between GNN and MLP as well as takes the complementation of different types of
model into consideration. This enhance the robustness of CoG when the graph
structure is heavily perturbed. Furthermore, the distillation process of CoG
is dynamic and non-differentiable, making it more difficult for attackers to
conduct adaptive attacks.

##  5 Experimental Results

In this section, we evaluate the performance of CoG on clean data and its
robustness against adversarial attacks. In particular, we focus on answering
the following questions:

Q1. How does CoG perform on clean data?

Q2. How does CoG perform under adversarial attacks compared to other state-of-
the-art defense methods?

Q3. How do model calibration, the class balancing technique and the hyper-
parameters affect CoG’s performance?

Q4. How does CoG perform against adaptive attacks?

Table 1: Statistics of the datasets.

| # of nodes | # of edges | Classes | Features
---|---|---|---|---
Cora | 2,485 | 5,069 | 7 | 1,433
Citeseer | 2,110 | 3,668 | 6 | 3,703
Pubmed | 19,717 | 44,338 | 3 | 500

###  5.1 Experiment Setup

####  5.1.1 Datasets

To obtain comparable results, we use three popular citation graphs, i.e., Cora
[21], Citeseer [25] and Pubmed [25] to evaluate our model. The basic
information of these three graphs is shown in Table 1. In these three
datasets, each node represents a document and edges are the citations between
documents. In terms of data splits, we randomly pick 10$\%$ of nodes for
training, 10$\%$ for validation and 80$\%$ for testing, following [15, 37,
38].

####  5.1.2 Baselines

We compare our model with the following baselines:

  * •

GCN [17]: As we describe in the previous sections, GCN is a representative
graph neural network proposed by Kipf and Welling.

  * •

GCN-SVD [9]: GCN-SVD is a defense method based on low-rank approximation.

  * •

GCN-Jaccard [33]: GCN-Jaccard is also a preprocessing-based defense, which
removes the edges between most dissimilar nodes to purify the graph.

  * •

SimP-GCN [14]: SimP-GCN adaptively combines the original graph and the kNN
graph to capture the similarity between nodes.

  * •

Pro-GNN [15]: Pro-GNN is a defense method that exploits three properties of
real-world graphs: low-rank, sparsity, and feature smoothness.

  * •

UM-GNN [26]: UM-GNN trains a feature-based model based on knowledge
transferred from GNN models.

####  5.1.3 Parameter Settings

We implement CoG under the framework of DeepRobust [19], which is a well-known
adversarial learning framework for GNNs. Similar to [15, 4], we run each
experiment 10 times and report the average performance. For GCN, we use the
settings of the original GCN [17], i.e., a two-layer structure with 16 hidden
units. We use a two-layer structure with 32 hidden units for the MLP model.
The number of nearest neighbors $k$ we set in kNN-GCN is 50. As for S-MLP, we
use the eigenvectors corresponding to the lowest 50 eigenvalues to distill
structure information. In the co-training process, each model adds 250 most
confident unlabeled data with their pseudo-labels in one iteration. For the
baseline models, we take the same experiment settings as in [15]. We set the
learning rate for all sub-models to 0.01, the weight decay to 5e-4, the
dropout rate to 0.5, and train for 200 epochs.

As for adversarial attacks, we use Metattack [38] which is an effective
poisoning attack method on graphs. It treats the graph as a hyperparameter and
modifies the graph to increase learning loss via meta-gradient. The edge
perturbation rate is set to $\\{5\%,10\%,15\%,20\%\\}$. We use the same random
seed as [15] to make fair comparisons with their reported results. When
comparing with UM-GNN, we also evaluate CoG under PGD attack [34]. The results
on PGD attack are in the supplementary material.

Table 2: Node classification accuracy ($\%$) on clean graphs and perturbed
graphs.

Datasets | Model | Perturb Rate$(\%)$
---|---|---
0 | 5 | 10 | 15 | 20
Cora | GCN | 83.50±0.44 | 76.55±0.79 | 70.39±1.28 | 65.10±0.71 | 59.56±2.72
GCN-SVD | 80.63±0.45 | 78.39±0.54 | 71.47±0.83 | 66.69±1.18 | 58.94±1.13
GCN-Jaccard | 82.05±0.51 | 79.13±0.59 | 75.16±0.76 | 71.03±0.64 | 65.71±0.89
SimP-GCN | 81.81±0.62 | 76.43±1.98 | 73.27±1.93 | 70.75±3.98 | 66.63±6.87
Pro-GNN | 82.98±0.23 | 82.27±0.45 | 79.03±0.59 | 76.40±1.27 | 73.32±1.56
GCN+F-MLP | 84.09±0.59 | 83.48±0.43 | 82.88±0.83 | 81.01±0.57 | 76.70±0.63
GCN+kNN-GCN | 84.27±0.31 | 83.16±0.28 | 82.54±0.29 | 80.57±0.40 | 76.86±0.75
S-MLP | 78.30±0.17 | 76.25±0.23 | 72.95±0.36 | 67.31±0.40 | 54.56±0.25
| S-MLP+F-MLP | 83.59±0.30 | 84.25±0.33 | 83.52±0.31 | 82.95±0.59 | 80.75±0.59
| S-MLP+kNN-GCN | 82.80±0.33 | 83.85±0.27 | 83.15±0.58 | 82.74±0.23 | 80.74±0.31
Citeseer | GCN | 71.96±0.55 | 70.88±0.62 | 67.55±0.89 | 64.52±1.11 | 62.03±3.49
GCN-SVD | 70.65±0.32 | 68.84±0.72 | 68.87±0.63 | 63.26±0.96 | 58.55±1.09
GCN-Jaccard | 72.10±0.63 | 70.51±0.97 | 69.54±0.56 | 65.95±0.94 | 59.30±1.40
SimP-GCN | 73.76.±0.78 | 73.12±0.85 | 72.38±0.67 | 71.75±1.54 | 69.37±1.50
Pro-GNN | 73.28±0.69 | 72.93±0.57 | 72.51±0.75 | 72.03±1.11 | 70.02±2.28
GCN+F-MLP | 75.14±0.54 | 74.83±0.58 | 73.70±0.60 | 73.68±0.83 | 71.91±0.93
GCN+kNN-GCN | 74.80±0.58 | 74.76±0.56 | 73.79±0.48 | 73.86±0.79 | 71.74±1.34
S-MLP | 69.31±0.34 | 68.90±0.40 | 62.41±0.65 | 62.11±0.38 | 54.99±0.31
| S-MLP+F-MLP | 74.31±0.47 | 74.12±0.38 | 74.82±0.22 | 74.60±0.58 | 69.30±0.69
| S-MLP+kNN-GCN | 74.94±0.41 | 75.04±0.23 | 74.53±0.24 | 75.01±0.22 | 72.80±0.74
Pubmed | GCN | 87.19±0.09 | 83.09±0.13 | 81.21±0.09 | 78.66±0.12 | 77.35±0.19
GCN-SVD | 83.44±0.21 | 83.41±0.15 | 83.27±0.21 | 83.10±0.18 | 83.01±0.22
GCN-Jaccard | 87.06±0.06 | 86.39±0.06 | 85.70±0.07 | 84.76±0.08 | 83.88±0.05
SimP-GCN | 87.59.±0.10 | 86.79±0.12 | 86.01±0.10 | 85.49±0.11 | 85.37±0.12
Pro-GNN | 87.26±0.23 | 87.23±0.13 | 87.21±0.13 | 87.20±0.15 | 87.15±0.15
GCN+F-MLP | 87.62±0.05 | 86.98±0.09 | 86.85±0.09 | 86.05±0.08 | 86.04±0.04
GCN+kNN-CCN | 84.92±0.14 | 83.74±0.15 | 82.97±0.15 | 81.79±0.16 | 81.35±0.08
S-MLP | 77.23±0.17 | 73.29±0.19 | 70.36±0.42 | 66.99±0.75 | 64.68±0.66
| S-MLP+F-MLP | 86.57±0.11 | 86.58±0.10 | 86.51±0.10 | 86.31±0.11 | 86.26±0.08
| S-MLP+kNN-GCN | 81.90±0.06 | 81.51±0.08 | 81.41±0.09 | 80.88±0.11 | 80.28±0.08

Figure 3: Accuracy improvement over GCN. As perturbation rate increases, CoG
(GCN+MLP) outperforms UM-GNN by a large margin.

###  5.2 Performance on clean graphs

The performance of the ensemble models on clean graphs is shown in the 3rd
column of Table 2. From the results, we make the following observations:

  * •

Our ensemble models outperform the baselines on clean graphs among three
datasets. Specifically, the GCN+MLP model outperforms GCN by $0.59\%$,
$3.18\%$, and $0.43\%$ for the three datasets respectively. For the other
three defense models, the best results are $-0.52\%$, $1.28\%$, and $0.07\%$,
respectively.

  * •

Our ensemble models perform better than any single sub-model from the
ensemble. For instance, the ensemble of GCN and kNN-GCN achieve an $84.27\%$
accuracy on Cora, while the accuracy is $83.50\%$ for GCN and $71.06\%$ for
kNN-GCN.

Results show that our method achieves competitive performance on clean data,
thus making it applicable in realistic settings where we have no idea if the
graphs are perturbed.

###  5.3 Performance under attacks

Results under attack are shown in the 4th to 7th columns in Table 2. As we can
see, the co-training framework effectively enhances the model’s robustness
against adversarial attacks. For example, the accuracy of GCN decreases
drastically from $83.5\%$ to $59.6\%$ as the perturbation rate increases from
0 to $20\%$ on Cora. In comparison, the ensemble of GCN and kNN-GCN still
achieves accuracy of $76.86\%$ even in the worst-perturbed case. Similar
results are obtained on the other two datasets.

Compared to the state-of-the-art defenses, our method improves the robustness
among three datasets, especially on Cora and Citeseer. Pro-GNN [15] is the
most robust model among the baselines. Therefore, we mainly compare our
results with those of Pro-GNN. Our ensemble models outperform Pro-GNN by a
large margin on Cora and Citeseer. Specifically, the combination of S-MLP and
kNN-GCN gains $0.61\%$, $2.01\%$, $4.49\%$, and $7.43\%$ improvements as
perturbation rate increases on Cora. The improvements on Citeseer are
$1.86\%$, $2.65\%$, $2.02\%$ and $2.78\%$, respectively. On Pubmed, we also
achieved comparable results. Moreover, it is worth noting that Pro-GNN trains
slow (over 100x longer than GCN), and requires lots of GPU memory (more than
16GB for Pubmed). Our approach is much faster (about 4x longer than GCN) and
less memory-demanding (same as GCN), making it more feasible in practical use.
We also compare CoG with SimP-GCN [14], which also focuses on exploiting the
feature information. Results show that CoG achieves better performance among
all three graphs and perturbation rates. Furthermore, the performance of SimP-
GCN possesses high variance, especially on Cora. One possible reason is that
SimP-GCN has unstable graph structure in its training stage.

UM-GNN [26] uses different data splits from the above methods. It follows the
original split in [17] and uses perturbation rate from $0\%$ to $10\%$.
Results in Fig. 3 show the accuracy improvement of CoG and UM-GNN over GCN. It
shows that although UM-GNN performs slightly better than CoG when the
perturbation rate is small, it is outperformed by a large margin as the
perturbation rate increases. This is caused by the one-directional knowledge
transferring scheme in UM-GNN. When the predictions of GNN is highly
inaccurate, it will mislead the MLP model. This defect is avoided in CoG since
CoG can transfer useful information in both direction. The performance of both
GNN and MLP improves during the co-training process.

###  5.4 Ablation Study

In this section, we conduct an ablation study on the model calibration and
label balancing techniques. We also evaluate how the number of iterations and
the number of nodes to add in each iteration affect the performance.

Figure 4: Confusion matrix and accuracy of the ensemble model’s prediction
during the co-training process (GCN+MLP with Cora). Without class balancing,
the co-training process will overfit to the dominant class, thus leading to
the decrease of accuracy as co-training continues. (Results on Citeseer can be
accessed in the supplementary material.) Figure 5: Reliability Diagram of GCN
w/o (left) and with (right) calibration. After calibrating, model’s output can
better reflect its ground truth correctness likelihood. Figure 6: Node
classification accuracy using different hyper-parameters. Different colors
denote the results under different perturbation rate. (Results on Citeseer can
be accessed in the supplementary material)

####  5.4.1 Model calibration and class balancing

Fig. 5 shows the outputs of GCN with/without model calibration. The
reliability diagram demonstrates that the vanilla GCN is under-confident.
After temperature scaling, model’s output can better reflect its the
correctness likelihood of ground truth more accurately. To validate the
benefit of the class balancing technique, we train an ensemble of GCN+MLP
with/without class balancing on Cora and present the results in Fig. 4. We
report the confusion matrix and the accuracy of the ensemble model in each
iteration. Results show that without class balancing, the co-training process
will overfit to the dominant class. As the co-training process continues, more
and more test data are labeled as the dominant class, thus impairing the
performance of the ensemble model.

####  5.4.2 Number of iterations and added nodes

There are two hyper-parameters to set in CoG, which are the number of
iterations and the number of nodes to add in each iteration. In our
experiments, we add 100, 250, 500 nodes in each iteration respectively and
evaluate CoG’s performance form 0 iteration (ensemble without co-training) to
max iterations (add until no test data left). Results are shown in Fig. 6. As
the co-training process continues, the performance of CoG improves quickly in
the beginning and stabilizes in the later iterations. Adding less nodes per
iteration can slightly improve the performance of CoG. However, it also takes
longer to train CoG, which means a trade-off between effectiveness and
efficiency. Furthermore, comparing to ensemble sub-models without co-training
(the starting points), CoG can improve the performance by a large margin,
especially on the perturbed data.

Figure 7: Node classification accuracy of GCN, MLP and CoG under different
ratios of feature and structure perturbations. (Results on Citeseer can be
accessed in the supplementary material)

###  5.5 Adaptive Attacks

Evaluating CoG against adaptive attacks is necessary since attackers may
attack CoG with both structure and feature perturbations in practice. However,
adaptive attacks on CoG is non-trivial since the co-training process works in
a non-differential manner. To evaluate CoG against adaptive attackers, we
adapt the Metattack to attack the MLP model, keep the total perturbation
budget as $20\%$, and evaluate CoG’s performance under different ratios
between feature and structure perturbations. Results are shown in Fig. 7. We
observe that CoG still achieves good robustness when both feature and
structure information are perturbed.

##  6 Conclusion

In this paper, we present a co-training framework, named CoG, to learn an
ensemble model that integrates the feature and the structure information of
graph data. The orthogonal nature of these two views diversifies the outputs
of sub-models and weakens the transferability of adversarial attacks between
them. Experiment results validate the effectiveness of CoG on both clean and
perturbed graph. In addition, CoG can still achieve good robustness under the
adaptive attack setting.

##  Appendix 0.A Algorithm

Here we demonstrate the overall algorithm of CoG in Algorithm 1. The co-
training process will stop once it reaches the preset number of iterations or
there is no unlabeled data left. In the inference phase, we average each sub-
model’s outputs to obtain the predictions of their ensemble.

Input: graph data $G=(X,A)$, nodes to add per iter $N^{add}$, Labeled dataset
$S$, Unlabel dataset $U$, number of classes $C$

Output: $f_{struct}$, $f_{feat}$

Choose a structure-dominant model $f_{struct}$;

Choose a feature-dominant model $f_{feat}$;

while  _stopping condition not met_ do

Train $f_{struct}$ under structure view of $S$;

Train $f_{feat}$ under feature view of $S$;

Calculate confidence score $s_{struct}$ on $U$;

Calculate confidence score $s_{feat}$ on $U$;

for  _$c=1$ to $C$_ do

Use $f_{struct}$ to label $N^{add}_{c}$ most confident nodes in $U$;

Use $f_{feat}$ to label $N^{add}_{c}$ most confident nodes in $U$;

Conducting temperature scaling for both sub-models;

if  _the same node is chosen_ then

Using the pseudo-label with higher confidence

endif

endfor for

Update $S\&U$;

end while

Algorithm 1 Two-view co-training for graph data

##  Appendix 0.B Ensemble of structure-dominant models

In the main paper, we show that a simple ensemble of GCN [17]+GAT [30] still
suffers from weak robustness against Metattack [38]. Fig.8 shows the
performance of GCN+GAT, GCN+S_MLP and GCN+MLP on Cora and Citeseer. From the
results we can conclude that the diversity of sub-models is crucial for the
robustness of their ensemble.

Figure 8: Performance of ensemble of GCN+GAT, GCN+S_MLP, GCN+MLP (from left to
right) under Metattack. Both the naive ensemble of GCN+GAT and the co-training
of GCN and S_MLP have little help for the robustness, compared to the co-
traning of GCN+MLP.

##  Appendix 0.C Robustness of CoG against PGD

In the main paper, we mainly report CoG’s performance against Metattack. Here
we report its performance (using GCN+MLP) against PGD [34] in Fig. 9. As the
results show, CoG also have good robustness under PGD attack.

Figure 9: The performance of CoG (GCN+MLP) under PGD attack (using the same
data splits with UM-GNN).

##  Appendix 0.D Comparing CoG and UM-GNN under PGD attack

In main paper, we compare CoG and UM-GNN [26] under Metattack. Here we report
the comparison under PGD attack in Fig. 10. Since the perturbed graph used in
UM-GNN is not publicly available, to make a fair comparison, we report the
performance improvement over GCN. For the comparison of absolute accuracy,
please refer to the results in Section 0.C and [26].

Figure 10: Accuracy improvement over GCN under PGD attack. Results keep
consistent with those under Metattack. As perturbation rate increases, CoG
(GCN+MLP) outperforms UM-GNN by a large margin.

##  Appendix 0.E Model calibration

In the main paper, we propose that model calibration is important for CoG. A
calibrated classification model means that its softmax outputs (confidence)
should accurately measure the likelihood of correct predictions. In the main
paper, we use the reliability diagram to demonstrate the miscalibration of
GCN. Reliability diagram is a visual representation of model calibration [7,
22, 11]. It splits the test data according to their output confidence, samples
from them and plots the expected sample accuracy as a function of confidence.
Here in Fig. 11, we show the performance improvement with model calibration.
Using temperature scaling to calibrate the model brings a $1\sim 2\%$
performance improvement.

Figure 11: Performance of CoG (GCN+MLP on Cora) with/without calibration. With
model calibration, CoG performs better, no matter the graph is perturbed or
not.

##  Appendix 0.F Class balancing

In the main paper, we compare CoG’s performance with/without classing
balancing on Cora. Results show that without class balancing, the ensemble
model will overfit to the dominant class. Here we report the confusion matrix
and ensemble accuracy on Citeseer in Fig. 12 and Fig. 13. Confusion matrix is
a specific layout for visualizing the performance of model’s predictions. Each
row of the matrix represents the instances in an actual class while each
column represents the instances in a predicted class. Instances in the
diagonal of the matrix are correctly classified. For more details about the
confusion matrix, we refer the reviewers to [24]. As we can see in Fig. 12,
class 5 is the dominant class in this dataset. Without class balancing, more
and more instances from other classes are misclassified as class 5 as the co-
training process continues, leading to an overfitting to class 5 and a
decrease in accuracy. We obtain the same results on the perturbed data. From
the confusion matrix in Fig. 13, we can also see how the co-training process
(with class balancing) corrects the misclassified instances.

Figure 12: Confusion matrix and accuracy of the ensemble model’s prediction
during the co-training process (GCN+MLP with unperturbed Citeseer data).
Figure 13: Confusion matrix and accuracy of the ensemble model’s prediction
during the co-training process (GCN+MLP with 20% perturbed Citeseer data).

##  Appendix 0.G Hyper-parameters

In main paper, we compare CoG’s performance with different hyper-parameters
(nodes to add per iteration & number of iterations) on Cora. Here we report
the results on Citeseer in Fig. 14. Results are consistent with those on Cora:
with less nodes to add per iterations, the final accuracy is slightly higher.
Note that it is also more time-demanding. As for co-training iterations,
performance improves quickly in the beginning and stabilizes in the final
stage.

Figure 14: Node classification accuracy using different hyper-parameters.
Generally, adding less nodes per iteration leads to better accuracy, although
it is more time-demanding.

##  Appendix 0.H Adaptive attack

In main paper, we demonstrate CoG’s performance under adaptive attacks on
Cora. Here we report the results on Citeseer in Fig. 15. CoG can still achieve
good robustness when the node features and graph structures are perturbed
simultaneously.

Figure 15: Node classification accuracy of GCN, MLP and CoG under different
ratios of feature and structure perturbations. The total perturbation budget
is 20%.

## References

  * [1] Blum, A., Mitchell, T.: Combining labeled and unlabeled data with co-training. In: Proceedings of the eleventh annual conference on Computational learning theory. pp. 92–100 (1998)
  * [2] Bojchevski, A., Günnemann, S.: Adversarial attacks on node embeddings via graph poisoning. In: International Conference on Machine Learning. pp. 695–704. PMLR (2019)
  * [3] Breiman, L.: Bagging predictors. Machine learning 24(2), 123–140 (1996)
  * [4] Chen, L., Li, X., Wu, D.: Enhancing robustness of graph convolutional networks via dropping graph connections. In: Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2020, Ghent, Belgium, September 14–18, 2020, Proceedings, Part III. pp. 412–428. Springer International Publishing (2021)
  * [5] Chung, F.R., Graham, F.C.: Spectral graph theory. No. 92, American Mathematical Soc. (1997)
  * [6] Dai, H., Li, H., Tian, T., Huang, X., Wang, L., Zhu, J., Song, L.: Adversarial attack on graph structured data. In: International conference on machine learning. pp. 1115–1124. PMLR (2018)
  * [7] DeGroot, M.H., Fienberg, S.E.: The comparison and evaluation of forecasters. Journal of the Royal Statistical Society: Series D (The Statistician) 32(1-2), 12–22 (1983)
  * [8] Dietterich, T.G.: Ensemble methods in machine learning. In: International workshop on multiple classifier systems. pp. 1–15. Springer (2000)
  * [9] Entezari, N., Al-Sayouri, S.A., Darvishzadeh, A., Papalexakis, E.E.: All you need is low (rank): Defending against adversarial attacks on graphs. In: Caverlee, J., Hu, X.B., Lalmas, M., Wang, W. (eds.) WSDM ’20: The Thirteenth ACM International Conference on Web Search and Data Mining, Houston, TX, USA, February 3-7, 2020. pp. 169–177. ACM (2020). https://doi.org/10.1145/3336191.3371789, https://doi.org/10.1145/3336191.3371789
  * [10] Fey, M., Lenssen, J.E.: Fast Graph Representation Learning with PyTorch Geometric (1), 1–9 (2019), http://arxiv.org/abs/1903.02428
  * [11] Guo, C., Pleiss, G., Sun, Y., Weinberger, K.Q.: On calibration of modern neural networks. In: International Conference on Machine Learning. pp. 1321–1330. PMLR (2017)
  * [12] Hamilton, W.L., Ying, Z., Leskovec, J.: Inductive representation learning on large graphs. In: Guyon, I., von Luxburg, U., Bengio, S., Wallach, H.M., Fergus, R., Vishwanathan, S.V.N., Garnett, R. (eds.) Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA. pp. 1024–1034 (2017), https://proceedings.neurips.cc/paper/2017/hash/5dd9db5e033da9c6fb5ba83c7a7ebea9-Abstract.html
  * [13] Hansen, L.K., Salamon, P.: Neural Network Ensembles. IEEE Transactions on Pattern Analysis and Machine Intelligence (1990). https://doi.org/10.1109/34.58871
  * [14] Jin, W., Derr, T., Wang, Y., Ma, Y., Liu, Z., Tang, J.: Node Similarity Preserving Graph Convolutional Networks (2020), http://arxiv.org/abs/2011.09643
  * [15] Jin, W., Ma, Y., Liu, X., Tang, X., Wang, S., Tang, J.: Graph structure learning for robust graph neural networks. In: Gupta, R., Liu, Y., Tang, J., Prakash, B.A. (eds.) KDD ’20: The 26th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, Virtual Event, CA, USA, August 23-27, 2020\. pp. 66–74. ACM (2020), https://dl.acm.org/doi/10.1145/3394486.3403049
  * [16] Kariyappa, S., Qureshi, M.K.: Improving adversarial robustness of ensembles with diversity training. arXiv (2019)
  * [17] Kipf, T.N., Welling, M.: Semi-supervised classification with graph convolutional networks. In: 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net (2017), https://openreview.net/forum?id=SJU4ayYgl
  * [18] Kuncheva, L.I., Whitaker, C.J.: Measures of diversity in classifier ensembles and their relationship with the ensemble accuracy. Machine Learning (2003). https://doi.org/10.1023/A:1022859003006
  * [19] Li, Y., Jin, W., Xu, H., Tang, J.: Deeprobust: A pytorch library for adversarial attacks and defenses (2020)
  * [20] Luo, D., Cheng, W., Yu, W., Zong, B., Ni, J., Chen, H., Zhang, X.: Learning to drop: Robust graph neural network via topological denoising. arXiv (2020)
  * [21] McCallum, A.K., Nigam, K., Rennie, J., Seymore, K.: Automating the construction of internet portals with machine learning. Information Retrieval 3(2), 127–163 (2000)
  * [22] Niculescu-Mizil, A., Caruana, R.: Predicting good probabilities with supervised learning. In: Proceedings of the 22nd international conference on Machine learning. pp. 625–632 (2005)
  * [23] Pang, T., Xu, K., Du, C., Chen, N., Zhu, J.: Improving adversarial robustness via promoting ensemble diversity. In: Chaudhuri, K., Salakhutdinov, R. (eds.) Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA. Proceedings of Machine Learning Research, vol. 97, pp. 4970–4979. PMLR (2019), http://proceedings.mlr.press/v97/pang19a.html
  * [24] Powers, D.M.: Evaluation: from precision, recall and f-measure to roc, informedness, markedness and correlation. arXiv preprint arXiv:2010.16061 (2020)
  * [25] Sen, P., Namata, G., Bilgic, M., Getoor, L., Galligher, B., Eliassi-Rad, T.: Collective classification in network data. AI magazine 29(3), 93–93 (2008)
  * [26] Shanthamallu, U.S., Thiagarajan, J.J., Spanias, A.: Uncertainty-matching graph neural networks to defend against poisoning attacks. In: Proceedings of the AAAI Conference on Artificial Intelligence. vol. 35, pp. 9524–9532 (2021)
  * [27] Stellargraph: Ensemble models for node classification (2020), https://stellargraph.readthedocs.io/en/stable/demos/ensembles/ensemble-node-classification-example.html
  * [28] Teixeira, L., Jalaian, B., Ribeiro, B.: Are graph neural networks miscalibrated? arXiv preprint arXiv:1905.02296 (2019)
  * [29] Tramèr, F., Papernot, N., Goodfellow, I., Boneh, D., McDaniel, P.: The space of transferable adversarial examples. arXiv preprint arXiv:1704.03453 (2017)
  * [30] Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., Bengio, Y.: Graph attention networks. In: 6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings. OpenReview.net (2018), https://openreview.net/forum?id=rJXMpikCZ
  * [31] Von Luxburg, U.: A tutorial on spectral clustering. Statistics and computing 17(4), 395–416 (2007)
  * [32] Wang, W., Zhou, Z.H.: Analyzing co-training style algorithms. In: European conference on machine learning. pp. 454–465. Springer (2007)
  * [33] Wu, H., Wang, C., Tyshetskiy, Y., Docherty, A., Lu, K., Zhu, L.: Adversarial examples for graph data: Deep insights into attack and defense. In: Kraus, S. (ed.) Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI 2019, Macao, China, August 10-16, 2019. pp. 4816–4823. ijcai.org (2019). https://doi.org/10.24963/ijcai.2019/669, https://doi.org/10.24963/ijcai.2019/669
  * [34] Xu, K., Chen, H., Liu, S., Chen, P.Y., Weng, T.W., Hong, M., Lin, X.: Topology attack and defense for graph neural networks: An optimization perspective. arXiv preprint arXiv:1906.04214 (2019)
  * [35] Yang, H., Zhang, J., Dong, H., Inkawhich, N., Gardner, A., Touchet, A., Wilkes, W., Berry, H., Li, H.: DVERGE: Diversifying vulnerabilities for enhanced robust generation of ensembles. arXiv (NeurIPS), 1–20 (2020)
  * [36] Zhang, X., Zitnik, M.: GNNGuard: Defending Graph Neural Networks against Adversarial Attacks (2020), http://arxiv.org/abs/2006.08149
  * [37] Zügner, D., Akbarnejad, A., Günnemann, S.: Adversarial attacks on neural networks for graph data. In: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. pp. 2847–2856 (2018)
  * [38] Zügner, D., Günnemann, S.: Adversarial attacks on graph neural networks via meta learning. In: 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net (2019), https://openreview.net/forum?id=Bylnx209YX