text:
# Soliton self-modulation of the turbulence amplitude and plasma rotation

## 1 Introduction

The ion temperature gradient turbulence, considered a major source of anomalous transport in the tokamak plasma, is characterized by the coexistence of irregular patterns (randomly fluctuating field) and intermittent robust cuasi-coherent structures. In closely related fluid models (for exemple, in the physics of atmosphere) modulational instabilities are known to produce solitary structures on the envelope of the fluctuating field. In the case of the tokamak, this can be particularly important since a poloidally nonuniform amplitude of the turbulence generates nonuniform transport rates. For a sufficiently high nonuniformity the torque arising via the mechanism initially mentioned by Stringer  may overcame the rotation damping due to the poloidal magnetic pumping. In this work we show that this can indeed be the case by proving that the uniform poloidal envelope of the ITG turbulence is an unstable state. It appears that the ITG turbulence has intrinsic resources to generate poloidal rotation via a combination of mechanisms which is not directly related to the Reynolds stress, inverse cascade, direct-ion loss, or other classical sources of rotation. We provide an essentially mathematical explanation of the instability of the solution consisting of poloidally uniform envelope of the turbulence. Based on the geometrico-algebraic method of solving integrable nonlinear differential equations on periodic domains we invoke an existing result, that any perturbation removes the degeneracies due to coincident eigenvalues in the main spectrum of the Lax operator, thus changing the topology of the hyperelliptic Riemann surface that provides the solution. The perturbed solution separates exponentially (in function space) from the initial one (uniform envelope) and this yields an exponential growth of the poloidal nonuniformity. The magnitude of the torque may be comparable to the poloidal magnetic damping.

There are many works in plasma physics related to the soliton dynamics , , , , , etc. In a recent paper the focusing solution has been used to study the formation of coherent motion or intermittent patterns (*streamers*) . Our work is basically a mathematical approach, using well developed technics related to the Inverse Scattering Transform method for periodic domains. However, more physical analysis may become possible combining the knowledge of the spectral properties of the ion turbulence with the nonlinear stability properties.

This work consists of several parts that may appear as developing separately: the derivation of the ion equation (barotropic equation), the multiple space-time scales analysis, the solution of the Nonlinear Schrodinger Equation and the stability of the solution, the torque arising from the Stringer mechanism and the possibility of the rotation. For being self-contained the paper also includes a review of the multiple space-time analysis (from the atmospheric physics applications) and of the geometric-algebraic method of integration of the Nonlinear Schrodinger Equation. Even if these parts could be found in basic works (see the references) they are included both for clarity and for their extreme importance for further applications related to our problem or independent of this.

## 2 The slab model of the ion mode instability

We consider cylindrical geometry with circular magnetic surfaces. Locally the model can be reduced to a slab geometry with $`(x,y)`$ cartezian coordinates replacing respectively the radius and poloidal angle coordinates $`(r,\theta )`$. At equilibrium the plasma parameters are constant on the magnetic surfaces. The effects of the toroidicity and of the particle drifts are not included instead the nonlinearity related to the ion polarization drift is fully retained. The plasma model consists of the continuity equations and the equations of motion for electrons and ions:

$$\frac{n_\alpha }{t}+\left(n_\alpha 𝐯_\alpha \right)=0$$
(1)
$$m_\alpha n_\alpha \left(\frac{𝐯_\alpha }{t}+\left(𝐯_\alpha \right)𝐯_\alpha \right)=p_\alpha \pi _\alpha +e_\alpha n_\alpha \left(𝐄+𝐯_\alpha \times 𝐁\right)+𝐑_\alpha $$
(2)
where $`\alpha =e,i`$. The friction forces $`𝐑_e=𝐑_i=n\left|e\right|𝐉_{}/\sigma _{}`$ , which are important for the parallel electron momentum balance, vanish for infinite plasma conductivity, which we will assume. The collisional viscosity $`\pi _{e,i}`$ will be neglected as well. However we will need to include it later when we will consider the balance of the forces contributing to the poloidal rotation. The electron and ion temperatures are considered constant along the magnetic lines $`_{}T_{e,i}=0`$. The equilibrium quantities are perturbed by the wave potential $`\varphi `$: $`n=n_0+\stackrel{~}{n}`$. A sheared poloidal plasma rotation is included, and we later will make explicit the corresponding part in the potential, $`\varphi _0`$.

The momentum conservation equations are used to determine the perpendicular velocities of the electrons and ions. The parallel momentum conservation equation for electrons, in the absence of *dissipation* or *drifts* gives the adiabatic distribution of the density fluctuation. The velocities are introduced in the continuity equations to find the dynamical equations for the density and electric potential.

From the equations of motion for the *ions* the velocities are obtained in the form:

$$𝐯_i=𝐯_i=𝐯_{dia,i}+𝐯_E+𝐯_{pol,i}$$

where the ion diamagnetic velocity is

$$𝐯_{dia,i}=\frac{1}{n_i}\frac{1}{m_i\mathrm{\Omega }_i}\widehat{𝐧}\times p_i$$

The versor of the magnetic field is $`\widehat{𝐧}`$; the versors along the transversal coordinate axis $`(x,y)`$ will be noted $`(\widehat{𝐞}_x,\widehat{𝐞}_y)`$. The ion-polarization velocity is:

$`𝐯_{pol,i}`$ $`=`$ $`{\displaystyle \frac{1}{\left|e\right|n_iB}}n_im_i\left({\displaystyle \frac{𝐯_E}{t}}+\left(𝐯_E\right)𝐯_E\right)\times \widehat{𝐧}`$
$`=`$ $`{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\left({\displaystyle \frac{}{t}}+\left(𝐯_E_{}\right)\right)_{}\varphi `$
Using the notation $`\frac{d}{dt}\frac{}{t}+\left(𝐯_E_{}\right)`$ the perpendicular ion velocity can be written

$`𝐯_{,i}`$ $`=`$ $`𝐯_{dia,i}+𝐯_{pol,i}+{\displaystyle \frac{\varphi \times \widehat{𝐧}}{B}}`$
$`=`$ $`{\displaystyle \frac{T_i}{\left|e\right|B}}{\displaystyle \frac{1}{n_i}}{\displaystyle \frac{dn_i}{dr}}\widehat{𝐞}_y{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}{\displaystyle \frac{d}{dt}}\left(_{}\varphi \right)+{\displaystyle \frac{\varphi \times \widehat{𝐧}}{B}}`$
The equation for the velocity of the *electrons* is

$$𝐯_e=𝐯_{,e}+𝐯_{,e}$$

$`𝐯_e`$ $`=`$ $`𝐯_E+𝐯_{dia,e}`$
$`=`$ $`{\displaystyle \frac{\varphi \times \widehat{𝐧}}{B}}+{\displaystyle \frac{1}{n}}{\displaystyle \frac{1}{m_e\mathrm{\Omega }_e}}\left(\widehat{𝐧}\times p\right)`$
$`=`$ $`{\displaystyle \frac{\varphi \times \widehat{𝐧}}{B}}+{\displaystyle \frac{T_e}{\left|e\right|B}}{\displaystyle \frac{1}{n_0}}{\displaystyle \frac{dn_0}{dr}}\widehat{𝐞}_y`$ (6)
We assume neutrality $`n_e=n_i=n`$ and introduce the expressions of the velocities in the continuity equations for ions and for electrons.

The ion continuity equation is

$$\frac{n}{t}+\left(n𝐯_{,i}\right)=0$$

The electron continuity equation is

$$\frac{n}{t}+_{}\left(n𝐯_{,e}\right)+_{}\left(nv_e\right)=0$$

The equations are substracted , to obtain

$`_{}\left(nv_e\right)`$
$`{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}_{}\left(n{\displaystyle \frac{d}{dt}}_{}\varphi \right)`$
$`+_{}\left(n𝐯_{dia,i}\right)_{}\left(n𝐯_{dia,e}\right)`$
$`=`$ $`0`$
From the last term in the left we get:

$`{\displaystyle \frac{dn_0}{dx}}\widehat{𝐞}_x𝐯_{dia,i}+\left(_{}\stackrel{~}{n}\right)𝐯_{dia,i}`$ $`=`$
$`=`$ $`{\displaystyle \frac{T_i}{T_e}}\left({\displaystyle \frac{T_e}{\left|e\right|B}}{\displaystyle \frac{1}{n_0}}{\displaystyle \frac{dn_0}{dx}}\right){\displaystyle \frac{\stackrel{~}{n}}{y}}`$
Including the similar term for the electrons, we obtain

$`_{}\left(nv_e\right)`$
$`{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}_{}\left(n{\displaystyle \frac{d}{dt}}_{}\varphi \right)`$
$`{\displaystyle \frac{T_i}{T_e}}\left({\displaystyle \frac{T_e}{\left|e\right|B}}{\displaystyle \frac{1}{n_0}}{\displaystyle \frac{dn_0}{dx}}\right){\displaystyle \frac{\stackrel{~}{n}}{y}}\left({\displaystyle \frac{T_e}{\left|e\right|B}}{\displaystyle \frac{1}{n_0}}{\displaystyle \frac{dn_0}{dx}}\right){\displaystyle \frac{\stackrel{~}{n}}{y}}`$
$`=`$ $`0`$
or:

$$\left(1+\frac{T_i}{T_e}\right)\left(\frac{T_e}{\left|e\right|B}\frac{1}{n_0}\frac{dn_0}{dx}\right)\frac{\stackrel{~}{n}}{y}+\frac{1}{B\mathrm{\Omega }_i}_{}\left(n\frac{d}{dt}_{}\varphi \right)+_{}\left(nv_e\right)=0$$
(7)
From the continuity equation for electrons

$$\frac{n}{t}+_{}\left(n𝐯_{}\right)+_{}\left(n𝐯_{}\right)=0$$
(8)
we obtain

$`{\displaystyle \frac{\stackrel{~}{n}}{t}}+{\displaystyle \frac{1}{B}}{\displaystyle \frac{\varphi }{y}}{\displaystyle \frac{dn_0}{dx}}++𝐯_{dia,e}_{}\stackrel{~}{n}+V_0{\displaystyle \frac{\stackrel{~}{n}}{y}}+`$
$`+n_0\left(_{}𝐯_{}\right)+v_{}_{}\stackrel{~}{n}`$
$`=`$ $`0`$
where the seed poloidal velocity is $`𝐕_0\left(x\right)=\frac{d\varphi _0\left(x\right)}{Bdx}\widehat{𝐞}_y`$. The parallel momenttum balance gives the parallel electron velocity

$$𝐯_e=\frac{\sigma _{}}{e^2}_{}\left(\left|e\right|\varphi +T_e\mathrm{ln}\left(n/n_0\right)\right)$$
(10)
In the absence of friction $`\left(\sigma \mathrm{}\right)`$ and of particle drifts the electron response is adiabatic

$$\frac{\stackrel{~}{n}}{n_0}=\frac{\left|e\right|\varphi }{T_e}$$

and the potential is determined from Eq.(7). To develop separately the ion-polarization drift term, we introduce the notation:

$`W`$ $``$ $`{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}_{}\left[\left(n_0+\stackrel{~}{n}\right)\left({\displaystyle \frac{}{t}}+𝐯_E_{}\right)_{}\varphi \right]`$
$`=`$ $`{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}_{}\left[\left(n_0+\stackrel{~}{n}\right)𝐈\right]`$
where we make explicit the electric potential $`\varphi _0`$ associated to the initial plasma poloidal rotation, $`V_0\widehat{𝐞}_y`$.

$$𝐈\left(\frac{}{t}+\left(V_0\widehat{𝐞}_y+\stackrel{~}{𝐯}\right)_{}\right)_{}\left(\varphi _0+\stackrel{~}{\varphi }\right)$$

We have

$`𝐈`$ $`=`$ $`{\displaystyle \frac{}{t}}\left(_{}\varphi _0\right)+`$
$`+{\displaystyle \frac{}{t}}\left(_{}\stackrel{~}{\varphi }\right)+`$
$`+\left[V_0{\displaystyle \frac{}{y}}+\left(\stackrel{~}{𝐯}_{}\right)\right]\left(_{}\varphi _0+_{}\stackrel{~}{\varphi }\right)`$
or

$`𝐈`$ $`=`$ $`{\displaystyle \frac{}{t}}\left(\widehat{𝐞}_xBV_0\right)+`$
$`+{\displaystyle \frac{}{t}}\left(_{}\stackrel{~}{\varphi }\right)+`$
$`+\left[V_0{\displaystyle \frac{}{y}}+\left(\stackrel{~}{𝐯}_{}\right)\right]\left(_{}\varphi _0+_{}\stackrel{~}{\varphi }\right)`$
$`=`$ $`{\displaystyle \frac{}{t}}\left(\widehat{𝐞}_xBV_0\right)+V_0B{\displaystyle \frac{V_0}{y}}\widehat{𝐞}_x+`$
$`+{\displaystyle \frac{}{t}}\left(_{}\stackrel{~}{\varphi }\right)+\left(\stackrel{~}{𝐯}_{}\right)BV_0\widehat{𝐞}_x+`$
$`+V_0{\displaystyle \frac{}{y}}\left(_{}\stackrel{~}{\varphi }\right)+\left(\stackrel{~}{𝐯}_{}\right)\left(_{}\stackrel{~}{\varphi }\right)`$
with the relations

$`_{}\stackrel{~}{\varphi }`$ $`=`$ $`B\left(\stackrel{~}{𝐯}\times \widehat{𝐧}\right)`$
$`=`$ $`\left(B\stackrel{~}{v}_y\right)\widehat{𝐞}_x+\left(B\stackrel{~}{v}_x\right)\widehat{𝐞}_y`$
$$\stackrel{~}{v}_y=\frac{1}{B}\frac{\stackrel{~}{\varphi }}{x}\text{and}\stackrel{~}{v}_x=\frac{1}{B}\frac{\stackrel{~}{\varphi }}{y}$$

After very simple calculations we obtain:

$`{\displaystyle \frac{1}{B}}I_x`$ $`=`$ $`{\displaystyle \frac{V_0}{t}}+\stackrel{~}{v}_y{\displaystyle \frac{V_0}{y}}+V_0{\displaystyle \frac{V_0}{y}}`$
$`+{\displaystyle \frac{\stackrel{~}{v}_y}{t}}+\stackrel{~}{v}_x{\displaystyle \frac{V_0}{x}}+V_0{\displaystyle \frac{\stackrel{~}{v}_y}{y}}`$
$`+\stackrel{~}{v}_x{\displaystyle \frac{\stackrel{~}{v}_y}{x}}+\stackrel{~}{v}_y{\displaystyle \frac{\stackrel{~}{v}_y}{y}}`$
and

$`{\displaystyle \frac{1}{B}}I_y`$ $`=`$ $`{\displaystyle \frac{\stackrel{~}{v}_x}{t}}V_0{\displaystyle \frac{\stackrel{~}{v}_x}{y}}`$
$`\stackrel{~}{v}_x{\displaystyle \frac{\stackrel{~}{v}_x}{x}}\stackrel{~}{v}_y{\displaystyle \frac{\stackrel{~}{v}_x}{y}}`$
It will be useful to calculate the derivaties of these quantities

$`{\displaystyle \frac{1}{B}}{\displaystyle \frac{I_x}{x}}`$ $`=`$ $`{\displaystyle \frac{}{x}}{\displaystyle \frac{V_0}{t}}+{\displaystyle \frac{\stackrel{~}{v}_y}{x}}{\displaystyle \frac{V_0}{y}}+\stackrel{~}{v}_y{\displaystyle \frac{^2V_0}{yx}}+{\displaystyle \frac{V_0}{x}}{\displaystyle \frac{V_0}{y}}+V_0{\displaystyle \frac{^2V_0}{xy}}`$
$`+{\displaystyle \frac{\stackrel{~}{v}_y}{t}}{\displaystyle \frac{V_0}{x}}+\stackrel{~}{v}_x{\displaystyle \frac{^2V_0}{x^2}}+{\displaystyle \frac{}{t}}{\displaystyle \frac{\stackrel{~}{v}_y}{x}}+`$
$`+{\displaystyle \frac{V_0}{x}}{\displaystyle \frac{\stackrel{~}{v}_y}{y}}+V_0{\displaystyle \frac{^2\stackrel{~}{v}_y}{xy}}+`$
$`+{\displaystyle \frac{\stackrel{~}{v}_x}{x}}{\displaystyle \frac{\stackrel{~}{v}_x}{y}}+\stackrel{~}{v}_x{\displaystyle \frac{^2\stackrel{~}{v}_y}{x^2}}+`$
$`+{\displaystyle \frac{\stackrel{~}{v}_y}{x}}{\displaystyle \frac{\stackrel{~}{v}_y}{y}}+\stackrel{~}{v}_y{\displaystyle \frac{^2\stackrel{~}{v}_y}{xy}}`$
and

$`{\displaystyle \frac{1}{B}}{\displaystyle \frac{I_y}{y}}`$ $`=`$ $`{\displaystyle \frac{V_0}{y}}{\displaystyle \frac{\stackrel{~}{v}_x}{y}}`$
$`{\displaystyle \frac{}{t}}{\displaystyle \frac{\stackrel{~}{v}_x}{y}}V_0{\displaystyle \frac{^2\stackrel{~}{v}_y}{y^2}}`$
$`{\displaystyle \frac{\stackrel{~}{v}_x}{y}}{\displaystyle \frac{\stackrel{~}{v}_x}{x}}\stackrel{~}{v}_x{\displaystyle \frac{^2\stackrel{~}{v}_x}{xy}}`$
$`{\displaystyle \frac{\stackrel{~}{v}_y}{y}}{\displaystyle \frac{\stackrel{~}{v}_x}{y}}\stackrel{~}{v}_y{\displaystyle \frac{^2\stackrel{~}{v}_x}{y^2}}`$
The quantity denoted by $`W`$ takes the form

$`W`$ $`=`$ $`{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\left({\displaystyle \frac{dn_0}{dx}}\right)I_x+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}n_0\left({\displaystyle \frac{I_x}{x}}\right)+`$
$`+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\left({\displaystyle \frac{\stackrel{~}{n}}{x}}\right)I_x+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\stackrel{~}{n}\left({\displaystyle \frac{I_x}{x}}\right)+`$
$`+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}n_0\left({\displaystyle \frac{I_y}{y}}\right)+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\left({\displaystyle \frac{\stackrel{~}{n}}{y}}\right)I_y+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\stackrel{~}{n}\left({\displaystyle \frac{I_y}{y}}\right)`$
### 2.1  <br>The mode evolution in a fixed plasma rotation profile

We will assume that the mode evolves initially without perturbing the equilibrium profiles, in particular the seed poloidal rotation. This allows us to simplify the expressions above, taking:

$`{\displaystyle \frac{V_0}{y}}`$ $`=`$ $`0`$
$`{\displaystyle \frac{V_0}{t}}`$ $`=`$ $`0`$
Then the first lines in the formulas Eqs.(2), (2), (2) are zero. Let us consider in the expression of $`W`$ the part $`W_0`$ which does not contain the fluctuating density $`\stackrel{~}{n}`$. Writting

$$W=W_0+\stackrel{~}{W}$$

we have

$`W_0`$ $``$ $`{\displaystyle \frac{1}{\mathrm{\Omega }_i}}\left({\displaystyle \frac{dn_0}{dx}}\right)\left[\stackrel{~}{v}_x{\displaystyle \frac{V_0}{x}}+V_0{\displaystyle \frac{\stackrel{~}{v}_y}{y}}+\stackrel{~}{v}_x{\displaystyle \frac{\stackrel{~}{v}_y}{x}}+\stackrel{~}{v}_y{\displaystyle \frac{\stackrel{~}{v}_y}{y}}+{\displaystyle \frac{\stackrel{~}{v}_y}{t}}\right]+`$
$`+{\displaystyle \frac{1}{\mathrm{\Omega }_i}}n_0_{}\left(I_x\widehat{𝐞}_x+I_y\widehat{𝐞}_y\right)`$
and

$$\stackrel{~}{W}=\frac{1}{B\mathrm{\Omega }_i}\left(\frac{\stackrel{~}{n}}{x}\right)I_x+\frac{1}{B\mathrm{\Omega }_i}\stackrel{~}{n}\left(\frac{I_x}{x}\right)+\frac{1}{B\mathrm{\Omega }_i}\left(\frac{\stackrel{~}{n}}{y}\right)I_y+\frac{1}{B\mathrm{\Omega }_i}\stackrel{~}{n}\left(\frac{I_y}{y}\right)$$

Replacing the perturbed velocity by the perturbed potential, writting all terms and summing, we get:

$`W_0`$ $`=`$ $`\left({\displaystyle \frac{}{t}}+V_0{\displaystyle \frac{}{y}}\right){\displaystyle \frac{1}{B}}{\displaystyle \frac{\stackrel{~}{\varphi }}{x}}+\left({\displaystyle \frac{1}{B}}{\displaystyle \frac{\stackrel{~}{\varphi }}{y}}\right){\displaystyle \frac{dV_0}{dx}}+\left({\displaystyle \frac{_{}\varphi \times \widehat{𝐧}}{B}}_{}\right){\displaystyle \frac{\stackrel{~}{\varphi }}{x}}+`$
$`+{\displaystyle \frac{1}{\mathrm{\Omega }_i}}n_0\{{\displaystyle \frac{1}{B}}V_0{\displaystyle \frac{}{y}}_{}^2\stackrel{~}{\varphi }+`$
$`+{\displaystyle \frac{1}{B}}\left({\displaystyle \frac{_{}\varphi \times \widehat{𝐧}}{B}}_{}\right)_{}^2\stackrel{~}{\varphi }`$
$`{\displaystyle \frac{1}{B}}{\displaystyle \frac{\stackrel{~}{\varphi }}{y}}{\displaystyle \frac{d^2V_0}{dx^2}}+`$
$`+{\displaystyle \frac{1}{B}}{\displaystyle \frac{}{t}}\left(_{}^2\stackrel{~}{\varphi }\right)\}`$
Collecting all what we have at this moment the ion continuity equation (7) becomes:

$`_{}\left(nv_e\right)+`$
$`\left(1+{\displaystyle \frac{T_i}{T_e}}\right)\left({\displaystyle \frac{T_e}{\left|e\right|B}}{\displaystyle \frac{1}{n_0}}{\displaystyle \frac{dn_0}{dx}}\right){\displaystyle \frac{\stackrel{~}{n}}{y}}+`$
$`+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}n_0\left\{\left({\displaystyle \frac{}{t}}+V_0{\displaystyle \frac{}{y}}\right)_{}^2\stackrel{~}{\varphi }V_0^{^{\prime \prime }}{\displaystyle \frac{\stackrel{~}{\varphi }}{y}}+\left({\displaystyle \frac{_{}\varphi \times \widehat{𝐧}}{B}}_{}\right)_{}^2\stackrel{~}{\varphi }\right\}`$
$`+{\displaystyle \frac{1}{B\mathrm{\Omega }_i}}\left({\displaystyle \frac{dn_0}{dx}}\right)\left[\left({\displaystyle \frac{}{t}}+V_0{\displaystyle \frac{}{y}}\right)\stackrel{~}{v}_y+\stackrel{~}{v}_x{\displaystyle \frac{\stackrel{~}{v}_y}{x}}+\stackrel{~}{v}_y{\displaystyle \frac{\stackrel{~}{v}_y}{y}}+\stackrel{~}{v}_x{\displaystyle \frac{V_0}{x}}\right]`$
$`+\stackrel{~}{W}+\left\{\text{terms of the first line in the expressions of }I_x\text{ and derivatives}\right\}`$
$`=`$ $`0`$
In the above equation (which is exact) we shall make the following approximations:

* neglect the term containing the parallel electron velocity, assuming infinite electric conductivity;
* neglect the term which contains $`dn_0/dx`$ since it is in the ratio $`k:1/L_n`$ with the other terms, and we consider

$$kL_n1$$
* neglect $`\stackrel{~}{W}`$ ; these terms are in the ratio $`\stackrel{~}{n}/n_0`$ with the terms which are retained;
* neglect of the terms in the first lines of the expressions for $`I_x`$ and $`dI_x/dx`$, $`dI_y/dy`$. (These are the terms in the curly brakets, the last line). As explained above, we assume that the mode evolves in a background of fixed rotation profile, $`V_0\left(x\right)`$.

The resulting equation is

$`B\mathrm{\Omega }_i\left(1+{\displaystyle \frac{T_i}{T_e}}\right)\left({\displaystyle \frac{T_e}{\left|e\right|B}}{\displaystyle \frac{1}{n_0}}{\displaystyle \frac{dn_0}{dx}}\right){\displaystyle \frac{}{y}}\left({\displaystyle \frac{\stackrel{~}{n}}{n_0}}\right)`$
$`+\left({\displaystyle \frac{}{t}}+V_0{\displaystyle \frac{}{y}}\right)_{}^2\stackrel{~}{\varphi }V_0^{^{\prime \prime }}{\displaystyle \frac{\stackrel{~}{\varphi }}{y}}+\left({\displaystyle \frac{_{}\varphi \times \widehat{𝐧}}{B}}_{}\right)_{}^2\stackrel{~}{\varphi }`$
$`=`$ $`0`$
and, replace the adiabatic form of the density perturbation

$$\frac{\stackrel{~}{n}}{n_0}=\frac{\left|e\right|\stackrel{~}{\varphi }}{T_e}$$

We define

$$\beta B\mathrm{\Omega }_i\left(1+\frac{T_i}{T_e}\right)\left(\frac{T_e}{\left|e\right|B}\frac{1}{n_0}\frac{dn_0}{dx}\right)\frac{\left|e\right|}{T_e}=\mathrm{\Omega }_i\left(1+\frac{T_i}{T_e}\right)\frac{1}{L_n}$$

and obtain

$$\beta \frac{\stackrel{~}{\varphi }}{y}+\left(\frac{}{t}+V_0\frac{}{y}\right)_{}^2\stackrel{~}{\varphi }V_0^{^{\prime \prime }}\frac{\stackrel{~}{\varphi }}{y}+\left(\frac{_{}\varphi \times \widehat{𝐧}}{B}_{}\right)_{}^2\stackrel{~}{\varphi }=0$$

which is the barotropic equation.

### 2.2 Nondimensional form of the equation

We consider that the ion mode extends in the spatial ($`x`$) direction over a length $`L`$. A typical value for the sheared poloidal rotation is noted $`U_0`$. We make the replacements

$$yLy$$

$$tt\frac{L}{U_0}$$

$$\stackrel{~}{\varphi }\varphi \frac{T_e}{\left|e\right|}$$

$$V_0UU_0$$

such that from now on $`y,t,\varphi `$ and $`U`$ are nondimensional quantities. We also change the radial coordinate into a dimensionless variable

$$xLx$$

and rewrite the equation

$`\left(\mathrm{\Omega }_i\left(1+{\displaystyle \frac{T_i}{T_e}}\right){\displaystyle \frac{1}{L_n}}{\displaystyle \frac{L^2}{U_0}}\right){\displaystyle \frac{\varphi }{y}}+`$
$`+\left({\displaystyle \frac{}{t}}+U{\displaystyle \frac{}{y}}\right)_{}^2\varphi `$
$`{\displaystyle \frac{d^2U}{dx^2}}{\displaystyle \frac{\varphi }{y}}+`$
$`+\left({\displaystyle \frac{1}{U_0}}{\displaystyle \frac{T_e}{\left|e\right|}}{\displaystyle \frac{1}{B}}{\displaystyle \frac{1}{L}}\right)\left[\left(_{}\varphi \times \widehat{𝐧}\right)_{}\right]_{}^2\varphi `$
$`=`$ $`0`$
The coefficients are

$$\beta ^{}\mathrm{\Omega }_i\left(1+\frac{T_i}{T_e}\right)\frac{1}{L_n}\frac{L^2}{U_0}\text{non-dimensional}$$

$$\epsilon \frac{1}{U_0}\frac{T_e}{\left|e\right|}\frac{1}{B}\frac{1}{L}\text{non-dimensional}$$

For an order of magnitude, $`\epsilon `$ is the ratio of the diamagnetic electron velocity to the rotation velocity $`U_0`$ multiplied by the ratio of the density gradient length to the length of the spatial domain. This quantity, $`\epsilon `$ is in general smaller than unity.

The quantity $`\beta ^{}`$ is the ratio of the ion cyclotron frequency to the inverse of the time required to cross the spatial domain with the typical flow velocity. Since the later $`\left(U_0/L\right)`$ involves macroscopic quantities this ratio can be large. It is multiplied by the ratio of the spatial length to the density gradient length (these quantities can be comparable and the ratio not too different of unity).

We change the notations eliminating the primes. The equation becomes

$$\beta \frac{\varphi }{y}+\left(\frac{}{t}+U\frac{}{y}\right)_{}^2\varphi \frac{d^2U}{dx^2}\frac{\varphi }{y}+\epsilon \left[\left(_{}\varphi \times \widehat{𝐧}\right)_{}\right]_{}^2\varphi =0$$

This is the *barotropic* equation, known in the physics of the atmosphere.

## 3 Multiple time and space scale analysis of the ion mode equation reduced to the barotropic form

The quasigeostrophic barotropic atmospheric model leads to the barotropic equation (see Horton  where the quasigeostrophic approximation is explained in relation with the Ertel’s theorem). This equation has been studied in the atmosphere science by means of the *multiple space and time scales method*. Our aim is similar, *i.e.* to obtain an equation for the envelope of the fluctuating field $`\varphi `$ in order to study the possible poloidally nonuniform profiles of the turbulence averaged amplitude. The connection between the original equation and the final slow-time and large-spatial scales equations is not simple : much numerical work is required in order to connect the input physical parameters of the plasma with the formal coefficients in the final equations. The detailed analytical work is presented in Ref() and the necessary formulas are reproduced here for convenience. To simplify the expressions the following notation is introduced $`[\varphi ,^2\varphi ]\left[\left(_{}\varphi \times \widehat{𝐧}\right)_{}\right]_{}^2\varphi `$ and the equation can be written

$$\left(\frac{}{t}+U\left(x\right)\frac{}{y}\right)^2\varphi +\left(\beta \frac{d^2U\left(x\right)}{dx^2}\right)\frac{\varphi }{y}=\epsilon [\varphi ,^2\varphi ]$$

The analysis on multiple scales starts by introducing the variables on scales separated by the parameter $`\epsilon `$:

$$T_1=\epsilon t,T_2=\epsilon ^2t$$

$$Y_1=\epsilon y,Y_2=\epsilon ^2y$$

This gives for the derivatives

$$\frac{}{t}\frac{}{t}+\epsilon \frac{}{T_1}+\epsilon ^2\frac{}{T_2}$$
(15)
$$\frac{}{y}\frac{}{y}+\epsilon \frac{}{Y_1}+\epsilon ^2\frac{}{Y_2}$$
(16)
The solution is adopted to be of the form

$$\varphi =\varphi ^{\left(1\right)}+\epsilon \varphi ^{\left(2\right)}+\epsilon ^2\varphi ^{\left(3\right)}+\mathrm{}$$
(17)
We denote the linear part of the operator in the equation by

$$L\left(\frac{}{t}+U\left(x\right)\frac{}{y}\right)^2+\left(\beta \frac{d^2U\left(x\right)}{dx^2}\right)\frac{}{y}$$

Substituting the Eq.(17) and the forms of the derivation operators Eqs.(15, 16) in the equation of motion we obtain from the equality between the coefficients of the powers of the variable $`\epsilon `$:

$$L\varphi ^{\left(1\right)}=0$$
(18)
$`L\varphi ^{\left(2\right)}`$ $`=`$
$`=`$ $`\left({\displaystyle \frac{}{T_1}}+U{\displaystyle \frac{}{Y_1}}\right)^2\varphi ^{\left(1\right)}+\left(\beta U^{\prime \prime }\right){\displaystyle \frac{\varphi ^{\left(1\right)}}{Y_1}}`$
$`2\left({\displaystyle \frac{}{t}}+U{\displaystyle \frac{}{y}}\right){\displaystyle \frac{^2\varphi ^{\left(1\right)}}{yY_1}}`$
$`[\varphi ^{\left(1\right)},^2\varphi ^{\left(1\right)}]`$
$`L\varphi ^{\left(3\right)}`$ $`=`$
$`=`$ $`\left({\displaystyle \frac{}{T_2}}+U{\displaystyle \frac{}{Y_2}}\right)^2\varphi ^{\left(1\right)}\left(\beta U^{\prime \prime }\right){\displaystyle \frac{\varphi ^{\left(1\right)}}{Y_2}}`$
$`2\left({\displaystyle \frac{}{t}}+U{\displaystyle \frac{}{y}}\right){\displaystyle \frac{^2\varphi ^{\left(1\right)}}{yY_2}}`$
$`\left({\displaystyle \frac{}{t}}+U{\displaystyle \frac{}{y}}\right){\displaystyle \frac{^2\varphi ^{\left(1\right)}}{Y_1^2}}`$
$`2\left({\displaystyle \frac{}{T_1}}+U{\displaystyle \frac{}{Y_1}}\right){\displaystyle \frac{^2\varphi ^{\left(1\right)}}{yY_1}}`$
$`\left({\displaystyle \frac{}{T_1}}+U{\displaystyle \frac{}{Y_1}}\right)^2\varphi ^{\left(2\right)}`$
$`2\left({\displaystyle \frac{}{t}}+U{\displaystyle \frac{}{y}}\right){\displaystyle \frac{^2\varphi ^{\left(2\right)}}{yY_1}}\left(\beta U^{\prime \prime }\right){\displaystyle \frac{\varphi ^{\left(2\right)}}{Y_1}}`$
$`{\displaystyle \frac{\varphi ^{\left(2\right)}}{Y_1}}{\displaystyle \frac{}{x}}^2\varphi ^{\left(1\right)}+{\displaystyle \frac{\varphi ^{\left(1\right)}}{x}}{\displaystyle \frac{}{Y_1}}^2\varphi ^{\left(1\right)}`$
$`2\left({\displaystyle \frac{\varphi ^{\left(1\right)}}{y}}{\displaystyle \frac{^3\varphi ^{\left(1\right)}}{yxY_1}}{\displaystyle \frac{\varphi ^{\left(1\right)}}{x}}{\displaystyle \frac{^3\varphi ^{\left(1\right)}}{y^2Y_1}}\right)`$
$`[\varphi ^{\left(1\right)},^2\varphi ^{\left(2\right)}][\varphi ^{\left(2\right)},^2\varphi ^{\left(1\right)}]`$
The solution is represented on the space and time scales.

### 3.1 The enveloppe equation

The solution of the linear part of the operator, Eq.(18) is adopted as the superposition of two propagating waves

$`\varphi ^{\left(1\right)}`$ $`=`$ $`A_1(T_1,T_2,Y_1,Y_2)\phi _1\left(x\right)\mathrm{exp}\left(ik_1yi\omega _1t\right)`$
$`+A_2(T_1,T_2,Y_1,Y_2)\phi _2\left(x\right)\mathrm{exp}\left(ik_2yi\omega _2t\right)`$
$`+c.c.`$
The two functions $`\phi _n`$, $`n=1,2`$, are solutions of the equation

$$\frac{d^2\phi _n}{dx^2}\left(k_n^2\frac{\beta U^{\prime \prime }}{Uc_n}\right)\phi _n=0$$
(22)
with

$$\phi _n\left(0\right)=\phi _n\left(L\right)=0$$

In this formula, $`(0,L)`$ are the limits of the spatial domain in the minor radius direction where the flow $`U\left(x\right)0`$ ; $`c_n=\omega _n/k_n`$ are the phase velocities. To find the solutions of this equation we have to use numerical methods. This is a Schrodinger-type equation where the potential depends on the energy, (since the phase velocity $`c_n`$ depends on the wavenumber $`k_n`$). We must use numerical methods suitable for Sturm-Liouville equations and in general the solution exists for only particular values of the parameters. We choose to fix the frequencies $`\omega _n`$ and consider $`k_n`$ as the parameter to be determined. The problem is complicated by the possible occurence of singularities in the potential due to resonances between the phase velocity and the fixed zonal flow. We will avoid any resonance since the physical content of the processes consisting of the direct energy transfer between the fluctuations and the flow represents a different physical problem and requires separate consideration. In the present work we examine the effect of the turbulence on the poloidal rotation as being due to the nonuniform diffusion rates and Stringer mechanism.

The parameters in the Eqs.(22) are $`\beta `$, $`\epsilon `$ and the sheared flow velocity $`U(x)`$. We take for $`U\left(x\right)`$ a symmetric Gaussian-like profile shifted in amplitude and retain its maximum, half-width and asymptotic value (*i.e.* at $`x=0`$ or $`L`$) as parameters. The two frequencies have fixed values during the eigenvalue search but they must be considered free parameters as well. This builds up a large space of parameters which should be sampled. For some point in this space the eigenvalue problems Eqs.(22) are solved and the phase velocities compared with the profile of $`U\left(x\right)`$. The presence of resonances renders the set of parameters useless. Since (as will result from the formuals below) there are also other possible resonances, the determination of a useful solution $`(k_1,k_2)`$ is difficult and requires many trials.

The numerical methods that have been used in solving Eqs.(22) belong to two classes: boundary value integration methods and respectively shooting methods. When the sign of the “potential” in the equation is negative everywhere on $`(0,L)`$ , which occurs for short wavelengths and smooth fluid velocity profile the boundary conditions cannot be fulfilled except for asymtotically, *i.e.* assuming arbitrary small nonzero values at these limits, since the solution decays exponentially. However most of the boundary value integrators we have tried gave amplitudes of magnitudes comparable to the boundary value, which essentially are homotopic to the trivial vanishing solution (the equation is homogeneous) and render useless this approach. We have to use a shooting method accepting the difficulty that the initial value of the parameter to be determined (eigenvalue) must be placed not far from the final value. We have used the specialized package SLEIGN2 interactively. For the investigation of ranges of parameters most of the calculations have been performed using the routine D02KEF of NAG. For identical conditions the two numerical codes gave the same results within the accepted tolerance.

Using the routine D02KEF we divide the integration range $`(0,L)`$ into four intervals corresponding roughly to a different behaviour of the potential. This should simplify the task of integrator. A single matching point is assumed, usually at the centre. The tolerence is $`10^8`$ which however does not exclude fake convergence in case of an initialisation of the eigenvalue very far from the true solution..

We notice from Eqs.(22) that the potential can easily be dominated by high shear $`U^{\prime \prime }`$ especially for smaller wavenumbers $`k_n`$. Actually we are more interested to investigate regimes where only mild effects of the seed sheared rotation arise during the process of self-modulation since, as explained before, the possible evolution to the nonuniformity in the poloidal direction will be the source of the torque on the plasma. Avoiding regimes of $`U^{\prime \prime }`$-dominated potentials we assume low shear and take the Gaussian profile with $`U_{\mathrm{max}}U_{\mathrm{min}}=2\times 10^4`$ ($`m/s`$) represented in Figure (1).

The following set of physical parameters has been used: minor radius $`a=1m`$, magnetic field $`B_T=3.5T`$, electron temperature $`T_e=1KeV`$, ion temperature $`T_i=1KeV`$, density $`n=10^{20}m^3`$, $`L_n=10m`$. The parameters in the barotropic equation are then: $`\beta =33.526`$ and $`\epsilon =0.28571`$. The two frequencies are $`\omega _1=85\times 10^2s^1`$ and $`\omega _2=85\times 10^4s^1`$. It is assumed that the space region of significant magnitude of $`U\left(x\right)`$ is $`L=0.1m`$. The following two eigenvalues are obtained, normalised to $`L^1`$: $`k_1=0.5525`$ and $`k_2=3.2029`$. The potentials in the Schrodinger-like equations and the eigenfunctions are shown in Figure (2).

The two eigenfunctions $`\phi _{1,2}`$ have an important role in the following steps of the calculation. Substituting the Eq.(3.1) in the barotropic equation and expanding the operators, we obtain:

$`L\varphi ^{\left(2\right)}`$ $`=`$ $`{\displaystyle \underset{n=1}{\overset{2}{}}}G_{1n}\mathrm{exp}\left(ik_nyi\omega _nt\right)`$
$`{\displaystyle \underset{n=1}{\overset{2}{}}}ig_{2n}A_n^2\mathrm{exp}\left[i2\left(k_ny\omega _nt\right)\right]`$
$`g_3A_1A_2\mathrm{exp}\left[i\left(k_{12}y\omega _{12}t\right)\right]`$
$`ig_4A_1A_2^{}\mathrm{exp}\left[i\left(\alpha _{12}y\sigma _{12}t\right)\right]`$
$`+c.c.`$
where

$`k_{12}`$ $`=`$ $`k_1+k_2`$
$`\omega _{12}`$ $`=`$ $`\omega _1+\omega _2`$
We have $`k_{12}=3.755`$, $`\omega =8.585`$.

$`\alpha _{12}`$ $`=`$ $`k_1k_2`$
$`\sigma _{12}`$ $`=`$ $`\omega _1\omega _2`$
or: $`\alpha _{12}=2.650`$, $`\sigma _{12}=8.415`$.

$`G_{1n}`$ $`=`$ $`\left({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n\right){\displaystyle \frac{A_n}{T_1}}+`$
$`+\left[U\left({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n\right)+\left(\beta U^{\prime \prime }\right)\phi _n2k_n\left(Uc_n\right)\right]{\displaystyle \frac{A_n}{Y_1}}`$
$$g_{2n}=k_n\left(\phi _n\frac{d}{dx}\frac{d\phi _n}{dx}\right)\frac{d^2\phi _n}{dx^2}$$

$`g_3`$ $`=`$ $`\left(k_1\phi _1{\displaystyle \frac{d}{dx}}k_2{\displaystyle \frac{d\phi _1}{dx}}\right)\left({\displaystyle \frac{d^2\phi _2}{dx^2}}k_2^2\phi _2\right)`$
$`+\left(k_2\phi _2{\displaystyle \frac{d}{dx}}k_1{\displaystyle \frac{d\phi _2}{dx}}\right)\left({\displaystyle \frac{d^2\phi _1}{dx^2}}k_1^2\phi _1\right)`$
$`g_4`$ $`=`$ $`\left(k_1\phi _1{\displaystyle \frac{d}{dx}}+k_2{\displaystyle \frac{d\phi _1}{dx}}\right)\left({\displaystyle \frac{d^2\phi _2}{dx^2}}k_2^2\phi _2\right)`$
$`+\left(k_2\phi _2{\displaystyle \frac{d}{dx}}+k_1{\displaystyle \frac{d\phi _2}{dx}}\right)\left({\displaystyle \frac{d^2\phi _1}{dx^2}}k_1^2\phi _1\right)`$
These functions can be obtained after calculating the eigenfunctions $`\phi _{1,2}`$ and are represented in Figure (3).

Various solutions to the equation (3.1) are possible, taking the space-time dependence similar to one of the terms which compose the right hand side.

The second term suggests solutions of the type

$$\varphi _{2n}^{\left(2\right)}=W_{2n}A_n^2\mathrm{exp}\left[i2\left(k_ny\omega _nt\right)\right]+c.c.$$

The third term suggets solutions of the type

$$\varphi _3^{\left(2\right)}=W_3A_1A_2\mathrm{exp}\left[i\left(k_{12}y\omega _{12}t\right)\right]+c.c.$$

The fourth term gives

$$\varphi _4^{\left(2\right)}=W_4A_1A_2^{}\mathrm{exp}\left[i\left(\alpha _{12}y\sigma _{12}t\right)\right]+c.c.$$

In the formulas above the functions $`W_{2n}\left(x\right)`$ are solutions of the equations:

$`{\displaystyle \frac{d^2W_{2n}}{dx^2}}\left(4k_n^2{\displaystyle \frac{\beta U^{\prime \prime }}{Uc_n}}\right)W_{2n}`$ $`=`$ $`{\displaystyle \frac{g_{2n}}{2\left(Uk_n\omega _n\right)}}`$
$`\text{with }W_{2n}\left(0\right)`$ $`=`$ $`W_{2n}\left(L\right)=0`$
$`{\displaystyle \frac{d^2W_3}{dx^2}}\left(k_{12}^2{\displaystyle \frac{k_{12}\left(\beta U^{\prime \prime }\right)}{Uk_{12}\omega _{12}}}\right)W_3`$ $`=`$ $`{\displaystyle \frac{g_3}{Uk_{12}\omega _{12}}}`$
$`\text{with }W_3\left(0\right)`$ $`=`$ $`W_3\left(L\right)=0`$
$`{\displaystyle \frac{d^2W_4}{dx^2}}\left(\alpha _{12}^2{\displaystyle \frac{\alpha _{12}\left(\beta U^{\prime \prime }\right)}{U\alpha _{12}\sigma _{12}}}\right)W_4`$ $`=`$ $`{\displaystyle \frac{g_4}{U\alpha _{12}\sigma _{12}}}`$
$`\text{with }W_4\left(0\right)`$ $`=`$ $`W_4\left(L\right)=0`$
These equations are integrated numerically with a boundary value integrator, a finite difference method. The results are shown in Figures (4), (5) and (6). In the following figures are explained also the functions $`W_{1n}`$ which will be introduced in Eq.(3.1).

The first term gives a solution of the following form

$$\varphi _{1n}^{\left(2\right)}=\phi _{1n}^{\left(2\right)}\mathrm{exp}\left(ik_nyi\omega _nt\right)+c.c.$$

where the function $`\phi _{1n}^{\left(2\right)}`$ satisfies the equation

$`{\displaystyle \frac{d^2\phi _{1n}^{\left(2\right)}}{dx^2}}\left(k_n^2{\displaystyle \frac{\beta U^{\prime \prime }}{Uc_n}}\right)\phi _{1n}^{\left(2\right)}`$ $`=`$
$`=`$ $`{\displaystyle \frac{i}{k_n\left(Uc_n\right)}}\{({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n){\displaystyle \frac{A_n}{T_1}}`$
$`+[U({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n)+(\beta U^{\prime \prime })\phi _n`$
$`2k_n^2(Uc_n)]{\displaystyle \frac{A_n}{Y_1}}\}`$
This equation serves to generate a condition on the two amplitudes $`A_{1,2}`$.

When the left hand side of the equation is multiplied by the function $`\phi _n`$ and integrated between the limits on the $`x`$ domain, it gives zero, due to the boundary conditions. The same must then be true for the right hand side. This gives a “solubility condition”, on the space time scales slower at order $`1`$:

$$\frac{A_n}{T_1}+c_{gn}\frac{A_n}{Y_1}=0$$

where

$$c_{gn}=c_n+2_0^Lk_n^2\phi _n^2𝑑x\mathrm{}_0^L\frac{\beta U^{\prime \prime }}{\left(Uc_n\right)^2}\phi _n^2𝑑x$$

Expressed in units of length $`L`$ and units of time $`t_0=10^4`$, the *group velocities* are : $`c_{g1}=0.30328`$ and $`c_{g2}=3.9908`$.

Then it is clear that the amplitude $`A_n`$ propagates at the speed $`c_{gn}`$ :

$$A_n=A_n(Y_1c_{gn}T_1,Y_2,T_2)$$

Now we introduce a new function

$$\phi _{1n}^{\left(2\right)}=iW_{1n}\frac{A_n}{Y_1}$$

The equation satisfied by $`W_{1n}`$ is

$`{\displaystyle \frac{d^2W_{1n}}{dx^2}}\left(k_n^2{\displaystyle \frac{\beta U^{\prime \prime }}{uc_n}}\right)W_{1n}`$
$`=`$ $`{\displaystyle \frac{1}{k_n\left(Uc_n\right)}}\left[\left(c_nc_{gn}\right)\left({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n\right)2k_n^2\left(Uc_n\right)\phi _n\right]`$
with the bounday conditions

$$W_{1n}\left(0\right)=W_{1n}\left(L\right)=0$$

At this order $`\left(2\right)`$ there are the following solutions

$`\varphi ^{\left(2\right)}`$ $`=`$ $`\varphi _{11}^{\left(2\right)}+\varphi _{12}^{\left(2\right)}+\varphi _{21}^{\left(2\right)}+\varphi _{22}^{\left(2\right)}`$
$`+\varphi _3^{\left(2\right)}+\varphi _4^{\left(2\right)}+`$
$`+\mathrm{\Phi }(x,T_1,Y_1)`$
The last function will be determined by the condition of solubility at the next order on the space-time scales ($`O\left(\epsilon ^2\right)`$). To obtain the evolution equations for $`\mathrm{\Phi }`$ and $`A_n`$, we use the expressions for $`\varphi ^{\left(1\right)}`$ Eq.(3.1) and $`\varphi ^{\left(2\right)}`$, Eq.(3.1) in the right hand side of the equation (3), and all inhomogeneous terms are obtained. In additions there are terms which are not dependent on $`y`$ and $`t`$, and these terms should be zero. This leads to the solubility condition relating the function $`\mathrm{\Phi }`$ (correction to the mean flow due to the presence of the wave packets) to the amplitude $`A_n`$.

$`\left({\displaystyle \frac{}{T_1}}+U{\displaystyle \frac{}{Y_1}}\right){\displaystyle \frac{^2\mathrm{\Phi }}{x^2}}+\left(\beta U^{\prime \prime }\right){\displaystyle \frac{\mathrm{\Phi }}{Y_1}}`$
$`=`$ $`{\displaystyle \underset{n=1}{\overset{2}{}}}\{{\displaystyle \frac{d^2}{dx^2}}\left[k_n(W_{1n}{\displaystyle \frac{d\phi _n}{dx}}\phi _n{\displaystyle \frac{dW_{1n}}{dx}})\right]4k_n^2\phi _n{\displaystyle \frac{d\phi _n}{dx}}`$
$`+(\phi _n{\displaystyle \frac{d}{dx}}{\displaystyle \frac{d\phi _n}{dx}})({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n)\left\}{\displaystyle \frac{}{Y_1}}\right|A_n|^2`$
This relation will be used later.

Returning to the equation (3), we examine the content of the inhomogeneous part (the right hand side). There are terms which have a space-time dependence

$$\mathrm{exp}\left[ik_nyi\omega _nt\right]$$
(27)
and in order to have solutions of this type, we need to impose solubility conditions, otherwise there will arise secular terms at the order $`\varphi ^{\left(3\right)}`$. To find these restrictions (solubility conditions) we multiply the equation with the function conjugated to (27)

$$\mathrm{exp}\left[ik_ny+i\omega _nt\right]$$

and integrate time on the interval $`(0,2\pi /\omega _n)`$ and space: $`y`$ on the interval $`(0,2\pi /k_n)`$ and $`x`$ on the interval $`(0,L)`$. On the left hand side we obtain zero, so this is what we must obtain also from the right hand side. The conditions which results are:

$$\left(\frac{}{T_2}+c_{g1}\frac{}{Y_2}\right)A_1i\alpha _1\frac{^2A_1}{Y_1^2}i\left(\rho _1\left|A_1\right|^2+\gamma _{12}\left|A_2\right|^2+\lambda _1\right)A_1=0$$

$$\left(\frac{}{T_2}+c_{g2}\frac{}{Y_2}\right)A_2i\alpha _2\frac{^2A_2}{Y_1^2}i\left(\rho _2\left|A_2\right|^2+\gamma _{21}\left|A_1\right|^2+\lambda _2\right)A_2=0$$

The notations which have been introduced are: $`\alpha _n=\frac{I_{n0}}{\mathrm{\Pi }_n}`$, $`\rho _n=\frac{I_{nn}}{\mathrm{\Pi }_n}`$, $`\gamma _{12}=\frac{I_{12}}{\mathrm{\Pi }_1}`$, $`\gamma _{21}=\frac{I_{21}}{\mathrm{\Pi }_2}`$, $`\lambda _n=\frac{I_{n3}}{\mathrm{\Pi }_n}`$ and

$$\mathrm{\Pi }_n=_0^L\frac{\phi _n}{Uc_n}f_n𝑑x,I_{n0}=_0^L\frac{\phi _n}{Uc_n}f_{n0}𝑑x$$

$$I_{nn}=_0^L\frac{\phi _n}{Uc_n}f_{nn}𝑑x,I_{12}=_0^L\frac{\phi _1}{Uc_1}f_{12}𝑑x$$

$$I_{21}=_0^L\frac{\phi _2}{Uc_2}f_{21}𝑑x,I_{n3}=_0^L\frac{\phi _n}{Uc_n}f_{n3}𝑑x$$

The equations satisfied by all these functions are given given in Ref.().

$$f_n=\left(\frac{d^2\phi _n}{dx^2}k_n^2\phi _n\right)$$

$$h_n=U\left(\frac{d^2\phi _n}{dx^2}k_n^2\phi _n\right)\left(\beta U^{\prime \prime }\right)\phi _n+2k_n^2\left(Uc_n\right)$$

$`f_{n0}`$ $`=`$ $`2k_n\left(Uc_{gn}\right)\phi _nk_n\left(Uc_n\right)`$
$`\left(Uc_{gn}\right)\left({\displaystyle \frac{d^2W_{1n}}{dy^2}}k_n^2W_{1n}\right)`$
$`+2k_n^2\left(Uc_n\right)W_{1n}\left(\beta U^{\prime \prime }\right)W_{1n}`$
$`f_{nn}`$ $`=`$ $`\left(2k_nW_{2n}{\displaystyle \frac{d}{dx}}+k_n{\displaystyle \frac{dW_{2n}}{dx}}\right)\left({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n\right)`$
$`+\left(k_n\phi _n{\displaystyle \frac{d}{dx}}+2k_n{\displaystyle \frac{d\phi _n}{dx}}\right)\left({\displaystyle \frac{d^2W_{2n}}{dx^2}}4k_n^2W_{2n}\right)`$
$`f_{12}`$ $`=`$ $`\left(k_{12}W_3{\displaystyle \frac{d}{dx}}+k_2{\displaystyle \frac{dW_3}{dx}}\right)\left({\displaystyle \frac{d^2\phi _2}{dx^2}}k_2^2\phi _2\right)`$
$`+\left(k_2\phi _2{\displaystyle \frac{d}{dx}}+k_{12}{\displaystyle \frac{d\phi _2}{dx}}\right)\left({\displaystyle \frac{d^2W_3}{dx^2}}k_{12}^2W_3\right)`$
$`+\left(k_2\phi _2{\displaystyle \frac{d}{dx}}\alpha _{12}{\displaystyle \frac{d\phi _2}{dx}}\right)\left({\displaystyle \frac{d^2W_4}{dx^2}}\alpha _{12}^2W_4\right)`$
$`\left(\alpha _{12}W_4{\displaystyle \frac{d}{dx}}k_2{\displaystyle \frac{dW_4}{dx}}\right)\left({\displaystyle \frac{d^2\phi _2}{dx^2}}k_2^2\phi _2\right)`$
$`f_{21}`$ $`=`$ $`\left(k_{12}W_3{\displaystyle \frac{d}{dx}}+k_1{\displaystyle \frac{dW_3}{dx}}\right)\left({\displaystyle \frac{d^2\phi _1}{dx^2}}k_1^2\phi _1\right)`$
$`\left(k_1\phi _1{\displaystyle \frac{d}{dx}}+\alpha _{12}{\displaystyle \frac{d\phi _1}{dx}}\right)\left({\displaystyle \frac{d^2W_4}{dx^2}}\alpha _{12}^2W_4\right)`$
$`+\left(k_1\phi _1{\displaystyle \frac{d}{dx}}+k_{12}{\displaystyle \frac{d\phi _2}{dx}}\right)\left({\displaystyle \frac{d^2W_3}{dx^2}}k_{12}^2W_3\right)`$
$`+\left(\alpha _{12}W_4{\displaystyle \frac{d}{dx}}+k_1{\displaystyle \frac{dW_4}{dx}}\right)\left({\displaystyle \frac{d^2\phi _1}{dx^2}}k_1^2\phi _1\right)`$
$$f_{n3}=k_n\left(\phi _n\frac{^3\mathrm{\Phi }}{x^3}\left(\frac{d^2\phi _n}{dx^2}k_n^2\phi _n\right)\frac{\mathrm{\Phi }}{x}\right)$$

We use a boundary value integrator for the inhomogeneous equations and (the NAG routine D02JAF). The results are given in the Figures (7) and (8).

After numerical integrations the follwing values are obtained for the constants: $`\alpha _1=0.39976`$, $`\alpha _2=0.45246`$, $`\rho _1=0.10182`$, $`\rho _2=1.8115`$, $`\gamma _{12}=0.091474`$, $`\gamma _{21}=5.3768`$.

The following transformation (Jeffrey and Kawahara) is performed in Ref.( ) :

$`T`$ $`=`$ $`T_2`$
$`X`$ $`=`$ $`{\displaystyle \frac{1}{\epsilon }}\left(Y_2c_{g1}T_2\right)=Y_1c_{g1}T_1`$
Then the equation of connection between $`\mathrm{\Phi }`$ and $`A_n`$ becomes

$`{\displaystyle \frac{}{Y}}\left[\left(Uc_{g1}\right){\displaystyle \frac{^2\mathrm{\Phi }}{x^2}}+\left(\beta U^{\prime \prime }\right)\mathrm{\Phi }\right]`$
$`=`$ $`{\displaystyle \underset{n=1}{\overset{2}{}}}[k_n{\displaystyle \frac{d^2}{dx^2}}(W_{1n}{\displaystyle \frac{d\phi _n}{dx}}\phi _n{\displaystyle \frac{dW_{1n}}{dx}})4k_n\phi _n{\displaystyle \frac{d\phi _n}{dx}}`$
$`(\phi _n{\displaystyle \frac{d}{dx}}{\displaystyle \frac{d\phi _n}{dx}})({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n)\left]{\displaystyle \frac{}{Y}}\right|A_n|^2`$
The solution of the equation (3.1) is considered to be of the form

$$\mathrm{\Phi }=\underset{n=1}{\overset{2}{}}H_n\left(x\right)\left|A_n\right|^2$$
(29)
and the equation for $`H_n\left(x\right)`$ results:

$`\left(Uc_{g1}\right){\displaystyle \frac{^2H_n}{x^2}}+\left(\beta U^{\prime \prime }\right)H_n`$
$`=`$ $`k_n{\displaystyle \frac{d^2}{dx^2}}\left(W_{1n}{\displaystyle \frac{d\phi _n}{dx}}\phi _n{\displaystyle \frac{dW_{1n}}{dx}}\right)4k_n\phi _n{\displaystyle \frac{d\phi _n}{dx}}`$
$`\left(\phi _n{\displaystyle \frac{d}{dx}}{\displaystyle \frac{d\phi _n}{dx}}\right)\left({\displaystyle \frac{d^2\phi _n}{dx^2}}k_n^2\phi _n\right)`$
The change of variables is performed in the equations for the amplitudes $`A_n`$ and, again, the Eq.(29) is used:

$$i\frac{A_1}{T}+\alpha _1\frac{^2A_1}{Y^2}+\left(\sigma _1\left|A_1\right|^2+\nu _{12}\left|A_2\right|^2\right)A_1=0$$

$$i\left(\frac{}{T}\delta \frac{}{Y}\right)A_2+\alpha _2\frac{^2A_2}{Y^2}+\left(\sigma _2\left|A_2\right|^2+\nu _{21}\left|A_1\right|^2\right)A_2=0$$

where

$$\delta =\frac{1}{\epsilon }\left(c_{g1}c_{g2}\right)$$

$$\sigma _n=\rho _n\frac{1}{\mathrm{\Pi }_n}_0^L\frac{k_n\phi _n}{Uc_n}\left[\phi _n\frac{d^3H_n}{dx^3}\left(\frac{d^2\phi _n}{dx^2}k_n^2\phi _n\right)\frac{dH_n}{dx}\right]𝑑x$$

$$\nu _{12}=\gamma _{12}\frac{1}{\mathrm{\Pi }_1}_0^L\frac{k_1\phi _1}{Uc_1}\left[\phi _1\frac{d^3H_2}{dx^3}\left(\frac{d^2\phi _1}{dx^2}k_1^2\phi _1\right)\frac{dH_2}{dx}\right]𝑑x$$

$$\nu _{21}=\gamma _{21}\frac{1}{\mathrm{\Pi }_2}_0^L\frac{k_2\phi _2}{Uc_2}\left[\phi _2\frac{d^3H_1}{dx^3}\left(\frac{d^2\phi _2}{dx^2}k_2^2\phi _2\right)\frac{dH_1}{dx}\right]𝑑x$$

The following transformation is convenient:

$$A_2=B\mathrm{exp}\left[i\frac{\delta }{2\alpha _2}\left(Y+\frac{\delta }{2}T\right)\right]$$

This leads to the following equations:

$$i\frac{A_1}{T}+\alpha _1\frac{^2A_1}{Y^2}+\left(\sigma _1\left|A_1\right|^2+\nu _{12}\left|B\right|^2\right)A_1=0$$
(30)
$$i\frac{B}{T}+\alpha _2\frac{^2B}{Y^2}+\left(\sigma _2\left|B\right|^2+\nu _{21}\left|A_1\right|^2\right)B=0$$
(31)
The constants are: $`\sigma _1=1.6324`$, $`\sigma _2=1.9331`$, $`\nu _{12}=0.53387`$, $`\nu _{21}=2.9014`$. The constants $`\alpha _{1,2}`$ are called the “dispersion”, $`\sigma _{1,2}`$ are called the nonlinearity (“Landau”) constants ; $`\nu _{12,21}`$ represent the coupling of the two amplitudes.

Numerical experience shows that the dispersion constants are sensitive to the precision of the integration scheme. We use for D02KEF the value for tolerence $`10^8`$.

The system of coupled Cubic Nonlinear Schrodinger Equations has been integrated numerically in relation with many applications , . One important conclusion is that, in the case of a single equation, the sign of the product of the dispersion and of the nonlinearity parameters establishes the possibility of a soliton formation.

### 3.2 Numerical search for classes of admissible envelope solutions

In principle an interactive search for eigenvalues (either with SLEIGN or D02KEF) is able to obtain the eigenvalues $`k_1,k_2`$ eventually after series of less-inspired initializations. However in many situations resonances occur and the sources and/or the potentials in the differential equations for the functions $`W_n`$, $`f_n`$ become discontinuous. Since these eigenvalues cannot be used we have to devise an automatic search over the space of parameters with checks of continuity of the functions and automatic discarding of wrong results. This is accomplished by the routine P01ABF of NAG.

The search reveals that the short and respectively long *radial* wavelength part of the ion wave spectrum have different behavior under the modulation instability. The most affected is the long wavelength part (small $`k_x`$) with $`\omega `$ of the order of the ion-diamagnetic frequency. This corresponds to the dominant part of the ion spectrum, because, as can be seen in numerous numerical simulations, the dominant structures are elongated in the radial direction, with various degrees of tilting compared to the radius. It is also important to note that finding the eigenvalues in the scan of the parameter space weas possible in the condition that the eigenfunctions $`\phi _{1,2}\left(x\right)`$ have low number of oscillations in the radial direction (one or two). Again this corresponds to the important part of the ITG turbulence spectrum.

The variation of the parameters of the system of NSE with the physical parameters of the problem can be studied by the automatic search of eigenvalues, as described before. Below are plotted four graphs showing this dependence.

## 4 Exact solutions of the Nonlinear Schrodinger Equation and the nonlinear stability problem

### 4.1 Introduction

The Nonlinear Schrodinger Equation can be solved *exactly* on an infinite spatial domain using the Inverse Scattering Transform. The first step is to represent the nonlinear equation as a condition of compatibility of two *linear* equations. This is achieved by finding a pair of linear operators (Lax operators) and noticing that for the first of them the eigenfunction problem is a Schrodinger equation leading to a quantum scattering problem. The unknown function of the NSE is the potential of the Schrodinger equation. Using the initial condition of the NSE the scattering data are determined at $`t=0`$. The scattering data have simple time dependence and so they can be evloved in time to the desired moment. Returning from this new set of scattering data (to the potential that has produced it) is achieved by solving an integral equation (Gelfand-Levitan-Marchenko) and provides the solution of the NSE. In general this solution consists of solitons and “radiation”.

The Inverse Scattering Transform on periodic domain , ,  is a more powerful procedure since it reveals and takes advantage of the deep geometric and topological nature of the problem. The admissible solutions of the Lax eigenvalue problem can now only be periodic functions. Since after one period the change of any solution can only be linear (via the *monodromy* matrix) the periodic or antiperiodic functions must be found as eigenfunctions of this monodromy operator (matrix), under the constrain that the eigenvalues are complex of unit absolute value. This singles out a set of complex values of the spectral parameter $`\lambda `$ (the formal eigenvalue in the Lax equation), the *main spectrum*. From this set, the particular values that make the two eigenfunctions to coincide are called *non-degenerate*. Here the squared Wronskian (which is a space and time invariant) has zeros on the complex $`\lambda `$ plane. Then the non-degenerate $`\lambda `$’s are singular points for the Wronskian. Removing the square-root indeterminancy one has to connect pairs of non-degenerate eigenvalues by cuts in the $`\lambda `$ plane. The Wronskian becomes a hyperelliptic Riemann surface. The evolution of the unknown NSE solutions (called hyperelliptic functions $`\mu _j(x,t)`$) on this surface is as complicated as the original NSE equation. However, via the *Abel map* the Riemann surface is mapped onto a torus and on this torus the motion magically becomes *linear*. After obtaining the space-time dependence on the torus we need to come back to the original framework. This is called the Jacobi inversion and can be done exactly in terms of Riemann *theta* functions, giving at the end the exact solution of the NSE.

Not only the geometric approach is more clear but it also allows a treatment of the stability problem for the solutions of the NSE since it allows to trace the changes of the main spectrum after a perturbation of the initial condition. The following subsections provide a more detailed discussion of the IST of a periodic domain. The next section will discuss the stability. This information is available from the abundant literature on the IST and is only mentioned here in order to understand the mechanism governing the stability of the envelope of the ion wave turbulence. For this reason we will focus on the determination of the *main spectrum* and its role in the construction of the solution. The effective steps to be undertaken to obtain the solution will only be briefly described. We strongly recommend the lecture of Ref.() on which this presentation is based.

### 4.2 The Lax operators and the main spectrum

The IST method starts by introducing a pair of linear operators, called the Lax pair (see ) which allow to express the nonlinear equation as a compatibility condition for a system of two linear equations. For the cubic Schrödinger equation

$$i\frac{u}{t}+\frac{^2u}{x^2}+2\left|u\right|^2u=0$$
(32)
the Lax pair of operators is

$$L\left(\begin{array}{cc}i_x& u(x,t)\\ u^{}(x,t)& i_x\end{array}\right)$$
(33)
$$A\left(\begin{array}{cc}i\left|u\right|^22i\lambda ^2& u_x+2i\lambda u\\ u_x^{}+2i\lambda u^{}& i\left|u\right|^2+2i\lambda ^2\end{array}\right)$$
(34)
The action of the operators on a two-component (column) wave function $`\varphi =\left(\begin{array}{c}\varphi _1\\ \varphi _2\end{array}\right)`$ is given by the equations

$$L\varphi =\lambda \varphi $$
(35)
$$\frac{\varphi }{t}=A\varphi $$
(36)
and the condition of compatibility of these two equations $`\varphi _{xt}=\varphi _{tx}`$ is precisely the cubic Nonlinear Schrödinger Equation.

### 4.3 Solving the eigenvalue equation for the operator $`L`$

As in the “infinite-domain” IST, we have to solve the eigenvalue equation Eq.(35) using the initial condition for the function$`u(x,t)`$, $`u(x,0)`$. However, for the particular case of periodic solutions of the NSE, we must have $`u(x+d,0)=u(x,0)`$ where $`d`$ is the length of the spatial period (*i.e.* actually $`L`$ in the previous notation; however we use $`d`$ in this and the next section, keeping $`L`$ for the *Lax* operator). Fixing an arbitrary base point $`x=x_0`$ one considers the two independent solutions of the equation (35) which take the following “initial” values at $`x=x_0`$:

$$\varphi \left(x_0\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)\text{ and }\stackrel{~}{\varphi }\left(x_0\right)=\left(\begin{array}{c}0\\ 1\end{array}\right)$$
(37)
The matrix $`\mathrm{\Phi }(x,x_0;\lambda )`$ of solutions is given by

$$\mathrm{\Phi }(x,x_0;\lambda )\left(\begin{array}{cc}\varphi _1(x,x_0;\lambda )& \stackrel{~}{\varphi }_1(x,x_0;\lambda )\\ \varphi _2(x,x_0;\lambda )& \stackrel{~}{\varphi }_2(x,x_0;\lambda )\end{array}\right)$$
(38)
*i.e.* this solution satisfies the equation

$$L\mathrm{\Phi }=\lambda \mathrm{\Phi }\text{ and }\mathrm{\Phi }(x_0,x_0;\lambda )=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)$$
(39)
It can be seen that the Wronskian of $`(\varphi ,\stackrel{~}{\varphi })`$ is the *determinant* of the matrix of fundamental solutions

$$W(\varphi ,\stackrel{~}{\varphi })=det\left(\mathrm{\Phi }\right)$$

It can be shown that $`\frac{W}{x}=0`$ *i.e.* the Wronskian is constant, which gives

$$det\left(\mathrm{\Phi }\left(x\right)\right)=det\left(\mathrm{\Phi }\left(x_0\right)\right)=1$$

Using the initial condition for $`\mathrm{\Phi }`$ one can solve the equation (39) and obtain $`\mathrm{\Phi }(x,x_0;\lambda )`$; in particular we can obtain it at $`x_0+d`$, *i.e.* after a period :

$$\mathrm{\Phi }(x_0+d,x_0;\lambda )=\left(\begin{array}{cc}\varphi _1(x_0+d,x_0;\lambda )& \stackrel{~}{\varphi }_1(x_0+d,x_0;\lambda )\\ \varphi _2(x_0+d,x_0;\lambda )& \stackrel{~}{\varphi }_2(x_0+d,x_0;\lambda )\end{array}\right)$$
(40)
This is called the monodromy matrix and is noted

$$M(x_0;\lambda )\mathrm{\Phi }(x_0+d,x_0;\lambda )$$
(41)
The monodromy matrix is the matrix of the fundamental solutions calculated for one spatial period. In general the *monodromy matrix* is the element of the monodromy group associated to a loop and to its homotopy class, for a marked point on a manifold (here the circle $`S^1[0,d`$). The linear monodromy group acts by replacing the elements of the vector column by *linear* combinations of the initial ones. For changes of the marked point $`x_0`$ the matrix has two invariants : the trace and the determinant:

$$\left[TrM\right](x_0;\lambda )=\left[TrM\right]\left(\lambda \right)\mathrm{\Delta }\left(\lambda \right)$$

$$detM=1$$

The function $`\mathrm{\Delta }\left(\lambda \right)`$ is called the discriminant and is independent of $`x_0`$.

##### Looking for Bloch (Floquet) solutions

The values the discriminant $`\mathrm{\Delta }\left(\lambda \right)`$ takes on the complex $`\lambda `$-plane control the monodromy and by consequence select the values of $`\lambda `$ for which admissible eigenfunction of the Lax problem exist. In other words $`\mathrm{\Delta }`$ governs the *spectral properties of the operator* $`L`$. To find under what conditions periodic solutions are possible, we construct the Bloch (or Floquet) solutions of the equation $`L\varphi =\lambda \varphi `$. The Bloch function has the property

$$\psi (x+d;\lambda )=e^{ip\left(\lambda \right)}\psi (x;\lambda )$$
(42)
where $`p\left(\lambda \right)`$ is the Floquet exponent. Like any other solution of the Lax eigenproblem $`\psi `$ can be expressed as a linear combination of the two fundamental solutions $`\varphi `$ and $`\stackrel{~}{\varphi }`$

$$\psi (x;\lambda )=A\varphi (x;\lambda )+B\stackrel{~}{\varphi }(x;\lambda )$$

For the particular point $`x_0`$ we have, taking into account the “initial” conditions (37)

$$\psi (x_0;\lambda )=A\varphi (x_0;\lambda )+B\stackrel{~}{\varphi }(x_0;\lambda )=\left(\begin{array}{c}A\\ B\end{array}\right)$$

After one period $`d`$ the function $`\psi `$ is linearly modified by the monodromy matrix and, according to Eq.(42) is multiplied by a complex number of unit absolute value. We write $`m\left(\lambda \right)`$ for this number and look for those $`\lambda `$’s where this is complex of modulus one. The regions of the $`\lambda `$ plane where $`m\left(\lambda \right)`$ has modulus different of one correspond to unstable functions $`\psi `$. We write

$$\psi (x_0+d;\lambda )=m\left(\lambda \right)\psi (x_0;\lambda )$$
(43)
$$M\left(\begin{array}{c}A\\ B\end{array}\right)=m\left(\lambda \right)\left(\begin{array}{c}A\\ B\end{array}\right)$$
(44)
The Bloch functions are invariant directions for the *monodromy operator* $`M`$ acting in the base formed by the fundamental solutions of the Lax eigenproblem. The monodromy operator preserves these directions and multiplies the Bloch vector with $`m\left(\lambda \right)=\mathrm{exp}\left(ip\left(\lambda \right)\right)`$.

$$det\left(Mm\right)=m^2\left(TrM\right)m+detM=m^2\mathrm{\Delta }\left(\lambda \right)m+1=0$$

This gives

$$m^\pm \left(\lambda \right)=\frac{\mathrm{\Delta }\left(\lambda \right)\pm \left(\mathrm{\Delta }^2\left(\lambda \right)4\right)^{1/2}}{2}$$

In general $`\mathrm{\Delta }\left(\lambda \right)`$ is an analytic function on the plane of the complex variable $`\lambda `$. The equation $`\mathrm{I}m\left[\mathrm{\Delta }\left(\lambda \right)\right]=0`$ is a single relation with two unknowns, $`(\mathrm{R}e\lambda ,\mathrm{I}m\lambda )`$ and calculating one of them leaves the dependence on the other as a free parameter. This gives a curve in the plane and the real-$`\lambda `$ axis is such a curve: for $`\mathrm{I}m\lambda =0`$ we have $`\mathrm{I}m\mathrm{\Delta }\left(\lambda \right)=0`$.

We look at the variation of $`\mathrm{\Delta }\left(\lambda \right)`$ on the complex plane to find the effect on $`m^\pm \left(\lambda \right)`$.

* Where the Floquet multiplier is a complex number of magnitude unity, the functions are only effected by a phase factor after one period. The Bloch functions (eigenfunctions of $`L`$) are stable under translations with $`d`$. The *real* $`\lambda `$ axis is a region of stability. The region where $`\mathrm{\Delta }`$ is real and

$$\mathrm{\Delta }\left(\lambda \right)<4$$

is the band of stability, the modulus of $`m`$ being $`1`$.
* When $`\lambda `$ is such that the discriminant ($`\mathrm{\Delta }\left(\lambda \right)Tr\left(M\right)`$ ) is equal to $`4`$, the eigenvalues of the monodromy matrix, $`m\left(\lambda \right)`$, are $`\pm 1`$, and the Bloch functions are *periodic* or *antiperiodic*. This is the Main Spectrum in $`\lambda `$ and is noted $`\left\{\lambda _j\right\}`$. The main spectrum consists of a set of discrete values $`\lambda _j`$.
* for those $`\lambda `$ which gives $`\left|m\left(\lambda \right)\right|1`$ the Bloch eigenfunctions are unstable.

Consider the case where $`\mathrm{\Delta }\left(\lambda \right)`$ is complex. For all values of the parameter $`\lambda `$ on the complex plane, \[with the exception of the Main Spectrum where $`\mathrm{\Delta }\left(\lambda \right)=\pm 2`$ \], the eigenvalues of the monodromy matrix are distinct which means that the eigenfunctions Bloch are distinct and Independent (The Wronskian is different of zero).

$`\psi ^+\left(x+d\right)`$ $`=`$ $`m^+\psi ^+\left(x\right)`$
$`\psi ^{}\left(x+d\right)`$ $`=`$ $`m^{}\psi ^{}\left(x\right)`$
However we can only be intersted in points of the *main spectrum*, as they can provide admissible (periodic) eigenfunctions of the Lax operator. (1) There are points $`\lambda _j`$ in the main spectrum where the two eigenfunctions Bloch are distinct and independent: these points are called *degenerate* (with the meaning that for two independent eigenfunctions there is only one eigenvalue). (2) There are points $`\lambda _j`$ where the two eigenfunctions Bloch are NOT independent: these points are called *nondegenerate*. For such values of $`\lambda `$ the Wronskian is zero.

The important class of initial conditions $`u\left(x,t=0\right)`$ with the property that there is only a *finite number of non-degenerate eigenvalues* $`\lambda _j`$ is called finite-band potential.

### 4.4 The “squared” eigenfunctions

Starting from the initial condition $`u(x,0)`$ and leting it evolve in time according to the NSE, $`u(x,0)\stackrel{NSE}{}u(x,t)`$, the discriminant remains invariant. All the spectral structure obtained from the discriminant $`\mathrm{\Delta }\left(\lambda \right)`$, *i.e.* the main spectrum, the stability bands, are invariant. The eigenvalues of the monodromy matrix $`\lambda _j`$ are invariant.

#### 4.4.1 The two Bloch functions for the operator $`L`$

Consider the two Bloch eigenfunctions of the operator $`L`$:

$$\psi ^+=\left(\begin{array}{c}\psi _1^+\\ \psi _2^+\end{array}\right)\text{ and }\psi ^{}=\left(\begin{array}{c}\psi _1^{}\\ \psi _2^{}\end{array}\right)$$

and define the following functions

$$f(x,t;\lambda )\frac{i}{2}\left(\psi _1^+\psi _2^{}+\psi _2^+\psi _1^{}\right)$$
(45)
$$g(x,t;\lambda )\psi _1^+\psi _1^{}$$
(46)
$$h(x,t;\lambda )\psi _2^+\psi _2^{}$$
(47)
These functions have the property that they are periodic, $`f(x+d,t;\lambda )=f(x,t;\lambda )`$for all values of $`\lambda `$. Since $`f`$ is composed of $`\psi ^+`$ and $`\psi ^{}`$ and these are Bloch solutions for the eigenvalue equation $`L\phi =\lambda \phi `$ we have

$$\frac{f}{x}=u^{}guh$$
(48)
$$\frac{g}{x}=2i\lambda g2uf$$
(49)
$$\frac{h}{x}=2i\lambda h+2u^{}f$$
(50)
and similarly we write the time derivatives taking account of the compatibility condition which involves the operator $`A`$.

$$\frac{f}{t}=\left(u_x^{}+2i\lambda u^{}\right)g+i\left(u_x+2i\lambda u\right)h$$
(51)
$$\frac{g}{t}=2i\left(\left|u\right|^22\lambda ^2\right)g+2i\left(u_x+2i\lambda u\right)f$$
(52)
$$\frac{h}{t}=2i\left(\left|u\right|^22\lambda ^2\right)h2i\left(u_x^{}+2i\lambda u^{}\right)f$$
(53)
It is shown in Ref.() that the condition on the initial function $`u(x,0)`$ for there to be only a finite number of nondegenerate points in the main spectrum is equivalent to the requirement that $`f,g`$ and $`h`$ be finite-order polynomials in the parameter $`\lambda `$, which we take of degree $`N+1`$: $`f(x,t;\lambda )=_{j=0}^{N+1}f_j(x,t)\lambda ^j`$, $`g(x,t;\lambda )=_{j=0}^{N+1}g_j(x,t)\lambda ^j`$, $`h(x,t;\lambda )=_{j=0}^{N+1}h_j(x,t)\lambda ^j`$. From the Eqs.(48 \- 50) it results however that $`g`$ and $`h`$ have degree $`N`$ in $`\lambda `$.

One can check that the following combination is invariant in time and space and we note that it is actually the square of the Wronskian:

$$\frac{}{x}\left(f^2gh\right)=\frac{}{t}\left(f^2gh\right)=0$$
(54)
$$f^2gh=\frac{1}{4}\left[W(\psi ^+,\psi ^{})\right]^2$$

The Wronskian only depends on the spactral parameter $`\lambda `$. We also know that for a subset of the *main spectrum*, the nondegenerate eigenvalues, the Wronskian is zero. Then the number of the nondegenerate eigenvalues is $`2N+2`$ (the degree of $`f^2gh`$ as a polynomial in $`\lambda `$) and we express the function $`f^2gh`$ as a polinomial with constant coefficients (not depending on $`x`$ and $`t`$).

$$\frac{1}{4}\left[W(\psi ^+,\psi ^{})\right]^2=f^2ghP\left(\lambda \right)=\underset{k=1}{\overset{2N+2}{}}P_k\lambda ^k=\underset{j=0}{\overset{2N+2}{}}\left(\lambda \lambda _j\right)$$
(55)
#### 4.4.2 Product expansion of $`g`$ and its zeros $`\mu _j(x,t)`$. Introduction of the *hyperelliptic* functions $`\mu _j`$

The functions $`g`$ and $`h`$ are both of order $`N`$ in $`\lambda `$. For $`g`$ we will note the $`N`$ roots by $`\mu _j(x,t)`$. The coefficient of $`\lambda ^{2N+2}`$ in $`f^2gh`$ is $`f_{N+1}^2`$. Due to Eq.(54) it is a constant which can be taken $`1`$. Now since $`f_{N+1}=1`$ the coefficient of $`\lambda ^N`$ in $`g`$ is $`iu(x,t)`$. Written as a product

$$g=iu(x,t)\underset{j=1}{\overset{N}{}}\left[\lambda \mu _j(x,t)\right]$$
(56)
By similar arguments we find that the coefficient of $`\lambda ^N`$ in $`h`$ is $`iu^{}(x,t)`$ and the zeros of $`h`$ are $`\mu _j^{}(x,t)`$. The functions $`\mu _j(x,t)`$ are called the *hyperelliptic functions*. It will be proved below that finding the *hyperelliptic* functions leads immediately to the solution $`u(x,t)`$.

Now we calculate the expression in Eq.(55) for $`\lambda `$ = a zero of the function $`g`$, *i.e.* $`\lambda `$ is equal to *hyperelliptic function* $`\mu _j(x,t)`$:

$$f^2\left(x,t;\lambda =\mu _m\right)gh=f^2(x,t;\mu _m)=P\left(\lambda =\mu _m\right)$$

or

$$f(x,t;\mu _m)=\sigma _m\sqrt{P\left(\mu _m\right)}$$
(57)
Here the factor $`\sigma _m`$ is a sheet index that indicates which sheet of the Riemann surface associated with $`\sqrt{P\left(\lambda \right)}`$ the complex $`\mu _m`$ lies on.

Let us calculate from Eq.(49) and Eq.(56) the derivative at $`x`$ of $`\mu _m`$:

$`{\displaystyle \frac{g(x,t;\lambda )}{x}}`$ $`=`$ $`i{\displaystyle \frac{u(x,t;\lambda )}{x}}{\displaystyle \underset{j=1}{\overset{N}{}}}\left(\lambda \mu _j(x,t)\right)`$
$`iu(x,t;\lambda ){\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{\mu _j(x,t)}{x}}{\displaystyle \underset{k=1,kj}{\overset{N}{}}}\left(\lambda \mu _k(x,t)\right)`$
Replace here a zero of $`g`$: $`\lambda =\mu _m`$

$$\frac{g(x,t;\mu _m)}{x}=iu(x,t;\mu _m)\frac{\mu _m(x,t)}{x}\underset{k=1,km}{\overset{N}{}}\left(\mu _m\mu _k(x,t)\right)$$

On the other hand we have, from Eq.(49)

$`{\displaystyle \frac{g}{x}}`$ $`=`$ $`2i\mu _mg\left(x,t;\lambda =\mu _m\right)2u\left(x,t;\lambda =\mu _m\right)f\left(x,t;\lambda =\mu _m\right)`$
$`=`$ $`2uf`$
then, since we have calculated $`f\left(x,t;\lambda =\mu _m\right)`$, Eq(57), we obtain

$$\frac{\mu _m(x,t)}{x}=\frac{2i\sigma _m\sqrt{P\left(\mu _m\right)}}{_{k=1,km}^N\left(\mu _m\mu _k(x,t)\right)}$$
(58)
for $`m=1,2,\mathrm{},N`$.

The coefficients of the term with $`\lambda ^N`$ in the equation (49) are matched and we obtain $`iu_x=2ig_{N1}2uf_N`$ .Using Eq.(56) we find

$$_x\mathrm{ln}u=2i\left(\underset{j=1}{\overset{N}{}}\mu _j+f_N\right)$$

From the Eq.(55) we obtain the coefficient of the term $`\lambda ^N`$, *i.e.* $`f_N`$

$$f_N=\frac{1}{2}\underset{k=1}{\overset{2N+2}{}}\lambda _k$$

Then from the preceding equation it results

$$_x\mathrm{ln}u=2i\left(\underset{j=1}{\overset{N}{}}\mu _j\frac{1}{2}\underset{k=1}{\overset{2N+2}{}}\lambda _k\right)$$
(59)
The same procedure is applied to the equation for $`g_t`$, (52) and we obtain

$$\frac{\mu _j(x,t)}{t}=2\left(\underset{mj}{}\mu _m\frac{1}{2}\underset{k=1}{\overset{2N+2}{}}\lambda _k\right)\frac{\mu _j(x,t)}{x}$$
(60)
$`_t\mathrm{ln}u`$ $`=`$ $`2i\left[{\displaystyle \underset{j>k}{}}\lambda _j\lambda _k{\displaystyle \frac{3}{4}}\left({\displaystyle \underset{k=1}{\overset{2N+2}{}}}\lambda _k\right)^2\right]`$
$`4i\left[\left({\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{2N+2}{}}}\lambda _k\right)\left({\displaystyle \underset{j=1}{\overset{N}{}}}\mu _j\right)+{\displaystyle \underset{j>k}{}}\mu _j\mu _k\right]`$
We can see that the problem has been reformulated: from the equation for the function $`u(x,t)`$ the Lax-operators formulation leads to the problem for Bloch functions, $`\psi ^+`$ and $`\psi ^{}`$ and their Wronskian; then, using the “squared” eigenfunctions $`f`$, $`g`$ and $`h`$ we arrive at a formulation for the *hyperelliptic* functions $`\mu _j(x,t)`$ . Finding $`\mu _j(x,t)`$ leads immediately to $`u(x,t)`$.

These operations have a significative geometric counterpart: every point of the complex plane of the spectral parameter $`\lambda `$ is mapped, via Eq.(55) on the two-sheeted Riemann surface $`\sqrt{\frac{1}{4}W(x,t;\lambda )}=\sqrt{P\left(\lambda \right)}`$ with singularities at the *nondegenerate* points of the main spectrum, $`\lambda _j`$. Removing the indeterminancy (by cutting and glueing the two sheets) we obtain a hyperelliptic Riemann surface of genus $`g=N`$. It results that we have to consider the variables $`\mu _j(x,t)`$, $`j=1,N`$ , as points on this surface but the motion of $`\mu _j(x,t)`$ is not simpler than the equation itself.

Now the next steps would be: (1) using the initial condition $`\mu _j(0,0)`$ on the function $`\mu _j(x,t)`$ we can solve the two equations (58) and (60) and find $`\mu _j(x,t)`$; (2) then using the initial condition $`\left|u(0,0)\right|`$ for the function $`u(x,t)`$ we can solve the equations (59) and (4.4.2). The procedure will be:

* Take the parameters $`\left\{\lambda _j\right|j=1,2,\mathrm{},2N+2\}`$ as known; these parameters are non-degenerate eigenvalues from the main spectrum, they can be found from the equation $`\mathrm{\Delta }\left(\lambda \right)=\pm 2`$ and $`\mathrm{\Delta }`$ is determined from $`u(x,0)`$, the initial condition.
* Choose initial conditions $`\mu _j(0,0)`$ and $`\left|u(0,0)\right|`$ (see below)
* solve the equations for $`\mu _j(x,t)`$ ; this means to find the *hyperelliptic functions*. Then solve the equations for $`u(x,t)`$.

### 4.5 Solution: the two-sheeted Riemann surface (hyperelliptic genus-$`N`$ Riemann surface)

The variable $`\mu _j`$’s are points lying on the two-sheeted Riemann surface associated with

$$\sqrt{P\left(\lambda \right)}=\left(\underset{k=1}{\overset{2N+2}{}}\left(\lambda \lambda _k\right)\right)^{1/2}R\left(\lambda \right)$$

which has branch cuts at each of the nondegenerate points $`\lambda _k`$. To eliminate the indeterminacy related to the square-root singularities (located at $`\lambda _j`$’s) the surfaces must be cut and reglued, obtaining a complex manifold of complex dimension one, a hyperelliptic Riemann surface. We need some constructions on this surface.

#### 4.5.1 Holomorphic differential forms and cycles on the Riemann surface

On this hyperelliptic Riemann surface (denoted $`M`$) it is possible to define $`N`$ linearly independent holomorphic (regular) differentials. The following is the canonical basis of differentials one-forms defined on the manifold $`M`$:

$$dU_1\frac{d\lambda }{R\left(\lambda \right)}$$

$$dU_2\frac{\lambda d\lambda }{R\left(\lambda \right)},\mathrm{}$$

$$dU_j\frac{\lambda ^jd\lambda }{R\left(\lambda \right)},\mathrm{}$$

$$dU_N\frac{\lambda ^{N1}d\lambda }{R\left(\lambda \right)}$$

On the surface $`M`$ there are $`2N`$ topologically distinct closed curves (cycles). On a torus ($`N=1`$, called *elliptic surface*) there are two curves which cannot be deformed one into another. A more general Riemann surface, with $`N>1`$ can be shown to be equivalent to a sphere with $`N`$ handles and is called *hyperelliptic Riemann surface*. The number $`N`$ is called the genus of the surface. We have

$$N=\text{genus of the surface}=\text{number of independent holomorphic differential forms}$$

There are $`2N`$ cycles which are split into two classes: $`a_j`$ cycles and $`b_j`$cycles. Each of these cycles has a specified direction (an arrow) attached to it. To construct the cycles, the rules to be applied are:

1. $`a_j`$ cycles do not cross any other $`a_j`$ cycle; $`b_j`$ cycles do not cross any other $`b_j`$ cycle;
2. the cycle $`a_k`$ intersects $`b_k`$ only once and intersects no other $`b`$ cycle;
3. the intersection is such that, at the point of intersection, the vector tangent to the cycle $`a_k`$, the vector tangent to the cycle $`b_k`$ and the normal to the tangent plane of $`M`$ represent a system of three vectors compatible with the *orientation* of $`M`$.

This can be represented as:

$$aa=0,bb=0,a_jb_k=\delta _{jk}$$

#### 4.5.2 Periods and matrices of periods

With the set of holomorphic differential forms and the set of cycles we can define the matrices of periods, $`A`$ and $`B`$

$$A_{kj}_{a_j}𝑑U_k$$

$$B_{kj}=_{b_j}𝑑U_k$$

and the matrices $`A`$ and $`B`$ are invertible. A change of basis of holomorphic differential forms is a linear transformation represented by the matrix $`C`$:

$$d\psi _j\underset{k=1}{\overset{N}{}}C_{jk}dU_k$$

We can choose the matrix $`C`$ as

$$C=A^1$$

Then the periods in the new basis becomes

$$_{a_n}𝑑\psi _j=\underset{k=1}{\overset{N}{}}C_{jk}_{a_n}𝑑U_k=\underset{k=1}{\overset{N}{}}C_{jk}A_{kn}=\delta _{jn}$$

$$_{b_n}𝑑\psi _j=\tau _{jn}$$

where

$$\tau =A^1B$$

It can be shown that the matrix $`\tau `$ is symmetric $`\tau _{jk}=\tau _{kj}`$ and has a positive-definite imaginary part $`\mathrm{I}m\tau >0`$.

#### 4.5.3 The Abel map

The Abel map is defined from the Riemann surface $`M`$ to the space $`C^N`$ and associates to a hyperelliptic Riemann surface (a sphere with $`N`$ hadles) a $`N`$-torus. Using the differentials $`d\psi _j`$ one constructs a change of variables by the following procedure:

* choose a base point on $`M`$ and call it $`p_0`$;
* define the variables $`W_j(x,t)`$ as integrals on the Riemann $`\lambda `$-surface of the *differential forms* from the arbitrary point $`p_0`$ to the point which is the hyperelliptic function $`\mu _k(x,t)`$:

$`W_j(x,t)`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle _{p_0}^{\mu _k(x,t)}}𝑑\psi _j`$
$`=`$ $`{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \underset{m=1}{\overset{N}{}}}C_{jm}{\displaystyle _{p_0}^{\mu _k(x,t)}}{\displaystyle \frac{\lambda ^{m1}d\lambda }{R\left(\lambda \right)}}`$
These variables have the remarkable property that their $`x`$ and $`t`$ dependence is trivial

$$\frac{d}{dx}W_j(x,t)=\underset{k=1}{\overset{N}{}}\underset{m=1}{\overset{N}{}}C_{jm}\frac{\mu _k^{m1}\frac{d\mu _k}{dx}}{\sigma _kR\left(\mu _k\right)}$$

and using

$$\frac{d\mu _k(x,t)}{dx}=2i\sigma _k\frac{R\left(\mu _k\right)}{_{nk}\left(\mu _k\mu _n\right)}$$

it results

$$\frac{d}{dx}W_j(x,t)=\underset{m=1}{\overset{N}{}}C_{jm}\underset{k=1}{\overset{N}{}}\frac{2i\mu _k^{m1}}{_{nk}\left(\mu _k\mu _n\right)}$$

To calculate the sum

$$\underset{k=1}{\overset{N}{}}\frac{\mu _k^{m1}}{_{nk}\left(\mu _k\mu _n\right)}$$

one can use the contour integral

$$I_m\frac{1}{2\pi i}_C\frac{\lambda ^{m1}d\lambda }{_{nk}\left(\mu _k\mu _n\right)}$$

where the contour $`C`$ encloses all of the poles $`\mu _n`$ counterclockwise. By the residue theorem

$$\underset{k=1}{\overset{N}{}}\frac{\mu _k^{m1}}{_{nk}\left(\mu _k\mu _n\right)}=I_m$$

One can evaluate this integral by compactifying the $`\lambda `$ plane into a sphere and noticing that the contour encloses the pole at $`z=\mathrm{}`$. Evaluating the residuu at $`z=\mathrm{}`$

$$\underset{k=1}{\overset{N}{}}\frac{\mu _k^{m1}}{_{nk}\left(\mu _k\mu _n\right)}=\delta _{m,N}$$

and from this we obtain

$$\frac{d}{dx}W_j(x,t)=2iC_{j,N}\frac{1}{2\pi }k_j$$

A similar calculation gives

$`{\displaystyle \frac{d}{dt}}W_j(x,t)`$ $`=`$ $`4i\left[C_{j,N1}+\left({\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{2N+2}{}}}\lambda _k\right)C_{j,N}\right]`$
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\mathrm{\Omega }_j`$
It results from these calculations

$$W_j(x,t)=\frac{1}{2\pi }\left(k_jx+\mathrm{\Omega }_jt+d_j\right)$$

where $`d_j`$ is a phase which is determined by the initial condition on $`\mu _k`$.

The Abel map has *linearized* the motion of the points $`\mu _k(x,t)`$. Now, in order to determine the function $`u(x,t)`$ we must invert the Abel map, returning from the variables $`W`$ to $`\mu `$. This is the *Jacobi inversion problem*.

#### 4.5.4 The Jacobi inversion problem

The return from the variables $`W_j(x,t)`$ of the $`N`$-torus back to the variables $`\mu _k(x,t)`$ of the Riemann surface $`M`$ (notice we will only need particular combinations of the functions $`\mu `$) can be done in an exact way using the Riemann $`\theta `$ function.

The argument of $`\theta `$ is the N-tuple of complex numbers $`\overline{z}C^N`$.

$$\theta \left(z|\tau \right)\underset{m_1=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{}\underset{m_N=\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{exp}\left(2\pi i\overline{m}\overline{z}+\pi i\overline{m}\tau \overline{m}\right)$$

with the notations $`\overline{m}\overline{z}=_{k=1}^Nm_kz_k`$ and $`\overline{m}\tau \overline{m}=_{k,j=1}^N\tau _{jk}m_jm_k`$. Adding to the vector $`z`$ a vector with a single nonzero element in the position $`k`$ we obtain from the definition:

$$\theta \left(z+e_k|\tau \right)=\theta \left(z|\tau \right)$$

which means that the $`\theta `$ -function has $`N`$ real periods. Adding the full $`k`$-th column of the matrix $`\tau `$, we obtain:

$$\theta \left(z+\tau _k|\tau \right)=\mathrm{exp}\left(2\pi iz_k\pi i\tau _{kk}\right)\theta \left(z|\tau \right)$$
(62)
where $`\tau _{kk}`$ is the diagonal element of the $`\underset{¯}{\tau }`$ matrix.

There are two quantities which must be determined in order to have the explicit form of the solution. The space and respectively time equations for $`u(x,t)`$ depend on the quantity $`_{j>k}\mu _j\mu _k=\frac{1}{2}\left[\left(_{m=1}^N\mu _m\right)^2_{m=1}^N\mu _m^2\right]`$. So we must determine the following combinations of the variables $`\mu _j`$:

$$\underset{m=1}{\overset{N}{}}\mu _m(x,t)\text{and}\underset{m=1}{\overset{N}{}}\mu _m^2(x,t)$$

In order to calculate these two quantities, we shall start by introducing a series of functions , $`\psi _j\left(p\right)`$ defined on the Riemann surface $`M`$:

$$\psi _j\left(p\right)=_{p_0}^p𝑑\psi _j$$

where $`p`$is a point on$`M`$ , $`d\psi _j`$ is the $`j`$ -th holomorphic differential form defined on $`M`$ and $`p_0`$ is an arbitrary fixed point on $`M`$. The functions $`\psi _j\left(p\right)`$ are multiple valued since the contour is defined up to addition of any combination of the cycles on $`M`$.

Now consider the function, $`F\left(p\right)`$:

$$F\left(p\right)\theta \left(\psi \left(p\right)K|\tau \right)$$

where $`\theta `$ is the $`N`$-dimensional Riemann *theta* function associated with the surface $`M`$; $`\psi \left(p\right)`$ is the $`N`$-dimensional column of complex functions $`\psi _j\left(p\right)`$; $`K`$ is an $`N`$-dimensional column of complex numbers independent of the point $`p`$. The function $`F\left(p\right)`$ is multivalued on the Riemann surface $`M`$ for the same reason as $`\psi `$: moving the point $`p`$ on the surface such as to turn around one of the $`b`$ cycle and returning to the initial position, the functions $`\psi _j`$ will add elements of the $`\tau `$ matrix. This will make the function $`F`$ to change as imposed by the properties of the $`\theta `$ function, shown in Eq.(62).

In order to render the function $`F`$ single-valued, we replace its domain $`M`$ by a new surface, obtained from $`M`$ by dissecting it in a canonical fashion, along a basis of cycles. This new surface, $`M^{}`$ is simply connected and has the bord composed of a number of arcs equal to $`4N`$.

This operation is necessary in order to render the function $`F\left(p\right)`$ entire and allow us define the inegral

$$I_0=\frac{1}{2\pi i}_M^{}d\mathrm{ln}F\left(p\right)$$
(63)
around the contour of the surface $`M^{}`$. Applying Cauchy theorem, the integral is the *number of zeroes* of $`F\left(p\right)`$ on the surface $`M^{}`$. This number is $`N`$ (Riemann).

We now impose the condition that the $`N`$ zeros of $`F\left(p\right)`$ coincide with the $`N`$ points $`(\mu _j(x,t),\sigma _j)`$. This fixes the values of the $`N`$ complex numbers $`K_j\left(j=1,\mathrm{},N\right)`$. By doing so the contour integral (63) calculated with the reziduum theorem will involve the variables $`\mu _j`$.

##### Calculation of the two combination of $`\mu `$’s

We must remember that the complex $`\lambda `$ plane is covered by the two-sheeted Riemann surface whose compact version is $`M`$. Further this is mapped by Abel map onto the Jacobi $`N`$-torus. To a variable on the Riemann surface $`M`$ (formally also on $`M^{}`$), say $`p`$, corresponds a certain $`\lambda \left(p\right)`$. One defines the following integrals, which are proved to be real constants :

$$I_1=\frac{1}{2\pi i}_M^{}\lambda \left(p\right)d\mathrm{ln}\left[F\left(p\right)\right]A_1$$

$$I_2=\frac{1}{2\pi i}_M^{}\lambda ^2\left(p\right)d\mathrm{ln}\left[F\left(p\right)\right]A_2$$

They can be evaluated by the residuu theorem. By the choice of the constants $`K`$’s, the zeroes of the of the function $`F\left(p\right)`$ are located at $`\mu _j`$. Then the residues are just the integrand ($`\lambda `$) calculated in $`\mu `$ plus the residuu at the infinite, $`\pm \mathrm{}`$:

$$I_1=\underset{m=1}{\overset{N}{}}\mu _m+\underset{\lambda \mathrm{}^+}{\mathrm{R}es}\left[\lambda \left(p\right)d\mathrm{ln}F\left(p\right)\right]+\underset{\lambda \mathrm{}^{}}{\mathrm{R}es}\left[\lambda \left(p\right)d\mathrm{ln}F\left(p\right)\right]$$

$$I_2=\underset{m=1}{\overset{N}{}}\mu _m^2+\underset{\lambda \mathrm{}^+}{\mathrm{R}es}\left[\lambda ^2\left(p\right)d\mathrm{ln}F\left(p\right)\right]+\underset{\lambda \mathrm{}^{}}{\mathrm{R}es}\left[\lambda ^2\left(p\right)d\mathrm{ln}F\left(p\right)\right]$$

The reason to write two residues at infinity is that $`\lambda =\mathrm{}`$ is not a branching point which means that there are two points on the manifold $`M`$ corresponding to $`\lambda =\mathrm{}`$, one on each sheet of the surface. We note $`r^\pm `$the value of $`\psi \left(p\right)`$ when $`p`$ is the point on $`M^{}`$ corresponding to $`\lambda \mathrm{}^\pm `$.

The result is

$$A_1=\underset{j=1}{\overset{N}{}}\mu _j+\frac{i}{2}\frac{}{x}\mathrm{ln}\left[\frac{F\left(r^{}K\right)}{F\left(r^+K\right)}\right]$$

or

$$A_1=\underset{j=1}{\overset{N}{}}\mu _j+\frac{i}{2}\frac{}{x}\mathrm{ln}\left[\frac{F\left(r^{}K\right)}{F\left(r^+K\right)}\right]$$

with $`\theta ^\pm =F\left(r^\pm K\right)`$. The expression of $`I_2`$

$$A_2=\underset{j=1}{\overset{N}{}}\mu _j^2\frac{1}{4}\frac{}{x}\mathrm{ln}\left[F\left(r^+K\right)F\left(r^{}K\right)\right]+\frac{i}{4}\frac{}{t}\mathrm{ln}\left[\frac{F\left(r^{}K\right)}{F\left(r^+K\right)}\right]$$

##### Return to the equations for the function $`u`$

With these expressions we come back to the two equations for the two partial derivatives of $`u`$.

$$\frac{}{x}\mathrm{ln}u=\frac{}{x}\mathrm{ln}\frac{F\left(r^+K\right)}{F\left(r^{}K\right)}+2iA_1i\underset{m=1}{\overset{2N+2}{}}\lambda _m$$

$$\frac{}{t}\mathrm{ln}u=\frac{}{t}\mathrm{ln}\frac{F\left(r^+K\right)}{F\left(r^{}K\right)}+i\text{const}$$

##### The solution

Let us note

$$\omega _0=Q$$

$$k_0=2A_1\underset{m=1}{\overset{2N+2}{}}\lambda _m$$

which are “external” frequency and wavelength.

The solution is

$$u(x,t)=u(0,0)\mathrm{exp}\left(ik_0i\omega _0t\right)\frac{\theta (W^{}|\tau )}{\theta \left(W^+|\tau \right)}$$

where

$$W_j^\pm =\frac{1}{2\pi }\left(k_jx+\mathrm{\Omega }_jt+\delta _j^\pm \right)$$

The phases $`\delta _j^\pm `$ are the part of $`r^\pm K_j`$ which is independent of $`(x,t)`$.

## 5 Stability of the envelope solutions

In this Section we return to the equations obtained by multiple space and time scales analysis. Then the meaning of the variables $`\varphi `$, $`y`$, $`t`$ is: first order (envelope) amplitude of the potential fluctuations in the ion turbulence and respectively large space and slow time variables. We focus on a single NSE in order to examine the properties of the stability of the solutions. The equation has the generic form

$$i\frac{\varphi }{t}+\alpha \frac{^2\varphi }{y^2}+2\sigma \left|\varphi \right|^2\varphi =0$$
(64)
and represents the equation for the envelope of a fast oscillation arising as solution of the barotropic equation. On infinite spatial domain this equation has soliton as solutions. On a periodic spatial domain (which is our case) the solutions can be cnoidal waves and they may be modulationally unstable. The instability can be examined in the neighbourhood of a solution by imposing a slight deviation and studying its time behaviour by classical perturbation expansion. This analysis and comparison with the geometric-algebraic results have been carried out in Ref.().

Consider $`\varphi (y,t)`$ is a solution and perturb it at the initial moment $`t=0`$ as

$$\stackrel{~}{\varphi }(y,t)=\varphi (y,0)+\epsilon h\left(y\right)$$

The question is to find the relative behaviour of the two functions $`\varphi (y,t)`$ and $`\stackrel{~}{\varphi }(y,t)`$. It they depart exponentially for time close to the initial moment, then $`\varphi (y,t)`$ is an unstable solution. Suppose the initial solution $`\varphi (y,t)`$ is the envelope of an exact plane wave *i.e.* it is constant in space. As a function of time it must have the form $`\varphi (y,t)=p\mathrm{exp}\left(2iqt\right)`$ and we find

$$\varphi (y,t)=p\mathrm{exp}\left(2i\sigma p^2t\right)$$
(65)
To this solution we add a perturbation

$$\stackrel{~}{\varphi }(y,t)=p\mathrm{exp}\left(2i\sigma p^2t\right)\left\{1+\epsilon \left[A_1\mathrm{exp}\left(ikyi\mathrm{\Omega }t\right)+A_2\mathrm{exp}\left(iky+i\mathrm{\Omega }t\right)\right]\right\}$$

where $`A_1`$ and $`A_2`$ are *real* coefficients. Linearising the equation for small $`\epsilon `$ we obtain:

$`\left\{\left(2p^3\sigma +p\mathrm{\Omega }p\alpha k^2\right)\mathrm{exp}\left(ikyi\mathrm{\Omega }t\right)+2p^3\sigma \mathrm{exp}\left(iky+i\mathrm{\Omega }t\right)\right\}A_1+`$
$`\left\{\left(2p^3\sigma p\mathrm{\Omega }^{}p\alpha k^2\right)\mathrm{exp}\left(iky+i\mathrm{\Omega }t\right)+2p^3\sigma \mathrm{exp}\left(ikyi\mathrm{\Omega }t\right)\right\}A_2`$
$`=`$ $`0`$
Taking the complex conjugate we obtain another equation for the two *real* coefficients $`A_1`$ and $`A_2`$. The condition of compatibility of this system of equations (the determinant equals zero) can only be fulfilled if the factor multiplying the function of $`y`$ and $`t`$ is allways zero. This gives us the dispersion relation

$$\mathrm{\Omega }=\pm \sqrt{\alpha }k\left(\alpha k^24p^2\sigma \right)^{1/2}$$
(66)
which shows that for

$$\left|k\right|<2|p|\sqrt{\frac{\sigma }{\alpha }}$$

(*i.e.* for long wavelengths) the two solutions $`\varphi (y,t)`$ and $`\stackrel{~}{\varphi }(y,t)`$ diverge exponentially in time and the function $`\varphi (y,t)`$ is modulationally unstable. From Eq.(66) we notice that when $`\sigma `$ (the Landau coefficient)$`<0`$ and $`\alpha `$ (the dispersion coefficient)$`>0`$ the frequency $`\mathrm{\Omega }`$ is real and there is no linear instability. This remark is useful in connection with the numerical calculations described above, which connect the physical parameters with the coefficients $`(\alpha ,\sigma )`$ of the NSE. Sampling the space of the physical parameters we must retain those points where the linear instability is possible. In the above formula $`\sigma `$ , $`\alpha `$ and $`p`$ are dimensionless constants and the variables of the equation are normalised : time to $`t_0`$ and length to $`L`$. The value of $`p`$ is connected with the frequency of oscillation of the uniform poloidal envelope $`\omega _{unif}^{phys}`$ ($`s^1`$) by

$$p=\sqrt{\frac{t_0}{2\sigma }\omega _{unif}^{phys}}$$
(67)
and the condition of linear instability is for the poloidal wavenumber $`\overline{m}`$

$$\overline{m}<\sqrt{2}\frac{a}{L}\sqrt{\frac{\sigma }{\alpha }}\left(t_0\omega _{unif}^{phys}\right)^{1/2}$$
(68)
In the case of the following parameters: $`\alpha =0.399`$, $`\sigma =1.452`$, $`d=2\pi a/L=2\pi \left(1/0.3\right)=\mathrm{20.\hspace{0.17em}944}`$, $`\omega _{unif}^{phys}=10^5`$ (sec<sup>-1</sup>) we have

$$\overline{m}8$$

For the Nonlinear Schrodinger Equation it is possible to follow the perturbed solution much further in time in an exact manner. In Ref.() it is shown that from a solution it is possible to return toward the initial time and to recover the results of the linear analysis. This is possible because the equation can be exactly solved for initial conditions starting arbitrary close to the envelope of the plane wave. The case of the plane wave (uniform envelope) is particular due to the simple form of the main spectrum, as will be discussed below.

### 5.1 Modulation instability of the envelope of a plane wave

#### 5.1.1 General method

The geometric-algebraic setting of Inverse Scattering Transform is the appropriate framework for understanding the nature of the stability of the exact solutions of the NSE. This is because the solution is strongly dependent on the topology of the hyperelliptic Riemann curve which in turn is dependent on the initial function. A perturbation of the initial condition may change the main spectrum turning degenerate eigenvalues into new nondegenerate pairs of eigenvalues. This means that the hyperelliptic curve will have a different genus $`N^{}>N`$ and new cycles and differential forms should be defined. The solution depends precisely on this topology.

We order to investigate the stability of an exact $`N`$-band solution of the NSE we have to develop in detail the construction of solution. In the particular case where the solution is the plane wave this analysis is simpler since all the spectral properties can be found analytically. It has been proved  that under the perturbation of the initial condition the degeneracies in the main spectrum will be removed and new degrees of freedom will become active. In the new solution some of the new degrees of freedom have little effect and can be neglected.

#### 5.1.2 Calculation of the main spectrum associated to the uniform amplitude solution

The unperturbed solution is the envelope of the plane wave, Eq.(65). To deteremine the main spectrum, the solution is considered initial condition (at $`t=0`$) and introduced in the Lax operator. The eigenvalue problem of the Lax operator has the form

$$\left(\begin{array}{cc}i\sqrt{\frac{\sigma }{\alpha }}\frac{}{y}& 1\\ 1& i\sqrt{\frac{\sigma }{\alpha }}\frac{}{y}\end{array}\right)\mathrm{\Phi }=\lambda \mathrm{\Phi }$$

where

$$\mathrm{\Phi }\left(\begin{array}{c}\varphi _1\\ \varphi _2\end{array}\right)$$

The equation is:

$$\frac{^2\varphi _1}{y^2}+\frac{\alpha }{\sigma }\left(1+\lambda ^2\right)\varphi _1=0$$

The general solution is written

$$\varphi _1=A\mathrm{cos}\left(\kappa y\right)+B\mathrm{sin}\left(\kappa y\right)$$

$$\kappa =\sqrt{\frac{\alpha }{\sigma }}\sqrt{1+\lambda ^2}$$

The second component will be also a combination of harmonic functions:

$$\varphi _2=C\mathrm{cos}\left(\kappa y\right)+D\mathrm{sin}\left(\kappa y\right)$$

The equation can have two independent one-column solutions. The one corresponding to the boundary conditions

$$\mathrm{\Phi }_1\left(y=0;\lambda \right)\left(\begin{array}{c}\varphi _1\\ \varphi _2\end{array}\right)_{y=0}=\left(\begin{array}{c}1\\ 0\end{array}\right)$$

is determined by $`A=1`$, $`B=\sqrt{\alpha /\sigma }\lambda /\left(i\kappa \right)`$ , $`C=0`$ and $`D=\sqrt{\alpha /\sigma }i/\kappa `$. The solution is

$$\mathrm{\Phi }_1(y;\lambda )=\left(\begin{array}{c}\mathrm{cos}\kappa y\lambda i\sqrt{\frac{\alpha }{\sigma }}\frac{\mathrm{sin}\kappa y}{\kappa }\\ i\sqrt{\frac{\alpha }{\sigma }}\frac{\mathrm{sin}\kappa y}{\kappa }\end{array}\right)$$
(69)
The second fundamental solution corresponds to the initial condition;

$$\mathrm{\Phi }_2\left(y=0;\lambda \right)=\left(\begin{array}{c}0\\ 1\end{array}\right)$$

and it results:

$$\mathrm{\Phi }_2(y;\lambda )=\left(\begin{array}{c}i\sqrt{\frac{\alpha }{\sigma }}\frac{\mathrm{sin}\kappa y}{\kappa }\\ \mathrm{cos}\kappa y+\lambda i\sqrt{\frac{\alpha }{\sigma }}\frac{\mathrm{sin}\kappa y}{\kappa }\end{array}\right)$$
(70)
The fundamental solution matrix $`\mathrm{\Phi }`$ is a $`2\times 2`$ matrix having as column the two above solutions. From the fundamental solution matrix we calculate the monodromy matrix

$$M\left(\lambda \right)\mathrm{\Phi }\left(y=d;\lambda \right)$$

The discriminant is the trace of the monodromy matrix.

$$\mathrm{\Delta }\left(\lambda \right)=2\mathrm{cos}\left(d\sqrt{\frac{\alpha }{\sigma }}\sqrt{1+\lambda ^2}\right)$$
(71)
We can proceed to the calculation of the *spectrum*.

### 5.2 The spectrum of the plane wave solution

#### 5.2.1 The spectrum of the unperturbed initial condition (case of plane wave)

We take a set of constants as derived from the numerical calculations that connects the physical conditions to the coefficients of the NSE, $`d=2\pi \left(\frac{a}{L}\right),\alpha =0.399,\sigma =1.452`$. With these constants we calculate the spectrum on the complex $`\lambda `$-plane. Although we can have analytical solutions the graphs below show the difficulty of resolving the details of the variation of the functions. To be noted that the function is intersected with a plane arbitrarly placed at the magnitude $`TrM=5`$, for easier 3d-representation.

This graph shows the presence of roots on the imaginary axis of $`\lambda `$ of the equation $`\left|TrM\left(\lambda \right)2\right|=0`$ with very narrow variation of the function. Actually it is easily to see that for the parameters adopted here, there are *three* roots of this equation, while in the graph we can hardly see two of them. A graphical better way of finding these roots is to look for the intersections of the contours of the two equations:

$`\mathrm{R}e\left(TrM\right)2`$ $`=`$ $`0`$
$`\mathrm{I}m\left(TrM\right)`$ $`=`$ $`0`$
where, again, the other two roots cannot be easily resolved. However the analytical expression allows to find the positions of the roots both on the $`\mathrm{R}e\lambda `$ axis

and on the $`\mathrm{I}m\lambda `$ axis:

By contrast, the variation along the real $`\lambda `$ axis easily evidences the zeroes of the corresponding to the *main spectrum*.

The stability or instability of the Bloch fuctions (*i.e.* the values of the quantity $`m\left(\lambda \right)`$ which is the eigenvalue of the monodromy matrix) is governed by the *discriminant* $`\mathrm{\Delta }\left(\lambda \right)`$. We have $`2\mathrm{\Delta }\left(\lambda \right)2`$ for all real $`\lambda `$ and for $`\lambda =i\alpha `$ with $`1\alpha 1`$. The *main spectrum* corresponds to the values of $`\lambda `$ for which $`\mathrm{\Delta }\left(\lambda \right)=\pm 2`$. For *real* $`\lambda `$, $`\lambda \lambda _{n,\mathrm{r}eal}^0`$ this means

$$\lambda _{n,\mathrm{r}eal}^0=\pm \left[\left(\frac{n\pi }{d\sqrt{\alpha /\sigma }}\right)^21\right]^{1/2}$$

for $`n3.`$ For pure imaginary $`\lambda `$, $`\lambda i\lambda _{n,imag}^0`$

$$\lambda _{n,imag}^0=\pm \left[1\left(\frac{n\pi }{d\sqrt{\alpha /\sigma }}\right)^2\right]^{1/2}$$

for $`n=0`$, $`1`$, $`2`$ (with the values of the parameters $`\alpha `$, $`\sigma `$ and $`d`$). These values are: $`i\lambda _{0,imag}^0=\pm i`$, $`i\lambda _{1,imag}^0=\pm i0.9422`$ and $`i\lambda _{2,imag}^0=\pm i0.7424`$. Since $`d`$ is finite only a finite number of degenerate pairs will appear on the imaginary axis. At $`n=0`$ there are two eigenvalues $`\lambda _0=\pm i`$ . These eigenvalues are *nondegenerate* in the sense that there is only one Bloch function for each of these values of $`\lambda `$. This can be shown directly, constructing the Bloch functions form the matrix of fundamental solutions (69) and (70) with $`\kappa 0`$ (because $`\mathrm{\Delta }\left(\lambda \right)=2`$ implies $`\kappa =0`$) and showing that only one Bloch function can be obtained.

#### 5.2.2 The spectrum of the perturbed initial condition (case of plane wave)

Starting from an initial condition slightly different of the *plane wave’s envelope*

$$u(y,0)u(y,0)+\epsilon h\left(y\right)$$
(72)
we shall have to repeat all steps in the determination of the spectrum. This must be done in a comparative manner, identifying the degenerate eigenvalues $`\lambda _j`$ of the main spectrum which becomes nondegenerate after perturbation. The time behaviour of the exact solution developing from the perturbed initial condition is dominated by the imaginary part of the variables $`W`$ which are obtained on the basis of the new (perturbed) quantities $`\mu _j`$ and this can be calculated.

Using the new initial condition Eq.(72) we shall calculate the new main spectrum, i.e. the new values of the parameter $`\lambda `$ for which the discriminant is $`\pm 2`$. All the spectrum of the simple initial condition $`u(y,0)`$ is perturbed by order-$`\epsilon `$ quantities. The main effect of this perturbation is that the degeneracies are broken. In order to construct the exact solution corresponding to the new position of the nondegenerate $`\lambda `$ values, we need to specify the cycles on the two-sheeted Riemann surface. The steps are as usual and consists of first choosing the cycles on the Riemann surface and the base of the holomorphic differentials; then, compute the matrices of periods (the matrix of $`a`$ periods, and the matrix of the $`b`$ periods); compute the $`\tau `$ matrix; construct the Abel map, which maps the points of the compact two-sheeted Riemann surface onto the Jacobi torus; find the wavenumbers and the frequencies in the Jacobi manifold; make the Jacobi inversion, using the theta functions;

##### Finding the base of the holomorphic differentials

One of the eigenvalue pair of the unperturbed system is $`(i,i)`$. Consider now that under the perturbation a degenerate eigenvalue $`\lambda _j^0`$ splits into the conjugated pair $`(\lambda _k^0\epsilon _k,\lambda _k^0+\epsilon _k)`$. In Ref.() it was constructed a new basis as a linear combination of the canonical basis:

$$dU_j=\frac{1}{2\pi i}\sqrt{1+\left(\lambda _j^0\right)^2}\frac{_{mj}\left(\lambda \lambda _m^0\right)d\lambda }{R\left(\lambda \right)},\text{for}j=1,2,\mathrm{},N$$

where

$$R\left(\lambda \right)=\sqrt{1+\lambda ^2}\underset{k=1}{\overset{N}{}}\left[\left(\lambda \lambda _k^0\epsilon _k\right)\left(\lambda \lambda _k^0+\epsilon _k\right)\right]^{1/2}$$

To understand the choice of basis it is useful to look at the limit $`\epsilon 0`$ where the basis of holomorphic differential becomes

$$\underset{\epsilon 0}{lim}dU_j=\frac{1}{2\pi i}\sqrt{1+\left(\lambda _j^0\right)^2}\frac{1}{\sqrt{1+\lambda ^2}}\frac{d\lambda }{\left(\lambda \lambda _j^0\right)}$$

This expression shows that, at the limit $`\epsilon 0`$, each differential “sees” only one pole, at $`\lambda _j`$ (which is a *double point* at this limit), and the square-root *branch* points $`\pm i`$.

The loop (cycle) $`a_j`$ encircles the pole $`\lambda _j^0`$ and the period is

$$\underset{\epsilon 0}{lim}A_{kj}=\underset{\epsilon 0}{lim}_{a_j}𝑑U_k=_{a_j}\underset{\epsilon 0}{lim}𝑑U_k=\delta _{kj}$$

This is because the cycle $`a_j`$ surrounds the pole $`\lambda _j^0`$ or it does not surrounds any other pole.

We find that, under the *limit operation* $`\epsilon 0`$, the matrix $`𝐀`$ goes over to the identity matrix to the first order in $`\epsilon `$, $`O\left(\epsilon \right)`$ and that the matrix $`𝐁`$ becomes

$$\underset{\epsilon 0}{lim}B_{jk}=\tau _{jk}+O\left(\epsilon \right)$$

To see what happens with the $`𝐁`$ matrix in this limit, we must examine the *off-diagonal* terms and the diagonal terms.

The *off-diagonal* terms are well-behaved and can be integrated analytically:

$$\underset{\epsilon 0}{lim}B_{jk}=_{b_j}\underset{\epsilon 0}{lim}𝑑U_k\text{for}jk$$

and

$$\tau _{kj}=\frac{1}{2\pi i}\sqrt{1+\left(\lambda _j^0\right)^2}_{b_j}\frac{1}{\sqrt{1+\lambda ^2}}\frac{d\lambda }{\left(\lambda \lambda _k^0\right)}+O\left(\epsilon \right)$$

The contour $`b_j`$ passes through the middle of the segment connecting the two new $`\lambda _j`$ and goes around the pole at $`\lambda =i`$, for example. If this contour does not surround the pole at $`\lambda _k^0`$ then the integral is elementary. If the contour intersecting the cut from $`i`$ to infinity encircles the pole at $`\lambda _k^0`$ then it may be chnaged such that to intersect the other cut, connecting the branch point $`+i`$ to infinity, and so it does not encircle any pole. In conclusion since the off-diagonal elements of the matrix $`\tau `$ depend (to $`0\left(\epsilon \right)`$) only of the positions of the original double points, which are determined only by the parameter $`d`$, it results that: the off-diagonal terms of the $`\tau `$ matrix does not carry any information on the *perturbed initial condition*. Only the diagonal terms of the $`\tau `$ matrix are affected by the initial condition.

The *diagonal* terms of the $`\tau `$ \- matrix are *singular* at the limit $`\epsilon 0`$. Consider

$$\tau _{kk}=\frac{1}{2\pi i}\sqrt{1+\left(\lambda _j^0\right)^2}_{b_k}\frac{1}{\sqrt{1+\lambda ^2}}\frac{d\lambda }{\sqrt{\left(\lambda \lambda _k^0\right)^2\epsilon _k^2}}+O\left(\epsilon \right)$$

Suppose that

$$\lambda _k^0\text{is on the }\text{imaginary axis}$$

Then we have

$$\tau _{kk}=\frac{1}{2\pi i}\sqrt{1+\left(\lambda _j^0\right)^2}\left(2_i^{\lambda _k^0}+\frac{1}{2}_{a_k}\right)\frac{1}{\sqrt{1+\lambda ^2}}\frac{d\lambda }{\sqrt{\left(\lambda \lambda _k^0\right)^2\epsilon _k^2}}+O\left(\epsilon \right)$$

$$\tau _{kk}=\frac{1}{2}+\frac{1}{\pi i}\sqrt{1+\left(\lambda _j^0\right)^2}_i^{\lambda _k^0}\frac{1}{\sqrt{1+\lambda ^2}}\frac{d\lambda }{\sqrt{\left(\lambda \lambda _k^0\right)^2\epsilon _k^2}}+O\left(\epsilon \right)$$

The integral can be made with *elliptic functions* but an approximation is (near $`\epsilon _k=0`$)

$$\tau _{kk}=\frac{1}{2}\frac{i}{\pi }\mathrm{ln}\left|\epsilon _k\right|+O\left(1\right)+O\left(\epsilon _k\right)$$

Consider now

$$\lambda _k^0\text{ is on }\text{real axis}$$

then the integral defining the diagonal elements of the matrix $`\tau `$ (the integral along the $`b_k`$ cycles) is

$$_i^{\lambda _k^0}=_i^0+_0^{\lambda _k^0}$$

where the second integral is along the real axis.

The first integral is done with the residuu theorem after taking safely the limit $`\epsilon _k=0`$ and gives a finite imaginary contribution which is well-behaved for this limit. The second can be done and the result is

$$\tau _{kk}=\frac{i}{\pi }\mathrm{ln}\left|\epsilon _k\right|+iO\left(1\right)+O\left(\epsilon _k\right)$$

We then conclude that for $`\epsilon _k0`$ the $`b`$ \- period matrix $`\tau `$ becomes logarithmically singular on the main diagonal.

The main conclusion is: (1) The real-$`\lambda `$ axis degenaracies have $`\mathrm{\Omega }_j`$ real so they are linearly stable. The imaginary axis degeneracies have $`\mathrm{\Omega }_j`$ imaginary so they are linearly unstable.

This conclusion can have important practical effects: in our case the degenerate eigenvalues on the *imaginary* $`\lambda `$-axis split and induce instability of the plane wave solution, *i.e.* of the envelope. At the same time the solutio evolves toward one of the stable (soliton-like) solutions. Still much work is required to examine the stability of that solution to additional perturbation. For example, the main spectrum of the solution of the type $`\mathrm{sec}h(y)`$ can be obtained from the discriminant whose form is very complicated, with many purely imaginary (dangerous) eigenvalues.

## 6 The growth of the perturbed solution as a source of poloidal asymmetry and spontaneous plasma spin-up

### 6.1 The minimal condition for instability of the poloidal rotation

We have found that the space-uniform solution (poloidal symmetry) is unstable and the turbulence will self-modulate according to one of the possible soliton-like solution. In a time picture, we can say that the poloidal symmetry will be left and the turbulence amplitude will develop a poloidal asymmetry which grows exponentially in time. The onset of the poloidal asymmetry will induce an asymmetric rate of transport on the poloidal direction . When the particle diffusion is poloidally asymmetric and the local particle confinement time is shorter than the damping time of the poloidal rotation, the poloidal rotation is unstable. The plasma can spontaneously spin-up. The growth rate of the asymmetry (which is given by the rate of departure from the uniform solution toward a soliton-like solution) determines the growth rate of the torque applied to plasma in the poloidal direction but the relation is not direct, depending essentially of the momentum balance in the poloidal direction. Any torque acting to generate a poloidal rotation must overcome the neoclassical magnetic pumping, which is an efficient damping mechanism. A balance of these forces depends on the comparison between the rates of (1) asymmetry growth, the radial variation of the asymmetry-induced terms, on one hand, and (2) rate of damping through magnetic pumping, on the other hand.

For the simplest case of the space-uniform envelope solution, we have the instability has the *linear* growth rate

$$i\gamma ^{env}=i\mathrm{I}m\mathrm{\Omega }=i\left|\sqrt{\alpha }k\left(\alpha k^24p^2\sigma \right)^{1/2}\right|$$

where $`k<2p\sqrt{\sigma /\alpha }`$ is here the wavenumber of the perturbation of the envelope. In the case that has served as example in the preceeding Section, we have $`p=\sqrt{t_0\omega _{unif}^{phys}/\left(2\sigma \right)}1.86`$ and $`k<7.08`$. Taking for example $`k=3`$ we obtain $`\gamma ^{env}=7.67`$ or $`\gamma _{phys}^{env}=7.67\times 10^4`$ ($`s^1`$).

### 6.2 The soliton developing from the initial perturbation

Since we know that the uniform state can be unstable we have to investigate the problem of the evolution of the system after a perturbation is applied and look for stable solutions. The solitons of this system can be realisable in real conditions since they are intrinsically more stable. We will look for soliton-like solution for the uncoupled equations, taking a single NSE like (64). The time dependence can be

$$\varphi =\mathrm{exp}\left(i\mu ^2\alpha t\right)w\left(x\right)$$

giving the reduced equation

$$\frac{^2w}{x^2}\mu ^2w+\frac{\sigma }{\alpha }w^3=0$$

where $`\mu `$ is a real constant. This equation has the following solution

$$w=\sqrt{\frac{2\alpha }{\sigma }}\mu \mathrm{sec}h\left(\mu x\right)$$

The constant $`\mu `$ can be related to the modulation wavenumber whose upper bound is (68). Taking $`\mu 3`$ and $`\alpha =0.399`$ we obtain the frequency of steadily oscillating soliton $`\omega _{osc}^{env}\mu ^2\alpha 1.2`$, which is slower than the growth rate of the deformation, $`\gamma ^{lin}=7.67`$ obtained above (all in absolute units). Since the rate of change of the poloidal velocity (which will be calculated below) has comparable magnitude with $`\omega _{osc}^{env}`$ we conclude that the rotation has the time to develop before the oscillation reverses the sense of the applied torque. The maximum amplitude of the deformation is $`\sqrt{2\alpha /\sigma }\mu =2.224`$. Large smooth deformations (small $`\mu `$) with less localised profile will have even slower steady oscillatory motion and however significant amplitude, and we can expect they are favored in real situations.

### 6.3 The poloidal asymmetry of the diffusion fluxes

Our approach, which is adapted to the investigation of the stability of the turbulence envelope, can only provide a *model* of the turbulent field, according to Eq.(3.1). This is composed of the primary modes propagating in the poloidal direction with amplitude modulated by the slowly varying envelopes $`A_{1,2}`$ obeying the system of coupled Nonlinear Schrödinger Equations (30) and (31). In the radial direction the field has a variation given by the functions $`\phi _{1,2}\left(x\right)`$ obtained by solving the equations (22). We will need to reconstitute the turbulent field in order to calculate the modified diffusion coefficient. We choose the following set of parameters:

$`a`$ $`=`$ $`1\left(m\right)`$
$`B_T`$ $`=`$ $`3.5\left(T\right)`$
$`T_e`$ $`=`$ $`1\left(KeV\right),T_i=1\left(KeV\right)`$
$`n`$ $`=`$ $`10^{20}\left(m^3\right)`$
$`L_n`$ $`=`$ $`0.1\left(m\right)`$
$$U_0=6\times 10^5\left(m/s\right)$$

The extension on the radial direction of the mode is taken $`L=0.1\left(m\right)`$ and this is also the normalising length. Other quantities are: normalising time $`t_0=0.16\times 10^6\left(s\right)`$ and $`\mathrm{\Omega }_i=0.335\times 10^9\left(s^1\right)`$, the ion-cyclotron frequency. From this it results the two parameters of the barotropic form of the ion mode, $`\beta =111.75`$ and $`\epsilon =0.476\times 10^2`$. We take the following frequencies for the primary poloidal modes: $`\omega _1=0.1451\times 10^6\left(s^1\right)`$ and $`\omega _2=0.1931\times 10^6\left(s^1\right)`$. From the solution of the equations (22) the two eigenvalues result: $`k_1=43.51\left(m^1\right)`$ and $`k_2=56.81\left(m^1\right)`$. Solving the series of differential equations described in detail in Section 3 we obtain the following constants in the coupled Nonlinear Schrödinger Equations: $`\alpha _1=0.1756\times 10^2`$ , $`\alpha _2=0.2278\times 10^2`$ and respectively $`\sigma _1=0.18272`$ , $`\sigma _2=0.5110`$. The frequency of the uniform oscillations, Eq.(67) is taken $`p=0.0675`$ giving the maximum poloidal number for the large scale instability, Eq.(68), $`\overline{m}=\left[5.89\right]=5.`$ The growth rate of the modulational instability is $`\gamma =4175\left(s^1\right)`$. The expression of the field becomes

$`\varphi (r,\theta )`$ $`=`$ $`\varphi _0[(1+A_1\left(\theta \right)\mathrm{sec}h(y))\mathrm{cos}(k_1y\omega _1t+\delta _1)\phi _1\left(r\right)`$
$`+(1+A_2\left(\theta \right)\mathrm{sec}h(y))\mathrm{cos}(k_2y\omega _2t+\delta _2)\phi _2\left(r\right)]`$
where $`\varphi _0`$ is explained below. The initial phases $`\delta _{1,2}`$ are arbitrary and in the numerical calculations are taken as random quantities with uniform distribution on the interval $`(0,1)`$.

We notice that the equations (22) fulfilled by $`\phi _{1,2}`$ are homogeneous and the amplitudes of the two functions are arbitrary. We fix the amplitude $`\varphi _0`$ with the condition that the *uniform* profile of the averaged square of the potential of the fluctuating field obtains a typical value of the diffusion coefficient, $`D_0\left|\varphi _0\right|^2`$, which is taken $`D_01\left(m^2/s\right)`$. The following two figures represent the fluctuating potential, resulting from the combination of the factors described above.

A simple way to take into account the effect of the turbulence amplitude modulation on the diffusion coefficient is to consider that locally the diffusion coefficient is proportional with the quantity $`\left|\varphi \right|^2`$, where the averages should be taken over a statistical ensemble of realizations of the fluctuating field. This requires a carefull examination when applied to numerical calculations. We have to remember that the diffusion is a statistical process where a test particle interacting with fluctuating background field performs random displacements whose mean square grows linearly in time. The random displacements have to belong to a Gaussian statistical ensemble and this requires that all kind of elementary experiences (i.e. displacements) to be represented with the corresponding probability. This is practically irrealisable in numerical calculations since the statistical ensemble of the stochastic field configurations is necessarly finite, i.e. incomplete as a Gaussian ensemble. Any average will be made on subensembles, more or less incomplete. The collection of field configurations are obtained from time *and* from space realizations of the field. The test particle should have the possibility to manifest the *linear* mean square displacement specific to the Gaussian ensembles and this rises the question on what is the minimal spatial extension on which the average should be done. We need in any case a space region larger than the correlation length of the fluctuating field which can be taken as the inverse of the characteristic wavenumber in the field (minimum of $`k_1`$, $`k_2`$).

The profile of the quantity $`D_0\left|\varphi \right|^2`$ is represented in the figure below. The exemple chosed for these graphical representation is characterized by the imposed constrain to obtain a rapid variation of the functions $`\phi _1`$ and $`\phi _2`$ on the radial $`r`$-interval, with the intention to simulate short wavelength ion turbulence, frequently observed in the numerical simulations of the Ion Temperature Gradient driven turbulence. This is achieved by asking the numerical procedures which solves the Schrödinger-type equations (22) to look for high order $`\kappa `$ of the eigenvalues. The pictures represented in this section correspond to $`\kappa =19`$ for both eigenfunctions. This generates a potential dominated by small spatial scales and consequently, the fast variation of the local value of the diffusion coefficient. We take, in the case $`\kappa =19`$ a space region containing at least three complete oscillations of the typical wavelength. This should allow at least several random steps in the typical displacement of the particle. When we do this (taking arbitrarly a factor of few units larger than the correlation length) we note that the diffusion coefficient exhibits space variations. These are natural and are inherited from the space dependence of the radial eigenfunctions $`\phi _{1,2}\left(r\right)`$. Since in real situations the radial eigenfunctions are superposed with random initial phases we need to enlarge the region of averaging over most of the radial interval considered, up to a level of saturation of these spatial variation which can be considered convenient for the estimation of the poloidal torque.

### 6.4 The poloidal torque applied on plasma

We use the equations describing the plasma poloidal and toroidal rotation derived in , :

$$\frac{n}{t}+\frac{1}{r}\frac{}{r}\left(rn\overline{v}_r\right)=0$$

$$\frac{}{t}\left(nV_\phi \right)+\frac{1}{r}\frac{}{r}\left[rn\left(V_\phi \overline{v}_rqV_\theta \stackrel{~}{v}_r\right)\right]=0$$

$$\frac{}{t}\left[V_\phi +\mathrm{\Theta }\left(1+2q^2\right)V_\theta \right]+\overline{v}_r\frac{V_\phi }{r}\stackrel{~}{v}_r\frac{}{r}\left(qV_\phi \right)+\text{magnetic pumping}=0$$

The definition of the variables is $`n=n\left(r\right)`$, $`nV_\phi \left(r\right)nv_\phi \frac{R}{R_0}`$ and $`V_\theta \left(r\right)=v_\theta \frac{R}{R_0}`$, where $`v_\phi `$ and $`v_\theta `$ are local toroidal and poloidal plasma velocities in the magnetic surface. The magnetic field is

$$𝐁=(0,\mathrm{\Theta }\left(r\right),1)\frac{B_0\left(r\right)}{R/R_0}$$

with $`R=R_0+r\mathrm{cos}\theta `$ and $`q=\frac{r}{R_0}\frac{1}{\mathrm{\Theta }}`$. The surface average is made according to the formula

$$f=\frac{d\theta }{2\pi }\frac{R}{R_0}f$$

In the collisional regime the *magnetic pumping* term determined by parallel viscosity is less effective and the spin up can be expected. The decay rate of the poloidal rotation by this mechanism is :

$$\gamma _{MP}\frac{3}{4}\left(1+\frac{1}{2q^2}\right)^1\left(\frac{l}{qR}\right)^2\nu _{ii}$$
(73)
where $`l`$ is the mean-free path. The velocities $`\stackrel{~}{v}_r`$ and $`\overline{v}_r`$ arise from the diffusive flux in the radial direction, to which is associated the local radial velocity $`V_r`$. We define the two velocities

$$\overline{v}_r=V_r$$

$$\stackrel{~}{v}_r=2\mathrm{cos}\theta V_r$$
(74)
If there is no poloidal variation of the local radial velocity $`V_r`$ the average in (74) is that of the trigonometric function $`\mathrm{cos}\theta `$ and the result is zero. The average will not be zero if there is a poloidal variation of the diffusion-generated particle radial velocity, $`V_r\left(\theta \right)`$. We can introduce the modified diffusion coefficient $`\stackrel{~}{D}`$ related to the nonsymmetry,

$$n\stackrel{~}{v}_r=\stackrel{~}{\mathrm{\Gamma }}_r=\stackrel{~}{D}\frac{n}{r}$$
(75)
with

$$\stackrel{~}{D}=D_0\left|\varphi \right|^2_{st}2\mathrm{cos}\theta _{pol}$$
(76)
where the first averaging is statistical and the second is on the poloidal surface. The equation for the plasma poloidal velocity is

$$\mathrm{\Theta }\left(1+2q^2\right)\left(\frac{V_\theta }{t}+\gamma _{MP}V_\theta \right)+qV_\theta \frac{1}{nr}\frac{}{r}\left(nr\stackrel{~}{v}_r\right)=0$$
(77)
The average velocity Eq.(75) still dependes on the radial coordinate via $`D_0`$, the density $`n\left(r\right)`$ and its derivatives and, more important, via the space dependence of the potential $`\phi \left(x\right)`$ left after the statistical averaging needed for the calculation of the diffusion coefficient (76). The following quantity must be numerically calculated

$$\gamma \frac{q}{\mathrm{\Theta }\left(1+2q^2\right)}\frac{1}{nr}\frac{}{r}\left(nr\stackrel{~}{v}_r\right)=\frac{q}{\mathrm{\Theta }\left(1+2q^2\right)}\frac{1}{nr}\frac{}{r}\left(rD_0\left|\varphi \right|^2_{st}2\mathrm{cos}\theta _{pol}\frac{n}{r}\right)$$

and this directly provides the rate of change of the poloidal velocity due to the asymmetry. This must be compared to the magnetic pumping term (73) and we expect rotation instability when

$$\gamma \gamma _{MP}>0$$

In this formula the soliton profile of the poloidal asymmetry is considered saturated and the motion of the solitons (which would give periodic variation of the torque) is neglected. As is shown in Ref. the high positive second derivative of the density profile, combined with a dependence of the diffusion coefficient $`D_0`$ on the shear of the poloidal velocity are the favorable circumstances for the asymmetry-induced torque to overcome the magnetic pumping term and a spontaneous spin up to occur. We do not take into account the variation of $`D_0`$ with the velocity shear, but notice that the space variation of (76) also contribute to obtain a positive value for $`\gamma `$ sufficient for being comparable to $`\gamma _{MP}`$ and even to be larger. The most important contribution comes however from the positive-valued second density gradient, which can be seen as an indication that the spin-up can be expected in such region on the minor radius.

The values represented in the figure have been calculated using a gaussian profile for the density. The decay due to magnetic pumping is of the order of $`\gamma _{MP}12\left(s^1\right)`$. We note that the torque applied to plasma shows a radial variation essentially inherited from the space variation of the statistically averaged potential fluctuations, which in general are not monotone functions of the radius on time scales of several eddy turn-over times. These spatial variations of the diffusion coefficient generate the profile of the torque which change the shear of the seed velocity and provide a characteristic which is expected for zonal flows. Even outside the domain of the parameters where the spin up is efficient compared to the magnetic pumping, the soliton asymmetry and the Stringer torque can contribute to the action of other mechanisms to generate poloidal rotation.

## 7 Discussion

We have examined the problem of stability of the envelope of the ion wave turbulence. Starting from a simple model for the ion dynamics and retaining the full nonlinear contribution of the ion polarization current we have derived an equation of the same form as the *barotropic* model of the atmosphere. The multiple space-time scale analysis shows that the envelope of the turbulence results from a set of coupled cubic nonlinear Schrodinger equations. The simplest solution is the space-constant envelope corresponding to a symmetric average amplitude of the turbulence on the poloidal circumference. In order to examine the stability of this solution we first review the method of inverse scattering transform in the geometric-algebraic framework, able to obtain the exact solution developing from a certain class of initial conditions (finite-band potentials). In this framework appears clear the origin of the sensitivity to the perturbations, connected to the topology of a hyperelliptic Riemann surface. The uniform solution is unstable and evolves to one of the soliton-type solution. The nonuniformity of the turbulence average profile generates a torque on the plasma, which may be sufficient to overcome the poloidal damping.

The simple physical model we have employed may be partly justified by the observation that the analysis is essentially independent of the other possible sources of poloidal asymmetry: the neoclassical variation of parameters on the magnetic surface, the toroidal effect on the nonlinearity represented as mode coupling, the effect of the particle drifts on the linear and renormalized propagators, etc. The cubic nonlinear Schrodinger equation has many solutions on the periodic domain and the relevance for the self-modulation of the envelope of the turbulence should be examined. Here we only considered the uniform envelope since its instability can be sufficient to induce plasma poloidal rotation. This rotation has an origin different of the usual sources of rotation: neoclassical, ion-loss, Reynolds stress or momentum injection. Our approach is essentially mathematical and so it operates at a more formal level. The physical origin of the instability seems to be the tendency of the system (averaged amplitude of the turbulence) to evolve to profiles characterised by the balance of the dispersion and the nonlinearity. What remains to be examined is the quantitative constraint on the minimal initial perturbation needed to produce a particular evolution: for a weak perturbation the solution may be simple “radiation” while for initially larger amplitude a soliton emerges ensuring higher robustness and a localised form. It is in this case that the asymmetry can induce plasma rotation. This question can be answered by examining the exact solution, with specific numerical tools for the determination of the main spectrum. It might be possible that the Riemann curve is not a hyperelliptic surface but a manifold close to a finite-genus curve.

Aknowledgements. The authors gretefully aknowledge the useful discussions and encouragements form J. H. Misguich and R. Balescu. The hospitality of the Département de Recherche sur la Fusion Controlée (Cadarche, France) is aknowledged. This work has been partly supported by the NATO Linkage Grant PST.CLG.977397.