text:
# Internal Modes and Magnon Scattering on Topological Solitons in 2d Easy–Axis Ferromagnets

## I Introduction

Nonlinear topologically nontrivial excitations (solitons) are well–known to play a special role in low–dimensional magnetic systems. For example, the presence of vortices in 2d easy–plane (EP) magnets gives rise to the Berezinskiĭ–Kosterlitz–Thouless (BKT) phase transition. Kinks in 1d magnets and localized Belavin–Polyakov (BP) solitons  in 2d isotropic magnets are responsible for the destruction of long–range order at finite temperature. The soliton signatures in dynamical response functions can be observed experimentally. Translational motion of solitons leads to the so–called soliton central peak. Another possibility to detect soliton signature is to look for magnon modes, localized at the soliton (local modes, LM); such modes are the cause of soliton magnetic resonance at the characteristic frequencies of “intrinsic” motion.

This can be explained within the scope of the so–called soliton phenomenology, where the magnet can be described as a two–component gas of elementary excitations: solitons and magnons. Such an approach was developed for 1d magnets. It was shown that the contribution of magnons and soliton–magnon interaction is important for this approach, as is obvious, for example, in the discussion of soliton magnetic resonance. But for 1d magnets the soliton–magnon scattering causes a change of the magnon density of states, which is necessary for a self–consistent calculation of the temperature dependence of the soliton density.

For the 2d case, the concept of a soliton (vortex) gas has been extended to describe EP magnets above the BKT vortex–antivortex unbinding transition, with a finite density of vortices. The concept of localized topological solitons ($`\pi _2`$–topological solitons, see details in Ref.) was used to explain the EPR line–width in easy–axis (EA) magnets. The classical example of a localized soliton is BP–solution, which exists in isotropic magnets. The scattering problem for such solitons allows an exact analytical solution because of the scale– and conformal– invariance of the model. For EA magnets it is necessary to take account of some important modifications due to the breaking of such symmetries. That is why, both for EA and EP 2d magnets, the soliton density has not been calculated, but has been used as an input parameter. The general features of the 2d soliton dynamics are not clear at present. In particular, the form of inertial terms in the dynamical equations for the soliton centre is unknown; so the study of localized states is an issue of current research. Note that the problem of the soliton dynamics and the problem of the existence of LM are intimately connected to each other, and to the problem of soliton–magnon scattering. For example, using numerical data for the scattering amplitude, a non–Newtonian effective equation of motion of the magnetic vortex was constructed.

Usually non–1d solitons and, especially, the problem of soliton–magnon scattering, are treated numerically. Only in the case of purely isotropic magnets, with the exact BP–solution, the analytical investigation of the scattering problem has been done recently. In particular, it was found that a soliton has a number of local modes of zero frequency; such modes determine the main features of the scattering picture. However, even a small uniaxial anisotropy leads to principal changes. Truly local modes (in contrast to quasi–local ones, below) with nonzero frequencies were found numerically in Ref.  for the soliton in EA ferromagnets (FM) for some special values of the azimuthal quantum number $`m`$; namely the mode with $`m=1`$ for the soliton with the topological charge $`q=1`$, and the mode with $`m=2`$ for $`q=2`$ were investigated. Note that local modes do not exist for the models of easy–plane and isotropic ferromagnets, where the magnon spectrum has a gapless dispersion law. The case of the 2d EP antiferromagnet is an exception; for this model a finite–frequency, truly localized internal mode with $`m=0`$ exists inside the continuum spectrum.

Recently there has also been renewed attention to the internal dynamics of the soliton. It is caused by the studying of small ferromagnetic particles in the so–called vortex state (magnetic dots). Magnetic dots are micron–sized magnetic samples, placed on a non–magnetic substrate; they have different shapes: circular, quadratic, etc. The magnetic dots are interesting for the usage in the high–density magnetic storage devices and from a fundamental point of view. As the dots are rather large (the size of magnetic dot exceeds the critical size of a single domain), their ground state could be inhomogeneous. Such a dot in the vortex state contains a magnetic singularity in the centre of the vortex, or the Belavin–Polyakov soliton type. The discrete magnon modes were observed for the uniformly magnetized dots by a resonance technique; their theory was constructed using the Landau–Lifshitz equations.

In this paper we consider the magnon modes which exist in 2d Heisenberg EA FM with a topological soliton. The model of an easy–axis ferromagnet is presented in Sec. II. An approximate solution for the soliton structure is examined. In Sec. III we formulate equations, describing magnons and solitons. In Sec. IV truly local magnon modes are calculated analytically by variational approach for the limiting cases of large and small radius solitons, and numerically for arbitrary radii. In Sec. V the problem of soliton–magnon scattering is formulated and treated numerically for different partial modes. The phenomenology of a pseudopotential is proposed in Sec. VI to analyze analytically both local states and the scattering problem. In Sec. VII we use results for local modes to derive the equations of soliton motion. Such equations are of Newtonian type (with the account of gyroscopic force), in contrast to both the case of the BP–soliton in the isotropic magnet and the case of the vortex in the EP magnet. In the framework of this equation the soliton effective mass is calculated. A discussion and concluding remarks are presented in Sec. VIII.

## II The Model and Elementary Excitations

The dynamics of the classical ferromagnet is described by the Landau–Lifshitz equations for the normalized magnetization $`𝒎`$. In angular variables, $`m_x+im_y=\mathrm{sin}\theta \mathrm{exp}(i\varphi )`$, these equations correspond to the Lagrangian

$$L=Ad^2x\left\{\frac{1}{D}\left(1\mathrm{cos}\theta \right)\frac{\varphi }{t}𝒲(\theta ,\varphi )\right\}$$
(1)
with the energy density $`A𝒲(\theta ,\varphi )`$,

$$𝒲(\theta ,\varphi )=\frac{1}{2}\left(\theta \right)^2+\frac{1}{2}\left(\varphi \right)^2\mathrm{sin}^2\theta +\frac{1}{2\mathrm{\Delta }_0^2}\mathrm{sin}^2\theta ,$$
(2)
where $`A=JS^2`$, $`J`$ is the exchange integral, $`S`$ the atomic spin, $`\mathrm{\Delta }_0=\sqrt{A/K}`$ the characteristic scale (“magnetic length”), $`K`$ the energy of the anisotropy, and $`D`$ the spin–wave stiffness.

The model of EA FM has well-known magnon excitations $`\theta =\text{const}1`$, $`\varphi =\omega t𝒌𝒓`$ above the homogeneous ground state. The dispersion law for the magnons has a finite activation frequency $`\omega _0=D/\mathrm{\Delta }_0^2`$,

$$\omega (𝒌)=\omega _0+D|𝒌|^2,$$
(3)
where $`𝒌`$ is the wave vector.

The simplest nonlinear excitation in the system is the dynamical (precessional) topological soliton

$$\varphi =\phi _0+q\chi +\mathrm{\Omega }t,\theta =\theta _0(r),\theta _0(0)=\pi ,\theta _0(\mathrm{})=0,$$
(4)
where $`r`$ and $`\chi `$ are polar coordinates in the plane of the magnet, the integer $`q`$ plays the role of a topological charge, and $`\mathrm{\Omega }`$ ($`0<\mathrm{\Omega }<\omega _0`$) is the frequency of the internal precession. The necessity to consider solitons with internal precession is caused by the fact that, according to the Derrick–Hobart theorem, static solitons in models like (2) are unstable. Formally, there is no function $`\theta _0(r)`$, which could provide a minimum of the energy

$$E=\pi A_0^{\mathrm{}}r𝑑r\left\{\left(\frac{d\theta _0}{dr}\right)^2+\mathrm{sin}^2\theta _0\left(\frac{1}{\mathrm{\Delta }_0^2}+\frac{q^2}{r^2}\right)\right\}.$$
(5)
The precessional soliton corresponds to the conditional minimum of the energy (5) for a given number $`N`$ of magnons bound in the soliton

$$N=\frac{N_2}{\mathrm{\Delta }_0^2}_0^{\mathrm{}}r𝑑r\left(1\mathrm{cos}\theta _0(r)\right),$$
(6)
which is conserved for the uniaxial case ($`N_2=2\pi A/\mathrm{}\omega _0`$ is the characteristic number of bound magnons in the 2d soliton), see also Refs. . This corresponds to the condition $`\delta \left[EN\mathrm{}\mathrm{\Omega }\right]=0`$. The form of the function $`\theta _0(r)`$ is defined by the ordinary differential equation,

$`{\displaystyle \frac{d^2\theta _0}{dr^2}}+{\displaystyle \frac{1}{r}}{\displaystyle \frac{d\theta _0}{dr}}\mathrm{sin}\theta _0\mathrm{cos}\theta _0\left({\displaystyle \frac{1}{\mathrm{\Delta }_0^2}}+{\displaystyle \frac{q^2}{r^2}}\right)`$ (7)
$`+{\displaystyle \frac{\mathrm{\Omega }}{\omega _0\mathrm{\Delta }_0^2}}\mathrm{sin}\theta _0=0,`$ (8)
with the boundary conditions (4). The precession frequency $`\mathrm{\Omega }`$ is fully determined by the number of bound magnons $`N`$. The dependence $`\mathrm{\Omega }(N)`$ was calculated numerically in Ref. ; large radius solitons ($`NN_2`$) have a very low precession frequency ($`\mathrm{\Omega }\omega _0`$); and $`\mathrm{\Omega }\omega _0`$ for small radius solitons ($`NN_2`$).

In the case of large radius solitons, $`R\mathrm{\Delta }_0`$, the approximate “domain wall” solution

$$\mathrm{cos}\theta _0(r)=\mathrm{tanh}\frac{rR}{\mathrm{\Delta }_0}$$
(9)
is applicable, which describes a 1d domain wall at the distance $`R`$ from the origin. One can obtain the relation between the soliton radius $`R`$ and the frequency of the soliton precession $`\mathrm{\Omega }`$ from variational condition, in which $`R`$ plays the role of a fitting parameter for the function (9). In the first approximation $`\mathrm{\Omega }(R)\omega _0\mathrm{\Delta }_0/R`$.

The opposite case of small radius solitons, $`R\mathrm{\Delta }_0`$, corresponds to the Belavin–Polyakov solution $`\mathrm{tan}(\theta _0/2)=(R/r)^{|q|}`$, which describes the isotropic limit of the model (1); it is valid for distances $`r\mathrm{\Delta }_0`$. At larger distances the algebraic decay is replaced by the general exponential form,

$$\theta _0(r)\mathrm{exp}\left(r/r_0\right),r_0=\sqrt{\frac{D}{\omega _0\mathrm{\Omega }}}.$$
(10)
The dependence $`\mathrm{\Omega }(R)`$ was calculated both analytically and numerically in Refs. . The frequency of the soliton precession satisfies $`\mathrm{\Omega }\omega _0`$ when $`R0`$, but the dependence $`\mathrm{\Omega }(N)`$ is not analytic, and $`d\mathrm{\Omega }/dN\mathrm{}`$ as $`N/N_20`$.

To describe approximately the soliton structure in the intermediate case of arbitrary $`R`$, we use a fitting method. We choose a trial function of the form

$$\mathrm{tan}\frac{\theta _0^{\text{(fit)}}(r)}{2}=\frac{R^{\text{(fit)}}}{r}\mathrm{exp}\left[\frac{\left(rR^{\text{(fit)}}\right)}{r_0}\right],$$
(11)
which has valid asymptotics for $`r0`$ and $`r\mathrm{}`$ (we discuss the case $`q=1`$ only). The parameter $`R^{\text{(fit)}}`$ can be found numerically by fitting the trial function (11) to the numerical solution of Eq. (8). Such an approximation shows that the trial function describes the soliton shape with an accuracy of about $`10^2`$, see Table I. From this table one can see that the values of the fitting parameter $`R^{\text{(fit)}}`$ are very close to the “exact” soliton radii; it should be understood by the fact that $`\theta _0\left(r=R^{\text{(fit)}}\right)=\pi /2`$. Numerical data for $`\theta _0(r)`$ and fitting results for $`\theta _0^{\text{(fit)}}(r)`$ are plotted in Fig. 1.

## III Magnon modes

To analyze magnons on the soliton background, it is convenient to introduce local coordinates $`\{𝒆_1,𝒆_2,𝒆_3\}`$, which describe the configuration of the magnetization unit vector $`𝒎`$ in the unperturbed soliton: $`𝒆_3`$ coincides with $`𝒎`$ in the soliton (4), $`𝒆_1=𝒆_y\mathrm{cos}\varphi _0𝒆_x\mathrm{sin}\varphi _0`$, and $`𝒆_2=𝒆_3\times 𝒆_1`$. Then the linear oscillations around the soliton can be described in terms of projections of the magnetization $`𝒎`$ on the local axes: $`\mu =𝒎𝒆_1`$ and $`\vartheta =𝒎𝒆_2`$. In the absence of the soliton, $`\theta _0=0`$, $`\varphi _0=\mathrm{\Omega }t`$, such oscillations correspond to free magnons in the rotating coordinate frame.

The linearized equations for $`\vartheta `$ and $`\mu `$ can be presented in the form:

$`\left[^2+𝒰_1(r)\right]\vartheta +{\displaystyle \frac{2q\mathrm{cos}\theta _0}{r^2}}{\displaystyle \frac{\mu }{\chi }}={\displaystyle \frac{1}{\omega _0\mathrm{\Delta }_0^2}}{\displaystyle \frac{\mu }{t}},`$ (13)
$`\left[^2+𝒰_2(r)\right]\mu {\displaystyle \frac{2q\mathrm{cos}\theta _0}{r^2}}{\displaystyle \frac{\vartheta }{\chi }}={\displaystyle \frac{1}{\omega _0\mathrm{\Delta }_0^2}}{\displaystyle \frac{\vartheta }{t}}`$ (14)
with the “potentials”

$`𝒰_1(r)`$ $`=`$ $`\left({\displaystyle \frac{1}{\mathrm{\Delta }_0^2}}+{\displaystyle \frac{q^2}{r^2}}\right)\mathrm{cos}2\theta _0(r){\displaystyle \frac{\mathrm{\Omega }}{\omega _0\mathrm{\Delta }_0^2}}\mathrm{cos}\theta _0(r),`$ (16)
$`𝒰_2(r)`$ $`=`$ $`\left({\displaystyle \frac{d\theta _0(r)}{dr}}\right)^2+\mathrm{cot}\theta _0(r)_r^2\theta _0`$ (17)
$`=`$ $`\left({\displaystyle \frac{d\theta _0(r)}{dr}}\right)^2+\left({\displaystyle \frac{1}{\mathrm{\Delta }_0^2}}+{\displaystyle \frac{q^2}{r^2}}\right)\mathrm{cos}^2\theta _0(r)`$ (19)
$`{\displaystyle \frac{\mathrm{\Omega }}{\omega _0\mathrm{\Delta }_0^2}}\mathrm{cos}\theta _0(r),`$
where $`_r^2`$ is the radial part of the Laplace operator.

It is convenient to represent the solutions of Eqs. (III) in terms of a partial–wave expansion for $`\vartheta `$ and $`\mu `$

$`\vartheta `$ $`=`$ $`{\displaystyle \underset{\alpha }{}}\left(u_\alpha (r)+v_\alpha (r)\right)\mathrm{cos}(m\chi +\stackrel{~}{\omega }_\alpha t)`$ (21)
$`\mu `$ $`=`$ $`{\displaystyle \underset{\alpha }{}}\left(u_\alpha (r)v_\alpha (r)\right)\mathrm{sin}(m\chi +\stackrel{~}{\omega }_\alpha t),`$ (22)
where $`\alpha =(k,m)`$ is a full set of eigennumbers, $`k`$ and $`m`$ being radial and azimuthal quantum numbers, respectively. The quantity $`\stackrel{~}{\omega }=\omega \mathrm{\Omega }`$ is the magnon frequency in the rotating frame; $`k`$ and $`\stackrel{~}{\omega }`$ are connected by the dispersion law, cf. Eq. (3),

$$\stackrel{~}{\omega }(𝒌)=\omega _{\text{gap}}+D|𝒌|^2,\omega _{\text{gap}}=\omega _0\mathrm{\Omega },$$
(23)
where $`\omega _{\text{gap}}`$ is the magnon gap frequency in the rotating frame. Below we will replace $`\stackrel{~}{\omega }`$ by $`\omega `$, almost everywhere without confusion, because we will work in the rotating frame mainly; the index $`\alpha `$ will be omitted, too. Thus we then finally obtain the following eigenvalue problem (EVP) for the functions $`u`$ and $`v`$:

$`\widehat{H}_1u`$ $``$ $`\left[\widehat{H}_0+{\displaystyle \frac{m^2}{r^2}}+V(r){\displaystyle \frac{\omega }{\omega _0\mathrm{\Delta }_0^2}}\right]u=W(r)v,`$ (24)
$`\widehat{H}_2u`$ $``$ $`\left[\widehat{H}_0+{\displaystyle \frac{m^2}{r^2}}V(r)+{\displaystyle \frac{\omega }{\omega _0\mathrm{\Delta }_0^2}}\right]v=W(r)u,`$ (25)
where $`\widehat{H}_0=_r^2+U_0(r)`$ is the 2d radial Schrödinger operator; the “potentials” are

$$U_0(r)=\frac{1}{2}\left(𝒰_1+𝒰_2\right),V(r)=\frac{2qm\mathrm{cos}\theta _0}{r^2}.$$
(27)
Note, that the Eqs. (23) are invariant under the conjugations $`\omega \omega `$, $`mm`$, and $`uv`$. In a classical theory we can choose any sign of the frequency; but in order to make contact with quantum theory, with a positive frequency and energy $`_k=\mathrm{}|\omega _k|`$, we will discuss the case $`\omega >0`$ only. Thus, there are two different sets of Eqs. (23) for $`m=|n|`$ and $`m=|n|`$. Besides there appears a difference between Eqs. (24) and (25). Since (25) has the asymptotically equivalent form, $`_r^2v+\omega v=0`$, the function $`v`$ has an exponential behaviour; at the same time Eq. (24) yields oscillating solutions. So, by choosing $`\omega >0`$, the variable $`v`$ becomes a slave variable in the EVP (23).

The “coupling potential” $`W`$,

$`W(r)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(𝒰_2𝒰_1\right)`$ (28)
$`=`$ $`{\displaystyle \frac{1}{2}}\left[\mathrm{sin}^2\theta _0\left({\displaystyle \frac{1}{\mathrm{\Delta }_0^2}}+{\displaystyle \frac{1}{r^2}}\right)\left({\displaystyle \frac{d\theta _0}{dr}}\right)^2\right],`$ (29)
is positive, it has a maximum value at a distance of about the soliton radius $`R`$ from the origin, vanishing both at $`r=0`$ and far from the soliton. Using the “domain wall” approximation (9), one obtains that the maximum value of $`W1/R^2`$, when $`R\mathrm{\Delta }_0`$. In the opposite case of small $`R`$, in the range $`rR`$ the soliton looks like the Belavin–Polyakov solution; here $`r(d\theta _0/dr)\mathrm{sin}\theta _0`$, and $`W=0`$ in the main exchange approximation. Therefore, one can suppose that $`W1/\mathrm{\Delta }_0^2`$ for any $`R`$. In fact, this is confirmed by numerical calculations, see Fig. 2.

However, one can see below in Sec. IV that the problem is impoverished when the “coupling potential” $`W(r)`$ is ignored; even in the limits of small (large) soliton radii, the main properties of local modes vanish, if we assume that $`W(r)=0`$. The reason is that a set of more than one coupled Schrödinger–like equations has properties, which are absent from a single equation. To explain this fact, remember that for a single 1d Schrödinger equation there exist two independent solutions; for the eigenfrequencies inside the continuous spectrum such solutions have oscillating asymptotics. For a set of 2, 3, … equations there exist 4, 6, … linearly independent solutions, respectively; some of them can have exponential asymptotics. Naturally, the real modes have an oscillational form here; we will use this fact below for the numerical analysis. But it is necessary to take into account the exponential solutions, too. In particular, as was shown in Ref. , such exponential asymptotics could correspond to truly localized states inside the continuous spectrum, which are forbidden for equations of the Schrödinger form. In the reduced EVP with two independent equations like (23) with $`W=0`$, there is a continuous spectrum $`\omega >\omega _{\text{gap}}`$ for the function $`u(r)`$ ($`u0`$, $`v=0`$) and $`\omega <\omega _{\text{gap}}`$ for the function $`v(r)`$ ($`u=0`$, $`v0`$). There are no limitations for the lowest value of the discrete levels; in particular, the frequencies of the discrete spectrum could be less then $`\omega _{\text{gap}}`$ for the function $`u`$. At the same time, as follows from the general set of coupled Eqs. (23), frequencies of local modes lie in the region $`\omega _{\text{gap}}<\omega <\omega _{\text{gap}}`$. Therefore, it is necessary to solve the general EVP (23) for an exact analysis of the problem.

We should note here that the model of an isotropic magnet is a special case, in which the “coupling” potential is exactly equal to zero. To study an EVP like (23) with $`W=0`$, it is sufficient to take into account one equation only (namely, the equation for the variable $`u`$); the second one, which has unphysical solutions for $`\omega >0`$, was ignored in Ref.. That is why the results for EA FM, obtained in this article, and the results for the isotropic magnet differ strongly even for small anisotropy. One difference is caused by the presence of truly local modes. Another one is due to the peculiar structure of the quasi–continuous spectrum with $`\omega >\omega _{\text{gap}}`$ for a finite geometry: the spectrum of the EA FM has a system of doublets with opposite signs of $`m`$, while there are no doublets for the isotropic magnet.

## IV Local modes

According to the dispersion law (23), free magnons only exist when $`\omega >\omega _{\text{gap}}`$. It is naturally to look for local modes in the range $`0<\omega <\omega _{\text{gap}}`$ (remember that we discuss the case $`\omega >0`$ ). In principle, such local modes could exist inside the continuous spectrum, as for the case of the nonlinear $`\sigma `$–model, describing the EP antiferromagnet. However, our study shows that this is not the case for the EA FM. Existence of local modes in such a magnet was found numerically in Ref.  for a soliton with some special parameters; namely the mode with $`m=1`$ for the soliton with the topological charge $`q=1`$, and the mode with $`m=2`$ for $`q=2`$ were calculated. Such modes appear in field theories, which include soliton–like solutions.

To study local modes, note that near the soliton centre, $`r\mathrm{\Delta }_0`$, one can obtain the asymptotics of the magnon amplitudes in the following form:

$`u_m(r)`$ $`=`$ $`𝒜_m\left({\displaystyle \frac{r}{\mathrm{\Delta }_0}}\right)^{|qm|}\left[1+O\left(r^2/\mathrm{\Delta }_0^2\right)\right],`$ (31)
$`v_m(r)`$ $`=`$ $`\epsilon 𝒜_m\left({\displaystyle \frac{r}{\mathrm{\Delta }_0}}\right)^{|q+m|}\left[1+O\left(r^2/\mathrm{\Delta }_0^2\right)\right].`$ (32)
Here the coefficient $`𝒜`$ can be determined by the normalization; the presence of a nonzero factor $`\epsilon `$ is caused by the “interaction” $`W(r)`$ between Eqs. (24) and (25); its value cannot be found through this asymptotic expansion.

In the opposite case of large distance, $`r\mathrm{\Delta }_0`$, the EVP (23) has the following general solution (here and below we consider the most interesting case of solitons with unit topological charge $`q=1`$, which has lowest energy)

$`u_m(r)`$ $``$ $`A_mK_{|m+1|}(\varkappa _1^{}r)+B_mI_{|m+1|}(\varkappa _1^{}r)`$ (34)
$`v_m(r)`$ $``$ $`C_mK_{|m1|}(\varkappa _1^+r)+D_mI_{|m1|}(\varkappa _1^+r),`$ (35)
where $`K_\nu `$ and $`I_\nu `$ are Macdonald and modified Bessel functions respectively, and $`\varkappa _1^\pm =\sqrt{\left(\omega _0\mathrm{\Omega }\pm \omega \right)/D}`$. Note, that at large distances, $`z1`$, the modified Bessel functions have exponential behaviour, $`I_\nu (z)\mathrm{exp}(z)/\sqrt{z}`$, $`K_\nu (z)\mathrm{exp}(z)/\sqrt{z}`$, so we use the asymptotical condition $`u_m,v_m\mathrm{exp}(\varkappa _1^\pm r)/\sqrt{r}`$ to obtain the localized solutions.

Firstly note that for the special cases $`m=0,+1`$ there exist “zero–modes” with $`\omega =0`$. Such modes are caused by the internal symmetry of the problem. They have the same form as for EP magnets. One of the zero–modes, the so–called translational mode with $`m=+1`$,

$$u_{+1}=\frac{d\theta _0}{dr}\frac{\mathrm{sin}\theta _0}{r},v_{+1}=\frac{d\theta _0}{dr}+\frac{\mathrm{sin}\theta _0}{r}$$
(37)
describes the position shift of the soliton. Another one with $`m=0`$

$$u_0=\mathrm{sin}\theta _0,v_0=\mathrm{sin}\theta _0$$
(38)
is a rotational mode, which corresponds to the presence of an arbitrary parameter, the phase $`\phi _0`$, in the soliton structure (4). For the isotropic magnet, as well as for the EP magnet, similar modes exist, too; but there they are quasi–local due to their slow power–law decay, $`1/r`$, or $`1/r^2`$. In the EA magnets the decay has an exponential form.

The remaining modes could be investigated numerically; analytical results were obtained for some limiting cases only, see below. To study the EVP numerically, we use the two–parameter shooting method for the integration of the set of Eqs. (23); a similar one–parameter shooting scheme was proposed in Ref.  to study the vortex–magnon scattering problem in 2d antiferromagnets. Numerically, we “kill” two growing exponents in (IV) by choosing two shooting parameters $`\epsilon `$ and $`\omega `$. As a result, we have found numerically that there are a number of local modes. Modes with opposite signs of $`m`$ bind to doublets; this picture is well–pronounced for $`|m|>1`$ in the limit of large $`R`$, see Fig. 3.

Such local modes exist only for not too small soliton radii (not for large $`\mathrm{\Omega }`$). Eigenfrequencies for such modes with arbitrary azimuthal numbers $`m`$ ($`|m|>1`$) increase rapidly with the decreasing of the soliton radius, and reach the boundary of the continuous spectrum. It looks as if all such modes “leave” the region of the gap in the spectrum, if the soliton radius is smaller than some critical value $`R_m^\text{c}`$, see Table II. It is natural to suppose that all modes with $`|m|>1`$ transform to the quasi–local ones with decreasing of $`R`$; this fact can be tested numerically due to the singularities of the scattering amplitudes on such quasi–local modes, see the next section.

Thus, in the range $`R1.52\mathrm{\Delta }_0`$ there is only one local mode with non–zero frequency, $`m=1`$, see Fig. 3. That is why only the mode with $`m=1`$ was observed in Ref.  for the soliton with $`q=1`$.

The eigenfunctions of bound states $`u`$ and $`v`$ are localized near the soliton core. The functions $`u`$ and $`v`$ are plotted in Fig. 4 for the mode with $`m=1`$; the picture of well–localized eigenfunctions is specially pronounced for the large radius solitons, see Fig. 4(a).

In order to analyze the spectrum of local modes analytically, let us reformulate the EVP (23) as a variational problem for the functional

$$F[u,v]=_0^{\mathrm{}}r𝑑r\left(u\widehat{H}_1u+v\widehat{H}_2v2Wuv\right).$$
(39)
We start with large radius solitons, $`R\mathrm{\Delta }_0`$, when the “domain wall” approximation (9) is valid. In the region of interest, $`rR\mathrm{\Delta }_0`$, the operator $`\widehat{H}_0`$ from the EVP (23) has in the main approximation for $`\mathrm{\Delta }_0/R`$ a simple reflectionless potential,

$$\widehat{H}_0\widehat{L}\frac{d^2}{dr^2}+\frac{1}{\mathrm{\Delta }_0^2}\left[1\frac{2}{\mathrm{cosh}^2\left((rR)/\mathrm{\Delta }_0\right)}\right]$$

with a single localized state $`\psi _0`$, corresponding to the zeroth eigenvalue

$$\widehat{L}\psi _0=0,\psi _0(r)=\left(2R\mathrm{\Delta }_0\right)^{1/2}/\mathrm{cosh}\left((rR)/\mathrm{\Delta }_0\right).$$

The other potentials in the EVP can be considered as corrections to $`\widehat{L}`$; this is the reason to put the trial functions in the form $`u=a\psi _0`$, $`v=b\psi _0`$ with trial parameters $`a`$ and $`b`$. The variational approach with such simple trial functions leads to eigenfrequencies for the different modes in the form

$`\omega _m=m\alpha +\sqrt{\left(m^2\beta \delta \right)^2\gamma ^2},`$ (40)
$`\alpha =\psi _0\left|{\displaystyle \frac{2\mathrm{cos}\theta _0(r)}{r^2}}\right|\psi _0,\beta =\psi _0\left|{\displaystyle \frac{1}{r^2}}\right|\psi _0,`$ (41)
$`\gamma =\psi _0\left|W(r)\right|\psi _0,\delta =\psi _0\left|{\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}+\mathrm{\Omega }\mathrm{cos}\theta _0(r)\right|\psi _0`$ (42)
where $`\psi _0\left|\mathrm{}\right|\psi _0=_0^{\mathrm{}}\psi _0(\mathrm{})\psi _0r𝑑r`$. However, there is a difference between the real solutions $`u`$, $`v`$, and the trial function $`\psi _0`$. This difference provides, e.g., the existence of “zero–modes”. We do not like to improve the form of the eigenfunctions, but only the coefficients in the Eqs. (40), which describe the eigenfrequencies. Using the additional condition for “zero–modes”,

$$\omega _m=0\text{for }\text{m=0,+1},$$
(43)
the eigenfrequencies $`\omega _m`$ can be expressed through two independent coefficients only. It is convenient to present the eigenfrequencies in the form

$$\omega _m=m\alpha +|m|\sqrt{\alpha ^2+\beta ^2(m^21)}.$$
(44)
Using trial functions $`\psi _0`$, we can calculate the parameters $`\alpha `$ and $`\beta `$ for the large radius solitons,

$$\alpha =\omega _0(\mathrm{\Delta }_0/R)^3\beta =\omega _0(\mathrm{\Delta }_0/R)^2.$$
(46)
In this case the spectrum has the form

$$\frac{\omega _m}{\omega _0}=|m|\sqrt{\left(\frac{\mathrm{\Delta }_0}{R}\right)^6+\frac{m^21}{(R/\mathrm{\Delta }_0)^4}}m\left(\frac{\mathrm{\Delta }_0}{R}\right)^3$$
(47)
This formula corresponds very well to the numerical data for $`R\mathrm{\Delta }_0`$, see Fig. 3. In the main approximation for $`\mathrm{\Delta }_0/R`$ there is a system of doublets (for all modes with $`|m|>1`$) with quadratic dependence of the mean frequency on $`\mathrm{\Delta }_0/R`$,

$$\overline{\omega }_{\pm m}\omega _0|m|\sqrt{m^21}\left(\frac{\mathrm{\Delta }_0}{R}\right)^2,$$
(48)
and the splitting is a small cubic correction, $`\mathrm{\Delta }\omega _{\pm m}m\left(\mathrm{\Delta }_0/R\right)^3`$. For the special case $`m=1`$ the dependence $`\omega (R)`$ has the form

$$\omega _{m=1}2\omega _0\left(\frac{\mathrm{\Delta }_0}{R}\right)^3.$$
(49)
Eq. (47) describes qualitatively all modes with arbitrary $`R`$. In particular, it describes the fact, mentioned above, that modes reach the boundary of the continuous spectrum at some finite frequency of the soliton precession $`\mathrm{\Omega }_m^\text{c}`$. Let us estimate the positions of the crossover points, where the modes “leave” the discrete spectrum. Using Eq. (47), one can obtain the condition for the critical values of $`x=\mathrm{\Omega }_m^\text{c}/\omega _0`$:

$$|m|\sqrt{x^6+(m^21)x^4}mx^3=1x.$$
(50)
Results are very close to the numerical data, see Table II. We expect that the local modes transform into the quasi–local ones inside the continuous spectrum; but their frequencies can still be described by the formula (47), see Sec. VI.

For the small radius solitons, $`R\mathrm{\Delta }_0`$, there only exist local modes with $`|m|=1`$. Thus, we will consider this case. Constructing the trial functions, one can suppose that they should be similar to the modes with $`m=\pm 1`$. Eigenfunctions for the translational “zero–mode”, $`u_{+1}`$ and $`v_{+1}`$, are well–known, see Eq. (37). To construct the trial functions it is convenient to use functions $`u_1`$ and $`v_1`$ in the form

$$u_1=r^2u_{+1},v_1=r^2v_{+1},$$

which coincide with the exact “zero–mode” eigenfunctions for the mode with $`m=1`$ in the limiting case of the BP–soliton. Now we can choose the trial functions $`u`$ and $`v`$ as follows,

$$u=au_{+1}+bu_1,v=cv_{+1}+dv_{1.}$$
(51)
The variational problem (39) can be solved analytically in the main approximation for $`R/\mathrm{\Delta }_0`$ using the explicit BP–solutions. Tedious calculations lead to the eigenfrequency for small $`R`$,

$$\omega _{m=1}C\left(\omega _0\mathrm{\Omega }\right).$$
(52)
Unfortunately, this method is not exact and leads to a value of the constant $`C`$ which is larger than the numerical value. This is caused by the fact that the trial functions, which are constructed from zero modes, have the wrong asymptotical behaviour (10) at large distance, $`\theta _0(r)`$, $`\theta _0^{}(r)`$ $`\mathrm{exp}(r/r_0)`$, instead of (IV), i.e. $`u\mathrm{exp}(\varkappa _1^{}r)`$, $`v\mathrm{exp}(\varkappa _1^+r)`$. To improve this result, one can change the asymptotical behaviour of the functions (51) by exponential factors:

$$uue^{r\left(1/r_0\varkappa _1^{}\right)},vve^{r\left(1/r_0\varkappa _1^+\right)}.$$

In this case the calculations can be done numerically only; the value of the constant $`C1`$. For the approximate description of the eigenfrequency of the mode with $`m=1`$ in the intermediate case of arbitrary $`R`$, we use the trial function

$$\omega _{m=1}^{(trial)}(\zeta )=\omega _0\frac{2\zeta ^3(1\zeta )}{1\zeta +2\zeta ^4},\zeta =\frac{\mathrm{\Omega }}{\omega _0}.$$
(53)
Providing the asymptotically correct behaviour, (49) and (52), this approximation reproduces to the numerical results with an accuracy of about $`210^3`$, see Fig. 3.

## V The scattering problem: numerical results

Let us describe the scattering of magnons by a soliton. Note that without a soliton, free magnon modes have the form

$$u_m(r)J_m(kr),v_m(r)=0,$$
(54)
where $`k`$ is a “radial wave number”, $`k`$ and $`\omega `$ are connected by the dispersion law (23), and $`J_m`$ are Bessel functions. The free–modes like $`u_m`$ play the role of the partial cylinder waves of a plane spin wave

$$\mathrm{exp}\left(i𝒌𝒓i\omega t\right)=\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}i^mJ_m(kr)e^{im\chi i\omega t}.$$
(55)
To describe magnon solutions in the presence of a soliton, one should note that far from the soliton, $`r\mathrm{\Delta }_0`$, the “coupling potential” $`W(r)`$ is exponentially small. Therefore there are two asymptotically independent EVP; the magnon amplitude $`v_m`$ has the same form (35) as for the local modes, but the function $`u_m`$ shows an oscillating behaviour like

$$u_m(r)J_{|m+1|}(kr)+\sigma _m(k)Y_{|m+1|}(kr),$$
(56)
where $`Y_{|m+1|}`$ are Neumann functions. The quantity $`\sigma (k)`$ derives from the soliton–magnon scattering; it can be interpreted as the scattering amplitude.

We use the one–parameter shooting method to study numerically the problem of soliton–magnon scattering, as described in Ref.  for vortices in the antiferromagnet. Choosing the shooting parameter $`\epsilon `$, see Eq. (IV), we “kill” the growing exponent in (35); as a result we have obtained a well–pronounced exponential decay for $`v_m`$; the corresponding oscillating solutions for $`u_m`$ differ from the asymptotical values (56) only in the vicinity of the soliton core, see Fig. 5. The scattering amplitude was found from these data by comparison with the asymptotics (56). The results are the following:

For all modes the scattering amplitude $`\sigma _m(k)`$ tends to zero as $`k0`$. In the long–wavelength limit the maximal scattering is related to the mode with $`m=1`$, for which the behaviour of $`d\sigma /dk`$ looks singular. However, other properties depend strongly on the soliton radius.

In the case of large $`R`$, the scattering amplitude $`\sigma _m(k)`$ is positive for all modes in the long–wavelength limit. There is a maximum of $`\sigma _m(k)`$ at about $`kR1`$ for modes with positive $`m`$; but there are two poles in the scattering amplitude for modes with negative $`m`$, and one pole for $`m=0`$. Naturally, there is no real divergence at such a pole: the physically observed phase of the scattering, or the phase shift $`\delta _m=\mathrm{tan}\sigma _m`$, varies monotonically, see Fig. 6. Thus, there are principally different pictures for positive and negative $`m`$: the total phase shift is

$$\delta _m^{tot}\delta (\mathrm{})\delta (0)=\{\begin{array}{cc}0,\hfill & m>0\hfill \\ \pi ,\hfill & m=0\hfill \\ 2\pi ,\hfill & m<0\hfill \end{array},R\mathrm{\Delta }_0.$$
(58)
The scattering picture is quite different for the case of small radius solitons. Each mode with $`|m|>1`$ has a single pole only, the position $`k_pR`$ of which increases when $`m`$ grows. For the special case of the mode with $`m=+1`$, the scattering amplitude is anomalously small; for the mode with $`m=1`$ there appears a second pole with anomalously small values of $`k_p`$. These results are presented in Fig. 7. The total phase shift is

$$\delta _m^{tot}=\{\begin{array}{cc}+\pi ,\hfill & m>1\hfill \\ 0,\hfill & m=1\hfill \\ 2\pi ,\hfill & m=1\hfill \\ \pi ,\hfill & m=0,m<1\hfill \end{array},R\mathrm{\Delta }_0.$$
(59)
There appears discontinuous change of the total phase shift from (58) to (59) for all modes with $`m0,1`$ for the cases of intermediate soliton radii. This transition occurs independently for each mode at some critical values $`k_m^\text{c}`$ of the wave–number.

## VI The scattering problem in the framework of the pseudopotential method

For the analytical description of the scattering problem, let us analyze the EVP qualitatively, replacing the real “potentials” (III) in the EVP (23) by simplified ones, which are amenable to an exact treatment. Such an approach is well–known in the quantum physics of solids as the pseudopotential method. Sometimes it is useful even when the potential is not small, see Ref.  for details.

Let us discuss the choice of the pseudopotential. The simplest way is to describe the soliton structure in the main approximation, using the ansatz

$$\mathrm{cos}\theta _0(r)\text{sign }\xi ,\xi =\frac{rR}{\mathrm{\Delta }_0},$$
(60)
which is a good approximation for both small and large distances. Replacing the soliton by this configuration, which does not have any in–plane structure ($`\mathrm{sin}\theta _0=0`$), reduces the complex EVP (23) to a very simple “centrifugal model” with $`v_m0`$:

$$\left(_r^2+k^2+\frac{|m+\text{sign }\xi |^2}{r^2}\right)u_m=0.$$
(61)
This model describes quasi–free magnons in each of the regions $`r<R`$ and $`r>R`$. The only effect of the soliton–magnon interaction is a shift of the mode indexes in comparison with the free magnons (54):

$$u_m(r)\{\begin{array}{cc}J_{|m1|}(kr),\hfill & r<R,\hfill \\ J_{|m+1|}(kr)+\sigma _m(k)Y_{|m+1|}(kr),\hfill & r>R,\hfill \end{array}$$
(62)
The usual matching condition for these solutions has the form

$$\left[\frac{u^{}}{u}\right]_R=0,$$
(63)
where $`\left[\mathrm{}\right]_R(\mathrm{})|_{R+0}(\mathrm{})|_{R0}`$. Calculations lead to the scattering amplitude in the form:

$$\sigma (\varkappa )=\frac{J_{|m1|}^{}(\varkappa )J_{|m+1|}(\varkappa )J_{|m+1|}^{}(\varkappa )J_{|m1|}(\varkappa )}{J_{|m1|}(\varkappa )Y_{|m+1|}^{}(\varkappa )J_{|m1|}^{}(\varkappa )Y_{|m+1|}(\varkappa )},$$
(64)
where $`\varkappa =kR`$, see Fig. 8. The total phase shift in this model

$$\delta _m^{tot}=\pi \text{sign }m.$$
(65)
### A Small radius solitons, $`𝑹\mathbf{}𝚫_\mathrm{𝟎}`$

The simple “centrifugal model” with the results (64) explains qualitatively the scattering data for the small radius solitons, except for modes with $`m=0,\pm 1`$, cf. (59). The cause for these exceptions is the existence of truly local modes (“zero–ones” for $`m=0,+1`$, and another one with non–zero frequency for $`m=1`$); while such modes are absent for the other $`m`$. The existence of a local mode diminishes the total phase shift by $`\pi `$ (and decreases the number of states in the continuum spectrum by $`1`$, due to the conserved total number of magnon states). Such a result is well–known in 1d, see e.g. Ref. . Except for these cases, one can expect for the scattering peculiarities at $`k1/R>>1/\mathrm{\Delta }_0`$. In this range the approximation by an isotropic magnet and BP–solitons is valid. The scattering problem for the isotropic magnets was studied in detail in Ref. . Over a wide range of wave–numbers, $`1/\mathrm{\Delta }_0k1/R`$, the asymptotics for the Belavin–Polyakov soliton are valid,

$`\sigma _{m1,0}(k)`$ $`=`$ $`{\displaystyle \frac{\pi (kR)^2}{2|m|(m+1)}},`$ (67)
$`\sigma _{m=0}(k)`$ $`=`$ $`{\displaystyle \frac{\pi }{2\mathrm{ln}\left(1/kR\right)}},`$ (68)
$`\sigma _{m=1}(k)`$ $`=`$ $`\pi (kR)^2\mathrm{ln}\left(1/kR\right).`$ (69)
For larger wave–numbers, $`k|m|/R`$, the scattering amplitude has a general, so–called, eikonal dependence, $`\sigma _m1/kR`$,

$$\sigma _m(k)=\frac{\pi (m1)}{kR}.$$
(70)
Our numerical data in this region are very close to the asymptotical behaviour (VI A), see Fig. 7. Note that the maximal scattering in the region $`1/\mathrm{\Delta }_0k1/R`$ occurs for the rotational mode with $`m=0`$, whose asymptotics (68) has a singularity for $`k0`$. There is one pole only for each mode; its position corresponds to the pole in the Belavin–Polyakov solution, $`k_p|m|/R`$.

Let us consider the local modes, for which the total phase shift is diminished by $`\pi `$ in comparison to the “centrifugal model” (65); it is equal to $`0`$, $`\pi `$, $`2\pi `$ for the cases $`m=+1`$, $`0`$, $`1`$, respectively.

There is one pole only for the simplest case of the mode with $`m=0`$. The scattering amplitude in the range $`1/\mathrm{\Delta }_0k1/R`$ is described by the formula (68). The position of the pole is determined by the BP–result $`k_p=1/R`$. The reason for such a nice agreement with the model of the isotropic magnet will be explained below in the text.

There is no pole for the translational mode with $`m=+1`$. This is due to a unique property of the isotropic model, where such a mode does not scatter at all; the magnon amplitude in this limit has the form

$$u_{m=+1}(r)=J_2(kr)\frac{2}{kr}\frac{J_1(kr)}{(r/R)^2+1},$$
(71)
which agrees with our numerical data, see Fig. 11. This unique property of reflectionlessness is caused by the high internal symmetry of the problem.

Fig. 7(a) shows that such a property is approximately valid for the small radius soliton in the anisotropic magnet, too: the scattering amplitude is anomalously small, so the total phase shift is absent. Another effect of this symmetry in the isotropic magnet is that the scattering problem for the BP–model can be reduced to a simplified one with effectively smaller scattering. Such a “reduced model” has no local (and quasi–local) modes; it is described by shifted indexes, $`mm1`$. It is natural to use the “centrifugal” approximation for such a reduced model (64). As a result, $`\sigma _{m=+1}=0`$ in the framework of the same “centrifugal model”; for the other modes, Eq. (65) leads to the total phase shift $`\delta _m^{tot}=\pi \text{sign}(m1)`$. Taking into account that $`\text{sign }0=0`$, one can see that such a simple dependence for $`\delta _m^{tot}`$ describes all modes, except of $`m=1`$. Its behaviour is most interesting.

For the mode with $`m=1`$ there are two poles, while it should be one pole only for the Belavin–Polyakov solitons. Indeed one of the poles lies in the region where the Belavin–Polyakov approximation is valid; its value, $`k_p2.79/R1/\mathrm{\Delta }_0`$, coincides exactly with the Belavin–Polyakov data. The second pole lies in the region $`k1/\mathrm{\Delta }_0`$, where the approximation by the isotropic magnet is not applicable. This pole is caused by the presence of the local mode, the typical wave–number of such a mechanism is $`k_p\varkappa _1^\pm 1/\mathrm{\Delta }_0`$. The use of “reduced” or “centrifugal” models is certainly wrong for this case. We will discuss the scattering amplitude in this region later, on the basis of a more general pseudopotential model, see Eq. (92).

### B Large radius solitons, $`𝑹\mathbf{}𝚫_\mathrm{𝟎}`$. Refinement of the pseudopotential model

Let us consider the case of large radius solitons, $`R\mathrm{\Delta }_0`$. Here exist local modes for each $`m`$. Therefore the simple “centrifugal model” as well as other simple models without a discrete spectrum seem to be wrong. This can be explained by the fact that in the very simple model (61) we omit the “coupling potential”; the real potential $`W`$ tends to zero at the origin and at infinity, but has maxima at some finite distances.

Thus we need a more general form of the pseudopotentials which accounts for the potentials $`U_0`$, $`V`$, and $`W`$. To describe the soliton shape, we can again use the expression (60) for the out–of–plane soliton configuration. However, we need to consider the in–plane signature, too. According to (60), it could be described as

$$\mathrm{sin}^2\theta _0(r)\delta (\xi ).$$

Therefore we can suppose that the “coupling potential” $`W(r)`$ is nonzero only at some distance around the soliton radius and approximately replace it by a $`\delta `$–potential as follows,

$$W(r)\frac{c_1}{\mathrm{\Delta }_0^2}\delta (\xi ).$$
(73)
The same assumptions lead to simple expressions for the other potentials in the EVP (23),

$`U_0(r)+{\displaystyle \frac{m^2}{r^2}}\pm V(r)`$ (74)
$`{\displaystyle \frac{1(\mathrm{\Omega }/\omega _0)\text{sign }\xi }{\mathrm{\Delta }_0^2}}+{\displaystyle \frac{|m\pm \text{sign}\xi |^2}{r^2}}+{\displaystyle \frac{c_2}{\mathrm{\Delta }_0^2}}\delta (\xi ).`$ (75)
Here $`c_1`$ and $`c_2`$ are trial parameters.

We choose this pseudopotential model (VI B) out of the following reasons: (i) it allows an exact analytical solution; (ii) it guarantees the correct asymptotic behaviour of the solutions; (iii) it offers a possibility to identify the parameters of the pseudopotential from the discrete spectrum (44).

Let us consider the discrete spectrum in the framework of the pseudopotential (VI B). Comparison with the eigenfrequencies (44) provides a possibility to determine all trial parameters $`c_k`$. For the discrete part of the spectrum, the EVP (23) with potentials (VI B) has solutions like

$`u_m(r)\{\begin{array}{cc}I_{|m1|}(\varkappa _2^{}r),\hfill & r<R\hfill \\ K_{|m+1|}(\varkappa _1^{}r),\hfill & r>R\hfill \end{array},`$ (79)
$`v_m(r)\{\begin{array}{cc}I_{|m+1|}(\varkappa _2^+r),\hfill & r<R\hfill \\ K_{|m1|}(\varkappa _1^+r),\hfill & r>R\hfill \end{array},`$ (82)
$`\varkappa _1^\pm =\sqrt{\left(\omega _0\mathrm{\Omega }\pm \omega \right)/D},\varkappa _2^\pm =\sqrt{\left(\omega _0+\mathrm{\Omega }\pm \omega \right)/D}.`$ (83)
Now it is possible to calculate the eigenspectrum by matching the solutions (79) and (82) at the distance $`R`$,

$`\left[u^{}\right]_R={\displaystyle \frac{c_2}{\mathrm{\Delta }_0}}u|_R{\displaystyle \frac{c_1}{\mathrm{\Delta }_0}}v|_R,`$ (85)
$`\left[v^{}\right]_R={\displaystyle \frac{c_2}{\mathrm{\Delta }_0}}v|_R{\displaystyle \frac{c_1}{\mathrm{\Delta }_0}}u|_R.`$ (86)
Calculations lead to the eigenfrequencies in the form (40). Using the condition for the “zero–modes” (43), the eigenfrequencies can be represented in the form (44); the constants $`c_1`$, $`c_2`$ can be identified by comparison with the spectrum (44),

$$c_1=\frac{1}{2}\left(\beta \frac{\alpha ^2}{\beta }\right),c_2=2+c_1+m\alpha +m^2\beta \frac{\mathrm{\Delta }_0^2}{4R^2}.$$
(87)
The parameters $`\alpha `$ and $`\beta `$ were calculated by the variational approach, $`\alpha 1/R^3`$, $`\beta 1/R^2`$, see Eq. (46).

Let us discuss the scattering problem. Solving the EVP (23) with pseudopotentials (VI B), we obtain for the continuum spectrum, $`\omega >\omega _{\text{gap}}`$, the following solutions. The function $`u`$ for $`r<R`$ has the form

$$u(r)\{\begin{array}{cc}I_{|m1|}(\varkappa _2^{}r),\hfill & \text{when }k<\sqrt{2\mathrm{\Omega }/D},\hfill \\ J_{|m1|}(k_2^{}r),\hfill & \text{when }k>\sqrt{2\mathrm{\Omega }/D},\hfill \end{array}$$
(88)
where $`k_2^{}=\sqrt{k^22\mathrm{\Omega }/D}`$. At large distances, $`r>R`$, it has the usual oscillating form (56). The function $`v`$ has the localized form (82). One can find the scattering amplitude by matching the solutions (56), (88), and (82) at the distance $`R`$, using the condition (VI B):

$$\sigma _m(k)=\frac{\left[_m(k)+𝔉_m(k)\right]J_{|m+1|}(kR)k\mathrm{\Delta }_0J_{|m+1|}^{}(kR)}{k\mathrm{\Delta }_0Y_{|m+1|}^{}(kR)\left[_m(k)+𝔉_m(k)\right]Y_{|m+1|}(kR)},$$
(89)
here

$`_m(k)`$ $`=`$ $`c_2c_1^2(c_2+\varkappa _2^+\mathrm{\Delta }_0{\displaystyle \frac{I_{|m+1|}^{}(\varkappa _2^+R)}{I_{|m+1|}(\varkappa _2^+R)}}`$
$``$ $`\varkappa _1^+\mathrm{\Delta }_0{\displaystyle \frac{K_{|m1|}^{}(\varkappa _1^+R)}{K_{|m1|}(\varkappa _1^+R)}})^1,`$
$`𝔉_m(k)`$ $`=`$ $`\{\begin{array}{cc}\varkappa _2^{}\mathrm{\Delta }_0{\displaystyle \frac{I_{|m1|}^{}(\varkappa _2^{}R)}{I_{|m1|}(\varkappa _2^{}R)}},\hfill & k<\sqrt{2\mathrm{\Omega }/D},\hfill \\ k_2^{}\mathrm{\Delta }_0{\displaystyle \frac{J_{|m1|}^{}(k_2^{}R)}{J_{|m1|}(k_2^{}R)}},\hfill & k>\sqrt{2\mathrm{\Omega }/D}.\hfill \end{array}`$
It is easy to analyze the expression for the scattering amplitude (89) in some limiting cases. In the long–wavelength limit, $`k1/R1/\mathrm{\Delta }_0`$, one can use the asymptotics for the cylindrical functions $`J(z)`$, $`Y(z)`$ at $`z1`$ and, at the same time, the asymptotics for $`I(z)`$, $`K(z)`$ for $`z1`$. Simple calculations show that the scattering intensity is maximal for the mode with $`m=1`$,

$$\sigma _{m=1}(k)=\frac{\pi }{2\mathrm{ln}\left(1/kR\right)},kR1,$$
(92)
which agrees with the numerical data, see Fig. 6(b). Note, that the Eq. (92) was obtained using one condition only, $`kR0`$. Therefore it is valid for the case of small $`R`$, too, explaining the existence of the second pole for the mode $`m=1`$ both for large and small radii.

For the other values of $`m`$, we can restore the general dependence of the amplitude of the magnon scattering on the large radius soliton, $`R\mathrm{\Delta }_0`$:

$$\sigma _{m1}(k)=\frac{\pi \left(kR/2\right)^{2|m+1|}}{|m+1|!(|m+1|1)!},kR1.$$
(93)
The scattering amplitude, according to Eq. (93), is positive for all modes in the long–wavelength limit.

In the range $`1/Rk1/\sqrt{R\mathrm{\Delta }_0}`$ the asymptotics for $`J(z)`$, $`Y(z)`$ for $`z1`$ are valid, which in the main approximation for $`\mathrm{\Delta }_0/R`$ leads to a linear dependence of the scattering phase shift on $`kR`$,

$$\delta _m(k)\frac{|m+1|\pi }{2}\frac{\pi }{4}kR,1kR\sqrt{\frac{R}{\mathrm{\Delta }_0}}.$$
(94)
For this range of $`k`$, all curves look like a family of parallel lines with equal distance of about $`\pi /2`$. Note that the real values of the trial parameters $`c_1`$ and $`c_2`$ lead to the corrections in Eqs. (VI B) like $`\mathrm{\Delta }_0/R`$.

Naturally the pseudopotential method cannot explain the scattering results for $`k1/\mathrm{\Delta }_0`$, where the influence of truly local modes and the differences between the real potential and the $`\delta `$–function are dominant. Nevertheless, the method gives quite good results even for the case $`k1/\mathrm{\Delta }_0`$, where the linear asymptotics does not work. Corresponding curves are plotted in Fig. 9 using the general dependence (89) with trial parameters $`c_1`$ and $`c_2`$. Fitting of these parameters gives a result in the zeroth–order approach for $`\mathrm{\Delta }_0/R`$: $`c_10`$, $`c_20.4`$ instead of $`c_10`$, $`c_22`$, according to (87). The difference is caused by the limitations of the model (VI B). For example, the above mentioned model leads to an increase of the scattering phase when $`k`$ increases; the total phase shift for this model is $`\delta _m^{tot}=\text{sign }m\pi `$, as for the “centrifugal model”, see (65), i.e. the qualitative properties of these models are equivalent. The appearance of a local mode diminishes the total phase shift by $`\pi `$. As a result, there is no total phase shift for all modes with positive $`m`$; but $`\delta ^{tot}=2\pi `$ for all modes with $`m<0`$.

The scattering picture in the short–wavelength limit, $`k|m|/R`$, gives the general, so–called eikonal dependence $`\sigma _mA_m/k`$. However, in contrast to the case of the isotropic magnet (70), the scattering amplitude has the same sign $`A_m<0`$ for all modes due to the presence of the effective attractive potential, which causes the existence of local modes. When the soliton radius decreases and the local mode disappears, the sign of $`\sigma _m`$ changes; the result becomes similar to (70).

Thus, we explain the particularities in the scattering picture of EA magnets, in comparison with the isotropic case, by the influence of local modes. To examine this influence, let us consider the scattering data $`\sigma _m`$ for the modes $`m=+2`$ and $`m=2`$ with different soliton radii $`R`$ as functions of the parameter $`kR`$, i.e. $`\sigma =\sigma _m(kR)`$; this gives a possibility to compare the scattering data with the isotropic model, see Fig. 10. For the small radii the scattering picture is similar to the BP–case (solid lines in Fig. 10). When the soliton radius increases, the scattering phase jumps up by $`\pi `$, decreasing the total phase shift. Such a bifurcation takes place for the mode with $`m=+2`$ at $`R=2.24\mathrm{\Delta }_0`$, and for the mode with $`m=1`$ at $`R=1.52\mathrm{\Delta }_0`$. These values agree with the crossover points $`R_m^\text{c}`$, see Table II. Thus, it is natural to suppose that the bifurcation is connected with the appearance of the quasi–local modes inside the continuum spectrum for $`R>R_m^\text{c}`$. Simple calculations show that the positions $`k_p`$ of the poles correspond to the frequencies of quasi–local modes, which can be estimated by Eq. (47).

## VII Translational modes and description of the soliton dynamics

The foregoing analysis of the spectrum for the modes with $`|m|=1`$ can be used to describe the translational motion of the soliton. Such an approach was first proposed in Ref. , where an effective mass of the vortex in a 2d EP magnet was calculated from numerical data for the eigenmode spectrum of the finite size magnet. An important progress for the vortices in FM was achieved in Ref. ; the specific structure of the spectrum of modes with $`|m|=1`$ (namely the existence of a Goldstone mode with $`\omega 1/L^2`$, and a number of doublets with mean frequency $`\overline{\omega }1/L`$ and a small splitting) leads to non–Newtonian equations of motion with $`3^{\text{rd}}`$, $`5^{\text{th}}`$, … time derivatives. Such equations were derived and discussed earlier in Ref. , using a phenomenological approach, for a review see Ref. . The analytical calculation verified their adequacy; besides, it allows to calculate the coefficients in the equations of motion with an error of only about $`0.8\%`$. Thus, it is possible to describe in this way such unusual properties as the appearance of the non–Newtonian equations of motion.

The spectrum of magnon modes with $`|m|=1`$ has a different structure for the isotropic FM. First of all, there are no doublets; for an approximate description of the dynamics only two modes are sufficient: the translational Goldstone mode with dependence $`\omega 1/L^2`$ and another mode with $`m=1`$. Therefore, the effective equations of motion have a Newtonian form with an effective mass, which diverges as $`L^2`$, see Ref. . Such behaviour agrees with the direct calculations done by Zaspel.

Let us note that all magnets mentioned above, i.e. EP and isotropic systems, have a gapless dispersion law and, as a result, a strong interaction of the soliton with the boundaries of the magnet; the interaction force is proportional to $`1/L`$ for the EP vortices and $`1/L^2`$ for the Belavin–Polyakov solitons. The absence of a gap manifests itself in the particularities of the dynamical coefficients (the effective mass $`M`$ for the BP–soliton in the isotropic FM and the coefficient $`G_3`$ in the term with $`\left[𝒆_z\times ^3𝑿/t^3\right]`$ for the vortex in the EP FM). This can be explained as follows. The existence of the mass $`M`$ and the gyroscopical force $`G\left[𝒆_z\times 𝑿/t\right]`$ leads to a finite frequency of the Larmor precession of the soliton. For a gapless dispersion law, this frequency lies in the continuous spectrum. Thus, the Larmor motion of the vortex leads to the generation of magnons. For a finite–size system without dissipation, such magnons are distributed through the whole magnet; the Larmor dynamics of the soliton centre is strongly coupled with a “magnon cloud”, having a scale of about the system size. Therefore, the position of the soliton centre, which can be determined as a point with $`m_z=1`$ (BP–soliton), or $`m_z=\pm 1`$ (vortex), in fact plays the role of a collective variable which governs the motion of the “magnon cloud”. Then it is not surprising that the corresponding equations of motion are nonlocal, leading to the divergence of the dynamic coefficients as $`L\mathrm{}`$.

The situation is quite different for the case of the EA FM. Our numerical analysis and analytical calculations show that the picture of doublets for local modes with $`|m|`$ and $`|m|`$ is valid for large $`R`$, see Sec. IV. On the other hand, the eigenfrequencies of these modes lie in the gap for the large soliton radii, $`\omega <\omega _{\text{gap}}`$, so one can expect the existence of a finite mass in the effective equations of the soliton motion. Therefore for the phenomenological description of the soliton dynamics we can use the $`2^{\text{nd}}`$–order differential equation

$$M\frac{^2𝑿}{t^2}G\left[𝒆_z\times \frac{𝑿}{t}\right]=𝑭_e\left(𝑿\right),$$
(95)
which corresponds to the picture with the lowest doublet of local modes. Note that further generalizations for the finite size magnet take into account the next doublets of quasi–continuous modes. Here the hierarchy of the effective equations of motion containing only even–order time derivatives manifests itself. In Eq. (95) $`𝑿`$ describes the position of the soliton, $`M`$ is the mass coefficient, $`G`$ the gyrocoefficient, $`𝒆_z`$ the unit vector along the easy ($`z`$) axis, and $`𝑭_e`$ the external force acting on the soliton due to the boundary or other solitons. Assuming that deviations of the soliton from the equilibrium position are small, $`𝑭_e\alpha 𝑿`$, the effective equation of motion (95) can be solved by a harmonic ansatz, which leads to

$$\omega _{0,1}=\frac{G}{2M}\pm \sqrt{\frac{G^2}{4M^2}+\frac{\alpha }{M}}.$$
(96)
For the localized solitons, the interaction between soliton and boundary (or between two solitons) has an exponential decay as $`\mathrm{exp}(L/r_0)`$, $`r_0=\sqrt{D/\omega _{\text{gap}}}`$, see Eq. (10). Thus $`\alpha /M`$ can be neglected for a large system size $`L`$, so the frequencies (96) have the form:

$$\omega _0=0,\omega _1=G/M.$$
(97)
The zero frequency $`\omega _0`$ corresponds to the Goldstone mode (position shift of the soliton), in which all spins rotate with the frequency of the soliton precession $`\mathrm{\Omega }`$ in the laboratory frame. The presence of the rotating frame frequencies in the phenomenological description of the soliton dynamics has a simple explanation:

The frequencies (97) describe, in essence, the motion of the soliton centre. Let us discuss their link with the frequencies of small oscillations, see (III). It is easy to write down expressions for the magnetization components in the laboratory coordinate frames:

$`m_x+im_y`$ $`=`$ $`[\mathrm{sin}\theta _0+\mathrm{cos}\theta _0{\displaystyle \underset{\alpha }{}}(u_\alpha +v_\alpha )\mathrm{cos}(m\chi +\stackrel{~}{\omega }_\alpha t)`$
$`+`$ $`i{\displaystyle \underset{\alpha }{}}(u_\alpha v_\alpha )\mathrm{sin}(m\chi +\stackrel{~}{\omega }_\alpha t)]e^{iq\chi +i\mathrm{\Omega }t},`$
$`m_z`$ $`=`$ $`\mathrm{cos}\theta _0\mathrm{sin}\theta _0{\displaystyle \underset{\alpha }{}}\left(u_\alpha +v_\alpha \right)\mathrm{cos}(m\chi +\stackrel{~}{\omega }_\alpha t).`$
One can see that the in–plane components have no simple frequency dependence, including the combined values $`\stackrel{~}{\omega }\pm \mathrm{\Omega }`$. At the same time, the out–of–plane components, $`m_z`$, depend on the frequency $`\stackrel{~}{\omega }`$ in the rotating frame only. In our phenomenological description only the $`z`$–components contain information about the position of the soliton centre $`𝑿`$, determined by the condition $`m_z=1`$. Thus, just the frequency $`\stackrel{~}{\omega }`$, which is calculated above and plotted in all figures, determines the soliton motion.

The solution of the effective equation (95) with the nonzero frequency $`\omega _1`$ allows us to calculate the effective mass of the soliton, $`M=G/\omega _1`$. The value of the gyroconstant is well–known, $`G=4\pi qA/D`$, see Refs. . As a result, the soliton mass is finite for every soliton radius (the local mode with $`\omega _1`$ exists for every $`R`$, see Fig. 3). This fact corresponds qualitatively to the calculations in Ref. . The soliton mass we have obtained numerically from the $`\omega _1`$–data is plotted in Fig. 12 as a function of $`R`$.

In the limit of a large radius soliton the formula (49) yields the dependence

$$M=M_0\left(\frac{R}{\mathrm{\Delta }}_0\right)^3,M_0=\frac{\mathrm{}N_2}{D},R\mathrm{\Delta }_0,$$
(99)
which agrees with the results obtained in Ref.  and our numerical data, see Fig. 12.

In the case of small radius, $`R\mathrm{\Delta }_0`$, the dependence (52) is valid, thus

$$M=M_0\frac{2\omega _0}{\omega _0\mathrm{\Omega }}=\frac{4\pi A}{D^2}r_0^2,R\mathrm{\Delta }_0,$$
(100)
where $`r_0`$ was introduced in Eq. (10). This result corresponds to the one for the Belavin–Polyakov soliton, $`M=4\pi AD^2L^2`$ , because here the characteristic length $`r_0`$ is to be replaced by the system size $`L`$.

## VIII Conclusion

We have studied the magnon–soliton system in the model of the 2d Heisenberg easy–axis ferromagnet. Combining numerical and analytical methods, we have obtained complete results about bound and scattering magnon states in the system. A rich spectrum of truly local modes was found along with bifurcations of these modes with the change of the soliton radius. For the modes with higher azimuthal numbers $`m`$, we have verified the picture of doublets with small splitting for large soliton radii $`R`$. The spectrum changes strongly with decreasing of $`R`$: all modes with $`|m|>1`$ “leave” the region of the discrete spectrum, transforming to quasi–local modes. Such modes could be observed experimentally by soliton magnetic resonance as was done for 1d solitons.

The influence of truly local modes is important for the scattering problem, leading to the bifurcations of the phase shift. As a result, when $`R<1.52\mathrm{\Delta }_0`$, the scattering picture is qualitatively the same as for the isotropic magnet which could be explained in the framework of a simple “centrifugal model”. There is one exception only, $`m=1`$, where the local mode exists for every soliton radius.

Our investigations can be applied to the description of the eigenmodes of magnetic dots. In the present paper we develop the theory of local modes on the soliton background; it could be a good guide for the study of the normal modes in the vortex–state magnetic dots. Our theory is constructed for the soliton in the model (1), where the soliton is stabilized by the internal precession. It is clear that such a model cannot guarantee the quantitative correspondence with the case of vortex–state magnetic dots, where the static soliton structure is stabilized by the magnetic–dipole interaction. We did not consider this type of interaction in the paper, as it is difficult to account for. Nevertheless we believe that the main features of the problem studied above are generic. First of all, we expect the appearance of modes with anomalously small frequencies, e.g. the mode of the translational oscillations of the vortex centre. The non–zero frequency of this mode is caused by the interaction with the boundary only. Secondly, doublets with $`m=\pm n`$, $`n>1`$, should appear, when the radius of the vortex is rather large.

We have used the results on the local modes for the most interesting case of the translational modes with $`|m|=1`$ to describe the soliton motion in the infinite–size magnet: it is possible to identify the soliton mass, which is finite due to the localized soliton structure. In contrast to both the Belavin–Polyakov soliton in the isotropic FM and the magnetic vortex in the easy–plane FM, the soliton motion in our case of the easy–axis FM is similar to the motion of a finite mass charged particle in a magnetic field.

Thus, to the best of our knowledge, the localized precessional soliton in the easy–axis FM is the first example of a 2d topological magnetic soliton with truly particle–like properties.

###### Acknowledgements.

D.D.Sh thanks the University of Bayreuth, where part of this work was performed as a post–doc of the Graduiertenkolleg “Nichtlineare Spektroskopie und Dynamik”, for kind hospitality. Work at Kiev was partially supported by INTAS–97-31311 Grant.


format: text

text:
# 1 Introduction

## 1 Introduction

There is long expectation to find the signal of quark-gluon plasma formation in the high energy experiments of nucleus-nucleus collision. Indeed, not a few of the experimental results suggest the existence of the quark phase , but it seems that the incident energy is a little lower to observe the clean signals even at CERN SPS. The recent experiment at BNL RHIC may give us cleaner and richer signals.

The purpose of the present paper is to give a manageable model for evolution of hadronic matter and quark-gluon plasma, produced in nucleus collisions, and to give some numerical predictions. There are two main ingredients in the model; the evolution equation of the matter and the equation of state for hadrons.

In the evolution equation, assumed is cylindrical symmetry with respect to the collision axis, and used are the averaged quantities over the transverse plane. The averaged evolution equation is already reported in  , and its benefits are in that the parameterization of the unknown initial-conditions becomes very simple and in that the concerned physical quantities, such as size and transverse velocity of the matter, are directly obtained, in contrast to the straightforward calculation of the relativistic fluid equation .

As for the equation of state for hadrons, the compressible bag model    is used, which has nice features as explained in the following. The simplest point-like pion gas model is often used but is lack of degree of freedom, and then it cannot reproduce the observed multiplicity of hadrons as long as natural parameters are employed. On the other hand, if the higher resonances are fully introduced in the point-like model, there appears the limiting temperature and the deconfinement transition to quark phase never takes place. The compressible bag model gives a solution for this dilemma, which gives sufficient degree of freedom and sufficiently high temperature because the mass of the hadrons becomes heavier as they are compressed by surrounding pressure. In this model the phase diagram with baryons is also calculable explicitly for whole region of parameters, and it coincides with the expected shape  .

In Fig. 1, shown are the three steps in the nucleus collisions: At the first step, gluons are exchanged between the nucleons, and the both nucleus are shifted to higher color-charge state. The initial energy density is essentially given at this stage. In the second step, the color flux is expanded in the longitudinal direction, and at proper time $`\tau _01\mathrm{fm}`$, the flux materializes into hadrons, quark-gluons or their mixed phase depending on the energy density. In the third step, the matter expands both in longitudinal and transverse directions, being cooled down to the freeze-out temperature, and then breaks into individual hadrons.

Since the present model treats the third step, in order to show explicit results of the model, we must take some assumptions for the initial energy density of the produced fluid. Therefore the details of the results depends on the assumption, but it is shown that the quark-gluon phase is achieved in RHIC experiments and that the signal of the phase transition, however, is difficult to observe in simple analyses of $`p_t`$ data. The latter is mainly because of the richness of freedom in hadron phase: The richness pushes up the entropy density in hadron phase, and its difference between hadron and quark phases become much smaller than that calculated in simple pion gas model. Therefore the period of mixed phase becomes also smaller, and the typical structure of the first order phase transition is smeared out.

In the next section, the equation of state in the compressible bag model is briefly reviewed on the basis of the results of  . In §3, evolution equations are presented for perfect fluid of hadron or quark-gluon or their mixture and it is revealed that the evolution equations have an expanding solution as well as a scaling solution. And relations to observables are also given. In §4, the results of numerical calculations are given and discussed on the basis of the expanding solution for $`NN`$-collisions and $`AA`$-collisions. Final section is devoted to conclusions and discussion.

## 2 Equation of states

In this section, we briefly summarize the results of  . Let us suppose a mixed gas of $`n`$-kinds of meson enclosed in volume $`V`$ at finite temperature $`T`$. In the compressible bag model, the free energy function $`\widehat{F}`$ of the gas is given by

$`\widehat{F}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}F_f(N_i,V^{},T,M_i(v_i)),`$ (2.1)
$`V^{}`$ $`=`$ $`Vb{\displaystyle \underset{j=1}{\overset{n}{}}}N_jv_j,`$ (2.2)
where $`N_i`$ is a number of $`i`$-th kind of meson and $`v_i`$ is its volume and $`M_i`$ its mass. The constant $`b`$ is a volume exclusion efficiency parameter. The function $`F_f`$ is a free energy function of free (point-like) boson gas. As for mass function $`M_i(v_i)`$, we assume that of MIT bag model:

$`M_i(v_i)=A_iv_i^{1/3}+Bv_i,`$ (2.3)
where $`B`$ is bag constant.

Under the approximation that the average of inverse Lorentz factor of $`i`$-th kind of meson equals to unity ($`\gamma _i^11`$) in the rest frame of the system, basic requirement of the compressible bag model that $`\widehat{F}/v_i=0`$ and the requirement that the the chemical potential of meson should vanish, $`\mu _i=\widehat{F}/N_i=0`$, determine the pressure $`p`$ of the system as a function of $`T`$ by the following equation:

$`p`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}g_iT{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\mathrm{ln}\left\{1\mathrm{exp}[(E_i\mu _i^{}(p,m_i))/T]\right\}},`$ (2.5)
$`E_i=\sqrt{𝒌^2+M_i(p,m_i)^2},`$
$`M_i(p,m_i)`$ $`=`$ $`m_i\left(1+{\displaystyle \frac{3bp}{4B}}\right)\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4},`$ (2.6)
$`\mu _i^{}(p,m_i)`$ $`=`$ $`bv_ip={\displaystyle \frac{bm_ip}{4B}}\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4},`$ (2.7)
where $`g_i`$ is a degeneracy factor of the $`i`$-th kind of meson and $`m_i=4(A_i/3)^{3/4}B^{1/4}`$ is its mass in the vacuum.

Now we take the limit $`n\mathrm{}`$ and replace the sum with the integral

$`{\displaystyle \underset{i=1}{\overset{n}{}}}g_i{\displaystyle 𝑑m\tau (m)},`$ (2.8)
where $`\tau (m)`$ is a level density function. Then (2.5) is written as

$`p`$ $`=`$ $`{\displaystyle 𝑑m\tau (m)T\frac{d^3𝒌}{(2\pi )^3}\mathrm{ln}\left[1e^{(E\mu ^{})/T}\right]},E=\sqrt{𝒌^2+M^2},`$ (2.9)
$`M`$ $`=`$ $`m\left(1+{\displaystyle \frac{3bp}{4B}}\right)\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4},`$ (2.10)
$`\mu ^{}`$ $`=`$ $`{\displaystyle \frac{bmp}{4B}}\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4}.`$ (2.11)
Thus we can calculate $`p`$ as functions of $`T`$ from (2.9) when an appropriate level density function is given. An entropy density $`s`$ is given by

$`s`$ $`=`$ $`{\displaystyle \frac{dp}{dT}}={\displaystyle \frac{C}{D}},`$ (2.12)
$`C`$ $`=`$ $`{\displaystyle \frac{p}{T}}+{\displaystyle 𝑑m\tau (m)\frac{d^3𝒌}{(2\pi )^3}\frac{(E\mu ^{})/T}{\mathrm{exp}[(E\mu ^{})/T]1}},`$ (2.13)
$`D`$ $`=`$ $`1+{\displaystyle 𝑑m\tau (m)\frac{d^3𝒌}{(2\pi )^3}\frac{1}{\mathrm{exp}[(E\mu ^{})/T]1}\frac{d(E\mu ^{})}{dp}},`$ (2.14)
$`{\displaystyle \frac{d(E\mu ^{})}{dp}}`$ $`=`$ $`{\displaystyle \frac{mb}{4B}}\left[1+\left(1+{\displaystyle \frac{3M}{E}}\right){\displaystyle \frac{bp}{4B}}\right]\left(1+{\displaystyle \frac{bp}{B}}\right)^{7/4}.`$ (2.15)
An energy density $`\epsilon `$ is calculated by the following thermodynamic relation

$`\epsilon =Tsp.`$ (2.16)
A number density $`\rho `$ for the hadron with vacuum mass $`m`$ is given by

$`\rho `$ $`=`$ $`{\displaystyle \frac{\rho ^{}}{1+b𝑑m\tau (m)\rho ^{}v}},`$ (2.17)
$`v`$ $`=`$ $`{\displaystyle \frac{m}{4B}}\left(1+{\displaystyle \frac{bp}{B}}\right)^{3/4},`$ (2.18)
$`\rho ^{}`$ $`=`$ $`{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{1}{\mathrm{exp}[(E\mu ^{})/T]1}}.`$ (2.19)
Precisely speaking, the above $`\rho `$ is not the number density but one should multiply it by $`\tau (m)`$ and integrate over some mass range, say $`[m\mathrm{\Delta }m,m+\mathrm{\Delta }m]`$, to get the number density of the hadron with vacuum mass around $`m`$. In practice, however, we use the above $`\rho `$ simply multiplied the degeneracy factor of the hadron, for example 3 for pion, instead of integrating, as the number density. This prescription is justified when narrow width approximation is applicable, which we will assume in this paper where only pion density is considered in numerical calculations.

As for the level density, we adopt the following form

$`\tau (m)=3\delta (mm_\pi )+a\theta (m2m_\pi )m^d\mathrm{exp}(m/T_M),`$ (2.20)
and whose parameters are

$`a=25.46T_M^2,d=3,T_M=(90B/16\pi ^2)^{1/4},`$ (2.21)
according to the Abelian bag model by Kapusta .

## 3 A simple fluid model

### 3.1 Evolution equations

In this subsection, we present evolution equations of perfect fluid of hadron or quark-gluon or their mixture. In the previous paper , hadrons had been assumed to be massless, but in the present paper, hadrons are assumed to obey the equation of states given in the previous section.

As shown in §3 of , transverse flow with cylindrical symmetry obeys the following equations:

$`{\displaystyle \frac{d\stackrel{~}{s}}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{s}}{\tau }},`$ (3.1)
$`{\displaystyle \frac{d\stackrel{~}{\epsilon }}{d\tau }}`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{\epsilon }+\stackrel{~}{p}}{\tau }},`$ (3.2)
where $`\tau `$ is defined by $`(t^2z^2)^{1/2}`$, and has the meaning of the proper time only on the collision axis (z-axis). $`\stackrel{~}{s}`$, $`\stackrel{~}{\epsilon }`$ and $`\stackrel{~}{p}`$ are quantities integrated over the transverse plane;

$`\stackrel{~}{s}`$ $`=`$ $`{\displaystyle d^2rs\gamma _t},\stackrel{~}{\epsilon }={\displaystyle d^2r[(\epsilon +p)\gamma _t^2p]},\stackrel{~}{p}={\displaystyle d^2rp},`$ (3.3)
where $`\gamma _t`$ is the transverse $`\gamma `$-factor of the fluid, and $`\epsilon `$, $`p`$ and $`s`$ are energy density, pressure and entropy density respectively.

Following , let us introduce a radius $`R`$ of the cylindrical fluid and assume the following kinematical relation;

$`{\displaystyle \frac{dR}{d\tau }}=v_t,`$ (3.4)
where $`v_t`$ is the transverse velocity and is related to $`\gamma _t`$ by the relation $`\gamma _t=(1v_t^2)^{1/2}`$. And further let us take a crude approximation,

$`\stackrel{~}{s}`$ $`=`$ $`\pi R^2s\gamma _t,\stackrel{~}{\epsilon }=\pi R^2[(\epsilon +p)\gamma _t^2p],\stackrel{~}{p}=\pi R^2p,`$ (3.5)
and hereafter regard $`s`$, $`\epsilon `$ and $`\gamma _t`$ as the averages over the transverse plane, or as the values at the representative point of $`r`$. All the quantities, then, only depend on $`\tau `$.

Integrating (3.1) and using (3.5), we get the following equation

$`\tau \stackrel{~}{s}=\pi R^2s\gamma _t\tau =\mathrm{const}..`$ (3.6)
Differentiation of this equation with respect to $`\tau `$ yields

$`0`$ $`=`$ $`{\displaystyle \frac{2\dot{R}}{R}}+{\displaystyle \frac{\dot{s}}{s}}+{\displaystyle \frac{\dot{\gamma }_t}{\gamma _t}}+{\displaystyle \frac{1}{\tau }}`$ (3.7)
$`=`$ $`{\displaystyle \frac{2v_t}{R}}+{\displaystyle \frac{\dot{s}}{s}}+{\displaystyle \frac{v_t}{1v_t^2}}\dot{v}_t+{\displaystyle \frac{1}{\tau }},`$ (3.8)
where dot represents differentiation with respect to $`\tau `$.

Before proceeding to derive an evolution equation, we should note the following relations obtained from the thermodynamic relations (2.12) and (2.16),

$`{\displaystyle \frac{\dot{\epsilon }+\dot{p}}{\epsilon +p}}={\displaystyle \frac{\dot{T}}{T}}+{\displaystyle \frac{\dot{s}}{s}},{\displaystyle \frac{\dot{p}}{\epsilon +p}}={\displaystyle \frac{\dot{T}}{T}},{\displaystyle \frac{\dot{\epsilon }}{\epsilon +p}}={\displaystyle \frac{\dot{s}}{s}}.`$ (3.9)
Then using (3.4)$``$(3.9) and (2.16), (3.2) is rewritten as follows;

$`0`$ $`=`$ $`{\displaystyle \frac{\dot{\stackrel{~}{\epsilon }}}{\stackrel{~}{\epsilon }+\stackrel{~}{p}}}+{\displaystyle \frac{1}{\tau }}`$ (3.10)
$`=`$ $`{\displaystyle \frac{2\dot{R}}{R}}+{\displaystyle \frac{\dot{\epsilon }+\dot{p}}{\epsilon +p}}+{\displaystyle \frac{2\dot{\gamma }_t}{\gamma _t}}{\displaystyle \frac{2\dot{R}p}{R(\epsilon +p)\gamma _t^2}}{\displaystyle \frac{\dot{p}}{(\epsilon +p)\gamma _t^2}}+{\displaystyle \frac{1}{\tau }}`$ (3.11)
$`=`$ $`{\displaystyle \frac{2v_t}{R}}+{\displaystyle \frac{\dot{T}}{T}}+{\displaystyle \frac{\dot{s}}{s}}+{\displaystyle \frac{2v_t}{1v_t^2}}\dot{v}_t{\displaystyle \frac{2v_tp}{R(\epsilon +p)\gamma _t^2}}{\displaystyle \frac{\dot{T}}{T\gamma _t^2}}+{\displaystyle \frac{1}{\tau }}`$ (3.12)
$`=`$ $`v_t\left\{v_t{\displaystyle \frac{\dot{T}}{T}}+{\displaystyle \frac{\dot{v}_t}{1v_t^2}}{\displaystyle \frac{2(1v_t^2)p}{R(\epsilon +p)}}\right\},`$ (3.13)
that is,

$`v_t=0,`$ (3.14)
or

$`v_t{\displaystyle \frac{\dot{T}}{T}}+{\displaystyle \frac{\dot{v}_t}{1v_t^2}}{\displaystyle \frac{2(1v_t^2)p}{R(\epsilon +p)}}=0.`$ (3.15)
The equation (3.14) has a scaling solution $`R=R_0`$, $`s=s_0\tau _0/\tau `$, where $`R_0`$ and $`s_0`$ are initial values at $`\tau =\tau _0`$. As for (3.15), it together with (3.8) is recasted as

$`\dot{v}_t`$ $`=`$ $`{\displaystyle \frac{1v_t^2}{1\alpha v_t^2}}\left[{\displaystyle \frac{\alpha v_t}{\tau }}+{\displaystyle \frac{1}{R}}\left\{2\alpha v_t^2+{\displaystyle \frac{2p(1v_t^2)}{\epsilon +p}}\right\}\right],`$ (3.16)
$`{\displaystyle \frac{\dot{T}}{T}}`$ $`=`$ $`{\displaystyle \frac{\alpha }{1\alpha v_t^2}}\left[{\displaystyle \frac{2v_t}{R}}\left\{1+{\displaystyle \frac{p(1v_t^2)}{\epsilon +p}}\right\}+{\displaystyle \frac{1}{\tau }}\right],`$ (3.17)
$`\alpha `$ $``$ $`{\displaystyle \frac{s}{Tds/dT}}.`$ (3.18)
The equations (3.4), (3.16) and (3.17) are a closed set of evolution equations of fluid, which is also applicable for the fluid in quark-gluon phase as well as hadron phase. These equations have an expanding solution if an initial value $`v_t(\tau _0)`$ is non-negative. In the next section, we will numerically analyze on the basis of this expanding solution.

For quark-gluon phase, equations of states are given by

$`p=K_QT^4B,s=4K_QT^3,\epsilon =3K_QT^4+B,`$ (3.19)
with

$`K_Q=\left({\displaystyle \frac{7N_f}{60}}+{\displaystyle \frac{8}{45}}\right)\pi ^2,`$ (3.20)
where $`N_f`$ is a number of quark flavor. The evolution equations (3.16) and (3.17) for quark-gluon fluid  read

$`\dot{v}_t`$ $`=`$ $`{\displaystyle \frac{1v_t^2}{3v_t^2}}\left\{{\displaystyle \frac{v_t}{\tau }}+{\displaystyle \frac{3+v_t^23(1v_t^2)B/K_QT^4}{2R}}\right\},`$ (3.21)
$`{\displaystyle \frac{\dot{T}}{T}}`$ $`=`$ $`{\displaystyle \frac{1}{3v_t^2}}\left\{{\displaystyle \frac{v_t(5v_t^2)}{2R}}+{\displaystyle \frac{1}{\tau }}\right\}.`$ (3.22)
For mixed phase, equations of states are given by

$`p`$ $`=`$ $`p_c=p_Q(T_c)=p_H(T_c),`$ (3.23)
$`s`$ $`=`$ $`xs_Q(T_c)+(1x)s_H(T_c),`$ (3.24)
$`\epsilon `$ $`=`$ $`x\epsilon _Q(T_c)+(1x)\epsilon _H(T_c),`$ (3.25)
where $`T_c`$ is critical temperature and the subscript $`H`$ ($`Q`$) denotes the quantity for hadron (quark-gluon), and $`x`$ is a fraction of quark-gluon. Corresponding equations to (3.16) and (3.17) are given by

$`\dot{v}_t`$ $`=`$ $`{\displaystyle \frac{2p}{\epsilon +p}}{\displaystyle \frac{(1v_t^2)^2}{R}},`$ (3.26)
$`\dot{x}`$ $`=`$ $`{\displaystyle \frac{\dot{s}}{s_Q(T_c)s_H(T_c)}}={\displaystyle \frac{s}{s_Q(T_c)s_H(T_c)}}\left[{\displaystyle \frac{2v_t}{R}}\left\{1+{\displaystyle \frac{p(1v_t^2)}{\epsilon +p}}\right\}+{\displaystyle \frac{1}{\tau }}\right].`$ (3.27)
Thus we have obtained evolution equations for all phases. Those are (3.4), (3.16) and (3.17) for hadron phase, (3.4), (3.21) and (3.22) for quark-gluon phase and (3.4), (3.26) and (3.27) for mixed phase. Throughout the evolution, $`\tau \stackrel{~}{s}`$ is conserved, as is seen in (3.6), so that it can be used as a check in solving the evolution equations numerically.

### 3.2 Relation to observables

Now that we have the evolution equations, we can solve them and get $`T`$ ($`x`$ for mixed phase), $`R`$ and $`\gamma _t`$ as functions of $`\tau `$ if some initial values $`T_0`$, ($`x_0`$ for mixed phase) and $`R_0`$ are given and the initial condition $`\gamma _t(\tau _0)=1`$ is assumed. The final condition of the evolution is set by the freeze-out of hadrons from hadronic fluid. If the freeze-out temperature $`T_f`$ is given, we can get $`\tau _f`$ from the solution,

$`T_f=T(\tau _f),`$ (3.28)
where the subscript $`f`$ means freeze-out. Then we get the values of the radius $`R`$ and $`\gamma _t`$ at the hadronization;

$`R_f=R(\tau _f),\gamma _{tf}=\gamma _t(\tau _f).`$ (3.29)
These quantities are related to experimental observables.

(i) Charged multiplicity

The produced hadron density in rapidity space is given by

$`{\displaystyle \frac{dN}{dy}}=\tau _f\pi R_f^2\rho _f\gamma _{tf},`$ (3.30)
where $`\rho _f`$ is the number density at $`\tau =\tau _f`$. If we assume that all the produced hadrons are pions and put $`N_{\mathrm{ch}}=(2/3)N`$, we get the charged particle density,

$`{\displaystyle \frac{dN_{\mathrm{ch}}}{dy}}={\displaystyle \frac{2}{3}}\tau _f\pi R_f^2\widehat{\rho }_{\pi f}\gamma _{tf},`$ (3.31)
where $`\widehat{\rho }_\pi `$ denotes effective pion number density defined by

$`\widehat{\rho }_\pi `$ $`=`$ $`\rho _\pi +N_M\rho _M,`$ (3.32)
$`\rho _\pi `$ $`=`$ $`={\displaystyle \frac{\rho _\pi ^{}}{1+𝑑m\tau (m)\rho ^{}v}},`$ (3.33)
$`\rho _M`$ $`=`$ $`{\displaystyle \frac{\rho _M^{}}{1+𝑑m\tau (m)\rho ^{}v}},`$ (3.34)
$`\rho _\pi ^{}`$ $`=`$ $`3{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{1}{\mathrm{exp}[(E_\pi \mu ^{})/T]1}},`$ (3.35)
$`\rho _M^{}`$ $`=`$ $`{\displaystyle _{2m_\pi }^{\mathrm{}}}𝑑m\tau (m){\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{1}{\mathrm{exp}[(E\mu ^{})/T]1}}.`$ (3.36)
In (3.32), the first term corresponds to pole term and the second term corresponds to continuum part of $`\tau (m)`$ in (2.20). The quantity $`N_M`$ is an average multiplicity to pions from the continuum of meson resonances.

(ii) Transverse momentum and transverse energy

The average transverse momentum of pion in the center of mass system is approximately given by

$`p_t=\sqrt{m_t^2\gamma _t^2m_\pi ^2},m_t^{}=\sqrt{p_t^2+m_\pi ^2},`$ (3.37)
where $`p_t^{}`$ ($`m_t^{}`$) is a transverse momentum (mass) of pion in the rest frame of the fluid. The average transverse mass of pion in the rest frame of the fluid is calculated as an average of the one for directly produced pion and the one for the contribution from pions as decay products of continuum of meson resonances;

$`m_t^2`$ $`=`$ $`{\displaystyle \frac{\rho _{\pi f}m_t^2(\pi _d)+N_M\rho _{Mf}m_t^2(\pi _c)}{\rho _{\pi f}+N_M\rho _{Mf}}}.`$ (3.38)
The quantity $`m_t^2(\pi _c)`$ is estimated on the basis of random walk approximation as

$`m_t^2(\pi _c)=m_t^2(M)/N_M.`$ (3.39)
Then it follows from (3.33), (3.34), (3.38) and (3.39)

$`m_t^2`$ $`=`$ $`{\displaystyle \frac{\rho _{\pi f}^{}m_t^2(\pi _d)+\rho _{Mf}^{}m_t^2(M)}{\rho _{\pi f}^{}+N_M\rho _{Mf}^{}}},`$ (3.40)
$`m_t^2(\pi _d)`$ $`=`$ $`{\displaystyle \frac{{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{k_t^2+M_\pi ^2}{\mathrm{exp}[(E_\pi \mu ^{})/T]1}}}{{\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{1}{\mathrm{exp}[(E_\pi \mu ^{})/T]1}}}},`$ (3.41)
$`m_t^2(M)`$ $`=`$ $`{\displaystyle \frac{{\displaystyle _{2m_\pi }^{\mathrm{}}}𝑑m\tau (m){\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{k_t^2+M^2}{\mathrm{exp}[(E_\pi \mu ^{})/T]1}}}{{\displaystyle _{2m_\pi }^{\mathrm{}}}𝑑m\tau (m){\displaystyle \frac{d^3𝒌}{(2\pi )^3}\frac{1}{\mathrm{exp}[(E_\pi \mu ^{})/T]1}}}},`$ (3.42)
where $`k_t`$ denotes a transverse momentum variable.

Transverse energy per unit rapidity interval is simply given by

$`{\displaystyle \frac{dE_t}{dy}}=\sqrt{m_t^2}\gamma _t{\displaystyle \frac{dN}{dy}}.`$ (3.43)
(iii) Rapidity interval and total multiplicity

Since the evolution of fluid does not take place yet at $`\tau =\tau _0`$, total energy of the fluid $`E_{\mathrm{tot}}`$ in the center of mass system is calculated as

$`E_{\mathrm{tot}}`$ $`=`$ $`{\displaystyle _{\eta _{\mathrm{max}}}^{\eta _{\mathrm{max}}}}\tau _0\epsilon _0\pi R_0^2\mathrm{cosh}\eta d\eta =2\tau _0\epsilon _0\pi R_0^2\mathrm{sinh}\eta _{\mathrm{max}}.`$ (3.44)
Then the rapidity interval of the fluid $`2\eta _{\mathrm{max}}`$ is determined when $`E_{\mathrm{tot}}`$ and initial temperature $`T_0`$ ($`x_0`$ for mixed-phase start) are given, that is,

$`2\eta _{\mathrm{max}}=2\mathrm{sinh}^1\left({\displaystyle \frac{E_{\mathrm{tot}}}{2\tau _0\epsilon _0\pi R_0^2}}\right).`$ (3.45)
Total multiplicity $`N_{\mathrm{ch}}`$ is approximately given by

$`N_{\mathrm{ch}}2\eta _{\mathrm{max}}{\displaystyle \frac{dN_{\mathrm{ch}}}{dy}}.`$ (3.46)
## 4 Numerical calculations

### 4.1 $`NN`$-collision

In order to start numerical calculations, we have to fix bag constant $`B`$, volume exclusion parameter $`b`$, the number of quark flavors $`N_f`$, average multiplicity of pions from continuum $`N_M`$, initial proper time $`\tau _0`$, initial temperature $`T_0`$ ($`x_0`$ for mixed-phase start), initial radius $`R_0`$ and freeze-out temperature $`T_f`$. First of all, we fix as

$`N_f=2.`$ (4.1)
The initial radius $`R_0`$ for $`NN`$-collision is set to the experimental value of proton charge radius

$`R_0(NN)=0.82\mathrm{fm}.`$ (4.2)
The parameter $`b`$ is determined by the relation

$`bv_N(p=0)={\displaystyle \frac{bm_N}{4B}}={\displaystyle \frac{4\pi }{3}}(0.82\mathrm{fm})^3,`$ (4.3)
when the bag constant is given. We adopt the following values

$`B^{1/4}=0.25\mathrm{GeV},b=5.00.`$ (4.4)
Then critical values at the deconfinement transition are given by

$`T_c=0.196\mathrm{GeV},\epsilon _H(T_c)=0.874\mathrm{GeV}/\mathrm{fm}^3,\epsilon _Q(T_c)=2.85\mathrm{GeV}/\mathrm{fm}^3.`$ (4.5)
For the purpose of determining $`T_f`$ and $`\tau _0`$, we use the fit to experimental data of $`NN`$-collisions as inputs ;

$`T_f`$ $`=`$ $`0.135\pm 0.003\mathrm{GeV},`$ (4.6)
$`\gamma _{tf}`$ $`=`$ $`1+(0.010\pm 0.003)\mathrm{log}(E_{\mathrm{cm}}/E_1)\theta (E_{\mathrm{cm}}E_1),`$ (4.8)
$`E_1=12\pm 8\mathrm{GeV},`$
$`{\displaystyle \frac{dN_{\mathrm{ch}}}{dy}}`$ $`=`$ $`(1.31\pm 0.02)+(0.20\pm 0.03)\mathrm{log}(E_{\mathrm{cm}}/E_2)`$
$`+(0.062\pm 0.007)\mathrm{log}^2(E_{\mathrm{cm}}/E_2),`$
$`E_2=12\mathrm{GeV},`$
where $`E_{\mathrm{cm}}`$ is the initial energy in the center of mass system. We choose the central value $`0.135\mathrm{GeV}`$ of (4.6) as $`T_f`$. Since the fluid immediately freeze-out at $`E_{\mathrm{cm}}=12\mathrm{GeV}`$ as seen from (4.8), the value of $`\tau _0`$ is fixed to reproduce the value $`dN_{\mathrm{ch}}/dy=1.31`$ at $`E_{\mathrm{cm}}=12\mathrm{GeV}`$ as seen from (4.1);

$`{\displaystyle \frac{dN_{\mathrm{ch}}}{dy}}={\displaystyle \frac{2}{3}}\tau _0\pi R_0^2\widehat{\rho }_{\pi f}=1.31.`$ (4.11)
Now varying $`T_0`$ (or $`x_0`$), we can calculate $`p_t`$ as a function of $`dN_{\mathrm{ch}}/dy`$ with giving the value of $`N_M`$. The results are shown in Fig. 2. The value of $`N_M`$ is searched for the range $`2.53.0`$. The value of $`p_t`$ is almost determined by $`m_t^2(\pi _c)=m_t^2(M)/N_M`$ since the value of $`N_M\rho _{Mf}`$ is almost an order of magnitude larger than that of $`\rho _{\pi f}`$. (See (3.37)$``$(3.39)) Therefore the calculated values of $`p_t`$ is sensitive to the value of $`N_M`$. Comparing with the experimental data  we choose as

$`N_M=2.7,\tau _0=1.29\mathrm{fm}.`$ (4.12)
The response to $`B`$ is shown in Fig. 3. For small bag constant, $`p_t`$ becomes large slightly fast as $`dN_{\mathrm{ch}}/dy`$ becomes large. This is recognized as a consequence of a larger acceleration for smaller $`B`$ as is seen from (3.21).

In Fig. 4, the initial temperature $`T_0`$ is shown as a function of $`dN_{\mathrm{ch}}/dy`$. For $`B^{1/4}=0.25\mathrm{GeV}`$ the critical temperature is $`T_c=0.196\mathrm{GeV}`$ and the region $`4.2<dN_{\mathrm{ch}}/dy<11`$ corresponds to the mixed-phase start. For a larger (smaller) value of $`B`$, $`T_c`$ becomes higher (lower) and the region of the mixed-phase start becomes wider (narrower). The value $`T_c=0.237\mathrm{GeV}`$ for $`B^{1/4}=0.30\mathrm{GeV}`$ is somewhat larger than the value suggested by lattice QCD calculations so that we have adopted $`B^{1/4}=0.25\mathrm{GeV}`$ in (4.4). Comparing with Fig. 3, we can see a slight structure of entrance to mixed phase while entrance to quark-gluon phase does not show a remarkable structure.

An initial energy density $`\epsilon _0`$ has a fundamental importance to prescribe an initial state of the fluid and determine its evolution in our picture as explained in §1. The values of $`\epsilon _0`$ can be obtained as a function of $`E_{\mathrm{cm}}`$ in such a way that the calculated values using (3.31) reproduce the predicted values by (4.1). However, (4.1) is only applicable for $`E_{\mathrm{cm}}1.8\mathrm{TeV}`$ since (4.1) was obtained as a result of fit to experimental data up to Tevatron region. In order to extrapolate above Tevatron region we assume the following form of $`\epsilon _0`$

$`\epsilon _0(NN)=\epsilon _{00}\left({\displaystyle \frac{E_{\mathrm{cm}}}{E_2}}\right)^\delta ,`$ (4.13)
and determine $`\epsilon _{00}`$ and $`\delta `$ to reproduce (4.1) below Tevatron region. The result is given by

$`\epsilon _{00}=0.184\mathrm{GeV}/\mathrm{fm}^3,\delta =0.294.`$ (4.14)
Using this result $`dN_{\mathrm{ch}}/dy`$ is shown as a function of $`E_{\mathrm{cm}}`$ in Fig. 5 where solid line show a consequence of (4.13) and (4.14) while dashed line show (4.1). We can see both lines are coincide below Tevatron region.

### 4.2 $`AA`$-collision

First we set the initial radius $`R_0`$ for $`AA`$-collision as

$`R_0(AA)=1.2A^{1/3}\mathrm{fm}.`$ (4.15)
In order to incorporate $`AA`$ collisions into numerical calculations, we need a relation to combine $`\epsilon _0(AA)`$ with $`\epsilon _0(NN)`$. For this purpose we assume that

$`\epsilon _0(AA)=A^\beta \epsilon _0(NN).`$ (4.16)
In the above, both sides should be compared at equal $`E_{\mathrm{cm}}`$ which means energy per nucleon in the center of mass system for $`AA`$-collision.

In Fig. 6 the average transverse energy per unit rapidity interval $`dE_t/dy`$ is shown as a function of $`A`$ for $`E_{\mathrm{cm}}=20\mathrm{GeV}/\mathrm{A}`$ and the values of $`\beta =0.400.45`$. Comparing with experimental data , where the initial energies are $`E_{\mathrm{cm}}=19.4\mathrm{GeV}/\mathrm{A}`$ ($`E_{\mathrm{Lab}}=200\mathrm{GeV}/\mathrm{A}`$) and $`E_{\mathrm{cm}}=17.3\mathrm{GeV}/\mathrm{A}`$ ($`E_{\mathrm{Lab}}=158\mathrm{GeV}/\mathrm{A}`$), we choose

$`\beta =0.43.`$ (4.17)
Having fixed all parameters we can now calculate various observables as functions of $`E_{\mathrm{cm}}`$ and $`A`$.

In Fig. 7 $`dE_t/dy`$ is shown as a function of $`A`$ for $`E_{\mathrm{cm}}=20\mathrm{GeV}/\mathrm{A}`$, $`200\mathrm{GeV}/\mathrm{A}`$, $`2\mathrm{TeV}/\mathrm{A}`$. $`dE_t/dy`$ grows almost linearly with $`A`$. In Fig. 8 the initial temperature $`T_0`$ is shown as a function of $`dN_{\mathrm{ch}}/dy`$ for $`A=1,30,200`$. These results correspond to the experiments for various incident energies. For the case of $`A=1`$ ($`R_0=R_0(NN)=0.82\mathrm{fm}`$) $`4.2<dN_{\mathrm{ch}}/dy<1.1\times 10`$ corresponds to mixed-phase start. For the case of $`A=30`$ ($`R_0=3.73\mathrm{fm}`$) $`8.6\times 10<dN_{\mathrm{ch}}/dy<2.3\times 10^2`$ corresponds to mixed-phase start. As for the case of $`A=200`$ ($`R_0=7.02\mathrm{fm}`$) $`dN_{\mathrm{ch}}/dy<8.2\times 10^2`$ corresponds to mixed-phase start. The behaviour of $`T_0`$ in the case of no phase transition is also shown by thin lines for comparison.

In Fig. 9 the average transverse momentum $`p_t`$ is shown as a function of $`dN_{\mathrm{ch}}/dy`$ for $`A=1,30,200`$. These results also correspond to the experiments for various incident energies. The behaviour of $`p_t`$ in the case of no phase transition is also shown by thin lines for comparison. It is clearly seen that the fluid is not accelerated so much in the mixed phase. The reason is that the acceleration in the mixed phase given by the right hand side of (3.26) has no $`1/\tau `$-term and an extra factor $`(1v_t^2)`$ in $`1/R`$-term so that it takes smaller values than the acceleration in hadron/quark-gluon phase given by (3.16)/(3.21). This is a common feature to the case of massless hadron, which is supposed to be pion, discussed in the previous paper . The difference lies in the fact that the interval of mixed-phase start in the present case is narrower than the one in the case of a single massless hadron so that the structure is less prominent especially for $`A=30,200`$. This is because the degree of freedom in the hadron phase is larger for the present case which takes into account the contribution from continuum of meson resonances. However it should be noticed that the entrance to quark-gluon phase can be seen in this figure plotted in the log scale of $`dN_{\mathrm{ch}}/dy`$ in contrast to Fig. 3 plotted in the linear scale.

In Fig. 10 $`p_t`$ is shown as a function of $`dN_{\mathrm{ch}}/dy`$ for $`E_{\mathrm{cm}}=20\mathrm{GeV},200\mathrm{GeV},2\mathrm{TeV}`$ with varying $`A`$ from 1 to 240. To this variation the response of $`p_t`$ is mild and does not become so large even for $`E_{\mathrm{cm}}=2\mathrm{TeV}`$. In Fig. 11 the final radius $`R_f`$ is shown as a function of $`dN_{\mathrm{ch}}/dy`$ for $`E_{\mathrm{cm}}=20\mathrm{GeV},200\mathrm{GeV},2\mathrm{TeV}`$ with varying $`A`$ from 1 to 240. The experimental data are taken from . The calculated values does not differ so much from the observed values. It should be noticed that $`R_f`$ is rather insensitive to the variation of $`E_{\mathrm{cm}}`$. As for the response to the variation of $`\beta `$, $`R_f`$ becomes larger for smaller $`\beta `$.

Finally in Fig. 12 the initial temperature is shown as a function of $`A`$ for $`E_{\mathrm{cm}}=20\mathrm{GeV},`$$`200\mathrm{GeV},2\mathrm{TeV}`$. For $`E_{\mathrm{cm}}=20\mathrm{GeV}`$ the entrance to the mixed phase occurs at $`A=27`$ and the entrance to the quark-gluon phase does not occur even at highest $`A`$. For $`E_{\mathrm{cm}}=200\mathrm{GeV}`$ the entrance to the mixed phase occurs at $`A=6`$ and the entrance to the quark-gluon phase occurs at $`A=87`$. As for $`E_{\mathrm{cm}}=2\mathrm{TeV}`$ the entrance to the mixed phase occurs at $`A=2`$ and the entrance to the quark-gluon phase occurs at $`A=18`$. Thus quark-gluon phase will be realized at BNL RHIC energies while it is not at CERN SPS energies.

## 5 Conclusions and discussion

We have shown that our simplified fluid model easily calculates observables in nucleon-nucleon and nucleus-nucleus collisions at high energy region. As for the equation of state for hadron phase, we have adopted the compressible bag model. This model can satisfy both requirements that one should take into account an infinite degrees of freedom of hadrons and that deconfinement transition should be take place, which usually contradict mutually in models other than compressible bag model. For an evolution of the fluid, we have assumed the cylindrical symmetry with respect to the collision axis and solved the averaged quantities over the transverse plane.

In the previous paper  of the averaged evolution equation, we used the equation of state of massless pion for hadron phase. Both calculations give the result that agrees well in terms of qualitative analysis. However, the present approach produces $`(\rho _{\pi f}+N_M\rho _{Mf})/\rho _{\pi f}`$ times of pions in comparison with the previous one. This gives a reasonable value for $`\tau _0`$ by the condition of (4.11). On the other hand, the ratio $`\epsilon _Q(T_c)/\epsilon _H(T_c)`$ becomes smaller because of the increase of hadronic degrees of freedom.

The increase of the average transverse momentum is brought about by a transverse flow of the fluid. As is seen from (3.16) the pressure slope $`p/R`$ in the transverse direction mainly causes the transverse flow. In the mixed phase the effect from the corresponding term became smaller as shown by equation (3.21). This gives a structure of an average transverse momentum as a function of $`dN_{\mathrm{ch}}/dy`$. In the hadron and the quark-gluon phase the ratio of the pressure to the energy density increases rapidly with the increase of the energy density. The ratio in both phases reaches $`1/3`$ in high energy density limit, so that the acceleration is generated similarly in both phases as shown in Fig. 9. In the case of the quark-gluon phase start the acceleration becomes weaker because of passing through the mixed phase. As for $`AA`$-collisions the acceleration becomes weaker with an increase of $`A`$ because of a decrease of the pressure slope in the transverse direction.

As shown in Figs. 8 and 9 $`dN_{\mathrm{ch}}/dy`$ is not a reasonable variable for plotting the observables in $`AA`$-collisions except for the transverse source radius. In our model the equations of (3.16) and (4.16) generate the $`A`$ dependence. We have to find a scaling variable in $`AA`$-collision.

We have assumed the relation between the initial energy densities in $`NN`$\- and $`AA`$\- collisions by equation (4.16). This power-law assumption is similar to the relations given by the dilute-gas and multiple-collision models . Finally we predict that the quark-gluon plasma can be observed for average multiplicity events in $`AA`$-collision with $`A87`$ at RHIC energies. The numerical calculation results are shown in Fig. 12.