# Topological Superconductivity and Superfluidity

Xiao-Liang Qi, Taylor L. Hughes, Srinivas Raghu and Shou-Cheng Zhang
Department of Physics, McCullough Building, Stanford University, Stanford, CA
94305-4045

###### Abstract

We construct time reversal invariant topological superconductors and
superfluids in two and three dimensions which are analogous to the recently
discovered quantum spin Hall and three-d $Z_{2}$ topological insulators
respectively. These states have a full pairing gap in the bulk, gapless
counter-propagating Majorana states at the boundary, and a pair of Majorana
zero modes associated with each vortex. We show that the time reversal
symmetry naturally emerges as a supersymmetry, which changes the parity of the
fermion number associated with each time-reversal invariant vortex. In the
presence of external T-breaking fields, non-local topological correlation is
established among these fields, which is an experimentally observable
manifestation of the emergent supersymmetry.

###### pacs:

74.20.Rp, 73.43.-f, 67.30.he, 74.45.+c

The search for topological states of quantum matter has become an active and
exciting pursuit in condensed matter physics. The quantum Hall (QH) effect
Prange and Girvin (1990) provides the first example of a topologically non-
trivial state of matter, where the quantized Hall conductance is a topological
invariantThouless et al. (1982). Recently, the quantum spin Hall (QSH) stateC.
L. Kane and E. J. Mele (2005a); B.A. Bernevig and S.C. Zhang (2006) has been
theoretically predicted B. A. Bernevig et al. (2006) and experimentally
observed in HgTe quantum well systemsKönig et al. (2007). The time reversal
invariant (TRI) QSH state is characterized by a bulk gap, a $Z_{2}$
topological number C. L. Kane and E. J. Mele (2005b), and gapless helical edge
states, where time-reversed partners counter-propagateC. Wu et al. (2006); C.
Xu and J. Moore (2006).

Chiral superconductors in a time reversal symmetry breaking (TRB)
$(p_{x}+ip_{y})$ pairing state in 2d have a sharp topological distinction
between the strong and weak pairing regimes Read and Green (2000). In the weak
pairing regime, the system has a full bulk gap and gapless chiral Majorana
states at the edge, which are topologically protected. Moreover, a Majorana
zero mode is trapped in each vortex coreRead and Green (2000), which leads to
a ground state degeneracy of $2^{n-1}$ in the presence of $2n$ vortices. When
the vortices wind around each other a non-Abelian Berry phase is generated in
the $2^{n-1}$ dimensional ground state manifold, which implies non-Abelian
statistics for the vorticesIvanov (2001). Chiral superconductors are analogous
to the QH state— they both break time reversal (TR) and have chiral edge
states with linear dispersion. However, the edge states of a chiral
superconductor have only half the degrees of freedom compared to the QH state,
since the negative energy quasi-particle operators on the edge of a chiral
superconductor describe the same excitations as the positive energy ones.

Given the analogy between the chiral superconducting state and the QH state,
and with the recent discovery of the TRI QSH state, it is natural to
generalize the chiral pairing state to the helical pairing state, where
fermions with up spins are paired in the $(p_{x}+ip_{y})$ state, while
fermions with down spins are paired in the $(p_{x}-ip_{y})$ state. Such a TRI
state have a full gap in the bulk, and counter-propagating helical Majorana
states at the edge (in contrast, the edge states of the TRI topological
insulator are helical Dirac fermions). Just as in the case of the QSH state, a
mass term for an odd number of pairs of helical Majorana states is forbidden
by TR symmetry, and therefore, a topologically protected superconducting or
superfluid state can exist in the presence of time-reversal symmetry.
Recently, a $Z_{2}$ classification of the topological superconductor has been
discussed in Refs Roy (a, 2008); Schnyder et al. (2008), by noting the
similarity between the Bogoliubov-de Gennes (BdG) superconductor Hamiltonian
and the QSH insulator Hamiltonian. The four types of topological states of
matter discussed here are summarized in Fig. 1. In this work, we give a
$Z_{2}$ classification of both the 2D and 3D cases which has a profound
physical implication. In two dimensions, we show that a time-reversal
invariant topological defect of a $Z_{2}$ non-trivial superconductor carries a
Kramers’ pair of Majorana fermions. Let $N_{F}$ be the operator which measures
the number of fermions of a general system, then the fermion-number parity
operator is given by $(-1)^{N_{F}}$. This operator is also referred to as the
Witten indexWitten (1982), which plays a crucial role in supersymmetric
theories. We prove the remarkable fact that in the presence of a topological
defect, the TR operator $\cal T$ changes the fermion number parity, ${\cal
T}^{-1}(-1)^{N_{F}}{\cal T}=-(-1)^{N_{F}}$ locally around the defect in the
$Z_{2}$ non-trivial state, while it preserves the fermion number parity,
${\cal T}^{-1}(-1)^{N_{F}}{\cal T}=(-1)^{N_{F}}$, in the $Z_{2}$ trivial
state. This fact gives a precise definition of the $Z_{2}$ topological
classification of any TRI superconductor state and is generally valid in the
presence of interactions and disorder. A supersymmetric operation can be
defined as an operation which changes the fermion number parity; therefore, in
this precise sense, we show that the TR symmetry emerges as a supersymmetry in
topological superconductors. Though supersymmetry has been studied extensively
in high energy physics, it has not yet been observed in Nature. Our proposal
offers the opportunity to experimentally observe supersymmetry in condensed
matter systems without any fine tuning of microscopic parameters. The physical
consequences of such a supersymmetry is also studied.

Figure 1: (Top row) Schematic comparison of $2d$ chiral superconductor and the
QH state. In both systems, TR symmetry is broken and the edge states carry a
definite chirality. (Bottom row) Schematic comparison of $2d$ TRI topological
superconductor and the QSH insulator. Both systems preserve TR symmetry and
have a helical pair of edge states, where opposite spin states counter-
propagate. The dashed lines show that the edge states of the superconductors
are Majorana fermions so that the $E<0$ part of the quasi-particle spectra are
redundant. In terms of the edge state degrees of freedom, we have
${\rm(QSH)}={\rm(QH)}^{2}={\rm(Helical~{}SC)}^{2}={\rm(Chiral~{}SC)}^{4}$.

As the starting point, we consider a TRI $p$-wave superconductor with spin
triplet pairing, which has the following $4\times 4$ BdG Hamiltonian:

$\displaystyle H=\frac{1}{2}\int
d^{2}x\Psi^{\dagger}(x)\left(\begin{array}[]{cc}\epsilon_{\bf
p}\mathbb{I}&i\sigma_{2}\sigma_{\alpha}\Delta^{\alpha j}p_{j}\\\
h.c.&-\epsilon_{\bf p}\mathbb{I}\end{array}\right)\Psi(x)$ (3)

with
$\Psi(x)=\left(c_{\uparrow}(x),c_{\downarrow}(x),c_{\uparrow}^{\dagger}(x),c_{\downarrow}^{\dagger}(x)\right)^{T}$,
$\epsilon_{\bf p}={\bf p}^{2}/2m-\mu$ the kinetic energy and chemical
potential terms and $h.c.\equiv(i\sigma_{2}\sigma_{\alpha}\Delta^{\alpha
j}p_{j})^{\dagger}$. The TR transformation is defined as
$c_{\uparrow}\rightarrow c_{\downarrow},~{}c_{\downarrow}\rightarrow-
c_{\uparrow}$. It can be shown that the Hamiltonian (3) is time-reversal
invariant if $\Delta_{\alpha j}$ is a real matrix. To show the existence of a
topological state, consider the TRI mean-field ansatz $\Delta^{\alpha
1}=\Delta(1,0,0),~{}\Delta^{\alpha 2}=\Delta(0,1,0)$. For such an ansatz the
Hamiltonian (3) is block diagonal with only equal spin pairing:

$\displaystyle H=\frac{1}{2}\int
d^{2}x\tilde{\Psi}^{\dagger}\left(\begin{array}[]{cccc}\epsilon_{\bf p}&\Delta
p_{+}&&\\\ \Delta p_{-}&-\epsilon_{\bf p}&&\\\ &&\epsilon_{\bf p}&-\Delta
p_{-}\\\ &&-\Delta p_{+}&-\epsilon_{\bf p}\end{array}\right)\tilde{\Psi}$ (8)

with
$\tilde{\Psi}(x)\equiv\left(c_{\uparrow}(x),c_{\uparrow}^{\dagger}(x),c_{\downarrow}(x),c_{\downarrow}^{\dagger}(x)\right)^{T}$,
and $p_{\pm}=p_{x}\pm ip_{y}$. From this Hamiltonian we see that the spin up
(down) electrons form $p_{x}+ip_{y}$ ($p_{x}-ip_{y}$) Cooper pairs,
respectively. In the weak pairing phase with $\mu>0$, the $(p_{x}+ip_{y})$
chiral superconductor is known to have chiral Majorana edge states propagating
on each boundary, described by the Hamiltonian $H_{\rm edge}=\sum_{k_{y}\geq
0}v_{F}k_{y}\psi_{-k_{y}}\psi_{k_{y}}$,where
$\psi_{-k_{y}}=\psi_{k_{y}}^{\dagger}$ is the quasiparticle creation operator
Read and Green (2000) and the boundary is taken to be parallel to the $y$
direction. Thus we know that the edge states of the TRI system described by
Hamiltonian (8) consist of spin up and spin down quasi-particles with opposite
chirality:

$\displaystyle H_{\rm edge}=\sum_{k_{y}\geq
0}v_{F}k_{y}\left(\psi_{-k_{y}\uparrow}\psi_{k_{y}\uparrow}-\psi_{-k_{y}\downarrow}\psi_{k_{y}\downarrow}\right).$
(9)

The quasi-particle operators $\psi_{k_{y}\uparrow},~{}\psi_{k_{y}\downarrow}$
can be expressed in terms of the eigenstates of the BdG Hamiltonian as

$\displaystyle\psi_{k_{y}\uparrow}$ $\displaystyle=$ $\displaystyle\int
d^{2}x\left(u_{k_{y}}(x)c_{\uparrow}(x)+v_{k_{y}}(x)c_{\uparrow}^{\dagger}(x)\right)$
$\displaystyle\psi_{k_{y}\downarrow}$ $\displaystyle=$ $\displaystyle\int
d^{2}x\left(u_{-k_{y}}^{*}(x)c_{\downarrow}(x)+v_{-k_{y}}^{*}(x)c_{\downarrow}^{\dagger}(x)\right)$
(10)

from which the time-reversal transformation of the quasiparticle operators can
be determined to be
${\cal{T}}^{-1}\psi_{k_{y}\uparrow}{\cal{T}}=\psi_{-k_{y}\downarrow},~{}{\cal{T}}^{-1}\psi_{k_{y}\downarrow}{\cal{T}}=-\psi_{-k_{y}\uparrow}$.
In other words, $(\psi_{k_{y}\uparrow},~{}\psi_{-k_{y}\downarrow})$ transforms
as a Kramers’ doublet, which forbids a gap in the edge states due to mixing of
the spin-up and spin-down modes when TR is preserved. To see this explicitly,
notice that the only $k_{y}$-independent term that can be added to the edge
Hamiltonian (9) is $im\sum_{k_{y}}\psi_{-k_{y}\uparrow}\psi_{k_{y}\downarrow}$
with $m\in\mathbb{R}$. However, such a term is odd under TR, which implies
that any back scattering between the quasi-particles is forbidden by TR
symmetry. The discussion above is exactly parallel to the $Z_{2}$ topological
characterization of the quantum spin Hall system. In fact, the Hamiltonian (8)
has exactly the same form as the four band effective Hamiltonian proposed in
Ref.B. A. Bernevig et al. (2006) to describe HgTe quantum wells with the QSH
effect. The edge states of the QSH insulators consist of an odd number of
Kramers’ pairs, which remain gapless under any small TR-invariant
perturbationC. Wu et al. (2006); C. Xu and J. Moore (2006). A no-go theorem
states that such a “helical liquid” with an odd number of Kramers’ pairs at
the Fermi energy can not be realized in any bulk 1d system, but can only
appear as an edge theory of a 2d QSH insulatorC. Wu et al. (2006). Similarly,
the edge state theory (9) can be called a “helical Majorana liquid”, which can
only exist on the boundary of a $Z_{2}$ topological superconductor. Once such
a topological phase is established, it is robust under any TRI perturbations.

The Hamiltonian (3) can be easily generalized to three dimensions, in which
case $\Delta^{\alpha j}$ becomes a $3\times 3$ matrix with $\alpha=1,2,3$ and
$j=x,y,z$. An example of such a Hamiltonian is given by the well-known 3He BW
phase, for which the order parameter $\Delta^{\alpha j}$ is determined by an
orthogonal matrix $\Delta^{\alpha j}=\Delta u^{\alpha j}$, $u\in{\rm
SO(3)}$Vollhardt and Wölfle (1990). Here and below we ignore the dipole-dipole
interaction term Leggett (1975) since it does not affect any essential
topological properties. By applying a spin rotation, $\Delta^{\alpha j}$ can
be diagonalized to $\Delta^{\alpha j}=\Delta\delta^{\alpha j}$, in which case
the Hamiltonian (3) has the same form as a 3d Dirac Hamiltonian with momentum
dependent mass $\epsilon({\bf p})={\bf p}^{2}/2m-\mu$. We know that a band
insulator described by the Dirac Hamiltonian is a 3d $Z_{2}$ topological
insulator for $\mu>0$Fu et al. (2007); Moore and Balents (2007); Roy (b), and
has nontrivial surface states. The corresponding superconductor Hamiltonian
describes a topological superconductor with 2d gapless Majorana surface
states. The surface theory can be written as

$\displaystyle H_{\rm surf}=\frac{1}{2}\sum_{\bf k}v_{F}\psi_{-\bf
k}^{T}\left(\sigma_{z}k_{x}+\sigma_{x}k_{y}\right)\psi_{\bf k}$ (11)

which remains gapless under any small TRI perturbation since the only
available mass term $m\sum_{\bf k}\psi_{-\bf k}^{T}\sigma_{y}\psi_{\bf k}$ is
time-reversal odd. We would like to mention that the surface Andreev bound
states in 3He-B phase have been observed experimentallyAoki et al. (2005).

To understand the physical consequences of the nontrivial topology we study
the TRI topological defects of the topological superconductors. We start by
considering the equal-spin pairing system with BdG Hamiltonian (8) in which
spin up and down electrons form $p_{x}+ip_{y}$ and $p_{x}-ip_{y}$ Cooper
pairs, respectively. A TRI topological defect can be defined as a vortex of
spin-up superfluid coexisting with an anti-vortex of spin-down superfluid at
the same position. In the generic Hamiltonian (3), such a vortex configuration
is written as $\Delta^{\alpha j}=[\exp\left(i\sigma_{2}\theta({\bf
r-r_{0}})\right)]^{\alpha j},~{}\alpha=1,2$ and $\Delta^{3j}=0$, where
$\theta({\bf r-r_{0}})$ is the angle of ${\bf r}$ with respect to the vortex
position ${\bf r}_{0}$. Since in the vortex core of a weak pairing
$p_{x}+ip_{y}$ superconductor there is a single Majorana zero modeRead and
Green (2000); Stone and Chung (2006), one immediately knows that a pair a
Majorana zero modes exist in the vortex core we study here. In terms of the
electron operators, the two Majorana fermion operators can be written as

$\displaystyle\gamma_{\uparrow}$ $\displaystyle=$ $\displaystyle\int
d^{2}x\left(u_{0}(x)c_{\uparrow}(x)+u_{0}^{*}(x)c_{\uparrow}^{\dagger}(x)\right)$
$\displaystyle\gamma_{\downarrow}$ $\displaystyle=$ $\displaystyle\int
d^{2}x\left(u_{0}^{*}(x)c_{\uparrow}(x)+u_{0}(x)c_{\uparrow}^{\dagger}(x)\right)$
(12)

where we have used the fact that the spin-down zero mode wave function can be
obtained from the time-reversal transformation of the spin-up one. The
Majorana operators satisfy the anti-commutation relation
$\left\\{\gamma_{\alpha},\gamma_{\beta}\right\\}=2\delta_{\alpha\beta}$. The
TR transformation of the Majorana fermions is

$\displaystyle{\cal{T}}^{-1}\gamma_{\uparrow}{\cal{T}}=\gamma_{\downarrow},~{}{\cal{T}}^{-1}\gamma_{\downarrow}{\cal{T}}=-\gamma_{\uparrow}.$
(13)

Similar to the case of the edge states studied earlier, the Majorana zero
modes are robust under any small TRI perturbation, since the only possible
term $im\gamma_{\uparrow}\gamma_{\downarrow}$ which can lift the zero modes to
finite energy is TR odd, i.e.,
${\cal{T}}^{-1}i\gamma_{\uparrow}\gamma_{\downarrow}{\cal{T}}=-i\gamma_{\uparrow}\gamma_{\downarrow}$.

The properties of such a topological defect appear identical to that of a
$\pi$-flux tube threading into a TRI topological insulatorQi and Zhang ; Ran
et al. , where a Kramers’ pair of complex fermions are trapped by the flux
tube. However, there is an essential difference. From the two Majorana zero
modes $\gamma_{\uparrow},\gamma_{\downarrow}$ a complex fermion operator can
be defined as $a=\left(\gamma_{\uparrow}+i\gamma_{\downarrow}\right)/2$, which
satisfies the fermion anticommutation relation
$\left\\{a,a^{\dagger}\right\\}=1$. Since
$\gamma_{\uparrow},\gamma_{\downarrow}$ are zero modes, we obtain
$\left[a,H\right]=0$ which implies that $a$ is the annihilation operator of a
zero-energy quasiparticle. Consequently, the ground state of the system is at
least two-fold degenerate, with two states $|G_{0}\rangle$ and
$|G_{1}\rangle=a^{\dagger}|G_{0}\rangle$ containing $0$ and $1$ $a$-fermions.
Since $a^{\dagger}a=\left(1+i\gamma_{\uparrow}\gamma_{\downarrow}\right)/2$,
the states $\left|G_{0(1)}\right\rangle$ are eigenstates of
$i\gamma_{\uparrow}\gamma_{\downarrow}$ with eigenvalues $-1(+1)$,
respectively. Thus from the oddness of $i\gamma_{\uparrow}\gamma_{\downarrow}$
under TR we know that $\left|G_{0}\right\rangle$ and
$\left|G_{1}\right\rangle$ are time-reversed partners. Note that
superconductivity breaks the charge $U(1)$ symmetry to $Z_{2}$, meaning that
the fermion number parity operator $(-1)^{N_{F}}$ is conserved. Thus, all the
eigenstates of the Hamiltonian can be classified by the value of
$(-1)^{N_{F}}$. If, say, $\left|G_{0}\right\rangle$ is a state with
$(-1)^{N_{F}}=1$, then
$\left|G_{1}\right\rangle=a^{\dagger}\left|G_{0}\right\rangle$ must satisfy
$(-1)^{N_{F}}=-1$. Since $\left|G_{0}\right\rangle$ and
$\left|G_{1}\right\rangle$ are time-reversal partners, we know that in the
Hilbert space of the zero-energy states the TR transformation changes the
fermion number parity:

$\displaystyle{\cal{T}}^{-1}(-1)^{N_{F}}{\cal{T}}=-(-1)^{N_{F}}.$ (14)

Eq. (14) is the central result of this paper. At a first glance it seems
contradict the fundamental fact that the electron number of the whole system
is invariant under TR. Such a paradox is resolved by noticing that there are
always an even number of topological defects in a closed system without
boundary. Under the TR transformation, the fermion number parity around each
vortex core is odd, but the total fermion number parity remains even as
expected. Once the anomalous transformation property (14) is established for a
topological defect in a TRI superconductor, it is robust under any TRI
perturbation as long as the bulk quasiparticle gap remains finite and other
topological defects are far away. Thus Eq. (14) is a generic definition of TRI
topological superconductors:

  * •

Definition I. A two-dimensional TRI superconductor is $Z_{2}$ nontrivial if
and only if fermion number parity around a TRI topological defect is odd under
TR.

A transformation changing fermion number by an odd number is a
“supersymmetry”; thus, the TR symmetry emerges as a discrete supersymmetry for
each TRI topological defect. The same analysis applies to the edge theory (9),
which shows that in the 1d helical Majorana liquid is a theory with TR
symmetry as a discrete supersymmetry.

All the conclusions above can be generalized to 3d topological
superconductors. In the 3He BW phase the Goldstone manifold of the order
parameter is $\Delta^{\alpha j}=\Delta u^{\alpha j}\in SO(3)\times
U(1)$Vollhardt and Wölfle (1990); Salomaa and Volovik (1987). A time-reversal
invariant configuration satisfies $\Delta^{\alpha j}\in\mathbb{R}$, which
restricts the order parameter to $SO(3)$. Since $\Pi_{1}(SO(3))=Z_{2}$, the
TRI topological defects are 1d “vortex” rings. By solving the BdG equations in
the presence of such vortex rings, it can be shown that there are linearly
dispersing quasiparticles propagating on each vortex ring, similar to the edge
states of the 2d topological superconductor. However, for a ring with finite
length the quasi-particle spectrum is discrete. Specifically, there may or may
not be a pair of Majorana modes at exactly zero energy. The existence of the
Majorana zero modes on the vortex rings turns out to be a topological property
determined by the linking number between different vortex rings. Due to the
length constraints of the present paper, we will write our conclusion and
leave the details for a separate work: There are a pair of Majorana fermion
zero modes confined on a vortex ring if and only if the ring is linked to an
odd number of other vortex rings. Such a condition is shown in Fig. 2.
Consequently, the generalization of Definition I to 3d is:

  * •

Definition II. A 3d TRI superconductor is $Z_{2}$ nontrivial if and only if
the fermion number parity around one of the two mutually-linked TRI vortex
rings is odd under TR.

Figure 2: Illustration of a 3d TRI topological superconductor with two TRI
vortex rings which are (a) linked or (b) unlinked. The $E-k_{\parallel}$
dispersion relations show schematically the quasiparticle levels confined on
the red vortex ring in both cases. “$\circ$” and “$\times$” stand for the
quasiparticle levels that are Kramers’ partners of each other. Only case (a)
has a pair of Majorana zero modes located on each vortex ring. Figure 3: (a)
Illustration of a 2d TRI topological superconductor with four TRI topological
defects coupled to a TR-breaking field. The arrows show the sign of the TR-
breaking field ${\rm sgn}(M({\bf r}_{s}))$ at each topological defect. In the
two configurations shown, only the field around vortex 1 is flipped, leading
to an opposite fermion number parity in the corresponding ground state
$|G\rangle$ and $|G^{\prime}\rangle$ (see text). (b) Illustration showing the
flow of the energy levels when the upper configuration in figure (a) is
deformed to the lower one. The flip of the TR-breaking field $M({\bf r}_{1})$
leads to a level crossing at $M({\bf r}_{1})=0$, where the fermion number
parity in the ground state changes sign.

Besides providing a generic definition of the $Z_{2}$ topological
superconductors, such an emergent supersymmetry also leads to physical
predictions. Consider the 2d topological superconductor coupled to a weak TR-
breaking field $M({\bf r})$, which is classical but can have thermal
fluctuations. This situation can be realized in an isolated superconductor
with vortices pinned to quenched weak magnetic impurities. The $n$-point
correlation function of $M({\bf r})$ can be obtained by

$\displaystyle\left\langle{\prod_{s=1}^{n}M({\bf
r}_{s})}\right\rangle\equiv\int\frac{D[M({\bf r})]}{Z}\prod_{s=1}^{n}M({\bf
r}_{s}){\rm Tr}\left(e^{-\beta H[M]}\right)_{\rm even}$ (15)

in which the trace is restricted to states with an even number of fermions.
For a closed system with $N$ vortices, the leading order effect of the TR-
breaking field is to lift the degeneracy between the two Majorana fermions in
each vortex core. Consequently, the perturbed Hamiltonian $H[M({\bf r})]$ to
first order can be written as

$\displaystyle H[M({\bf r})]=\sum_{s=1}^{N}iM({\bf
r}_{s})a_{s}\gamma_{s\uparrow}\gamma_{s\downarrow}$ (16)

in which $\gamma_{s\uparrow(\downarrow)}$ are the Majorana fermion operators,
and $a_{s}\in\mathbb{R}$ depend on the details of the perturbation. The
important fact is that the mass term induced is linear in $M({\bf r})$ at the
defect position ${\bf r}_{s}$, since $i\gamma_{s\uparrow}\gamma_{s\downarrow}$
is TR odd.

Since the superconductor has a full gap, naively one would expect all the
correlations of $M({\bf r})$ field to be short ranged. However, for a system
with $N$ topological defects the $N$-point correlation function has a long
range order when ${\bf r}_{s},s=1,2,..N$ are chosen to be the coordinates of
the topological defects. In other words, the correlation function

$\displaystyle\lim_{|{\bf r}_{i}-{\bf r}_{j}|\rightarrow\infty,~{}\forall
i,j}\left\langle{\prod_{s=1}^{N}M({\bf r}_{s})}\right\rangle\neq 0,$ (17)

though all the $n$ point correlations in Eq. (15) with $n<N$ remain short
ranged. Physically, such a non-local correlation can be understood by
comparing two states $|G\rangle$ and $|G^{\prime}\rangle$, which are the
ground states of the systems with the field configurations
$\mathcal{M}\equiv(M({\bf r}_{1}),M({\bf r}_{2}),...,M({\bf r}_{N}))$ and
$\mathcal{M}^{\prime}\equiv(-M({\bf r}_{1}),M({\bf r}_{2}),...,M({\bf
r}_{N}))$, respectively. From Hamiltonian (16) it can be seen that $|G\rangle$
and $|G^{\prime}\rangle$ have opposite fermion number parity, since the
fermion number parity around the first topological defect
$i\gamma_{1\uparrow}\gamma_{1\downarrow}$ is reversed while that of all the
other topological defects remains invariant, as shown in Fig. 3. Without loss
of generality, we can assume $(-1)^{N_{F}}$ is even for $|G\rangle$ and odd
for $|G^{\prime}\rangle$. Since the whole system is required to have an even
number of fermions, the lowest energy state in the Hilbert space for the field
configuration $\mathcal{M}^{\prime}$ is not $|G^{\prime}\rangle$, but the
lowest quasiparticle excitation $a_{\rm min}^{\dagger}|G^{\prime}\rangle$.
Thus, the two field configurations $\mathcal{M}$ and $\mathcal{M}^{\prime}$
have different free energies, which leads to the non-vanishing correlation
function in Eq. (17). Even when the topological defects are arbitrarily far
away, the energy difference between the two configurations remains finite,
which shows the non-local topological correlation. Similar non-local
correlations can also be obtained for a 3d TRI topological superconductor with
linked vortex rings.

Acknowledgement.—We acknowledge helpful discussions with S. B. Chung, A. L.
Fetter, L. Fu, R. Roy and S. Ryu. This work is supported by the NSF under
grant numbers DMR-0342832, the US Department of Energy, Office of Basic Energy
Sciences under contract DE-AC03-76SF00515, and the Stanford Institute for
Theretical Physics (S.R.).

## References

  * Prange and Girvin (1990) R. E. Prange and S. M. Girvin, eds., _The Quantum Hall effect_ (Springer-Verlag, USA, 1990).
  * Thouless et al. (1982) D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Phys. Rev. Lett. 49, 405 (1982).
  * C. L. Kane and E. J. Mele (2005a) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005a).
  * B.A. Bernevig and S.C. Zhang (2006) B.A. Bernevig and S.C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
  * B. A. Bernevig et al. (2006) B. A. Bernevig, T. L. Hughes, and S.C. Zhang, Science 314, 1757 (2006).
  * König et al. (2007) M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * C. L. Kane and E. J. Mele (2005b) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005b).
  * C. Xu and J. Moore (2006) C. Xu and J. Moore, Phys. Rev. B 73, 045322 (2006).
  * C. Wu et al. (2006) C. Wu, B.A. Bernevig, and S.C. Zhang, Phys. Rev. Lett. 96, 106401 (2006).
  * Read and Green (2000) N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * Ivanov (2001) D. A. Ivanov, Phys. Rev. Lett. 86, 268 (2001).
  * Roy (a) R. Roy, arxiv: cond-mat/0608064.
  * Roy (2008) R. Roy, arxiv: 0803.2868 (2008).
  * Schnyder et al. (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, e-print arXiv: 0803.2786 (2008).
  * Witten (1982) E. Witten, Nucl. Phys. B 202, 253 (1982).
  * Vollhardt and Wölfle (1990) D. Vollhardt and P. Wölfle, _The Superfluid Phases of Helium 3_ (Taylor and Francis, USA, 1990).
  * Leggett (1975) A. J. Leggett, Rev. Mod. Phys. 47, 331 (1975).
  * Fu et al. (2007) L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
  * Moore and Balents (2007) J. E. Moore and L. Balents, Phys. Rev. B 75, 121306 (2007).
  * Roy (b) R. Roy, arxiv: cond-mat/0607531.
  * Aoki et al. (2005) Y. Aoki, Y. Wada, M. Saitoh, R. Nomura, Y. Okuda, Y. Nagato, M. Yamamoto, S. Higashitani, and K. Nagai, Phys. Rev. Lett. 95, 075301 (2005).
  * Stone and Chung (2006) M. Stone and S.-B. Chung, Phys. Rev. B 73, 014505 (2006).
  * (23) X. L. Qi and S. C. Zhang, arxiv: cond-mat/0801.0252.
  * (24) Y. Ran, A. Vishwanath, and D.-H. Lee, arxiv: cond-mat/0801.0627.
  * Salomaa and Volovik (1987) M. M. Salomaa and G. E. Volovik, Rev. Mod. Phys. 59, 533 (1987).


#  A disk instability model for the quasi-periodic eruptions of GSN 069

Xin Pan Key Laboratory for Research in Galaxies and Cosmology, Shanghai
Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan Road,
Shanghai 200030, People’s Republic of China  University of Chinese Academy of
Sciences, 19A Yuquan Road, 100049, Beijing, People’s Republic of China
Shuang-Liang Li Key Laboratory for Research in Galaxies and Cosmology,
Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80 Nandan
Road, Shanghai 200030, People’s Republic of China  University of Chinese
Academy of Sciences, 19A Yuquan Road, 100049, Beijing, People’s Republic of
China  Xinwu Cao Zhejiang Institute of Modern Physics, Department of Physics,
Zhejiang University, 38 Zheda Road, Hangzhou 310027, People’s Republic of
China  Shanghai Astronomical Observatory, Chinese Academy of Sciences, 80
Nandan Road, Shanghai, 200030, People’s Republic of China  Giovanni Miniutti
Centro de Astrobiología (CSIC-INTA), Camino Bajo del Castillo s/n, Villanueva
de la Cañada, E-28692 Madrid, Spain  Minfeng Gu Key Laboratory for Research in
Galaxies and Cosmology, Shanghai Astronomical Observatory, Chinese Academy of
Sciences, 80 Nandan Road, Shanghai 200030, People’s Republic of China Xin Pan,
Shuang-Liang Li, Xinwu Cao<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

GSN 069 is a recently discovered QPE (Quasi-periodic eruptions) source
recurring about every 9 hours. The mechanism for the QPEs of GSN 069 is still
unclear so far. In this work, a disk instability model is constructed to
explain GSN 069 based on Pan et al. (2021) (PLC21), where the authors proposed
a toy model for the repeating changing-look (CL) active galactic nuclei (AGN).
We improve the work of PLC21 by including a non-zero viscous torque condition
on the inner boundary of disk and adopting a general form for the viscous
stress torque in Kerr metric. It is found that the 0.4-2 keV light curves, the
light curves at different energy bands and the phase-resolved X-ray spectrum
of GSN 069 can all be qualitatively reproduced by our model. Furthermore, the
profiles of light curve in QPEs can be significantly changed by the parameter
$\mu$ in viscous torque equation, which implies that our model may also be
applied to other QPEs.

accretion disks — instabilities: active — galaxies: Seyfert — quasars:
magnetic fields

##  1 Introduction

Quasi-periodic eruptions (QPEs) are new phenomena exhibiting quasi-periodic
rapid and high-amplitude bursts in soft X-ray, which were reported by Miniutti
et al. (2019) firstly in a low-mass Seyfert 2 galaxy GSN 069. The burst
duration and period of GSN 069 are about 1 and 9 hours, respectively. Except
for GSN 069, four other QPE sources, i.e., RX J1301.9+2747 (Sun et al., 2013;
Giustini et al., 2020), 2MASS 02314715-1020112, 2MASX J02344872-4419325
(Arcodia et al., 2021) and XMMSL1 J024916.6-041244 (Chakraborty et al., 2021),
have been discovered recently. Compared with others, more physical properties
of GSN 069 have been inferred, such as the central black hole mass ($M_{\rm
BH}\sim 4\times 10^{5}\rm M_{\odot}$ with a factor of few uncertainty), the
X-ray spectrum in various stages during outbursts and the nearly constant
eruption period with alternating long/short QPE time-separation as well as
alternating strong/weak QPEs. All of these are very important for us to study
the physical origin of QPEs.

Several models have been proposed to explain the physics of QPEs: radiation
pressure instabilities in the transition zone between inner advection-
dominated accretion flow (ADAF) and outer thin disk (Sniegowska et al. 2020;
Pan et al. 2021, hereinafter PLC21); Roche lobe overflow from one/two stars
orbiting the central black hole Zhao et al. 2021; Metzger et al. 2021; mass
overflow at pericentre from a white dwarf around a super-massive black hole (a
near miss tidal disruption, King 2020); warped disk tearing (Raj & Nixon,
2021); self-lensing of binary massive black hole (Ingram et al., 2021); and
star-disk collisions (Suková et al., 2021; Xian et al., 2021). These models
can partly explain the outbursts of QPEs, but mainly focused on their short
timescale and regular period. So far, there is no model can fit both the
outburst period and the X-ray spectrum of GSN 069 simultaneously. The observed
X-ray spectra at various burst stages in GSN 069 have been successfully
reproduced by using a constant disk blackbody plus a variable blackbody
component (Miniutti et al., 2019). Together with the quasi-periodic outbursts,
we notice that all these characteristics may be explained by the model of
Sniegowska et al. (2020) and PLC21.

The radiation pressure dominated inner region of a standard thin disk is both
thermally and viscously unstable (Shakura & Sunyaev, 1973, 1976). Adopting
some special values, this unstable region can be shrunken to a narrow zone
between inner ADAF and outer thin disk (Sniegowska et al., 2020). This model
has been successfully applied to the repeating changing-look (CL) active
galactic nuclei (AGN). Due to the important role of magnetic field playing on
the formation of jet/winds, PLC21 investigated how the large-scale magnetic
field affects the outbursts of unstable region and found that both the period
and amplitude of outbursts can be significantly changed. However, the mass
accretion rate of GSN 069 should be very high due to its high luminosity
($\sim 10^{42}\rm{ergs^{-1}}$) and a relatively small black hole mass ($\sim
4\times 10^{5}\rm M_{\odot}$, Miniutti et al. 2019). In this case, the inner
ADAF will disappear and the thin disk can extend to the innermost stable
circular orbit (ISCO, Narayan et al. 1998). The high mass accretion rate will
also result on a larger unstable inner region in a thin disk and a much longer
outburst period, which is inconsistent with the behaviour of QPEs. However, as
shown in PLC21, this point can be resolved if we consider the effects of
large-scale magnetic field, which can take away both the energy and angular
momentum from disk and decrease its temperature significantly (Cao & Spruit,
2013; Li & Begelman, 2014; Li & Cao, 2019). In this work, we construct a disk
instability model base on PLC21 to explain the quasi-periodic eruptions of GSN
069.

##  2 Model

###  2.1 steady outer disk

Figure 1: The schematic picture of our model. The outer blue region represents
a stable thin disk dominated by gas pressure and the inner orange region is
the unstable zone dominated by radiation pressure.

The inner ADAF will disappear when the mass accretion rate is high enough.
Therefore, the disk contains only two regions: the inner unstable and outer
stable regions (see Figure 1). On the presence of large-scale magnetic field,
the unstable region can be constrained to a narrow zone just outside the ISCO.
In order to improve the precision of our model, we adopt the general
relativistic correction factors defined by Novikov & Thorne (1973) to modify
our equations:

$\mathscr{A}=1+\frac{a_{*}^{2}}{4\hat{r}^{2}}+\frac{a_{*}^{2}}{4\hat{r}^{3}},$
(1)

$\mathscr{B}=1+a_{*}\sqrt{\frac{1}{8\hat{r}^{3}}},$ (2)

$\mathscr{C}=1-\frac{3}{2\hat{r}}+2a_{*}\sqrt{\frac{1}{8\hat{r}^{3}}},$ (3)

$\mathscr{D}=1-\frac{1}{\hat{r}}+\frac{a_{*}^{2}}{4\hat{r}^{2}},$ (4)

$\mathscr{E}=1+\frac{a_{*}^{2}}{\hat{r}^{2}}-\frac{a_{*}^{2}}{2\hat{r}^{3}}+\frac{3a_{*}^{4}}{16\hat{r}^{4}},$
(5)

$\mathscr{F}=1-2a_{*}\sqrt{\frac{1}{8\hat{r}^{3}}}+\frac{a_{*}^{2}}{4\hat{r}^{2}},$
(6)

$\Omega_{\rm k}=\sqrt{\frac{GM}{R^{3}}}\frac{1}{\mathscr{B}},$ (7)

$l_{\rm k}=\sqrt{GMR}\frac{\mathscr{F}}{\mathscr{C}^{1/2}},$ (8)

where $a_{*}=cJ/GM^{2}$ is the dimensionless spin parameter of black hole,
$\hat{r}(=R/R_{\rm s})$, $J$, $R_{\rm s}=2GM/c^{2}$ and $M$ are the
dimensionless radius, angular momentum of black hole, Schwarzschild radius and
black hole mass, respectively.

As PLC21, the general relativistic continuity equation for the outer thin disk
with winds driven by large-scale magnetic field can be written by:

$\frac{\mathrm{d}\dot{M}}{\mathrm{d}R}+4\pi R\dot{m}_{\rm w}=0,$ (9)

where $\dot{m}_{\rm w}$ is the mass-loss rate per unit surface area of disk
(see PLC21 for the details).

The angular momentum equation can be rewritten as

$-\frac{1}{2\pi}\frac{\mathrm{d}(\dot{M}l_{\rm
k})}{\mathrm{d}R}-\frac{\mathrm{d}}{\mathrm{d}R}(R^{2}\mathscr{B}\mathscr{C}^{-1/2}\mathscr{D}T_{r\phi})+T_{\rm
m}R=0,$ (10)

where $T_{r\phi}$ and $T_{\rm m}$ are the viscous and the magnetic torques
exerted on the accretion disk, respectively (see PLC21 for the details).

While the viscosity in accretion disk is believed to be related to the
turbulence induced by the magnetorotational instability (MRI, Balbus & Hawley
1998), the viscous torque has been suggested to be either proportional to the
total pressure (Shakura & Sunyaev, 1973), or proportional to the gas pressure
(Sakimoto & Coroniti, 1981), or both (Czerny et al., 2009). In this work, we
adopt a general form for the viscous torque:

$T_{r\phi}=-2\alpha P_{\rm tot}^{1-\mu}P_{\rm gas}^{\mu}H,$ (11)

where $P_{\rm tot}$ and $P_{\rm gas}$ are the total and gas pressure,
respectively (Horiuchi & Kato, 1990; Done & Davis, 2008). This choice is
directly motivated by the requirement that the outer radius of unstable zone
should be close enough to ISCO with the parameters presented in Section 3, in
order to get a narrow unstable zone. The widely used $\alpha$-viscosity is
$\nu=-\alpha C_{\rm s}H(\mathrm{dln}\Omega_{\rm k}/\mathrm{dln}R)^{-1}$ (and
$T_{r\phi}=\nu\Sigma R\mathrm{d}\Omega_{\rm k}/\mathrm{d}R$), where $C_{s}$ is
sound speed of local disk and $H=(P_{\rm tot}/\rho)^{1/2}\Omega_{\rm k}^{-1}$
is the half thickness of disk. With $C_{\rm s}=(P_{\rm tot}/\rho)^{1/2}$, we
have $\nu\propto\alpha P_{\rm tot}$ (Shakura & Sunyaev, 1973), corresponding
to $\mu=0$ in equation 11. Taam & Lin (1984) proposed that the turbulent
velocity may be related to sound speed $C_{s}=(P_{\rm gas}/\rho)^{1/2}$
instead of $C_{s}=(P_{\rm tot}/\rho)^{1/2}$, which leads to
$\nu\propto\alpha(P_{\rm gas}P_{\rm tot})^{1/2}$, i.e., $\mu=0.5$. In this
case, the thermal instability of the disk is mostly suppressed. The parameter
of $\mu$ in this work is employed to describe to what extent the turbulent
velocity depends on the gas pressure. The exact value of $\mu$ is currently
unknown, as is its physical interpretation. Stability criteria applied to
X-ray binary systems suggest $\mu\sim 0.56$ (see Done & Davis (2008) for
details). However, in this work we are not attempting to model the general
behavior of accreting black holes, but rather to seek a possible explanation
of a relatively rare phenomenon (QPEs) in terms of disk instabilities. As
such, as discussed in Section 3, we adopt a value of $\mu=0.27$ that allows
for a small instability region close to the ISCO without affecting the
stability of the outer disk.

The total pressure reads:

$P_{\rm tot}=(1+\frac{1}{\beta_{1}})(P_{\rm gas}+P_{\rm rad}),$ (12)

where $\beta_{1}=(P_{\rm gas}+P_{\rm rad})/(B^{2}/8\pi)$, $B$ and $P_{\rm
rad}$ are magnetic strength and radiation pressure, respectively. In the
standard thin disk model, the zero viscous torque condition is usually adopted
at ISCO. However, an additional magnetic torque would be exerted on the ISCO
when large scale magnetic fields are present (e.g., Krolik 1999). The non-zero
torque on the inner boundary will increase the radiative efficiency by
producing additional dissipation (e.g., Agol & Krolik 2000). Following Agol &
Krolik (2000), we define a parameter $f$ in our model to describe the effect
of non-zero torque condition at ISCO in Kerr metric:

$T_{r\phi}\bigg{|}_{R=R_{\rm ISCO}}=-(1-f)\left[\frac{\dot{M}l_{\rm
k}\mathscr{C}^{1/2}}{2\pi\mathscr{BD}R^{2}}\right]_{R=R_{\rm ISCO}},$ (13)

where the parameter $f$ can be related to the parameter $f_{\rm ms}$ in Agol &
Krolik (2000) as:

$f_{\rm ms}=\frac{3(1-f)}{2x_{\rm ms}^{2}\mathscr{C}^{1/2}_{\rm ms}}.$ (14)

Our $f$ parameter can take any value between $-\infty$ and 1, while $f_{ms}$
is limited to vary between 0 and $\infty$. The zero-torque condition
corresponds to $f=1$ in our model instead of $f_{ms}=0$.

Without an inner ADAF region, the energy equation will return to the classical
form as:

$-\frac{3}{2}\Omega_{\rm
k}T_{r\phi}\frac{\mathscr{BD}}{\mathscr{C}}=\frac{8acT_{\rm c}^{4}}{3\tau}.$
(15)

The optical depth $\tau$ is given by $\tau=\bar{\kappa}\Sigma/2$, where
$\bar{\kappa}$ and $\Sigma=2\rho H$ are the opacity and surface density,
respectively.

###  2.2 variable inner unstable region

The instability of a standard thin accretion disk in Kerr metric had been
investigated by Xue et al. (2011), where the continuity equation is given as :

$u^{t}\frac{\partial\Sigma}{\partial t}=-\frac{1}{r}\frac{\partial}{\partial
r}(\Sigma ru^{r})-\Sigma\frac{\partial u^{t}}{\partial t},$ (16)

where $u^{t}$ and $u^{r}$ are the time and radial components of four-velocity,
respectively. Similar with PLC21, the equation (16) can be rewritten as:

$\begin{split}&\left[u^{t}-\frac{C_{\rm
H}H\left(1-\beta_{2}\right)}{\Sigma\left(1+\beta_{2}\right)}\right]\frac{\mathrm{d}\Sigma}{\mathrm{d}t}+\frac{C_{\rm
H}H\left(4-3\beta_{2}\right)}{T\left(1+\beta_{2}\right)}\frac{\mathrm{d}T}{\mathrm{d}t}\\\
&-\frac{\dot{M}_{0}-\dot{M}-4\pi R\dot{m}_{\rm w}\Delta R}{2\pi R\Delta
R}=0,\end{split}$ (17)

where $\dot{M}_{0}$ and $\Delta R$ are the inflow rate and the width of the
unstable zone, respectively. $\beta_{2}=P_{\rm gas}/(P_{\rm gas}+P_{\rm rad})$
is the ratio of gas pressure to the sum of gas pressure and radiation
pressure. In order to simplify the form of equation, we define a coefficient
$C_{\rm H}$ as:

$C_{\rm H}=-\frac{u^{t}u^{r}\Sigma\dot{M}_{0}}{2\pi
RH_{0}\Sigma_{0}\left[c^{2}\mathscr{D}+(u^{r})^{2}\right]},$ (18)

where $H_{\rm 0}$ and $\Sigma_{\rm 0}$ are the scale-height and surface
density on the outer boundary of unstable zone, respectively.

The energy conservation equation is (Xue et al., 2011):

$\Sigma T\left[u^{t}\frac{\partial S}{\partial t}+u^{r}\frac{\partial
S}{\partial r}\right]=Q^{+}-Q^{-},$ (19)

where $S$, $Q^{+}$ and $Q^{-}$ are the entropy of gas, the viscous heating
rate and the cooling rate, respectively. The second term in the left side of
Equation (19) represents the advection rate of energy $Q_{\rm adv}$, which can
be written as (see PLC21):

$Q_{\rm adv}=\frac{\dot{M}PH}{2\pi R\Delta R\Sigma\mathscr{D}^{\frac{1}{2}}}.$
(20)

Therefore, the evolution equation of temperature can be rewritten as:

$\displaystyle\frac{\mathrm{d}T}{\mathrm{d}t}=$
$\displaystyle\frac{T(Q^{+}-Q^{-}-Q_{\rm
adv})(1+\frac{1}{\beta_{1}})(1+\beta_{2})}{2PHu^{t}(28-22.5\beta_{2}-1.5\beta_{2}^{2}+\frac{12-9\beta_{2}}{\beta_{1}})}$
(21)
$\displaystyle+2\frac{T\mathrm{d}\Sigma}{\Sigma\mathrm{d}t}\frac{4-3\beta_{2}+\frac{2-\beta_{2}}{\beta_{1}}}{28-22.5\beta_{2}-1.5\beta_{2}^{2}+\frac{12-9\beta_{2}}{\beta_{1}}}.$

##  3 Numerical Results

With equations (17) and (21), we can study the limit-cycle behaviour of thin
disk. The overall model is complex and depends on several free parameters.
Exploring the full parameter space to compare the model with the available
X-ray data would require to generate a very large number of simulations that
would take decades to run. Moreover, the large number of free parameters that
is inherent to any disk instability model implies that general model
predictions can be obtained, but that different details can be produced
easily. This suggests to look for qualitative comparisons with the data rather
than detailed fits, at least in this initial stage. We adopt here a fiducial
model with a black hole mass of $2\times 10^{5}\rm M_{\odot}$ and a high black
hole spin value of 0.98. A discussion on the choice of parameters for the
fiducial model is given in Section 4. This value of mass is 2 times smaller
than that given by Miniutti et al. (2019), but is still range within the
estimated black hole mass. Other parameters are conveniently chosen as
$f=0.9$, $\dot{m}=0.08\dot{M}_{\rm Edd}$ ($\dot{M}_{\rm Edd}=1.5\times
10^{18}M/M_{\odot}\rm{gs^{-1}}$), $\alpha=0.15$, $\beta_{1}=38$, $\mu=0.27$
and $\Delta R=0.12R_{\rm s}$ (see Group 1 in Table 1). Here the outer radius
of $\Delta R$ is given self-consistently by the disk instability criteria (see
section 4 for details).

GSN 069 is observed to display QPEs about every 9 hours. During the outbursts,
the X-ray flux can increase by about two order of magnitudes (Miniutti et al.,
2019). We compare the results of our model with the observed quasi-periodic
0.4-2 keV light curve of GSN 069 in Figure 2. It is found that our theoretical
results can qualitatively fit both the period and the outburst duration,
although the alternating strong/weak peaks and short/long recurrence times are
not well reproduced. Furthermore, in Figure 3, we present the normalized light
curves for the three different energy bands given by Figure 2b of Miniutti et
al. (2019). The narrower width found in high energy bands is also
qualitatively consistent with observational results of GSN 069. On the other
hand, QPEs in higher energy bands peak earlier in GSN 069, and our model does
not reproduce this behavior.

Figure 2: The 0.4-2 keV light curve, where the black line represents the light
curve given by our model and the colored lines are the seven phase segments of
eruption to fit the X-ray spectrum in Figure 4. The grey dots are the data
observed by XMM-Newton Director Discretionary Time (DDT) performed on 2019
January 16/17. Figure 3: The light curves normalized by the peak flux in three
different energy bands. The red, green and blue lines correspond to the
0.2-0.3 $\rm keV$, 0.4-0.5 $\rm keV$ and 0.8-1 $\rm keV$ light curves,
respectively. Figure 4: Phase-resolved Spectral analysis of the GSN 069. The
colored points with various shapes are the observational data after doing
corrections of instrumental effective area and absorption.  Figure 5: The
0.4-2 keV light curves and phase-resolved Spectral analysis under different
parameter groups. Figure 5(1), 5(2), and 5(3) correspond to $M=4\times
10^{5}M_{\odot}$, $f=1$(zero torque condition), and $a_{*}=0.9$, respectively.
The best-fitting parameters of Figure 5(1), 5(2), and 5(3) correspond to Group
2, Group 3, and Group 4 in Table 1, respectively.

The timing properties of GSN 069 has been reproduced by several works (e.g.,
King 2020; Raj & Nixon 2021; Xian et al. 2021). However, besides the timing
properties, the X-ray spectrum during outburst have also been well resolved by
XMM-Newton (see Figure 4 in Miniutti et al. 2019, where a whole outburst is
divided into seven phases). We adopt seven phase segments of outburst
($\rm{Q_{pre}}$, $\rm{R_{1}}$, $\rm{R_{2}}$, $\rm{P}$, $\rm{D_{1}}$,
$\rm{D_{2}}$, $\rm{Q_{post}}$, see the colored lines in Figure 2) to fit the
observed X-ray spectrum. It is found that all the X-ray spectrum during
outburst can be qualitatively reproduced by our model. The average effective
temperature of unstable zone corresponding to the seven segments are $kT_{\rm
Q_{pre}}=57\rm eV$, $kT_{\rm R_{1}}=67\rm eV$, $kT_{\rm R_{2}}=81\rm eV$,
$kT_{\rm P}=91\rm eV$, $kT_{\rm D_{1}}=82\rm eV$, $kT_{\rm D_{2}}=68\rm eV$
and $kT_{\rm Q_{post}}=59\rm eV$, respectively. The phase-resolved X-ray data
in Figure 4 are taken from Figure 3 in Miniutti et al. (2019) after correcting
the instrumental effective area and galaxy absorption, in which the cross-
sections given by Morrison & McCammon (1983) and the equivalent hydrogen
column $n_{\rm H}$ ($=3\times 10^{20}\rm cm^{-2})$ are adopted. Some of the
spectra show a high-energy excess with respect to our model prediction, i.e.
the agreement is qualitative but there is still room for improvement, as
discussed in the following section.

##  4 Choice of parameters in our model

In addition to the width of unstable region $\Delta R$, our model include 7
free parameters, i.e., the black hole mass $M$, black hole spin $a_{*}$, $f$,
mass accretion rate $\dot{m}$, $\alpha$, $\beta_{1}$ and $\mu$. $\Delta R$ is
self-consistently given by the disk instability criteria given in Kato et al.
(1998), which indicates that the disk is unstable when
$4-10\beta_{2}-7\mu+7\beta_{2}\mu>0$ , where $\beta_{2}$ is determined by the
mass accretion rate and radius. There are five main characteristics required
to fit the data points in Figures 2 and 4 well, i.e., 1) the eruption period;
2) the eruption duration and shape; 3) the amplitude of eruption; 4) the
constant temperature component during spectral evolution (representing the
constant outer disk emission); 5) a variable temperature component during
spectral evolution. As mentioned, our goal here is to provide a possible
framework for QPE production, obtaining results that are close enough, both in
timing and spectral properties, to the observed data to suggest that this type
of modeling is worth exploring further in the future.

The maximal disk temperature inferred from observation is $kT_{\rm max}\sim
50\rm eV$. However, for a thin accretion disk surrounding a black hole with
mass of $4\times 10^{5}M_{\odot}$, the maximal disk temperature is only
$kT_{\rm max}\sim 30\rm eV$
($kT_{\rm{max}}=11.5\left(M/10^{8}M_{\odot}\right)^{-1/4}\dot{m}^{1/4}(\rm
eV)$). In order to match the disk temperature before outburst (the constant
temperature component), a small black hole mass $2\times 10^{5}M_{\odot}$, a
high black hole spin ($a=0.98$) and a non-zero viscous torque are all
simultaneously necessary at the first step. We notice that a spectral
hardening factor $\kappa$ ($\sim 1.7$) between the effective temperature of
disk $T_{\rm eff}$ and observed disk temperature $T_{\rm max}$ ($T_{\rm
max}\sim\kappa T_{\rm eff}$) is adopted in Miniutti et al. (2019) when fitting
the black hole mass. With the parameter $\kappa$, a disk temperature $T_{\rm
eff}\sim 30\rm eV$ is enough to fit the observed spectra. In this case, a
black hole with mass of $4\times 10^{5}M_{\odot}$ should be able to fit the
observational data well. Our model ignores any spectral hardening, and
therefore the matching with the observed constant temperature can only be
obtained with a smaller black hole mass. In Figure 5(1), we investigate the
effects of black hole mass $M$ on the 0.4-2 kev light curves and phased
resolved spectral analysis. During the process of fitting, the eruption period
and duration are the first two characteristics should be matched. It is found
that the parameters with $M=4\times 10^{5}M_{\odot}$ (the detailed parameters
adopted are listed in Table 1.) can reproduce a similar evolution trend for
both the light curves and the spectral analysis, but provide a much worse
match to the data. The worse fit of the spectral evolution (see Figure 5(1)b
and (1)c) is driven by a too low constant temperature of the thermal component
that is present in all fits, and could likely be cured by introduced a
hardening factor, as explained above. On the other hand, QPEs have the right
recurrence time but are too long (broad), as shown in Figure 5(1)a. However,
if the mismatch of the spectral evolution fits can be ignored due to the
hardening factor, the QPEs may be made narrower (shorter) by playing with a
small value of $\mu$ (see Figure 6). Next, we adjust the other four parameters
($\alpha$, $\beta_{1}$, $\mu$ and $\dot{m}$) to reproduce the other four
characteristics observed. The parameters $\mu$ can change the ratio of
eruption duration to the period (see Figure 6). How the last three parameters
affect the outburst behaviour had been reported by PLC21. It is found that the
stronger magnetic field strength (smaller $\beta_{1}$) will result on both the
decreasing eruption period and higher eruption amplitude. However, increasing
$\alpha$ can only increase the eruption period (see PLC21 for details). The
increase of mass accretion rate will improve the maximal disk temperature
($T_{\rm max}\sim\dot{m}^{1/4}$) and the width of unstable zone, resulting on
the enhancement of all the five characteristics.

The non-zero torque in equation 13 will emerge when considering the magnetic
torque exerted on the gas at ISCO by the in-falling gas inside the ISCO
(Krolik, 1999). However, this torque should be small due to the small scale-
height and the weak magnetic field in a thin disk (Nixon & Pringle, 2021),
which is also suggested by GRMHD simulations (Noble et al., 2010). In our
model, a small non-zero torque $f=0.9$ ($\sim 10\%$ of the angular momentum at
ISCO) is adopted, which is consistent with the results of simulations ($f=1$
corresponding to the zero torque condition). We give the results with zero
torque in Figure 5(2). It is found that the results with non-zero torque are
much better than that with zero torque because the effective gas temperature
close to ISCO can be greatly increased with non-zero torque (the effective
temperature at ISCO is zero for zero torque).

The spin of black hole is adopted as 0.98 in this work. As discussed by Thorne
(1974), the maximal black hole spin can reach 0.998. However, the maximal
value of spin will be smaller than 0.998 when taking the non-zero torque into
account, which seems to be inconsistent with our assumption. However, the
torque at ISCO is small as suggested by the MHD simulations and thus it will
not affect the maximal value of spin significantly. Furthermore, QPEs are
likely transient phenomena in the accretion history of a given source (no QPEs
are observed in GSN 069 after June 2020), and if they correspond to non-zero
torque transient phases, the black hole spin could well have been established
during past long-lived zero-torque phases. Our model only allows a slightly
lower BH spin in order to fit the observational data well (see Figure 5(3),
the fitting results with $a_{*}=0.9$ are a little worse but still acceptable).
Finally, although the width of the unstable region is small ($\Delta
R=0.12R_{\rm s}$), it can still significantly change the 0.4-2 kev light curve
(Figure 2). However, the variation of bolometric luminosity is less than
$20\%$ for the reason that the light curves are almost constant at optical-UV
bands, which is consistent with observational results (Miniutti et al., 2019).

Table 1: Detailed parameter of our fitting Number | $M$ | $a_{*}$ | $f$ | $\dot{m}$ | $\alpha$ | $\beta_{1}$ | $\mu$
---|---|---|---|---|---|---|---
1 | $2\times 10^{5}$ | 0.98 | 0.9 | 0.08 | 0.15 | 38 | 0.27
2 | $4\times 10^{5}$ | 0.98 | 0.9 | 0.08 | 0.15 | 22 | 0.3
3 | $2\times 10^{5}$ | 0.98 | 1 | 0.08 | 0.15 | 50 | 0.24
4 | $2\times 10^{5}$ | 0.9 | 0.9 | 0.08 | 0.15 | 35 | 0.25

##  5 Conclusions and Discussion

Figure 6: The 0.4-2 keV light curve, where $L_{\rm const}$ is the luminosity
of outer constant disk and the value of $\mu$ in the top, middle and bottom
panels are $0.2$, $0.3$ and $0.4$, respectively. Other parameters adopted in
this figure are: $f=0.9$, $\dot{m}=0.1\dot{M}_{\rm Edd}$, $\alpha=0.15$,
$\beta_{1}=50$ and $\Delta R=0.1R_{\rm s}$.

The mechanism of QPEs is still unclear. In this work, we adopt a disk
instability model to explore the physical origin of QPEs in GSN 069. The model
of PLC21 is improved in that: 1) importing the non-zero torque on the inner
boundary of thin disk; 2) adopting a general form of viscous stress tensor
($T_{r\phi}=-2\alpha P_{\rm tot}^{1-\mu}P_{\rm gas}^{\mu}H$); 3) adopting a
Kerr metric. Our model can qualitatively fit the 0.4-2keV light curves, the
light curves at different energy bands and the phase-resolved X-ray spectrum
(Figures 2, 3, 4). However, we emphasize that the parameters adopted above may
not represent the real best-fitting set. Exploring the full parameter space
is, however, beyond the scope of this work as a full set of simulations would
require several decades of computing time (by personal computer). A more
quantitative comparison with data could be possible by fixing some/several
free parameters to fiducial values, and by considering only a limited set of
variable ones. We defer this study to future work.

The period, outburst duration and light-curve profile of QPEs are source-
dependent (see e.g., Miniutti et al. 2019; Giustini et al. 2020; Arcodia et
al. 2021). The period and amplitude of outburst can be changed significantly
by the presence of large-scale magnetic field (PLC21). Furthermore, we suggest
that the $\mu$ parameter in the general form of viscous torque can be an
useful tool to fit the light-curve profile of QPEs (Figure 6). It is found
that the ratio of outburst duration to the period increases with increasing
$\mu$, which means that our model may also be applicable to other QPEs.
Although our model can qualitatively fit the characteristics of light curves
in GSN 069, some observational details are still problematic for our toy
model, e.g., the asymmetric eruptions (a faster rise and slower decay is seen
in some of the QPE sources); the peak delay at different energy band (the flux
of higher energy bands peak at earlier times); the hard X-ray excess ($\sim
1\rm keV$) that is not well reproduced by our model during the $\rm{R_{1}}$,
$\rm{R_{2}}$ and $\rm{P}$ phases; and the alternating long/short and
strong/weak QPEs. The first two issues may be resolved if we allow the
unstable zone to propagate outwards instead of remaining constant (Li et al.,
2007), the radiative area of unstable zone will increase after eruption and
result on the the asymmetric eruptions and the peak delay at different energy
band. The hard X-ray excess at $\rm{R_{1}}$, $\rm{R_{2}}$ and $\rm{P}$ phases
may be related to the appearance of a warm comptonizing corona, as generally
observed in AGN (Petrucci et al., 2020). The alternating strong/weak QPEs and
long/short recurrence times may be partly solved if the unstable zone is
variable instead of a constant due to the slight variation of magnetic field,
perhpas responding to the QPEs themselves in a sort of feedback loop. We can
expect a longer recurring time when the eruptions are stronger, which will
drag more gas from accretion disk into the black hole. We will further
investigate these issues in the future work.

## ACKNOWLEDGEMENTS

We thank the referee for his/her helpful suggestions and comments. This work
is supported by the NSFC (grants 11773056, 11773050, 11833007, 11873073, and
12073023) and partly funded by Project No. MDM-2017-0737 Unidad de Excelencia
"María de Maeztu" - Centro de Astrobiología (CSIC-INTA) by
MCIN/AEI/10.13039/501100011033. Support for this work was also provided by the
science research grants from the China Manned Space Project with NO. CMS-
CSST-2021-A06 and Shanghai Pilot Program for Basic Research - Chinese Academy
of Science, Shanghai Branch (JCYJ-SHFY-2021-013).

## References

  * Agol & Krolik (2000) Agol, E., & Krolik, J. H. 2000, ApJ, 528, 161, doi: 10.1086/308177
  * Arcodia et al. (2021) Arcodia, R., Merloni, A., Nandra, K., et al. 2021, Nature, 592, 704–707, doi: 10.1038/s41586-021-03394-6
  * Balbus & Hawley (1998) Balbus, S. A., & Hawley, J. F. 1998, Reviews of Modern Physics, 70, 1, doi: 10.1103/RevModPhys.70.1
  * Cao & Spruit (2013) Cao, X., & Spruit, H. C. 2013, ApJ, 765, 149, doi: 10.1088/0004-637X/765/2/149
  * Chakraborty et al. (2021) Chakraborty, J., Kara, E., Masterson, M., et al. 2021, ApJ, 921, L40, doi: 10.3847/2041-8213/ac313b
  * Czerny et al. (2009) Czerny, B., Siemiginowska, A., Janiuk, A., Nikiel-Wroczyński, B., & Stawarz, Ł. 2009, ApJ, 698, 840, doi: 10.1088/0004-637X/698/1/840
  * Done & Davis (2008) Done, C., & Davis, S. W. 2008, ApJ, 683, 389, doi: 10.1086/589647
  * Giustini et al. (2020) Giustini, M., Miniutti, G., & Saxton, R. D. 2020, A&A, 636, L2, doi: 10.1051/0004-6361/202037610
  * Horiuchi & Kato (1990) Horiuchi, T., & Kato, S. 1990, PASJ, 42, 661
  * Ingram et al. (2021) Ingram, A., Motta, S. E., Aigrain, S., & Karastergiou, A. 2021, MNRAS, 503, 1703, doi: 10.1093/mnras/stab609
  * Kato et al. (1998) Kato, S., Fukue, J., & Mineshige, S. 1998, in Black-hole Accretion Disks Edited by Shoji Kato, Jun Fukue, & Sin Mineshige Publisher: Kyoto, Japan: Kyoto University Press, Isbn
  * King (2020) King, A. 2020, MNRAS, 493, L120, doi: 10.1093/mnrasl/slaa020
  * Krolik (1999) Krolik, J. H. 1999, ApJ, 515, L73, doi: 10.1086/311979
  * Li & Cao (2019) Li, J., & Cao, X. 2019, ApJ, 872, 149, doi: 10.3847/1538-4357/ab0207
  * Li & Begelman (2014) Li, S.-L., & Begelman, M. C. 2014, ApJ, 786, 6, doi: 10.1088/0004-637X/786/1/6
  * Li et al. (2007) Li, S.-L., Xue, L., & Lu, J.-F. 2007, ApJ, 666, 368, doi: 10.1086/519834
  * Metzger et al. (2021) Metzger, B. D., Stone, N. C., & Gilbaum, S. 2021, arXiv e-prints, arXiv:2107.13015.  https://arxiv.org/abs/2107.13015
  * Miniutti et al. (2019) Miniutti, G., Saxton, R. D., Giustini, M., et al. 2019, Nature, 573, 381, doi: 10.1038/s41586-019-1556-x
  * Morrison & McCammon (1983) Morrison, R., & McCammon, D. 1983, ApJ, 270, 119, doi: 10.1086/161102
  * Narayan et al. (1998) Narayan, R., Mahadevan, R., & Quataert, E. 1998, in Theory of Black Hole Accretion Disks, ed. M. A. Abramowicz, G. Björnsson, & J. E. Pringle, 148–182.  https://arxiv.org/abs/astro-ph/9803141
  * Nixon & Pringle (2021) Nixon, C. J., & Pringle, J. E. 2021, New A, 85, 101493, doi: 10.1016/j.newast.2020.101493
  * Noble et al. (2010) Noble, S. C., Krolik, J. H., & Hawley, J. F. 2010, ApJ, 711, 959, doi: 10.1088/0004-637X/711/2/959
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), 343–450
  * Pan et al. (2021) Pan, X., Li, S.-L., & Cao, X. 2021, ApJ, 910, 97, doi: 10.3847/1538-4357/abe766
  * Petrucci et al. (2020) Petrucci, P. O., Gronkiewicz, D., Rozanska, A., et al. 2020, A&A, 634, A85, doi: 10.1051/0004-6361/201937011
  * Raj & Nixon (2021) Raj, A., & Nixon, C. J. 2021, ApJ, 909, 82, doi: 10.3847/1538-4357/abdc25
  * Sakimoto & Coroniti (1981) Sakimoto, P. J., & Coroniti, F. V. 1981, ApJ, 247, 19, doi: 10.1086/159005
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astronomy and Astrophysics, 24, 337.  http://adsabs.harvard.edu/abs/1973A%26A....24..337S
  * Shakura & Sunyaev (1976) Shakura, N. I., & Sunyaev, R. A. 1976, MNRAS, 175, 613, doi: 10.1093/mnras/175.3.613
  * Sniegowska et al. (2020) Sniegowska, M., Czerny, B., Bon, E., & Bon, N. 2020, A&A, 641, A167, doi: 10.1051/0004-6361/202038575
  * Suková et al. (2021) Suková, P., Zajaček, M., Witzany, V., & Karas, V. 2021, ApJ, 917, 43, doi: 10.3847/1538-4357/ac05c6
  * Sun et al. (2013) Sun, L., Shu, X., & Wang, T. 2013, ApJ, 768, 167, doi: 10.1088/0004-637X/768/2/167
  * Taam & Lin (1984) Taam, R. E., & Lin, D. N. C. 1984, ApJ, 287, 761, doi: 10.1086/162734
  * Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507, doi: 10.1086/152991
  * Xian et al. (2021) Xian, J., Zhang, F., Dou, L., He, J., & Shu, X. 2021, ApJ, 921, L32, doi: 10.3847/2041-8213/ac31aa
  * Xue et al. (2011) Xue, L., Sądowski, A., Abramowicz, M. A., & Lu, J.-F. 2011, ApJS, 195, 7, doi: 10.1088/0067-0049/195/1/7
  * Zhao et al. (2021) Zhao, Z. Y., Wang, Y. Y., Zou, Y. C., Wang, F. Y., & Dai, Z. G. 2021, arXiv e-prints, arXiv:2109.03471.  https://arxiv.org/abs/2109.03471