text:
# Relaxation schemes for normal modes of magnetic vortices: Application to S-matrix

## I Introduction: Vortex Normal Modes

The spectrum of small-amplitude vibrations in the presence of a nonlinear excitation, such as a vortex, soliton, or other inhomogeneous magnetization, contains information about the scattering properties, the translational properties, and the internal modes of vibration and instabilities of that object. Any modes that acquire zero frequency, for example, as a parameter is changed, signal a change of symmetry or phase transition, or perhaps a process such as magnetization reversal. Other modes may be coupled to the translation of the center of the object: such translation modes contain information about the effective mass of the excitation. If the modes are calculated by diagonalization methods, the size of the matrix and cpu time necessary to solve it very easily exceed the practical limits of any computer even for fairly small systems. Thus it is interesting to consider other methods for the calculation of at least the lowest frequency modes in the presence of an inhomogeneous state.

In particular, classical magnetic models are known to support localized or at least partially localized nonlinear excitations, one of the most well-known examples being the vortices of two-dimensional (2D) models with easy-plane (XY) anisotropy:

$`H=`$ $``$ $`{\displaystyle \frac{J}{2}}{\displaystyle \underset{𝐧,𝐚}{}}(S_𝐧^xS_{𝐧+𝐚}^x+S_𝐧^yS_{𝐧+𝐚}^y+\lambda S_𝐧^zS_{𝐧+𝐚}^z)`$ (1)
$``$ $`\stackrel{}{h}{\displaystyle \underset{𝐧}{}}\stackrel{}{S}_𝐧,`$ (2)
where subscripts $`𝐧`$ and $`𝐚`$ label the lattice sites and displacements to the nearest neighbors, and $`0\lambda <1`$ determines the anisotropy, and we have included an applied magnetic field $`\stackrel{}{h}`$. Quite generally, linearization of

the associated equations of motion about a nonlinear excitation such as a vortex will lead to a normal mode or spinwave problem, with an operator $``$ producing the time evolution of a spinwave wavefunction $`\mathrm{\Psi }`$:

$$\frac{d}{dt}\mathrm{\Psi }=\mathrm{\Psi }.$$
(3)
The calculation of the spectrum of spinwaves in the presence of a vortex by diagonalization was accomplished for small circular systems up to a radius $`R20a`$, where $`a`$ is the lattice constant. The method was used to determine vortices’ spinwave spectrum for ferromagnets (FM, $`J>0`$) and antiferromagnets (AFM, $`J<0`$) on a lattice. In addition to usual scattering states, it was found that FM and AFM vortices possess an internal mode whose frequency goes to zero at a critical anisotropy value $`\lambda _c`$, which is lattice dependent ($`0.7034`$ for square lattice). The mode is associated with the instability of out-of-plane vortices (those with $`S^z0`$) to change to in-plane vortices $`S^z=0`$ everywhere) at strong easy-plane anisotropy ($`\lambda <\lambda _c`$). Other modes found, which tend to have high intensity at the vortex core, are associated with the translation of the vortex center of mass, and can be related to a collective coordinate description for vortex dynamics.

Analytic calculations of the vortex normal modes, using a continuum limit, have partially described the scattering spectrum. In the earliest calculations, the phase shifts of spinwaves scattered by FM in-plane vortices were calculated, using a Born approximation, for the exchange anisotropy model (1) considered here  and for a similar FM model with site-anisotropy. The weak point of these calculations is that they rely on the choice of a short distance cutoff for the Born approximation integrals. For the AFM in-plane vortices, Pereira et al. were able to make an exact calculation of the scattering phase shifts, for the optical spinwave branch, without the need for a cutoff. In the continuum limit equations, AFM vortices do not scatter the acoustic branch spinwaves; the associated phase shifts vanish.

The validity of the continuum limit and other approximations, however, can only be tested by comparison with numerical calculations. For example, Ivanov et al.  found by analytic means that the vortices of the continuum AFM model for $`\lambda `$ near 1 should possess a true internal local mode. It was shown only through numerical diagonalization that this mode is identical to the vortex instability mode mentioned above. Ivanov et al. were also able to calculate the scattering (S-matrix) by a numerical shooting solution for the continuum theory equations of motion. The S-matrix was not calculated from the numerical diagonalization data, however, because the system size that could be solved was too small.

For $`\lambda _c<\lambda 1`$, the nonzero $`S^z`$ component of the out-of-plane vortex leads to a smooth structure decaying over the length scale

$$r_v=\sqrt{\frac{\lambda }{1\lambda }},$$
(4)
with no singularity at the vortex center. Thus, the continuum description used in Ref.  is quite accurate. At the other limit, $`\lambda <\lambda _c`$, the in-plane vortex spin field is singular at the vortex center. In this case, the continuum limit cannot describe very well the spin dynamical motions near the vortex core, without the introduction of some kind of lower radius cutoff. The length scale in the problem approaches zero, and derivatives of the static vortex spin structure on a lattice are not small in the core region, violating the usual continuum limit assumptions about small spatial derivatives across one lattice constant. Nevertheless, even for the in-plane AFM vortices, a true local mode of the vortex is present, and can be approximately described in continuum theory, but only if a cutoff is introduced. Therefore, for both FM and AFM in-plane vortices, it is interesting to consider the calculation of the S-matrix (or equivalently, phase shifts) for the lattice system, and compare with continuum scattering theory.

It has been realized that the translational modes that come from the spinwave spectra have unusual properties and can be used to calculate a vortex effective mass. Although only a limited size of system could be studied ($`R<20a`$) it was found that the FM vortex mass increases faster than $`\mathrm{ln}R`$, for in-plane as well as out-of-plane vortices. AFM in-plane vortices also were found to have a mass increasing faster than $`\mathrm{ln}R`$, however, it appeared that the mass of AFM out-of-plane vortices may actually reach a finite limit for large $`R`$. It is also possible that AFM vortex mass tends to zero for $`\lambda 1`$, the isotropic limit.

More recently, a newer theory for the FM out-of-plane vortex mass and related collective coordinate description has been developed by Mertens et al. There, a vortex mass approximately independent of system size was determined, both by including effects due to an image vortex outside the studied system, and also, by using a dynamical equation of motion higher than second order in time. This newer theory appears to describe well the individual vortex dynamics by assigning a mass, gyrovector, and a new higher order gyro-tensor to the vortex, considered as a particle-like object. A careful application of the theory was made by Ivanov et al, using spinwave translation mode data obtained from another numerical scheme for finding some of the lowest modes. In Ref. , Wielandt’s version of inverse iteration procedure was used: an operator $`(\omega I)^1`$ was applied to an initial randomly chosen vector, which after many iterations converges to an eigenvector of $``$, provided the frequency $`\omega `$ is updated appropriately. It is also a type of relaxation procedure similar to what we describe here, with great speed and memory advantages over the diagonalization method. This vortex mass theory, however, does not apply directly to the AFM vortices, nor to the in-plane vortices for both FM and AFM models.

In this paper we present and apply some numerical methods for finding the normal modes of a magnetic vortex on a lattice, which are applicable to systems larger than those solvable by diagonalization schemes. We use these methods to analyze the FM and AFM spinwave-vortex scattering problems for in-plane vortices. In addition to the calculation of the S-matrix and phase shifts, we present results on the translational modes and make estimates of in-plane vortex mass. Analysis of out-of-plane vortices will be considered elsewhere.

The paper is organized as follows. First, in Sec. II, a short overview of relaxation schemes is given. In Sec. III, we describe the magnetic spinwave equations of motion in the presence of a nonuniformly magnetized configuration. This is followed in Sec. IV by a description of the relaxation schemes using time evolution and Gauss-Seidel iteration. In Sec. VI the schemes are used to calculate the vortex-spinwave S-matrix, for a single vortex in a circular system. After a summary of the translational dynamics of magnetic vortices, in Sec. VII we present and analyze the results obtained for the FM and AFM in-plane vortex mass.

## II Relaxation Methods for Eigenfunctions

The basic idea of these relaxation schemes is the following: The eigenfrequencies of a linear spinwave problem, such as Eq. (3), always come in positive/negative pairs, corresponding essentially to creation and annihilation operators, respectively. Thus, its time evolution leads to oscillatory behavior, starting from an arbitrary initial condition. It is interesting, however, to consider how the evolution can be made diffusive, such that high (absolute magnitude) frequency modes decay away faster than low frequency modes. For example, a Schrödinger equation with the replacement $`it\tau `$, becomes a diffusion equation in imaginary time $`\tau `$ (because it is a parabolic equation). It is well known that its imaginary time evolution will lead to the lowest frequency wavefunction, i.e., the ground state of the associated Hamiltonian. In the equations of motion (3), such a replacement $`it\tau `$, unfortunately, does not lead to a diffusion-like equation, due to the presence of the two fundamental solutions (creation/annihilation operators) behaving as $`e^{\pm i\omega _kt}e^{\pm \omega _k\tau }`$ (because the spinwave equations are essentially hyperbolic). There is a fundamental solution growing with $`\tau `$ for every mode $`k`$, leading to numerical instability in imaginary time. However, if one uses the square of the spinwave Hamiltonian to evolve forward in real time, then both fundamental solutions associated with mode $`k`$ decay as $`e^{\omega _k^2t}`$. Evolution forward in time will lead to the lowest frequency eigenmode. By combining with Gram-Schmidt orthogonalization, it is then possible to iterate this procedure and generate a small set of the lowest frequency eigenmodes. The method offers a drastic savings in memory needed compared to numerical diagonalization, at the expense of rather long cpu time due to its diffusive behavior. Also note, it is not essential that the evolution be performed in time. We have also found that the eigenmodes can be found even more quickly by using a Gauss-Seidel iteration of the spinwave equations, where the unknown frequency $`\omega `$ within those equations is determined self-consistently during the iteration.

## III 2D Easy-Plane Magnetic Model and Spinwave Equations

We review briefly the derivation of the linearized spinwave problem for normal modes on a vortex, or other nonuniform state. The nonlinear equations of motion from (1) are simple:

$$\dot{𝐒}_𝐧=𝐒_𝐧\times \left[\stackrel{}{h}+\underset{𝐧^{}=𝐧+𝐚}{}\stackrel{~}{J}𝐒_𝐧^{}\right],$$
(5)
where $`\stackrel{~}{J}`$ is a constant diagonal exchange matrix:

$$\stackrel{~}{J}=J\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& \lambda \end{array}\right)$$
(6)
and the spins are considered as column vectors. A vortex is a particular static configuration $`\{𝐒_𝐧^0\}`$ that solves this equation with a zero time derivative. The structure of in-plane and out-of-plane vortices has been described elsewhere using continuum theory and in the presence of a lattice. Because the spin length is a conserved quantity, there are effectively only two degrees of freedom per spin, and it is usually convenient to describe the static structure using an in-plane angle $`\varphi _𝐧^0`$ and an out-of-plane angle $`\theta _𝐧^0`$, in planar spherical coordinates:

$$𝐒_𝐧^0=S\left(\begin{array}{c}\mathrm{cos}\theta _𝐧^0\mathrm{cos}\varphi _𝐧^0\\ \mathrm{cos}\theta _𝐧^0\mathrm{sin}\varphi _𝐧^0\\ \mathrm{sin}\theta _𝐧^0\end{array}\right).$$
(7)
Then one could obtain the spinwave effective equations of motion by assuming a perturbation to this structure in the form,

$$\varphi _𝐧=\varphi _𝐧^0+\phi _𝐧,$$
(9)
$$\theta _𝐧=\theta _𝐧^0+\vartheta _𝐧,$$
(10)
where the equations of motion are now linearized in terms of $`\phi _𝐧`$ and $`\vartheta _𝐧`$.

In Ref.  a slightly different but equivalent description was used, which we follow here. The unperturbed spins of the static vortex structure are considered to define local quantization axes $`\stackrel{~}{z}_𝐧`$, specifically,

$$𝐒_𝐧^0=S\stackrel{~}{z}_𝐧.$$
(11)
Then the perturbation of this structure involves fluctuations orthogonal to the $`\stackrel{~}{z}_𝐧`$-axis, along two new local axes $`\stackrel{~}{x}_𝐧`$ and $`\stackrel{~}{y}_𝐧`$. The $`\stackrel{~}{x}_𝐧`$-axis is taken to be along the direction defined by the cross product $`\stackrel{~}{x}_𝐧=z_𝐧\times \stackrel{~}{z}_𝐧`$, which is within the original xy (easy) plane. Quantities $`x_𝐧`$, $`y_𝐧`$, and $`z_𝐧`$ here refer to the original coordinates of Hamiltonian (1). Then $`\stackrel{~}{y}_𝐧=\stackrel{~}{z}_𝐧\times \stackrel{~}{x}_𝐧`$ makes the last member of the local orthogonal set. The perturbation of the static vortex structure can then be expressed in terms of its spin components along these new local axes:

$$𝐒_𝐧=S_𝐧^{\stackrel{~}{z}}\stackrel{~}{z}_𝐧+S_𝐧^{\stackrel{~}{x}}\stackrel{~}{x}_𝐧+S_𝐧^{\stackrel{~}{y}}\stackrel{~}{y}_𝐧.$$
(12)
A short calculation shows that these are related to the angular perturbation coordinates by

$$S_𝐧^{\stackrel{~}{x}}=S\phi _𝐧\mathrm{cos}\theta _𝐧^0,S_𝐧^{\stackrel{~}{y}}=S\vartheta _𝐧.$$
(13)
The variables $`S_𝐧^{\stackrel{~}{x}}`$ and $`\phi _𝐧`$ relate to purely in-plane spin motions, while $`S_𝐧^{\stackrel{~}{y}}`$ and $`\vartheta _𝐧`$ measure the out-of-easy-plane tilting, relative to the unperturbed vortex.

In Refs.  and  it is shown that under the assumptions $`S_𝐧^{\stackrel{~}{x}}1`$, $`S_𝐧^{\stackrel{~}{y}}1`$, $`S_𝐧^{\stackrel{~}{z}}S`$, this notation leads to a usual linear time evolution problem, which we write here in the form:

$$\frac{d}{dt}\left(\begin{array}{c}S_𝐧^{\stackrel{~}{x}}\\ S_𝐧^{\stackrel{~}{y}}\end{array}\right)=\underset{𝐦=𝐧,𝐧^{}}{}\left(\begin{array}{cc}M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{x}\stackrel{~}{x}}& M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{x}\stackrel{~}{y}}\\ M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{y}\stackrel{~}{x}}& M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{y}\stackrel{~}{y}}\end{array}\right)\left(\begin{array}{c}S_𝐦^{\stackrel{~}{x}}\\ S_𝐦^{\stackrel{~}{y}}\end{array}\right)$$
(14)
or equivalently in components,

$$\dot{S}_𝐧^\alpha =\underset{𝐦=𝐧,𝐧^{}}{}\underset{\beta }{}M_{\mathrm{𝐧𝐦}}^{\alpha \beta }S_𝐦^\beta $$
(15)
where $`\alpha `$ and $`\beta `$ range over the set $`\{\stackrel{~}{x},\stackrel{~}{y}\}`$, and the site $`𝐦`$ must either be equal to site $`𝐧`$ or one of its neighbors, $`𝐧^{}=𝐧+𝐚`$. The matrix elements $`M_{\mathrm{𝐧𝐦}}^{\alpha \beta }`$ are:

$`M`$ $`{}_{\mathrm{𝐧𝐧}}{}^{\stackrel{~}{x}\stackrel{~}{x}}=M_{\mathrm{𝐧𝐧}}^{\stackrel{~}{y}\stackrel{~}{y}}=0,`$ (17)
$`M`$ $`{}_{\mathrm{𝐧𝐧}}{}^{\stackrel{~}{x}\stackrel{~}{y}}=M_{\mathrm{𝐧𝐧}}^{\stackrel{~}{y}\stackrel{~}{x}}=(h_x\mathrm{cos}\varphi _𝐧^0+h_y\mathrm{sin}\varphi _𝐧^0)p_𝐧^0+h_zm_𝐧^0`$ (18)
$`+`$ $`JS{\displaystyle \underset{𝐧^{}=𝐧+𝐚}{}}\left[p_𝐧^0p_𝐧^{}^0\mathrm{cos}(\varphi _𝐧^0\varphi _𝐧^{}^0)+\lambda m_𝐧^0m_𝐧^{}^0\right],`$ (19)
$`M`$ $`{}_{𝐧,𝐧^{}}{}^{\stackrel{~}{x}\stackrel{~}{x}}=JSm_𝐧^0\mathrm{sin}(\varphi _𝐧^0\varphi _𝐧^{}^0),`$ (20)
$`M`$ $`{}_{𝐧,𝐧^{}}{}^{\stackrel{~}{y}\stackrel{~}{y}}=JSm_𝐧^{}^0\mathrm{sin}(\varphi _𝐧^0\varphi _𝐧^{}^0),`$ (21)
$`M`$ $`{}_{𝐧,𝐧^{}}{}^{\stackrel{~}{x}\stackrel{~}{y}}=JS[m_𝐧^0m_𝐧^{}^0\mathrm{cos}(\varphi _𝐧^0\varphi _𝐧^{}^0)+\lambda p_𝐧^0p_𝐧^{}^0]`$ (22)
$`M`$ $`{}_{𝐧,𝐧^{}}{}^{\stackrel{~}{y}\stackrel{~}{x}}=JS\mathrm{cos}(\varphi _𝐧^0\varphi _𝐧^{}^0).`$ (23)
For simplification, we use the definitions of in-plane projection and out-of-plane magnetization,

$`p_𝐧^0`$ $``$ $`\sqrt{1(S_𝐧^{z0}/S)^2}=\mathrm{cos}\theta _𝐧^0,`$ (25)
$`m_𝐧^0`$ $``$ $`S_𝐧^{z0}/S=\mathrm{sin}\theta _𝐧^0.`$ (26)
The operator on the right hand side of the linear equations of motion, Eq. (15), is not Hermitian. (For example, consider the case $`\lambda =0`$, with all $`\theta _𝐧^0=0`$.) This led to complications in the previous numerical diagonalization calculations. Its eigenvalues are pure imaginary, in complex conjugate pairs; this results in the usual spinwave solutions oscillatory in time. However, this implies that its square (which need not be Hermitian) has pure real eigenvalues, a fact we will exploit below in Sec. IV.

In Refs.  and , the variables $`S^{\stackrel{~}{x}}`$ and $`S^{\stackrel{~}{y}}`$ were considered as quantum operators, in a semiclassical picture. Then an eigenmode labeled by index $`k`$, has a creation operator $`B_k^{}`$ that is an appropriate linear combination, with complex expansion coefficients, $`w_{k,𝐧}^\mathit{1},w_{k,𝐧}^\mathit{2}`$,

$$B_k^{}=\underset{𝐧}{}\left[w_{k,𝐧}^\mathit{1}S_𝐧^{\stackrel{~}{x}}+w_{k,𝐧}^\mathit{2}S_𝐧^{\stackrel{~}{y}}\right],$$
(27)
such that its time evolution is simple, with eigenfrequency $`\omega _k`$:

$$\dot{B}_k^{}=i\omega _kB_k^{}.$$
(28)
The corresponding annihilation operator is $`B_k=(B_k^{})^{}`$, since the spin operators are real. The local spin fluctuations are given from the inverse relationship,

$$S_𝐧^{\stackrel{~}{x}}=i\mathrm{}S\underset{k}{}\left(w_{k,𝐧}^\mathit{2}B_kw_{k,𝐧}^\mathit{2}B_k^{}\right),$$
(30)
$$S_𝐧^{\stackrel{~}{y}}=i\mathrm{}S\underset{k}{}\left(w_{k,𝐧}^\mathit{1}B_kw_{k,𝐧}^\mathit{1}B_k^{}\right),$$
(31)
Note that the italic superscripts on the $`w_𝐧^\alpha ,\alpha =\mathit{1},\mathit{2}`$ in this article should not be confused with powers, because this is a linear problem.

The collection of coefficients describes the wavefunction, and satisfies

$$i\omega _k\left(\begin{array}{c}w_{k,𝐧}^\mathit{1}\\ w_{k,𝐧}^\mathit{2}\end{array}\right)=\underset{𝐦=𝐧,𝐧^{}}{}\left(\begin{array}{cc}M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{x}\stackrel{~}{x}}& M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{y}\stackrel{~}{x}}\\ M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{x}\stackrel{~}{y}}& M_{\mathrm{𝐧𝐦}}^{\stackrel{~}{y}\stackrel{~}{y}}\end{array}\right)\left(\begin{array}{c}w_{k,𝐦}^\mathit{1}\\ w_{k,𝐦}^\mathit{2}\end{array}\right)$$
(32)
The matrix on the RHS is the transpose of that in Eq. (14). Thus the two systems clearly have the same eigenvalues, however, we prefer to use the $`w^\mathit{1},w^\mathit{2}`$ notation, together with its underlying creation and annihilation operators, because the overall wavefunction orthogonalizations are determined by related commutators. A wavefunction can be considered as a column vector of these coefficients:

$$\mathrm{\Psi }_k=\left(\begin{array}{c}\psi _{k,\mathrm{𝟏}}\\ \psi _{k,\mathrm{𝟐}}\\ .\\ .\\ \psi _{k,𝐍}\end{array}\right),$$
(33)
where some labeling from 1 to $`N`$ has been given to the $`N`$ sites, each of which has a local two-component wavefunction holding the $`S^{\stackrel{~}{y}}`$ ($`w^\mathit{1}`$) and $`S^{\stackrel{~}{x}}`$ ($`w^\mathit{2}`$) fluctuations:

$$\psi _{k,𝐧}=\left(\begin{array}{c}w_{k,𝐧}^\mathit{1}\\ w_{k,𝐧}^\mathit{2}\end{array}\right).$$
(34)
The dynamics of these individual wavefunctions can be written in terms of the matrices appearing in Eq. (32), which we denote here for convenience as $`_{\mathrm{𝐧𝐦}}`$:

$$i\omega _k\psi _{k,𝐧}=\underset{𝐦=𝐧,𝐧^{}}{}_{\mathrm{𝐧𝐦}}\psi _{k,𝐦}.$$
(35)
Then we can see that there exists an operator $``$ composed of the individual $`_{\mathrm{𝐧𝐦}}`$ matrices, with real elements, that evolves the entire wavefunction:

$$i\omega _k\mathrm{\Psi }_k=\mathrm{\Psi }_k.$$
(36)
This is the linear problem to be solved.

We have previously solved this non-Hermitian problem, Eq. (36), through numerical diagonalization, using the EISPACK routine RG() for diagonalizing real-general matrices (no symmetry assumed). In this way the full spectrum could be found, however, the matrix $``$ is sparse but of size $`2N\times 2N`$, where $`N`$ is the number of lattice sites. A circular square lattice system with radius $`R=20a`$ has $`N=1264`$. This leads to a matrix with 6390784 elements, which when stored in double precision (8 bytes per number), requires 49 MB RAM. In practice, it is also necessary to store another array of equal size that contains the calculated eigenmodes, thus the required storage for a calculation at one value of $`\lambda `$ is twice this number. This gets incremented by another 49 MB if one wants to compare eigenmodes at nearby values of $`\lambda `$ for tracking their changes with anisotropy. For these reasons, on a machine with 128 MB RAM, the upper limit for this calculation was system radius $`R=`$19–20, depending on the desired results.

Next, we consider how the eigenmodes can be found through squaring the operator of the RHS of Eq. (36), combined with introducing a fictitious time evolution, resulting in considerable memory savings.

## IV Relaxation Method

### A Time Evolution

A mode $`\mathrm{\Psi }_k`$ that solves Eq. (36) is a fundamental solution of

$$\frac{d}{dt}\mathrm{\Psi }=\mathrm{\Psi }.$$
(37)
with time dependence $`e^{i\omega _kt}`$, determined by the corresponding eigenvalue of $``$, $`i\omega _k`$. Applying the $``$ operation again leads to an equivalent equation,

$$\frac{d^2}{dt^2}\mathrm{\Psi }=^2\mathrm{\Psi }.$$
(38)
Both of these equations have fundamental oscillatory solutions varying as $`\mathrm{\Psi }_k(t)=\mathrm{\Psi }_k(0)e^{i\omega _kt}`$, i.e., the usual spinwaves. The advantage is that operator $`^2`$ has purely real, negative eigenvalues, in contrast to the purely imaginary eigenvalues (both positive and negative imaginary parts) of operator $``$.

Now consider that we invent a fictitious time evolution, by rather arbitrarily changing the LHS of Eq. (38) to a first time derivative:

$$\frac{d}{dt}\mathrm{\Psi }=^2\mathrm{\Psi }.$$
(39)
The fundamental solutions here are determined by exponentials involving the eigenvalues of $`^2`$, which are $`\omega _k^2`$. Thus, the general solution from an arbitrary initial state is

$$\mathrm{\Psi }(t)=\underset{k}{}\left[c_k^+\mathrm{\Psi }_k^++c_k^{}\mathrm{\Psi }_k^{}\right]e^{\omega _k^2t},$$
(40)
where the coefficients $`c_k^\pm `$ are determined by the initial state, and the $`\mathrm{\Psi }_k^\pm `$ are normalized eigenfunctions of $``$, with respective eigenvalues $`\pm i\omega _k`$, corresponding to the creation and annihilation operators for mode $`k`$. (Of course, $`\mathrm{\Psi }_k^\pm `$ have the same eigenvalues, $`\omega _k^2`$, with respect to the operator $`^2`$.) All modes decay away in time, leaving the lowest frequency modes present in the initial state to dominate at large time. In practice, the time evolution will be started from a randomly chosen initial state (i.e., randomly assigned values of the set of $`\{w_𝐧^\mathit{1},w_𝐧^\mathit{2}\}`$ coefficients), to ensure a nonzero overlap with the lowest frequency mode. The total wavefunction $`\mathrm{\Psi }(t)`$ will be re-normalized to unity periodically to avoid numerical underflow, by enforcing the overlap appropriate for a non-Hermitian operator described below (Sec. IV B).

$$\mathrm{\Psi }|\mathrm{\Psi }=1.$$
(41)
We have used a second order Runge-Kutta time integration scheme to solve Eq. (38), and found that a time step of $`0.02/JS`$ is adequate to provide stability and insure convergence to an eigenstate.

At large time, the system will be relaxed to a linear combination of creation and annihilation wavefunctions:

$$\mathrm{\Psi }=c_{k_0}^+\mathrm{\Psi }_{k_0}^++c_{k_0}^{}\mathrm{\Psi }_{k_0}^{},$$
(42)
where $`k_0`$ labels the lowest mode that was present in the intial state. The frequency for this combination is evaluated in the usual way,

$$\omega _{k_0}^{}{}_{}{}^{2}=\frac{\mathrm{\Psi }|^2|\mathrm{\Psi }}{\mathrm{\Psi }|\mathrm{\Psi }}.$$
(43)
The creation and annihilation wavefunctions can be separated by an application of operator $``$:

$$\mathrm{\Psi }=i\omega _{k_0}\left(c_{k_0}^+\mathrm{\Psi }_{k_0}^+c_{k_0}^{}\mathrm{\Psi }_{k_0}^{}\right),$$
(44)
Combining the two equations (42) and (44) leads to

$$c_{k_0}^\pm \mathrm{\Psi }_{k_0}^\pm =\frac{i\omega _{k_0}\mathrm{\Psi }\pm \mathrm{\Psi }}{2i\omega _{k_0}}.$$
(45)
This finally will be followed by a normalization appropriate to the non-Hermitian operator $``$, to be described below \[Eq. (47)\].

### B Orthogonality of Modes and Gram-Schmidt Process

To obtain higher modes, it is necessary to evolve forward in time together with a Gram-Schmidt orthogonalization to any previously found lower frequency modes, denoted here by index $`\mathrm{}`$. We have found that the necessary way to do this is by making the current wavefunction $`\mathrm{\Psi }(t)`$ orthogonal to both $`\mathrm{\Psi }_{\mathrm{}}^+`$ and $`\mathrm{\Psi }_{\mathrm{}}^{}`$, rather than to a somewhat arbitrary combination such as in Eq. (42). Equivalently, both constants $`c_{\mathrm{}}^+`$ and $`c_{\mathrm{}}^{}`$ must be enforced to be zero in the current wavefunction, for all modes $`\mathrm{}`$ already determined. To obtain these two conditions, it is necessary to use the overlap appropriate to the non-Hermitian operator $``$, since $`\mathrm{\Psi }_{\mathrm{}}^+`$ and $`\mathrm{\Psi }_{\mathrm{}}^{}`$ are distinct eigenfunctions of $``$ corresponding to its distinct eigenvalues $`\pm i\omega _{\mathrm{}}`$.

The overlap of two eigenmodes of $``$ can be defined in one of two equivalent ways. The first is to realize that an annihilation ($`B_j`$) and creation ($`B_k^{}`$) operator, with associated wavefunctions $`\mathrm{\Psi }_j^{}`$ and $`\mathrm{\Psi }_k^+`$, must have a commutator,

$$[B_j,B_k^{}]=\delta _{j,k}.$$
(46)
From the definition, Eq. (27), with the help of the basic spin commutators $`[S_𝐧^{\stackrel{~}{x}},S_𝐧^{}^{\stackrel{~}{y}}]=i\delta _{𝐧,𝐧^{}}S_𝐧^{\stackrel{~}{z}}i\delta _{𝐧,𝐧^{}}S`$, this is equivalent to:

$$[B_j,B_k^{}]=iS\underset{𝐧}{}\left(w_{j,𝐧}^\mathit{1}w_{k,𝐧}^\mathit{2}w_{j,𝐧}^\mathit{2}w_{k,𝐧}^\mathit{1}\right)=\mathrm{\Psi }_j|\mathrm{\Psi }_k.$$
(47)
In general, this is the definition of the overlap used between two arbitrary states $`j`$ and $`k`$.

The second way to define the overlap, equivalent to this, is to make use of the left and right eigenvectors of $``$. The solutions $`\mathrm{\Psi }_k`$ of Eq. (36), where $``$ acts to the right, are the right eigenvectors. For each eigenvalue $`i\omega _k`$, there also exists a corresponding left eigenvector, $`{}_{k}{}^{}\mathrm{\Psi }`$ satisfying

$$i\omega _k{}_{k}{}^{}\mathrm{\Psi }={}_{k}{}^{}\mathrm{\Psi }.$$
(48)
Each creation operator as well as each annihilation operator has a distinct left eigenvector. For a Hermitian problem, the left eigenvectors are obtained by the familiar relationship, $`{}_{k}{}^{}\mathrm{\Psi }=\mathrm{\Psi }_k^{}`$, and then overlaps are in the usual form $`\mathrm{\Psi }_j|\mathrm{\Psi }_k={}_{j}{}^{}\mathrm{\Psi }\mathrm{\Psi }_k=\mathrm{\Psi }_j^{}\mathrm{\Psi }_k`$.

For this problem, although $``$ is not Hermitian, there exists a real matrix $`𝒜=𝒜^1`$, that transforms $``$ to $`^{}`$:

$$𝒜𝒜^1=^{}$$
(49)
where $`𝒜`$ is diagonal in $`2\times 2`$ submatrices:

$$𝒜=\left(\begin{array}{cccc}\alpha & 0& 0& \mathrm{}\\ 0& \alpha & 0& \mathrm{}\\ 0& 0& \alpha & \mathrm{}\\ .& .& .& \mathrm{}\\ .& .& .& \mathrm{}\end{array}\right),\alpha =\left(\begin{array}{cc}\hfill 0& \hfill 1\\ \hfill 1& \hfill 0\end{array}\right).$$
(50)
Then using Eq. (49) in Eq. (36), rearranging and comparing with Eq. (48) determines a relationship between left and right eigenvectors,

$${}_{k}{}^{}\mathrm{\Psi }=i\left(𝒜\mathrm{\Psi }_k\right)^{},$$
(51)
where the $``$ operation is the usual Hermitian conjugate (transpose of complex conjugate), and the factor of $`i`$ has been added for consistency with Eq. (47). Then the overlap between two states $`\mathrm{\Psi }_j`$ and $`\mathrm{\Psi }_k`$ is defined by the scalar product of a left with a right eigenvector:

$$\mathrm{\Psi }_j|\mathrm{\Psi }_k={}_{j}{}^{}\mathrm{\Psi }\mathrm{\Psi }_k=i\left(𝒜\mathrm{\Psi }_j\right)^{}\mathrm{\Psi }_k=i\mathrm{\Psi }_j^{}𝒜^{}\mathrm{\Psi }_k.$$
(52)
which is seen to give exactly the expression in (47).

The application of Gram-Schmidt orthogonalization to remove the previously found modes (labeled by $`\mathrm{}`$) proceeds according to replacing the current $`\mathrm{\Psi }(t)`$ by

$$\mathrm{\Psi }(t)\mathrm{\Psi }(t)\underset{\mathrm{}}{}\left[\mathrm{\Psi }_{\mathrm{}}^+|\mathrm{\Psi }(t)\mathrm{\Psi }_{\mathrm{}}^++\mathrm{\Psi }_{\mathrm{}}^{}|\mathrm{\Psi }(t)\mathrm{\Psi }_{\mathrm{}}^{}\right].$$
(53)
Combined with evolution forward in time and periodic renormalization according to Eq. (41), the wavefunction will evolve to the next low-frequency mode, of the form Eq. (42). Continuing to iterate the above procedures, a small set of the lowest frequency modes (about 5-10) can be determined with double precision accuracy in reasonable time for systems of size $`100\times 100`$.

### C Self-Consistent Gauss-Seidel Scheme

The eigenmodes of Eq. (36) alternatively can be found by a Gauss-Seidel type of iteration, not applied to that equation, but rather, to the equation with the squared operator,

$$=^2,$$
(55)
$$\mathrm{\Psi }_k=\omega _k^2\mathrm{\Psi }_k.$$
(56)
This Gauss-Seidel scheme to be described is about twice as fast (in total cpu time) as the time evolution using second order Runge Kutta described above. It is based on using the matrix elements of $``$, and in particular, one needs to single out the matrix elements that couple a site to itself ($`_{𝐧,𝐧}`$) from the intersite couplings ($`W_𝐧`$). We write this equation as a set of equations with the form,

$`\omega ^2w_𝐧^\mathit{1}`$ $`=`$ $`_{𝐧,𝐧}^{\mathit{11}}w_𝐧^\mathit{1}+_{𝐧,𝐧}^{\mathit{12}}w_𝐧^\mathit{2}+W_𝐧^\mathit{1}`$ (58)
$`\omega ^2w_𝐧^\mathit{2}`$ $`=`$ $`_{𝐧,𝐧}^{\mathit{21}}w_𝐧^\mathit{1}+_{𝐧,𝐧}^{\mathit{22}}w_𝐧^\mathit{2}+W_𝐧^\mathit{2}`$ (59)
$`W_𝐧^\alpha `$ $`=`$ $`{\displaystyle \underset{𝐧^{\prime \prime }}{}}{\displaystyle \underset{\beta =1,2}{}}_{𝐧,𝐧^{\prime \prime }}^{\alpha \beta }w_{𝐧^{\prime \prime }}^\beta ,`$ (60)
where the sites $`𝐧^{\prime \prime }`$ include all first and second nearest neighbors of site $`𝐧`$, but not the site $`𝐧`$. The second nearest neighbor interations result directly from squaring the nearest neighbor operator, $``$. These matrix elements of $``$ are sums of quadratic products,

$$_{𝐧,𝐧^{\prime \prime }}=\underset{𝐦=𝐧,𝐧^{},𝐧^{\prime \prime }}{}_{𝐧,𝐦}_{𝐦,𝐧^{\prime \prime }},$$
(61)
where the summation is over any sites $`𝐦`$ that are either first neighbors of sites $`𝐧`$ and $`𝐧^{\prime \prime }`$, or, these sites themselves. Due to the structure of the matrix elements of $``$, the on-site, off-diagonal elements $`_{𝐧,𝐧}^{\mathit{12}},_{𝐧,𝐧}^{\mathit{21}}`$, are both zero, while the corresponding diagonal elements are equal: $`_{𝐧,𝐧}^{\mathit{11}}=_{𝐧,𝐧}^{\mathit{22}}`$.

Now in the spirit of Gauss-Seidel iteration, we can try to “solve” for the new values, $`\psi _{𝐧}^{}{}_{}{}^{}`$, with some certain freedom of choice. That is, the terms involving $`\omega ^2`$ can be taken as explicit or implicit, depending on whether we consider them to involve the wavefunction at the previous step or the wavefunction at the current step being solved. We have found for this problem, that the iteration is more stable if the $`\omega ^2`$ terms are explicit, and the equation is re-arranged for each component $`\alpha =\mathit{1},\mathit{2}`$ as

$$w_{𝐧}^{\alpha }{}_{}{}^{}=\frac{1}{_{𝐧,𝐧}^{\alpha \alpha }}\left[\omega ^2w_𝐧^\alpha +W_𝐧^\alpha \right].$$
(62)
For each step of the iteration, the frequency $`\omega `$ must be evaluated according to an expectation value of the form of Eq. (43). Also, in the interest of numerical efficiency, we actually store the matrix elements of $``$ and only the diagonal elements of $`=^2`$. Then it is useful to realize that an equivalent expression for $`W_𝐧`$ is in terms of the entire $``$ operation with the diagonal part subtracted out:

$$W_𝐧^\alpha =\left[\left(\mathrm{\Psi }\right)_𝐧^\alpha _{𝐧,𝐧}^{\alpha \alpha }w_𝐧^\alpha \right]$$
(63)
This means Eq. (62) can be written as

$$w_{𝐧}^{\alpha }{}_{}{}^{}=w_𝐧^\alpha \frac{1}{_{𝐧,𝐧}^{\alpha \alpha }}\left[\omega ^2w_𝐧^\alpha +\left(\mathrm{\Psi }\right)_𝐧^\alpha \right].$$
(64)
The stability of this process is further enhanced if we mix a fraction $`f<1`$ of the updated configuration $`\mathrm{\Psi }^{}`$ \[Eq. (64)\] with a fraction $`(1f)<1`$ of the old configuration $`\mathrm{\Psi }`$. This then finally defines the iteration procedure:

$$w_{𝐧}^{\alpha }{}_{}{}^{}=w_𝐧^\alpha \frac{f}{_{𝐧,𝐧}^{\alpha \alpha }}\left[\omega ^2w_𝐧^\alpha +\left(\mathrm{\Psi }\right)_𝐧^\alpha \right].$$
(65)
The equation is to be iterated, together with Eq. (43) for the frequency and the usual normalization, Eq. (41), until a self-consistent convergence is achieved. We should note also that the application of Eq. (65) can be performed either synchronously or asynchronously. In a synchronous step, the old values $`w_𝐧^\alpha `$ are updated into new ones $`w_{𝐧}^{\alpha }{}_{}{}^{}`$ simultaneously, or in parallel. In an asynchronous step, which is the usual Gauss-Seidel scheme, once a site is updated into the new value $`w_{𝐧}^{\alpha }{}_{}{}^{}`$, that new value is recycled immediately into the right hand side of Eq. (43), and this new value then modifies the results for its neighboring sites. Either way, the value of $`\omega ^2`$ is corrected only after the entire lattice has been updated. Both ways work well, although the asynchronous step generally produces faster convergence. Just as in the time evolution scheme, a random initial configuration of $`\mathrm{\Psi }`$ is used for these iterations. In our actual calculations, we have used fractions $`f`$ from $`0.6`$ to as high as $`1.9`$ (an over-relaxation). Values of $`f`$ below $`1.0`$ provide reasonably fast convergence with reliable stability. Fractions $`f`$ that are closer to zero result in slower convergence, while fractions closer to one give faster convergence at the expense of greater instability, where the instability causes convergence to the highest frequency eigenmode rather than the desired lowest frequency mode.

The above procedure will tend to the lowest frequency mode. As for the time evolution scheme, to get higher modes it is only necessary to combine the Gauss-Seidel iteration with the Gram-Schmidt orthogonalization procedure described in Sec. IV B.

### D Boundary Conditions

To completely define the eigenvalue problem, boundary conditions must be specified. The simplest choices are Dirichlet and free boundary conditions. For Dirichlet boundary conditions, we imagine placing fictitious spin vectors $`𝐒_{𝐧_𝐛}`$, fixed to the static vortex directions, at sites $`𝐧_𝐛`$ just outside the system. The sites just inside the system interact with these, leading to on-site terms in the matrix $`M`$ due to the boundary, see Eq. (15b),

$`M_{\mathrm{𝐧𝐧}}^{\stackrel{~}{x}\stackrel{~}{y}}`$ $`=`$ $`M_{\mathrm{𝐧𝐧}}^{\stackrel{~}{y}\stackrel{~}{x}}`$ (66)
$`=`$ $`JS{\displaystyle \underset{𝐧_𝐛}{}}\left[p_𝐧^0p_{𝐧_𝐛}^0\mathrm{cos}(\varphi _𝐧^0\varphi _{𝐧_𝐛}^0)+\lambda m_𝐧^0m_{𝐧_𝐛}^0\right].`$ (67)
These are the only effects of the boundary.

For free boundary conditions, the system is open at the boundary, and no fictitious spins are needed. Therefore, the above terms of Eq. (66) are absent. This is the only difference in the matrix between Dirichlet and free boundary conditions.

However, in the absence of an in-plane magnetic field, free boundary conditions lead to a zero frequency mode, due to the in-plane rotational invariance of the model. This mode is unusual in that it cannot be made orthogonal to the other modes by the overlap integrals discussed above. Therefore, it needs to be removed from the evolving mode solution by a different tactic.

Consider its structure for $`\lambda <\lambda _c`$, i.e., in-plane vortices, where $`m_𝐧^0=0`$, $`p_𝐧^0=1`$ on all sites. Eqns. (32) become fairly simple,

$$JS\underset{𝐦=𝐧,𝐧^{}}{}\mathrm{cos}(\varphi _𝐧^0\varphi _𝐦^0)\left(w_𝐧^\mathit{2}w_𝐦^\mathit{2}\right)=i\omega _kw_𝐧^\mathit{1},$$
(69)
$$JS\underset{𝐦=𝐧,𝐧^{}}{}\left\{\mathrm{cos}(\varphi _𝐧^0\varphi _𝐦^0)w_𝐧^\mathit{1}\lambda w_𝐦^\mathit{1}\right\}=i\omega _kw_𝐧^\mathit{2}.$$
(70)
For the zero frequency mode, there is the solution,

$$\psi _{0,𝐧}=\left(\begin{array}{c}w_{k,𝐧}^\mathit{1}\\ w_{k,𝐧}^\mathit{2}\end{array}\right)=\left(\begin{array}{c}0\\ \phi _0\end{array}\right),$$
(71)
where $`\phi _0`$ is an arbitrary constant. Obviously this mode is not normalizable by (47), and represents the fact that a uniform in-plane rotation costs no energy.

For $`\lambda >\lambda _c`$, i.e., out-of-plane vortices, there is no great simplification of (32), however, a short calculation shows that the zero frequency mode is

$$\psi _{0,𝐧}=\left(\begin{array}{c}w_{k,𝐧}^\mathit{1}\\ w_{k,𝐧}^\mathit{2}\end{array}\right)=\left(\begin{array}{c}0\\ \phi _0\mathrm{cos}\theta _𝐧^0\end{array}\right),$$
(72)
where $`\phi _0`$ again is an arbitrary constant, and $`\theta _𝐧^0`$ is the static out-of-plane structure of the vortex. This wavefunction is slightly diminished in amplitude near the vortex core, where the spins point out of the easy plane. One can think of the $`w^\mathit{2}`$ component as a “coordinate” ($`w_𝐧^\mathit{2}=\phi _𝐧\mathrm{cos}\theta _𝐧^0`$ in-plane angle), which has the specified value, and $`w^\mathit{1}`$ as the corresponding conjugate momentum (part of $`S^z`$ due to the spinwave), which has been set to zero.

In the evolution to find other modes, any component of the current solution proportional to $`\psi _0`$ must be removed. For $`\lambda <\lambda _c`$, it is fairly obvious that we only need to enforce the constraints, $`w_𝐧^\mathit{1}=0`$, $`w_𝐧^\mathit{2}=0`$, where $``$ indicates the spatial average. These are essentially the conditions that the total conjugate momentum is zero, and that modes have no uniform in-plane drift component.

For $`\lambda >\lambda _c`$, physically, the higher modes must also have a total conjugate momentum equal to zero, and, clearly, must satisfy $`\varphi _𝐧=0`$. This leads to the general constraints we applied for free boundary conditions,

$$\delta S_𝐧^z=w_𝐧^\mathit{1}\mathrm{cos}\theta _𝐧^0=0,$$
(74)
$$\phi _𝐧=w_𝐧^\mathit{2}/\mathrm{cos}\theta _𝐧^0=0,$$
(75)
where we used $`S_𝐧^z=S_𝐧^{0z}+\delta S_𝐧^z=\mathrm{sin}\theta _𝐧^0+w_𝐧^\mathit{1}\mathrm{cos}\theta _𝐧^0`$ to define the dynamical conjugate momentum, $`\delta S_𝐧^z`$.

### E Memory and CPU Considerations

These relaxation methods use considerably less memory than the diagonalization scheme, provided that we are satisfied to obtain only a few of the lowest frequency modes, which is the case here. In order to obtain good stability for systems up to radius $`R=100a`$, it is necessary to use double precision.

If we want to calculate the lowest $`K`$ modes of a system $`N`$ sites, then we need to save $`4N`$ complex numbers for each state to be calculated (creation and annihilation operators), or $`8NK`$ doubles for all states. The matrix elements of $``$ are real and occupy $`(4\times 4+2)\times N=18N`$ memory locations, where the factor accounts for not storing the elements that are zero. For either type of relaxation scheme, it is convenient to have two complex arrays each with $`2N`$ elements to hold $`\mathrm{\Psi }`$ and $`^2\mathrm{\Psi }`$ ($`8N`$ doubles). For the Runge-Kutta time evolution, another two complex arrays of size $`2N`$ are needed. For the Gauss-Seidel iteration, on the other hand, the real terms $`\frac{f}{_{𝐧,𝐧}^{\alpha \alpha }}`$ are saved, requiring $`N`$ memory locations. For the asynchronous Gauss-Seidel, cpu time is reduced by saving $`_{𝐧,𝐧^{\prime \prime }}`$, which requires $`4\times 12\times N`$ doubles (there is a total of 12 nearest and next-nearest neighbors on the square lattice). Thus the total program size required for a calculation at a single value of anisotropy $`\lambda `$ scales as $`(34+8K)N`$ for Runge-Kutta and as $`(27+8K)N`$ (synchronous) or as $`(75+8K)N`$ (asynchronous) for Gauss-Seidel.

For a circular system of radius $`R=20a`$, with $`N=1264`$, and choosing $`K=15`$, these come to 1.52 MB for Runge-Kutta and 1.45 MB or 1.93 MB for Gauss-Seidel, considerably smaller than the approximately 100 MB necessary by diagonalization. The required memory increases if we want to track the changes in the modes with changing anisotropy, which requires storage of the set of modes at two different values of $`\lambda `$. The primary drawback of these methods is that the cpu time may be quite long, because the solution method is essentially diffusive. We find that the cpu time increases as $`R^4`$, where one factor of $`R^2`$ is the increase in system area, and the other factor of $`R^2`$ is due to the diffusion time increasing, $`t_{\mathrm{diffuse}}R^2`$. Often we are interested to calculate for a set of 20 to 40 values of $`\lambda `$ between zero and one; the large cpu time has prevented us from going beyond system radius $`R=100a`$, where each mode may take on the order of a day on typical 100 MHz processors.

It is also clear that the rate at which each mode converges diminishes with increasing system size, because the frequency differences between modes decreases as $`1/R`$. Furthermore, it is especially difficult to converge to any mode that is barely split from the next higher mode. However, the doubly degenerate modes present in this magnetic spinwave problem do not cause difficulties. Any arbitrary linear combination of such a pair will solve Eq. (36), and then the Gram-Schmidt process during the relaxation to the next higher mode will always produce another equal frequency mode, orthogonal to the first, spanning the degenerate subspace.

## V Relation to Continuum Theories

For a full theoretical description of vortex-spinwave interactions we need to mention how to compare our results for the wavefunctions with the spin fields determined in continuum limit theories. This is straightforward for the FM model, for which we will only make a few comments. For the AFM model, however, we are usually analyzing the staggered magnetization $`\stackrel{}{\mathrm{}}`$. Below we describe how $`\stackrel{}{\mathrm{}}`$ is related to the $`w_𝐧^\alpha `$ fields.

### A FM Continuum Coordinates

For the calculations on a lattice described so far, it was most convenient to use the Cartesian spin components. To make some comparison to theory we used planar spherical coordinates \[Eq. (7)\], with $`\theta _{\mathrm{planar}}`$ measured from the $`xy`$-plane, such that $`\theta _{\mathrm{planar}}`$ is conveniently zero in the ground state. Of course, in much of the literature it is common to use polar spherical coordinates, with $`\theta _{\mathrm{polar}}`$ measured from the positive $`z`$-axis. Clearly, these are related by $`\theta _{\mathrm{polar}}=\pi /2\theta _{\mathrm{planar}}`$. Then the relation that replaces (13), between our local Cartesian components and the angular fluctuations, using polar spherical coordinates, is

$$S_𝐧^{\stackrel{~}{x}}=S\phi _𝐧\mathrm{sin}\theta _𝐧^0,S_𝐧^{\stackrel{~}{y}}=S\vartheta _𝐧.$$
(76)
The reason we mention this is that the quantity $`\mu =\phi (𝐫)\mathrm{sin}\theta ^0(𝐫)`$ appears naturally in the continuum theory for FM vortex-spinwave scattering, and it is equivalent as well, via Eq. (28), to our $`w_𝐧^\mathit{2}`$ field. Similarly, our $`w_𝐧^\mathit{1}`$ field represents the $`\vartheta `$-fluctuations discussed in continuum theory.

### B AFM Continuum Coordinates

For the AFM scattering problem, there have been two primary ways to approach the continuum description, where we consider only lattices composed of two sublattices.

For example, in Refs. , the theory is developed directly in terms of $`\stackrel{}{\mathrm{}}`$, defined by the difference of spins on the two sublattices,

$$\stackrel{}{\mathrm{}}=\frac{𝐒_𝐧𝐒_𝐧^{}}{2S}=\stackrel{~}{z}_𝐧+\frac{S_𝐧^{\stackrel{~}{x}}+S_𝐧^{}^{\stackrel{~}{x}}}{2S}\stackrel{~}{x}_𝐧+\frac{S_𝐧^{\stackrel{~}{y}}S_𝐧^{}^{\stackrel{~}{y}}}{2S}\stackrel{~}{y}_𝐧,$$
(77)
where the local axes are related by $`\stackrel{~}{z}_𝐧^{}=\stackrel{~}{z}_𝐧`$, $`\stackrel{~}{x}_𝐧^{}=\stackrel{~}{x}_𝐧`$, $`\stackrel{~}{y}_𝐧^{}=\stackrel{~}{y}_𝐧`$, and $`𝐧,𝐧^{}`$ are neighboring sites, on the two different sublattices. Assuming $`\stackrel{}{\mathrm{}}`$ has perturbations and is described by planar spherical angles $`\varphi _{\mathrm{}}=\varphi _{\mathrm{}}^0+\phi _{\mathrm{}}`$, $`\theta _{\mathrm{}}=\theta _{\mathrm{}}^0+\vartheta _{\mathrm{}}`$, and using Eq. (28), we then obtain the equivalences,

$$\mathrm{}^{\stackrel{~}{x}}=\frac{S_𝐧^{\stackrel{~}{x}}+S_𝐧^{}^{\stackrel{~}{x}}}{2S}=\frac{w_𝐧^\mathit{2}+w_𝐧^{}^\mathit{2}}{2S}=\phi _{\mathrm{}}\mathrm{cos}\theta _{\mathrm{}}^0,$$
(79)
$$\mathrm{}^{\stackrel{~}{y}}=\frac{S_𝐧^{\stackrel{~}{y}}S_𝐧^{}^{\stackrel{~}{y}}}{2S}=\frac{w_𝐧^\mathit{1}+w_𝐧^{}^\mathit{1}}{2S}=\vartheta _{\mathrm{}}.$$
(80)
Here the first equation gives the in-plane fluctuations while the second gives the out-of-plane tilting fluctuations. The mapping to the $`w_𝐧^\alpha `$ fields is correct up to an arbitrary phase. If polar spherical coordinates are used, then change $`\mathrm{cos}\theta _{\mathrm{}}^0\mathrm{sin}\theta _{\mathrm{}}^0`$ and $`\vartheta _{\mathrm{}}\vartheta _{\mathrm{}}`$. Once again, it is found that $`\mathrm{}^{\stackrel{~}{x}}`$ is equivalent to the variable $`\mu _{\mathrm{}}`$ introduced by Ivanov et al, that appears together with $`\vartheta _{\mathrm{}}`$ in symmetrical Schrodinger-like equations.

Alternatively, in other vortex-magnon scattering articles, the spins on the two sublattices have been represented in polar spherical coordinates introduced by Mikeska, using four angular fields, as

$$𝐒=\pm S\left(\begin{array}{c}\mathrm{sin}(\mathrm{\Theta }\pm \theta )\mathrm{cos}(\mathrm{\Phi }\pm \varphi )\\ \mathrm{sin}(\mathrm{\Theta }\pm \theta )\mathrm{sin}(\mathrm{\Phi }\pm \varphi )\\ \mathrm{cos}(\mathrm{\Theta }\pm \theta )\end{array}\right),$$
(81)
where for simplicity we suppressed the site indices, and upper/lower signs refer to the two sublattices. The “small” angles $`\theta `$ and $`\varphi `$ are slave to the “large” angles $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$, via $`\theta =\dot{\mathrm{\Phi }}\mathrm{sin}\mathrm{\Theta }/8JS`$, and $`\varphi =\dot{\mathrm{\Theta }}/(8JS\mathrm{sin}\mathrm{\Theta })`$, showing that both are truly small for low-frequency modes. Then assuming that the large angles have small perturbations away from a vortex configuration $`\mathrm{\Theta }^0,\mathrm{\Phi }^0`$, $`\mathrm{\Theta }=\mathrm{\Theta }^0+\vartheta `$, $`\mathrm{\Phi }=\mathrm{\Phi }^0+\phi `$, a short calculation show that the staggered magnetization will have components, $`\mathrm{}^{\stackrel{~}{x}}=\phi \mathrm{sin}\mathrm{\Theta }^0`$, $`\mathrm{}^{\stackrel{~}{y}}=\vartheta `$. We see that these are equivalent to the perturbation variables in Eq. (77), once the change from polar to planar coordinates is taken into account.

## VI Application: Vortex-Spinwave Spectra and S-Matrix

The methods described above were tested by comparing their results with numerical diagonalization results for small systems (up to radius $`R20a`$). At these small system sizes, it was even possible to get very precise results (more than 6 digits for eigenvalues) for as many as the lowest 50 modes. Typically, for $`\lambda 0`$, the mode frequencies for the AFM model are higher than those for the FM model, which means that the calculations tend to converge faster for the AFM model.

As an example of the utility of these methods, we consider the spinwave spectrum (lowest 30 modes) for a single vortex in a circular system, with Dirichlet boundary conditions, and use the results to calculate the scattering S-matrix. We used a square lattice of spins. For most of the data presented, the synchronous self-consistent Gauss-Seidel scheme with fraction $`f=0.7`$ was used. The FM model for $`\lambda `$ very close to the critical value is an exception; there we used the asynchronous sheme, with $`f=0.9`$, which has greater stability at the expense of greater CPU time.

Typical mode spectra for some small systems were given in Ref. , as functions of the anisotropy parameter $`\lambda `$. In the present work we can consider system radii as large as $`R=100a`$, however, the CPU time to get the complete dependence versus $`\lambda `$ at such a large radius is prohibitive, and for most calculations we used $`R50a`$. We concentrate on the spectra at a few isolated values of $`\lambda <\lambda _c`$, corresponding to in-plane vortices, and focus more on the variations with system radius $`R`$.

### A FM/AFM In-plane vortex S-Matrix

Due to azimuthal symmetry, the modes generally are classified according to azimuthal quantum number $`m`$, corresponding to an expected $`e^{im\chi }`$ dependence, where a point in the system has polar coordinates $`(r,\chi )`$, measured from the center. Similarly, a principal quantum number $`n`$ can be defined, that is the number of nodes in a radial direction in the wavefunction. There are only minor differences in calculating a scattering amplitude $`\rho _m(k)`$ \[defined below, see Eq. (VI A)\] for the FM model compared to the AFM model. For FM out-of-plane vortices, the rotational symmetry around the $`S^z`$ axis is broken by the nonzero topological charge and nonzero $`S^z`$ spin components, and generally, $`m`$ and $`m`$ modes are not degenerate. For the in-plane case we consider here, however, the topological charge is zero, rotational symmetry holds, and $`\pm m`$ modes are degenerate. In the continuum limit, the zero-field AFM model has complete rotational and inversion symmetry around the $`S^z`$ axis, with no topological charge, so that the modes $`\pm m`$ are degenerate at any value of $`\lambda `$. For both FM and AFM in-plane vortices, the lack of $`S^z`$ spin components insures degeneracy of $`\pm m`$ modes. The degeneracy is partially broken on a lattice; for even $`m`$, the eigenfunctions generally are mixtures of $`+m`$ and $`m`$ components.

We can analyze the scattering S-matrix both for FM and AFM models as follows. In the continuum description, the lowest modes, which belong to the acoustic spinwave branch, have the asymptotic form,

$$w^\mathit{1}(r,\chi )0,$$
(83)
$$w^\mathit{2}(r,\chi )\left[J_m(kr)+\rho _m(k)Y_m(kr)\right]e^{im\chi },$$
(84)
where $`J_m`$ and $`Y_m`$ are Bessel and Neumann functions, and $`k`$ is the wavevector as determined from the frequency of the mode. For the AFM model, the free spinwave dispersion relation is

$$\omega _𝐤=4JS\sqrt{(1\gamma _𝐤)(1\pm \lambda \gamma _𝐤)},$$
(85)
where upper/lower signs are for the acoustic/optical branch, and

$$\gamma _𝐤=\frac{1}{2}(\mathrm{cos}k_x+\mathrm{cos}k_y).$$
(86)
For the FM model, we have only an acoustic branch with dispersion relation

$$\omega _𝐤=4JS\sqrt{(1\gamma _𝐤)(1\lambda \gamma _𝐤)}.$$
(87)
Note that the AFM acoustic branch dispersion can be obtained from the FM dispersion simply by the change $`\lambda \lambda `$. After the relaxation scheme is used to obtain a mode frequency $`\omega `$, then we invert the dispersion relation to get the corresponding wavevector needed in the fitting expression (84). For the FM dispersion, assuming $`𝐤=(k,0)`$, we have $`k=\mathrm{cos}^1(2\gamma _𝐤1)`$, where

$$\gamma _𝐤=(2\lambda )^1\left\{(1+\lambda )\sqrt{(1+\lambda )^24\lambda [1(\omega /4JS)^2]}\right\}.$$
(88)
The same formula applies to the AFM acoustic branch when the change $`\lambda \lambda `$ is made.

In Eq. (VI A), the coefficient $`\rho _m(k)`$ is a measure of the scattering; in the absence of scattering it vanishes, and the free spinwave is described by the $`J_m(kr)`$ function alone. Alternatively, the $`w^\mathit{2}`$ wavefunction can be expressed using the S-matrix scattering function, $`S_m(k)=e^{2i\mathrm{\Delta }_m(k)}`$, where $`\mathrm{\Delta }_m(k)`$ is the phase shift, in terms of incoming and outgoing waves, as

$$w^\mathit{2}\sqrt{\frac{1}{2\pi kr}}\left\{e^{i(kr+\varphi _m)}+S_m(k)e^{i(kr+\varphi _m)}\right\},$$
(89)
where the phase angle is $`\varphi _m=\frac{\pi }{4}(2m+1)`$. In the absence of scattering, $`S_m(k)1`$, or $`\mathrm{\Delta }_m(k)0`$, and the asymptotic form of $`J_m(kr)`$ results. Expressions (VI A) and (89) are equivalent up to a constant, provided

$$S_m(k)=\frac{1i\rho _m(k)}{1+i\rho _m(k)}.$$
(90)
Then the phase shift is related to $`\rho _m(k)`$ by

$$\mathrm{tan}\mathrm{\Delta }_m(k)=\rho _m(k).$$
(91)
The scattering amplitude $`\rho _m(k)`$ can be found by making a least squares fitting of a wavefunction to the above asymptotic form (84). However, in the case $`m>0`$, when $`m`$ and $`m`$ are degenerate, we actually fit to the expression,

$`w^\mathit{2}(r,\chi )`$ $`=`$ $`\left[a_1J_m(kr)+b_1Y_m(kr)\right]e^{im\chi }`$ (92)
$`+`$ $`\left[a_2J_m(kr)+b_2Y_m(kr)\right]e^{im\chi },`$ (93)
where $`a_1`$, $`a_2`$, $`b_1`$, and $`b_2`$ are complex fitting constants. Then, $`\rho _m(k)=b_1/a_1`$ and $`\rho _m(k)=b_2/a_2`$ are extracted, and should be expected to be comparable. For in-plane vortices, whose static structure has no $`\lambda `$-dependence, we considered the asymptotic region to be $`r>8a`$.

### B The Scattering Results: In-Plane Vortices

We considered in-plane FM and AFM vortices for $`\lambda =0.0,0.5,0.7`$. The static in-plane vortex spin structure is confined to the xy plane, resulting in a singular point at the vortex core, and making continuum theory difficult. For any $`\lambda <\lambda _c`$, the static structure of in-plane FM vortices is the same as the in-plane AFM vortex structure on one sublattice. However, the dynamics are different, in particular, the FM mode frequencies tend to be much lower than those for the AFM, although the appearances of their wavefunctions are very similar. The only exception to this is at $`\lambda =0`$, where the FM and AFM models have the same mode frequencies. Thus it is interesting to compare the FM/AFM scattering properties and also investigate the anisotropy dependence.

Wavefunctions of some of the lowest modes were presented in Refs. , and more complete wavefunction and other data, from the new methods described here, can be found at http://www.phys.ksu.edu/~wysin/vortexmodes/ . From those diagrams it is easy to identify the principal and azimuthal quantum numbers $`n`$ and $`m`$ for each mode. In Fig. 1 we show the dependence of some of the lowest eigenfrequencies on system radius $`R`$, for $`\lambda =0`$ (identical results for FM and AFM here). Differences appear in these spectra as $`\lambda `$ approaches the critical value, as seen in Fig. 2, where the spectra are shown for $`\lambda =0.70`$, which is just below $`\lambda _c0.7034`$ . Additional data for $`\lambda =0.5`$ are presented at the www page listed above. With the exception of the local mode in the AFM, most modes have frequencies diminishing asymptotically as $`1/R`$, as expected in a continuum theory. The local mode in the AFM model has a frequency independent of the system size, but dependent on $`\lambda `$, and can be thought of as a state taken out of the optical spinwave branch. The lowest mode in the FM model is only quasi-local: its wavefunction is local in the out-of-plane fluctuations but extended in the in-plane fluctuations, and its frequency decreases slowly with increasing system size.

We calculated the spectrum of the lowest 30 modes for different system sizes $`R`$, ranging from $`R=13a`$ to $`R=100a`$. Inverting the $`\omega _k`$ relationship (87) as described above \[via Eq. (88)\], and by using these different systems sizes, we could observe a particular mode labeled by $`(n,m)`$ at different $`k`$. In this way we obtained the $`k`$-dependence of the scattering.

#### 1 Scattering by FM Vortices

Some results obtained for the scattering of spinwaves by FM vortices are shown in Fig. 3 and Fig. 4, where the scattering amplitude $`\rho _m(k)`$ is displayed for $`\lambda =0.0,0.7`$. The scattering tends to be much weaker as $`m`$ increases above 2. As $`k0`$, $`\rho _m(k)`$ tends to zero, a physically reasonable limit. An interesting result, especially as $`\lambda \lambda _c`$, is that $`\rho _m(k)`$ has points where it is singular and changes sign. Indeed, the $`m=0`$ scattering curve at $`\lambda =0.7`$ exhibits two singularities, one near $`ka=0.035`$ and another near $`ka=0.65`$. Note that the scattering results obtained for $`\lambda =0.0,0.5`$, probe down to $`ka=0.024`$; down to this point they do not display the lower singularity.

Using Eq. (91), we see that these singular points in $`\rho _m(k)`$ correspond to the phaseshift $`\mathrm{\Delta }_m(k)`$ passing through the value $`\pi /2`$ near $`ka=0.035`$ and near $`ka=0.65`$. In Figs. 5 and 6 we show the resulting phaseshifts $`\mathrm{\Delta }_m(k)`$ for $`\lambda =0.0,0.7`$; the phaseshifts are smooth functions of $`k`$. It is impressive that the $`m=0`$ phase shift for $`\lambda =0.7`$ rapidly surpasses $`\pi /2`$, reaches a maximum and then falls back through $`\pi /2`$. Additional data for other parameters is found at the www page. Note that as long as $`\rho _m(k)1`$ there is little difference between it and the phase shift, except for the reverse sign.

Costa et al. calculated the $`m=0`$ and $`m=1`$ phase shifts for the FM model, only for $`\lambda =0`$, in the continuum limit via a Born approximation. Our results here are somewhat in contrast to theirs; from Fig. 3 and Eq. (91) one sees that we obtain the opposite sign, $`\mathrm{\Delta }_m(k)>0`$ for $`m=0,1,2`$, and over the range, $`0ka<0.5`$, there are no singularities in $`\rho _m(k)`$. Thus these phase shifts do not reach $`\pi /2`$ in this interval. For $`m=0`$, the Costa et al. result is similar to ours (except for the opposite sign), whereas for $`m=1`$ their result for the phase shift $`\mathrm{\Delta }_1(k)`$ has reached $`\pi `$ at $`ka0.3`$. Note, however, that the analytic calculations require the (somewhat arbitrary) choice of a short distance cutoff on the integrals appearing in the Born approximation, at a distance on the order of a lattice constant ($`a/2`$ was used in Ref. ). It is clear that a different cutoff could lead to drastically modified results, as was shown, for example, in a calculation of the frequencies of the local mode of AFM in-plane vortices. The choice of how to smoothly convert the singular field in the core of the vortex into a smooth continuum field is not clear, and always makes it difficult to write out an unambiguous continuum theory for the in-plane vortices. Therefore it may not be so fruitful to make any detailed comparison of the discrete lattice results here to continuum results for in-plane vortices, unless no cutoff is used. (This is not a problem for the out-of-plane vortices, which are not singular at their core.)

#### 2 Scattering by AFM Vortices

At $`\lambda =0`$, there is no difference in the scattering for AFM vortices compared to FM vortices; Figs. 3 and 5 apply. (Only the acoustic AFM branch is being considered; the optical branch is too high in the spectrum for these relaxation techniques to be useful.) Similarly, results for the scattering by AFM vortices at $`\lambda =0.7`$ are shown in Fig. 7. Results obtained at $`\lambda =0.5`$ are very similar to these, and shown on the web page indicated above. For the $`m=1`$ scattering, there is a singular point somewhere around $`ka0.55`$; thus $`\mathrm{\Delta }_1(k)`$ passes through $`\pi /2`$ at this point, as seen in Fig. 8. Although the AFM frequencies are higher than the corresponding frequencies for the FM model, the data obtained probe to similar low values of $`ka0.024`$, which is essentially determined by the largest system size used, $`R=100a`$; no singularity for $`m=0`$ was found down to this wavevector limit. Therefore we can say that the $`m=0`$ scattering by AFM in-plane vortices, just below $`\lambda _c`$, is considerably weaker than that for FM in-plane vortices. On the other hand, the $`m=1`$ scattering in this case is much stronger by the AFM in-plane vortices than by the FM in-plane vortices.

Continuum results for the AFM acoustic branch spinwave scattering by in-plane vortices are intriguing. Pereira et al. have calculated the scattering phase shifts for the optical branch, which excite out-of-plane spin fluctuations (our $`w^\mathit{1}`$ field), even without the need for a Born approximation. However, in the leading order continuum theory, the potential due to the vortex appears only in the perturbation equation for the out-of-plane fluctuations (angle $`\vartheta _{\mathrm{}}`$); it is absent from the equation for the fluctuations of the in-plane components (angle $`\phi _{\mathrm{}}`$ or our $`w^\mathit{2}`$ field). Indeed, the dynamical equation for $`\phi _{\mathrm{}}`$ is simply a wave equation for free acoustic spinwaves. It means that the usual continuum limit predicts no scattering of acoustic spinwaves (which excite in-plane fluctuations) by in-plane AFM vortices. We might comment, however, that for out-of-plane vortices, the potential due to the vortex appears in both the equations, one for $`\vartheta _{\mathrm{}}`$, the other for the quantity $`\mu _{\mathrm{}}=\phi _{\mathrm{}}\mathrm{cos}\theta _{\mathrm{}}^0`$ \[Eq. (77)\]. So the scattering from in-plane vortices is a special case.

Of course, with our numerical calculations on the lattice for AFM in-plane vortices, we did obtain some relatively weak scattering for $`m=0,2`$, and rather stronger scattering for $`m=1`$. It is possible that this unexpected scattering is caused primarily by the vortex core, where the static spin directions rotate by $`90^o`$ for nearest neighbor sites. This gives large gradients in the core region, which cannot be described in the usual lowest order continuum limits. Some theoretical analysis of this kind of discrete effect, which always occurs around the vortex core, is needed.

## VII Vortex Mass

As another application, information about certain normal modes can be used to estimate a vortex mass. The results here apply only to in-plane vortices. The modes with $`(n=0,m=1)`$, which are doubly degenerate for in-plane vortices, are most strongly associated with the translational motion of the vortex center. If only these modes are added to the original static vortex structure, the resulting configuration produces either a circular or linear simple harmonic oscillatory motion of the vortex center. These motions are the same as those predicted from a simple Newtonian dynamical equation of motion for the vortex center, involving a force $`\stackrel{}{F}`$ and effective mass $`M`$, i.e.,

$$\stackrel{}{F}=M\dot{\stackrel{}{V}},$$
(94)
where $`\stackrel{}{V}`$ is the velocity of the vortex center $`\stackrel{}{X}(t)`$. The translational modes can be considered to be driven by the force caused by the system boundary and the lattice itself, which can be evaluated independently of the spinwave spectrum. In the simplest approximation, we assume a linear restoring force, $`\stackrel{}{F}=K\stackrel{}{X}`$, where the spring constant $`K`$ has been found by a relaxation scheme applied to a static vortex in Ref. . For in-plane vortices, it was found that $`K4JS^2/a^2`$, independent of the value of $`\lambda `$. Then, for this simple harmonic motion, the mass is found as

$$M=K/\omega _{0,1}^2,$$
(95)
where $`\omega _{0,1}`$ is the translation mode frequency. The results of this calculation are shown in Fig. 9, for both the FM and AFM models at various anisotropies. The fact that the translational mode frequencies diminish as the reciprocal system size is reflected there; the effective mass increases as the square of the system radius $`R`$. The masses are identical for FM and AFM vortices at $`\lambda =0`$, where their spinwave spectra are identical. The FM vortex mass increases with increasing $`\lambda `$, while the AFM vortex mass decreases. Thus we have the result that in some real dynamical sense AFM vortices will be easier to move than FM vortices. The fact that the mass depends on the system size should not be unexpected, because a vortex is not a localized object. Its spin configuration extends to the limits of the system, causing its inertia against motion inside the system to depend on the system size. The mass defined here really is more a property of the entire system.

## VIII Conclusions

We have discussed some relaxation methods for determining the lower part of the eigenspectrum for magnetic models in the presence of an individual magnetic nonlinear excitation, such as a vortex. The methods are based on the fact that the effective time evolution of the squared Hamiltonian is diffusive, and therefore evolution from an arbitrary intial state will lead to the lowest mode present. By combining with Gram-Schmidt orthogonalization, a reasonably large set of the lowest modes, including their eigenfunctions and eigenfrequencies, can be determined. Special considerations also were described for applying free boundary conditions, in contrast to Dirichlet boundary conditions; for the former the zero-frequency mode must be eliminated from the evolving solution.

As an application of the schemes, the vortex-spinwave scattering and S-matrix were determined numerically for both FM and AFM in-plane vortices on square lattice circular systems. The results obtained exhibit some interesting singular points in the scattering amplitude $`\rho _m(k)`$, especially as $`\lambda `$ approaches the critical value from below. A singularity in the $`m=1`$ scattering amplitude by FM or AFM vortices ($`\lambda =0.0`$) occurs near $`ka0.65`$. A singularity in the $`m=0`$ scattering amplitude for the FM at $`\lambda =0.7`$ occurs at very low wavevector ($`ka0.035`$) or over long distances, and it may be interesting to ask whether it could be described in a continuum theory. These singular points occur where the phase shift passes smoothly through values like $`\pm \pi /2`$.

At $`\lambda =0`$, the scattering results for FM and AFM vortices are identical, as a result of identical mode frequencies. For the AFM model, although theoretically there should be no scattering of acoustic spinwaves by a vortex, we did find nonzero scattering amplitude, especially for $`m=1`$, and even a singular point in this amplitude for $`\lambda =0.7`$ . This unexpected strong scattering is probably due to the large gradients in the spin field around the vortex core, which cannot accurately be described in the usual continuum limit theory.

The lowest vortex translational modes (lowest mode most closely coupled to motion of the vortex center) were used to estimate a vortex mass, for in-plane vortices, in a very simple phenomenolgical theory. The theory is based on the idea that the vortex in a small system vibrates in its translational mode in response to the force on it due to the lattice and system boundary, with a frequency determined by the effective force constant and the vortex mass. Using the observed translational mode frequencies, this simple theory leads to a vortex mass increasing as the square of the system size, for both FM and AFM vortices. AFM vortices have lower effective mass and thus are expected to be easier to move. For the in-plane vortices, FM vortices at $`\lambda \lambda _c`$ have the highest mass, while AFM vortices at $`\lambda \lambda _c`$ have the lowest mass.

With appropriate modifications to the Hamiltonian, etc., in general, the relaxation schemes discussed here would apply to other kinds of stable nonuniformly magnetized systems, such as, for example, small magnetic particles with a frozen-in magnetic configuration. Then the methods might be very useful for determining any interesting response properties of such materials.

Acknowledgements.—The author gratefully acknowledges the support of NSF Grant DMR-9412300, NSF/CNPq International Grant INT-9502781, and a FAPEMIG Grant for Visiting Researchers while at UFMG. Discussions with M. E. Gouvêa, A. S. T. Pires, and B. A. Ivanov are greatly appreciated.


format: text

text:
# PHENOMENOLOGY OF THE ELECTRON STRUCTURE FUNCTION22footnote 2Work supported by the Polish State Committee for Scientific Research (grant No. 2P03B06116).

## Abstract

Advantages of introducing the electron structure function (ESF) in electron induced processes are demonstrated. Contrary to the photon structure function it is directly measured in such processes. At present energies a simultaneous analysis of both the electron and the photon structure functions gives an important test of the experimentally applied methods. Estimates of the ESF at currently measured momenta are given. At very high momenta contributions from $`W`$ and $`Z`$ bosons together with $`\gamma `$-$`Z`$ interference can be observed. Predictions for next generation of experiments are given.

In a series of papers  we have presented the construction of the electron structure function — a useful notion in the QCD analysis of electron induced hadron production. The $`Q^2`$ evolution equations have been constructed and asymptotic solutions found for the quark and gluon content of the electron in the leading logarithmic approximation. We included contributions from all intermediate bosons, in particular we found that the $`\gamma `$-$`Z`$ interference is important at very high energies. This spoils the usual probabilistic interpretation of separate $`\gamma /Z/W`$ structure functions. We also found that in certain experimental situations the commonly used convolution of the photon structure function and the photon flux, used to describe the electron induced processes, is incorrect.

In this letter we study phenomenology of the electron structure function comparing it with the well known approach which makes use of the photon structure function. Theoretical framework which allows to calculate the photon structure is known since long . It appears as perturbative QCD contribution, in addition to the modelled Vector Meson Dominance term. To measure this photonic structure, experiments  use the electron (or positron) beam as a source of photons. Despite precise measurement the data are not easy to extract. The problem is displayed in Fig. 1a. The tagged (upper) electron emits a probing photon whereas the untagged (lower) one goes nearly along the beam, emitting the target photon. (The situation where at very high energies the probing boson can be also a $`Z`$ boson is considered below.) The large scale $`Q^2`$ is determined by the tagged electron:

$$Q^2=(kk^{})^2=2EE_{\mathrm{tag}}(1\mathrm{cos}\theta _{\mathrm{tag}}),$$
(1)
where $`E`$ is the initial electron energy and $`E_{\mathrm{tag}}`$ and $`\theta _{\mathrm{tag}}`$ are the energy and polar angle of the measured electron. The antitag condition (if present) requires the virtuality of the target photon $`P^2`$ to be less than certain $`P_{\mathrm{max}}^2`$:

$$(pp^{})^2P^2P_{\mathrm{max}}^2.$$
(2)
This photon is clearly not a beam particle and has the energy diffused according to the equivalent photon spectrum:

$$f_\gamma ^e(y_\gamma ,P^2)=\frac{\alpha }{2\pi P^2}\left[\frac{1+(1y_\gamma )^2}{y_\gamma }\frac{2y_\gamma m_e^2}{P^2}\right]$$
(3)
where $`y_\gamma `$ is the photon momentum fraction, $`\alpha `$ is the QED structure constant, and $`m_e`$ is the electron mass. The measured cross-section for the production of a hadronic system $`X`$, expressed in terms of the photon structure functions $`F_2^\gamma `$ and $`F_L^\gamma `$ reads:

$`{\displaystyle \frac{d^3\sigma _{eeeX}}{dzdQ^2dx}}`$ $`=`$ $`{\displaystyle \frac{2\pi \alpha ^2}{x^2Q^4}}{\displaystyle \underset{P_{\mathrm{min}}^2}{\overset{P_{\mathrm{max}}^2}{}}}[(1+(1y)^2)F_2^\gamma (x,Q^2,P^2)`$ (4)
$`y^2F_L^\gamma (x,Q^2,P^2)]f_\gamma ^e(y_\gamma ,P^2)dP^2`$
where

$`y=1(E_{\mathrm{tag}}/E)\mathrm{cos}^2(\theta _{\mathrm{tag}}/2),`$
$`P_{\mathrm{min}}^2=m_e^2y_\gamma ^2/(1y_\gamma )`$ (5)
and $`x`$ ($`z`$) are fractions of the parton momentum with respect to the photon (electron). They are related to the photon momentum fraction with respect to the electron by

$$z=xy_\gamma $$
(6)
The integral over the photon virtuality is usually performed by assuming $`P^2=0`$ in the photon structure functions (Weizsäcker-Williams approximation ) which leads to:

$`{\displaystyle \frac{d^3\sigma _{eeeX}}{dzdQ^2dx}}={\displaystyle \frac{2\pi \alpha ^2}{x^2Q^4}}`$
$`\times [(1+(1y)^2)F_2^\gamma (x,Q^2,0)y^2F_L^\gamma (x,Q^2,0)]`$
$`\times f_\gamma ^{\mathrm{WW}}(z/x,P_{\mathrm{max}}^2)`$ (7)
where

$`f_\gamma ^{\mathrm{WW}}(y_\gamma ,P_{\mathrm{max}}^2)`$ $`=`$ $`{\displaystyle \frac{\alpha }{2\pi }}[{\displaystyle \frac{1+(1y_\gamma )^2}{y_\gamma }}\mathrm{ln}{\displaystyle \frac{P_{\mathrm{max}}^2(1y_\gamma )}{m_e^2y_\gamma ^2}}`$ (8)
$`2{\displaystyle \frac{1y_\gamma }{y_\gamma }}]`$
This is how the real photon structure functions $`F_2`$, $`F_L`$ appear in lepton-lepton scattering.

Three remarks are important for further considerations. First, the splitting of the process into a distribution of photons inside electron $`f_\gamma ^{\mathrm{WW}}`$ and that of partons inside the photon $`F_2^\gamma `$ is an approximation. The optimal form of the equivalent photon formula is still being discussed . Even if most of experimental groups choose the same formula, one should keep in mind that the photon structure function depends on this convention.

Second, the target photon is always off-shell and its virtuality is experimentally not measured. Although the equivalent photon distribution is peaked at minimum (nearly zero) virtuality, treating the photon as real is another approximation. One should keep in mind that the measured photon structure function depends on $`x`$, $`Q^2`$ and $`P^2`$, and the, usually neglected, $`P^2`$ dependence can be quite strong . A hint that we are not measuring the real photon structure function comes also from the analysis of the QED structure function of the photon where analytical solutions (and thus the $`P^2`$ dependence) are known . The data (coming from $`e^+e^{}e^+e^{}\mu ^+\mu ^{}`$ scattering ) agree with theory only if non-zero virtuality is taken into account. In the case of QCD the situation is more difficult — experiments yield the photon structure function at virtuality $`P^2`$ which is not measured and theory does not predict this $`P^2`$ dependence.

Third, in order to fix $`x`$, one is forced to measure — in addition to the tagged electron — the hadronic momenta. In fact,

$$x=\frac{Q^2}{Q^2+P^2+W^2}\frac{Q^2}{Q^2+W^2},$$
(9)
where $`W`$ is the invariant mass of the produced hadronic system $`X`$. Its determination is more difficult than of other (tagged electron) variables since substantial part of hadrons is lost in the beam pipe. The uncertainty in the determination of the $`x`$ variable is the source of large uncertainties in the analysis (unfolding procedure). The data are indirectly biased by theoretical assumptions.

Many of the above problems can be avoided when we introduce the structure function of the electron (Fig. 1b). To see how it works let us first write the cross-section in terms of the electron structure functions $`F_2^e`$ and $`F_L^e`$:

$`{\displaystyle \frac{d^2\sigma _{eeeX}}{dzdQ^2}}`$ $`=`$ $`{\displaystyle \frac{2\pi \alpha ^2}{zQ^4}}[(1+(1y)^2)F_2^e(z,Q^2,P_{\mathrm{max}}^2)`$ (10)
$`y^2F_L^e(z,Q^2,P_{\mathrm{max}}^2)].`$
The structure function $`F_2^e(z,Q^2,P_{\mathrm{max}}^2)`$ which dominates the cross-section at small $`y`$, has simple partonic interpretation:

$$F_2^e(z,Q^2,P_{\mathrm{max}}^2)=z\underset{i}{}e_{q_i}^2q_i(z,Q^2,P_{\mathrm{max}}^2),$$
(11)
where $`e_{q_i}`$ and $`q_i`$ are the $`i`$-th quark fractional charge and density. This is a standard deep inelastic scattering process where the cross-section is related to the structure function via simple kinematical factors. More precisely the above defined electron structure function corresponds exactly to the proton structure function if no antitag condition is imposed and $`P_{\mathrm{max}}^2`$ goes up to its kinematical limit of the order of $`Q^2`$. (This is the ‘inclusive case’ in our terminology , the electron structure function depends then on $`z`$ and $`Q^2`$ only.)

The argument $`z`$ — the parton momentum fraction with respect to the electron — is measured, as in the standard deep inelastic scattering, by means of the tagged electron variables only:

$$z=\frac{Q^2}{2pq}=\frac{\mathrm{sin}^2(\theta _{\mathrm{tag}}/2)}{E/E_{\mathrm{tag}}\mathrm{cos}^2(\theta _{\mathrm{tag}}/2)}.$$
(12)
There is no need a priori to reconstruct the hadronic mass $`W`$. In present experimental analyses one introduces in both the photon and electron structure analyses lower limit on $`W`$ because the reconstruction of the hadronic mass is unreliable below this limit. But even so, there is an important difference between the case when $`W`$ (in small $`W`$ region) is used only as a selection cut to reduce background (electron structure) and the case when in addition it is used to determine the $`x`$ variable (photon structure). In particular a lower cut on $`W`$ is much less sensitive to the unfolding procedure than the value of $`W`$ itself. The quantitative analysis of $`W`$ influence is discussed below.

All these features cause that the same experiment can produce more precise and analysis independent data when looking at the electron structure. What is most important — the electron structure function contains the same information about QCD as the photon one and it is known theoretically with at least the same accuracy. Moreover, it allows to avoid problems which arise in the photon structure function at very high energies.

At present energies, where the $`W`$ and $`Z`$ bosons contributions are negligible, one can reanalyse the existing data in terms of the electron structure function. This can be treated also as a consistency check of both photon and electron structure. Phenomenologically, having a parametrisation of the photon structure function which describes well the existing data, we can predict the electron structure function for $`Q^2P_{\mathrm{max}}^2`$ by taking the convolution of this parametrisation with the equivalent photon spectrum:

$`F_2^e(z,Q^2,P_{\mathrm{max}}^2)`$ $`=`$ $`{\displaystyle \underset{z}{\overset{1}{}}}𝑑y_\gamma {\displaystyle \underset{P_{\mathrm{min}}^2}{\overset{P_{\mathrm{max}}^2}{}}}𝑑P^2f_\gamma ^e(y_\gamma ,P^2)`$ (13)
$`\times F_2^\gamma ({\displaystyle \frac{z}{y_\gamma }},Q^2,P^2)`$
Such convolution is correct when the experiments use antitag condition (‘exclusive case’ in our terminology ). The curves with momenta corresponding to LEP2 and TESLA/NLC/JLC experiments, resulting from some popular parametrisations  of real photon structure — i.e. $`P^2=0`$ in $`F_2^\gamma `$ of Eq. (13) — are shown in Fig. 2. To test how significant is the (neglected above) $`P^2`$ dependence we also plot the resulting curve of SaS-1D parametrisation with the $`P^2`$ dependence built in . In the latter case the integration over $`P^2`$ in Eq.( 13) is performed numerically. Our previous arguments are confirmed: already at LEP2 the photon virtuality is non-negligible, it can produce effects of the same order as differences between various parametrisations. The electron variables $`z,Q^2`$ and $`P_{\mathrm{max}}^2`$ are well known, the photon data do not contain one crucial information — the photon virtuality, at which the structure is measured.

As noticed in Ref. , the shape of the electron structure function is strongly influenced by the QED part (“photon flux”). All $`F_2^e`$ functions resulting from various photon structure parametrisations are falling functions of $`z`$ and look “similar”. In addition the contributions at given electron momentum fraction $`z=xy_\gamma `$ come both from large $`x`$ and small $`y_\gamma `$ as well as large $`y_\gamma `$ and small $`x`$. Both these features can be regarded as drawbacks of the electron structure function. The above arguments are only partly true. First, one should keep in mind that the variables $`x`$, $`y_\gamma `$ and $`z`$ are not independent and their interplay under the integral Eq. (13) is specific. To see the problem in more detail let us take the Eq. (13) within the Weizsäcker-Williams approximation ($`P_{\mathrm{max}}^2,P^2,Q^2`$ suppressed)

$$F_2^e(z)=\underset{z}{\overset{1}{}}𝑑y_\gamma f_\gamma ^{\mathrm{WW}}(y_\gamma )F_2^\gamma \left(\frac{z}{y_\gamma }\right)$$
(14)
Noting that essentially $`f_\gamma ^{\mathrm{WW}}(y_\gamma )1/y_\gamma `$ we get from Eq. (14):

$$F_2^e(z)\underset{z}{\overset{1}{}}\frac{dy_\gamma }{y_\gamma }F_2^\gamma \left(\frac{z}{y_\gamma }\right)=\underset{z}{\overset{1}{}}\frac{dx}{x}F_2^\gamma (x).$$
(15)
One sees that in this approximation the $`z`$ dependence comes via the lower limit of integration. Moreover the importance of the small $`x`$ region under the integral (14) is enhanced by the $`1/x`$ factor. Second, due to the same kinematics the data points are generally shifted towards lower $`z`$ (as compared to $`x`$). Therefore the experimental results are more accurate at small $`z`$, a feature mostly welcome in the region where new effects are to be expected. In addition, we remind that the photon data points at low $`x`$ are strongest influenced by the unfolding procedure.

As already mentioned, the present way of data analysis introduces a lower cut on the hadronic mass $`W`$. In the photon case it causes that we do not measure the photon structure function above certain $`x`$ (see Eq. (9)). In the measurement of the electron structure function $`W`$ is not needed at all (as it is not used in the analysis of the proton structure function). A cut on $`W`$ imposed in present experiments lowers the cross-section in the whole $`z`$ range, see Fig. 3. We checked that e.g. for $`Q^2=17.8\mathrm{GeV}^2`$ the effect of the condition $`W1.7`$ GeV (with SaS-1D parametrisation in Eq. (13)) is less than 5% for $`z0.01`$, and less than 9% for $`z0.1`$. This suppression gets smaller with growing $`Q^2`$ (e.g. at $`Q^2=120\mathrm{GeV}^2`$ it is below 1% for $`z0.1`$).

The concept of the electron structure function introduces new interesting effects at momenta much higher than presently available. One has then to take into account not only the photon flux contributing to the electron structure but also those resulting from $`Z`$ and $`W`$ bosons. As shown in Ref.  $`\gamma `$-$`Z`$ interference comes into play and is comparable to the $`Z`$ contribution itself. This means that a notion of separate gauge boson structure functions ($`\gamma ,Z`$ or $`W`$) looses sense and only the electron structure function preserves probabilistic interpretation. The question is, can we observe these effects in the next generation of experiments? In Fig. 4 we give quantitative estimate in the case of single tag $`e^+e^{}`$ scattering at CLIC  momenta, choosing $`Q^2=10000\mathrm{GeV}^2`$ and $`P_{\mathrm{max}}=1000\mathrm{GeV}^2`$. In this case the picture of Fig. 1 gets modified. The upper (tagged) electron emits now both the photon and the $`Z`$ boson. In the calculation of the electron structure one has to take into account contributions from the photon, $`Z`$ boson, their interference (with the antitag condition fixed by $`P_{\mathrm{max}}^2)`$ and the $`W`$ boson (no antitag condition). The presented curves are asymptotic solutions of our electron evolution equations . One sees that the effect of $`Z`$ and $`\gamma `$-$`Z`$ terms is of the order of $`515\%`$. We checked that it can be enhanced to $`2025\%`$ in a double tag experiment.

Let us add a few final remarks. First concerns the study of the virtual photon structure  (double tag experiment). The analysis can be reformulated in terms of the $`P_{\mathrm{max}}^2`$ dependence of the electron structure function. Studying a real, convention independent object is first advantage. Another one is the fact that at very high virtualities the $`Z`$ admixture and the $`\gamma `$-$`Z`$ interference are properly taken into account.

Second is a comment on the QED structure function of the photon. It is obtained from the process $`e^+e^{}e^+e^{}\mu ^+\mu ^{}`$ by dividing out the (approximate) equivalent photon distribution and assuming some effective photon virtuality. The use of the QED electron structure function avoids the above approximations. The exactly known (in given order of $`\alpha `$) electron structure function can be compared directly with the electron data.

Finally, the photon structure has been also measured  in dijet production at HERA. Again the extraction of the $`x`$ variable is difficult. In addition to jets, one has to measure essentially the whole hadronic system in order to obtain the photon energy. The data, when presented in terms of the electron structure, require only measurement of the two jets. Practically the new approach means plotting the dijet cross-section as a function of $`z`$. A more ambitious program would be to extract the parton densities in the electron, as it has been done in the case of the photon  or to construct a parametrisation of parton densities inside the electron by a direct fit to the HERA data.

To summarise, we propose to look at the electron as surrounded by a QCD cloud of quarks and gluons (in order $`\alpha ^2`$), very much like it is surrounded by a QED cloud of equivalent photons (in order $`\alpha `$). We argue that the use of the electron structure function in electron induced processes has some advantages over the photon one. Experimentally it leads to more precise, convention independent data. Theoretically it allows for more careful treatment of all variables. It also takes into account all electroweak gauge boson contributions, including their interference, which will be important in the next generation of $`e^+e^{}`$ colliders. At present energies it should certainly be used as a cross-check of the photon structure analysis.

The authors would like to thank Danuta Kisielewska, Maria Krawczyk, Aharon Levy, Bogdan Muryn, Mariusz Przybycień and Jacek Turnau for discussions.


text:
# Holomorphic Disks and the Chord Conjecture

## Arnold’s Chord Conjecture

In this paper we prove the following theorem which was conjectured by Arnold :

###### Theorem 1.

For every closed, compact Legendrian submanifold in $`S^{2n1}`$ with the standard contact structure and any contact form for this structure, there is a Reeb chord, i.e. an integral curve of the Reeb vector field which begins and ends on the Legendrian submanifold.

Theorem 1 will follow as a corollary from the main result of this paper, Theorem 4. In fact it can be applied to a more general situation:

###### Theorem 2.

Let $`(M,\xi )`$ be a contact structure which arises as smooth boundary of a compact subcritical Stein manifold (see  for a definition). Then for any Legendrian submanifold and any contact one form corresponding to $`\xi `$ there is a Reeb chord.

Our results include the existence of chords for Legendrians in the standard contact structure on $`P^{2n1}`$ proved by Ginzburg and Givental , although it does not provide their statement of linear growth. They cover results by Abbas  and Cieliebak  who treat subcases of the problem on the sphere and on boundaries of subcritical Stein manifolds.

## Lagrangian out of Legendrian embeddings

Consider a closed Legendrian submanifold $`lM^{2n1}`$ in a contact manifold. Given a contact one form $`\alpha `$ we will construct Lagrangian embeddings of the torus $`l\times S^1`$ into the symplectization $`(M\times ,d(e^s\alpha ))`$ of $`(M,\xi )`$ and study them. Denote by $`\varphi `$ the flow of the Reeb vector field $`R=R_\alpha `$.

Assume $`l`$ has no Reeb chords of length at most $`T>0`$. Then there is an embedding $`\mathrm{\Phi }:l\times (\times [0,T])M\times `$ such that $`\mathrm{\Phi }^{}(d(e^s\alpha ))=e^sdsdt`$, $`(s,t)\times [0,T]`$ being the coordinates of the infinite strip. The map is, of course, constructed using the Reeb flow: $`\mathrm{\Phi }(x,s,t):=(\varphi ^t(x),s)`$.

The upshot of all this is, that any smooth embedded loop in the rectangle $`[S,0]\times [0,T]`$ defines a smooth Lagrangian embedding of $`l\times S^1`$ into $`M\times [S,0]`$. Here $`S<0`$ is any negative parameter. The integral of the primitive $`e^s\alpha `$ over closed curves in such a Lagrangian depends only on the homology class of that curve. It vanishes on classes arising as closed loops in $`l`$, since $`\alpha |_l0`$. On $`\{x\}\times l`$ it is given by the area enclosed by the loop in the rectangle with respect to the volume form $`e^sdsdt`$. Approximating the boundary of $`[S,0]\times [0,T]`$ by smooth embedded loops we obtain the following

###### Proposition 3.

Assume the Legendrian $`lM`$ in a contact manifold has no Reeb chord of length smaller or equal to $`T`$ with respect to a contact one form $`\alpha `$. Then there exists a smooth Lagrangian embedding of $`l\times S^1`$ into the symplectic manifold $`(M\times [S,0],d(e^s\alpha ))`$ such that the integral of the primitive $`e^s\alpha `$ over a smooth closed curve in the image is an integer multiple of some constant $`C>(1e^S)T`$.

## Proof of Arnold’s Chord Conjecture

The results described above will follow from the following

###### Theorem 4.

Assume $`W`$ is a symplectic manifold of bounded geometry, i.e. there is an almost complex structure $`J`$ together with a hermitian metric $`g`$ whose sectional curvature is globally bounded from above and injectivity radius bounded from below, such that $`\omega `$ tames $`J`$ uniformly w.r.t. $`g`$: $`\omega (X,JX)\text{const.}X_g^2`$ for any tangent vector $`XTW`$. Suppose $`W`$ is monotone, i.e. any sphere in $`W`$ has vanishing symplectic area. Let $`(M\times [S,0],d(e^s\alpha ))(W,\omega )`$ be a symplectic embedding of a finite cylinder in the symplectization of a contact manifold into $`W`$, such that $`\pi _1(M)`$ injects into $`\pi _1(W)`$. Suppose the image of the embedding can be displaced by a (time-dependent) Hamiltonian flow with compact support and oscillation

$$H:=_0^1\left(\underset{xW}{\mathrm{max}}H(t,x)\underset{xW}{\mathrm{min}}H(t,x)\right)𝑑t.$$

Then any closed Legendrian $`lM`$ admits a Reeb chord of length not bigger than $`H/(1e^S)`$.

###### Proof.

Assume that there is a Legendrian embedding into $`M`$ which admits no Reeb chord of length not bigger than $`T>0`$. Consider a Lagrangian embedding as constructed in the previous section. It will be displaced from itself by $`H`$. Fix any almost complex structure $`J`$ as in the theorem. Due to monotonicity of $`W`$ there are no $`J`$–holomorphic spheres in $`W`$. Then an argument given by Chekanov in  produces a non–constant $`J`$–holomorphic disk in $`W`$ with boundary on $`L`$ and (symplectic) area smaller than or equal to $`H`$. But due to the injectivity of $`\pi _1(M)\pi _1(W)`$ there is another disk completely lying in $`M`$. Since $`W`$ is monotone the symplectic areas of both disks agree. The latter is an integer multiple of $`C>T(1e^S)`$. Non-constant holomorphic disks have positive symplectic area, hence $`TH/(1e^S)`$. ∎

Remarks. (1) Arnold’s conjecture follows from this result since every contact one form corresponding to the standard contact structure can be realized as the restriction of the primitive $`\omega (z,.)`$ of the standard symplectic structure in $`^n`$ to the smooth boundary of a star-shaped domain.
(2) Biran and Cieliebak observed in  that any compact set in a (complete) subcritical Stein manifold can be displaced from itself by a Hamiltonian isotopy with compact support. In particular, one can apply this to the strongly pseudo–convex domain in the Stein manifold within its completion to obtain Theorem 2.
(3) One could prove the result using Gromov’s Fredholm alternative for maps of the disk into $`W`$ with boundary on a Lagrangian satisfying a $`\overline{}`$–equation with a non-zero right-hand side (see  and ). It should be an amusing and (non-trivial) exercise to produce exactly the same estimate. For our convenience we chose to use Chekanov’s result but feel obliged to point out that the result is a consequence of Gromov’s original ideas.

Problem. We have not answered Arnold’s original question yet, whether there always exists a Reeb orbit intersecting the Legendrian submanifold at least twice. We proved that there is a chord whose length is estimated by $`S`$ and the displacement energy. But this chord could be a closed orbit which intersects the Legendrian once. Notice that this question is bound to the sphere. E.g. pick the contact form $`x_1d\theta +\frac{1}{2}(x_2dx_3x_3dx_2)`$ on $`S^1\times S^2`$ with coordinate $`\theta S^1`$ and $`S^2=\{x=(x_1,x_2,x_3)^3x=1\}`$. Then $`S^1\times \{0\}`$ is a Legendrian loop which intersects each of the (closed) Reeb orbits $`\{\theta \}\times \{x_1=0\}`$ only once.

Acknowledgements. I would like to thank Leonid Polterovich, Kai Cieliebak, Yasha Eliashberg, Tobias Ekholm and John Etnyre for their interest and comments, and the Department of Mathematics of Stanford for the warm hospitality during my visit.