format: text

text:
# Intermediate phase in the spiral antiferromagnet Ba2CuGe2O7

## I Introduction

A recent experimental investigation  of the magnetic properties of Ba2CuGe2O7 in its low-temperature phase ($`T<T_N=3.2`$ K) established the occurrence of spiral antiferromagnetic order due to a Dzyaloshinskii-Moriya (DM) anisotropy . A schematic illustration of the spiral abstracted from experiment may be found in Fig. 5 of Ref. 1. It was further demonstrated that a Dzyaloshinskii-type  commensurate-incommensurate (CI) phase transition is induced by a magnetic field $`H`$ applied along the $`c`$-axis. As the field approaches a critical value $`H_c2`$ T, the spiral is highly distorted while its period (pitch) grows to infinity. For $`H>H_c`$ the ground-state configuration is thought to degenerate into a uniform spin-flop state. This phase transition is similar to the cholesteric-nematic transition induced by an external magnetic field in liquid crystals .

It is of obvious interest to describe theoretically the magnon excitations measured by inelastic neutron scattering , but progress has been hindered by the great formal complexity of the calculation. Here we explore a new approach in which the original discrete system is replaced by a continuum field theory. We are thus able for the first time to carry out a complete calculation of the low-energy excitation spectrum for any strength of the applied field and any direction of spin-wave propagation. In addition, our analysis reveals the existence of a new intermediate phase whose properties we examine and compare with experiment.

In Sec. II the low-energy dynamics is described in terms of a nonlinear $`\sigma `$ model that is compatible with symmetry. In Sec. III we present a brief demonstration of the conventional CI transition which will provide the basis for all subsequent work. The complete field theory is first applied in Sec. IV for an analytical calculation of the field dependence of the magnon spectrum in the high-field commensurate phase. Interestingly, the uniform spin-flop state is shown to be locally stable only for $`H>H_2>H_c`$ where the new critical field $`H_2`$ is predicted to be equal to $`2.9`$ T. A first contact with the measured spectrum is also made in Sec. IV.

The main thrust of our calculation is presented in Sec. V where the determination of the magnon spectrum in the low-field spiral phase is reduced to a quasi-one-dimensional band structure problem that is solved numerically. While an earlier calculation  of the spectrum at $`H=0`$ is confirmed, we are also in a position to analyze existing experimental data at nonzero field and to predict the results of possible future experiments. A byproduct of this analysis is yet another critical field $`H_1=1.7`$ T $`<H_c`$ beyond which the flat spiral ceases to be locally stable. Therefore, the combined results of Secs. IV and V suggest the existence of an intermediate phase in the field region $`H_1<H<H_2`$ whose nature is studied in Sec. VI where we show that a nonflat spiral becomes energetically favorable. The main results are summarized in the concluding Sec. VII, while discussion of some technical issues is relegated to two Appendices.

## II Low-energy dynamics

The unit cell of Ba2CuGe2O7 is partially illustrated in Fig. 1 where we display only the magnetic Cu sites. The lattice constants are $`a=b=8.466`$ Å and $`c=5.445`$ Å. Since the Cu atoms form a perfect square lattice within each plane, with lattice constant $`d=a/\sqrt{2}6`$ Å, it is also useful to consider the orthogonal axes $`x,y`$ and $`z`$ obtained from the original crystal axes $`a,b`$ and $`c`$ by a $`45^{}`$ azimuthal rotation. The complete magnetic lattice is formally divided into two sublattices labeled by $`A`$ and $`B`$ because the major spin interaction between nearest in-plane neighbors is antiferromagnetic. In contrast, the interaction between out-of-plane neighbors is ferromagnetic and weak . Therefore, the interlayer coupling is not crucial for our purposes and is thus ignored in the following discussion which concentrates on the 2D spin dynamics within each layer.

The space group of this crystal is $`D_{2d}^3`$ or $`P\overline{4}2_1m`$ and imposes significant restrictions on the possible types of spin interactions. Such symmetry constraints underlie most of the earlier work  but were not spelled out in sufficient detail. We have thus found it necessary to carry out afresh a complete symmetry analysis, including both nearest-neighbor (nn) and next-nearest-neighbor (nnn) couplings. For the moment, we restrict attention to nn interactions and write the 2D spin Hamiltonian as the sum of four terms :

$$W=W_E+W_{DM}+W_A+W_Z.$$
(1)
Here

$$W_E=\underset{<kl>}{}J_{kl}(𝐒_k𝐒_l)$$
(2)
describes the isotropic exchange over nn in-plane bonds, denoted by $`<kl>`$ , with $`J_{kl}=J`$ for all such bonds. Similarly,

$$W_{DM}=\underset{<kl>}{}𝐃_{kl}(𝐒_k\times 𝐒_l)$$
(3)
stands for antisymmetric DM anisotropy where the vectors $`𝐃_{kl}`$ assume four distinct values:

$`𝐃_I`$ $`=`$ $`D𝐞_2+D^{}𝐞_3,𝐃_{II}=D𝐞_1+D^{}𝐞_3,`$ (4)
$`𝐃_{III}`$ $`=`$ $`D𝐞_2D^{}𝐞_3,𝐃_{IV}=D𝐞_1D^{}𝐞_3,`$ (5)
which are distributed over the 2D lattice as shown in Fig. 2 where nn bonds are accordingly labeled by I, II, III or IV.

Here $`D`$ and $`D^{}`$ are two independent scalar constants, while $`𝐞_1,𝐞_2`$ and $`𝐞_3`$ are unit vectors along the $`x,y`$ and $`z`$ axes of Fig. 1. It should be noted that the $`z`$-components of the DM vectors alternate in sign on opposite bonds, a feature that could lead to weak ferromagnetism. No such alternation occurs for the in-plane components of the DM vectors (4) which are responsible for the observed spiral magnetic order or helimagnetism.

The third term in Eq. (1) contains all “symmetric” anisotropies. Since single-ion anisotropy is not possible in this spin $`s=\frac{1}{2}`$ system, the most general form of $`W_A`$ is

$$W_A=\frac{1}{2}\underset{<kl>}{}\underset{i,j}{}G_{kl}^{ij}\left(S_k^iS_l^j+S_k^jS_l^i\right),$$
(6)
where the indices $`i`$ and $`j`$ are summed over three values corresponding to the Cartesian components of the spin vectors along the axes $`x,y`$ and $`z`$. Accordingly, $`G_{kl}=(G_{kl}^{ij})`$ are $`3\times 3`$ symmetric matrices, one for each bond $`<kl>`$. Again, there exist four distinct such matrices:

$`G_I=\left(\begin{array}{ccc}K_1& 0& 0\\ 0& K_2& K_4\\ 0& K_4& K_3\end{array}\right),`$ $`G_{II}=\left(\begin{array}{ccc}K_2& 0& K_4\\ 0& K_1& 0\\ K_4& 0& K_3\end{array}\right),`$ (13)
$`G_{III}=\left(\begin{array}{ccc}K_1& 0& 0\\ 0& K_2& K_4\\ 0& K_4& K_3\end{array}\right),`$ $`G_{IV}=\left(\begin{array}{ccc}K_2& 0& K_4\\ 0& K_1& 0\\ K_4& 0& K_3\end{array}\right),`$ (20)
which are all expressed in terms of the four scalar parameters $`K_1,K_2,K_3`$ and $`K_4`$. The latter may be further restricted by the trace condition $`K_1+K_2+K_3=0`$ because the isotropic component of the exchange interaction is already accounted for by Eq. (2). Finally,

$$W_Z=\underset{l}{}(g\mu _B𝐇𝐒_l)$$
(22)
describes the usual Zeeman interaction with an external field $`𝐇`$.

The discrete Hamiltonian could be employed to analyze this system by standard spin-wave techniques, but the calculational burden is rather significant and has so far prevented a complete determination of the magnon spectrum . Nevertheless, the relevant low-energy dynamics can be efficiently calculated in terms of a continuum field theory which provides a reasonable approximation for Ba2CuGe2O7 because the period of the observed spiral is equal to about 37 lattice constants along the $`x`$-direction. A similar approach is often invoked in the related subject of weak ferromagnetism  and can be implemented by a straightforward step-by-step procedure starting from the original discrete Hamiltonian .

The first step is to group spins into dimers as shown in Fig. 2. Each dimer contains a pair of spins denoted by $`𝐀`$ and $`𝐁`$ and labeled by a common set of sublattice indices $`\alpha `$ and $`\beta `$ that advance along the crystal axes $`a`$ and $`b`$. A more convenient set of variables is given by the “magnetization” $`𝐦`$ and the “staggered magnetization” $`𝐧`$ which are defined as

$$𝐦=\frac{1}{2s}(𝐀+𝐁),𝐧=\frac{1}{2s}(𝐀𝐁),$$
(23)
and satisfy the classical constraints $`𝐦𝐧=0`$ and $`𝐦^2+𝐧^2=1`$. We also introduce space-time variables according to

$$\eta =\sqrt{2}\epsilon \alpha ,\xi =\sqrt{2}\epsilon \beta ,\tau =2s\sqrt{2}\epsilon Jt,$$
(24)
where $`\epsilon `$ is a dimensionless scale whose significance will become apparent as the discussion progresses. The final result will be stated in terms of the coordinates

$$x=\frac{\xi +\eta }{\sqrt{2}},y=\frac{\xi \eta }{\sqrt{2}},$$
(25)
along the $`x`$ and $`y`$ axes of Fig. 1. One should keep in mind that actual distances are given by $`xd/\epsilon `$ and $`yd/\epsilon `$ where $`d=a/\sqrt{2}`$ is the lattice constant of the square lattice formed by the Cu atoms. Finally, we introduce rescaled anisotropy constants and magnetic field as

$`\lambda ={\displaystyle \frac{D}{\epsilon J}},\lambda ^{}={\displaystyle \frac{\sqrt{2}D^{}}{\epsilon J}},𝐡={\displaystyle \frac{g\mu _B𝐇}{2s\sqrt{2}\epsilon J}}`$ , (26)
$`\kappa _0={\displaystyle \frac{2}{\epsilon ^2J}}\left(K_1+K_22K_3\right),`$ (28)
where we display only those combinations of constants that survive in the effective low-energy dynamics. In particular, the constant $`K_4`$ does not appear to leading order. The further notational abbreviations

$$\kappa =\kappa _0\lambda ^2+\lambda _{}^{}{}_{}{}^{2},𝐝_z=\lambda ^{}𝐞_3$$
(29)
will prove convenient in all subsequent calculations.

Now, a consistent low-energy expansion is obtained by treating $`𝐦`$ as a quantity of order $`\epsilon `$ while $`𝐧`$ is of order unity. To leading order, the classical constraints reduce to

$$𝐦𝐧=0,𝐧^2=1,$$
(30)
$`𝐦`$ is expressed entirely in terms of $`𝐧`$ by

$$𝐦=\frac{\epsilon }{2\sqrt{2}}\left[𝐧\times \left(\dot{𝐧}+𝐝_z𝐧\times 𝐡\right)\right]\frac{\epsilon }{2}_1𝐧,$$
(31)
and the $`T=0`$ dynamics of the staggered magnetization $`𝐧`$ is governed by the Lagrangian density $`=_0V`$ where

$`_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{𝐧}^2+𝐡(𝐧\times \dot{𝐧}),`$ (32)
$`V`$ $`=`$ $`{\displaystyle \frac{1}{2}}(_1𝐧\lambda 𝐞_2\times 𝐧)^2+{\displaystyle \frac{1}{2}}(_2𝐧\lambda 𝐞_1\times 𝐧)^2`$ (34)
$`+{\displaystyle \frac{1}{2}}\kappa n_3^2+{\displaystyle \frac{1}{2}}(𝐧𝐡)^2+(𝐡\times 𝐝_z)𝐧.`$
The dot denotes differentiation with respect to the time variable $`\tau `$, $`_1`$ and $`_2`$ are partial derivatives with respect to $`x`$ and $`y`$, and ($`n_1,n_2,n_3`$) are the Cartesian components of $`𝐧`$ along the axes $`xyz`$ of Fig. 1. Consistency requires that all physical predictions derived from Eqs. (31) and (32) must be independent of the specific choice of the scale parameter $`\epsilon `$. This fact will be explicitly demonstrated or used to advantage in the continuation of the paper.

We have further examined possible modifications of the low-energy dynamics due to nnn spin interactions along the diagonals of the Cu plaquettes. Our symmetry analysis revealed that both antisymmetric (DM) and symmetric anisotropies are present over nnn bonds and introduce a new set of parameters. Nevertheless, in the continuum limit, all new parameters merge with those already present in the Lagrangian (32). The implied remarkable rigidity of the effective low-energy spin dynamics is obviously due to the special crystal structure of Ba2CuGe2O7.

In the remainder of this section we make contact with the static energy functional derived by Zheludev et al. , restricted to $`T=0`$, which appears to differ in some respects from the potential $`V`$ of Eq. (32). First, we note that we have omitted from the potential some additive field-dependent constants which play no role except to relate the energy to the magnetization. The latter will be obtained in Sec. III by a direct application of Eq. (31). A more interesting point concerns the special choice of exchange anisotropy made in Ref. 5, which was suggested by the work of Kaplan , Shekhtman, Entin-Wohlman and Aharony , and is referred to as the KSEA anisotropy. If the original perturbative derivation of the antisymmetric DM interaction  is carried to second order , a symmetric anisotropy results that is described by a special case of the matrices (13) with

$`K_1=0,K_2={\displaystyle \frac{D^2}{2J}},`$ (35)
(36)
$`K_3={\displaystyle \frac{D^2}{2J}},K_4={\displaystyle \frac{DD^{}}{2J}},`$ (37)
in addition to a simple renormalization of the exchange constant $`J`$. The parameter $`\kappa _0`$ of Eq. (26) is then given by $`\kappa _0=\lambda ^2\lambda ^2`$ and the parameter $`\kappa `$ of Eq. (29) vanishes. Since a nonzero $`\kappa `$ is allowed by symmetry, we shall keep it throughout our theoretical development. However, our numerical demonstrations will also be restricted to the KSEA limit ($`\kappa =0`$).

Finally, the term $`(𝐡\times 𝐝_z)𝐧`$ in the potential $`V`$ of Eq. (32) is absent from the energy functional of Zheludev et al. . A contribution of that nature is present in the early work of Andreev and Marchenko  and plays a significant role in various aspects of weak ferromagnetism . This term vanishes when the field is applied along the $`c`$-axis ($`𝐡\times 𝐝_z=0`$) and thus does not affect the analysis of the CI transition. However, such a term is important in the case of an in-plane magnetic field which is also of experimental interest  and is briefly discussed in the concluding paragraph of Sec. III.

## III Ground state

An important first step in the calculation of the $`T=0`$ dynamics is the search for the classical spin configuration that minimizes the static energy

$$W=V𝑑x𝑑y,$$
(38)
where $`V`$ is the potential of Eq. (32). For a field applied along the $`c`$-axis, $`𝐡=(0,0,h)`$, the potential is given by

$`V`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\left(_1𝐧\right)^2+\left(_2𝐧\right)^2+\gamma ^2n_3^2+\lambda ^2\right]`$ (40)
$`\lambda \left[\left(_1n_1_2n_2\right)n_3\left(n_1_1n_2_2\right)n_3\right],`$
which depends only on the parameter $`\lambda `$ that measures the strength of the in-plane component of the DM anisotropy, and the combination of parameters

$$\gamma ^2=\kappa +\lambda ^2+h^2$$
(41)
that includes the external field $`h`$. A notable feature of the potential (40) is its invariance under the simultaneous transformations

$`x+iy(x+iy)e^{i\psi _0},`$ (42)
$`n_1+in_2(n_1+in_2)e^{i\psi _0}.`$ (43)
This is a peculiar realization of $`U(1)`$ symmetry in that the usual 2D rotation of spatial coordinates with an angle $`\psi _0`$ is followed by an azimuthal rotation of the staggered magnetization with an angle $`\psi _0`$.

The minimization problem was extensively studied in the earlier work . Here we briefly describe a slightly simplified version of the obtained solution in order to establish convenient notation for our subsequent dynamical calculations. If we invoke the usual spherical parametrization of the unit vector $`𝐧`$ defined from

$$n_1+in_2=\mathrm{sin}\mathrm{\Theta }e^{i\mathrm{\Phi }},n_3=\mathrm{cos}\mathrm{\Theta },$$
(44)
the minimum of the energy is sought after in the form of the one-dimensional (1D) Ansatz

$$\mathrm{\Theta }=\theta (x),\mathrm{\Phi }=0,$$
(45)
which assumes that the staggered magnetization is confined in the $`xz`$-plane and depends only on the spatial coordinate $`x`$, modulo a $`U(1)`$ transformation given by Eq. (42). The potential (40) then simplifies to

$$V=\frac{1}{2}\left[\left(\theta ^{}\lambda \right)^2+\gamma ^2\mathrm{cos}^2\theta \right],$$
(46)
where the prime denotes differentiation with respect to $`x`$, and stationary points of the energy (38) satisfy the ordinary differential equation $`\theta ^{\prime \prime }+\gamma ^2\mathrm{cos}\theta \mathrm{sin}\theta =0`$ whose distinct feature is that it does not depend on $`\lambda `$. A first integral of this equation is given by $`\theta ^2\gamma ^2\mathrm{cos}^2\theta =C=\delta ^2`$, where we anticipate the fact that the minimum of the energy is achieved at positive integration constant $`C`$. Thus the desired solution $`\mathrm{\Theta }=\theta (x)`$ is given by the implicit equation

$$x=_0^\theta \frac{d\vartheta }{\sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\vartheta }},$$
(47)
and is a monotonically increasing function of $`x`$. The corresponding spin structure repeats itself when $`\theta `$ is changed by an amount $`2\pi `$ ; i.e., when $`x`$ advances by a distance

$$L=4_0^{\frac{\pi }{2}}\frac{d\theta }{\sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }},$$
(48)
which will be called the period of the spiral. The free parameter $`\delta `$ is determined by the requirement that the average energy density $`w=\frac{1}{L}_0^LV𝑑x`$ is a minimum, where $`V`$ is the potential (46) calculated for the specific configuration (47). A direct computation shows that $`\delta `$ must satisfy the algebraic equation

$$\frac{2}{\pi }_0^{\frac{\pi }{2}}𝑑\theta \sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }=\lambda ,$$
(49)
and the corresponding energy density is

$$w=\frac{1}{2}(\lambda ^2\delta ^2).$$
(50)
The configuration described above will be referred to as the flat spiral because the staggered magnetization is confined in the $`xz`$-plane.

It is clear that the root $`\delta `$ of Eq. (49) decreases with increasing $`\gamma `$. In fact, $`\delta `$ vanishes at a critical value of $`\gamma `$ which is easily calculated by setting $`\delta =0`$ in Eq. (49) to obtain $`\gamma =\gamma _c=\lambda \pi /2`$. In view of Eq. (41), the corresponding critical field is given by

$$h_c=\left[\left(\frac{\pi }{4}^21\right)\lambda ^2\kappa \right]^{1/2},$$
(51)
and a spiral state is possible only for $`h<h_c`$. At the critical point, the energy density (50) becomes $`w=\lambda ^2/2`$ and is equal to the energy of the uniform spin-flop state $`𝐧=(1,0,0)`$. The latter is a stationary point of the energy functional for any strength of the applied field and is thought to be the absolute minimum for $`h>h_c`$. The actual stability of the spin-flop state for $`h>h_c`$, and of the spiral state for $`h<h_c`$, will be addressed more carefully in Secs. IV and V.

Next we calculate the $`T=0`$ magnetization $`𝐦=(m_1,m_2,m_3)`$ which can be obtained from Eq. (31) applied for the static configuration $`𝐧=(\mathrm{sin}\theta ,0,\mathrm{cos}\theta )`$ and averaged over the period $`L`$ of the spiral. The only term that survives in the average is

$$m_3=\frac{\epsilon h}{2\sqrt{2}}\frac{1}{L}_0^L(1\mathrm{cos}^2\theta )𝑑x,$$
(52)
and can be expressed in terms of quantities already considered, namely

$$m_3=\frac{\epsilon h}{2\sqrt{2}\gamma ^2}\left(\gamma ^2+\delta ^2\frac{2\pi \lambda }{L}\right),h<h_c.$$
(53)
For $`h>h_c`$, the spin-flop state $`𝐧=(1,0,0)`$ is inserted in Eq. (31) to yield after a trivial computation

$$m_3=\frac{\epsilon h}{2\sqrt{2}},h>h_c,$$
(54)
while $`m_1=0`$ and $`m_2=\epsilon \lambda ^{}/2\sqrt{2}`$. The latter formula is the only place where the oscillating component of the DM anisotropy appears and produces a field-independent weak ferromagnetic moment along the $`y`$-axis.

In order to make definite quantitative predictions we use as input  the spin value $`s=1/2`$, an exchange constant $`J=0.96`$ meV, and a gyromagnetic ratio $`g=g_c=2.474`$ for a field applied along the $`c`$-axis. Concerning anisotropy, we adopt the KSEA limit ($`\kappa =0`$) and thus the only relevant parameter is $`\lambda `$ which may be estimated from the observed spin rotation by an angle $`\mathrm{\Delta }\theta 2\pi \zeta `$ over a distance $`d=a/\sqrt{2}`$ along the $`x`$-axis. The incommensurability parameter $`\zeta `$ is related to the period $`L`$ of Eq. (48) by $`\zeta =\epsilon /L`$, where $`\epsilon `$ is the scale parameter introduced in Eq. (24). One may actually choose the free parameter $`\epsilon `$ as $`\epsilon =D/J`$ and thus $`\lambda 1`$ and $`\gamma ^2=\kappa +\lambda ^2+h^2=1+h^2`$. At zero field, Eq. (49) is applied for $`\lambda =1=\gamma `$ to yield $`\delta ^2=0.53189772`$ and the period is calculated from Eq. (48) as $`L=6.49945169`$. Hence, $`\epsilon =\zeta L=0.1774`$, where we have also used the value $`\zeta =0.0273`$ measured at zero field . To summarize, our final choice of constants is

$`\kappa `$ $`=`$ $`0,\lambda 1,`$ (55)
$`\gamma ^2`$ $`=`$ $`1+h^2,\epsilon =D/J=0.1774,`$ (56)
and should be completed with the stipulation that the unit of field ($`h=1`$) corresponds to $`2s\sqrt{2}\epsilon J/g_c\mu _B=1.682`$ T, while the unit of frequency (energy) is $`2s\sqrt{2}\epsilon J=0.241`$ meV. The magnetization per Cu atom is given by Eqs. (53) and (54) in units of $`sg_c\mu _B=1.237\mu _B`$. Distance is measured in units of $`d/\epsilon =33.75`$ Å.

The constants (55) are inserted in Eq. (51) to yield a critical field $`h_c=1.21`$ in rationalized units or $`H_c=2.04`$ T in physical units. This theoretical prediction is consistent with experiment and is thought to be a good indication that the KSEA limit ($`\kappa =0`$) may provide an accurate description of anisotropy . Now, Eqs. (49) and (48) are applied with $`\lambda =1`$ and $`\gamma ^2=1+h^2`$ to yield the root $`\delta =\delta (h)`$ and the period $`L=L(h)`$ at field $`h`$. The field dependence of the energy density computed from Eq. (50) for $`h<h_c`$, and $`w=1/2`$ for $`h>h_c`$, is depicted by a solid line in Fig. 3a. Similarly, the field dependence of the incommensurability parameter $`\zeta =\zeta (h)`$ is calculated from

$$\frac{\zeta (h)}{\zeta (0)}=\frac{L(0)}{L(h)},$$
(57)
where $`\zeta (0)`$ and $`L(0)`$ are the zero-field parameters already discussed, and is depicted by a solid line in Fig. 3b. The results of Fig. 3 will be completed and further discussed in Sec. VI. The same numerical data may be employed in Eqs. (53) and (54) to calculate the field dependence of the magnetization and the corresponding susceptibility.

Finally, we return to the $`U(1)`$ transformation (42) which may be applied to the special solution (45) to yield a family of degenerate ground-state configurations:

$$\mathrm{\Theta }=\theta (x\mathrm{cos}\psi _0+y\mathrm{sin}\psi _0),\mathrm{\Phi }=\psi _0,$$
(58)
where $`\psi _0`$ is an arbitrary angle. The propagation vector of the resulting spiral forms an angle $`\psi _0`$ with the $`x`$-axis, while the normal to the spin plane forms an angle $`\pi /2\psi _0`$ with the same axis. For the special rotation $`\psi _0=\pi /4`$, the magnetic propagation vector and the normal to the spin plane are parallel (screw-type spiral). This symmetry operation is the basis for the bisection rule discovered by Zheludev et al.  when the external field is applied in a direction *perpendicular* to the $`c`$-axis, at an angle $`\chi _0`$ with respect to the $`x`$-axis. The normal to the spin plane rotates almost freely to align with the external field, and thus $`\chi _0=\pi /2\psi _0`$, in order to minimize (eliminate) the positive term $`(𝐧𝐡)^2`$ in the potential (32). The new term $`(𝐡\times 𝐝_z)𝐧`$ in the above potential does not affect the bisection rule but it does modify the profile of the spiral. For example, when the field is applied along the $`y`$-axis ($`\chi _0=\pi /2`$, $`\psi _0=0`$, $`𝐡\times 𝐝_z=\lambda ^{}h𝐞_1`$) the staggered magnetization is again confined in the $`xz`$-plane but the potential (46) becomes

$$V=\frac{1}{2}\left[\left(\theta ^{}\lambda \right)^2+\gamma ^2\mathrm{cos}^2\theta \right]+\lambda ^{}h\mathrm{sin}\theta ,$$
(59)
where $`\gamma ^2=\kappa +\lambda ^2`$ is now field independent. Nevertheless, the external field reappears in a different form and requires a new calculation of the spiral based on Eq. (59). Such a calculation might actually explain the observed (weak) field dependence of the magnitude of the magnetic propagation vector  and provide an estimate for the strength $`\lambda ^{}`$ (or $`D^{}`$) of the oscillating component of the DM anisotropy.

## IV Spin-flop phase

We now begin to address questions of dynamics based on the complete Lagrangian $`=_0V`$ of Eq. (32) applied for a field $`𝐡=(0,0,h)`$. If we also insert the spherical parameters (44) we find that

$$_0=\frac{1}{2}\left(\dot{\mathrm{\Theta }}^2+\mathrm{sin}^2\mathrm{\Theta }\dot{\mathrm{\Phi }}^2\right)+h\mathrm{sin}^2\mathrm{\Theta }\dot{\mathrm{\Phi }},$$
(60)
and

$`V=`$ $`{\displaystyle \frac{1}{2}}\left[\left(\mathrm{\Theta }\right)^2+\mathrm{sin}^2\mathrm{\Theta }\left(\mathrm{\Phi }\right)^2+\gamma ^2\mathrm{cos}^2\mathrm{\Theta }+\lambda ^2\right]`$ (63)
$`+\lambda [\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\mathrm{\Theta }(\mathrm{sin}\mathrm{\Phi }_1\mathrm{\Phi }+\mathrm{cos}\mathrm{\Phi }_2\mathrm{\Phi })`$
$`\mathrm{cos}\mathrm{\Phi }_1\mathrm{\Theta }+\mathrm{sin}\mathrm{\Phi }_2\mathrm{\Theta }],`$
where $`=(_1,_2)`$ is the usual 2D gradient operator, while the Laplacian will be denoted in the following by $`\mathrm{\Delta }=_1^2+_2^2`$.

We first study the high-field commensurate phase ($`h>h_c`$) where the absolute minimum of the classical energy is thought to be the uniform spin-flop state $`𝐧=(1,0,0)`$ or $`\mathrm{\Theta }=\pi /2`$ and $`\mathrm{\Phi }=0`$. Small fluctuations around this state are calculated by introducing $`\mathrm{\Theta }=\pi /2+f`$ and $`\mathrm{\Phi }=g`$ in Eqs. (60) and (63) and keeping terms that are at most quadratic in the small amplitudes $`f=f(x,y,\tau )`$ and $`g=g(x,y,\tau )`$. Linear terms do not appear because we are expanding around a stationary point of the energy functional, whereas constants and total derivatives can be omitted because they do not contribute to the equations of motion. Thus the corresponding linearized equations are found to be

$$\ddot{f}\mathrm{\Delta }f+\gamma ^2f=2\lambda _2g,\ddot{g}\mathrm{\Delta }g=2\lambda _2f.$$
(64)
Performing the usual Fourier transformation with frequency $`\omega `$ and wave vector $`𝐪=(q_1,q_2)`$ one obtains a homogeneous system whose solution requires that the corresponding determinant vanish. This condition leads to two branches of eigenfrequencies

$$\omega _\pm (𝐪)=\left[q_1^2+q_2^2+\frac{1}{2}\left(\gamma ^2\pm \sqrt{\gamma ^4+16\lambda ^2q_2^2}\right)\right]^{1/2},$$
(65)
which will be referred to as the optical or acoustical mode, corresponding to the plus or minus sign, respectively.

A notable feature of the calculated dispersions is their strong anisotropy. In particular, the low-$`𝐪`$ acoustical branch reads

$$\omega _{}(𝐪)\left[q_1^2+\left(14\lambda ^2/\gamma ^2\right)q_2^2\right]^{1/2}$$
(66)
and demonstrates that the spin-wave velocity depends on the direction of propagation. It also makes it clear that an instability arises when $`\gamma ^2<4\lambda ^2`$. In fact, the complete acoustical frequency of Eq. (65) becomes purely imaginary over a nontrivial region in $`𝐪`$-space when $`\gamma ^2=\kappa +\lambda ^2+h^2<4\lambda ^2`$. Therefore, the uniform spin-flop state is unstable for $`h<h_2`$ where

$$h_2=\sqrt{3\lambda ^2\kappa }$$
(67)
is a new critical field. For our choice of parameters (55) $`h_2=\sqrt{3}`$ or $`H_2=1.682\sqrt{3}=2.91`$ T. The important conclusion is that the spin-flop state is locally stable only for $`H>H_2>H_c`$.

It is also interesting to examine the gap of the optical branch at $`𝐪=0`$ where $`\omega _+(𝐪=0)=\gamma =(\kappa +\lambda ^2+h^2)^{1/2}`$. This result may be used to illustrate our earlier claim concerning the role of the scale parameter $`\epsilon `$. If we recall the definition of the rescaled parameters (26) and also include the factor $`2s\sqrt{2}\epsilon J`$ to account for the physical unit of frequency, the calculated gap is independent of $`\epsilon `$ and is expressed entirely in terms of constants that appear in the original discrete Hamiltonian of Sec. II. Hence, in the KSEA limit, we find that

$$\omega _+(𝐪=0)=\left[\left(2s\sqrt{2}D\right)^2+\left(g_c\mu _BH\right)^2\right]^{1/2},$$
(68)
in agreement with the magnon gap given in Ref. 5. Incidentally, this special result is the only feature of the spectrum actually calculated in the above reference for nonzero field.

The complete dispersions are illustrated in Fig. 4 for $`H=3`$ T, and for spin-wave propagation along the $`x`$\- or the $`y`$-axis.

The anisotropy of the spectrum is made especially apparent by the fact that the dispersion of the acoustical mode is strictly linear in the $`x`$-direction, but almost ferromagnetic-like in the $`y`$-direction because the chosen field is only slightly greater than the critical field $`H_22.9`$ T. The numerical data for Fig. 4 were obtained from Eq. (65) applied for our choice of units and constants given in Eq. (55). Thus we set $`\lambda =1`$ and $`\gamma ^2=1+h^2`$, with $`h=3/1.682=1.784`$, and also include an overall factor $`0.241`$ meV to account for the physical unit of energy. Finally, $`𝐐=\epsilon 𝐪`$ is the wave vector defined on the complete square lattice formed by the Cu atoms within each layer, while relative units are defined from $`Q[r.l.u.]=Q/2\pi =(\epsilon /2\pi )q`$. Therefore, Eq. (65) is applied with $`q=(2\pi /\epsilon )Q[r.l.u.]=35.418Q[r.l.u.].`$

Unfortunately, there seem to exist no experimental data in the field region $`H3`$ T. In fact, the only published data  were obtained for $`H=2.5`$ T $`<H_2`$ and spin-wave propagation along the $`x`$-axis. For this special direction ($`q_2=0`$) the theoretical dispersions (65) do not “see” the instability. One may then deliberately apply them for $`H=2.5`$ T and compare the results to the actual data, as is done in Fig. 5 where a systematic disagreement is apparent in both dispersions. In particular, the numerical fits to the data represented by dashed lines indicate a significant 20% reduction in the measured spin-wave velocity; as was already noted in Ref. 5.

Of course, our earlier discussion makes it clear that the dispersions (65) cannot be applied for $`H=2.5`$ T because the corresponding ground state is predicted to be unstable. At best, the fully polarized spin-flop state $`𝐧=(1,0,0)`$ survives in the field region $`H<H_2`$ as a metastable state thanks to some small tetragonal anisotropy that may be present in the discrete system  but drops out of the leading continuum approximation. An appealing scenario suggested by our calculation is that the system actually enters a different (intermediate) phase for $`H<H_2`$ which consists of some sort of mixed domains with no definite axis of polarization. Such a picture could explain the effective reduction of the spin-wave velocity, also taking into account the anisotropy of the acoustical mode.

As mentioned already, the continuum model does not contain anisotropies that would necessarily polarize the staggered magnetization along the $`x`$ (or the $`y`$) axis. Instead, there is a family of degenerate spin-flop states $`𝐧=(\mathrm{cos}\mathrm{\Phi }_0,\mathrm{sin}\mathrm{\Phi }_0,0)`$ with the same energy for any constant angle $`\mathrm{\Phi }_0`$. The corresponding small fluctuations are now studied by introducing $`\mathrm{\Theta }=\pi /2+f`$ and $`\mathrm{\Phi }=\mathrm{\Phi }_0+g`$ in Eqs. (60) and (63). A short calculation similar to the one presented for $`\mathrm{\Phi }_0=0`$ leads to the magnon dispersions

$$\omega _\pm (𝐪)=\left[𝐪^2+\frac{1}{2}\left(\gamma ^2\pm \sqrt{\gamma ^4+16\lambda ^2(𝐞𝐪)^2}\right)\right]^{1/2},$$
(69)
where $`𝐞=\mathrm{sin}\mathrm{\Phi }_0𝐞_1+\mathrm{cos}\mathrm{\Phi }_0𝐞_2`$ is the unit vector obtained by rotating $`𝐞_2`$ with an angle $`\mathrm{\Phi }_0`$. The emerging picture is yet another manifestation of the peculiar nature of the $`U(1)`$ symmetry (42), in some respects similar to the bisection rule discussed in the concluding paragraph of Sec. III. In any case, the main conclusion of the present section persists; namely, the acoustical mode develops maximum instability along the direction $`𝐞`$ and leads to the same critical field given earlier in Eq. (67).

The nature of the intermediate phase will be discussed in Sec. VI. The present section is concluded with a word of caution concerning the validity of the continuum approximation at nonzero field, which roughly requires that $`g_c\mu _BHJ`$. This strong inequality becomes increasingly marginal for field strengths in the region $`HH_2`$.

## V Spiral phase

The calculation of the low-energy magnon spectrum in the spiral phase ($`h<h_c`$) is significantly more complicated, but the general strategy is identical to that followed in Sec. IV. Hence we introduce new fields according to

$$\mathrm{\Theta }=\theta +f,\mathrm{\Phi }=\frac{g}{\mathrm{sin}\theta },$$
(70)
where $`\theta =\theta (x)`$ is the profile of the ground-state spiral given by Eq. (47) while $`f=f(x,y,\tau )`$ and $`g=g(x,y,\tau )`$ account for small fluctuations. The special rescaling chosen in the second equation is equivalent to working in a *rotating frame*  whose third axis is everywhere parallel to the direction of the background staggered magnetization $`𝐧=(\mathrm{sin}\theta ,0,\mathrm{cos}\theta )`$.

The new fields (70) are introduced in the complete Lagrangian given by Eqs. (60) and (63) which is then expanded to second order in $`f`$ and $`g`$. The required algebra is lengthy but the final result for the linearized equations is sufficiently simple:

$`\ddot{f}\mathrm{\Delta }f+U_1f`$ $`=`$ $`2h\mathrm{cos}\theta \dot{g}+2\lambda \mathrm{sin}\theta _2g,`$ (71)
$`\ddot{g}\mathrm{\Delta }g+U_2g`$ $`=`$ $`2h\mathrm{cos}\theta \dot{f}2\lambda \mathrm{sin}\theta _2f,`$ (72)
where

$`U_1`$ $`=`$ $`\gamma ^2\mathrm{cos}(2\theta ),`$ (73)
$`U_2`$ $`=`$ $`2\lambda \sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }2\gamma ^2\mathrm{cos}^2\theta \delta ^2,`$ (74)
are effective potentials that can be calculated for any desired set of parameters, as explained in Sec. III. The general idea that the calculation of the spectrum in a spiral antiferromagnet can be reduced to a Schrödinger-like problem in a periodic potential is not new , but the specific structure of Eqs. (71) requires special attention.

We found it instructive to consider first the special case of spin-wave propagation along the $`x`$-axis ($`_2f=0=_2g`$) at zero external field ($`h=0`$). This is actually the only case for which the low-energy spectrum was previously calculated starting from the discrete Hamiltonian . If we further perform the temporal Fourier transformation with frequency $`\omega `$, Eqs. (71) reduce to

$$f^{\prime \prime }+U_1f=\omega ^2f,g^{\prime \prime }+U_2g=\omega ^2g,$$
(75)
where the prime denotes differentiation with respect to $`x`$. Therefore, in this special case, the eigenvalue problem is reduced to two decoupled 1D Schrödinger equations of the standard type with potentials $`U_1`$ and $`U_2`$ calculated at zero field. Also note that both potentials are periodic functions of $`2\theta `$ and thus their period is actually $`L/2`$ where $`L`$ is the period of the background spiral.

The eigenvalue problems (75) are solved in Appendix A. The numerical procedure yields eigenfrequencies $`\omega =\omega (q_1)`$ as functions of Bloch momentum $`q_1`$. The latter can be restricted to the zone $`[2\pi /L,2\pi /L]`$, because the period of the potentials is $`L/2`$, or to the zone $`[\zeta ,\zeta ]`$ in relative units defined as in Sec. IV. Several low-lying eigenvalues are illustrated in Fig. 6a using a reduced-zone scheme. Solid and dashed lines correspond to the first and second eigenvalue problem in Eq. (75) and are superimposed in the same graph for convenience. We also find it convenient to refer to the two types of modes as acoustical and optical. In either case, there is only one discernable gap that occurs between the first and the second band at the zone boundary. The calculated boundary gaps are $`0.123`$ meV and $`0.049`$ meV, respectively, while the absolute gap of the optical mode at the zone center is $`0.170`$ meV. All of the above theoretical predictions agree with those obtained in Ref. 5 by a different method. They also agree with experiment, except for the small ($`0.049`$ meV) gap that has not yet been resolved at zero field.

The same results are depicted in Fig. 6b using an extended-zone scheme. In fact, this figure displays two replicas of the acoustical mode centered at $`\pm \zeta `$. The need for two replicas follows from the structure of dynamic correlation functions in the laboratory frame, rather than in the rotating frame actually used in the calculation of the magnon spectrum . Our results in Fig. 6b are obviously consistent with both the experimental and the theoretical results obtained in the above reference at zero field.

We are now in a position to extend the calculation to the general case of nonzero field and arbitrary direction of spin-wave propagation. The external field enters Eqs. (71) in two distinct ways. First, it affects the structure of the potentials $`U_1`$ and $`U_2`$ because the background spiral is further distorted. Second, the field induces first-order time derivatives which originate in the “nonrelativistic” term of Eq. (60) and couple the two linear equations (71). Additional coupling between the two equations appears in the case of arbitrary direction of propagation because $`_2f`$ and $`_2g`$ no longer vanish. Altogether we are faced with a nonstandard eigenvalue problem that is also solved in Appendix A.

Here we present explicit results for four typical values of the rationalized field $`h=0`$, $`0.3`$, $`0.6`$ and $`0.9`$ which will be quoted from now on by their rounded physical values $`H=0`$, $`0.5`$, $`1`$ and $`1.5`$ T. In Fig. 7 we illustrate the calculated spectrum for spin-wave propagation along the $`x`$-axis ($`q_2=0`$) using a highly reduced zone scheme. An important check of consistency is provided by the fact that the $`H=0`$ results of Fig. 7 agree with those presented earlier in Fig. 6a, except that the zone is now reduced down to \[$`\zeta /2,\zeta /2`$\] for reasons explained in Appendix A. Furthermore, we no longer employ solid and dashed lines to distinguish between acoustical and optical modes. Such a distinction is not a priori possible in the current algorithm because of the coupling (hybridization) of the two types of modes at nonzero field.

One should keep in mind that the extent of the zone \[$`\zeta /2,\zeta /2`$\] slides with the applied field, a feature that is not apparent in Fig. 7 because the scale of the abscissa is adjusted accordingly. The incommensurability parameter $`\zeta =0.0273`$ measured at $`H=0`$ is used as input in our calculation. The calculated values for $`H=0.5`$, $`1`$ and $`1.5`$ T are $`\zeta =0.0271`$, $`0.0264`$ and $`0.0245`$.

At first sight, it would seem difficult to extract useful information from the highly convoluted spectra shown in Fig. 7. Nevertheless, the most vital information concerning the low-energy dynamics is easily abstracted from Fig. 7 because the low-lying bands are clearly segregated. In particular, it is still possible to distinguish between the acoustical and the optical mode, at least in an operational sense. Thus we unfold the first six branches back to the zone \[$`\zeta ,\zeta `$\] and then proceed to the extended-zone scheme of Fig. 6b including two replicas of the acoustical mode centered at $`\pm \zeta `$. The resulting low-energy spectra are shown in Fig. 8.

The $`H=0`$ entry of Fig. 8 is but a magnified version of the lower-central portion of Fig. 6b, as expected. This version reveals a certain “anomaly” that is not conspicuous in Fig. 6b, namely, a relative crossing between the two modes in a narrow region around the zone center. The calculated maximum splitting of $`0.005`$ meV is within the error margin of the continuum approximation and, in any case, beyond experimental detection. But the resolution of this theoretical curiosity is interesting: when the direction of spin-wave propagation departs slightly from the $`x`$-axis ($`q_20`$) and/or a finite field is turned on, the crossing points become avoided crossings. Therefore, strictly speaking, the solid and dashed lines must be interchanged in the narrow region between the two crossing points. This explains the apparent slight inconsistency in the labeling of the five characteristic points of the spectrum denoted by 1, 2, 3, 4 and 5 in Fig. 8. The calculated magnon energies at those points are summarized in Table I.

We now concentrate on the optical mode. The gap $`E_2=0.176`$ meV calculated at zero field agrees with the measured $`0.18`$($`1`$) meV. Our calculation further shows that the above gap evolves quickly with increasing field to reach the asymptotic value $`0.26`$ meV around which it oscillates mildly. The complete optical mode evolves into a snake-like dispersion with energy values in the range $`0.25`$ meV $`<E<0.29`$ meV. These predictions are generally consistent with experiment . However, some of the finer details deserve closer attention. The calculated energy at point 5 in the spectrum remains practically constant at $`E_50.31`$ meV for $`H1`$ T, while a steep crossover takes place for higher field values which leads to $`E_50.35`$ meV for $`H=1.5`$ T. These predictions are also in agreement with experiment . But the calculated splittings of the optical dispersion $`E_5E_3=0.02`$ meV and $`0.06`$ meV, for $`H=1`$ and $`1.5`$ T, disagree with the measured 0.05 meV and 0.11 meV. It appears that the observed splittings are better described by $`E_5E_2=0.04`$ meV and 0.09 meV. In fact, the above identification may not be completely arbitrary. For instance, the lowest branch in the optical dispersion measured for $`H=1.5`$ T shows a clear local maximum of 0.28 meV at the zone center, which agrees with the calculated maximum $`E_3=0.285`$ meV at the zone boundaries $`\pm \zeta `$ rather than the gap $`E_2=0.255`$ meV at the zone center. It seems that the lowest branch in the observed optical dispersion for $`H=1.5`$ T is composed of two replicas of the calculated dispersion centered at $`\pm \zeta `$. On the other hand, experimental data  at higher energies not shown in Fig. 8 indicate the appearance of two replicas centered at $`\pm 2\zeta `$. Unfortunately, we cannot resolve this issue of proper replication of the basic modes because our current formalism does not directly address the relevant dynamic correlation functions.

Next we discuss the acoustical mode. Our calculation shows that the energy at point 4 in the spectra of Fig. 8 remains remarkably stable at $`E_40.30`$ meV for all field values considered. This feature is also in agreement with experiment which indicates only a mild decline from the above value with increasing field. Nevertheless, a clear disagreement occurs in the lowest branch of the acoustical mode. Although explicit data points are not given for this branch by Zheludev et al. , the solid lines in their Figs. 6 and 7, and the corresponding wording in their text, suggest that the lowest branch in the measured spectrum is also largely insensitive to the applied field. In contrast, our calculation predicts a robust reduction of the energy gap $`E_1`$ with increasing field (see Fig. 8 and Table I). The calculated spin-wave velocity is also reduced, albeit at a slower rate.

The preceding apparent disagreement with experiment is especially important because it is directly related to the issue of local stability of the spiral phase. Indeed, a careful numerical investigation reveals that the gap $`E_1`$ vanishes at the critical field $`h_11.01`$, or $`H_11.70`$ T, while an unstable mode develops for $`H>H_1`$. This mode is first detected by the appearance of a *real* eigenvalue in the matrix $`M`$ of Eq. (A5), when $`H`$ crosses $`H_1`$, which corresponds to purely imaginary frequency. As the field increases beyond $`H_1`$ the instability occurs over a nontrivial region in $`𝐪`$-space. Therefore, the flat spin spiral constructed in Sec. III is predicted to be locally stable only for $`H<H_1<H_c`$.

It is interesting that the experimental work  already provided evidence for the existence of a critical field $`H_1=1.7`$ T that coincides with our theoretical prediction. However, one should also contemplate the possibility that such a coincidence may be fortuitous, in view of the apparent contradiction between experimental and theoretical predictions for the gap $`E_1`$. In any case, our current result together with the discussion of Sec. IV clearly suggest the existence of an intermediate phase in the field region $`1.7`$ T $`<H<2.9`$ T. The nature of the intermediate phase is discussed in Sec. VI.

In the remainder of this section we take a different view of the low-energy magnon spectrum by considering spin-wave propagation along the normal to the plane of the flat spiral. Our algorithm is adapted to this case simply by setting the Bloch wave number $`q_1=0`$ and calculating frequencies as functions of the wave number $`q_2`$ in the $`y`$-direction. It is interesting that no theoretical or experimental results exist in this case even at zero field. Our results are illustrated in Fig. 9 for the same set of field values employed in the preceding discussion.

The most stable feature of Fig. 9 is its lowest branch which exhibits *quadratic* dependence on $`q_2`$ near the origin. Clearly this branch is the extension of the acoustical dispersion in the $`y`$-direction originating at its points where $`E=0`$. Therefore, the complete acoustical mode is Goldstone-like in the $`x`$-direction but ferromagnetic-like in the $`y`$-direction. Such a characteristic anisotropy is in some respects similar to the situation encountered in the spin-flop phase discussed in Sec. IV.

Higher branches labeled as 1, 2, 3, 4 and 5 in Fig. 9 also possess a simple interpretation, for they are the extensions in the $`y`$-direction of the special spectral points numbered accordingly in our earlier Fig. 8. In contrast to the fundamental ferromagnetic-like branch, higher branches evolve vigorously with the applied field. In particular, branch 1 in Fig. 9 is quickly depressed with increasing field to become degenerate with the fundamental branch at the critical field $`H_1=1.70`$ T not included in the figure. For $`H>H_1`$ this mode becomes unstable over a nontrivial region of wave numbers around the origin. Of course, this is the instability described earlier in the text viewed from a different perspective.

We have thus provided a fairly complete theoretical picture of the low-energy magnon spectrum, including predictions for which there exist no experimental data at present. It is interesting to see whether or not future experiments could resolve the apparent discrepancy in the field dependence of the magnon gap $`E_1`$, and thus illuminate the important issue of local stability of the spiral phase, as well as confirm the predicted characteristic anisotropy in the low-energy spectrum.

## VI Intermediate phase

We now focus on the predicted intermediate phase and examine its nature through a direct numerical minimization of the complete energy functional $`W`$ of Eqs. (38) and (40). The method of calculation is a relaxation algorithm formulated on the basis of a discretized form of the energy functional defined on a square grid. After long experimentation with 2D simulations, it progressively became apparent that the optimal configuration for $`h>h_1`$ is actually a 1D nonflat spiral characterized by a staggered magnetization whose three components are all different than zero.

Therefore, an accurate calculation of the nonflat spiral was eventually obtained by a relaxation algorithm applied directly to a 1D restriction of the energy functional whose variation leads to the coupled stationary equations:

$`\mathrm{\Theta }^{\prime \prime }+(\gamma ^2\mathrm{\Phi }^2)\mathrm{cos}\mathrm{\Theta }\mathrm{sin}\mathrm{\Theta }=2\lambda \mathrm{sin}^2\mathrm{\Theta }\mathrm{sin}\mathrm{\Phi }\mathrm{\Phi }^{}`$ (76)
$`(\mathrm{sin}^2\mathrm{\Theta }\mathrm{\Phi }^{})^{}=2\lambda \mathrm{sin}^2\mathrm{\Theta }\mathrm{sin}\mathrm{\Phi }\mathrm{\Theta }^{}`$ (77)
These are ordinary differential equations because both angular variables $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ are assumed to be functions of the single coordinate $`x`$, while the prime again denotes differentiation with respect to $`x`$. Nevertheless, it does not seem possible to obtain analytical solutions of Eqs. (76), except for the case of the flat spiral ($`\mathrm{\Phi }=0`$) discussed in Sec. III. A significant obstacle is the fact that the period of the nonflat spiral is not known a priori. Hence our numerical solution was carried out on a periodic 1D grid with specified length $`L`$, until a relaxed configuration was obtained with energy density $`w=w(L)`$. We then varied $`L`$ to achieve the least possible energy for each field $`h`$, and the corresponding optimal period $`L=L(h)`$.

An important check of consistency is that the above algorithm reproduces the results for the flat spiral obtained more directly in Sec. III, but only when $`h<h_1=1.01`$. Instead, a nonflat spiral emerges as the optimal solution for $`h>h_1`$. The calculated configuration is illustrated in Fig. 10 for a field value $`h=1.21`$ deliberately chosen to be equal to the critical field $`h_c`$ of the conventional CI transition.

The energy of the nonflat spiral depicted by a dashed line in Fig. 3a is smaller than the energy of both the flat spiral and the uniform spin-flop state throughout the intermediate region $`h_1<h<h_2`$. One should also stress that the nonflat spiral is here predicted to occur for a field applied strictly along the c-axis, and is not due to sample misalignment  or the presence of a transverse magnetic field .

In a sense, the predicted intermediate phase smooths out the original sharp CI transition. This smoothing is also apparent in the calculated field dependence of the period $`L=L(h)`$ which is inserted in Eq. (3.17) to yield the results for the incommensurability parameter shown by a dashed line in Fig. 3b. The same figure displays experimental data taken from Ref. 4 where they were analyzed in terms of the conventional CI transition based solely on a flat spiral. It should be noted that both the measured zero-field incommensurability parameter $`\zeta (0)=0.0273`$ and an experimental critical field $`H_c=2.15`$ T were used as adjustable parameters in the theoretical analysis of Refs. 4 and 5 to obtain a reasonable overall fit. Yet the experimental data indicate some smoothing of the CI transition near the critical field. This fact is made apparent in our Fig. 3b where theoretical results for both the flat spiral (solid line) and the nonflat spiral (dashed line) are calculated using as input only the zero-field parameters given earlier in Eq. (3.16).

Nevertheless, the results of Fig. 3b cannot be interpreted as unambiguous evidence for the existence of an intermediate phase, especially because the experimental data were taken at the relatively high temperature $`T=2.4`$ K. It is feasible that the $`T=0`$ theoretical predictions could be further focused by invoking deviation from KSEA anisotropy that is allowed by symmetry; i.e., by repeating the calculation for nonzero values of the free parameter $`\kappa `$. One should also keep in mind that a completely accurate description of the CI transition may not be attainable within the classical approximation.

The nonflat spiral exists as a stationary point of the energy functional throughout the intermediate phase and degenerates into a uniform spin-flop state polarized along the $`y`$-axis near the upper critical field $`h_2=1.73`$. Actually, our calculation was not pushed all the way to the critical field $`h_2`$ because of numerical difficulties that occur as the period grows to infinity. The theoretical analysis should be completed with a detailed study of the stability and dynamics of the nonflat spiral within the full 2D context, in a manner analogous to our treatment of the flat spiral in Sec. V. The required computational effort is too great to be included in the present paper, especially because the profile of the nonflat spiral is obtained numerically through the relaxation algorithm. A future analysis could, in principle, reveal the existence of yet another critical field within the intermediate region, beyond which the nonflat spiral may cease to be locally stable. It is thus important to also examine the nature of instability at the upper critical field $`h_2`$, as discussed further in Appendix B.

The configuration of Fig. 10 may be viewed as a conical spiral that nutates around the $`y`$-axis. It is interesting that a simple conical spiral without nutation had been discussed theoretically in connection with the cholesteric-nematic transition in liquid crystals  but has not yet been observed experimentally because its realization requires an anomalously small bend modulus . In contrast, the parameters of Ba2CuGe2O7 favor the occurrence of the currently predicted intermediate phase.

## VII Conclusion

We have presented a field theoretical description of the low-energy dynamics in the spiral antiferromagnet Ba2CuGe2O7. We have thus been able to calculate the low-energy magnon spectrum for any strength of the applied field and any direction of spin-wave propagation. In this respect, the present work significantly extends the results of Ref. 5 where the spectrum was calculated only at zero field and for propagation along the direction of the spiral. Therefore, our theoretical results are relevant for the analysis of experimental data obtained for nonzero field, which were previously analyzed mostly in terms of empirical formulas.

An interesting byproduct of this detailed spin-wave analysis is the identification of the two new critical fields $`H_1`$ and $`H_2`$, and a corresponding prediction of an intermediate phase that does not seem to be inconsistent with available experimental data. The apparent discrepancy in the field dependence of the magnon gap $`E_1`$ pointed out in Sec. V needs to be clarified, but could be due to poor experimental resolution at this rather low energy scale ($`0.1`$ meV or less). The field dependence of the incommensurability parameter discussed in Sec. VI could be rectified by invoking a slight deviation from the KSEA limit that is allowed by symmetry. Susceptibility data  taken at $`T=2`$ K display a rounded maximum which could be explained as a finite-temperature effect but does not a priori exclude an intermediate phase. Furthermore, the set of data for the magnon dispersion discussed in connection with Fig. 5 is too limited to provide a clear picture. Therefore, a clear identification or disproof of the intermediate phase may require additional experimental work guided by the theoretical predictions of the present paper.

On the other hand, it is desirable to carry out a complete theoretical analysis of the stability and dynamics of the intermediate phase along the lines outlined in Sec. VI. A related project is to extend our approach to the case of a field applied in a direction perpendicular to the $`c`$-axis . The field-dependent modifications of the spiral can be computed on the basis of Eq. (59), and a corresponding calculation of the low-energy magnon spectrum can be carried out by a straightforward extension of the methods developed in Sec. V.

Finally, we must comment on the two basic approximations made in the present work. The adopted classical approach is equivalent to the usual semiclassical approximation obtained by the $`1/s`$ expansion restricted to leading order. The omitted quantum (anharmonic) corrections are not negligible in this 2D problem but are offset in part by the fact that the input parameters are consistently estimated within the classical approximation . One should also question the validity of the continuum approximation whose relative accuracy can be roughly estimated from $`\epsilon ^20.03`$ at zero field, but may deteriorate in the presence of a strong external magnetic field. Incidentally, the corresponding parameter $`\epsilon `$ in a typical weak ferromagnet such as an orthoferrite (YFeO3) or a high-Tc superconductor (La2CuO4) is at least one order of magnitude smaller. In any case, the physical picture derived is sufficiently complete to provide a basis for a meaningful discussion of further refinements.

###### Acknowledgements.

We thank A. Bogdanov for bringing Refs. 20 and 21 to our attention, S. Trachanas for valuable suggestions concerning the eigenvalue problems studied in the present paper, and M. Marder for a careful reading of the manuscript. The work was supported in part by a Marie Curie Fellowship (HPMT-GH-00-00177-03), by a TMR program (ERBFMRXCT-960085), and by VEGA 1/7473/20.

## A Eigenvalue problems

The eigenvalue problems (75) were solved numerically, as explained here for the first equation. Taking into account that the period of the potential is $`L/2`$, the Bloch representation of the wave function reads

$$f(x)=e^{iq_1x}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}f_nexp(i4n\pi x/L)$$
(A1)
and the wave equation becomes

$$(q_1+4n\pi /L)^2f_n+\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}U_{1,nm}f_m=\omega ^2f_n,$$
(A2)
where the Fourier coefficients of the potential are given by

$`U_{1,n}`$ $`=`$ $`{\displaystyle \frac{2}{L}}{\displaystyle _{L/4}^{L/4}}exp(i4n\pi x/L)U_1[\theta (x)]𝑑x`$ (A3)
$`=`$ $`{\displaystyle \frac{4}{L}}{\displaystyle _0^{\pi /2}}\mathrm{cos}\left[{\displaystyle \frac{4n\pi }{L}}x(\theta )\right]{\displaystyle \frac{U_1(\theta )d\theta }{\sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }}}.`$ (A4)
Here we use the fact that $`U_1`$ is an even function of $`\theta `$ or $`x`$, and $`x=x(\theta )`$ is given by the integral (47). Thus the last step of Eq. (A3) is in effect a double integral that is computed by an adaptive Newton-Cotes algorithm. The eigenvalue equation (A2) is then solved by diagonalizing the finite matrix that results from a restriction of the indices $`m`$ and $`n`$ to the interval $`[N,N]`$ where $`N`$ can be as low as 20. To be sure, only the first few Fourier coefficients of the potential $`U_1`$ and $`U_2`$ are important, as demonstrated in Table II using as input the zero-field parameters quoted in Sec. III. The numerical procedure just described yields eigenfrequencies $`\omega =\omega (q_1)`$ as functions of Bloch momentum $`q_1`$ that can be restricted to the zone $`[2\pi /L,2\pi /L]`$, or to $`[\zeta ,\zeta ]`$ in relative units defined as in Sec. IV.

We now return to the general case of nonzero field and arbitrary direction of spin-wave propagation. We first rewrite Eqs. (71) in a form that contains only first-order time derivatives. Hence we treat $`u=\dot{f}`$ and $`v=\dot{g}`$ as independent fields and introduce the four-component spinor $`𝒳`$ defined from $`𝒳^T=(u,v,f,g)`$. Then Eqs. (71) read

$$\dot{𝒳}=M𝒳,$$
(A5)
where $`M`$ is the differential operator

$`M`$ $`=`$ $`\left[\begin{array}{cccc}0& 0& I& 0\\ 0& 0& 0& I\\ D_1& D_3& 0& D_4\\ D_3& D_2& D_4& 0\end{array}\right].`$ (A10)
Here $`D_1=\mathrm{\Delta }+U_1`$, $`D_2=\mathrm{\Delta }+U_2`$, $`D_3=2\lambda \mathrm{sin}\theta _2`$, $`D_4=2h\mathrm{cos}\theta `$, and $`I`$ is the unit operator. The chief advantage of $`M`$ is that it does not contain time derivatives. A superficial disadvantage is that $`M`$ is not a hermitian operator. In fact, Eq. (A5) suggests that the eigenvalues of $`M`$ are purely imaginary and come in pairs $`\pm i\omega `$ where $`\omega `$ is the desired physical frequency. A real eigenvalue in $`M`$ would correspond to purely imaginary physical frequency and thus indicate instability of the ground-state spiral. All of these features are explicitly realized in the following numerical calculation.

Our task is then to construct a matrix representation of the differential operator $`M`$. Attention should be paid to the fact that the Bloch theorem must now be applied with the full period $`L`$ of the spiral because of those terms in Eq. (A10) that are proportional to $`\mathrm{cos}\theta `$ and $`\mathrm{sin}\theta `$. Hence the operator $`D_1=\mathrm{\Delta }+U_1`$ is replaced by a matrix ($`D_{1,nm}`$) with elements

$$D_{1,nm}=\left[\left(q_1+2n\pi /L\right)^2+q_2^2\right]\delta _{nm}+U_{1,nm},$$
(A11)
where $`q_1`$ is now restricted to the zone $`[\pi /L,\pi /L]`$, or $`[\zeta /2,\zeta /2]`$ in relative units, while $`q_2`$ is unrestricted because the spiral depends only on $`x`$. Accordingly, the Fourier coefficients of the potential are given by

$$U_{1,n}=\frac{2}{L}_0^\pi \mathrm{cos}\left[\frac{2n\pi }{L}x(\theta )\right]\frac{U_1(\theta )d\theta }{\sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }},$$
(A12)
which differs from Eq. (A3) only in that the full period $`L`$, instead of $`L/2`$, is employed. As a result, odd coefficients in Eq. (A12) vanish, while the collection of even coefficients coincides with that obtained from Eq. (A3). The operator $`D_2`$ is treated in exactly the same way replacing $`U_1`$ with $`U_2`$. On the other hand, the operator $`D_3=2\lambda \mathrm{sin}\theta _2`$ in Eq. (A10) is replaced by $`2\lambda q_2S`$ where $`S`$ is an antisymmetric matrix whose n-th codiagonal has all its elements equal to

$$S_n=\frac{2}{L}_0^\pi \mathrm{sin}\left[\frac{2n\pi }{L}x(\theta )\right]\frac{\mathrm{sin}\theta d\theta }{\sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }},$$
(A13)
and $`D_4=2h\mathrm{cos}\theta `$ is replaced by $`2hC`$ where $`C`$ is a symmetric matrix whose n-th codiagonal has all its elements equal to

$$C_n=\frac{2}{L}_0^\pi \mathrm{cos}\left[\frac{2n\pi }{L}x(\theta )\right]\frac{\mathrm{cos}\theta d\theta }{\sqrt{\delta ^2+\gamma ^2\mathrm{cos}^2\theta }}.$$
(A14)
An interesting fact is that both $`S_n`$ and $`C_n`$ vanish for even $`n`$. The most important terms are those with $`n=\pm 1`$, whereas higher-order terms account for distortion of the spiral from its ideal shape $`\theta =\lambda x`$. Such a distortion occurs even at zero field in the presence of KSEA anisotropy.

A finite-matrix representation of the differential operator $`M`$ is then obtained by restricting the indices $`m`$ and $`n`$ to the finite interval $`[N,N]`$ where $`N`$ may again be chosen as low as 20. The resulting nonsymmetric $`4(2N+1)\times 4(2N+1)`$ matrix is diagonalized numerically to yield eigenvalues that are indeed purely imaginary and come in pairs $`\pm i\omega `$ where $`\omega =\omega (q_1,q_2)`$ is the sought after physical frequency. We have thus obtained a number of results using as input the spiral parameters $`\lambda =1`$, $`\gamma ^2=1+h^2`$, $`\delta =\delta (h)`$ and $`L=L(h)`$ calculated for each field $`h`$ as explained in Sec. III. The numerical burden is insignificant and can be carried out interactively. Explicit results are discussed in Section V.

## B Vortex states

In the original picture of the CI transition  the high-field commensurate phase is rendered unstable through domain-wall nucleation at the critical field $`h_c`$ to become a spiral phase for $`h<h_c`$. The instability at the higher field $`h_2>h_c`$ suggested by the spin-wave analysis of Sec. IV is clearly caused by 2D fluctuations. Therefore, it is conceivable that the uniform spin-flop phase is actually destabilized by nucleation of 2D vortices rather than 1D domain walls, as advocated by Bogdanov et al.  in a number of related models.

We thus search for genuinely 2D stationary points of the static energy that are compatible with $`U(1)`$ symmetry. First, we introduce the usual polar coordinates ($`r,\psi `$) from

$$x=\frac{r}{\gamma }\mathrm{cos}\psi ,y=\frac{r}{\gamma }\mathrm{sin}\psi ,$$
(B1)
where the overall rescalling by the constant $`\gamma `$ will simplify subsequent calculations. A configuration that is strictly invariant under the $`U(1)`$ transformation (42) reads

$$\mathrm{\Theta }=\theta (r),\mathrm{\Phi }=\psi ,$$
(B2)
where the minus sign in the second equation is again due to the peculiar nature of $`U(1)`$ symmetry in the present problem. Under normal circumstances, e.g., an isotropic antiferromagnet in an external field , both choices $`\mathrm{\Phi }=\psi `$ and $`\mathrm{\Phi }=\psi `$ are compatible with axial symmetry and are referred to as vortex and antivortex. Here only antivortices are possible within the axially-symmetric ansatz but will be called vortices for brevity.

When the Ansatz (B2) is introduced in the potential $`V`$ of Eq. (63) the corresponding total energy $`W=V𝑑x𝑑y`$ reads

$`W`$ $`=`$ $`\pi {\displaystyle _0^{\mathrm{}}}rdr[\left({\displaystyle \frac{d\theta }{dr}}\right)^2+{\displaystyle \frac{\mathrm{sin}^2\theta }{r^2}}+\mathrm{cos}^2\theta `$ (B4)
$`\nu ({\displaystyle \frac{d\theta }{dr}}+{\displaystyle \frac{\mathrm{cos}\theta \mathrm{sin}\theta }{r}})],`$
where $`\nu =2\lambda /\gamma `$ is the only relevant parameter in this static calculation. Also note that we have droped the additive constant term $`\lambda ^2/2`$ from the potential (63) and thus the energy of the uniform spin-flop state is set equal to zero. Variation of the energy functional (B4) with respect to the unknown amplitude $`\theta (r)`$ leads to the ordinary differential equation

$$r\frac{d^2\theta }{dr^2}+\frac{d\theta }{dr}+\left(r\frac{1}{r}\right)\mathrm{cos}\theta \mathrm{sin}\theta =\nu \mathrm{sin}^2\theta ,$$
(B5)
which reduces to the familiar equation for ordinary spin vortices in the extreme limit $`\nu =0`$. For $`\nu 0`$, solutions of Eq. (B5) exhibit slow decay at large distances, namely

$$\theta (r)\frac{\pi }{2}\frac{\nu }{r}+\mathrm{},$$
(B6)
which turns into exponential decay for $`\nu =0`$. Explicit solutions were obtained by a straightforward relaxation algorithm and are illustrated in Fig. 11 for three characteristic values of the parameter $`\nu `$= 0, 1 and 2.

One may restrict the integral in Eq. (B4) to the finite range $`0<r<R`$ and examine its behavior for large $`R`$. A short calculation taking into account the asymptotic expansion (B6) leads to

$$W=\pi \left(1\nu ^2\right)\mathrm{ln}R+finiteterms,$$
(B7)
and thus the energy exhibits the familiar logarithmic divergence. This asymptotic result demonstrates the crucial role played by the parameter $`\nu `$. For $`\nu <1`$, the energy of a single vortex is greater than the energy of the uniform spin-flop state by a logarithmically divergent quantity. This is the usual situation encountered in the case of ordinary vortices ($`\nu =0`$). The vortex energy is finite for $`\nu =1`$ and becomes again logarithmically divergent but *negative* for $`\nu >1`$. The special point $`\nu =2\lambda /\gamma =1`$ leads to the same critical field $`h_2`$ given earlier in Eq. (67).

Therefore, for $`h<h_2`$, the energy of the uniform spin-flop state can be lowered by vortex nucleation. Because of the logarithmic dependence of the energy on the size of the system, it is clear that a single vortex cannot by itself produce a thermodynamically significant effect. Instead, one should expect that a large number of vortices is created for $`h<h_2`$, probably in the form of a vortex lattice . We have actually performed several numerical experiments using the full 2D relaxation algorithm described in the beginning of Sec. VI. Although we have already obtained some “spectacular” pictures indicating the formation of a vortex lattice, we have not yet been able to lower its energy below that of the nonflat spiral. It appears that the complete (2D) energy functional displays glassy behavior in the intermediate region, which may lead to several nearly degenerate local minima.


text:
# 1 Introduction

## 1 Introduction

In this paper we shall review some recent results on analytical solutions of the Mikheyev–Smirnov–Wolfenstein (MSW) effect  for the propagation of three neutrino flavours. For two model densities which are relevant to the neutrino propagation in the sun, such results have been obtained. These are the exponential density  and the linear density . (The case of a constant density, which is of some relevance for propagation through the Earth, has also been studied recently .)

For the exponential density, the solutions for the three neutrino wave functions  can be expressed in terms of generalized hypergeometric functions, $`{}_{2}{}^{}F_{2}^{}`$ and $`{}_{3}{}^{}F_{1}^{}`$. For the linear density, the solutions can be expressed as a Fourier transform of a rather simple expression, which, in the case $`N=2`$ (two flavours) reduces to the well-known parabolic cylinder functions or confluent hypergeometric functions .

We also briefly discuss the application of the results to the Lehmann–Newton–Wu (LNW) mass matrix .

## 2 Exponential density

Neutrino propagation through a medium where the electron neutrino (denoted $`\varphi _1(r)`$) interacts differently from the others, is governed by the equation

$$\begin{array}{c}i\frac{\text{d}}{\text{d}r}\left[\begin{array}{c}\varphi _1(r)\\ \varphi _2(r)\\ \varphi _3(r)\end{array}\right]=(\left[\begin{array}{ccc}D(r)& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\hfill \\ \hfill +\frac{1}{2p}\left[\begin{array}{ccc}M_{11}^2& M_{12}^2& M_{13}^2\\ M_{21}^2& M_{22}^2& M_{23}^2\\ M_{31}^2& M_{32}^2& M_{33}^2\end{array}\right])\left[\begin{array}{c}\varphi _1(r)\\ \varphi _2(r)\\ \varphi _3(r)\end{array}\right],\end{array}$$
(1)
where the mass matrix is real and symmetric, $`M_{ji}^2=M_{ij}^2(M^2)_{ij}`$. Here, $`D(r)=\sqrt{2}G_\mathrm{F}N_e(r)`$, with $`G_\mathrm{F}`$ the Fermi constant and $`N_e(r)`$ the solar electron density.

For the sun, the density  is well approximated by an exponential,

$$N_e(r)=N_e(0)e^{r/r_0},r_00.1\times R_{}.$$
(2)
It is convenient to introduce a new radial variable: $`u=r/r_0+u_0`$, and perform a rotation on the second and third components,

$$\begin{array}{c}\left[\begin{array}{cc}\mathrm{cos}\theta _0& \mathrm{sin}\theta _0\\ \mathrm{sin}\theta _0& \mathrm{cos}\theta _0\end{array}\right]\frac{r_0}{2p}\left[\begin{array}{cc}M_{22}^2& M_{23}^2\\ M_{32}^2& M_{33}^2\end{array}\right]\hfill \\ \hfill \times \left[\begin{array}{cc}\mathrm{cos}\theta _0& \mathrm{sin}\theta _0\\ \mathrm{sin}\theta _0& \mathrm{cos}\theta _0\end{array}\right]=\left[\begin{array}{cc}\omega _2& 0\\ 0& \omega _3\end{array}\right].\end{array}$$
(3)
Eq. (1) then takes the form

$$i\frac{\text{d}}{\text{d}u}\left[\begin{array}{c}\psi _1(u)\\ \psi _2(u)\\ \psi _3(u)\end{array}\right]=\left[\begin{array}{ccc}\omega _1+e^u& \chi _2& \chi _3\\ \chi _2& \omega _2& 0\\ \chi _3& 0& \omega _3\end{array}\right]\left[\begin{array}{c}\psi _1(u)\\ \psi _2(u)\\ \psi _3(u)\end{array}\right].$$
(4)
Let $`\mu _1`$, $`\mu _2`$ and $`\mu _3`$ be the eigenvalues of the $`3\times 3`$ matrix

$$\left[\begin{array}{ccc}\omega _1& \chi _2& \chi _3\\ \chi _2& \omega _2& 0\\ \chi _3& 0& \omega _3\end{array}\right].$$
(5)
These $`\mu _j`$ are the squares of the neutrino masses multiplied by $`r_0/(2p)`$. Together with $`\omega _1`$ and $`\omega _2`$ they control the evolution of the $`\psi _i`$.

We again introduce a new variable, $`z=ie^u`$. Then, the solutions to Eq. (4) can be expressed in terms of solutions to

$$\begin{array}{c}[(z\frac{\text{d}}{\text{d}z}i\mu _1)(z\frac{\text{d}}{\text{d}z}i\mu _2)(z\frac{\text{d}}{\text{d}z}i\mu _3)\hfill \\ \hfill z(z\frac{\text{d}}{\text{d}z}i\omega _2)(z\frac{\text{d}}{\text{d}z}i\omega _3)]\psi =0,\end{array}$$
(6)
namely generalized hypergeometric functions :

$`\psi ^{(1)}=e^{i\mu _1u}`$
$`\times {}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}i(\omega _2\mu _1),& i(\omega _3\mu _1)\\ 1i(\mu _2\mu _1),& 1i(\mu _3\mu _1)\end{array}|ie^u\right]`$
$`\psi ^{(2)}=e^{i\mu _2u}`$
$`\times {}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}i(\omega _2\mu _2),& i(\omega _3\mu _2)\\ 1i(\mu _1\mu _2),& 1i(\mu _3\mu _2)\end{array}|ie^u\right]`$
$`\psi ^{(3)}=e^{i\mu _3u}`$
$`\times {}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}i(\omega _2\mu _3),& i(\omega _3\mu _3)\\ 1i(\mu _1\mu _3),& 1i(\mu _2\mu _3)\end{array}|ie^u\right]`$ (7)
The $`{}_{2}{}^{}F_{2}^{}`$ can be defined in terms of the series expansions

$${}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}\alpha _1,& \alpha _2\\ \rho _1,& \rho _2\end{array}|z\right]=\underset{k=0}{\overset{\mathrm{}}{}}\frac{(\alpha _1)_k(\alpha _2)_k}{(\rho _1)_k(\rho _2)_k}\frac{z^k}{k!}$$
(8)
where $`(\alpha )_k`$ is a Pochhammer symbol, $`(\alpha )_k=\alpha (\alpha +1)\mathrm{}(\alpha +k1)`$. The solutions to Eq. (4) are thus

$$\psi _i=C_1\psi _i^{(1)}+C_2\psi _i^{(2)}+C_3\psi _i^{(3)},$$
(9)
where the constants $`C_j`$ are determined by the boundary conditions: $`\psi _1(r=0)=1`$, $`\psi _2(r=0)=0`$, $`\psi _3(r=0)=0`$.

For neutrino masses and energies of physical interest, the parameters and arguments of the $`{}_{2}{}^{}F_{2}^{}`$ become too large for the series expansion to be useful. A practical procedure is to adopt a stationary-phase approximation for the $`{}_{2}{}^{}F_{2}^{}`$ of $`\psi _i^{(3)}`$, and express the others by $`{}_{3}{}^{}F_{1}^{}`$ functions. This procedure is outlined below.

Consider the ordinary differential equation for $`{}_{2}{}^{}F_{2}^{}`$ in the form

$$\begin{array}{c}[(z\frac{d}{dz}+\beta _1)(z\frac{d}{dz}+\beta _2)(z\frac{d}{dz}+\beta _3)\hfill \\ \hfill z(z\frac{d}{dz}+\alpha _1)(z\frac{d}{dz}+\alpha _2)]f=0.\end{array}$$
(10)
This can be converted to an ordinary differential equation for $`{}_{3}{}^{}F_{1}^{}`$ by substituting $`\widehat{z}=z^1`$:

$$\begin{array}{c}[(\widehat{z}\frac{d}{d\widehat{z}}\alpha _1)(\widehat{z}\frac{d}{d\widehat{z}}\alpha _2)\hfill \\ \hfill +\widehat{z}\left(\widehat{z}\frac{d}{d\widehat{z}}\beta _1\right)\left(\widehat{z}\frac{d}{d\widehat{z}}\beta _2\right)\\ \hfill (\widehat{z}\frac{d}{d\widehat{z}}\beta _3)]f=0.\end{array}$$
(11)
Two solutions are of the form

$$z^{\mathrm{power}}{}_{3}{}^{}F_{1}^{}\left[\begin{array}{ccc}a_1,& a_2,& a_3\\ b& & \end{array}|z^1\right].$$
(12)
The full solutions (9) can then be constructed schematically (leaving out powers) as follows:

$$\psi _i=A_i{}_{2}{}^{}F_{2}^{}(1)+B_i{}_{3}{}^{}F_{1}^{}(1)+C_i{}_{3}{}^{}F_{1}^{}(2).$$
(13)
The series expansion for $`{}_{3}{}^{}F_{1}^{}`$ has zero radius of convergence. However, it can be expressed in terms of an integral involving the familiar hypergeometric function $`{}_{2}{}^{}F_{1}^{}`$,

$$\begin{array}{c}{}_{3}{}^{}F_{1}^{}\left[\begin{array}{ccc}a_1,& a_2,& a_3\\ b& & \end{array}|x^1\right]\hfill \\ \hfill =\frac{1}{\mathrm{\Gamma }(a_1)}x^{a_1}_0^{\mathrm{}}𝑑te^{xt}t^{a_11}\\ \hfill \times {}_{2}{}^{}F_{1}^{}(a_2,a_3;b;t)\end{array}$$
(14)
For large parameters and argument, the Pochhammer contour $`𝒫`$ is useful :

$$\begin{array}{c}{}_{2}{}^{}F_{1}^{}(a,b;c;t)\hfill \\ \hfill =\frac{\mathrm{\Gamma }(c)e^{i\pi c}}{4\mathrm{\Gamma }(b)\mathrm{\Gamma }(cb)\mathrm{sin}\pi b\mathrm{sin}\pi (cb)}\\ \hfill \times _𝒫s^{b1}(1s)^{cb1}(1+ts)^a𝑑s\end{array}$$
(15)
This approach leads to an accurate and efficient evaluation of the three neutrino wave functions in terms of stationary phase approximations to the $`{}_{2}{}^{}F_{2}^{}`$ and $`{}_{3}{}^{}F_{1}^{}`$ functions.

## 3 Linear density

The terminology “linear electron density” is used to mean that $`N_e(x)`$ is a linear function of $`x`$.

### 3.1 Two generations

The case of two states and a linear potential has been studied extensively, starting with Landau and Zener in the 1930s and applied to neutrino mixing in the 1980s . After a suitable scaling and shift of the variable, one has

$$i\frac{d}{dt}\left[\begin{array}{c}\psi _1(t)\\ \psi _2(t)\end{array}\right]=\left[\begin{array}{cc}t& a_2\\ a_2& 0\end{array}\right]\left[\begin{array}{c}\psi _1(t)\\ \psi _2(t)\end{array}\right],$$
(16)
or more explicitly:

$`i{\displaystyle \frac{d}{dt}}\psi _1(t)`$ $`=`$ $`t\psi _1(t)+a_2\psi _2(t),`$
$`i{\displaystyle \frac{d}{dt}}\psi _2(t)`$ $`=`$ $`a_2\psi _1(t).`$ (17)
Elimination of $`\psi _2`$ gives

$$\frac{d^2\psi _1(t)}{dt^2}it\frac{d\psi _1(t)}{dt}+(a_2^2i)\psi _1(t)=0.$$
(18)
The first-derivative term can be removed by taking

$$\psi _1(t)=e^{it^2/4}\varphi _1(t).$$
(19)
Then the equation for $`\varphi _1(t)`$ is

$$\frac{d^2\varphi _1(t)}{dt^2}+(\frac{1}{4}t^2+a_2^2\frac{1}{2}i)\varphi _1(t)=0.$$
(20)
Two linearly independent solutions of this equation are the parabolic cylinder functions

$$D_\rho (\pm e^{i\pi /4}t)$$
(21)
where $`\rho =ia_2^21`$.

Parabolic cylinder functions are special cases of the confluent hypergeometric function,

$$D_\rho (z)=2^{(\rho 1)/2}e^{z^2/4}z\mathrm{\Psi }(\frac{1}{2}\frac{1}{2}\rho ,\frac{3}{2};\frac{1}{2}z^2)$$
(22)
In terms of confluent hypergeometric functions $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$:

$`\psi _1(t)`$ $`=`$ $`t[C\mathrm{\Phi }(1+\frac{1}{2}ia_2^2,\frac{3}{2};\frac{1}{2}it^2)`$ (23)
$`+C^{}\mathrm{\Psi }(1+\frac{1}{2}ia_2^2,\frac{3}{2};\frac{1}{2}it^2)]`$
Unfortunately, it is not clear how to generalize this approach to $`N3`$.

Let us therefore consider an alternative solution to the $`N=2`$ case . We start by writing

$$F(\zeta )=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑te^{i\zeta t}\psi _1(t).$$
(24)
Then it follows from Eq. (18) that $`F(\zeta )`$ satisfies the first-order differential equation

$$\zeta ^2F(\zeta )\frac{d}{d\zeta }[i\zeta F(\zeta )]+(a_2^2i)F(\zeta )=0$$
(25)
or

$$\frac{1}{F(\zeta )}\frac{dF(\zeta )}{d\zeta }=\frac{i}{\zeta }(a_2^2\zeta ^2).$$
(26)
Integrating over $`\zeta `$, one finds

$$F(\zeta )=\text{const.}e^{i\zeta ^2/2}\zeta ^{ia_2^2}.$$
(27)
With the notation \[cf. Eq. (32)\]

$$b_1=\mathrm{}\text{and}b_{N+1}=+\mathrm{},$$
(28)
the two solutions can be written as

$$\psi _1^{(n)}(t)=_{b_n}^{b_{n+1}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}|\zeta |^{ia_2^2}$$
(29)
for $`n=1,2`$. It can be shown that they are confluent hypergeometric functions of the correct parameters and argument \[i.e., identical to Eq. (23)\].

### 3.2 General $`N`$

This second approach outlined above has the advantage that it can be generalized to an arbitrary number of neutrino flavors . We start out by writing the equation analogous to (1) in dimensionless standard form

$$i\frac{d}{dt}\psi (t)=A(t)\psi (x),$$
(30)
where

$$\psi (t)=\left[\begin{array}{c}\psi _1(t)\\ \psi _2(t)\\ \psi _3(t)\\ \mathrm{}\\ \psi _N(t)\end{array}\right]$$
(31)
and

$$A(t)=\left[\begin{array}{ccccc}t& a_2& a_3& \mathrm{}& a_N\\ a_2& b_2& 0& \mathrm{}& 0\\ a_3& 0& b_3& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& & \mathrm{}\\ a_N& 0& 0& \mathrm{}& b_N\end{array}\right].$$
(32)
There are two kinds of equations:

$$i\frac{d\psi _1(t)}{dt}=t\psi _1(t)+\underset{j=2}{\overset{N}{}}a_j\psi _j(t)$$
(33)
and, for $`k=2,3,4\mathrm{}N`$:

$$\left(i\frac{d}{dt}b_k\right)\psi _k(t)=a_k\psi _1(t).$$
(34)
One finds

$$\psi (t)=\underset{n=1}{\overset{N}{}}C_n\psi ^{(n)}(t),$$
(35)
where

$$\begin{array}{c}\psi ^{(n)}(t)=_{b_n}^{b_{n+1}}𝑑\zeta e^{i\zeta t}e^{i\zeta ^2/2}\hfill \\ \hfill \times \left(\underset{j=2}{\overset{N}{}}|\zeta b_j|^{ia_j^2}\right)\left[\begin{array}{c}1\\ \frac{a_2}{\zeta b_2}\\ \frac{a_3}{\zeta b_3}\\ \mathrm{}\\ \frac{a_{N1}}{\zeta b_{N1}}\\ \frac{a_N}{\zeta b_N}\end{array}\right]\end{array}$$
(36)
While these solutions are reasonably simple, a numerical Fourier transform is required. As $`t\pm \mathrm{}`$, explicit expressions can be written out for the different $`\psi ^{(n)}`$. However, these can not directly be used for imposing the boundary conditions, since they would correspond to negative density at $`t\mathrm{}`$.

## 4 The LNW mass matrix

For quark mixing, it was found  that a particular, simple texture for the $`d`$ ($`d`$, $`s`$, $`b`$) and $`u`$ ($`u`$, $`c`$, $`t`$) quark mass matrices leads to an acceptable CKM matrix . This same mass matrix has been applied to the case of three neutrinos , and rather good fits to the atmospheric  and solar \[??\] neutrino data have been obtained.

The mass matrix is assumed to have the form

$$M=\left(\begin{array}{ccc}0& d& 0\\ d& c& b\\ 0& b& a\end{array}\right)$$
(37)
with $`b^2=8c^2`$. The eigenvalues are given by $`m_1`$, $`m_2`$, and $`m_3`$, with $`m_1m_3`$.

In order to outline the diagonalization, whereby $`M=RM_{\mathrm{diag}}R^\mathrm{T}`$, let us introduce the notation

$`S_1`$ $``$ $`m_3m_2+m_1,`$
$`=`$ $`a+c`$
$`S_2`$ $``$ $`m_3m_2m_3m_1+m_2m_1,`$
$`=`$ $`8c^2+d^2ac`$
$`S_3`$ $``$ $`m_1m_2m_3`$ (38)
$`=`$ $`ad^2.`$
Then, a cubic equation for the parameter $`a`$ can be written as

$$9a^317S_1a^2+(8S_1^2+S_2)aS_3=0.$$
(39)
A physical solution requires $`a`$ real and positive. This is equivalent to having three real solutions for $`a`$. One of these is negative and two are positive. At any point inside the allowed domain in the $`m_1/m_3`$$`m_2/m_3`$ plane, there are thus two allowed solutions, denoted Solutions 1 and 2.

Let us consider first the atmospheric neutrino data. The Super-Kamiokande results  give $`\mathrm{\Delta }m^2(23)\times 10^3`$ eV, with $`\mathrm{sin}^2(2\theta )1`$. The survival of muon neutrinos is given by

$$\begin{array}{c}P_{\nu _\mu \nu _\mu }(t)=14[U_{\mu 1}^2U_{\mu 2}^2\mathrm{sin}^2\left(\frac{\mathrm{\Delta }m_{21}^2t}{4p}\right)\hfill \\ \hfill +U_{\mu 1}^2U_{\mu 3}^2\mathrm{sin}^2\left(\frac{\mathrm{\Delta }m_{31}^2t}{4p}\right)\\ \hfill +U_{\mu 2}^2U_{\mu 3}^2\mathrm{sin}^2\left(\frac{\mathrm{\Delta }m_{32}^2t}{4p}\right)],\end{array}$$
(40)
where $`U`$ is the neutrino mixing matrix. In the limit of $`\mathrm{\Delta }m_{21}^2t/4p1`$ this simplifies, and invoking further unitarity, one finds

$$P_{\nu _\mu \nu _\tau }(t)4U_{\mu 3}^2U_{\tau 3}^2\mathrm{sin}^2\left(\frac{\mathrm{\Delta }m_{32}^2t}{4p}\right),$$
(41)
which suggests that one needs $`|U_{\mu 3}U_{\tau 3}|=𝒪(1)`$. This can be achieved within the model (for both Solutions 1 and 2), for $`m_1m_3`$, with also $`m_2`$ small compared with $`m_3`$. Furthermore, the data suggest that the scale $`m_3`$ must be such that $`m_3^2(23)\times 10^3`$ eV.

Fits to atmospheric data confirm this qualitative analysis. Invoking also the solar Cl, Ga and Super-Kamiokande neutrino data \[??\], one finds that both Solutions 1 and 2 give good fits for $`m_1m_3`$, with $`m_2`$ also small as compared with $`m_3`$. Forming a $`\chi ^2`$ from these different atmospheric and solar survival probabilities, we found good fits , with $`m_3`$ of the order of 0.05 eV, $`m_2`$ about 0.01 eV, and $`m_1`$ 0.001–0.004 eV. In terms of the more conventional two-flavour analyses for the solar-neutrino sector, these fits roughly correspond to the large-mixing-angle solution.

## 5 Summary

We have reviewed analytic work on the solutions to the MSW equations for three neutrino flavours. Such results are very valuable for a fast scanning over the parameters of some given model for the mass matrix.

Also, we have more briefly reviewed the LNW mass matrix, as applied to the neutrino data. This is a very constrained model that in the quark sector describes the CKM matrix, and in the neutrino sector gives the mixing in terms of the mass eigenvalues.

The solar neutrino data has also been studied within the same model, using numerical integration methods (no $`{}_{2}{}^{}F_{2}^{}`$’s) . An additional fit was then found at $`m_12.8\times 10^6`$ eV, corresponding to the small-mixing-angle solutions. However, this point is disfavoured by the atmospheric neutrino data, and by the electron recoil spectrum.

Acknowledgments

This work is supported in part by the United States Department of Energy under Grant No. DE-FG02-84ER40158, and by the Research Council of Norway.