text:
\section{Introduction}
\label{sec:intro}

The best known examples of topological magnetic solitons are magnetic bubbles
or skyrmions observed in abundance in ferromagnetic films with easy-axis
anisotropy \cite{malozemoff,hubert}.
The experimental situation is less clear in the case of ferromagnets with easy-plane
anisotropy. The relevant topological structures are theoretically predicted
to be half skyrmions or vortices, with a logarithmically divergent energy
that may inhibit production of an isolated vortex on an infinite film.
Thus the early studies of magnetic vortex dynamics
had been mostly theoretical \cite{huber82,mertens99}
drawing on various analogies with related work
on ferromagnetic bubbles \cite{thiele74},
with vortex dynamics in classical fluids and superfluids
\cite{batchelor,saffman,donnelly},
as well as with the dynamics of interacting electric charges in a uniform
magnetic field.

The situation has changed dramatically in recent years.
It has been realized that a disc-shaped magnetic element, with a diameter
of a few hundred of nanometers, provides an excellent geometry for the
realization of a magnetic vortex configuration.
In particular, the exchange energy is finite on a finite element while the
magnetostatic field vanishes everywhere except at the vortex core.
As a result, the vortex is actually the lowest energy magnetic state in a
disc-shaped element. In other words, interest in the vortex stems from
the fact that it is a nontrivial magnetic state which can, nevertheless,
be spontaneously created in magnetic elements \cite{raabe00}.

It is then natural to ask whether nontrivial magnetic states other than
the single vortex may play an important role in the dynamics
of magnetic elements \cite{shigeto02,castano03}.
An answer to this question comes from a somewhat unlikely direction.
Recent experiments have shown a peculiar dynamical behavior of vortices
and magnetic domain walls when these are probed by external magnetic fields.
Vortices may switch their polarity under the influence of a very weak
external magnetic field of the order of a few mT \cite{neudert05,waeyenberge06}.
The same switching phenomenon was observed by passing an a.c. electrical current
through a magnetic disc \cite{yamada07}.
Since the polarity of the vortex contributes to its topological structure,
the switching process clearly implies a discontinuous (topologically forbidden)
change of the magnetic configuration.
This is certainly a surprise especially because the external field is rather weak.
The key to this phenomenon is the appearance of vortex pairs which are
spontaneously created in the vicinity of existing vortices \cite{waeyenberge06,hertel07}.
The creation of topological excitations (vortex pairs) by alternating external fields
had been anticipated by an early study based on collective coordinates
\cite{pokrovskii85}.

In this paper we study vortex-antivortex pairs (VA pairs) which are nontrivial
magnetic states that play an important role in the dynamics
of magnetic elements.
However, unlike a single vortex, a VA pair is a localized object whose energy
remains finite even on an infinite film. It is then reasonable to expect that
the essential features of the dynamics of VA pairs can be understood in
the infinite-film approximation which is adopted in the following.
A brief summary of the relevant dynamical equations and related topological
structures is given in Section~\ref{sec:model}.
In Section~\ref{sec:kelvin} we study a VA pair in which the vortex and the antivortex
carry the same polarity. Such a pair is shown to undergo translational
Kelvin motion analogous to that observed in fluid dynamics \cite{papanicolaou99}.
In Section~\ref{sec:rotate} we study a VA pair in which the vortex and the antivortex
carry opposite polarities. Such a pair is shown to behave as a rotating
vortex dipole \cite{komineas07} because its topological structure
is substantially different than that of the pair in Kelvin motion.

The preceding results are combined in Section~\ref{sec:switching}
to demonstrate that a rotating vortex dipole may be annihilated by a quasi-continuous
process in spite of its nontrivial topological structure.
In particular, no energy barrier has to be overcome in contrast to the usual
expectations for topological solitons.
This opens the possibility for switching mechanisms between topologically
distinct states in ferromagnets.
The possibility to change the topological structure leads to a dramatic change
in the magnetization dynamics as a VA pair is created or annihilated.
Such a pair annihilation process lies in the heart of the counter-intuitive
vortex polarity switching event that was observed in magnetic elements
\cite{neudert05,waeyenberge06,yamada07}.

Some concluding remarks are summarized in Section~\ref{sec:conclusion}.
Finally, an Appendix is devoted to a brief description of the dynamics
of interacting electric charges in a uniform magnetic field,
which also exhibits most of the peculiar features of the dynamics of VA pairs.

\section{The model}
\label{sec:model}

A ferromagnet is characterized by the magnetization $\bm{m} = (m_1, m_2, m_3)$
measured in units of the constant saturation magnetization $M_s$.
Hence $\bm{m}$ is a vector field of unit length, $\bm{m}^2 = m_1^2+m_2^2+m_3^2 = 1$,
but is otherwise a nontrivial function of position and time
$\bm{m} = \bm{m}(\bm{r},t)$ that satisfies the rationalized Landau-Lifshitz (LL)
equation
\begin{equation}  \label{eq:lle}
\frac{\partial\bm{m}}{\partial t} = \bm{m} \times \bm{f}, \quad
\bm{f} \equiv \Delta\bm{m} - q\, m_3\, \bm{\hat{e}}_{\rm 3}, \quad
\bm{m}^2 = 1.
\end{equation}
Here distances are measured in units of the exchange length
$\ell_{\rm ex} = \sqrt{A/2\pi M_s^2}$, where $A$ is the exchange constant,
and the unit of time is $\tau_0 \equiv 1/(4\pi\gamma M_s)$ where $\gamma$ is the
gyromagnetic ratio.
Typical values are $\ell_{\rm ex} \sim 5{\rm nm}$
and $\tau_0 \sim  10{\rm ps}$ which set the scales for the phenomena
described by Eq.~(\ref{eq:lle}).
To complete the description of the LL equation we note that we consider
ferromagnetic materials with uniaxial anisotropy.
Then $\bm{\hat{e}}_{\rm 3}$ in Eq.~(\ref{eq:lle}) is a unit vector along
the symmetry axis and the dimensionless parameter $q \equiv K/2\pi M_s^2$,
where $K$ is an anisotropy constant, measures the strength of anisotropy.
In particular, $q$ is taken to be positive throughout this paper,
a choice that corresponds to easy-plane ferromagnets.

An important omission in Eq.~(\ref{eq:lle}) is the demagnetizing field produced
by the magnetization itself \cite{malozemoff,hubert}. However, in the limit
of a very thin film, the effect of the demagnetizing field is thought to
amount to a simple additive renormalization of the anisotropy constant
\cite{gioia97}.
Also note that we may perform the rescalings $\sqrt{q} \bm{r} \to \bm{r}$ and
$q t \to t$ which further renormalize the units of space and time
discussed earlier and lead to a completely rationalized LL equation
where we may set $q=1$ without loss of generality.
With this understanding, all calculations presented in this paper are based
on a two-dimensional (2D) restriction of Eq.~(\ref{eq:lle}),
i.e., $\bm{r} = (x,y)$ and $\Delta = \partial^2/\partial x^2 + \partial^2/\partial y^2$, while
$q$ is set equal to unity without further notice.

The effective field $\bm{f}$ in Eq.~(\ref{eq:lle}) may be derived from
a variational argument:
\begin{equation}  \label{eq:energy}
\bm{f} = -\frac{\delta E}{\delta \bm{m}}\,; \quad
E = \frac{1}{2} \int{\left[ (\bm{\nabla}\bm{m})^2 + m_3^2 \right]\, dx dy}
\end{equation}
where $E$ is the conserved energy functional.
A standard Hamiltonian form is obtained by resolving the constraint $\bm{m}^2=1$
through, say, the spherical parameterization
\begin{equation}  \label{eq:thetaphi}
m_1 + i\, m_2 = \sin\Theta\;e^{i\Phi}, \quad m_3 = \cos\Theta.
\end{equation}
The LL equation is then written as
\begin{equation}  \label{eq:llesymplectic}
\frac{\partial\Phi}{\partial t} = \frac{\delta E}{\delta \Pi}\,, \quad
\frac{\partial\Pi}{\partial t} = -\frac{\delta E}{\delta \Phi}
\end{equation}
where $\Pi\equiv \cos\Theta$ is the canonical momentum conjugate to the
azimuthal angle $\Phi$. Taking into account the specific form
of the energy in Eq.~(\ref{eq:energy}) or
\begin{equation}  \label{eq:energy1}
E = \frac{1}{2} \int \left[ (\bm{\nabla}\Theta)^2 + \sin^2\Theta\, (\bm{\nabla}\Phi)^2
    + \cos^2\Theta\right]\, dx dy
\end{equation}
the Hamilton equations (\ref{eq:thetaphi}) yield
\begin{eqnarray}  \label{eq:llethetaphi}
\sin\Theta\, \frac{\partial \Phi}{\partial t} & = &
   \Delta\Theta + [1 - (\bm{\nabla}\Phi)^2]\,\cos\Theta\,\sin\Theta\,, \nonumber \\
\sin\Theta\, \frac{\partial \Theta}{\partial t} & = &
                  - \bm{\nabla}\cdot (\sin^2\Theta\, \bm{\nabla}\Phi).
\end{eqnarray}
Another useful parameterization is obtained through the stereographic variable $\Omega$:
\begin{equation}  \label{eq:omega}
\Omega = \frac{m_1 + i\,m_2}{1+m_3}\,;\quad
m_1 + i\,m_2 = \frac{2\Omega}{1+\overline{\Omega}\Omega}\,,\quad
m_3 = \frac{1-\overline{\Omega}\Omega}{1+\overline{\Omega}\Omega}
\end{equation}
in terms of which the LL equation reads
\begin{equation}  \label{eq:lleomega}
i\,\frac{\partial \Omega}{\partial t} + \Delta\Omega + \frac{1-\overline{\Omega}\Omega}{1+\overline{\Omega}\Omega}\,\Omega =
\frac{2\overline{\Omega}}{1+\overline{\Omega}\Omega}\,(\bm{\nabla}\Omega\cdot\bm{\nabla}\Omega)
\end{equation}
where $\overline{\Omega}$ is the complex conjugate of $\Omega$.
Equations (\ref{eq:lle}),(\ref{eq:llethetaphi}), and (\ref{eq:lleomega})
are three equivalent versions of the LL equation and may be used at convenience
depending on the specific calculation considered.

A key quantity for describing both topological and dynamical properties of the 2D
LL equation is the local topological vorticity $\gamma = \gamma(x,y,t)$ defined
from \cite{papanicolaou91,komineas96}:
\begin{equation}  \label{eq:vorticity1}
\gamma = \epsilon_{\alpha\beta}\, \partial_\alpha\Pi\partial_\beta\Phi
       = \epsilon_{\alpha\beta}\, \sin\Theta\, \partial_\beta\Theta \partial_\alpha\Phi
       = \frac{1}{2} \epsilon_{\alpha\beta}\, (\partial_\beta\bm{m} \times \partial_\alpha\bm{m}) \cdot \bm{m},
\end{equation}
where the usual summation convention is invoked for the repeated indices $\alpha$
and $\beta$, which take over two distinct values corresponding to the two spatial
coordinates $x$ and $y$, and $\epsilon_{\alpha\beta}$ is the 2D antisymmetric tensor.
In particular, one may consider the total topological vorticity $\Gamma$
and the Pontryagin index or skyrmion number ${\cal N}$ defined from
\begin{equation}  \label{eq:skyrmionnumber}
\Gamma = \int{\gamma\, dx dy},\quad  {\cal N} = \frac{\Gamma}{4 \pi}.
\end{equation}
A naive partial integration using Eq.~(\ref{eq:vorticity1})
yields $\Gamma=0={\cal N}$ for all magnetic configurations for which
such an integration is permissible.
However, nonvanishing values for $\Gamma$ and ${\cal N}$ are possible
and are topologically quantized.
Specifically, for field configurations that approach a constant (uniform)
magnetization at spatial infinity, the skyrmion number ${\cal N}$ is quantized
according to ${\cal N} = 0,\pm 1,\pm 2, \ldots$. Half integer values are also
possible in the case of field configurations with more complicated structure
at infinity such as half skyrmions or vortices (see below).

The local topological vorticity $\gamma$ is also important for an unambiguous
definition of conservation laws in the LL equation.
Hence the linear momentum (impulse) $\bm{P} = (P_x, P_y)$ is defined from
\begin{equation}  \label{eq:linmomentum}
P_x = - \int{y\gamma\, dx dy},\quad P_y = \int{x\gamma\, dx dy},
\end{equation}
while the angular momentum (impulse) is given by
\begin{equation}  \label{eq:angmomentum}
L = \frac{1}{2}\int{\rho^2\gamma\, dx dy}
\end{equation}
where $\rho^2 = x^2+y^2$. Since detailed discussions of these conservation laws have
already appeared in the literature \cite{papanicolaou99,papanicolaou91,komineas96}
we simply note here that analogous conservation laws were defined as moments
of ordinary vorticity in fluid dynamics \cite{batchelor,saffman}.

\begin{figure}
\epsfig{file=vortex.eps,width=0.5\linewidth}
   \caption{Profile of a single vortex calculated numerically.
The complete solution is given by Eq.~(\ref{eq:vortex}).
   }
 \label{fig:vortex}
\end{figure}

We first search for static (time independent) solutions of the LL equation
which may be obtained by omitting time derivatives in Eq.~(\ref{eq:llethetaphi})
and further introducing the axially symmetric ansatz $\Theta=\theta(\rho)$
and $\Phi= \kappa(\phi-\phi_0)$,
where $\rho$ and $\phi$ are the usual cylindrical coordinates
$(x =\rho\cos\phi,\, y =\rho\sin\phi)$, $\kappa=\pm 1$
will be referred to as the vortex number,
and $\phi_0$ is an arbitrary constant phase reflecting the
azimuthal invariance. The resulting ordinary differential equation
for the amplitude $\theta=\theta(\rho)$ is solved numerically
with standard boundary condition $\theta(\rho \to \infty) = \pi/2$ and the result
is shown in Fig.~\ref{fig:vortex}.
The corresponding magnetization is then given by
\begin{equation}  \label{eq:vortex}
m_1 + i\, m_2 = \sin\theta\, e^{i\kappa(\phi-\phi_0)},
\quad m_3 = \lambda\cos\theta,
\end{equation}
where $\lambda=\pm 1$ will be called the polarity.
The total energy is accordingly reduced to
\begin{equation}  \label{eq:energy2}
E = \frac{1}{2} \int\limits_0^\infty\left[\left(\frac{\partial\theta}{\partial\rho}\right)^2
 + \frac{\sin^2\theta}{\rho^2} + \cos^2\theta \right]\, (2\pi\rho d\rho).
\end{equation}
where the centrifugal (second) term  is logarithmically divergent
for the assumed boundary condition at spatial infinity.
However, the anisotropy energy given by the last term is finite and is
actually predicted to be
\begin{equation}  \label{eq:virial1}
E_{\rm a} = \frac{1}{2} \int\limits_0^\infty \cos^2\theta\, (2\pi\rho d\rho) = \frac{\pi}{2}
\end{equation}
by a careful derivation of a suitable virial relation \cite{komineas98}.
Finally, the total vorticity $\Gamma$ and the skyrmion number ${\cal N}$
are calculated from Eq.~(\ref{eq:skyrmionnumber}) to be
\begin{equation}  \label{eq:vortextopo}
\Gamma = -2\pi\kappa\lambda\,,\quad {\cal N} = -\frac{1}{2}\,\kappa\lambda\,,
\end{equation}
where the vortex number $\kappa=\pm 1$ and the polarity $\lambda=\pm 1$
may be taken in any combination.
Thus, we must consider four possibilities; namely, a vortex that comes
in two varieties ($\kappa=1, \lambda=\pm 1$) and thus ${\cal N}= \mp 1/2$,
and an antivortex which also comes
in two varieties ($\kappa=-1, \lambda=\pm 1$) and thus ${\cal N}= \pm 1/2$.
In all cases the calculated static solution is a topological soliton
with half integer skyrmion number, in contrast to ordinary skyrmions
(such as magnetic bubbles) which carry integer ${\cal N}$.

Our main aim in the continuation of this paper is to search for nontrivial
solutions that may combine a vortex and an antivortex (VA pair)
in a way that the total energy is finite.
We leave aside for the moment the LL equation and construct model VA pairs
in terms of the basic single-vortex configuration described above.
This is easily accomplished by invoking the stereographic variable $\Omega$
of Eq.~(\ref{eq:omega}) to write a single ($\kappa,\lambda$) vortex located at the
origin of coordinates as
\begin{equation}  \label{eq:omega-thetaphi}
\Omega = \frac{\sin\theta}{1+\lambda \cos\theta}\, e^{i\kappa(\phi-\phi_0)}
\end{equation}
where $\theta=\theta(\rho)$ is the vortex profile taken from Fig.~\ref{fig:vortex}
or simply the model profile defined from $\cos\theta=1/\cosh\rho$ and
$\sin\theta=\tanh\rho$.
We now produce two replicas of the basic vortex to describe a pair of vortices
by the product ansatz
\begin{equation}  \label{eq:productomega}
\Omega = \Omega_1\,\Omega_2
\end{equation}
where $\Omega_1$ is configuration (\ref{eq:omega-thetaphi}) applied for
$(\kappa,\lambda) = (\kappa_1,\lambda_1)$ and the origin displaced to, say,
$(x,y) = (-d/2,0)$
while $\Omega_2$ is a $(\kappa_2,\lambda_2)$ vortex located around $(x,y) = (d/2,0)$.
The skyrmion number of this configuration is given by
\begin{equation}  \label{eq:addS}
{\cal N} = -\frac{1}{2}\, (\kappa_1\lambda_1 + \kappa_2\lambda_2)
\end{equation}
and its energy is finite only if we restrict attention to vortex-antivortex (VA) pairs
$(\kappa_1=-\kappa_2)$.
For definiteness we choose $\kappa_1=1$ and $\kappa_2=-1$ and thus the skyrmion number
\begin{equation}  \label{eq:subtractS}
{\cal N} = -\frac{1}{2}\, (\lambda_1-\lambda_2)
\end{equation}
depends on the polarities $\lambda_1$ and $\lambda_2$.
Now a VA pair with equal polarities $(\lambda_1=\lambda_2 = \pm 1)$ is topologically
trivial $({\cal N} = 0)$. In contrast, a VA pair with opposite polarities
is topologically equivalent to a skyrmion
(${\cal N} = 1$, for $\lambda_1 = -1$ and $\lambda_2=1$)
or an antiskyrmion (${\cal N} = -1$, for $\lambda_1 = 1$ and $\lambda_2=-1$).
In all three cases configuration (\ref{eq:productomega}) carries finite energy
because its overall phase cancels out at spatial infinity where the magnetization
approaches the uniform configuration $\bm{m}=(1,0,0)$
modulo an overall azimuthal rotation which depends on the choice
of individual  phases $\phi_0$ in the ansatz (\ref{eq:productomega}).
This explains, in particular, why VA pairs are characterized by an integer
skyrmion number (${\cal N}=0,\pm 1$).

Needless to say, the VA pairs constructed above are not solutions of the
LL equation. Yet an interesting picture arises when configuration (\ref{eq:productomega})
is used as an initial condition in the complete equation (\ref{eq:lleomega}).
A topologically trivial (${\cal N}=0$) VA pair undergoes Kelvin motion
in which the vortex and the antivortex initially located at a relative distance $d$
along the x axis move in parallel along the y axis with nearly constant velocity.
In contrast, a topologically nontrivial (${\cal N}=\pm 1$) VA pair
undergoes rotational motion around a fixed guiding center at nearly constant
angular velocity. In both cases the main trajectories are decorated by Larmor-type
oscillations \cite{volker94} which are tamed when the relative distance between the
vortex and the antivortex is large.
The effect of the polarity of vortex pairs has been studied experimentally
in patterned ferromagnetic ellipses \cite{buchanan05}. Two vortices were created an ellipse
and it was found that their dynamics depended on their relative polarity,
as is indicated by Eq.~(\ref{eq:subtractS}).

In the following two sections we examine the two cases in turn.
In particular, we aim at constructing true steady-state solitary waves that
describe VA pairs in pure Kelvin motion for ${\cal N}=0$ and pure rotational motion
for ${\cal N}=\pm 1$.

\section{Kelvin motion}
\label{sec:kelvin}

As the title of this section suggests, VA pairs in Kelvin motion were
originally studied in the context of ordinary fluid dynamics \cite{batchelor,saffman}.
A further analogy exists with the 2D motion of an electron-positron pair
interacting via the Coulomb potential and placed in a uniform magnetic field
perpendicular to the plane.
If the electron and the positron are initially at rest their guiding centers
will move along two parallel straight lines while the actual trajectories
will display the familiar Larmor oscillations. However, when the electron and the
positron are given a common initial velocity such that the Coulomb force is
exactly balanced by the magnetic force,
both the guiding centers and the actual positions of the
charges will move steadily along parallel lines, even though the two sets
of trajectories do not coincide.
The resulting special configuration may be thought of as a peculiar
electron-positron bound state in steady translational motion (see our Appendix).

It is thus reasonable to expect that a solitary wave exists in a 2D easy-plane
ferromagnet which describes a VA pair that proceeds rigidly (without Larmor oscillations)
in a direction perpendicular to the line connecting the vortex and the antivortex,
probably because the mutual force is exactly balanced by a topological
``Magnus force''.
The actual existence of such a solitary wave was established in Ref.~\cite{papanicolaou99}
whose main result is briefly reviewed in the remainder of this section.

As it turns out, when the relative distance is large, the sought after solitary wave
resembles in its gross features  the model VA pair of Eq.~(\ref{eq:productomega})
applied for, say, $\kappa_1=-\kappa_2=1$ and $\lambda_1=\lambda_2=1$
(thus ${\cal N}=0$). We may then invoke this model to motivate some important
asymptotic results valid for large $d$.
For instance, the local topological vorticity $\gamma$ is then peaked around
the two points $(x,y)=(-d/2,0)$ and $(x,y)=(d/2,0)$ with weights $-2\pi$ and $2\pi$
respectively, corresponding to the total vorticities $\Gamma$ of the individual
vortex and antivortex.
Then the impulse defined from Eq.~(\ref{eq:linmomentum}) yields $P_x=0$,
thanks to reflexion symmetry, while $P = P_y \sim 2\pi d$. One may also
invoke the Derrick-like scaling relation applied to the extended energy functional
$F=E- v P$:
\begin{equation}  \label{eq:derrick1}
v P = \int m_3^2\, dx dy = 2\,E_{\rm a}
\end{equation}
where $E_{\rm a}$ is the total anisotropy energy of the VA pair.
For large $d$ this energy approaches the sum of the anisotropy energies
of the individual  vortex and antivortex, each given by $\pi/2$ in view of the
virial relation (\ref{eq:virial1}). Hence,
$\int{m_3^2\, dx dy} \sim 2 (\pi/2+\pi/2) = 2\pi$ and $vP \sim 2\pi$.
To summarize,
\begin{equation}  \label{eq:estimates1}
P \sim  2\pi d, \quad v P \sim 2\pi, \quad v \sim\frac{1}{d}.
\end{equation}
One may further consider the familiar group-velocity relation
\begin{equation}
v = \frac{dE}{dP}
\end{equation}
in which we insert the estimate $v\sim 2\pi/P$ to obtain an elementary
differential equation for $E$ whose integral is
\begin{equation}  \label{eq:kelvinpairenergy}
E \approx 2\pi\,\ln(P/P_0)
\end{equation}
where $P_0$ is an integration constant that cannot be fixed by the present leading-order
argument. Nevertheless, Eq.~(\ref{eq:kelvinpairenergy}) provides the essence
of the energy-momentum dispersion at large relative distance $d$
or small velocity $v\sim 1/d$ and hence large momentum $P\sim 2\pi d$.

To obtain an accurate numerical solution it is convenient to work
with the stereographic variable $\Omega$ of Eq.~(\ref{eq:omega}).
Then the LL Eq.~(\ref{eq:lleomega}) restricted to a solitary wave in
rigid motion with constant velocity $v$ along, say, the y axis reads
\begin{equation}  \label{eq:lleomega1}
-iv\,\frac{\partial \Omega}{\partial y} + \Delta\Omega
  + \frac{1-\overline{\Omega}\Omega}{1+\overline{\Omega}\Omega}\,\Omega =
\frac{2\overline{\Omega}}{1+\overline{\Omega}\Omega}\,(\bm{\nabla}\Omega\cdot\bm{\nabla}\Omega)
\end{equation}
and is supplemented by the boundary condition $\Omega \to 1$ at spatial infinity
where the magnetization approaches a uniform configuration.
Once a solution $\Omega = \Omega(x,y;v)$ of Eq.~(\ref{eq:lleomega1}) is obtained
for a specific value of the velocity $v$, the sought after solitary wave
is given by $\Omega(x,y-vt;v)$.

\begin{figure}
\epsfig{file=kelvin02.eps,width=0.9\linewidth}
   \caption{Snapshot of a topologically trivial (${\cal N}=0$)
VA pair in Kelvin motion, illustrated through the ($m_1,m_2$) projection
of the magnetization $\bm{m} = (m_1,m_2,m_3)$ on the plane of the film
(left panel) as well as level contours of $m_3$ (right panel).
The contour levels $m_3=0.1,0.3,0.5,0.7,0.9$ are shown. The outer curve corresponds
to $m_3=0.1$, while the two smallest circles are the contours for $m_3=0.9$
and surround the two vortex centers.
The pair moves along the $y$ axis with velocity $v=0.2$ at a
calculated relative distance $d=4.97$ (along the x axis), energy $E=20$,
and impulse $P=30$.
   }  \label{fig:kelvin02}
\end{figure}

\begin{figure}
\epsfig{file=kelvin095.eps,width=0.9\linewidth}
   \caption{Solitary wave in Kelvin motion with $v=0.95$ along
the y axis (conventions as in Fig.~\ref{fig:kelvin02}).
Note that the vortex-antivortex character is lost ($d=0$) and the wave
is a lump with no apparent topological features.
The calculated energy and impulse are $E=18$ and $P=17$.
   }  \label{fig:kelvin095}
\end{figure}

Equation (\ref{eq:lleomega1}) was solved numerically via a Newton-Raphson
iterative algorithm for velocities in the range $0.1 \leq v < 0.99$.
Note that $v=1$ is the familiar magnon velocity (in rationalized units)
and provides an upper bound for the existence of a solitary wave
in rigid motion. On the contrary, there is no lower bound for the velocity -- the
restriction $v \geq 0.1$ was dictated only by numerical expedience.
Now, the calculated solitary wave is illustrated in Fig.~\ref{fig:kelvin02}
for the relatively low velocity $v=0.2$ and does indeed describe a VA pair
(with ${\cal N} = 0$) at a relative distance $d=4.97\sim 1/v$,
as anticipated by the discussion of model VA pairs in Section~\ref{sec:model}
and the heuristic asymptotic analysis earlier in this section. But we are now in
a position to carry out an accurate calculation practically throughout
the allowed range $0< v < 1$, to discover that there exists a characteristic velocity
$v_0 \sim 0.78$ above which the vortex-antivortex character is lost ($d\sim 0$)
and the solitary wave becomes a lump with no apparent topological features,
as illustrated in Fig.~\ref{fig:kelvin095} for $v=0.95$.

\begin{figure}
\epsfig{file=figure04.eps,width=0.80\linewidth}
   \caption{Energy $E$ and impulse $P$ as functions of velocity $v$
(left panel) and $E$ vs $P$ dispersion (right panel)
for a solitary wave in Kelvin motion.
The dotted line was calculated from the asymptotic dispersion
(\ref{eq:kelvinpairenergy2}) applied for $Q_0=4.7$.
   }  \label{fig:kelvin_energy}
\end{figure}

The existence of a characteristic velocity $v_0$ becomes
apparent when we calculate energy $E$ and impulse $P$ as functions of
$v$, as shown in Fig.~\ref{fig:kelvin_energy}.
Note that both $E$ and $P$ develop a minimum at a common velocity $v=v_0=0.78$.
As a result, the energy vs impulse dispersion shown in
Fig.~\ref{fig:kelvin_energy}
develops a cusp at a point $(E_0,P_0)$ that corresponds to the
values of $E$ and $P$ at $v=v_0$.
Thus the calculated family of solitary waves consists of two branches.
Branch I consists of VA pairs that propagate with
velocities in the range  $v < v_0$.
The corresponding branch in the dispersion of Fig.~\ref{fig:kelvin_energy}
approaches the asymptotic dispersion (\ref{eq:kelvinpairenergy})
for large $P$ (or $v\to 0$). Indeed, an excellent fit of the data
is obtained by Eq.~(\ref{eq:kelvinpairenergy}) if we choose the subleading
constant according to $\ln P_0 \sim 1/4$.
Branch II consists of lumps with no apparent topological features
which propagate with velocities in the range $v_0 < v < 1$.
The corresponding branch in the dispersion of Fig.~\ref{fig:kelvin_energy}
is accurately described in the asymptotic ($v\to 1$) region, where
$P$ again becomes large, by
\begin{equation}  \label{eq:kelvinpairenergy2}
E = P\, \left(1 + \frac{Q_0^2}{2P^2} + \ldots\right)
\end{equation}
with $Q_0\approx 4.7$. Actually this asymptotic dispersion can be derived by
showing that in the limit $v \to 1$ the 2D Landau-Lifshitz equation
reduces to what is a modified Kadomtsev-Petviashvili (KP)
equation \cite{papanicolaou99}.

A cusp in the energy vs impulse dispersion
occurred previously in a calculation of vortex rings
in a model superfluid by Jones and Roberts \cite{jones82}.
The same authors together with Putterman later argued that the solitary waves
that correspond to branch II are actually  unstable \cite{jones86}.
Interestingly, a dispersion with a cusp appears also in the much simpler
problem of electron-positron motion in a uniform magnetic field,
where the motion that corresponds to branch II is also unstable
(see our Appendix).
Hence, while a stability analysis has not yet been carried out for the solitary
waves described in this section, it is reasonable to expect that the lumps
of branch II may be unstable.
But there is every reason to believe that the Kelvin motion of the VA pairs
of branch I is indeed stable.

\section{Rotational motion}
\label{sec:rotate}

The possibility of topologically nontrivial (${\cal N}=\pm 1$) VA pairs
in steady rotational motion was recently examined by one of us \cite{komineas07}.
Again, one may invoke the model VA pair of Section~\ref{sec:model}
to understand some important features of the rotational motion
at large distance $d$.
For definiteness we consider a VA pair defined by Eq.~(\ref{eq:productomega})
with $\kappa_1=-\kappa_2=1$ and $\lambda_1=-\lambda_2=-1$, thus ${\cal N} = 1$.
The local topological vorticity $\gamma$ is then peaked around the positions of the
vortex and the antivortex, now with weight equal to $2\pi$ in both cases.
Then, for large $d$, the angular momentum defined by Eq.~(\ref{eq:angmomentum})
is estimated to be $L \sim \frac{1}{2} 2(2\pi)(\frac{d}{2})^2=\frac{\pi}{2} d^2$.
One may also consider the Derrick-like scaling relation applied to the extended
energy functional $F=E-\omega L$:
\begin{equation}  \label{eq:derrick2}
 \omega L = \frac{1}{2}\int{m_3^2\, dx dy} = E_{\rm a}
\end{equation}
where $E_{\rm a}$ is the total anisotropy energy of a VA pair which approaches
asymptotically $E_{\rm a} \sim (\frac{\pi}{2}+\frac{\pi}{2}) = \pi$, and hence $\omega L = \pi$,
because the anisotropy energy of a single vortex is equal to  $\pi/2$
according to Eq.~(\ref{eq:virial1}).
To summarize,
\begin{equation}  \label{eq:estimates2}
L \sim \frac{\pi}{2}\, d^2\,,\quad \omega L \sim \pi\,,\quad \omega \sim \frac{2}{d^2}.
\end{equation}
One may further employ the familiar relation
\begin{equation}  \label{eq:dedl}
\omega = \frac{dE}{dL}
\end{equation}
in which we insert the estimate $\omega\sim\pi/L$ to obtain an elementary differential
equation for $E$ whose integral is
\begin{equation}  \label{eq:rotatepairenergy}
E \approx \pi\ln(L/L_0)
\end{equation}
where $L_0$ is an integration constant that cannot be fixed by the present
leading-order argument. Nevertheless, Eq.~(\ref{eq:rotatepairenergy})
provides the essence of the energy vs angular momentum
dispersion for large relative distance $d$ or small angular frequency
$\omega \sim 2/d^2$ and hence large angular momentum $L\sim\frac{\pi}{2}\,d^2$.

While the preceding heuristic asymptotic analysis is very useful for
understanding some basic aspects of a rotating VA pair, it does not give us any clue
concerning  the fate of the pair at small vortex-antivortex separation.
In principle, such a question could be settled by solving numerically
the analog of Eq.~(\ref{eq:lleomega1}) for a solitary wave rotating
at constant angular frequency $\omega$:
\begin{equation}  \label{eq:llerotateomega}
i\omega\,\epsilon_{\alpha\beta}\,x_\alpha\partial_\beta\Omega + \Delta\Omega
  + \frac{1-\overline{\Omega}\Omega}{1+\overline{\Omega}\Omega}\,\Omega =
\frac{2\overline{\Omega}}{1+\overline{\Omega}\Omega}\,(\bm{\nabla}\Omega\cdot\bm{\nabla}\Omega).
\end{equation}
A numerical solution could again be attempted by an iterative Newton-Raphson
algorithm. Actually, in this case, we found it more convenient to employ
a relaxation algorithm to derive approximate numerical solutions as stationary
points of the extended energy functional $F=E-\omega L$ \cite{komineas07}.

\begin{figure}
\epsfig{file=rotate700.eps,width=0.9\linewidth}
  \caption{Snapshot of a topologically nontrivial (${\cal N}=1$)
VA pair in rotational motion (conventions as in Fig.~\ref{fig:kelvin02},
with solid lines in the right panel corresponding to positive values
of $m_3$ and dashed lines to negative ones).
The pair rotates around a fixed guiding center taken at the origin
of coordinates, with angular velocity $\omega=0.06$ and calculated
relative distance $d=5.3$, energy $E=21$, and angular impulse $L=64$.
}  \label{fig:rotate700}
\end{figure}

\begin{figure}
\epsfig{file=rotate200.eps,width=0.9\linewidth}
  \caption{${\cal N}=1$ VA pair in rotational motion with
angular velocity $\omega=0.18$ and calculated
relative distance $d=1.7$, energy $E=15$, and angular impulse $L=11$.
The only difference from Fig.~\ref{fig:rotate700} is that the
overall size of the pair is now reduced.
}  \label{fig:rotate200}
\end{figure}

The calculated configuration for $\omega=0.06$ is shown in Fig.~\ref{fig:rotate700}
and does indeed correspond to a topologically nontrivial (${\cal N}=1$)
rotating VA pair consisting of a vortex with negative polarity $(\kappa,\lambda)=(1,-1)$
and an antivortex with positive polarity $(\kappa,\lambda)=(-1,1)$,
as anticipated by the general discussion of Section~\ref{sec:model}.
In particular, for this relatively small value of $\omega$,
the calculated distance $d$ in the rotating pair is relatively large
$(d=5.3)$, so is the angular momentum $(L=64)$;
in rough agreement with the asymptotic estimates of Eq.~(\ref{eq:estimates2}).

The calculation was repeated for a larger value of angular velocity
$(\omega=0.18)$ to find that both the relative distance $(d=1.7)$
and the angular momentum $(L=11)$ are reduced to smaller values.
But the general structure of the solution shown in Fig.~\ref{fig:rotate200}
for $\omega=0.18$ remains basically the same as that for $\omega=0.06$,
except that the overall size of the VA pair is reduced.
For yet larger values of $\omega$ the angular momentum $L$ tends to vanish.
This trend is apparent in Fig.~\ref{fig:rotate_energy} which illustrates the
dependence of $E$ and angular momentum $L$ as functions of angular velocity $\omega$,
as well as the $E$ vs $L$ dispersion. For large values of $L$ the above
dispersion exhibits logarithmic dependence, as predicted by the asymptotic
result of Eq.~(\ref{eq:rotatepairenergy}). But the most important feature
of the calculated dispersion is that energy approaches the finite value
$E=4\pi$ as $L\to 0$, while the corresponding rotating VA pair becomes vanishingly
small. This result is of central importance for our main argument and will
be analyzed in some detail in the remainder of this section.

\begin{figure}
\epsfig{file=figure07.eps,width=0.80\linewidth}
   \caption{Energy $E$ and angular impulse $L$ as functions of
angular velocity $\omega$ (left panel) and $E$ vs $L$ dispersion (right panel)
for an ${\cal N}=1$ VA pair in rotational motion.
   }  \label{fig:rotate_energy}
\end{figure}

The best way to describe a vanishing VA pair, in the limit $L\to 0$,
is to invoke yet another model configuration through the stereographic
variable
\begin{equation}  \label{eq:VAskyrmion}
\Omega = \frac{\overline{\zeta}-\frac{d}{2}}{\overline{\zeta}+\frac{d}{2}}\,,\quad
\zeta = x+i y\,,\;\; \overline{\zeta} = x - i y\,,
\end{equation}
where the constant $d$ is taken to be real for simplicity.
Now, configuration (\ref{eq:VAskyrmion}) reaches a finite value $\Omega=1$
at spatial infinity and the corresponding magnetization is uniform:
\begin{equation}  \label{eq:bc}
\bm{m} = (1,0,0)\quad as \quad |\zeta| \to \infty
\end{equation}
which is an appropriate boundary value (modulo a constant azimuthal rotation)
for easy-plane anisotropy discussed here.
Because of (\ref{eq:bc}) the skyrmion number is expected to be an integer
and is actually computed to be ${\cal N}=1$ by a direct application
of Eq.~(\ref{eq:skyrmionnumber}). Furthermore, the spin configuration derived from
(\ref{eq:VAskyrmion}) is an exact solution of the LL equation, if we neglect
anisotropy, with exchange energy
\begin{equation}  \label{eq:skyrmionenergy}
E=E_{\rm e} = \frac{1}{2} \int \left(\nabla\bm{m}\cdot\nabla\bm{m}\right)\, dx dy = 4\pi\,,
\end{equation}
for any $d$, as discussed long time ago by Belavin and Polyakov \cite{belavin75}.

We now return to the main line of argument by noting that configuration
(\ref{eq:VAskyrmion}) may also be thought of as a topologically nontrivial
(${\cal N}=1$) model VA pair consisting of a
$(\kappa,\lambda)=(1,-1)$ vortex and a $(\kappa,\lambda)=(-1,1)$ antivortex
located at a distance $d$ apart, in close analogy with the model VA pair
constructed in Section~\ref{sec:model}.
Although (\ref{eq:VAskyrmion}) is not an exact solution in the presence of
anisotropy, it provides a good model for the behavior of a vanishing VA pair
in the limit $L\to 0$. Anisotropy sets a distance scale $R\sim 1/\sqrt{q} = 1$
(in rationalized units) beyond which a physically acceptable configuration
must reach the uniform value (\ref{eq:bc}).
In a sense, this property is shared by configuration (\ref{eq:VAskyrmion})
because it becomes uniform almost everywhere when $d \ll 1$ while it retains
its topological structure as long as $d$ remains finite.
The strict limit $d\to 0$ is not uniform.
If taken naively, all topological structure appears to be lost and both
energy $E$ and skyrmion number ${\cal N}$ appear to vanish.
However, if integrals are performed before taking the $d\to 0$ limit,
$E$ approaches $4\pi$ without encountering an energy barrier
while ${\cal N}=1$ for all $d$.
Clearly the limit $d\to 0$ creates a singular point which hides all
topological structure. The main point is the claim that a similar
situation arises in the calculated rotating VA pair in the limit $L\to 0$,
as discussed further in Section~\ref{sec:switching}.

This section is completed with a comment concerning the manner in
which the $L\to 0$ limit is reached. Our numerical data as well as
virial relation (\ref{eq:derrick2}) are consistent with a linear dispersion
\begin{equation}  \label{eq:rotatedispersion}
E \approx 4\pi + \frac{1}{2}\,L
\end{equation}
in the limit $L\to 0$, which implies a finite value of the angular velocity
$\omega = dE/dL=1/2$ in the limit of a vanishing VA pair,
and thus an upper limit $\omega_{\rm max} \approx  1/2$.

\section{Vortex core switching}
\label{sec:switching}

Applied for long time intervals, our relaxation algorithm
revealed some tendency for instability of a rotating VA pair,
probably due to radiation effects analogous to those expected
for a pair of rotating electric charges discussed in the Appendix.
However, the basic features of rotating VA pairs studied in the
preceding section persist over sufficiently long time intervals
and are thus relevant for practical applications.
In fact, in a realistic ferromagnet, some dissipation is always
present and can be modeled by introducing Gilbert damping
in the LL equation through the replacement
\begin{equation}  \label{eq:gilbert}
\frac{\partial \bm{m}}{\partial t} \to \frac{\partial \bm{m}}{\partial t} +
 \alpha \left(\bm{m}\times \frac{\partial \bm{m}}{\partial t}\right)
\quad \rm{or}\quad
i\,\frac{\partial \Omega}{\partial t} \to (i-\alpha)\,\frac{\partial \Omega}{\partial t}
\end{equation}
in Eq.~(\ref{eq:lle}) or Eq.~(\ref{eq:lleomega}), respectively,
where $\alpha$ is a dissipation constant.

The dynamics of a topologically nontrivial (${\cal N}=\pm 1$)
VA pair may thus be summarized as follows. The vortex and the antivortex rotate
around each other, while the pair  shrinks due to dissipation.
The energy of the pair follows approximately the curve
of the right panel of Fig.~\ref{fig:rotate_energy} as its size (and its angular momentum)
decreases. At vanishing size a singular point of the type discussed
in the preceding section would be created and the total energy
would reach the finite value $E=4\pi$ (in rationalized units).
However, the discreteness of the lattice actually interrupts the process
when the size of the pair becomes comparable to the lattice spacing.
The VA pair disappears (i.e., the skyrmion number changes
abruptly from ${\cal N}=\pm 1$ to ${\cal N}=0$) and a burst
of energy equal to $4\pi$ is released into the system,
probably in the form of spin waves.
In physical units, the amount of energy released is given by $E=8\pi t A$
where $t$ is the film thickness and $A$ the exchange constant.
For typical values $t=10$ nm and $A=10^{-11}$ J/m we obtain the estimate
$E \sim 2.5\times 10^{-18}$ J which is apparently in rough
agreement with numerical simulations \cite{hertel06,tretiakov07}.

The scenario described above explains how a topologically forbidden
$(\Delta{\cal N} = 1)$ transition can take place in a real ferromagnet
but does not by itself account for the experimentally observed
vortex core switching. The complete scenario involves two distinct steps.
First, application of a short burst of an alternating magnetic field
creates a VA pair in the vicinity of a preexisting single vortex.
Second, a three-body collision takes place during which a suitable
pair of vortices is annihilated through a $\Delta{\cal N} = 1$
transition of the type described above, and the final product is a
single vortex with polarity opposite to that of the original vortex
\cite{waeyenberge06}.
We also note that a system with three vortices
(two vortices and an antivortex which form a cross-tie wall)
was observed in a rectangular platelet in \cite{kuepper07}.
Their spectrum of eigenmodes was studied experimantally.

   Here we do not address the question of how a VA pair is actually
created. Rather we concentrate on step 2 of the process and explain in some
detail how vortex core switching may occur in a three-body collision.
Specifically, let us assume that a single $(1,-1)$=C vortex
is initially at
rest at some specified point which is taken to be the origin of the
coordinate system. Let us further assume that a topologically trivial $({\cal N} =0)$
VA pair consisting of a $(1,1)$=A vortex and a $(-1,1)$=B antivortex is somehow
created in the neighborhood of the original vortex. Once created the AB pair
will undergo Kelvin motion of the type described in Section \ref{sec:kelvin} and eventually
collide with the single vortex C.

\begin{figure}
\epsfig{file=figure08.eps,width=0.75\linewidth}
   \caption{Three snapshots for the collision of a VA pair in Kelvin motion
(the AB pair), initially located at $(0,-15)$ and propagating
with velocity $v=0.1$, against a target vortex C initially located
at the origin.
During collision, antivortex B rotates around vortex C before rejoining
its original partner A to form a new VA pair that scatters off at an angle
in the third quadrant.
The target vortex C is shifted to a
new position in the fourth quadrant thanks to transmutation
of VA pair momentum to vortex position.
   }  \label{fig:switch_v01}
\end{figure}

   The process was simulated by a numerical solution of the
corresponding initial-value problem in the LL equation. Figure \ref{fig:switch_v01}
provides an illustration with three characteristic snapshots in the case of a relatively
slow Kelvin pair initially moving along the y axis with velocity $v=0.1$ for
which the vortex and the antivortex are separated by a distance
$d \approx 1/v=10$. As the pair approaches, the original C=$(1,-1)$ vortex teams
up with the B=$(-1,1)$ partner of the AB pair to form a new, topologically
nontrivial $({\cal N} =1)$ VA pair in quasi-rotational motion. In fact, B rotates
almost a full circle around C before rejoining its original partner A. The new
AB pair is again a topologically trivial $({\cal N} =0)$ VA pair in Kelvin motion that
moves away from the target vortex, having suffered a total scattering angle
that is greater than $\pi/2$ from its original direction.
The scattering is
inelastic in the sense that the outgoing AB pair moves out with greater
velocity $(v=0.15)$. And, most remarkably, the target vortex C moves away from
the origin and comes to rest at a new location in the fourth quadrant of the
xy plane.

   This rather unusual behavior is explained by the unusual nature of
the conservation laws (11) and (12) which allow for a transmutation between
position and impulse in the case of topologically nontrivial systems, such
as the three-vortex system considered here (a more detailed discussion
will be given in a future publication).

\begin{figure}
\epsfig{file=figure09.eps,width=0.75\linewidth}
   \caption{
Same as Fig.~\ref{fig:switch_v01} but for a larger initial velocity $v=0.5$ of
the AB pair. During collision, antivortex B begins
to rotate around vortex C but the rotating BC pair is
eventually annihilated leaving behind vortex A (with
polarity opposite to that of the target vortex C) and
a burst of spin waves that propagate away from the scattering region.
   }  \label{fig:switch_v05}
\end{figure}

   The preceding numerical experiment was repeated for a Kelvin pair
with relatively large velocity $v=0.5$ for which the vortex and the antivortex
are tightly bound at a relative distance $d=2.6$ \cite{papanicolaou99}. The process is
again illustrated by three characteristic snapshots in Figure \ref{fig:switch_v05}.
While the initial stages of the process are similar to those encountered in the case
of slow Kelvin motion (Figure \ref{fig:switch_v01}) a substantial departure occurs when the
pair now approaches the target vortex. In particular, as soon as antivortex
B=$(-1,1)$ begins to rotate around the target vortex C=$(1,-1)$ they collide and
undergo a spectacular $\Delta{\cal N} =1$ transition (annihilation) leaving behind the
A=$(1,1)$ vortex which may be thought of as the target vortex C=$(1,-1)$ with
polarity flipped from $-1$ to 1 (vortex core switching) and a burst of spin
waves propagating away from the scattering region.

  A detailed numerical investigation of the three-vortex process for
Kelvin waves with velocities in the allowed range $0<v<1$ suggests the
existence of the three characteristic regions separated by two critical
velocities $v_1=0.3$ and $v_2=0.9$ (such that $0<v_1<v_0<v_2<1$, with $v_0=0.78$ being the
critical velocity discussed in Section \ref{sec:kelvin}). For $0<v<v_1$, the Kelvin pair
undergoes nearly elastic scattering of the type depicted in Figure \ref{fig:switch_v01}.
For $v_1<v<v_2$, the process leads to a topologically forbidden $\Delta{\cal N} =1$
transition of the type illustrated in Figure \ref{fig:switch_v05}.
There is also some evidence that fast
Kelvin waves with velocities in the narrow range $v_2<v<1$ undergo a nearly
elastic scattering without inversion of the polarity of the target vortex.

\section{Conclusion}
\label{sec:conclusion}

The VA pairs analyzed in this paper are special examples of solitary waves
whose dynamics is closely related to their topological structure.
For example, the VA pairs studied in Section~\ref{sec:kelvin}
can undergo free translational motion because their topological charge
vanishes (${\cal N} = 0$).
In contrast, a VA pair with nonvanishing ${\cal N}$ performs
rotational motion around a fixed guiding center and is thus spontaneously
pinned within the ferromagnetic medium, as discussed
in Section~\ref{sec:rotate}.
Such a peculiar dynamical behavior would have been surprising had it not
occurred previously in the case of interacting electric charges in the
presence of a magnetic field.

It is then interesting to ascertain conditions under which
a certain field theory would exhibit a similar link between
topology and dynamics. A simple criterion was introduced in Ref.~\cite{komineas98}
and is briefly summarized as follows.
We restrict attention to 2D Hamiltonian systems described in terms of
$\Lambda$ pairs of canonically conjugate variables $(\Pi_i,\Phi_i)$
with $i=1,2,\ldots \Lambda$. Then one may define the local vorticity
\begin{equation}  \label{eq:vorticity2}
\gamma = \sum_{i=1}^\Lambda \epsilon_{\alpha\beta}\, \partial_\alpha\Pi_i\, \partial_\beta\Phi_i
\end{equation}
which is a simple generalization of the first step of Eq.~(\ref{eq:vorticity1}),
the remaining two steps being special to the specific example
considered in the present paper (ferromagnets).
Now, in general, the total vorticity $\Gamma=\int \gamma\, dx dy$ is expected
to vanish by a trivial partial integration using Eq.~(\ref{eq:vorticity2}).
On the other hand, a nonzero $\Gamma$ would signal a special topological structure
of the field theory under consideration and may lead to peculiar dynamics.
The theory analyzed in the present paper is an example (with $\Lambda=1$) which yields
a nonzero $\Gamma$ that may be identified with the Pontryagin index
${\cal N} = \Gamma/4\pi$. An example with $\Lambda=2$ is provided by
a 2D antiferromagnet where $\Gamma = 0$ except when an external field
is present which may lead to $\Gamma \neq 0$ and an interesting link
between topology and dynamics \cite{komineas98}.

Implicit in the preceding general argument is the fact that a topological charge
${\cal N}$ is conserved. Also taking into account that ${\cal N}$
is quantized, one would expect that a topological (${\cal N} \neq 0$) soliton
cannot be annihilated in a continuous manner.
Nevertheless, a quasi-continuous process was described in Ref.~\cite{komineas07}
and in the present paper according to which a rotating VA pair with ${\cal N} = 1$
may be reduced to a singular point and thereby be eliminated by lattice
discreteness without encountering an energy barrier.

A mechanism for changing the topological number of a magnetic
configuration makes it possible to obtain controlled switching
between topologically distinct (and thus robust) magnetic states.
This was achieved in the experiments of Refs.~\cite{waeyenberge06,yamada07}.
The dynamics underlying both experiments involves a three-vortex process
\cite{hertel07} initiated by the production of a topologically trivial
VA pair in the vicinity of a preexisting vortex by an
alternating magnetic field.
It has been predicted that the same phenomenon would occur
if one uses a rotating external field \cite{kravchuk07}.
The resulting three-vortex system carries nonzero topological charge
and is thus by itself a rotating object spontaneously pinned
in the magnet.
Also due to dissipation, a quasi-continuous process takes place
that changes the topological number by one unit,
leaving behind a burst of energy in the form of spin waves
and a single vortex with polarity opposite to that of the original vortex.

Although our strictly 2D treatment provides a detailed scenario for the
three-vortex process that leads to polarity switching, it does not
account for the initial production of a topologically trivial VA pair.
This is partly due to our approximation of a thin film with infinite extent.
Implicit in this approximation is the assumption that the demagnetizing
field amounts to a simple (additive) renormalization of easy-plane
anisotropy. While this assumption appears to be firmly established
for static magnetic states \cite{gioia97} we do not know of a
corresponding mathematical derivation for dynamical processes of the
type discussed here.

It is also important to visualize how the formation of a singular
point discussed in the present strictly 2D context appears within
a realistic magnetic element of finite extent.
Numerical simulations \cite{hertel07} show that, in an element
of finite thickness, a singular point is first created at one of the
surfaces of the element. The VA pair then vanishes by formation
and subsequent annihilation of a singular point at successive
levels away from the surface. At the stage when a singular point has
been formed and annihilated, say, near the top surface, while the
VA pair is still present in the bulk of the element,
a Bloch Point (BP) is created in the element.
This is a somewhat simplified realization of the BP studied
in Ref.~\cite{doering68}.
Needless to say, the BP created near the top surface is eventually
annihilated when the VA pair exits the system through the lower surface.
It is important to emphasize the unusual fact that during
creation and annihilation of the BP the system does not have to overcome
an energy barrier, unlike the case discussed in \cite{slonczewski75},
essentially for the same reasons explained for the strictly
2D VA pairs studied in the main text.

\begin{acknowledgments}
N.P. is grateful for hospitality at the Max-Planck Institute
for the Physics of Complex Systems (Dresden) where this work was completed.
\end{acknowledgments}

\begin{appendix}
\section{Electric charges in a magnetic field}

Most of the distinct features of the dynamics of VA pairs occur also
in the dynamics of electric charges in the presence of a uniform magnetic field $\bm{B}$.
Although we shall mainly be interested in 2D motion in a plane perpendicular to $\bm{B}$,
it is convenient to keep for the moment 3D notation and write the equations
of motion for two interacting charges $e_1$ and $e_2$:
\begin{equation}  \label{eq:A1}
m\, \frac{d\bm{v}_1}{d t} = \bm{F}_1 + e_1 (\bm{v}_1\times\bm{B})\,,\quad
m\, \frac{d\bm{v}_2}{d t} = \bm{F}_2 + e_2 (\bm{v}_2\times\bm{B})\,,
\end{equation}
where
\begin{equation}  \label{eq:A2}
\bm{F}_1 = -\bm{F}_2 = -\frac{\bm{r}_1-\bm{r}_2}{|\bm{r}_1-\bm{r}_2|}\;
                                          V'(|\bm{r}_1-\bm{r}_2|)
\end{equation}
is the mutual force derived from a potential energy $V=V(|\bm{r}_1-\bm{r}_2|)$
and $V'$ denotes derivative with respect to the argument $|\bm{r}_1-\bm{r}_2|$.
The conserved energy functional is then given by
\begin{equation}  \label{eq:A3}
E = \frac{1}{2} m (v_1^2 + v_2^2) + V(|\bm{r}_1-\bm{r}_2|),
\end{equation}
which does not depend explicitly on the magnetic field, while the conserved
linear momentum (impulse) is now given by
\begin{equation}  \label{eq:A4}
\bm{P} = m (\bm{v}_1+\bm{v}_2) - (e_1\bm{r}_1 + e_2 \bm{r}_2)\times \bm{B}
\end{equation}
and differs from the usual mechanical definition by an important field dependent term.
This (second) term actually indicates a rather profound influence of the magnetic field
on the dynamics of electric charges.
For instance, note that a shift of the origin of coordinates by a constant vector
$\bm{c}$, thus $\bm{r}_1 \to \bm{r}_1 - \bm{c}$
and $\bm{r}_2 \to \bm{r}_2 - \bm{c}$, induces a nontrivial change on the impulse
$\bm{P}$ of Eq.~(\ref{eq:A4}) given by
\begin{equation}  \label{eq:A5}
\bm{P} \to \bm{P} + (e_1+e_2) (\bm{c}\times \bm{B}).
\end{equation}
This unusual behavior is analogous to that of the impulse defined
by Eq.~(\ref{eq:linmomentum}) in the case of field configurations
with nonvanishing total topological vorticity $\Gamma$ or skyrmion number ${\cal N}$.
Here we may abstract from Eq.~(\ref{eq:A5}) the analog of the topological vorticity
in the present problem:
\begin{equation}  \label{eq:A6}
\Gamma \sim (e_1+e_2)\,B.
\end{equation}
An electron-positron pair ($\Gamma=0$) may undergo Kelvin motion,
while two like charges ($\Gamma \neq 0$) perform rotational motion around a fixed guiding
center, as determined by the explicit solutions
constructed and briefly analyzed in the remainder of this appendix.

Consider first the case of an electron-positron pair ($e_1=-e_2=e$) for which
a special 2D solution of Eqs.~(\ref{eq:A1}) is given by
\begin{eqnarray}  \label{eq:A6a}
 x_1 = \frac{d}{2}\,,\quad y_1 & = & v t\,,\quad z_1=0\,,  \nonumber \\
 x_2 = -\frac{d}{2}\,,\quad y_2 & = & v t\,,\quad z_2=0\,,
\end{eqnarray}
where the electron and the positron are located at a constant
relative distance $d$ along the x axis and move in formation
along the y axis with constant velocity
\begin{equation}  \label{eq:A7}
 v = \frac{V'(d)}{e\,B}\,,
\end{equation}
in close analogy with the Kelvin motion of the VA pair discussed
in Section~\ref{sec:kelvin}.
We further calculate energy $E$ from Eq.~(\ref{eq:A3}) and impulse
$\bm{P}=(0,P,0)$ from Eq.~(\ref{eq:A4}) to find
\begin{eqnarray}  \label{eq:A8}
 E & = & m v^2 + V(d) = m \left[\frac{V'(d)}{eB}\right]^2 + V(d) \nonumber \\
 P & = & 2 m v + eBd  = 2m \frac{V'(d)}{eB} + eBd.
\end{eqnarray}
Hence all relevant quantities are given in parametric form as functions
of relative distance $d$ once the potential energy $V=V(d)$ has been specified.
But there are some generic properties of this solution that are practically
independent of the choice of $V(d)$.
We may take derivatives with respect to $d$ of both sides
of Eq.~(\ref{eq:A8}) to find $E' = V'(1+2 m V''/e^2B^2)$ and
$P' = eB(1+2 m V''/e^2B^2)$.
An immediate consequence of these relations is the group-velocity
relation $v=dE/dP$. We also note that both $E$ and $P$ may acquire an extremum
at a common value of $d$ (or $v$) determined from
\begin{equation}  \label{eq:A9}
 1 + \frac{2m V''(d)}{e^2B^2} = 0.
\end{equation}
This is actually the reason for the appearance of a cusp in the $E$ vs $P$ dispersion
analogous to that encountered in the Kelvin motion of VA pairs, which now appears
to be a generic feature of a wide class of physical systems.

\begin{figure}
\epsfig{file=figure10.eps,width=0.80\linewidth}
   \caption{Energy $E$ and impulse $P$ as functions of velocity $v$
(left panel) and $E$ vs $P$ dispersion (right panel)
for an electron-positron pair ($e_1=-e_2$) in Kelvin motion.
   }  \label{fig:pn_energy}
\end{figure}

For an explicit demonstration we make the special choice of potential energy
\begin{equation}  \label{eq:A10}
 V = 2\pi\, \ln|\bm{r}_1-\bm{r}_2|
\end{equation}
in order to model the behavior of VA pairs at large relative distance
\cite{pokrovskii85,kovalev02}.
For a graphical illustration we also make the special choice
of constants $m=1$ and $eB=2\pi$, as suggested by Eq.~(\ref{eq:A6}),
to write
\begin{equation}  \label{eq:A11}
 v = \frac{1}{d}\,,\quad E = \frac{1}{d^2} + 2\pi\, \ln d\,,\quad
P = \frac{2}{d} + 2\pi\,d
\end{equation}
where we note that both $E$ and $P$ acquire a minimum at a distance
$d=d_0=1/\sqrt{\pi}$ or velocity $v=v_0=\sqrt{\pi}$.
The dependence of $E$ and $P$ on the velocity $v$ as well as
the $E$ vs $P$ dispersion are shown in Fig.~\ref{fig:pn_energy}
and are found to be closely analogous to the results of
Fig.~\ref{fig:kelvin_energy}
pertaining to Kelvin motion of VA pairs.
In particular, for a widely separated pair (branch I)
we find from Eq.~(\ref{eq:A11}) that $P\sim 2\pi\,d,\; vP \sim 2\pi,\; v=1/d$,
and $E=2\pi\,\ln (P/P_0)$ with $P_0=2\pi$,
in close analogy with the asymptotic results of
Eqs.~(\ref{eq:estimates1})--(\ref{eq:kelvinpairenergy}).
The appearance of a cusp and consequently of branch II in the
spectrum is also notable.
But the details of branch II are different than those
of Fig.~\ref{fig:kelvin_energy} and Eq.~(\ref{eq:kelvinpairenergy2}).
There is now no upper limit in the velocity $v$.
In fact, all $v, E$ and $P$ in Eq.~(\ref{eq:A11}) diverge in the limit of small $d$ and
$E\sim P^2/4$, which coincides with the dispersion $P^2/2M$
of a free particle with
mass equal to the total mass of the pair $(M=2 m = 2)$.

We have also carried out a stability analysis of the special Kelvin-like
solution (\ref{eq:A6a}) to find that the motion is marginally stable along branch I
but becomes unstable along branch II.
This conclusion is in agreement with a similar result obtained in Ref.~\cite{jones86}
in the case of a vortex ring in a superfluid, as is further discussed in the concluding
remarks of our Section~\ref{sec:kelvin}.

As a last example we consider the case of 2D motion of two like charges $e_1=e_2=e$.
Then a special solution of Eq.~(\ref{eq:A1}) is given by
\begin{equation}  \label{eq:A12}
x_1 = - x_2 = R\,\cos\omega t\,,\quad y_1 = -y_2 = R\,\sin\omega t\,,\quad z_1=z_2=0
\end{equation}
which describes a pair rotating at constant radius $R=d/2$ and angular
frequency $\omega=v/R$ where the velocity $v$ is determined from the algebraic equation
\begin{equation}  \label{eq:A13}
 \frac{m v^2}{R} + e B v - V'(d) = 0
\end{equation}
that expresses the exact balance of the centrifugal, the magnetic, and the mutual force.
The conserved energy is still calculated from Eq.~(\ref{eq:A3}) with
$v_1^2=v_2^2=v^2$, but the conservation of the impulse $\bm{P}$ of Eq.~(\ref{eq:A4})
is simply equivalent to the statement that rotation takes place around a
fixed guiding center. More relevant is now the conserved angular momentum
(impulse) which is given by
\begin{equation}  \label{eq:A14}
 L = m (x_1 \dot{y}_1 - y_1 \dot{x}_1) + \frac{e_1 B}{2} (x_1^2 + y_1^2)
   + m (x_2 \dot{y}_2 - y_2 \dot{x}_2) + \frac{e_2 B}{2} (x_2^2 + y_2^2),
\end{equation}
where the overdot denotes time derivative. Again, the angular momentum
differs from its standard mechanical expression by important
field-dependent terms.

\begin{figure}
\epsfig{file=figure11.eps,width=0.80\linewidth}
   \caption{Energy $E$ and angular impulse $L$ as functions of
angular velocity $\omega$
(left panel) and $E$ vs $L$ dispersion (right panel)
for a pair of like charges ($e_1=e_2$) rotating around
a fixed guiding center.
   }  \label{fig:pp_energy}
\end{figure}

Now, for the specific choice of potential energy given by Eq.~(\ref{eq:A10}),
and constants $m=1, e_1=e_2=e$ and $eB=2\pi$,
the algebraic equation (\ref{eq:A13}) yields
\begin{equation}  \label{eq:A15}
 v = \pi \left(\sqrt{R^2+\frac{1}{\pi}} - R\right)\,,
\quad R \equiv \frac{d}{2}\,,
\end{equation}
while the angular velocity $\omega$, the energy $E$, and the angular momentum $L$,
read
\begin{equation}  \label{eq:A16}
 \omega = \frac{v}{R}\,,\quad E = v^2 + 2\pi\,\ln(2R)\,,\quad
 L = 2 R v + 2\pi\,R^2.
\end{equation}
In view of Eq.~(\ref{eq:A15}) all three quantities in (\ref{eq:A16})
are expressed
in terms of a single parameter $R=d/2$. As a check of consistency one may
verify the relation $\omega=dE/dL$ using Eqs.~(\ref{eq:A15}-\ref{eq:A16}).

The dependence of $E$ and $L$ on angular velocity $\omega$ as well as the
$E$ vs $L$ dispersion are shown in Fig.~\ref{fig:pp_energy}.
Again there exists a close analogy with the results of Section~\ref{sec:rotate}
on rotating VA pairs.
In particular, for large diameter $d$, Eqs.~(\ref{eq:A16}) yield
$L\sim \frac{\pi}{2} d^2,\; \omega L \sim \pi,\; \omega\sim 2/d^2$,
and $E = \pi \ln(L/L_0)$ with $L_0=\pi/2$,
which should be compared with the asymptotic results for VA pairs
given in Eqs.~(\ref{eq:estimates2})--(\ref{eq:rotatepairenergy}).
On the other hand, some quantitative differences arise at small $d$, where the
angular momentum vanishes as expected ($L\sim \sqrt{\pi} d$) but the
angular frequency diverges ($\omega \sim 2\sqrt{\pi}/d$).
Furthermore, the energy $E$ does not reach a finite value at $L=0$
(as was the case for rotating VA pairs) but diverges logarithmically
to minus infinity.

Finally, an analysis of mechanical stability \cite{dialynas97} shows that
circular motion of two like charges in a magnetic field is marginally stable
 for all values of $d$, in contrast to the Kelvin motion discussed earlier
in this section which becomes unstable at small $d$.
However, unlike Kelvin motion which proceeds with no acceleration,
a rotating pair is expected to radiate when full electrodynamics is turned on.
Surely, this is also a source of instability and may indicate a similar instability
for the rotating VA pairs discussed in Section~\ref{sec:rotate}.

\end{appendix}


text:
\section{Introduction}
In 2006 in the {\it On-Line Encyclopedia of Integer Sequences} (OEIS)
\cite{OEIS},
sequence \texttt{A119326}, Paul Barry
introduced a modified Pascal triangle, defined for integers
$0 \leq k \leq n$, as follows:
\[T(n,k)=\sum_{\substack{0\leq j\leq n-k\\2\mid j}}\binom{k}{j}\binom{n-k}{j}.\]
The first few rows of this triangle are as follows:
\begin{gather*}
1 \\
1 \qquad 1 \\
1 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 1 \qquad 1  \\
1 \qquad 1 \qquad 2 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 4 \qquad 4 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 7 \qquad 10 \qquad 7 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 11 \qquad 19 \qquad 19 \qquad 11 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 16 \qquad 31 \qquad 38 \qquad 31 \qquad 16 \qquad 1 \qquad 1
\end{gather*}
\newpage
Similarly, one can consider $T(n,k) \bmod 2$, whose terms are given by
sequence \texttt{A114213}:
\begin{gather*}
1 \\
1 \qquad 1 \\
1 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 1 \qquad 1  \\
1 \qquad 1 \qquad 0 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 0 \qquad 0 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 1 \qquad 0 \qquad 1 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 1 \qquad 1 \qquad 1 \qquad 1 \qquad 1 \qquad 1 \\
1 \qquad 1 \qquad 0 \qquad 1 \qquad 0 \qquad 1 \qquad 0 \qquad 1 \qquad 1
\end{gather*}

Sequences \texttt{A114212} and \texttt{A114214},
respectively, are the row sums and
diagonal sums of this latter triangle.  We denote them by
$r(n)$ and $d(n)$, respectively:
\begin{align*}
r(n) &= \sum_{k=0}^n \left( T(n,k) \bmod 2 \right) \\
d(n)&=\sum_{k=0}^{\lfloor n/2 \rfloor} \left( T(n-k,k) \bmod 2 \right).
\end{align*}

In May 2016, the first author observed, empirically, a connection between
$d(n)$ and the binary representation of $n$.  In this note we prove this
connection, and also prove a formula for $r(n)$.  The connection
involves Stern's ``diatomic sequence'' $s(n)$, defined by
$s(0) = 0$, $s(1) = 1$, $s(2n) = s(n)$, and
$s(2n+1) = s(n) + s(n+1)$; see \cite{Stern:1858}.
\section{The diagonal sums}

Let the binary representation of $n$ be denoted by
$\sum_{i=0}^j \varepsilon_i (n) 2^i$.
We show that the diagonal sum $d(n)$ can be expressed in terms of this
representation.
Given a string $s$ of $0$'s and $1$'s,
we consider its {\it run lengths}:
the lengths of maximal blocks of consecutive identical elements.  For example,
if $s = 111000011111$, then the run lengths of $s$ are $(3,4,5)$.

If $m$ is a sequence of positive integers, we may associate an
integer with it via the continued fraction expansion: if $m=(m_0,\ldots,m_k)$, we say that the \emph{continuant} of $m$ is the numerator of the
continued fraction $[m_0;m_1,\ldots,m_k]$ (see \cite[Ch.~34, \S 4]{Chrystal}).
\begin{theorem}\label{mainthm}
Let $n\geq 0$ be an integer and let $m$ be the sequence of
run lengths of the binary representation of $n$.
Then $d(n)$ equals the continuant of $m$.
\end{theorem}

We will use Lucas' famous congruence for binomial
coefficients \cite[p.~230]{Lucas:1878}:
if $p$ is a prime number and $n=(n_\nu\cdots n_0)_p$ and $k=(k_\nu\cdots k_0)_p$, then
\[\binom nk\equiv \binom{n_\nu}{k_\nu}\cdots \binom{n_0}{k_0}\pmod p . \]
This implies that $\binom nk$ is not divisible by $p$ if and only if $k_i\leq n_i$ for all $i$.
Moreover, it follows that the number of odd binomial coefficients $\binom nk$ equals $2^{s_2(n)}$, where $s_2$ is the binary sum-of-digits function
\cite{Glaisher:1899}.

We prove the following statement, which reduces the problem to divisibility by $2$ of binomial coefficients.
We will derive Theorem~\ref{mainthm} from it in a moment.
\begin{proposition}\label{prp_reduction}
Let $n$ and $k$ be nonnegative integers such that $k\leq n$.
If $2\mid n+k$, then $T(n,k)\equiv \modd{\binom nk} {2}$.
Otherwise, $T(n,k)\equiv \modd{\binom{n-1}k} {2}$.
\end{proposition}
\begin{proof}
We prove the first statement. By replacing $n$ with $n+k$ we get
the equivalent assertion that if $2\mid n$ or $2\mid k$, then
\begin{equation}\label{eqn_sufficient}
\sum_{\substack{0 \leq j \leq n\\2\mid j}}
\binom nj\binom kj\equiv \binom {n+k}k\pmod 2 .
\end{equation}
By Lucas' congruence the left-hand side is congruent to
\[\sum_{j=0}^{n}\binom nj\binom kj\equiv\sum_{j=0}^{n}\binom {n\wedge k}j
\equiv 2^{s_2(n\wedge k)} \pmod 2,
\]
where $n\wedge k$ is the integer whose binary digits satisfy $\varepsilon_i(n\wedge k)=\min(\varepsilon_i(n),\varepsilon_i(k))$.
This expression is odd if and only if $s_2(n\wedge k)=0$, which is the case if and only if the binary representations of $n$ and $k$ are disjoint.
To handle the right-hand side of Eq.~\eqref{eqn_sufficient}, we note that $\binom{n+k}{k}$ is odd if $n\wedge k=1$.
On the other hand,
if the binary representations of $n$ and $k$ are not disjoint,
then the condition $\varepsilon_i(k)\leq \varepsilon_i(n+k)$ is
violated for
$i=\min\{j:\varepsilon_j(n)=1,
\varepsilon_j(k)=1\}$; therefore $\binom{n+k}{k}$ is even.
This proves the first assertion.

For the second assertion, we use Lucas' congruence again: for $2\mid j$ and $2\mid m$ we have $\binom mj\equiv \binom{m+1}j\pmod 2$. Since $2\nmid n-k$, we obtain $\binom{n-k}j\equiv\binom{n-1-k}j\pmod 2$.
Moreover, by $2\nmid n-k$ the ranges of summation in $T(n,k)$ and $T(n-1,k)$ are the same.
\end{proof}
From this proposition we obtain in particular the identity
\begin{equation}\label{eqn_double}
d(2n)=d(2n+1).
\end{equation}

Carlitz \cite{Carlitz:1962} proved
that Stern's diatomic sequence $s(n)$ satisfies
$s(n+1)=\sum_{k=0}^{\lfloor n/2\rfloor} \left( \binom{n-k}k \bmod 2 \right)$.
By Proposition~\ref{prp_reduction} and Eq.~\eqref{eqn_double}
we therefore have
\begin{equation}\label{eqn_stern_relation}
d(2n)=d(2n+1)=s(2n+1).
\end{equation}

It is well-known~\cite{L1929, L1969} that
if $m=(m_0,\ldots,m_k)$ is the sequence of run-lengths of the binary
representation of $n$ and $n$ is odd,
then $s(n)$ is the continuant of $m$.
Therefore $d(n)$ is the continuant of $m$.
In order to complete the proof of the conjecture,
we have to show that the same is true for even $n$.
By Eq.~\eqref{eqn_stern_relation} it is sufficient to prove the following lemma.

\begin{lemma}
If $n$ is even, then the continuant of the sequence of run-lengths of the binary representation of $n$ is equal to the continuant corresponding to $n+1$.
\end{lemma}
\begin{proof}
Let $n=\tL^{m_0}\tO^{m_1}\cdots \tL^{m_{k-1}}\tO^{m_k}$.
We distinguish between two cases. If $m_k=1$, then
$n+1=(\tL^{m_0}\tO^{m_1}\cdots \tO^{m_{k-2}}\tL^{m_{k-1}+1})$
and the statement follows from the identity $[m_0;m_1,\ldots,m_{k-1},1]=[m_0;m_1,\ldots,m_{k-1}+1]$. If $m_k\geq 2$, then
$n+1=(\tL^{m_0}\tO^{m_1}\cdots \tO^{m_{k-2}}\tL^{m_{k-1}}\tO^{m_k-1}1)$ and the statement follows from
$[m_0;m_1,\ldots,m_k]=[m_0;m_1,\ldots,m_{k-1},m_k-1,1]$.
\end{proof}

\begin{remark}
The sequence $(d(n))_{n\geq 0}$ is a $2$-regular sequence
\cite{Allouche&Shallit:1992}, as it satisfies the equalities
\begin{align*}
d(2n+1) &= d(2n) \\
d(4n+2) &= 3d(2n) - d(4n) \\
d(8n) &= -d(2n) + 2 d(4n) \\
d(8n+4) &= 4d(2n) - d(4n) .
\end{align*}
\end{remark}
\section{The row sums}

We will prove

\begin{theorem}
$$r(n)= \begin{cases}
	2^{s_2 (n)}, & \text{if $n$ odd}; \\
	2^{s_2(n)} +2^{s_2(n-2)}, & \text{if $n$ even}.
\end{cases}
$$
\end{theorem}

A similar characterization was stated, without proof or attribution,
in the notes to \texttt{A114212} of the OEIS.

\begin{proof}
From Proposition~\ref{prp_reduction} we get,
for integers $n \geq k \geq 0$, that
\begin{align*}
T(2n,2k) &\equiv T(2n+1,2k) \equiv T(2n+1,2k+1) \equiv \modd{{n \choose k}} {2} ; \\
T(2n,2k+1)&\equiv \modd{ { {n-1} \choose k}} {2} .
\end{align*}
Then
\begin{align*}
r(2m) &= \sum_{k=0}^{2m} \left( T(2m,k) \bmod 2 \right) \\
&=  \sum_{k=0}^m \left( T(2m,2k) \bmod 2 \right)
\quad + \quad
\sum_{k=0}^{m-1} \left( T(2m, 2k+1) \bmod 2 \right)  \\
&=  \sum_{k=0}^m \left( {m \choose k} \bmod 2 \right)
\quad + \quad
\sum_{k=0}^{m-1} \left( {{m-1} \choose k} \bmod 2 \right)  \\
&= 2^{s_2 (m)} + 2^{s_2 (m-1)}  \\
&= 2^{s_2 (2m)} + 2^{s_2 (2m-2)}.
\end{align*}
Similarly,
\begin{align*}
r(2m+1) &= \sum_{k=0}^{2m+1} \left( T(2m+1,k) \bmod 2 \right) \\
&=  \sum_{k=0}^m \left( T(2m+1,2k) \bmod 2 \right)
\quad + \quad
 \sum_{k=0}^m \left( T(2m+1, 2k+1) \bmod 2 \right)  \\
&=  \sum_{k=0}^m \left( {m \choose k} \bmod 2 \right)
\quad + \quad
 \sum_{k=0}^m \left( {m \choose k} \bmod 2 \right)   \\
&= 2^{s_2 (m)} + 2^{s_2 (m)} \\
&= 2^{s_2 (2m+1)} .
\end{align*}
\end{proof}
\bibliographystyle{siam}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
A long-standing goal of the solar interior physics is to reproduce active regions, composed mainly of the sunspots, self-consistently from
magnetic fluxes generated in the solar interior. We now approach the subject from two different theoretical perspectives: one focuses on
emergence and organization processes of the magnetic flux in the uppermost part of the convection zone (CZ), and the other explores
the flux generation and maintenance processes, i.e., the dynamo process, operating deeper down.

Several leading-edge numerical studies, which focus on the uppermost part of the solar CZ, have succeeded to simulate spontaneous
formations of concentrated magnetic structures reminiscent of active regions  \citep[e.g.,][]{cheung+10,stein+12,rempel+14,kapyla+15}.
In these studies, the solar surface convection and its nonlinear interaction with the magnetic field were simulated in a more or less realistic
manner with the steep density gradient just below the photosphere and/or the radiative transfer with the ionization in Cartesian domains.
However, since some sort of the large-scale seed magnetic field has been assumed inconsistently as initial or boundary condition, the dynamo
mechanism and its connection to the formation process of the active region were beyond the scope of these studies.

A growing body of evidence is accumulating to demonstrate that solar-like cyclic large-scale magnetic field is organized in global
spherical-shell convections \citep[e.g.,][]{ghizaru+10,kapyla+12,masada+13,yadav+15,augustson+15}. Despite some differences in the
numerical setup and method, there is a common outcome of the convective dynamo in these studies: diffuse magnetic flux extending over
the CZ and/or the tachocline instead of magnetic flux tubes expected in the standard solar dynamo paradigm \citep[e.g.,][and references
therein]{charbonneau10}. Although the flux emergence like event from distributed magnetic flux has been occasionally observed in some
models \citep[][]{nelson+13,fan+14}, its universality or feasibility in the Sun is still a matter of considerable debate.

There is still a large gap between the dynamo in the interior and the active region formation at the surface. Our study in this Letter would
be a first step aiming to bridge the gap between them. By advancing our previous works of weakly-stratified MHD convection
\citep[][hereafter MS14a,b]{masada+14a,masada+14b}, we perform convective dynamo simulation in a strongly-stratified atmosphere
resembling the solar interior in Cartesian geometry. The spontaneous formation of large-scale magnetic structures in the CZ surface
self-consistently from the large-scale convective dynamo is reported.
\section{Numerical setup}
Convective dynamo system is solved numerically in Cartesian domain. Our model covers only the CZ of depth $d_{\rm cz}$
($0 \le z \le d_{\rm cz}$) with omitting a stably-stratified layer below it, where $x$- and $y$-axes are taken to be horizontal and
$z$-axis is pointing downward. We set the width of the domain to be $W = 4d_{\rm cz}$.

We solve the fully-compressible MHD equations in the rotating frame of reference with a constant angular velocity of
${\bm \Omega} = -\Omega_0 {\bm e}_z$,
\begin{eqnarray}
\frac{\partial \rho}{\partial t} & = & - \nabla\cdot (\rho {\bm u})  \;, \label{eq1} \\
\frac{{\mathcal{D} }{\bm u}}{\mathcal{D} t}  & = & - \frac{\nabla P}{\rho}
+ \frac{{\bm J}\times{\bm B} }{\rho} - 2{\bm \Omega}\times{\bm u} + \frac{\nabla \cdot {\bm \Pi}}{\rho} + {\bm g}  \;, \ \ \ \ \label{eq2} \\
\frac{\mathcal{D}\epsilon }{\mathcal{D} t} & = & - \frac{P\nabla\cdot {\bm u}}{\rho}  + \mathcal{Q}_{\rm heat} \;, \label{eq3} \\
\frac{\partial {\bm B} }{\partial t} & = & \nabla \times ({\bm u} \times {\bm B} - \eta_0 {\bm J})\;, \label{eq4}
\end{eqnarray}
with the viscous stress ${\bm \Pi}$ and the heating term $\mathcal{Q}_{\rm heat}$ of
\begin{eqnarray}
 \Pi_{ij} &=& 2\rho\nu_0S_{ij} = \rho\nu_0\left( \frac{\partial u_i}{\partial x_j} +  \frac{\partial u_j}{\partial x_i} - \frac{2}{3}\delta_{ij}\frac{\partial u_i}{\partial x_i} \right) \;,  \label{eq5} \\
\mathcal{Q}_{\rm heat} &=&   \frac{\gamma\nabla\cdot(\kappa_0 \nabla\epsilon)}{\rho} + 2\nu_0 \bm{S}^2 + \frac{\mu_0\eta_0\bm{J}^2}{\rho}\;, \label{eq6}
\end{eqnarray}
where $\mathcal{D}/\mathcal{D}t$ is the total derivative, $\epsilon = c_{\rm V} T$ is the specific internal energy,
${\bm J} = \nabla \times {\bm B}/\mu_0$ is the current density with the vacuum permeability $\mu_0$, and ${\bm g} = g_0 {\bm e}_z$ is the
gravity of the constant $g_0$. The viscosity, magnetic diffusivity, and thermal conductivity are represented by $\nu_0$, $\eta_0$, and
$\kappa_0$, respectively. A perfect gas law $P = (\gamma -1)\rho \epsilon$ with $\gamma = 5/3$ is assumed.

The initial hydrostatic state is polytropic stratification given by
\begin{equation}
\epsilon = \epsilon_0 + \frac{g_0z}{(\gamma-1)(m + 1)}\;, \ \ \ {\rm and} \ \ \ \rho = \rho_0 (\epsilon/\epsilon_0)^m \;,
\end{equation}
with the initial surface internal energy $\epsilon_0$, surface density $\rho_0$, and the polytropic index $m = 1.49$,
providing the super-adiabaticity of $\delta\ \equiv \nabla - \nabla_{\rm ad}\ = 1.6 \times 10^{-3}$, where
$\nabla_{\rm ad} = 1-1/\gamma$ and $\nabla = (\partial \ln T/\partial \ln P)$.

Normalization quantities are defined by setting $d_{\rm cz}/2 = g_0= \rho_0 = \mu_0 = c_p = 1$.
The normalized pressure scale-height at the surface, defined by $\xi = H_p/d_{\rm cz} = (\gamma-1)\epsilon_0/(g_0d_{\rm cz})$, controls
the stratification level and is chosen here as $\xi = 0.01$, yielding a strong stratification with the density contrast between top and
bottom CZs about $700$.

Figure 1a shows the initial profiles of the density (solid) and temperature (dashed) of our model. The
density profile in the range $0.71 \le r/R_\odot \le 0.991$ of the standard solar model is also shown in
such a way as to fit the computational domain (dash-dotted) \citep[Model S:][]{CD+96}.
Our model has a stratification almost encompassing the solar CZ except its uppermost part.

All the variables are assumed to be horizontally periodic. Stress-free boundary conditions are used in the vertical
direction for the velocity. Perfect conductor and vertical field conditions are used for the magnetic field at the bottom and top
boundaries. A constant energy flux which drives the thermal convection is imposed on the bottom boundary, while the specific
internal energy is fixed at the top boundary.

The fundamental equations are solved by the second-order Godunov-type finite-difference scheme that employs an approximate MHD
Riemann solver \citep{sano+98}. The magnetic field evolves with the Consistent MoC-CT method \citep{evans+88,clarke96}.
Non-dimensional parameters ${\rm Pr} =20$, ${\rm Pm} =2$, ${\rm Ra} = 6\times10^7$, angular velocity of $\Omega_0 = 0.5$ and
the spatial resolution of ($N_x, N_y, N_z$) $=$ $(256,256,256)$  are adopted, where the Prandtl, magnetic Prandtl, and Rayleigh
numbers are defined by
\begin{equation}
{\rm Pr}  =  \frac{\nu_0}{(\kappa_0/\rho c_p )} \;,\ {\rm Pm} = \frac{\nu_0}{\eta_0}\;, \ {\rm Ra}  =  \frac{g_0 d_{\rm cz}^4}{\chi_0 \nu_0}\frac{\delta}{H_{p}}  \;,
\end{equation}
where $\rho$, $\delta$, and $H_p$ are evaluated at $z = d_{\rm cz}/2$.

In the following, the volume-, horizontal-, $x$- and $y$- averages are denoted by single angular brackets with subscripts ``\rm v", ``\rm h", ``\rm x"
and ``\rm y", respectively. The time-average of each spatial mean is denoted by additional angular brackets. The
relative importance of convection to the rotation is measured by the Rossby number ${\rm Ro} = (u_{\rm cv}k_f)/(2\Omega_0)$, where
$u_{\rm cv} \equiv \sqrt{\langle \langle u_z^2 \rangle_{\rm v} \rangle}$ is the mean convective velocity, and $k_f = 2\pi/d_{\rm cz}$.
The global convective turn-over time and the equipartition field strength are defined by
$\tau_{\rm cv} \equiv 1/(u_{\rm cv}k_f)$ and $B_{\rm eq} (z) \equiv \sqrt{\langle \mu_0 \rho {\bm u}^2 \rangle_{\rm h}} $.
Note that $B_{\rm eq}(z)$ is evaluated from the local convective energy and thus has a depth-dependence.

Since the sound speed in the deep CZ becomes very large in the strongly-stratified model and imposes a strict limit on the time-step,
a long thermal relaxation time is required in our fully-compressible simulation. To alleviate it, we first construct a progenitor model, in
which the convection reaches a fully-developed state and the system becomes thermally-relaxed, by evolving a non-rotating
hydrodynamic run for $800\tau_{\rm cv}$.

Shown in Figure 1b and 1c are the vertical profile of $\langle\langle u_z^2 \rangle\rangle_{\rm h}^{1/2}$ (dashed) and the distributions
of $u_z$ in the $x$-$y$ plane at $z/d_{\rm cz} = 0.04,\ 0.23,\ 0.78$ in the equilibrated state of the progenitor model. The black (orange)
tone denotes down-(up-)flows. The multi-scale convection with the strong up-down asymmetry, i.e., the slower and broader upflow cell
surrounded by networks of the faster and narrower downflow lanes, is developed in the progenitor model \citep[e.g.,][]{spruit+90,miesch05}. The
dynamo run is started by adding the rotation and a seed weak horizontal field to the progenitor model.

\begin{figure}[tbp]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=8.5cm,clip]{f1a.pdf} \\
\includegraphics[width=8.5cm,clip]{f1b.pdf} \\
\includegraphics[width=8.5cm,clip]{f1c.pdf}
\end{tabular}
\caption{(a) Vertical profiles of the initial density (solid) and temperature (dashed) of the simulation model, and the density profile of the
CZ of the standard solar model (dash-dotted). Normalization units are their surface values. (b) Vertical profiles of $\langle\langle u_z^2
\rangle\rangle_{\rm h}^{1/2}$ for the progenitor (dashed) and dynamo (solid) runs. (c) The horizontal distributions of the $u_z$ at
$z/d_{\rm cz} = 0.04,\ 0.23,\ 0.78$ of the progenitor (upper) and dynamo (lower) runs. }
\label{fig1}
\end{center}
\end{figure}
\begin{figure*}[htbp]
\begin{minipage}{0.35\hsize}
\begin{center}
\includegraphics[width=5.5cm,clip]{f2a.pdf}
\end{center}
\end{minipage}
\begin{minipage}{0.35\hsize}
\begin{center}
\includegraphics[width=10.5cm,clip]{f2b.pdf}
\end{center}
\end{minipage}
\caption{(a) The temporal evolution of $\epsilon_{\rm K}$ (solid), $\epsilon_{\rm Mt}$ (dashed) and $\epsilon_{\rm Mm}$ (dash-dotted).
(b) Time-depth diagrams of $\langle B_x \rangle_{\rm h}/B_{\rm eq}(z)$ and $\langle B_y \rangle_{\rm h}/B_{\rm eq}(z)$. }
\label{fig2}
\end{figure*}
\section{Simulation Results}
\subsection{Basic Properties of Convection and Dynamo}
The temporal evolutions of $\epsilon_{\rm K} \equiv \langle \rho\bm{u}^2/2 \rangle_{\rm v}$,
$\epsilon_{\rm Mt} \equiv \langle \bm{B}^2/2\mu_0 \rangle_{\rm v}$ and
$\epsilon_{\rm Mm} \equiv (\langle B_x \rangle_{\rm v}^2 + \langle B_y \rangle_{\rm v}^2)/2\mu_0$ (the energy of mean magnetic
components) are shown by solid, dashed and dash-dotted lines in Figure 2a. The evolution of the $\epsilon_{\rm K}$ of the progenitor run
before starting the dynamo run is also shown. Note that, from the horizontal symmetry and ${\rm div} \bm{B} = 0$,
$\langle B_z\rangle_{\rm h}$ and $\langle B_z\rangle_{\rm v}$ are zero independent of time.

After a short relaxation time, the convective kinetic energy reaches a quasi-steady state at $t \simeq 50\tau_{\rm cv}$. The magnetic
energy is gradually amplified by the convection and is saturated at $t \simeq 120\tau_{\rm cv}$. The mean values evaluated there are
$u_{\rm cv} = 5.9\times 10^{-3}$ and $B_{\rm eq,v} \equiv \sqrt{\langle\langle \mu_0 \rho {\bm u}^2 \rangle\rangle_{\rm v}}
= 7.2\times 10^{-2}$, providing $\tau_{\rm cv} = 54$ and ${\rm Ro} = 0.02$.
The vertical profile of $\langle\langle u_z^2 \rangle\rangle_{\rm h}^{1/2}$ (solid) and the distributions of the $u_z$ on the horizontal
planes in the dynamo-saturated stage are also shown in Figure 1b and 1c. Since the rotation gives rise to the Coriolis force
acting on the convective motion, the convective cell shrinks and thus the scale separation becomes larger in the rotating system.
Since the mean kinetic helicity, which is a prerequisite for exciting the large-scale dynamo, arises as a natural consequence
of the rotation, the dynamo-generated magnetic field also affects the convective motion shown in Figure 1c.

Figure 2b shows the time-depth diagrams of $\langle B_x \rangle_{\rm h}$ and $\langle B_y \rangle_{\rm h}$ normalized by
$B_{\rm eq} (z)$. Note that the turbulent magnetic component is eliminated by taking horizontal average. It is found that the oscillatory large-scale
horizontal magnetic component is spontaneously organized in the bulk of the CZ. It has a peak with the super-equipartition strength in the mid-part
of the CZ and propagates from there to the top and base of the CZ. Since there exists a phase difference of $\pi/2$ between
$\langle B_x \rangle_{\rm h}$ and $\langle B_y \rangle_{\rm h}$, the mean horizontal magnetic flux, defined by
$B_h \equiv \sqrt{\langle B_x \rangle_{\rm h}^2 + \langle B_y \rangle_{\rm h}^2}$, has a quasi-steady vertical profile.

The large-scale dynamo observed here in the strongly-stratified model has physical properties similar to those in the weakly-stratified
convective dynamo simulations \citep[e.g.,][MS14a]{kapyla+13}. Because of the horizontal symmetry and thus no differential rotation
in our system, the turbulent electro-motive force would be solely responsible for the dynamo (see MS14b for a mean-field $\alpha^2$-dynamo
model which can quantitatively reproduce the DNS results). Our intriguing finding in this Letter, which have not been observed in the
weakly-stratified model with the similar boundary conditions, is spontaneous formation of large-scale magnetic structures
in the CZ surface, which will be reported in the following.
\subsection{Spontaneous Formation of Surface Magnetic Structure}
\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=17cm,clip]{f3.pdf}
\caption{(a) A series of snapshots for the horizontal distribution of the $B_z$ at $z/d_{\rm cz} = 0.04$ and vertical distribution of the
$\langle B_z \rangle_{\rm y}/B_{\rm eq}(z)$ or $\langle B_z \rangle_x/B_{\rm eq}(z)$. (b) The temporal evolution of 2D Fourier
spectrum of the $B_z^2$ in the upper CZ, where $k_c = 2\pi/W$.}
\label{fig3}
\end{center}
\end{figure*}
A series of snapshots where the distribution of the $B_z$ at different time on the horizontal cutting plane at $z/d_{\rm cz} = 0.04$ is shown
in the top panel of Figure 3a. The darker (lighter) tone denotes positive (negative) $B_z$. While the $B_z$ has a small-scale tangled structure
with the typical size comparable to the convective cell in the initial evolutionary stage [(a1)--(a2)], it evolves as time passes to organize the large-scale
structure with the spatial-scale much larger than the convective cell [(a3)--(a4)]. The surface magnetic structure has the dynamically-important
strength comparable to $B_{\rm eq}(z)$ and is recurrently-appeared in the dynamo-saturated stage [(a5)], implying that it should be a
consequence of the strongly-stratified MHD convection rather than a transiently-formed structure.

The one-sided horizontal average ($x$- or $y$-average) of the $B_z$ may be helpful to clearly demonstrate its
depth-dependence because the horizontal average of it, i.e., $\langle B_z \rangle_{\rm h}$, is zero in our setup. Shown in
lower panel of Figure 3a is a corresponding time series of snapshots for the distributions of the $\langle B_z \rangle_{\rm y}$
[(a1)--(a4)] or $\langle B_z \rangle_{\rm x}$ [(a5)] in units of $B_{\rm eq} (z)$. The blue (red) tone represents the positive (negative) value.
The initial small-scale vertical magnetic structure gradually evolves to the larger scale in the upper CZ. While it prevails in the mid-part of the
CZ in the dynamo-saturated stage, the core of the large-scale $B_z$ structure stays in around the CZ surface.

Figure 3b shows the temporal evolution of two-dimensional (2D) Fourier spectrum of the $B_z^2$ in the upper CZ. Here the spectrum at each
depth is projected onto a one-dimensional wavenumber $k^2 = k_x^2 + k_y^2$ and then is averaged over the normalized depth from $0.0$
to $0.25$. The energy-containing scale of the $B_z$ becomes larger with the time and finally reaches the possible largest scale in the
upper CZ, confirming the existence of the large-scale $B_z$ structure in the upper CZ. Since the large-scale vertical magnetic
component cannot be generated solely by the $\alpha^2$-type dynamo only with the vertical helicity variation (see, MS14b), some additional
mechanism should play a role in organizing it in this simulation.
\subsection{What Makes Surface Magnetic Structure Organized ?}
\begin{figure*}[tbp]
\begin{center}
\begin{minipage}{0.45\hsize}
\begin{center}
\includegraphics[width=7.5cm,clip]{f4a.pdf}
\end{center}
\end{minipage}
\begin{minipage}{0.45\hsize}
\begin{center}
\includegraphics[width=7.9cm,clip]{f4b.pdf}
\end{center}
\end{minipage}
\caption{(a) The growth rate of the Parker instability evaluated from the simulation result. (b) The comparison between the growth time of
the fastest-growing mode of the Parker instability ($\tau_{\rm pk,m}$) and the local convective turn-over time ($\tau_{\rm cv,l}$)
as a function of the depth.}
\label{fig5}
\end{center}
\end{figure*}
The so-called negative magnetic pressure instability (NEMPI) has been proposed as a mechanism for the self-assembly of
the magnetic flux near the CZ surface \citep[e.g.,][]{kleeorin+96,brandenburg+10}. Although its presence has been confirmed numerically
in the forced MHD turbulence \citep[][]{brandenburg+11,brandenburg+13,mitra+14,jabbari+14,warnecke+13}, it would not play a
significant role in organizing the surface magnetic structure in our simulation. This is because the relatively-rapid rotation with ${\rm Ro} =0.02$
is assumed here, but according to \citet{lasoda+12}, ${\rm Ro} \gtrsim 5 $ is required for exciting the NEMPI. Here possible two alternatives
are discussed as a cause of the magnetic structure formation: one is the Parker instability \citep[][]{parker79} and the other is the flux
expulsion accompanied with the strong downflow \citep[e.g.,][]{weiss66, kitiashvili+10, stein+12, kapyla+15}.

Figure 4a shows the growth rate of the Parker instability obtained from the WKB dispersion equation
\begin{eqnarray}
(c_\star^2 + v_{A}^2 )\gamma^4 & + & v_{A}^2\left[ 2(c_\star^2 + v_A^2)k_h^2 + g_0 \frac{D}{Dz}\ln\left(\frac{B_h}{\rho}\right)\right] \gamma^2  \\ \nonumber
& + & k_h^2v_A^4\left( k_h^2 c_\star^2 + g_0\frac{D}{Dz}\ln B_h \right)  = 0 \;,
\end{eqnarray}
\citep[e.g.,][]{gilman70}, where $\gamma$ is the growth rate, $k_h$ is the horizontal wavenumber, $B_h$ is the horizontal magnetic flux,
$c_\star \equiv \sqrt{P/\rho} $, and $v_{A} \equiv B_h/\sqrt{\mu_0 \rho}$. For deriving the depth-dependent growth rate, we adopt
$P = \langle P \rangle_{\rm h}$, $\rho = \langle \rho \rangle_{\rm h}$ and $B_h = \sqrt{\langle B_x \rangle_{\rm h}^2 + \langle B_y \rangle_{\rm h}^2}$
evaluated from the simulation model. The different line-type denotes the growth rate at the different depth. The vertical and horizontal axes are
normalized by $1/\tau_{\rm cv}$ and $k_c = 2\pi/W$. Note that the instability is inhibited in the range $k/k_c < 1$ due to the box-width constraint.

The dynamo-maintained magnetic flux is unstable to the Parker instability in the span $0 \le z/d_{\rm cz} \lesssim 0.4$. Since the
typical growth time of it is comparable to or a bit smaller than $\tau_{\rm cv}$, it has a sufficient time to grow during the simulation.
However, the most unstable mode has a smaller wavelength in comparison with the box-width, implying the difficulty to explain the dominance of the
box-sized surface magnetic structure in our simulation.

Additionally to the mismatch of the magnetic spatial-scales, the growth of the Parker instability itself may be inhibited by vigorous
convective motions. In Figure 4b, the growth time of the most unstable mode of the Parker instability (solid: $\tau_{\rm pk,max}$) and the local
convective turn-over time defined by $\tau_{\rm cv, l} (z) = H_{\rm \rho} (z)/\langle u_{z}^2 \rangle_{\rm h}^{1/2}$ (dashed) are compared as a
function of the depth, where $H_{\rho} \equiv {\rm d}z/{\rm d}\ln\langle\rho\rangle_{\rm h}$.  In the upper CZ, the condition
$\tau_{\rm cv, l} \ll \tau_{\rm pk,max}$ is always satisfied. Since, in such a situation, the small-scale convective motion violently disturbs the
coherency of the magnetic flux, we would have to say that the Parker instability would not be responsible for the large-scale structure
formation observed in our simulation.

Next, the large-scale flow and its association with the surface magnetic structure are analyzed. For casting light on the
large-scale pattern, the small-scale structures with $k/k_c \gtrsim 8$ are eliminated by applying Fourier filtering \citep[e.g.,][]{warnecke+15,jabbari+16}.
A series of snapshots where $\langle B_z \rangle_{k_8}$ and $\langle u_z \rangle_{k_8}$ on the horizontal plane at $z/d_{\rm cz} = 0.04$ are shown in
Figure 5a, where $\langle \cdot \rangle_{k_8}$ denotes Fourier filtering. The over-plotted arrows are the velocity vectors composed of
$\langle u_x \rangle_{k_8}$ and $\langle u_y \rangle_{k_8}$. Additionally, 2D spectra of $B_h^2$, $B_z^2$, $\rho v_h^2$ and $\rho v_z^2$
in the upper CZ are also shown in Figure 5b. The spectrum at each depth is spatially-averaged over the normalized depth from $0.0$ to
$0.25$ and is temporally-averaged over $10\tau_{\rm cv}$ around the corresponding reference time.

It is prominent that, in the dynamo-saturated stage, the bipolar ``band-like" structure elongated along the direction of the
horizontal magnetic flux is predominant (see Fig.2b). Although the faster horizontal flow and stronger downflow can be found in/around the
region with the stronger $B_z$ before the dynamo-saturation (left column), large-scale flow pattern is not necessarily associated with the
magnetic structure in the dynamo-saturated stage (middle and right columns). In addition, we can find from the spectra that the energy
contained in the large-scale magnetic components is much larger than that of the large-scale flow in the upper CZ. This would suggest
that the large-scale flows are not the cause, but a consequence of the large-scale magnetic structures in the upper CZ.

Overall our analyses indicate that there should be an as-yet-unknown mechanism for the self-organization of large-scale
magnetic structures, which would be inherent in the strongly-stratified atmosphere. The band-like magnetic structure observed in our
simulation is similar to that observed in the large-scale dynamo by the forced turbulence in a strongly-stratified atmosphere
\citep[][]{mitra+14,jabbari+16}. This may imply that the surface magnetic structure formation is a common universal
feature of the strongly-stratified MHD turbulence regardless of its details.
\begin{figure*}[tbp]
\begin{center}
\includegraphics[width=18cm,clip]{f5.pdf}
\caption{
(a) A series of snapshots of $\langle B_z \rangle_{k_8}$ and $\langle u_z \rangle_{k_8}$ on the horizontal plane at $z/d_{\rm cz} = 0.04$.
The over-plotted arrows are the velocity vectors. (b) 2D spectra of $B_h^2$, $B_z^2$, $\rho v_h^2$ and $\rho v_z^2$ in the upper CZ
at around the reference time $t_r$. }
\label{fig5}
\end{center}
\end{figure*}
\section{Summary}
In this Letter, we studied numerically MHD convections in the strongly-stratified atmosphere resembling the solar interior.
The large-scale dynamo observed in our simulation had physical properties similar to those observed in the weakly-stratified model (see
MS14ab): the oscillatory large-scale horizontal magnetic component with the dynamically-important strength was spontaneously organized
in the bulk of the CZ. Its spatio-temporal evolution strongly suggests that the $\alpha^2$-type mechanism would be responsible for the MHD
dynamo in our system.

Our intriguing finding, which have not been observed in the weakly-stratified model, was the spontaneous formation of the large-scale
$B_z$ structure in the CZ surface. Small-scale tangled components of the $B_z$ which was seen in the earlier evolutionary stage gradually
evolved to the large-scale organized structure with the size much larger than convective cell as time passes in the upper CZ. In the
dynamo-saturated stage, the bipolar ``band-like" $B_z$ structure was predominant and was recurrently-appeared. Since the possible
candidates, such as the NEMPI, Parker instability, and flux expulsion due to the strong downflow, had difficulties to explain the surface
magnetic structure formation, our results may suggest the existence of an as-yet-unrecognized mechanism inducing the spontaneous
formation of the large-scale magnetic structure in the strongly-stratified convective atmosphere.

The large-scale dynamo observed here maintains the quasi-steady magnetic flux and thus is different from the actual solar dynamo with a
quasi-periodic flux modulation. Furthermore, since we adopt the faster rotation than that achieved in the actual Sun and ignore the global
effects, such as differential rotation and meridional flows, our model still remains a long way from the solar interior.  However, despite a lack
of some solar elements, it would be interesting that the relatively shallow root of the surface magnetic structure, implied from our simulation,
seems to be compatible with some observations of the magnetic patches on the solar surface \citep[e.g.,][]{brandenburg05,hara09}.
Understanding the self-organization mechanism of the surface magnetic structure observed in our study should deepen our knowledge
of the solar magnetism and is a priority target of our future work.
\acknowledgments
We acknowledge an anonymous referee for constructive comments. Computations were carried out on XC30 at NAOJ. This work
was supported by JSPS KAKENHI Grant number 15K17611 and the joint research project of the Institute of Laser Engineering,
Osaka University.