数理解析研究所講究録
第 1606 巻 2008 年 1-17

1Massive point-vortex simulation using
special-purpose computer, MDGRAPE-3
YATSUYANAGI, Yuichi
Faculty of Education, Shizuoka University*

1

Introduction

Our final goal is to understand characteristics of negative temperature state appearing
in dynamics of two-dimensional point vortex system [1].
We began a simulation of the point vortex system to understand results of nonneutral
plasma experiments obtained by Kiwamoto group at Kyoto University. It can be shownthat the electron motion in the plane perpendicular to the strong, uniform magnetic fieldis given by two-dimensional Euler equation. In the same sense, the governing equation for
the two-dimensional guiding-center plasma that is a collection of charged rods is identicalto that of two-dimensional, inviscid, incompressible fluid.
One of the most remarkable features of the point vortex system is the negative temperature state first introduced by Onsager [2]. Two-dimensional point vortex system confinedin a finite area has limited phase space volume, so that in the limit of infinite system energy, the density of state approaches zero and has at least one peak at some energy.
Following the statistical definition of temperature, $1/T=dS/dE=k_{B}d\ln W/dE$ , whereis temperature, is entropy, $E$ is energy,
is Boltzmann constant, and $W$ is densityof state, the sign of temperature changes if system energy is larger than
. Note that
there is no concept of negative absolute temperature in the thermodynamic sense. However, as is mentioned by Onsager, there are many phenomena that can be well-describedby assuming negative statistical temperature, even if the thermodynamic temperature
is positive [2]. We consider that characteristic features that originate in the negative
temperature may appear in the dynamics of the point vortex system.
Wide variety of research effort has been devoted to understand the negative temperature state in the context of two-dimensional turbulence. The condensation of thesame-sign vortices in the rectangular domain is first demonstrated numerically by Joyce$E_{0}$$T$$S$$k_{B}$$E_{0}$’eyyatsu@ipc.shizuoka.ac.jp

2and Montgomery [3]. The equilibrium state in the negative temperature is discussed by
Pointin and Lundgren [4]. The density of states is obtained by analytically by Seyler [5]and numerically by Buehler 6] and Johnson [7]. However, these simulation results arevery restricted ones in the number of particles and resolution due to the limited calculation power. Thus, we started the point vortex simulation on the negative temperaturestate using recent high performance computers.
In the point vortex simulation, time evolution of the point vortices is traced by thevorticity equation. To update the position of the point vortex, the flow velocity at eachpoint vortex must be calculated. This calculation is called the Biot-Savart integral, andits calculation cost is proportional to the square number of the total point vortices. There
are some different solutions for reducing the calculation cost:$[$$\bullet$$\bullet$$\bullet$Use a fast algorithm.
Use a PC cluster.
Use a supercomputer.

For PC-based simulations, the most practical solution may be to use a fast algorithm,
for example, tree code and Vortex-In-Cell approximation. However, the integral precision
is also usually reduced. For our case, this is a severe problem as the target system innegative temperature is high energy one and to trace the trajectory of the motion ofvortices, the velocity at each point vortex should be evaluated exactly at the position of
the point vortex.
The second candidate, PC cluster, has become a very popular solution for massive
simulations. Even a university laboratory can bear the cost and the space to maintaina small cluster. However, in our case, total calculation speed necessary is approximately
100 GFlops over and more than 50 PCs are required to achieve this speed. So, we decided
to choose the third solution, using a supercomputer called MDGRAPE series. It was
originally designed for molecular dynamics simulation. By using a capability of arbitrarychoice of Green’s function, Biot-Savart integral also can be calculated by it. It provides
350 GFlops calculation speed. This means approximately 100 times faster than a usualPC.
By the massive simulations using MDGRAPE series, two-dimensional point vortex
system is examined and characterized by temperature. At first, we examine the densityof states if it has a peak. Equilibrium distributions are obtained by time asymptotically.
We find that slope of the energy spectrum in the intermediate wave number depends oncorrelation function is used for detailed analysis for thethe system energy.
equilibrium distribution.
The organization of this paper is as follows. In \S 2, the prelude of this research isreviewed. In \S 3, the negative temperature state in the point-vortex system is reviewed.
$Tw\mathfrak{c}\succ body$

3

Figure 1: Schematic of confinement geometry is shown. Strong axial magnetic field andedge electrostatic potential confines electrons. They rotate in $ExB$ direction.
In \S 4, special-purpose supercomputer, MDGRAPE-2 andare introduced. In \S 5, weshow simulation results on the point vortex system characterized by the temperature.
$- 3$

2

Prelude of this research: Nonneutral plasma experiment

We started numerical simulations of the point vortex system to understand the resultsof nonneutral plasma experiments carried out by Kiwamoto group at Kyoto University.
Electrons are confined in a circular cylinder axially by an electrostatic potential andradially by a strong uniform magnetic field. Schematic configuration of confinement isshown in Fig. 1. Typical time scales of the electron motion satisfy$\tau_{c}\ll\tau_{ExB}\ll\tau_{bounce}$

,

(1)

where
denotes a period for an electron cyclotron motion,
is a time scale for$ExB$ motion, andis a time scale for an electron to bounce between the ends by
the electrostatic potential. As long as the radius of the cyclotron motion is small enough,the electron motion can be represented by the motion of the guiding center. In this limit,
the equation of the electron motion in the two-dimensional plane perpendicular to themagnetic field is identical to the Euler equation for inviscid, incompressible fluid.
Let the direction of the strong magnetic field be ,$\tau_{ExB}$$\tau_{c}$$\tau_{bounce}$$z$$B=B_{0}\hat{z}$.

(2)Equation of motion of the guiding center is given by$m_{e} \frac{dv}{dt}=-e(E+v\cross B)$.

(3)

By time averaging Eq. (3) over a cycle, the drift velocity of the guiding center is givenby the usual $E\cross B$ drift$v= \frac{E\cross B}{|B|^{2}}$

.

(4)

4

Incorporating electrostatic potential induced by the electrons themselves(5)$E=-\nabla\phi$and Eq. (2) into Eq. (4), we obtain$v= \frac{1}{B_{0}}\hat{z}\cross\nabla\phi$.

(6)Thus, the self electrostatic potential behaves like a stream function for a two-dimensional
flow. Rotation of Eq. (6) gives the vorticity$\omega_{z}\hat{z}=\nabla\cross v=\frac{\hat{z}}{B_{0}}\nabla^{2}\phi=\frac{en}{\epsilon_{0}B_{0}}\hat{z}$

where

$n$

,

(7)

is the number density of the electron. In the last term, Poisson equation is used$\nabla^{2}\phi=\frac{en}{\epsilon_{0}}$.

(8)Equation (7) tells that the vorticity is proportional to the number density of the electron.
By these relations, the governing equation of guiding-center plasma is found to be identicalto the Euler equation. More complete description for the electron confinement is available
in Ref. [8].

3Negative temperature state appearing in two dimensional point vortex system

Large-scale, long-lived vortices are commonly observed in nearly two-dimensional flow.
It can also be seen in the nonneutral plasma experiment. Initially doughnut-shaped
electron distribution is deformed into several clumps by the diocotron instability, or morefamous words, Kelvin-Helmholtz instability. The number of clumps produced initially isunderstood by the linear theory 9 . Shortly, merger continuously occurs and the numberof clumps gradually decreases in time. In this process, each clump size increases, whichindicates the energy transfer from a short scale to a large scale, i.e., the inverse energy
cascade.
For decades, various investigations have been made to understand the inverse energy
cascade [10, 11]. Onsager introduced a concept of “negative temperature” for the twodimensional point vortex system in 1949 to understand the large-scale vortex formation[2]. When temperature $T$ is positive, vortices of opposite sign will tend to approach eachother. On the other hand, when $T$ is negative, vortices of the same sign will tend tocluster and preferably form stronger vortex. This statistical tendency of the same sign
to cluster in the negative temperature region is clear from a description by a canonical
$[$

$]$

5

energy distribution proportional to $\exp(-\beta E)$ . Negative corresponds to reversing thesign of the interaction, where the same-sign vortices statistically attract and opposite ones$\beta$repel 12].
$[$

Here we explain a target point vortex system, including the basic concept of negative
temperature. We consider a point vortex system consisting of $N$ positive and $N$ negative
point vortices with circulation, respectively, bounded by a circularconst. andwall of radius . Position vector and circulation of the i-th point vortex is given by
and , respectively. Constants of motion are Hamiltonian and inertia:$\Gamma_{0}=$$-\Gamma_{0}$$R$$r_{i}$$\Gamma_{i}$

$H$

$=$

$- \frac{1}{4\pi}\sum_{i}^{N}\sum_{j\neq i}^{N}\Gamma_{i}\Gamma_{j}\ln|r_{i}-r_{j}|+\frac{1}{4\pi}\sum_{i}^{N}\sum_{j}^{N}\Gamma_{i}\Gamma_{j}\ln|r_{i}-\overline{r}_{j}|-\frac{1}{4\pi}\sum_{i}^{N}\sum_{j}^{N}\Gamma_{i}\Gamma_{j}\ln\frac{R_{/}}{|r_{j}\}}9)$

$I$

$=$

$\sum_{i}^{N}\Gamma_{i}|r_{i}|^{2}$(10)

The last term in Eq. (9) is introduced to make the value of stream function at the circularwall zero. Equations of motion for each point vortex are given by$\Gamma_{i}\frac{dx_{i}}{dt}=\frac{\partial H}{\partial y_{i}}$,$\Gamma_{i}\frac{dy_{i}}{dt}=-\frac{\partial H}{\partial x_{i}}$

,

(11)

or explicitly.

(12)$\frac{dr_{i}}{dt}=-\frac{1}{2\pi}\sum_{j\neq i}^{2N}\Gamma_{j}\frac{(r_{i}-r_{j})x\hat{z}}{|r_{i}-r_{j}|^{2}}+\frac{1}{2\pi}\sum_{j}^{2N}\Gamma_{j}\frac{(r_{i}-\overline{r}_{j})\cross\hat{z}}{|r_{i}-\overline{r}_{j}|^{2}}$The wall effect is introduced by the image vortex located at
The inverse temperature is statistically defined by$\overline{r}_{i}=R^{2}r_{i}/|r_{i}|^{2}$.

$\beta$$\beta=\frac{dS}{dE}=k_{B}\frac{d\ln W(E)}{dE}$

.

(13)

In usual cases, $W(E)$ increases as $E$ increases, and the slope $dS/dE$ never changes thesign, i.e., always positive (Fig. 2). On the other hand, if the total phase space volume
is finite, the total number of states is limited, and the asymptotic value of the density ofstate at infinite energy must be zero. Then the density of state should have at least apeak at some energy. So that, in the region of $E>E_{0}$ , temperature is negative.
Onsager considered that theand coordinates of each point vortex are canonicalconjugates due to the similarity of the equations of motion of the point vortices (11) tothe common Hamilton equation. Namely, the phase space of the point vortex system is
identical to the configuration space. Assuming that the vortices are confined in an area
, total phase space volume is given by$E_{0}$$x_{i}$$y_{i}$$A$$/d\Omega=/dx_{1}dy_{1}\cdots dx_{2N}dy_{2N}=(/dxdy)^{2N}=A^{2N}$

.

(14)

6

Figure 2: Usual relation between the density of state and the system energy is shown.

Figure 3: The density of state has at least a peak when the total number of states islimited. In this case, negative temperature state appears when $E>E_{0}$ .

7Figure 4: MDGRAPE-3 is a PCI-X card.
Thus, Onsager concluded that the negative temperature state appears in the two-dimensional
point vortex system confined in a finite area. Since then, many simulation results havebeen presented to understand the two-dimensional point vortex system in statistical mechanical framework. However, numerical simulation on point vortex system demands
relatively large calculation power. We noticed that there were few large scale simulationon the point vortex system in negative temperature state, so that we decided to investigate
it by on-hand supercomputer, MDGRAPE-2 and MDGRAPE-3.

4

An accelerator for vortex simulations, MDGRAPE3

To speed up the calculation in central processing unit (CPU), it is very efficient to
restrict the types of calculation and completely implement them as a wired-logic device[13]. An extreme example is a special-purpose computer that can calculate usually only
one kind of calculation. One of the most famous special-purpose computers may be
GRAPE hardware that was at first developed at University of Tokyo [14]. The type of
calculation is limited to the gravitational force forstars. So that, GRAPE can be
regarded as a fast Poisson solver for a three-dimensional system. As a branch of GRAPE,MDGRAPE series was developed. It has a capability to integrate Poisson equation withdifferent Green function from the three-dimensional one, for example, Van der Waals
force, two-dimensional Coulomb force, Biot-Savart integral.
Calculation time for each time step is measured in the point vortex simulation withtime development by 4th order Runge-Kutta method. The simulations are executed as(a) a single thread on dual core CPU, (b) two threads on dual core CPU, (c) four threads
on quad core CPU, (d) a single thread with MDGRAPE-2 and (e) a single thread withMDGRAPE-3. For the two cases using MDGRAPE hardware, the Biot-Savart integral
is accelerated by them. The result is shown in Fig. 5. Detailed specifications of the
hardware used here are shown in Table 1. For 8036 point vortices case, calculation time$n$

8$\overline{\frac{y}{\propto\Phi}}$$\tilde{n}$$\Sigma^{\sim}0\varpi\underline{\in\Phi}Q)$\={o}$\sigma E\Phi$$\frac{\approx oc\varpi}{\frac{8}{0\alpha}}$Number of point vortices$(x10^{4})$Figure 5: Calculation time per one time step is plotted against the number of point
vortices without image vortices.

Table 1: Specifications of CPUs and memories.

9

Figure 6: Density of state is plotted against system energy $E$ and inertia .
$I$

is (a) 6.01 sec., (b) 3.01 sec., (c) 1.67 sec., (d) 0.49 sec., (e) 0.068 sec.. The speed ratio to(a) is (a) 1.00, (b) 1.99, (c) 3.59, (d) 12.2, (e) 88.5. In other words, MDGRAPE-3 finishes
simulation for a day that takes 89 days with PC without MDGRAPE-3.

5

Characterization of point vortex system by temperature

In this section, we demonstrate the simulation results for the statistical understandingof the two-dimensional point vortex system [1].

5.1

Density of state

To check the existence of the negative temperature state in the target point vortex
system, the density of state is obtained by a random sampling of states following the
micro canonical statistics. Each state is characterized by energy $E$ and inertia which
are determined by randomly generated distribution of point vortices. Density of state
consisting ofstates is shown in Figs. 6 and 7. The density of state has a peak at$E=29.1$ and $I=0.O$ , which gives an evidence for the negative temperature state. The
ridge extends on $I=0$ plane that is a symmetric plane of the density of state. As the$I$$10^{8}$

10Energy $E(x10)$

Figure 7: Density of state on $I=0$ plane is shown.
number of vortices increases, the peak becomes steep one and peak position approaches$E=0$ .

5.2

Equilibrium distribution

It may be likely that equilibrium distribution of the vortices changes with the temperature. So, equilibrium distributions of the vortices at various values of temperature
are obtained time-asymptotically by time development simulations. The temperature iscontrolled by the initial distribution of the vortices that determines the system energy.
The results are shown in Fig. 8. In positive temperature, both sign of vortices mix witheach other and spread over the circular area uniformly. On the other hand, when thesign of the temperature changes with energy increase, the same-sign vortices tend to form
small clumps. In larger energy case, the clumps size becomes gradually large and theconfiguration finally reaches a dipole one. Note that the background vortices enable theclump formation that needs large energy in the energy conserving system. The energy
belonging to the background vortices is relatively low compared with the vortices in theclumps.

5.3

Two-body correlation function

Here two-body correlation function is defined by a distribution of the distances for
all combinations of two vortices. The value of the correlation function is normalized by

11

Figure 8: Equilibrium distributions of the vortices are shown. In the right edge distribution, the upper and lower clump exclusively consists of the positive and negative vortices,
respectively.
the total number of the combinations of two vortices. In Fig. 9, the correlation function
for the negative temperature case whose vortex distribution is given in the right edge in
Fig. 8 is shown. In Fig. 10, the correlation function for the positive temperature casewhose distribution is given in the left edge in Fig. 8 is shown. In each Fig., two lines
are indicated. The first one corresponds the distance distribution of the positive vortices
around the positive vortices, and the second one corresponds the distance distributionof the negative vortices around the positive vortices. The distance distribution of thenegative vortices around the positive vortices and the one of the negative vortices aroundthe negative vortices are omitted here, because the distributions of positive and negative
vortices are symmetric in the circular area, and the omitted correlation functions arealmost identical to the two cases indicated in Figs. 9 and 10. For the negative temperaturecase, the left peak corresponds the clump distribution of the same sign, and the rightpeak corresponds the clump distribution of the other sign. It clearly shows the clumping

distribution.
In positive temperature case, two lines overlap with each other perfectly, so that we
can see only a line. It means the distributions of the negative and positive vortices are
uniform in the circular area.

5.4

Energy $(k-)$ spectrum

Energy spectrum of the point vortex system bounded by a circular wall is obtained
analytically by [15, 16],$E(k)$$=$$\frac{1}{4\pi k}\sum_{i}^{N}\Gamma_{i}^{2}+\frac{1}{4\pi k}\sum_{i}^{N}\sum_{j\neq i}^{N}\Gamma_{i}\Gamma_{j}J_{0}(k|r_{i}-r_{j}|)$$- \frac{1}{2\pi k}\sum_{i}^{N}\sum_{j}^{N}\Gamma_{i}\Gamma_{j}\sum_{l=0}^{\infty}\epsilon_{\ell}(\frac{|r_{j}|}{R})^{\ell}\cross J_{\ell}(kR)J_{\ell}(k|r_{i}|)\cos(\ell(\varphi_{i}-\varphi_{j}))$12$A-\dot{\overline{D}}oe-\llcorner\Phi\circ>$

Distance between two vortices (R)Figure 9: The correlation function for the negative temperature case whose vortex distribution is given in the right edge in Fig. 8 is shown..

$\delta-$$\overline{\dot{\overline{A}}}$$o\dot{e}ro\varpi$

Distance between two vortices (R)Figure 10: The correlation function for the positive temperature case whose distribution
is given in the left edge in Fig. 8 is shown.

13

$k$

Figure 11: Energy spectra are plotted at various temperatures.
controlled by the initial configuration of the vortices.

The temperature is$+ \frac{1}{2\pi k}\sum_{i}^{N}\sum_{j}^{N}\Gamma_{i}\Gamma_{j}\sum_{\ell=0}^{\infty}\epsilon_{\ell}(\frac{|r_{j}|}{R})^{p}xJ_{\ell}^{2}(kR)\cos(l(\varphi_{i}-\varphi_{j}))$$\epsilon_{\ell}$$x_{i}$

$=$

$=$

(16)

$\{\begin{array}{l}1 P=0,2 \ell\geq 1,\end{array}$$r_{i}|\cos(\varphi_{i})$

,

(15)

$y_{i}=|r_{i}|\sin(\varphi_{i})$

(17)

.

The first and the second terms obtained by Novikov gives the spectrum for unboundedpoint vortex system, and the rest represents the effect of the circular boundary. By
using the formula (15), we obtain the spectra for various temperatures shown in Fig.
11. As is readily seen, the slope of the intermediate decreases as the energy increases.
Spectrum range is limited by the radius of the boundary in the small region and by theminimum distance between the vortices in the large region. In Fig. 12, the slope of
the intermediate is plotted against the system energy. Up to $E\approx 2x10^{4}$ , the negative
slope decreases with the energy increase. In larger energy, the slope approaches to a
constant value, and the asymptotic value may be $-2.0$ . However, it is possible that theslope cannot develop sharply due to a lack of simulation accuracy, i.e., integral precision is
insufficient. Thus more detailed investigation is needed to conclude the slope approaches$-2.0$ in the high energy region.
In the positive temperature case, the slope should be $-1$ due to the cancellation of
the terms in Eq. (15) by the uniform distribution of the positive and negative vortices.
$k$

$k$

$k$

$k$

14

Energy $(x10^{4})$Figure 12: The slope in the intermediate$k$is plotted against the system energy.

However, there is a case where the slope is not $-1$ as shown in the lowest line in Fig. 11.
Equilibrium distributions of the two cases where the slope is equal to-l and unequal to-l

are shown in Fig. 13. One may understand that there is no remarkable difference between
the two cases. However, the difference can be elucidated by the two body correlation
function. The correlation functions are shown in Figs. 14 and 15 for the distributions in
Fig. 13. The difference appears in the distribution probability of positive vortices aroundpositive vortices. In Fig. 14 where the slope is equal to-l, the value approaches zero in
the limit of zero distance. On the other hand, in Fig. 15, the value does not approachzero in the limit of zero distance. Thus the distribution probability of positive vortices

Figure 13: Equilibrium distributions of the two cases where the slope is equal to-l andunequal to-l are shown.

15$\hat{\uparrow 1\not\in\frac{.g}{\triangleleft^{-}.}\frac{o}{s}}$Figure 14: Two-body correlation function where the slope is equal to-l is plotted.

.

$\uparrow^{\wedge}\frac{g}{\dot{\overline{B}}}n^{n}\frac{o}{\dot{\check}}$$Ae$

Figure 15: Two-body correlation function where the slope in unequal to-l is plotted.

16(a) Slope $=-1$(b)

Slope#-lFigure 16: Two equilibrium distributions are illustrated. The signs
the positive and negative vortices, respectively.

$+$

and

$\cross$

represent

around positive vortices and that of negative vortices around positive vortices does notmatch, and the cancellation of the terms in Eq. (15) is imperfect, so that the slope is not$-1$ . The correlation functions also suggest that there are at least two equilibrium states
as shown in Fig. 16.

6

Discussion

The massive numerical simulation results using the special-purpose supercomputer,
MDGRAPE series have been presented to understand the statistical mechanical featureof the two-dimensional point-vortex system.
Recently, computational power has grown rapidly and PC cluster enables massive numerical simulations in university laboratory. However, simulations that can be accelerated
by distributed computing is restricted, which depends mainly the physical model itself.
Physical problems that need integrals over whole space, including Poisson equation, BiotSavart integral, may be inappropriate for the distributed computing. One of the ways
to solve the above problem is a special-purpose computer. Nowadays, integrated circuitsthat can be programmed by a user are widely used in industry, which is called Field
Programmable Gate Array (FPGA). I believe the FPGA technology breaks the usability
and the performance of the classic vector-type supercomputer.
I did not mention here, but an intriguing aspect of two-dimensional turbulence is avortex crystal that is understood as a metastable state. Several ideas to understand themetastable state have been proposed. However these are applicable for only limited cases
and no fundamental theory is established.

17

Acknowledgement
The author thanks Professor Yasuhito Kiwamoto, Professor Hiroyuki Tomita and Dr.
Mitsusada M. Sano for valuable discussion. The author thanks Dr. Toshikazu Ebisuzaki,
Dr. Makoto Taiji and Dr. Tetsu Narumi for a support of MDGRAPE-2 and MDGRAPE-

3.

References[1] Y. Yatsuyanagi, Y. Kiwamoto, H. Tomita, M. M. Sano, T. Yoshida, and T. Ebisuzaki.
Phys. Rev. Lett., Vol. 94, p. 054502, 2005.
[2] L. Onsager. Nuovo Cimento Suppl., Vol. 6, p. 279, 1949.
[3] G. Joyce and D. Montgomery. J. Plasma Phys., Vol. 10, p. 107, 1973.

[4] Y. B. Pointin and T. S. Lundgren. Phys. Fluids, Vol. 19, p. 1459, 1976.
[5] C. E. Seyler, Jr. Phys. Fluids, Vol. 19, p. 1336, 1976.
[6] O. B\"uhler. Phys. Fluids, Vol. 14, p. 2139, 2002.
[7 D. J. Johnson. Phys. Fluids, Vol. 31, p. 1856, 1988.
$|$[8] D. H. E. Dubin and T. M. O’Neil. Rev. Mod. Phys., Vol. 71, p. 87, 1999.
[9] R. C. Davidson. Physics of nonneutral plasmas. Imperial College, London, 2001.

[10] J. Sommeria. Two-dimensional turbulence. In M. Lesieur, A. Yaglom, and F. David,
editors, New Trends in turbulence. Springer, 2000.
[11] P. Tabeling. Phys. Rep., Vol. 362, p. 1, 2002.

[12] G. L. Eyink and K. R. Sreenivasan. Rev.. Mod. Phys., Vol. 78, p. 87, 2006.
[13] D. Patterson and J. Hennessy. Computer Organization and Design, 3rd ed. Morgan
Kaufmann, San Fransisco, 2004.
[14] T. Ito, K. Makino, T. Ebisuzaki, and D. Sugimoto. A special-purpose n-body machine
grape-l. Comput. Phys. Comm., Vol. 60, p. 187, 1990.

[15] T. Yoshida and M. M. Sano. J. Phys. Soc. Jpn., Vol. 74, p. 587, 2005.

[16] E. A. Novikov. $Sov$ . Phys. JETP, Vol. 41, p. 937, 1975.