#  Critical Drift in a Neuro-Inspired Adaptive Network
Supplemental material

Silja Sormunen  Department of Computer Science, Aalto University, 00076 Espoo,
Finland Thilo Gross  Helmholtz Institute for Functional Marine Biodiversity at
the University of Oldenburg (HIFMB),26129 Oldenburg, Germany Alfred-Wegener
Institute, Helmholtz Centre for Marine and Polar Research, 27570 Bremerhaven,
Germany Institute for Chemistry and Biology of the Marine Environment (ICBM),
Carl-von-Ossietzky University, 26129 Oldenburg, Germany Jari Saramäki
Department of Computer Science, Aalto University, 00076 Espoo, Finland

##  I Critical values of $\langle k\rangle$ and $\lambda_{1}$

Fig 1(a) shows the empirical estimates for the critical mean degree $\langle
k\rangle^{*}_{N,\ \rm static}$ in static ER networks as a function of $1/N$.
The critical values are determined by first finding the mean degree for which
the size/duration of finite avalanches in 1000 successive runs is maximal,
then repeating this procedure for 30 network realizations and taking the
average of the found values. We fit a quadratic polynomial to the points to
find the intersection with the $y$-axis. This intersection lies around
$\langle k\rangle^{*}_{N\to\infty,\ \rm static}=2.27$, slightly above the
theoretical estimate $\langle\hat{k}\rangle^{*}_{\rm static}=2.21$.

In the evolved ER networks of size $N=10^{5}$, the leading eigenvalue
$\lambda_{1}$ self-organizes to a value close to $\langle
k\rangle^{*}_{N=10^{5},\ \rm static}$ during the drift phase, and this value
moves closer to the theoretical estimate
$\hat{\lambda}_{1}^{*}=\langle\hat{k}\rangle^{*}_{\rm static}$ as the network
size increases [Fig. 1(b)].

Figure 1: (a) Critical value $\langle k\rangle^{*}_{N,\ \rm static}$ in static
ER networks as a function of $1/N$ calculated using avalanche sizes (circles)
and durations (triangles). The markers overlay the 95% confidence intervals.
The dashed line marks the theoretical estimate $\langle\hat{k}\rangle^{*}_{\rm
static}=2.21$. (b) The mean value of $\lambda_{1,N}$ in evolved networks
during the drift phase as a function of $1/N$. The initial degree is 3 and
$\mu$ is set to $1/(100N)$. The extrapolated value
$\lambda^{*}_{1,N\to\infty}$ lies close to the value $\langle
k\rangle^{*}_{N\to\infty,\ \rm static}$.

##  II Effect of network size on the drift

We verify that the increase in the mean degree $\langle k\rangle$ during the
drift phase does not arise from some finite size effect. To this end, we first
calculate the average value of $\lambda_{1}$ after it has clearly stabilized,
and subsequently define the drift to start when $\lambda_{1}$ first crosses
this value. Next, we subtract the value of $\langle k\rangle$ at the start of
the drift from the subsequent values of $\langle k\rangle$ during the drift.
After this normalization, it becomes evident that increasing $N$ has very
little effect on the drift for networks larger than $N=50000$ [Fig. 2].

Figure 2: Increase in $\langle k\rangle$ during the drift phase. Initial mean
degree is 1 and the spontaneous firing rate equals $\frac{1}{100N}$.

##  III Parameter considerations

Figure 3: Time evolution of $\langle k\rangle$ and $\lambda_{1}$ normalized by
the stabilized value of $\lambda_{1}$ in a network of size $N=10^{5}$ for
different parameter values of $\eta$, $l$ and $\epsilon$. The default
parameter values are $\eta=10^{-7},l=0.001$ and $\epsilon=0.001$. As the
parameters $l$ and $\epsilon$ control the timescale of the drift, we have
additionally scaled the time series in figures (b) and (c) to correspond to
the pace of the drift for the default values. The scaling is done by
multiplying the $x$-axis by the factor by which the parameter value is larger
than 0.001.

###  A Spontaneous firing rate $\eta$

As discussed in [2], the spontaneous firing rate $\eta$ plays a role in how
close to the critical threshold $\lambda_{1,N}^{*}$ the network self-organizes
to; when $\eta$ is increased, the stabilized value of $\lambda_{1,N}$
decreases. This happens because spontaneous activity contributes to the
overall level of activity, which the self-organization process strives to keep
at a constant level. However, as demonstrated in [2], the effect of the exact
value of $\eta$ on the distance to criticality becomes negligible for small
enough $\eta$ and large enough networks. We observe that for a network of size
$N=10^{5}$, the value of $\eta$ has little effect also on the drift of
$\langle k\rangle$ when $\eta<10^{-4}$ [Fig. 3(a)].

In networks of finite size, the value of $\eta$ affects how smoothly the mean
degree changes. If $\eta$ is low, intervals between subsequent bursts of
activity lengthen, and the network is likely to cross well into the
supercritical regime before a new avalanche takes place and the link removal
mechanism is activated. Consequently, the fluctuations between sub- and
supercritical states become more pronounced. However, this effect subsides as
the network size increases.

###  B Network evolution parameters

In general, self-organized criticality is possible only if the system’s
dynamics can be separated into two parts; the dynamics that becomes critical,
and the controlling dynamics that steers the former to criticality.
Dynamically speaking, this division is justified only if a timescale
separation exists between the two. Consequently, the firing dynamics need to
have a faster timescale than the network evolution, _i.e._
$\beta,\delta,\gamma\gg l,g$. In addition, the rate $l$ at which nodes lose
links needs to be clearly higher than the rate $g$ at which they gain links.
As we are interested in the magnitude of $g$ compared to $l$, we define
$g=\epsilon l$, where $\epsilon$ is a scaling parameter. This second time-
scale separation stems from the fact that the plasticity rules act locally
instead of controlling the average activity in a more centralized manner.
Since firing nodes characterize only the active phase but inactive nodes are
common in both the quiescent and the active phase, local excitability has to
be increased gradually for inactive nodes and decreased quickly for active
nodes. As shown in [2], the network self-organizes to criticality in the limit
$l\to 0,\epsilon\to 0$.

The exact values of $l$ and $\epsilon$ have little effect on the magnitude of
increase in $\langle k\rangle$ during the drift phase [see Fig. 3(b),(c)].
Both parameters control the timescale of the drift; if one of the parameters
is decreased by an order of magnitude, the duration of the drift
correspondingly increases by an order of magnitude.

##  IV Scaling of signatures of criticality

In Figure 2(e) in the main text, we observed that the divergence peaks of
avalanche sizes and durations did not always align with the point where the
average activity $\langle[F]\rangle$ abruptly increased. Here, we show that
this discrepancy is a finite-size effect, and that the different ways of
measuring criticality converge as the network size increases.

To show this, we repeat the procedure of Fig. 2(e) in the main text for
networks of different sizes. We then determine the critical degree $\langle
k\rangle^{*}_{N}$ in three different ways; either based on where the curve for
maximum size or duration of finite avalanches reaches its maximal value, or by
looking at where the average activity $\langle[F]\rangle$ changes the most
between two consecutive mean degrees. We observe that these different ways to
determine $\langle k\rangle^{*}_{N}$ converge as the network size $N$
increases [Fig. 4].

Looking at Fig. 4, we observe that as $N$ increases, the mean degree of the
evolved networks seems to converge to a value slightly above $\langle
k\rangle^{*}_{N}$. In fact, this is to be expected; as discussed in the
previous section, the network self-organizes to criticality in the limit
$l,\epsilon\to 0$. Hence, the evolved mean degree is expected to lie slightly
above the critical value for finite values of $l$ and $\epsilon$.

[2] derive a theoretical estimate for this deviation for specific values of
$l$ and $\epsilon$ under the assumption that the network is an ER network (see
Eq. (11g) of their paper). In Fig. 4, the dashed line marks this theoretical
deviation from the theoretical critical value $\langle k\rangle^{*}=2.21$. We
observe that the mean degree of the simulated networks would seem to deviate
from the critical value less than suggested by Eq. (11g). This may be affected
by the evolved networks violating the assumption of the network being an ER
network. In addition, since the theoretical deviation depends on $\langle
k\rangle^{*}_{N}$ for $l,\epsilon\neq 0$, the true distance to criticality may
vary slightly as $\langle k\rangle^{*}_{N}$ increases during the drift.
Finally, the mean degree is expected to converge to the theoretical value from
below only as the spontaneous activity goes to zero.

Figure 4: Distance to criticality at different times during the network
evolution as a function of inverse network size $1/N$. Distance to criticality
is measured based on the average activity (red circles), maximum size (blue
triangles) and duration (green diamonds) of finite avalanches. The initial
mean degree is 2 and the spontaneous firing rate is set to $1/(100N)$. For
each $N$ and each timepoint, we first determine the estimates of $\langle
k\rangle^{*}_{N}$ for one time series by averaging over 10 rounds of degree
manipulations. The figures display the average and the 95% confidence
intervals over 30 time series. Note that the resolution of the $y$-axis is
0.01 before any averaging.

##  V PDFs of avalanche sizes and durations

Figure 5 shows the PDFs of the sizes and durations of finite avalanches in
static networks with different mean degrees. When the mean degree equals the
experimentally extrapolated value $\langle k\rangle^{*}_{\rm
static,N\to\infty}=2.27$, the distributions resemble power law distributions
with exponents matching the theoretical predictions for critical SIS-like
systems derived in [4]. When the mean degree lies clearly below or above the
critical value, the distributions’ tails decay markedly faster.

We observe that the tail of avalanche sizes is in general heavier than that of
avalanche durations, which aligns with the results of [1] and [4]. The
beginning of the distribution of avalanche durations (right-hand column in
Fig. 5) is furthermore affected by the exponential distribution of the times
that nodes spend in the firing state.

To further verify that our self-organizing system resides at criticality
during the drift phase, we plot the PDFs of avalanche sizes and durations at
different points during the drift [Fig. 6]d. We observe that these
distributions resemble power law distributions with exponents matching the
theoretical predictions at criticality. Most of the distributions have a small
bump in their tail, which might indicate that the systems are slightly
supercritical. This result would align with the theoretical prediction derived
in [2] already discussed in sections III and IV, namely that our model is
strictly critical in the limit $l,\epsilon\to\infty$ and deviates slightly to
the supercritical regime for finite parameter values. When the network size
increases, this bump moves further towards the tail [see Fig. 7].

Figure 5: PDFs of the sizes and durations of finite avalanches in static
random networks with different mean degrees and varying network sizes. The
mean degrees of the top and bottom rows are chosen to be below and above
criticality, respectively, while the middle row corresponds to the
experimentally extrapolated critical value $\langle k\rangle^{*}_{\rm
static,N\to\infty}=2.27$. The number of runs for each mean degree is
originally 10000, but only avalanches lasting less than $t_{\rm max}=50000$
are considered. The green dashed lines mark the slopes derived for critical
SIS-like systems in [4]. The x-axis in the right-hand column starts at $0.95$,
which is the expected value of the time a node spends in the firing state
before transitioning to refractory.  Figure 6: PDFs of the sizes and durations
of finite avalanches in the evolved networks at different stages during the
network evolution. The colors correspond to the distinct timepoints shown in
Fig. 2 in the main text. Before the drift (light orange), the distributions
decay faster, while during the drift (red shades) the distributions are heavy-
tailed. The green dashed lines mark the slopes derived for critical SIS-like
systems in [4]. The network size is $N=500000$ and the spontaneous firing rate
is equal to $1/(100N)$, while all other parameters controlling the dynamics
and the network evolution are identical to those used in Fig. 2 in the main
text. Number of runs for each mean degree is originally 10000, but only
avalanches lasting less than $t_{\rm max}=50000$ are considered.

Figure 7: PDFs of evolved networks at time $t=3900000$ for different network
sizes. The networks have been evolved with spontaneous firing rate equal to
$1/(100N)$.

##  VI Illustration of the critical drift

Figure 8 illustrates that the value of the critical mean degree $\langle
k\rangle^{*}_{N=10^{5}}$ increases during the drift phase. In addition, we
observe that during the initial phase, the mean degree of the evolved network
lies below $\langle k\rangle^{*}_{N=10^{5}}$ while during the drift phase, the
system resides at the onset of activity.

Figure 8: The mean degree $\langle k\rangle$ (black circles) and the critical
$\langle k\rangle^{*}_{N=10^{5}}$ (red crosses) in a network evolved with the
plasticity rules. The initial degree is 1 and the parameters are set to
$N=10^{5},l=10^{-3}$, $g=10^{-6}$ and $\eta=\frac{1}{100N}$. The values of
$\langle k\rangle^{*}_{N=10^{5}}$ are obtained by first manipulating $\langle
k\rangle$ of the evolved network as in Fig. 2(e) in the main text and then
finding $\langle k\rangle^{*}_{N=10^{5}}$ based on avalanche durations as in
Fig. 1 of SM.

##  VII Leading eigenvalue as an indicator of criticality

The leading eigenvalue is known to reflect the critical threshold in epidemic
models such as SIS and SIRS (identical to our static IFRI model) on undirected
networks. However, previous proofs for the SIRS/IFRI model have mostly relied
on the assumption that the states of two neighboring nodes are independent.
This assumption leads to a threshold condition $\lambda_{1}=\delta/\beta$ (see
e.g. [6]), which does not accurately predict the critical threshold in our
model.

To improve the accuracy, we no longer assume that the probability of a link
connecting a node in state $X$ to a node in state $Y$ is simply given by the
product of the states’ probabilities. Instead, we derive evolution equations
for the probability of node $i$ being in state $X$ and node $j$ being in state
$Y$, which we denote by $[X_{i}Y_{j}]$. Often, these probabilities depend on
the probabilities of two-link structures, say, the probability of having a
$FI$-link followed by a $II$-link or two $F$-nodes pointing to the same
$I$-node. To avoid rendering the system overly complex, we approximate these
probabilities using the link and node probabilities. For example, given a path
$k\to j\to i$, the probability of the first link being an $FI$-link and the
second link being an $II$-link would be given by
$[F_{k}I_{j}][I_{j}I_{i}]/I_{j}$, where $I_{j}$ denotes the probability of
node $j$ being in state $I$. This type of approximation – called the pair
approximation – is used in [5] to study the critical threshold for
$\lambda_{1}$ in the undirected SIS model, and in [2] to derive a critical
value for the mean degree in the SIRS/IFRI model when the exact network
structure is not known. Here, we assume that the adjacency matrix $\mathbf{A}$
is known and derive a critical threshold for the leading eigenvalue.

With the pair approximation, we obtain the following differential equations:

$\displaystyle\dot{F}_{i}$ $\displaystyle=-\delta
F_{i}+\beta\sum_{j=1}^{N}A_{ji}[F_{j}I_{i}]$ (1) $\displaystyle\dot{R}_{i}$
$\displaystyle=\delta F_{i}-\gamma R_{i}$ (2)
$\displaystyle\dot{[F_{j}I_{i}]}$
$\displaystyle=A_{ji}\Big{(}\beta\frac{[I_{j}I_{i}]}{I_{j}}\sum_{k=1}^{N}A_{kj}[F_{k}I_{j}]+\gamma[F_{j}R_{i}]-\beta\frac{[F_{j}I_{i}]}{I_{i}}\sum_{k\neq
j}^{N}A_{ki}[F_{k}I_{i}]-(\beta+\delta)[F_{j}I_{i}]\Big{)}$ (3)
$\displaystyle\dot{[I_{j}I_{i}]}$
$\displaystyle=A_{ji}\Big{(}\gamma[I_{j}R_{i}]+\gamma[R_{j}I_{i}]-\beta\frac{[I_{j}I_{i}]}{I_{i}}\sum_{k\neq
j}^{N}A_{ki}[F_{k}I_{i}]-\beta\frac{[I_{j}I_{i}]}{I_{j}}\sum_{k=1}^{N}A_{kj}[F_{k}I_{j}]\Big{)}$
(4) $\displaystyle\dot{[F_{j}R_{i}]}$
$\displaystyle=A_{ji}\Big{(}\delta[F_{j}F_{i}]+\beta\frac{[I_{j}R_{i}]}{I_{j}}\sum_{k=1}^{N}A_{kj}[F_{k}I_{j}]-(\gamma+\delta)[F_{j}R_{i}]\Big{)}$
(5) $\displaystyle\dot{[I_{j}R_{i}]}$
$\displaystyle=A_{ji}\Big{(}\gamma[R_{j}R_{i}]+\delta[I_{j}F_{i}]-\gamma[I_{j}R_{i}]-\beta\frac{[I_{j}R_{i}]}{I_{j}}\sum_{k=1}^{N}A_{kj}[F_{k}I_{j}]\Big{)}$
(6) $\displaystyle\dot{[R_{j}I_{i}]}$
$\displaystyle=A_{ji}\Big{(}\delta[F_{j}I_{i}]+\gamma[R_{j}R_{i}]-\beta\frac{[R_{j}I_{i}]}{I_{i}}\sum_{k\neq
j}^{N}A_{ki}[F_{k}I_{i}]-\gamma[R_{j}I_{i}]\Big{)}$ (7)
$\displaystyle\dot{[F_{j}F_{i}]}$
$\displaystyle=A_{ji}\Big{(}\beta\frac{[I_{j}F_{i}]}{I_{j}}\sum_{k=1}^{N}A_{kj}[F_{k}I_{j}]+\beta\frac{[F_{j}I_{i}]}{I_{i}}\sum_{k\neq
j}^{N}A_{ki}[F_{k}I_{i}]+\beta[F_{j}I_{i}]-2\delta[F_{j}F_{i}]\Big{)}$ (8)
$\displaystyle\dot{[R_{j}R_{i}]}$
$\displaystyle=A_{ji}\Big{(}\delta[F_{j}R_{i}]+\delta[R_{j}F_{i}]-2\gamma[R_{j}R_{i}]\Big{)}$
(9) $\displaystyle\dot{[R_{j}F_{i}]}$
$\displaystyle=A_{ji}\Big{(}\beta\frac{[R_{j}I_{i}]}{I_{i}}\sum_{k\neq
j}^{N}A_{ki}[F_{k}I_{i}]+\delta[F_{j}F_{i}]-(\gamma+\delta)[R_{j}F_{i}]\Big{)},$
(10)

where we use the conservation laws

$\displaystyle I_{i}$ $\displaystyle=1-F_{i}-R_{i}$ (11)
$\displaystyle[I_{j}F_{i}]$
$\displaystyle=1-[F_{j}F_{i}]-[F_{j}I_{i}]-[F_{j}R_{i}]-[I_{j}I_{i}]-[I_{j}R_{i}]-[R_{j}F_{i}]-[R_{j}I_{i}]-[R_{j}R_{i}].$
(12)

The onset of activity occurs when the so-called trivial fixed point
$(I_{i}=[I_{j}I_{i}]=1$ for all $i,j$) becomes unstable. From dynamical
systems theory we know that the fixed point loses stability when the leading
eigenvalue of the associated Jacobian matrix crosses from negative to
positive, _i.e._ $\lambda_{1}(\mathbf{J})=0$. The Jacobian associated with
equations 1-10 can be interpreted as a block matrix of size
$(2N+8N^{2})\times(2N+8N^{2})$ consisting of four blocks of size $N\times N$,
16 blocks of size $N\times N^{2}$, 16 blocks of of size $N^{2}\times N$ and 64
blocks of size $N^{2}\times N^{2}$. At the trivial steady state, all blocks
except for $\mathbf{X}$, $\mathbf{Y}$ and $\mathbf{Z}$ are diagonal matrices.
The diagonal matrices corresponding to variables $F-R$ are identity matrices
multiplied by the constant given below. The other diagonal matrices are given
by $\mathrm{diag}(\mathrm{vec}(\mathbf{A}))$ multiplied by the constant given
below:

$\mathbf{J}=\setcounter{MaxMatrixCols}{11}\begin{matrix}&F&R&FF&FI&FR&II&IR&RF&RI&RR\\\
F&-\delta&&&\mathbf{X}\\\ R&\delta&-\gamma\\\ FF&&&-2\delta&\beta\\\
FI&&&&\mathbf{Y}&\gamma\\\ FR&&&\delta&&-\gamma-\delta\\\
II&&&&\mathbf{Z}&&&\gamma&&\gamma\\\
IR&&&-\delta&-\delta&-\delta&-\delta&-\delta-\gamma&-\delta&-\delta&-\delta+\gamma\\\
RF&&&\delta&&&&&-\gamma-\delta&&\\\ RI&&&&\delta&&&&&-\gamma&\gamma\\\
RR&&&&&\delta&&&\delta&&-2\gamma\\\ \end{matrix}$ (13)

To find the stability condition, it suffices to examine a smaller part of the
Jacobian. To see this, we divide the Jacobian into four blocks, where
$\mathbf{B_{1}}$ corresponds to variables $F-[FR]$ and $\mathbf{B_{3}}$ to
$[II]-[RR]$. The eigenvalue equation is then given by

$\begin{bmatrix}\mathbf{B_{1}}&0\\\
\mathbf{B_{2}}&\mathbf{B_{3}}\end{bmatrix}\begin{bmatrix}\mathbf{e_{1}}\\\
\mathbf{e_{2}}\end{bmatrix}=\lambda(\mathbf{J})\begin{bmatrix}\mathbf{e_{1}}\\\
\mathbf{e_{2}}\end{bmatrix},$ (14)

where the eigenvector $\mathbf{e}$ has been split into two parts
$\mathbf{e_{1}}$ and $\mathbf{e_{2}}$. The multiplication results in equations

$\displaystyle\mathbf{B_{1}}\mathbf{e_{1}}$
$\displaystyle=\lambda(\mathbf{J})\mathbf{e_{1}}$ (15)
$\displaystyle\mathbf{B_{2}}\mathbf{e_{1}}+\mathbf{B_{3}}\mathbf{e_{2}}$
$\displaystyle=\lambda(\mathbf{J})\mathbf{e_{2}}.$ (16)

The first equation shows that $\lambda(\mathbf{J})$ must be an eigenvalue of
$\mathbf{B_{1}}$ if $\mathbf{e_{1}}$ is not a zero vector. On the other hand,
if $\mathbf{e_{1}}$ equals zero, $\lambda(\mathbf{J})$ must be an eigenvalue
of $\mathbf{B_{3}}$. Consequently, the eigenvalues of $\lambda(\mathbf{J})$
are given by the eigenvalues of $\mathbf{B_{1}}$ and $\mathbf{B_{3}}$.
However, since the eigenvalues of $\mathbf{B_{3}}$ are non-positive constants,
it suffices to examine the eigenvalues of $\mathbf{B_{1}}$.

With similar logic, we can divide $\mathbf{B_{1}}$ into four blocks where one
block is a zero matrix. Consequently, the eigenvalues of $\mathbf{B_{1}}$ are
given by the blocks formed by variables $F-R$ and $[FF]-[FR]$. The eigenvalues
of the first block are non-positive constants, meaning that the stability
conditions can be found by examining the matrix

$\mathbf{C}=\begin{bmatrix}-2\delta&\beta&0\\\ 0&\mathbf{Y}&\gamma\\\
\delta&0&-\gamma-\delta\end{bmatrix},$ (17)

where
$\mathbf{Y}=\beta\mathbf{M}-(\beta+\delta)\mathrm{diag}(\mathrm{vec}(\mathbf{A}))$
and $\mathbf{M}$ is of the form

$\begin{matrix}&[F_{1}I_{1}]&[F_{2}I_{1}]&...&[F_{1}I_{2}]&[F_{2}I_{2}]&....\\\
[F_{1}I_{1}]&\mathbf{A}_{11}\mathbf{A}_{11}&\mathbf{A}_{21}\mathbf{A}_{11}&...\\\
[F_{2}I_{1}]&&&&\mathbf{A}_{12}\mathbf{A}_{21}&\mathbf{A}_{22}\mathbf{A}_{21}&...\\\
...\\\
[F_{1}I_{2}]&\mathbf{A}_{11}\mathbf{A}_{12}&\mathbf{A}_{21}\mathbf{A}_{12}&...\\\
[F_{2}I_{2}]&&&&\mathbf{A}_{12}\mathbf{A}_{22}&\mathbf{A}_{22}\mathbf{A}_{22}&...\\\
...\\\ \end{matrix}$ (18)

In words, the entry of $\mathbf{M}$ corresponding to row $[F_{i}I_{j}]$ and
column $[F_{k}I_{l}]$ equals one if there is a path $k\to l=j\to i$.

We observe that the nonzero eigenvalues of the matrix $\mathbf{M}$ are equal
to those of the edge adjacency matrix $\mathbf{E}$. The rows and columns of
$\mathbf{E}$ correspond to the existing edges of the network, and
$\mathbf{E}_{ij}=1$ if edges $j$ and $i$ form a directed path. The eigenvalues
are equal because of the symmetric structure of $\mathbf{M}$; if the $i$th
column consists of zeros because $A_{kl}=0$, the same must be true for the
$i$th row. Consequently, $\mathbf{M}$ corresponds to a linear transformation
where the $i$th unit vector (corresponding to the $i$th column) collapses to
the origo. As the $i$th component of all other columns is zero, the
transformation can be depicted in a lower-dimensional space (_i.e._ the $i$th
column and row can be removed) without affecting the nonzero eigenvalues.
Finally, we note that the nonzero eigenvalues of the edge adjacency matrix are
equal to the eigenvalues of the adjacency matrix $\mathbf{A}$. Consequently,
$\lambda_{i}(\mathbf{Y})=\beta\lambda_{i}(\mathbf{A})-\beta-\delta$.

Given the eigenvalues $\lambda_{i}(\mathbf{Y})$, the eigenvalues of the matrix
$\mathbf{C}$ can subsequently be found by setting each eigenvalue
$\lambda_{i}(\mathbf{Y})$ to the diagonal of the middle block in $\mathbf{C}$
and solving the eigenvalues of this modified matrix $\mathbf{C_{D}}$. This is
possible because the part of any eigenvector $\mathbf{e}(\mathbf{C})$ that
$\mathbf{Y}$ modifies must be an eigenvector of $\mathbf{Y}$. To see this, we
look at the eigenvalue equation
$\mathbf{C}\mathbf{e}=\lambda(\mathbf{C})\mathbf{e}$ where
$\mathbf{e}=[\mathbf{e_{1}}\mathbf{e_{2}}\mathbf{e_{3}}]^{T}$. The
multiplication results in equations

$\displaystyle-2\delta\mathbf{e_{1}}+\beta\mathbf{e_{2}}$
$\displaystyle=\lambda\mathbf{e_{1}}$ (19)
$\displaystyle\mathbf{Y}\mathbf{e_{2}}+\gamma\mathbf{e_{3}}$
$\displaystyle=\lambda\mathbf{e_{2}}$ (20)
$\displaystyle\delta\mathbf{e_{1}}+(-\gamma-\delta)\mathbf{e_{3}}$
$\displaystyle=\lambda\mathbf{e_{3}}.$ (21)

Substituting equations 19 and 21 into equation 20, we obtain

$\mathbf{Y}\mathbf{e_{2}}=\Big{(}\lambda(\mathbf{C})+\frac{\beta\delta\gamma}{(\lambda(\mathbf{C})+\gamma+\delta)(-2\delta-\lambda(\mathbf{C}))}\Big{)}\mathbf{e_{2}}.$
(22)

As the first term on the right-hand side is a constant, it must be the case
that $\mathbf{e_{2}}$ is an eigenvector of $\mathbf{Y}$.

The modified matrix $\mathbf{C_{D}}$ can be expressed as a Kronecker product

$\begin{bmatrix}-2\delta&\beta&0\\\
0&\beta\lambda_{i}(\mathbf{A})-\beta-\delta&\gamma\\\
\delta&0&-\gamma-\delta\end{bmatrix}\mathbf{I}^{N^{2}\times N^{2}}.$ (23)

Consequently, the eigenvalues of $\mathbf{C}$ are given by the eigenvalues of
the small $3\times 3$ matrix above. This matrix has the same form as the
matrix in [2], and the signs of its eigenvalues can be solved similarly using
the Hurwitz theorem. Solving for the critical value of $\lambda_{1}$, we
obtain that the trivial steady state loses stability when

$\lambda_{1}=\frac{\delta}{\beta}+\frac{\delta+\gamma/2}{\delta+\gamma}.$ (24)

##  VIII Critical drift in a Watts-Strogatz network

As mentioned in the main text, the leading eigenvalue $\lambda_{1}$ of the
network’s adjacency matrix can change during the critical drift in some highly
structured networks as well as in networks with a significant amount of loops.
In such cases, the value of $\lambda_{1}$ alone does not determine whether or
not the system resides at criticality. To illustrate this, we create a
directed Watts-Strogatz network, a directed ring lattice where each node is
connected to the next node as well as to the node directly after. When this
network is evolved according to the adaptation rules, the drift is clearly
visible in all tracked parameters. Contrary to the initialization with ER
networks, however, $\lambda_{1}$ does not stay constant during the drift [Fig.
9]. Instead, its value stabilizes only after the link removals and random
additions have erased some of the original structure and increased randomness
in the network.

Figure 9: Drift in a network initialized as a directed Watts-Strogatz model.
The graphs are obtained similarly to those in Fig. 2 in the main text, only
the initial network structure differs.

##  IX Correlation of the eigenvector centrality and the probability of firing

We show that a node’s eigenvector centrality correlates with its probability
of being in the firing state $F$. This has previously been shown for several
different epidemic models on undirected networks, including the SIS model
([3]). The proof for our directed IFRI model follows a similar pattern. Let us
denote the probability of node $i$ being in state $F$ at time $t$ by $F_{i}$
and the probability of being in state $R$ by $R_{i}$. For simplicity, the
corresponding vectors are denoted by $F$ and $R$. Assuming that the states of
all nodes are independent of each other, the time evolution of the state
vectors obey the differential equations

$\displaystyle\dot{F}$ $\displaystyle=-\delta F+\beta I\circ F\mathbf{A}$ (25)
$\displaystyle\dot{R}$ $\displaystyle=\delta F-\gamma R,$ (26)

where $\mathbf{A}$ denotes the network’s adjacency matrix and $\circ$ denotes
the Hadamard product. In addition, we know that

$\displaystyle I$ $\displaystyle=1-F-R.$ (27)

At criticality, the derivatives equal zero as we are at a steady state.
Substituting expressions 26 and 27 into Eq. 25, we obtain

$\displaystyle F$
$\displaystyle=\frac{\beta}{\delta}(1-F-\frac{\delta}{\gamma}F)\circ
F\mathbf{A}$ (28)

Assuming that the effect of terms involving $F\circ F$ is negligible, this can
be expressed as

$\displaystyle F$ $\displaystyle=\frac{\beta}{\delta}F\mathbf{A}.$ (29)

For this equation to apply, it must be the case that $F$ is a left eigenvector
of $\mathbf{A}$. If we assume that the network is strongly connected, Perron-
Frobenius theorem for non-negative matrices guarantees that the principal
eigenvector is the only eigenvector where all elements are non-negative.
Hence, as we require $F$ to be non-negative, it must correspond to the
principal eigenvector, which gives the eigenvector centralities of individual
nodes. If the network is not strongly connected, we can apply Perron-Frobenius
theorem individually to each of the network’s strongly connected components.

##  X Driving force of the drift

As discussed in the main text, the drift of the mean degree $\langle k\rangle$
results from the targeted link removals decreasing the leading eigenvalue
$\lambda_{1}$ more efficiently than the random link additions increase it on
average. We verify this by freezing a network at two different times during
the critical drift, first at the beginning of the drift and second after the
value of $\langle k\rangle$ has stabilized. In these frozen networks, we
repeatedly remove or add a link to the network and measure the resulting
change in $\lambda_{1}$ [Fig. 10]. Since a node’s eigenvector centrality
correlates with its probability of being in the firing state, we simulate the
targeted link removals from firing nodes by first choosing a node with
probability proportional to its eigenvector centrality and then randomly
removing one of its incoming links.

We observe that at the beginning of the drift, the targeted link removals
decrease $\lambda_{1}$ more efficiently than random link additions increase it
on average, while at a later stage of the drift, the effects of link removal
and addition have evened out. At this later stage, link removals have a weaker
effect on $\lambda_{1}$ than in ER networks with the same $\langle k\rangle$,
as expected. In addition, the network topology has changed in a way that the
link additions increase $\lambda_{1}$ more than expected in an ER network with
the same mean degree.

Figure 10: Average absolute change in the leading eigenvalue $\lambda_{1}$
after a random link addition (red circles) or targeted link removal (blue
circles) in the evolved networks at two different times during the critical
drift. The triangles show the effect of a link addition/removal in static ER
networks with mean degree equal to that of the evolved network at the
corresponding time point. The markers show the mean effect of 10000 repeated
link additions/removals with the bars displaying the 95% confidence intervals.
Initial $\langle k\rangle$ of the evolved network is 2, and the parameters are
identical to those of Fig. 2 in the main text.

##  XI Dependence of the drift on parameters $\beta$, $\delta$ and $\gamma$

It has been shown in [2] that our model self-organizes to criticality
irrespective of the specific choices of $\beta,\delta$ and $\gamma$
controlling the dynamics. In this section, we investigate how these parameters
affect the critical drift of $\langle k\rangle$.

As argued in the previous section, the increase in $\langle k\rangle$ during
the drift originates from the imbalanced effect of targeted link removals and
random additions on the leading eigenvalue $\lambda_{1}$. If the network
structure is such that the effects of a link addition and removal are already
balanced when the system reaches criticality, we do not expect to observe this
drift. Hence, as $\beta$, $\delta$ and $\gamma$ affect the critical value of
$\lambda_{1}$ that the network self-organizes towards (see Eq. 1 in the main
text), they indirectly influence the network topology at the beginning of the
critical drift, which in turn determines whether or not $\langle k\rangle$
will increase during the drift.

To understand the effect of parameters $\beta,\delta$ and $\gamma$ on the
drift, we first examine how the average effects of link addition and removal
change when the mean degree of a static ER graph increases. As shown in Fig.
11, increasing $\langle k\rangle$ causes the distribution of eigenvector
centralities (corresponding to the elements of the left principal eigenvector)
to concentrate heavily around one value. This means that all nodes have an
increasingly equal probability of firing, which causes the targeted link
removals to reduce to random link removals. Consequently, the average value of
a random link addition and a targeted removal converge to the same value,
$1/N$ [Fig. 12].

Since the drift of $\langle k\rangle$ is driven by the difference in the
expected effect of a link removal and addition, we expect to observe this
drift in ER networks only if the value of $\langle k\rangle$ is relatively low
when the system first reaches criticality. This value can be approximated by
the critical value $\lambda_{1}^{*}$ given by Eq. 1 in the main text (the
approximation works especially well when the initial $\langle k\rangle$ at the
start of the topological evolution is smaller than $\lambda_{1}^{*}$). As
$\lambda_{1}^{*}$ is determined by the parameters $\beta,\gamma$ and $\delta$
controlling the dynamics, the combination of these parameters controls whether
we observe the critical drift of $\langle k\rangle$ in ER networks. Note,
however, that this reasoning applies directly only to ER networks, and a
network with a more varied topology and a more refined local structure may
behave differently.

Figure 11: Histogram of the left principal eigenvector elements for static ER
networks with different mean degrees. The size of the networks is $N=10^{5}$.
The distributions of right principal eigenvector elements behave in a similar
manner. Figure 12: Change in $\lambda_{1}$ after a random link addition or
targeted link removal for ER networks with different mean degrees. The
targeted link removals are simulated by repeatedly choosing a node with
probability proportional to its eigenvector centrality and randomly removing
one of its incoming links. The markers show the mean of 10000 trials and
largely overlay the 95% confidence intervals.

To verify the previous reasoning with simulations, we evolve networks with
different combinations of parameters $\beta$, $\delta$ and $\gamma$ [Fig. 13].
We observe that the slope in which $\langle k\rangle$ increases during the
drift seems to be unaffected by the specific parameter combination, and – as
expected – this slope decreases as $\lambda_{1}^{*}$ increases. We also
observe that the parameter combination has an influence on how much the
stabilized value of $\lambda_{1}$ deviates from the theoretical critical value
$\lambda_{1}^{*}$. This latter observation, however, can be explained by the
network evolution parameters $l$ and $\epsilon$ being finite. As already
discussed in section IIIb, the system self-organizes to criticality in the
limit $l,\epsilon\to 0$, and for finite values of $l$ and $\epsilon$, the mean
degree of the evolved networks is expected to deviate slightly from its
critical value. [2] have shown that for networks with a random ER structure,
the magnitude of this deviation depends on $\delta,\gamma$ and $\langle
k\rangle^{*}$ when $l$ and $\epsilon$ are finite. While this proof does not
directly apply to our case, we expect to observe a similar dependency. Indeed,
as shown in Fig. 14, decreasing the value of $\epsilon$ causes the evolved
$\lambda_{1}$ to decrease towards its theoretical value or even below it,
which – to our understanding – is a result of the finite spontaneous firing
rate pushing the system slightly towards the subcritical regime.

Figure 13: Time evolution of $\lambda_{1}$ and $\langle k\rangle$ for
different combinations of parameters $\beta,\delta$ and $\gamma$. On each row,
the parameter combinations have been chosen to produce a specific value of
$\lambda_{1}^{*}$ (dashed lines) according to Eq. 1 of the main text. We have
first chosen (with no specific logic) the values of $\delta$ and $\gamma$,
after which the parameter $\beta$ is chosen so that Eq. 1 gives the desired
value of $\lambda_{1}^{*}$. Other parameters are identical to those in Fig. 2
of the main text. Figure 14: Time evolution of $\lambda_{1}$ and $\langle
k\rangle$ for different values of $\epsilon$. Each column corresponds to one
parameter combination resulting in $\lambda_{1}^{*}=10$ (dashed line)
according to Eq. 1 of the main text.

##  XII Degree distribution of the evolved networks

Figure 15 shows the distributions of in- and out-degrees at three different
points during the drift.

Figure 15: Distribution of in- and out-degrees at three different times during
the drift

##  XIII Time series of firing rates

Figure 16 shows a few examples of firing rate time series of individual nodes
during a burst of activity along the critical drift. The nodes are chosen
randomly.

Figure 16: Firing rate time series of individual nodes. The parameters are
identical to those in Fig. 2 of the main text.

## References

  * Bak, Tang, and Wiesenfeld [1987] Bak, P., Tang, C., and Wiesenfeld, K., “Self-organized criticality: An explanation of 1/f noise,” Physical Review Letters 59 (1987), 10.1103/PhysRevLett.59.381.
  * Droste, Do, and Gross [2013] Droste, F., Do, A., and Gross, T., “Analytical investigation of self-organized criticality in neural networks,” J. Roy. Soc. Interface 10 (2013), 10.1098/rsif.2012.0558.
  * Goltsev  _et al._ [2012] Goltsev, A., Dorogovtsev, S., Oliveira, J., and Mendes, J. F., “Localization and spreading of diseases in complex networks,” Physical Review Letters 109 (2012), 10.1103/PhysRevLett.109.128702.
  * Larremore  _et al._ [2012] Larremore, D., Carpenter, M., Ott, E., and Restrepo, J., “Statistical properties of avalanches in networks,” Physical Review E 85 (2012), 10.1103/PhysRevE.85.066131.
  * Mata and Ferreira [2013] Mata, A.and Ferreira, S., “Pair quenched mean-field theory for the susceptible-infected-susceptible model on complex networks,” EPL (Europhysics Letters) 103 (2013), 10.1209/0295-5075/103/48003.
  * Prakash  _et al._ [2012] Prakash, B., Chakrabarti, D., Faloutsos, M., Valler, N., and Faloutsos, C., “Threshold conditions for arbitrary cascade models on arbitrary networks,” Knowledge and Information Systems 33 (2012), 10.1109/ICDM.2011.145.


# Topological dynamics and dynamical scaling behavior of vortices in a two-
dimensional XY model

Wei-Kai Qi  Institute of Theoretical Physics, Lanzhou University, Lanzhou
$730000$, China  Yong Chen  Email<EMAIL_ADDRESS>Institute of Theoretical
Physics, Lanzhou University, Lanzhou $730000$, China

###### Abstract

By using topological current theory we study the inner topological structure
of vortices a two-dimensional (2D) XY model and find the topological current
relating to the order parameter field. A scalar field, $\psi$, is introduced
through the topological current theory. By solving the scalar field, the
interaction energy of vortices in a 2D XY model is revisited. We study the
dynamical evolution of vortices and present the branch conditions for
generating, annihilating, crossing, splitting and merging of vortices. During
the growth or annihilation of vortices, the dynamical scaling law of relevant
length in a 2D XY model, $\xi(t)\propto(t-t^{*})^{1/z}$, is obtained in the
neighborhood of the limit point, given the dynamic exponent $z=2$. This
dynamical scaling behavior is consistent with renormalization group theory,
numerical simulations, and experimental results. Furthermore, it is found that
during the crossing, splitting and merging of vortices, the dynamical scaling
law of relevant length is $\xi(t)\propto(t-t^{*})$. However, if vortices are
at rest during splitting or merging, the dynamical scaling law of relevant
length is a constat.

###### pacs:

47.32.C-, 61.72.Cc, 64.70.qj

##  I INTRODUCTION

Vortices play an important role in understanding a variety of problems in
physics. In the 1970s, Kosterlitz and Thouless constructed a detailed and
complete theory of 2D systems Kt . The KT phase transition theory predicts
that vortices pair unbinding will lead to a second-order transition in 2D
systems, such as superfluid films, superconductors, and XY models two ; Dn .
In a 2D XY model, there exist meta-stable states corresponding to vortices,
which are closely bound in pairs below a critical temperature. Above the
critical temperature, the paired vortices unbind and become free. Vortices in
a 2D XY model disrupt the spin alignments even at large distances and
correspond to singularities of the order-parameter field Cl .

In 2D XY model, vortices are topologically stable configurations. It is found
that the high temperature disorder phase with an exponential correlation is a
result of the formation of vortices. The critical temperature at which the KT
transition occurs is, in fact, that at which vortex generation becomes
thermodynamically favorable. At temperatures below this, the system has a
power law correlation. In a low temperature phase, vortices appear in a small
density of tightly bound dipole pairs. With an increase in temperature, vortex
pairs dissociate and become free in the disorder phase. The evolution of
vortices plays an important role in the KT transition of the 2D XY model.

There has been progress on the study of the defects associated with an
n-component vector order parameter field, $\phi(\vec{r},t)$. For the scalar
case (n=1), the defects are domain walls, which are points of the spatial
dimensionality d=1, lines for d=2, planes for d=3, etc. More generally, for
n=d, one has point defects. This leads to vortices for n=d=2. It is
interesting to consider the appropriate form for the point defect densities
when expressed in terms of the vector order parameter field. This has been
carried out by Halperin Ha , and exploited by Liu and Mazenko LM ; however,
their analyses are incomplete Duan00 . In a 2D system, a gauge field-
theoretical formalism has been developed by Kleinert HK . Furthermore, the
gauge theory of topological quantum melting in a $2+1$ dimensional Bose system
was developed by Nussinov et al, and the superfluidity and superconductivity
can arise in a strict quantum field-theoretical setting Z1 .

A topological field theory for topological defects has been developed by Duan
et al Duan . By using the $\phi$-mapping method and topological current
theory, the evolution of the topological defect that relates to singularities
of the order-parameter field, such as the vortex in BEC Duan04  and
superconductivity Duan05 , was studied. In this paper, we will discuss the
topological quantization and evolution of vortex in two-dimensional XY model.
We introduce a scalar field, $\psi$, through the topological current theory.
By solving the scalar field, the Hamiltonian of the interaction of vortices in
the 2D XY model is revisited.

More recently, the dynamical behavior of the 2D XY model following a quench of
the temperature to below the KT critical temperature has been studied Br01 .
During the quench of a dynamic 2D XY model, aging phenomena and dynamic
scaling behavior become extremely important Bo . The assumption of dynamical
scaling for the predicted asymptotic growth law the characteristic length,
$\xi(t)\propto(t/lnt)^{1/2}$ Br02 , which is also characteristic of the
spacing between defects. However, scaling violations were reported. According
to an expansion in $\epsilon=4-d$ using standard field-theoretic
renormalization group theory, it can be shown that the growing length is given
as $\xi(t)\propto t^{1/z}$ for large t, where z is the critical exponent for
equilibrium critical dynamics Jan . This standard theoretical approach also
shows that the result $\xi(t)\propto t^{1/z}$ is independent of the initial
conditions. However, this renormalization group method does not involve the
effects of topological defects, such as vortices in the 2D XY model. In the
recent work, it was shown that, for the specific case of the 2D XY model,
$\xi(t)\propto(t/lnt)^{1/2}$ if free vortices are present, while
$\xi(t)\propto t^{1/2}$ if there are no free vortices present in the initial
state Brprl .

Through our topological current theory for the 2D XY model, the dynamic of
vortices is studied, and the branch conditions for generating, annihilating,
crossing, splitting and merging of vortices are given. During the growth or
annihilation of vortices, the dynamical scaling law of relevant length in the
2D XY model, which is $\xi(t)\propto(t-t^{*})^{1/2}$, can be obtained in the
neighborhood of the limit point. This indicates that vortices in the 2D XY
model are the source of the scaling violations. This dynamical scaling
behavior is consistent with numerical simulations LD  and experimental results
Pa . Furthermore, we have also found that during the crossing, splitting, and
merging of vortices, the dynamical scaling law of relevant length in the 2D XY
model, which is $\xi(t)\propto(t-t^{*})$. It has been shown that why the
relevant length of vortices scales with time as $t^{-1/2}$ for small values of
times and then deviate toward a linear dependence in nematic liquid-crystal
experimental observations Pa2 . However, if vortices are at rest during
splitting and merging, the dynamical scaling law of relevant length is
$\xi=const$. Moreover, it is worthwhile to note that the dynamical scaling law
of vortices, which was deduced from the topological current theory, only
depends on topological properties of the order parameter field.

The organization of this paper is as follows. In Sec. II, we describe the 2D
XY model and the $\phi$-mapping method, and topological current theory is
discussed. The Hamiltonian of vortices in the 2D XY model is revisited in
section III. The evolution of vortices and dynamical scaling discussed in
section IV. Finally, in section V, we summarize our results.

##  II Topological Current in the 2D XY model

The 2D XY model is a system of spins confined to rotate in the plane of the
lattice. The Hamiltonian of the system is give as:

$H_{0}=-J\sum_{<ij>}\textbf{S}_{i}\textbf{S}_{j}=-J\sum_{<ij>}\cos(\theta_{i}-\theta_{j}),$
(1)

where J is the strength of the nearest-neighbor interaction. $\theta_{i}$
denotes the angle of the spin on site i with respect to arbitrary polar
direction in the 2D vector space containing spins.

By using the continuum limit, we can approximate $\cos(\theta_{i}-\theta_{j})$
by the first two terms $1-(\theta_{i}-\theta_{j})^{2}/2$ of the Taylor
expansion. We can express partial derivatives through
$\theta_{i}-\theta_{j}=\partial\theta$ for the two sites i and j. This leads
to the continuum Hamiltonian:

$H=H_{0}-E_{0}=\frac{J}{2}\int d\textbf{r}(\nabla\theta)^{2},$ (2)

where $E_{0}=2JN$ is the energy of the completely aligned ground state of $N$
spins.

Kosterlitz and Thouless suggested that the disordering is caused by
topological defects, such as vortices in two-dimensional XY model, which are
characterized by a mapping from some loop $\Gamma$ in real space onto the
order parameter space. For two-dimensional XY model, this implies

$\oint d\theta=\oint\nabla\theta d\vec{s}=2\pi W,W=0,\pm 1,\ldots,n,$ (3)

where $W$ is the winding number of the vortex. By introducing a 2D XY model,
the local order parameter is defined as:

$\Psi=\Psi_{0}e^{i\theta}=\phi^{1}+i\phi^{2}.$ (4)

Quantized vortices are topological objects associated with topological
properties of the order parameter $\Psi$. It is worth noting that the phase of
the order parameter is undefined at the vortex core. In other words, vortices
correspond to singularities of the order-parameter field. We define the unit
vector field $\vec{n}$ as

$n^{a}=\frac{\phi^{a}}{||\phi||},~{}||\phi||=\sqrt{\phi^{a}\phi^{a}},~{}a=1,2.$
(5)

where $n^{a}n^{a}=1$. From the unit vector, $\vec{n}(x)$, we can construct a
topological current of the order parameter field in the 2D XY model, which
carries the topological information of $\vec{\phi}(x)$:

$j^{k}=\frac{1}{2\pi}\epsilon^{ijk}\epsilon_{ab}\partial_{i}n^{a}\partial_{j}n^{b},\quad
i,j,k=0,1,2.$ (6)

We will see that this current does not vanish at the zero point of
$\vec{\phi}(x)$ or the singularities of the unit field $\vec{n}(x)$. By using
the 2D topological current theorem, Eq. (6) can be rewritten in the compact
form;

$j^{k}=\delta^{2}(\vec{\phi})J^{k}\left(\frac{\phi}{x}\right).$ (7)

It is the important relation between the $\delta-$like topological current of
vortices and the order parameter, $\Psi$, where $J^{k}(\phi/x)$ is the vector
Jacobian of $\vec{\phi}$:

$J^{k}\left(\frac{\phi}{x}\right)=\frac{1}{2}\epsilon^{ijk}\epsilon_{ab}\partial_{i}\phi^{a}\partial_{j}\phi^{b},$
(8)

or

$\displaystyle J^{0}\left(\frac{\phi}{x}\right)$ $\displaystyle=$
$\displaystyle\det\left(\begin{array}[]{cc}\partial_{x}\phi^{1}&\partial_{y}\phi^{1}\\\
\partial_{x}\phi^{2}&\partial_{y}\phi^{2}\end{array}\right),$ (11)
$\displaystyle J^{1}\left(\frac{\phi}{x}\right)$ $\displaystyle=$
$\displaystyle\det\left(\begin{array}[]{cc}\partial_{y}\phi^{1}&\partial_{t}\phi^{1}\\\
\partial_{y}\phi^{2}&\partial_{t}\phi^{2}\end{array}\right),$ (14)
$\displaystyle J^{2}\left(\frac{\phi}{x}\right)$ $\displaystyle=$
$\displaystyle\det\left(\begin{array}[]{cc}\partial_{t}\phi^{1}&\partial_{x}\phi^{1}\\\
\partial_{t}\phi^{2}&\partial_{x}\phi^{2}\end{array}\right).$ (17)

According to the implicit function theorem, the Jacobian’s determinant can be
given as:

$J^{0}\left(\frac{\phi}{x}\right)=J\left(\frac{\phi}{x}\right)\neq 0,$ (18)

The solutions of the zero point of $\vec{\phi}(x)$ can be generally expressed
as:

$x=x_{l}(t),\quad y=y_{l}(t),\quad l=1,2,\ldots,N,$ (19)

which represent $N$ zero points, $\vec{z_{l}}(t)$ (l=1,2,…,N), or a world line
of $N$ vertices in space-time.

With the $\delta$-function theory, $\delta^{2}(\phi)$, can be expanded as:

$\delta^{2}(\phi)=\sum_{l=1}^{N}\frac{\beta_{l}}{|J(\phi/x)_{z_{l}}|}\delta^{2}\left(\vec{r}-\vec{z_{l}}(t)\right)$
(20)

where the positive integer, $\beta_{l}$, is called the Hopf index of map
$x\rightarrow\vec{\phi}$. The meaning of $\beta_{l}$ is that when the point
$\vec{r}$ covers the neighborhood of the zero, $\vec{z_{l}}$, once the vector
field, $\vec{\phi}$, covers the corresponding region for $\beta_{l}$ times.
Using the implicit function theorem and the definition of the vector Jacobian
(Eq. 8), we can find the velocity of the $l$-th defect,

$\displaystyle\vec{v_{l}}$ $\displaystyle=$
$\displaystyle\frac{d\vec{z_{l}}}{dt}=\left[\frac{\vec{J}(\phi/x)}{J(\phi/x)}\right]_{\vec{z_{l}}}$
$\displaystyle\vec{J}\left(\frac{\phi}{x}\right)$ $\displaystyle=$
$\displaystyle\left[J^{1}\big{(}\phi/x\big{)},J^{2}\left(\phi/x\right)\right]$
(21)

The spatial and temporal components of the defect current, $j^{k}$, can be
written as the form of the current and the density of a system of $N$
classical point particles moving in a $(2+1)$-dimensional space-time,

$\vec{j}=\sum_{l=1}^{N}\beta_{l}\eta_{l}\vec{v_{l}}\delta^{2}(\vec{r}-\vec{z_{l}}(t))$
(22)
$j^{0}=\rho=\sum_{l=1}^{N}\beta_{l}\eta_{l}\delta^{2}(\vec{r}-\vec{z_{l}}(t))$
(23)

where $\eta_{l}$ is the Brouwer degree,

$\eta_{l}=\frac{J(\phi/x)}{|J(\phi/x)|}\bigg{|}_{\vec{z_{l}}}=\pm 1.$ (24)

It can clearly be seen that Eq. (22) shows the movement of vertices. The
topological charge of vertices in the 2D XY model are conserved:

$\frac{\partial\rho}{\partial t}+\nabla\vec{j}=0.$ (25)

In addition, there is a constraint of charge neutrality:

$\int\rho d^{2}x=\frac{1}{2\pi}\sum_{l=1}^{N}\beta_{l}\eta_{l}=0,$ (26)

which indicates that vortices in the 2D XY model appear in pairs.

##  III The Hamiltonian of vortices in the 2D XY model revisited

In analogy to the velocity field in a superfluid, the distortion field,
$\vec{u}=\nabla\theta$, carries the topological information of vortices in the
2D XY model. By using Eq. (4) and (5), we can prove that

$\vec{u}=\epsilon_{ab}n^{a}\nabla n^{b},$

and the vorticity is given as:

$\nabla\times\vec{u}=\textbf{e}_{i}(\epsilon^{ijk}\epsilon_{ab}\partial_{j}n^{a}\partial_{k}n^{b})$
(27)

where $\textbf{e}_{k}(k=1,2,3)$ are the base vectors in the Cartesian
coordinate system. Comparing Eq. (27) and Eq. (6), we conclude that

$\nabla\times\vec{u}=2\pi j^{0}\textbf{z}.$ (28)

Therefore, in the 2D XY model, the vorticity of distortion field,
$\nabla\times\vec{u}$, can be expressed in terms of the topological current of
the order parameter field. By comparing the $\delta-$like topological current
Eq. (7) and Eq. (28), we have the important relation between vorticity and the
order parameter field in the 2D XY model:

$\nabla\times\vec{u}=2\pi\delta^{2}(\phi)J\left(\frac{\phi}{x}\right)\textbf{z},$
(29)

From Eq. (29), we see that the vorticity of distortion field,
$\nabla\times\vec{u}$, does not vanish at the zero points of $\Psi$. The
location and direction of the $i$th vortex are determined by the $i$th
singular point, $\vec{z_{l}}(t)$, and the vector Jacobian, $J(\phi/x)$, on
$\vec{z_{l}}(t)$, respectively. In the absence of vorticity, there are no zero
values of the order parameter field, $\delta^{2}(\phi)$ is zero, and Eq. (29)
becomes:

$\nabla\times\vec{u}_{0}=0,$ (30)

which is the condition of irrotationality. Thus, Eq. (29) describes both the
vortex-state and the irrotationality-state. To describe a collection of
vortices at locations $\vec{z_{l}}(t)$, we substitute Eq. (23) into Eq. (28),
leading to:

$\nabla\times\vec{u}=2\pi\sum_{l=1}^{N}\beta_{l}\eta_{l}\delta^{2}(\vec{r}-\vec{z_{l}}(t))\textbf{z}.$
(31)

It can be see that Eq. (31) represents N vortices that are charged with the
topological charge, $W$. $W$ is $\beta_{l}\eta_{l}$, which describes the inner
topological structure of the vortex. To find the solution for Eq. (31), we can
introduce a harmonious scalar field, $\psi$ defining by:

$\displaystyle\nabla^{2}\psi$
$\displaystyle=\partial_{i}\partial_{j}\theta-\partial_{j}\partial_{i}\theta$
(32)
$\displaystyle=\epsilon^{ij}\epsilon_{ab}\partial_{i}n^{a}\partial_{j}n^{b}.$

In the absence of a vortex in the 2D XY model, we can see that the condition,
$\partial_{i}\partial_{j}\theta=\partial_{j}\partial_{i}\theta$, leads to
$\nabla^{2}\psi=0$. Thus, the 2D distortion can be written as,
$\vec{u}=\vec{u}_{0}-\nabla\times(\textbf{z}\psi)$, which contains two parts:
the first part, $\vec{u}_{0}=\nabla\phi$, describe spin waves, while the
second part is associated with vortices. By using 2D topological current
theorem, it can be shown that:

$\nabla^{2}\psi=2\pi\sum_{l=1}^{N}\beta_{l}\eta_{l}\delta^{2}(\vec{r}-\vec{z_{l}}(t)).$
(33)

The scalar field, $\psi$, behaves like the potential due to a set of charged
particles. The solution of Eq. (33) is given as:

$\psi(x)=\sum_{l=1}\beta_{l}\eta_{l}ln(|\vec{r}-\vec{z_{l}}(t)|).$ (34)

The continuum Hamiltonian of the vortices in Eq. (2) can be rewritten as:

$H-\sum_{i}E_{ci}=\frac{J}{2}\int(\nabla\psi)^{2}d^{2}\textbf{x}=-\frac{J}{2}\int\psi\nabla^{2}\psi,d^{2}\textbf{x}$
(35)

where $E_{ci}$ is the core energy of the vortex. Substituting Eq. (33) and
(34) into the Hamiltonian (35), one can obtain:

$\beta H=\beta E_{c}\sum_{i}W_{i}^{2}-\pi
K\sum_{ij}W_{i}W_{j}ln(\vec{z_{i}}(t)-\vec{z_{j}}(t)),$ (36)

which is identical to the Hamiltonian of a 2D Coulomb gas with point charges
of charge $W_{i}=\beta_{i}\eta_{i}$. The constraint given by Eq. (26) is thus
the constraint of charge neutrality.

##  IV The evolution of vortices in the 2D XY model

The zero point of the order parameter field, $\Psi$, which is associated with
the locations of the core of vortices, plays an important role in describing
the evolution of vortices in the 2D XY model. If the Jacobian determinant,
$J^{0}(\phi/x)\neq 0$, we will have the isolated solution of the zeros of the
order parameter field in $2+1$ dimensional space-time. But, when
$J^{0}(\phi/x)=0$, the above results will change in some way and will lead to
the branch process of vortices. We denote one of the vectors Jacobians at the
zero points as $(t^{*},\vec{z_{l}})$. According to the values of the vector
Jacobian at the zero points of the order parameter, there are usually two
kinds of branch points, namely, the limit points and bifurcation points KU .
Each kind corresponds to different cases of branch processes.

Figure 1: Generating and annihilating of vortices pairs. (a) The origin of two
vortices. (b) Two vortices annihilate in collision at the limit point.

###  IV.1 The generation and annihilation of vortices

We explore what will happen to vortices at the limit point
$(t^{*},\vec{z_{l}})$. The limit points are determined by

$\displaystyle
J^{0}\left.\left(\frac{\phi}{x}\right)\right|_{t^{*},\vec{z_{l}}}=0,\quad
J^{1}\left.\left(\frac{\phi}{x}\right)\right|_{t^{*},\vec{z_{l}}}\neq 0$ (37)
$\displaystyle
J^{0}\left.\left(\frac{\phi}{x}\right)\right|_{t^{*},\vec{z_{l}}}=0,\quad
J^{2}\left.\left(\frac{\phi}{x}\right)\right|_{t^{*},\vec{z_{l}}}\neq 0.$ (38)

Considering the condition given by Eq. (37) and making use of the implicit
function theorem, the solution of the zero points of $\phi(x)$ in the
neighborhood of the point ($t^{*},\vec{z_{l}}$) is given as:

$t=t(x),\quad y=y(x),$ (39)

where $t^{*}=t(z_{l}^{1})$. In this case, one can see that:

$\left.\frac{dx}{dt}\right|_{(t^{*},\vec{z_{l}})}=\left.\frac{J^{1}(\phi/x)}{J(\phi/x)}\right|_{(t^{*},\vec{z_{l}})}=\infty,$
(40)

or

$\left.\frac{dt}{dx}\right|_{(t^{*},\vec{z_{l}}}=0.$ (41)

The Taylor expansion of $t=t(x)$ at the limit points $(t^{*},\vec{z_{l}})$ is
given as:

$t-t^{*}=\left.\frac{1}{2}\frac{d^{2}t}{dx^{2}}\right|_{t^{*},\vec{z_{l}}}(x-z_{l}^{1})^{2},$
(42)

which is a parabola in the x-t plane. From this equation, we can obtain two
solutions $x_{1}(t)$ and $x_{2}(t)$, which give two branch solutions (world
lines of the vortices). If

$\left.\frac{d^{2}t}{dx^{2}}\right|_{(t^{*},\vec{z_{l}})}>0,$

we have the branch solutions for $t>t^{*}$ which are related to the origin of
a dipole pair. Otherwise, we have the branch solutions for $t<t^{*}$, which
related to the annihilation of a dipole pair (Fig. 1).

Figure 2: Two vortices collide with different directions of motion at the
bifurcation point in $(2+1)$-dimensional space-time. Two world lines intersect
with different directions at the bifurcation point; i.e., two vortices
encounter at the bifurcation point.

Since the topological current is identically conserved, the topological
charges of these two generated or annihilated vortices must be opposite of one
another at the limit points, such as:

$\beta_{1}\eta_{1}+\beta_{2}\eta_{2}=0.$ (43)

This indicates that the vortices always generate and annihilate in pairs. From
Eq. (42), one also obtains that the velocity of the vortices is infinite when
they are annihilating, which agrees with the result introduced by Bray Bray .
Furthermore, a new result can be obtained that shows that the velocity of the
vortices is infinite when they are generating, which is gained only from the
topology of the order parameter field. For a limit point, it is required that,
$J^{1}(\phi/x)|_{t^{*},\vec{z_{l}}}\neq 0$. A bifurcation point, on the other
hand, must satisfy a more complex condition. This case will be discussed in
the following.

###  IV.2 Encountering, Splitting and Merging of vortices

Now, let us turn to consider the case in which the restrictions on the zero
point $(t^{*},\vec{z_{l}})$ are given by:

$J^{k}\left.\left(\frac{\phi}{x}\right)\right|_{(t^{*},\vec{z_{l}})}=0,\quad
k=0,1,2,$ (44)

which imply an important fact that the function relationship between t and x
or y is not unique in the neighborhood of the bifurcation point
$(t^{*},\vec{z_{l}})$. This fact is easily seen from

$\frac{dx}{dt}=\left.\frac{J^{1}(\phi/x)}{J(\phi/x)}\right|_{t^{*},\vec{z_{l}}},~{}~{}\frac{dy}{dt}=\left.\frac{J^{2}(\phi/x)}{J(\phi/x)}\right|_{t^{*},\vec{z_{l}}},$
(45)

which, under Eq. (44), directly shows the indefiniteness of the direction of
the integral curve of Eq. (45) at $(t^{*},\vec{z_{l}})$. For this reason, the
point $(t^{*},\vec{z_{l}})$ is called a bifurcation point of the orientation
order parameter.

Figure 3: Vortices with the same direction of motion. (a) Two world lines
tangentially contact; i.e., two vortices tangentially encounter at the
bifurcation point. (b) Two world lines merge into one world line; i.e., two
vortices merge into one vortex at the bifurcation point. (c) One world line
resolves into two world lines; i.e., one vortex splits into two vortices at
the bifurcation point.

As we know, at the bifurcation point, $(t^{*},\vec{z_{l}})$, the rank of the
Jacobian matrix, $\left[\partial\phi/\partial\textbf{x}\right]$, is $1$ in the
2D vector order parameter. With the aim of finding the different directions of
all branch curves at the bifurcation point, we assume:

$\left.\frac{\partial\phi^{1}}{\partial y}\right|_{t^{*},\vec{z_{l}}}\neq 0.$
(46)

From the implicit function theorem, there is one function relation:

$y=f(x,t).$ (47)

According to the $\phi$-mapping theory, the Taylor expansion of the solution
of the zeros of the order parameter field in the neighborhood of
$(t^{*},\vec{z_{l}})$ can be expressed as: Duan

$A(x-z_{l}^{1})^{2}+2B(x-z_{l}^{1})(t-t^{*})+C(t-t^{*})^{2}+\cdots=0,$ (48)

which leads to

$A\left(\frac{dx}{dt}\right)^{2}+2B\frac{dx}{dt}+C=0,$ (49)

and

$C\left(\frac{dt}{dx}\right)^{2}+2B\frac{dt}{dx}+A=0,$ (50)

where A, B, and C are constants determined by the order parameter. The
solutions of Eq. (49) or Eq. (50) give different directions for the branch
curves (world line of the vortices) at the bifurcation point. There are four
possible cases which demonstrate the physical meaning of the bifurcation
points.

Figure 4: Two important cases of Eq. (53). (a) One world line resolves into
three world lines; i.e., one vortex resolves into three vortices at the
bifurcation point. (b) Three world lines merge into one world line; i.e.,
three vortices merge into one vortex at the bifurcation point.

Case 1 $(A\neq 0)$: For $\Delta=4(B^{2}-AC)>0$, from Eq. (49) we get two
different motion directions of the core of the vortex given as:

$\left.\frac{dx}{dt}\right|_{1,2}=\frac{-B\pm\sqrt{B^{2}-AC}}{A},$ (51)

which is shown in Fig. 2, where two world lines of two vortices intersect with
different directions at the bifurcation point. This shows that two vortices
encounter and then depart at the bifurcation point.

Case 2 $(A\neq 0)$: For $\Delta=4(B^{2}-AC)=0$, form Eq. (49), we obtain only
one motion direction of the core of the vortex given as:

$\left.\frac{dx}{dt}\right|_{1,2}=-\frac{B}{A},$ (52)

which includes three important cases (Fig. 3). These three cases are: (i) Two
world lines tangentially contact; i.e., two vortices tangentially encounter at
the bifurcation point. (ii) Two world lines merge into one world line; i.e.,
two vortices merge into one vortex at the bifurcation point. (iii) One world
line resolves into two world lines; i.e., one vortex splits into two vortices
at the bifurcation point.

Case 3 $(A=0,C\neq 0)$: For $\Delta=4(B^{2}-AC)=0$, from Eq.(49) we have:

$\left.\frac{dt}{dx}\right|_{1,2}=\frac{-B\pm\sqrt{B^{2}-AC}}{C}=0,~{}-\frac{2B}{C}.$
(53)

There are two important cases (Fig. 4): (i) One world line resolves into three
world lines; i.e., one vortex splits into three vortices at the bifurcation
point. (ii) Three world line merge into one world line; i.e., three vortices
merge into one vortex at the bifurcation point.

Case 4 (A=C=0): Equations (49) and (50) give respectively:

$\frac{dx}{dt}=0,\quad\frac{dt}{dx}=0.$ (54)

This case shows that two world lines intersect normally at the bifurcation
point, which is similar to case 3. It is no surprise that both parts of Eq.
(54) are correct because they give the slope coefficients of two different
curves at the same point $(t^{*},\vec{z_{l}})$.

The remaining components of $dy/dt$ can be calculated from a function relation
of Eq. (47):

$\frac{dy}{dt}=\frac{\partial f}{\partial x}\frac{\partial x}{\partial
t}+\frac{\partial f}{\partial t}.$ (55)

The above solutions reveal the evolution of the vortices. The topological
structure of the vortices is detailed in the neighborhood of the bifurcation
points of the order parameter field. Besides the encountering of vortices,
i.e., two vortices encountered at and then depart from the bifurcation point
along different branch curves (Fig. 2 and Fig. 3(a)), splitting and merging of
vortices are also included. When multicharged vortices pass the bifurcation
point, it may split into several vortices along different branch curves (Fig.
3(c), Fig. 4(a)). On the other hand, several vortices can merge into one
vortex at the bifurcation point (Fig. 3(b), Fig. 4(c)). As before, since the
topological current of the vortices is identically conserved, the sum of
topological charges of the final vortices must be equal to that of the initial
vortices at the bifurcation point, which is given as (for fixed l):

$\sum_{f}\beta_{lf}\eta_{lf}=\sum_{i}\beta_{li}\eta_{li}.$ (56)

This indicates that vortices with a higher value of Burgers vector can evolve
to the lower value of Burgers vector, or that vortices with a lower value of
Burgers vector can evolve to a the higher value of Burgers vector through the
bifurcation process. Furthermore, we see that the generation, annihilation,
and bifurcation of vortices are not gradual changes, but that they start at a
critical value of a parameter, i.e., a sudden change. It is important to note
that further bifurcations are possible during the evolution of vortices
besides the cases studied in this work. It is necessary need to assume that
the terms in the Taylor series considered above vanish and that the expansion
of the field is dominated by higher order terms.

###  IV.3 Dynamical scaling law of vortices in the 2D XY model

In the neighborhood of the limit point, we denote the scale length $l=\nabla
x$. The growth velocity or annihilation velocity of vortices is $v=l/\nabla t$
Xu . From Eq. 42), one can obtain the scaling law as:

$v\varpropto(t-t^{*})^{-1/2}.$ (57)

It can be seen that $E_{k}\varpropto(t-t^{*})^{-1}$. In the 2D XY model, the
growth or annihilation is parameterized in terms of a relevant characteristic
length $\xi(t)$, which is also characteristic of the mean vortex-antivortex
separation distance. From Eq. (57), we find that the relevant length,
$\xi(t)$, obeys the following:

$\xi(t)\sim(t-t^{*})^{1/2}.$ (58)

Equation (58) is the dynamic scaling law of the vortex-antivortex pairs. The
low temperature equilibrium phase has essentially no vortex pairs, and
$\xi(t)$ is infinite. The relationship between the relevant length and the
metastable vortex density below the critical temperature is given as:

$\xi(t)=\sqrt{\frac{1}{\rho_{v}}}.$ (59)

then the number of vortices that satisfies the power law $N\propto t^{-1}$.
The dynamic scaling law of the vortex-antivortex pairs is consistent with
renormalization group theory and results from recent numerical simulations
using Langevin equations and Monte Carlo methods Brprl ; LD . The specially
prepared nematic liquid-crystal system Pa , which is developed, exhibits 2D XY
behavior and can be utilized to obtains the dynamical scaling behavior of
vortices (Fig. 5). This result is consistent with the expansion in
$\epsilon=4-d$ using standard field-theoretic renormalization group theory.
This standard theory shows that the growing length, $\xi(t)\propto t^{1/z}$
for large $t$, where $z$ is the critical exponent for equilibrium critical
dynamics Jan . This standard theoretical approach also shows that the result
$\xi(t)\propto t^{1/z}$ is independent of the initial conditions. However,
this renormalization group method does not involve the effects of topological
defects, such as vortices in the 2D XY model.

Figure 5: Number of defects versus time. Experimental and simulation results
are shown. The results of simulation were the average of 20 quenches LD , for
which $T^{*}=5.0$ initially and is suddenly reduced to $T^{*}=0.8$ ($\square$)
and $T^{*}=0.9$ ($\triangle$). The experimental results, which are obtained by
Pargellis et al Pa , were averaged over nine experimental runs ($\bullet$).
Solid lines are the expected $N\propto t^{-1}$ scaling.

From the assumptions of dynamical scaling law, it has been predicted that the
asymptotic growth law of the characteristic length is given by
$\xi(t)\propto(t/lnt)^{1/2}$ Br02 ; however, scaling violations were reported.
In a recent work, it was shown that the specific case of the 2D XY model is
given as $\xi(t)\propto(t/lnt)^{1/2}$ if free vortices are present, while
$\xi(t)\propto t^{1/2}$ if there are no free vortices present in the initial
state Brprl . Since the topological current is identically conserved, there
are no free vortices during the growth or annihilation process.

In the neighborhood of the bifurcation point, we denote scale length $\nabla
x=l$. From Eqs. (51)-(53), we can then obtain the growth or annihilation
velocity of the vortices, which is given as

$v\propto const.$ (60)

The approximation asymptotic relation of is then

$\xi(t)\propto(t-t^{*}).$ (61)

It can be seen that this is the reason why the relevant length of the vortices
scales with time as $t^{-1/2}$ for small values of $t$ and then deviates
towards a linear dependence in nematic liquid-crystal experimental
observations Pa2 .

From Eq. (54), one can obtain

$\xi(t)=const,\quad v=0.$ (62)

It can be seen that that vortices are relatively at rest when $\xi(t)=const$.

During the growth or annihilation of vortices, the dynamical scaling law of
the relevant length in the 2D XY model is given as
$\xi(t)\propto(t-t^{*})^{1/2}$. This indicates that the vortices in the 2D XY
model are the source of the scaling violations. During the crossing, splitting
and merging of vortices, the dynamical scaling law of relevant length in the
2D XY model is given as $\xi(t)\propto(t-t^{*})$. However, if vortices are at
rest during merging or splitting, the dynamical scaling law of relevant length
is $\xi(t)=const$. Moreover, it is worth noting that the dynamical scaling law
of vortices depends only on topological properties of the order parameter
field.

##  V Conclusions

In summary, we have studied the inner structure and evolution of vortices in
the 2D XY model by making using of the $\phi$-mapping topological current
theory. A scalar field, $\psi$, is introduced through the topological current
theory. By solving the scalar field, the interaction energy of vortices in the
2D XY model is revisited.

By using the $\phi$-mapping current theory, the densities of vortices in terms
of the order parameter field in the 2D XY model are obtained directly from the
definition of the topological charges of the vortices. The inner topological
structure of the charge of the vortices, which is characterized by the Hopf
index and the Brouwer degree, was obtained. By using the 2D topological
current theorem, the Hamiltonian of the interaction vortices in the 2D XY
model was revisited. This result is identical to the Hamiltonian of a 2D
Coulomb gas with point charges of charge $W_{i}=\beta_{i}\eta_{i}$.

Furthermore, we have studied the evolution of vortices in the 2D XY model and
concluded that there are crucial cases of branch processes in the evolution of
vortices when $J^{0}(\phi/x)=0$ and $\phi=0$, i.e., $\eta_{l}$ is indefinite.
It is important to note that according to the Landau-Ginzburg approach, there
is instability in the 2D XY model at the zero point of the order parameter
fields in the low temperature phase. This indicates that vortices are unstable
in the low temperature phase. In the high temperature phase, the zero point of
the order parameter field is stable. Due to the branch condition, the vortices
generate or annihilate at the limit points and encounter, split, or merge at
the bifurcation points of the order parameter field. This result also
indicates that the velocity of the vortices is infinite when they are being
annihilated or generated, which is obtained only from the topological
properties of the order parameter field in the 2D XY model. The scaling law of
relevant length, which is given as $R(t)\propto(t-t^{*})^{1/2}$, can be
obtained in the neighborhood of the bifurcation point during the growth or
annihilation of vortices. During the crossing, splitting and merging of
vortices, the dynamical scaling law of relevant length is
$\xi(t)\propto(t-t^{*})$. However, if vortices are at rest during merging or
splitting, the dynamical scaling law of relevant length is $\xi=const$. The
dynamical scaling law of vortices only depends on topological properties of
the order parameter field.

##  Appendix A Topological Current theory

We study a two-component vector order parameter,
$\vec{\phi}=(\phi^{1},\phi^{2})$, over the base manifold, M (in this paper
$M=R^{2}\bigotimes R$):

$\Psi=\phi^{1}+i\phi^{2},\quad\phi^{a}=\phi^{a}(x,y,t),\quad a=1,2.$ (63)

Quantized defects are topological objects associated with topological
properties of the order parameter $\Psi$. It is worth noting that the phase of
the order parameter is undefined at the vortex core. In other words, vortices
correspond to singularities of the order-parameter field. Let us define the
unit vector field, $\vec{n}$, as:

$n^{a}=\frac{\phi^{a}}{||\phi||},\quad||\phi||=\sqrt{\phi^{a}\phi^{a}},\quad
a=1,2.$ (64)

where, $n^{a}n^{a}=1$. From the unit vector, $\vec{n}(x)$, we can construct a
topological current of the order parameter field in the 2D XY model, which
carries the topological information of $\vec{\phi}(x)$:

$j^{k}=\frac{1}{2\pi}\epsilon^{ijk}\epsilon_{ab}\partial_{i}n^{a}\partial_{j}n^{b},\quad
i,j,k=0,1,2.$ (65)

It can be seen that this current does not vanish only at the zero point of
$\vec{\phi}(x)$, or the singularities of the unit field $\vec{n}(x)$.
Obviously, the current given by Eq.(65) is identically conserved:

$\partial_{k}j^{k}=0.$ (66)

Suppose there is a vortex located at $z_{i}$. The topological charge of the
defect is defined by the Gauss map, $n:\partial\sum_{i}\rightarrow S^{1}$:

$W(\phi,z_{i})=\frac{1}{2\pi}\oint_{\partial\sum_{i}}\epsilon_{ab}n^{a}dn^{b}.$
(67)

Using Stokes’ theorem in the exterior differential form, one can deduce that:

$W(\phi,z_{i})=\frac{1}{2\pi}\oint_{\sum_{i}}\epsilon_{ab}\epsilon^{ij}\partial_{i}n^{a}\partial_{j}n^{b}d^{2}x.$
(68)

##  Appendix B Two-dimensional topological current Theorem

By using the $\phi$-mapping method, we can prove that the topological current
given by Eq. (65) is a $\delta$-like current:

$j^{k}=\delta(\vec{\phi})J^{k}\left(\frac{\phi}{x}\right).$ (69)

From the above formula, it can be seen that the topological current does not
vanish only at the zero points of the vector field $\vec{\phi}(x)$.

Substituting Eq. (64) into Eq. (65) and considering that:

$\partial_{k}n^{a}=\frac{\partial_{k}\phi^{a}}{||\phi||}-\phi^{a}\partial_{k}\left(\frac{1}{||\phi||}\right).$
(70)

we have

$j^{k}=\frac{1}{2\pi}\epsilon^{ijk}\epsilon_{ab}\partial_{i}\phi^{a}\partial_{j}^{b}\frac{\partial}{\partial\phi^{a}}\frac{\partial}{\partial\phi^{b}}(\ln(||\phi||)),$
(71)

If we define the Jacobian as:

$\epsilon^{ab}J^{k}\left(\frac{\phi}{x}\right)=\epsilon^{ijk}\partial^{j}\phi^{a}\partial^{k}\phi^{b},$
(72)

and, by virtue of the Laplacian relation in $\phi$ space GS

$\triangle_{\phi}\left(\ln(||\phi||)\right)=2\pi\delta^{2}(\vec{\phi}),$ (73)

where

$\triangle_{\phi}=\frac{\partial}{\partial\phi^{a}}\frac{\partial}{\partial\phi^{a}}$

is the two-dimensional Laplacian operator in $\phi$ space, we obtain a
$\delta$-function like current. This current is given by:

$j^{k}=\delta(\vec{\phi})J^{k}\left(\frac{\phi}{x}\right).$ (74)

From the above formula, one can see that the topological current does not
vanish only at the zero points of the vector field $\vec{\phi}(x)$. Therefore,
it is essential to investigate the solutions of $\vec{\phi}(x)=0$.

###### Acknowledgements.

We would like to thank Y. X. Liu, T. Zhu, and Y. S. Duan for their helpful
discussion. Y.C. was supported by the SRF for ROCS, SEM, and by the
Fundamental Research Fund for Physics and Mathematics of Lanzhou University.

## References

  * (1) N. D. Mermin, Phys. Rev 176, 250 (1968); N. D. Mermin and H. Wagner, Phys. Rev. Lett. 22, 1133 (1966); V. L. Berezinskii, Sov. Phys.-JETP, 32, 493 (1970); H. E. Stanley and T. A. Kaplan, Phys. Rev. Lett 17, 913 (1966); J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1181 (1973).
  * (2) K. J. Strandburg, Rev. Mod. Phys. 60, 161 (1988).
  * (3) D. R. Nelson, Defects and Geometry in Condensed Matter Physics, (Cambridge University Press, Cambridge, 2002).
  * (4) P. M. Chaikin and T. C. Lubensky, Principles of condensed matter physics, (Cambridge University Press, Cambridge, 1995).
  * (5) B. I. Halperin, Physics of Defects, edited by R. Balian et al. (North-Holland, Amsterdam, 1981).
  * (6) F. Liu and G. F. Mazenko, Phys. Rev. B 46, 5963 (1992).
  * (7) Y. S. Duan, H. Zhang, and L. B. Fu, Phys. Rev. E 59, 528 (1999).
  * (8) H. Kleinert, Gauge Theory in Condensed Matter, (World Scientific, Singapore, 1989).
  * (9) J. Zaanen, Z. Nussinov and S. I. Mukhin, Ann. Phys(New York) 310, 181(2004); V. Cvetkovic, Z. Nussinov, S. Mukhin, and J. Zaanen, Europhys. Lett 81, 27001 (2007).
  * (10) Y. S. Duan and S. L. Zhang, Int. J. Eng. Sci. 29, 1593 (1991); G. H. Yang and Y. S. Duan, Int. J. Theor. Phys. 37, 2371 (1998); Y. S. Duan and H. Zhang, Phys. Rev. E 60, 2568 (1999).
  * (11) Y. S. Duan and H. Zhang, Eur. Phys. J. D 5, 47 (1999).
  * (12) Y. S. Duan, H. Zhang, and S. Li, Phys. Rev. B 58, 125 (1998).
  * (13) A. J. Bray, Adv. Phys. 43, 357 (1994).
  * (14) B. Zheng, Int. J. Mod. Phys. B 12, 1419 (1998).
  * (15) A. J. Bary and A. D. Rutenberg, Phys. Rev. E 49, R27 (1994).
  * (16) H. K. Janssen, B. Schaub, and B. Schmittmann, Z. Phys.B 73, 539 (1989).
  * (17) A. J. Bray, A. J. Briant, and D. K. Jervis. Phys. Rev. Lett 84, 1503 (2000).
  * (18) R. Loft and T. A. DeGrand, Phys. Rev. B 35, 8528 (1987).
  * (19) A. N. Pargellis, S. Green, and B. Yurke, Phys. Rev. E 49, 4250 (1994).
  * (20) A. N. Pargellis, P. Finn, J. W. Goodby, P. Panizza, B. Yurke and P. E. Cladis, Phys. Rev. A 46, 7765 (1992).
  * (21) M. Kubicek and M. Marek, Computational Methods in Bifurcation Theory and Dissipative Structures, (Springer, New York, 1983).
  * (22) A. J. Bray, Phys. Rev. E 55, 5297 (1997).
  * (23) T. Xu, Phys. Rev. E 72, 036303 (2005).
  * (24) I. M. Gel’fand and G. E. Shilov, Generalized Function, (Academic, New York, 1964) Vol. 1.