The magneto-rotational instability near threshold: spatio-temporal amplitude equation and saturation

0.23 MB·English

by  O. Regev

Title : will be set by the publisher
Editors : will be set by the publisher
EAS Publications Series, Vol. ?, 2008
8
0
THE MAGNETO-ROTATIONAL INSTABILITY NEAR
0
THRESHOLD: SPATIO-TEMPORAL AMPLITUDE EQUATION
2
AND SATURATION
n
a
J Oded Regev1,2

0
1
Abstract. Weshow,bymeansofaperturbativeweaklynonlinearanal-
]
ysis, that the axisymmetric magneto-rotational instability (MRI) in
h
p a magnetic Taylor-Couette (mTC) flow in a thin-gap gives rise, for
- verysmallmagneticPrandtlnumbers(Pm),toarealGinzburg-Landau
o
equation for the disturbance amplitude. The saturation amplitude A
s
tr is found to scale in this regime as Pmδ, with 1/2 < δ < 2/3 (depend-
s ing on the boundary conditions adopted). The asymptotic results are
a
[ shown to comply with numerical calculations performed by using a
spectral code. They suggest that thetransport due to the nonlinearly

1 developed MRI may bevanishingly small for Pm ≪1.
v
7
3
1 Introduction
6
1
. 1.1 General
1
0 Asymptotic approaches to nonlinear stellar pulsation were pioneered by Robert
8
Buchler and his collaborators (see this volume). Amplitude equations, which
0
capture the essential dynamics near the instability threshold, result from such
:
v
approaches and may usually be derived using singular perturbation theory (see
i
X Marie-Jo Goupil, this volume). In this contribution we report on the results of
r applying a perturbative asymptotic analysis to an instability that has acquired,
a in the past 15 years or so, a paramount importance in the quest to understand
angular momentum transport in accretion disks, namely the magneto-rotational
instability (MRI). Even though the system studied here is rather different from
pulsating stars, the ideas and techniques are quite similar. One important differ-
enceisthatwhiletheamplitudeequationsforstellarpulsationaresetsofordinary
differentialequations,herewemusttakerecoursetopartialdifferentialequation(s)
1 DepartmentofPhysics,Technion-IsraelInstitute ofTechnology
e-mail: [email protected]
2 DepartmentofAstronomy,ColumbiaUniversity
(cid:13)c EDPSciences 2008
DOI:(willbeinsertedlater)
2 Title : will be set by the publisher
(PDE), in which the amplitude (usually called in this case envelope) is a function
of time (as in the ordinary amplitude equations) and space. The application of
singularperturbationtheoryrequiresthustheuseofslowvariationintime aswell
asinspace. Inthiswayspatio-temporalslowdynamics (i.e.,patterns)iscaptured
by a generic PDE - the real Ginzburg-Landau Equation, rGLE, in our case.
The rudiments of singular perturbation theory are well explained in the book
by Bender & Orszag (1999), while a detailed account on problems dealing with
the dynamics on the center manifold can be found in Manneville (1990). It seems
that these powerful analytical and semi-analytical techniques have not yet found
their way to enough astrophysical applications. In the early 1980s Robert Buch-
ler started an ambitious program of introducing such a new approach to stellar
pulsation. I was fortunate to be his postdoc then, but have moved on to other
topics. Robert and his postdocs (most of which contributed to this volume) have
pursued this program and during this conference one could learn just how much
has been achieved in the understanding of stellar pulsation, following asymptotic
approaches, combined with numerical and phenomenological ones. It is quite un-
fortunate that ”driving forces” of Robert’s magnitude have unfortunately been
quite rare in other branches of astrophysics.
1.2 The MRI - background
ThelinearMRIhasbeenknownforalmost50years(Velikhov1959,Chandrasekhar
1960): RayleighstablerotatingCouette flowsofconductingfluidsaredestabilized
in the presence of a vertical magnetic field if dΩ2/dr < 0 (angular velocity de-
creasing outward). However, the MRI only acquired importance to astrophysics
after the influential work of Balbus & Hawley (1991), who demonstrated its vi-
ability to cylindrical accretion disks (for any sufficiently weak field). The linear
analysis was later supplemented by nonlinear numerical simulations (albeit of a
small, shearing-box (SB), segment of the disk). Enhanced transport (conceivably
turbulent) is necessary for accretion to proceed in disks, found in a variety of as-
trophysical settings - from protostars to active galactic nuclei. As is well known,
Keplerian rotation law is hydrodynamically linearly stable and thus the MRI has
been widely accepted as an attractivesolution for enhanced transport(see the re-
viewsbyBalbus&Hawley1998andBalbus2003),eventhoughsomequestionson
the nature ofMRI-driven turbulence still remain(see, e.g., Branderburg2005), in
particularin view ofsome recentnumericalstudies (Pessah,Chan& Psaltis 2007,
Fromang & Papaloizou 2007, Fromang et al. 2007). Because of the MRI’s im-
portance and some of its outstanding unresolvedissues, severalgroups have quite
recently undertaken projects to investigate the instability under laboratory con-
ditions. Additionally, numerical simulations specially designed for experimental
setups have been conducted. It appears, however, quite impossible to make defi-
nite deductions from these experiments and simulations (see, e.g., Ji et al. 2001,
Noguchiet al. 2002, Sisan et al. 2004, Stefani et al. 2006, andreferencestherein).
We shall report here on the result of a weakly nonlinear analysis of the MRI
nearthreshold,foradissipativemTCflow. Wehavedonetheanalysisfortheabove
Regev: MRI - amplitude equation & saturation 3
simplified(relativelytoanaccretiondisk)problemfortwotypesofboundarycon-
ditions(seeUmurhan,Menou&Regev2007,hereafterUMRandUmurhan,Regev
& Menou 2007, hereafter URM, for a detailed account). This kind of approach
is important because the viability of the MRI as the driver of turbulence and
angular-momentum transport relies on understanding its nonlinear development
andsaturation. Bycomplementingnumericalsimulationswithanalyticalmethods
useful physical insight can be expected. To facilitate an analytical approach we
made a number of simplifying assumptions so as to make the system amenable
to well-known asymptotic methods (see, e.g., Cross and Hohenberg 2003, Regev
2006) for the derivation of a nonlinear envelope equation, valid near the linear
instability threshold.
2 Linear theory
2.1 Equations
Thehydromagneticequationsincylindricalcoordinates(Chandrasekhar1961)are
applied to the neighborhood of a representative radial point (r0) in a mTC setup
with an imposed constant background vertical magnetic field. The steady base
flow has only a constant vertical magnetic field, B = B0ˆz, and a velocity of
the form V = U(x)yˆ. In this base state the velocity has a linear shear profile
U(x)= qΩ0x,representinganazimuthalflowaboutapointr0,thatrotateswith
arateΩ0−,definedfromthedifferentialrotationlawΩ(r) Ω0(r/r0)−q (q =3/2for
∝
Keplerianrotation). The totalpressure in the base state (divided by the constant
density), Π≡ρ0−1 P + B8π02 , is a constant and thus its gradient is zero.
This base flow i(cid:16)s disturb(cid:17)ed by 3-D perturbations on the magnetic field b =
(b ,b ,b ),aswellasonthe velocity-u=(u ,u ,u ), andonthetotalpressure-
x y z x y z
̟. We consideronlyaxisymmetricdisturbances,i.e. perturbationswithstructure
only in the x and z directions. This results, after non-dimensionalization, in the
following set of equations, given here in the rotating frame:
du 1
2Ω0ˆz u qΩ0uxyˆ b b B0∂zb = ̟+ 2u, (2.1)
dt − × − −C ·∇ −C −∇ ∇
R
db 1
b u+qΩ0bxyˆ B0∂zu = 2b, (2.2)
dt − ·∇ − ∇
m
R
u ∂ u +∂ u =0, b ∂ b +∂ b =0. (2.3)
x x z z x x z z
∇· ≡ ∇· ≡
The Cartesiancoordinatesx,y,z correspondto the radial(shear-wise),azimuthal
(stream-wise) and vertical directions, respectively, and since axisymmetry is as-
sumed xˆ∂ +ˆz∂ and 2 ∂2+∂2. Lengths have been non-dimensionalized
byL( ∇th≡egapxsize)z,time∇byt≡helxocalrzotationrateΩ˜0 (tildesdenotedimensional
≈
quantities). Because the dimensional rotation rate of the box (about the central
object)isΩ˜0 =Ω˜0ˆz,thenon-dimensionalquantityΩ0 issimplyequalto1,butwe
keepittoflagtheCoriolisterms. VelocitieshavebeenscaledbyΩ˜0Landthemag-
netic fieldbyitsbackgroundvalueB˜0. ThusB0 1aswell,butagain,weleaveit
≡
4 Title : will be set by the publisher
in the equation set for later convenience. The hydrodynamicpressure is scaledby
ρ˜0L2Ω˜20 and the magnetic one by B˜02/(8π). The non-dimensional perturbation ̟
ofthetotalpressuredividedbythedensity(whichisequalto1innon-dimensional
units), which survives the spatial derivatives, is thus given by ̟ = p+ 1 b2,
C2| |
where p is the hydrodynamicpressureperturbation. The non-dimensionalparam-
eter C ≡B˜02/(4πρ˜0Ω˜20L2)=V˜A2/V˜2 is the Cowling number, measuring the relative
importance of the magnetic pressure to the hydrodynamical one. It is equal to
the inverse square of the typical Alfv´en number (V˜ is the typical Alfv´en speed).
A
The Cowling number appears in the non-linear equations, together with the two
Reynolds numbers: Ω˜0L2/ν˜ and m Ω˜0L2/η˜, where ν˜ and η˜ are, respec-
R ≡ R ≡
tively,themicroscopicviscosityandmagneticresistivityofthefluid. Weshallalso
see that the magnetic Prandtl number, given as / , plays an important
m m
P ≡R R
role in the nonlinear evolution of this system.
Itisusefultorewritetheaboveequationsofmotionintermsofmoreconvenient
dependent variables. Because the flow is incompressible and y-independent, the
radial and vertical velocities can be expressed in terms of the streamfunction, Ψ,
that is, (u ,u ) = (∂ Ψ, ∂ Ψ). Also, since the magnetic field is source free,
x z z x
−
one can similarly express its vertical and radial components in terms of the flux
function, Φ, that is, (b ,b )=(∂ Φ, ∂ Φ) (see URM for details). The system is
x z z x
−
supplementbyappropriateboundaryconditiononthethincylindricalgap”walls”,
as detailed in UMR and URM and we also point out that the stress relevant for
angular momentum transport, Σ, is
Σ=ΣR+ΣM, ΣR uxuy, ΣM bxby.
≡ ≡−C
whichconsistsofΣR andΣM,theReynolds(hydrodynamic)andMaxwellstresses,
expressing the velocity and magnetic field disturbance correlations, respectively.
2.2 Linear stability
Linearizationof (2.1-2.3), for perturbations of the form est+ikxx+ikzz, gives rise
∝
to the dispersion relation,
(s;kx,kz, m, , ,q)=a0s4+a1s3+a2s2+a3s+a4 =0,
D P S C
where all the a are functions of k , k and the other parameters, but we shall
i x z
write explicitly only the coefficient that will be used in what follows,
a4 = C4 kT2C(CkT4Pm+kz2S2)2+κ2S2CkT4kz2−2qS4kz4 , (2.4)
S h i
where the notation κ2 2(2 q), k2 k2+k2 is introduced.
≡ − T ≡ x z
For given values of the parameters (call them p ) there will be four distinct
j
modes. The linear theory in various limits for this problem has been discussed in
numerous publications and we shall not elaborate upon it any further. Rather,
we focus on situations where the most unstable mode (of the four) is marginal
(critical) for some k at a given value of k = K. We fix K in order to focus on
z x
Regev: MRI - amplitude equation & saturation 5
(a) (b)
R = 5.3 5.4
0.02 m m
R5.3 Unstable
e(s) 0 s #: 5.2
R d
e: −0.02 nol5.1
owth rat−0.04 Rm = 4.9 etic Rey4.59 Rm(crit)
gr−0.06 Rm = 4.6 Magn4.8
k=Q
−0.08 4.7
0.5 1 1.5 0 0.5 1 1.5
vertical wavenumber: k vertical wavenumber: k
Fig. 1. Summary of linear theory. This example is for C =0.08, P =0.001, q =3/2,
m
and thefundamentalmode. (a) Growth rates, Re(s),as afunction ofwavenumber k for
three values of R . (b) Solid line depicts those values of R and k where Re(s) = 0.
m m
The shaded region shows unstable modes. The locations of k = k ≡ Q and R =
crit m
R (crit)≡R are shown.
m m
the marginal vertical dynamics within the thin gap (actually a rotating channel,
see below). Marginality(s=0) implies the vanishing ofthe realcoefficient a4(kz)
(as expressed above) and its derivative with respect to k at some k =Q:
z z
∂a4
a4(kz =Q;K,pj)=0, (kz =Q;K,pj)=0. (2.5)
∂k
z
The second condition and (2.4) yield
Q Q2 4( Q2 4q)+(K2+Q2)2Q2 2( 2+κ2)+
S C − CS S
(Kh2+Q2)2 2(6 mQ2+κ2)+
CS CP
2(K2+Q2)3 m 2 2+5(K2+Q2)4 m2 3 =0. (2.6)
P C S P C
i
This equation, together with a4(kz = Q,K;pi) = 0, can be solved for and Q.
S
Thegeneralexpressionsfor (K, m, ,q)andQ(K, m, ,q)arelengthybuttheir
S P C P C
asymptotic forms, to ( m) (for m 1), are simple,
O P P ≪
16 q(2 q)K 2q K2
= C − , Q2 =K2 −C (2.7)
S (2q K2) 2q+ K2
p −C C
If K2 > 2q the solutions are not physically meaningful, while K2 = 2q corre-
C C
sponds to the ideal MRI limit.
Our system is a rotating thin channel, whose walls are at x = 0,π/K. All
quantitiesareverticallyperiodiconascaleL commensuratewithintegermultiples
z
6 Title : will be set by the publisher
of 2π/Q. As long as L 1/Q, the limit of a vertically extended system (and
z
≫
thus a continuous spectrum of vertical modes) may be effected. Linear theory, as
discussed above, is summarized in Figure 1
3 Weakly nonlinear analysis
We follow the fluid into instability by tuning the vertical background field away
from the steady state, i.e. B0 1 ǫ2λ˜, where ǫ 1 and λ˜ is an (1) control
→ − ≪ O
parameter. This means that we are now in a position to apply procedures of
singular perturbation theory to this problem, by employing the multiple-scale (in
z and t) method (e.g, Bender & Orszag1999). It facilitates, by imposing suitable
solvabilityconditionsateachexpansionorderofthecalculationinordertoprevent
abreakdowninthesolutions,aderivationofanenvelopeequationfortheunstable
mode. Thus for any dependent fluid quantity F(x,z,t) we assume
F(x,z,t)=ǫF1(x,z,t)+ǫ2F2(x,z,t)+ . (3.1)
···
The fact that the x and z components of the velocity and magnetic field pertur-
bations can be derived from a streamfunction, Ψ(x,z) and magnetic flux function
Φ(x,z) reduces the number of relevant dependent variables F to four (u and b
y y
are the additional two). These four dependent variables will also be used in the
spectral numerical calculation (see below).
For the lowestǫ order of the equations, resulting from substituting the expan-
sions into the originalPDE and collecting same order terms, we make the Ansatz
F1(x,z,t) = Fˆ1A˜(ǫz,ǫ2t)eiQzsinKx+c.c., where Fˆ1 is a constant (according to
the variable in question), and where the envelope function A˜ (an arbitrary con-
stant amplitude in linear theory) is now allowed to have weak space (on scale
Z ǫz) and time (on scale T ǫ2t) dependencies. Because this system is tenth
≡ ≡
order in x-derivatives, a sufficient number of conditions must be specified at the
edges. andthe Ansatz hasto obeythem. InUMRwe chosemathematicallyexpe-
dient boundary conditions, which allowed fully analytical treatment for the limit
m 1 andhere we shallreport,in some detail, only onthe resultsofthat work.
P ≪
Wehavealsotriedsomedifferentboundaryconditions(seeURM)andinthatcase
the rGLE coefficients had to be calculated numerically. The qualitative behavior
of the saturation amplitude was found to be quite similar (see below).
3.1 Real Ginzubrg-Landau equation
The end result of the asymptotic procedure procedure is the well-known real
Ginzburg-Landau Equation (rGLE) which, for m 1, is
P ≪
1
∂ A=λA AA2+D∂2A. (3.2)
T − m | | Z
P C
We shall consider only the magnitude (real part) of A. In a 1D gradient system
like this one phase dynamics can only give rise to wavelength modulations and
Regev: MRI - amplitude equation & saturation 7
consequently our main results are in no way affected by it. The real envelope is
A √ξA˜, λ ζλ˜, and the coefficients, for q =3/2 are
≡ ≡
3 5 4 18 2 32+2( 2+16)√ 2+1
ξ = S − S − S S , (3.3)
4 · ( 2+1)(4√ 2+1 3)
S S S −
2+2 2√ 2+1 ( 2+1)
D = 6 S − S S , (3.4)
3(4√ 2+1 3)
(cid:0) (cid:1)
S S −
3 √ 2+1 6
ζ = S S − S , (3.5)
4 2+1+√ 2+1
S S
where (K) is asgivenin(2.7). For 1,(i.e. as oneapproachesthe idealMRI
S S ≫
limit), these simplify toξ =15/16, ζ =3/4, D =3/2;andingeneraltheyremain
(1) quantities for all reasonable values of .
O S
3.2 Scaling of angular momentum transport
Contributionsto the totalangularmomentum transport(J˙ =J˙H+J˙B) due to the
hydrodynamic (J˙H) and magnetic correlations (J˙B) are, to leading order,
9ǫ2 2+ 2 2√1+ 2
J˙H = S 1−+ 2 S !A2+O ǫ3,ǫ2Pm ,
S S
(cid:0) (cid:1)
3ǫ2 1
J˙B = 1 A2+ ǫ3,ǫ2 m . (3.6)
− √1+ 2 O P
S (cid:18) S (cid:19)
(cid:0) (cid:1)
The envelope A is found by solving the rGLE, which is a well-studied system
(see, e.g., Regev 2006, for a summary and references). It has 3 steady uniform
solutions in 1D (here, the vertical): an unstable state A = 0, and 2 stable states
A=±As, whereA2s =|ζ(S)|PmC (note thatζ(S) is anO(1)quantity). As is also
the saturation amplitude, because the system will develop towards it.
SettingA As in(3.6)andintheexpressionforEV (thetotaldisturbanceen-
→
ergy,seeUMR),followedbysomealgebra,revealsthattheangularmomentumflux
in the saturated state is, to leading order in m and ǫ, J˙ =ǫ2 ζ( )γ( ) m −1,
P | S | S P CS
while in the expression for the energy one term is independent of m, EV =
ǫ4 3ζ2( )β2( ) + ( m). The key results of this analysis (see URMP) are the
C S S O P
following scalings (with the magnetic Prandtl number)
As √ m EV E0, J˙ m (3.7)
∼ P −→ ∼ ∼P
for m 1, where E0 is a constant, independent of m 1. For fixed resistivity
thisPimp≪lies J˙ −1. P
∼R
1FortheboundaryconditionsusedinURMwenumericallygotAs∝Pm2/3
8 Title : will be set by the publisher
(a) (b)
100 10−2
y10−1 sport10−3 R = 5x103 R = 104
g n
al ener10−2 R = 104 R = 5x104 um tra10−4 R = 5x104
: totV10−3 ment
E10−4 R = 5x103 R = 2x105 . J: mo10−5
10−5 10−6 R = 2x105
100 101 102 103 100 102
Time (in box orbits) Time (in box orbits)
Fig. 2. EV (panel a) and J˙ (panel b) as a function of time from numerical calculation.
The different lines are labeled by the value of the Reynolds number R The diamonds
in panel (b) show the scaling predicted by our asymptotic analysis, predicting also a
constantfinalvalueofthedisturbanceenergy,asapparentinpanel(a). One”boxorbit”
=2π in nondimensional units.
We havealsoperformedfullynumericalcalculations,usinga2-Dspectralcode
to solve the originalnonlinear equations, in the streamfunctionand magnetic flux
function formulation, near MRI threshold. The code implemented a Fourier-
Galerkin expansion in 64 64 modes in each of the four independent physical
variables,i.e. F= n,mFn×m(t)sinKnxeiQmz+c.c.,whereF=(ψ,uy,Φ,by)T and
where Fnm(t) is the time-dependent amplitude of the particular Fourier-Galerkin
P
mode in question (denoted by the indices n,m). We typically started with white
noise initial conditions on Fnm(0) at a level of 0.1 the energy of the background
shear. Because of space limitations we display only a representative plot of runs
made with = 5.0, = 0.05, (i.e. m fixed), q = 3/2 and ǫ2 = 0.2 for a few
S C R
successively increasing values of (Figure 2). With this value of ǫ the fastest
R
growing linear mode has a growth rate of 0.065, and thus fully developed ideal
∼
MRIcannotbe expected. EV saturatesataconstantindependentof (forsmall
enough m), while J˙ saturates at values that scale as −1. Thus theRasymptotic
P R
analysis fits very well the fully numerical results after the system saturates.
4 Summary
The numerical and asymptotic solutions developed here show that for m 1,
it is the azimuthal velocity perturbation, arising from the ǫ2 term,Ptha≪t be-
O
comes dominant in the saturated state. It appears to be the primary agent in
(cid:0) (cid:1)
the nonlinearsaturationofthe MRI inthe channel,actinganisotropicallyso asto
modify the shear profile and results in a non-diagonal stress component (relevant
Regev: MRI - amplitude equation & saturation 9
for angular momentum transport). Our analysis is complementary to the study
of Knobloch and Julien (2005), who have performed an asymptotic MRI analysis,
but for a developed state, far from marginality.
The trends predicted by this simplified model are not qualitatively altered by
different boundary conditions, the value of δ changing only by a small amount.
Qualitatively similar scalings, but with somewhat different values of δ > 0, have
very recently been reported by Lesur & Longaretti (2007) and Fromang et al.
(2007),howeverinthoseSBsimulations m wasnottakentobevanishinglysmall.
P
It should be remarkedthat conventionalaccretiondisks do have m 1, but the
P ≪
numericalSBsimulationscannotfaithfullytreatverysmallsuchPrandtlnumbers,
for numerical reasons. Further analytical investigations of the kind reported here
should contribute to better physical understanding of the MRI saturation..
References
Balbus, S.A.2003, ARA&A,41, 555
Balbus, S.A.and Hawley, J.F. 1991, ApJ,376, 214
Balbus, S.A.and Hawley, J.F. 1998, Rev.Mod. Phys., 70, 1
Bender,C.M.andOrszag,S.A.1999,Advanced Mathematical Methods forScientists and
Engineers, Springer, New York
Brandenburg,A. 2005, Astron.Nachr., 326, 787.
Chandrasekhar, S. 1960, Proc. Natl. Acad. Sci. USA,46, 253.
Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability, Oxford University
Press, Oxford
Cross, M.C. and Hohenberg,P.C. 2003, Rev.Mod. Phys. 65, 851
Fromang, S. and Papaloizou, J. 2007, A&A,476, 1113
Fromang, S., Papaloizou, J., Lesur, G. and Heinemann, T. 2007, A&A,476, 1123
Ji, H., Goodman, J. and Kageyama, A.2001, MNRAS,325, L1
Knobloch, E. & Julien, K. 2005, Phys. Fluids, 17, 094106
Lesur, G., & Longaretti, P.-Y.2007, MNRAS,378, 1471
Manneville, P. 1990, Dissipative Structures and Weak Turbulence, Academic Press,
Boston
Noguchi, K., Pariev, V.A., Colgate, S.A., Beckley, H.F. and Nordhaus, J. 2002, ApJ,
575, 1151
Pessah, M.E., Chan, C-k. & Psaltis, D.2007, ApJ, 668, L51
Regev, O. 2006, Chaos and Complexity in Astrophysics, Cambridge University Press,
Cambridge
Sisan, D.R., Mujica, N., Tillotson, W.A., Huang, Y., Dorland, W., Hassam, A.B., An-
tonsen, T.M. and Lathrop, D.P., 2004, Phys. Rev.Lett., 93, 114502
Stefani,F.,Gundrum,T.,Gerbeth,G.,Rdiger,G.,Schultz,M.,Szklarski,J.andHoller-
bach, R., 2006, Phys. Rev.Lett., 97, 184502
Umurhan,O.M. and Regev, O. 2004, A&A,427, 855.
Umurhan,O.M., Menou, K.& Regev,O. 2007, Phys. Rev.Lett., 98, 034501
Umurhan,O.M., Regev, O. & Menou, K. 2007, Phys. Rev.E, 76, 036310
Velikhov,E.P. 1959, Sov.Phys. JETP, 9, 959.