text:
# Two-band second moment model and an interatomic potential for caesium

## I Introduction

Semiempirical models for metallic binding have had a long and successful history in computer modelling. The most significant development came in the mid eighties with the implementation of ‘embedded atom’ potentials db  based loosely on density functional theory HKS  and ‘N-body’ potentials FS  based on the tight binding second moment approximation duc .

Although the rationale for these potentials suggests applicability to free-electron and transition metal systems respectively, the implied functional form is similar in each case. A large number of successful parameterisations of the functional forms suggested by these methods have been used, and they have established themselves as the standard method for modelling metallic systems.

In the second moment approximation to tight binding, the cohesive energy is proportional to the square root of the bandwidth, which can be approximated as a sum of pairwise potentials representing squared hopping integrals. Assuming atomic charge neutrality, this argument can be extended to all band occupancies and shapes avf . For simplicity, consider a rectangular $`d`$-band of full width W centred on the atomic energy level $`E_0`$. The energy for this band relative to the free atom (the bond energy) is given by:

$$U_{bond}=_{W/2}^{E_f=\left(\frac{n}{N}\frac{1}{2}\right)W}E\frac{N}{W}𝑑E=\frac{W}{2N}n\left(nN\right)$$
(1)
where $`n`$ is the occupation of the band and $`N`$ the capacity (for $`d`$-bands N=10, for $`s`$-bands N=2). The band width is defined as the square root of the second moment of the density of states. If we project the density of states onto each atom, and assume that each atom is charge neutral then the energy becomes:

$$U_{bond}=\underset{i}{}_{W_i/2}^{E_f=(\frac{n}{N}\frac{1}{2})W_i}E\frac{N}{W_i}𝑑E=\underset{i}{}\frac{W_i}{2N}n(Nn)$$
(2)
and the second moment of the density of states can be calculated as the sum of the squares of the hopping integrals to nearest neighbours. This latter operation can be written as a sum of pair potentials,

$$W_i=\sqrt{\underset{j}{}\varphi (r_{ij})}$$
(3)
The enormous success of Finnis-Sinclair (and embedded atom-type) potentials arises from their extreme computational efficiency: essentially they are no more computationally expensive than a conventional pair potential. This allows them to be applied to extremely large scale simulations, addressing complex geometries and phenomena which are intractability with more accurate quantum mechanical models.

The computational simplicity follows from the formal division of the energy into a sum of energies per atom, which can in turn be evaluated locally. Here we follow a similar philosophy, seeking to incorporate as much of the relevant physics as possible, without increasing computational cost.

The alkali and alkaline earth metals appear at first glance to be close packed metals, forming fcc, hcp or bcc structures at ambient pressures. However, compared with transition metals they are easily compressible, and at high pressures adopt more complex “open” structures (with smaller interatomic distances). The simple picture of the physics here is of a transfer of electrons from an $`s`$ to a $`d`$ band takemura4 ; abd ; macmahan , the $`d`$-band being more compact but higher in energy. Hence, at the price of increasing their energy, U, atoms can reduce their volume, V, and since the stable structure at 0K is determined by minimum enthalpy, H=U+PV at high pressure this transfer becomes energetically favourable. The net result is a metal-metal phase transformation characterised by a large reduction in volume (and often also in conductivity), the crystal structure itself is not the primary order parameter, and in some cases the transition may even be isostructural.

The mechanism of the phase transition is unknown. Although it is obvious that an isostructural phase transition is accompanied by an instability of the bulk modulus, the onset of this instability may occur after that in the shear modulus. Thus the mechanism may involve shearing rather than isostructural collapse, particularly if a continuous interface between the two phases exists, as in a shockwave shock .

## II The Two Band Model

### II.1 Energies

The Finnis Sinclair formalism has been successfully used to study a large variety of systems. However, we are interested in systems in which electrons change from one type of orbit to another and in particular, from an $`s`$-type orbital to a $`d`$-type orbital as the sample is pressurised.

Consider therefore two rectangular bands of widths $`W_1`$ and $`W_2`$ as shown in figure 1

with widths evaluated using equation 3. The bond energy of an atom may be written as the sum of the bond energies of the two bands on that atom as in equation 2, and a third term giving the energy of promotion from band 1 to band 2 (see equation 6):

$$U_{bond}=\underset{i}{}\frac{W_{i1}}{2N_1}n_{i1}(n_{i1}N_1)+\frac{W_{i2}}{2N_2}n_{i2}(n_{i2}N_2)+E_{prom}$$
(4)
where $`N_1`$ and $`N_2`$ are the capacities of the bands and $`n_{i1}`$ and $`n_{i2}`$ are the numbers of electrons in each band on the $`i^{\text{th}}`$ atom.

For an element, enforcing charge neutrality, the sum of the numbers of electrons is necessarily constrained to be equal to the total number of electrons on an atom,

$$n_{i1}+n_{i2}=T.$$
(5)
The difference between the energies of the band centres $`\alpha _1`$ and $`\alpha _2`$ is assumed to be fixed. The values of $`\alpha `$ correspond to the appropriate energy levels in the isolated atom. Thus, $`\alpha _2\alpha _1`$ is the excitation energy from one level to another. For alkali and alkaline earth metals, the free atom occupies only $`s`$-orbitals, so the promotion energy term is therefore simply

$$E_{prom}=n_2(\alpha _2\alpha _1)=n_2E_0,$$
(6)
where $`E_0=\alpha _2\alpha _1`$.

Thus the band energy can be written as a function of $`n_{i1}`$, $`n_{i2}`$ and the bandwidths (evaluated at each atom as a sum of pair potentials, within the second moment approximation). Defining:

$$\eta _i=n_{i1}n_{i2}$$
(7)
and using equation 5 we can write:

$$U_{bond}=\underset{i}{}\frac{\eta _i}{4}\left(W_{i1}W_{i2}\right)\frac{T}{4}\left(W_{i1}+W_{i2}\right)+\frac{\eta _i^2+T^2}{8}\left(\frac{W_{i1}}{N_1}+\frac{W_{i2}}{N_2}\right)+\frac{\eta _iT}{4}\left(\frac{W_{i1}}{N_1}\frac{W_{i2}}{N_2}\right)+\frac{T\eta _i}{2}E_0,$$
(8)
Although this expression looks unwieldy, it is in fact computationally efficient, requiring for its evaluation only two sums of pair potentials for $`W`$ (see equation 3) and a minimisation at each site independently with respect to $`\eta _i`$. Since they are local variables, it is possible to write closed form expressions for $`\eta _i`$ which minimise the total energy (see equation 14).

In addition to the bonding term, the Finnis-Sinclair form includes a pairwise repulsion between the ions, which is due primarily to the screened ionic charge and orthogonalisation of the valence electrons. In the present case, this pair potential should be a function of $`\eta _i`$. In keeping with maintaining locality of the energy we write the pairwise contribution to the energy in the intuitive form, as the sum of two terms, one from each “band”, proportional to the number of electrons in that band

$$V(r_{ij})=(n_{i1}+n_{j1})V_1(r_{ij})+(n_{i2}+n_{j2})V_2(r_{ij})$$
(9)
We rearrange this to give the energy as a sum over atoms

$$U_{pair}=\underset{i}{}\left[n_{i1}\underset{ji}{}V_1(r_{ij})+n_{i2}\underset{ji}{}V_2(r_{ij})\right]$$
(10)
The total energy is now simply

$$U_{tot}=U_{pair}+U_{bond}.$$
(11)
This depends on $`\eta _i`$, but as described above the $`\eta _i`$ take the values which minimise the energy and one can solve:

$$\frac{U_{tot}}{\eta _i}=0$$
(12)
explicitly for $`\eta _{i_0}`$ independently at each atom, whence

$$\eta _{i_0}=\frac{N_1N_2}{W_{i1}N_2+W_{i2}N_1}\left[W_{i1}W_{i2}T\left(\frac{W_{i1}}{N_1}\frac{W_{i2}}{N_2}\right)+2E_02U_{i1,pair}+2U_{i2,pair}\right].$$
(13)
Depending on the number of electrons in the system it may not be possible to realise this. Charge neutrality requires that $`|\eta _i|`$ cannot be greater than the total number of electrons $`T`$ per atom. The fixed capacities of the bands ($`N_1`$ and $`N_2`$) can also prohibit the realisation of $`\eta _{i_0}`$. It is therefore necessary to limit the values which $`\eta _i`$ may have:

$$\eta _i=\{\begin{array}{cc}\mathrm{min}(T,2N_1T),\hfill & \text{if }\eta _{i_0}>\mathrm{min}(T,2N_1T)\hfill \\ \mathrm{max}(T,T2N_2),\hfill & \text{if }\eta _{i_0}<\mathrm{max}(T,T2N_2)\hfill \\ \eta _{i_0},\hfill & \text{otherwise}\hfill \end{array}$$
(14)
where $`\eta _{i_0}`$ is given by equation 13. The expressions for $`\eta _i`$ involves only constants and sums of pair potentials, and can be evaluated independently at each atom at a similar computational cost to a standard many-body type potential.

### II.2 Forces

Since the energy is variational in $`\eta _i`$ within this model, we have

$$\frac{U_{tot}}{\eta _i}=0$$
(15)
and the force on the $`i^{\mathrm{th}}`$ atom is

$`𝐟_i`$ $`=`$ $`{\displaystyle \frac{dU_{tot}}{d𝐫_i}}`$ (16)
$`=`$ $`{\displaystyle \frac{U_{tot}}{𝐫_i}}|_\eta {\displaystyle \frac{U_{tot}}{\eta }}{\displaystyle \frac{\eta }{d𝐫_i}}`$
$`=`$ $`{\displaystyle \frac{U_{tot}}{𝐫_i}}|_\eta `$
Hence the force is simply the derivative of the energy at fixed $`\eta `$. Basically, this is the Hellmann-Feynman theorem p:Feynman  which arises here because $`\eta `$ is essentially a single parameter representation of the electronic structure.

This result means that, like the energy, the force can be evaluated by summing pairwise potentials. Hence the model is well suited for large scale molecular dynamics.

The force derivation is somewhat tedious (see appendix B), the result being:

$`𝐟_𝐢={\displaystyle \underset{j}{}}`$ $`[(\stackrel{~}{W_{i1}}+\stackrel{~}{W_{j1}})\varphi _1^{}(r_{ij})+(\stackrel{~}{W_{i2}}+\stackrel{~}{W_{j2}})\varphi _2^{}(r_{ij})`$
$`(n_{i1}+n_{j1})𝐕_{}^{}{}_{1}{}^{}(r_{ij})(n_{i2}+n_{j2})𝐕_{}^{}{}_{2}{}^{}(r_{ij})]\widehat{𝐫}_{ij}`$
with

$$\stackrel{~}{W_{ib}}=\frac{n_{ib}(N_bn_{ib})}{4N_bW_{ib}}$$

and where the numbers of electrons in each band, $`n_{ib}`$, have been calculated analytically using equations 13, 14 and 7. The contribution to the force from each neighbour acts along the direction of the vector between the atoms.

## III A Model Potential for Caesium

One application of the model is to simulate the pressure-induced isostructural phase transition in materials such as caesium. Here the transformation arises from electronic transition from the $`s`$ to the $`d`$-band. Caesium adopts the bcc structure (Cs I) at ambient pressure, transforming under slight additional pressure to fcc (Cs II) and under further increase transforming isostructurally to fcc (Cs III) p:HallMerrillBarnett  111There is some recent evidence that the Cs III phase may be a more complex structure p:MalcomCsIII .. Here we seek to represent the volume collapse Cs II $``$ Cs III transition.

### III.1 Parametrisation

To make a usable potential, the functional forms of $`\varphi `$ and $`V`$ must be chosen. Although this is somewhat arbitrary, the physical picture of hopping integral and screened ion-ion potential suggests that both should be short ranged, continuous and reasonably smooth.

In the present case, we are interested in the gross features of the model, so we adopt simple forms fitted to the cohesive energies and volume of the high and low pressure polytypes of caesium. For more complex applications, it may be necessary to fit other properties.

Following Finnis and Sinclair, we choose for the hopping integral

$$\varphi _b(r_{ij})=\{\begin{array}{cc}C_b(d_br_{ij})^3,\hfill & \text{if }r_{ij}d_b\hfill \\ 0,\hfill & \text{if }r_{ij}>d_b\hfill \end{array}$$
(18)
and for the pairwise part, following Lennard-Jones

$$V_b(r_{ij})=\underset{j}{}\frac{A_b}{r_{ij}^{12}}$$
(19)
where $`b`$ labels the band.

We select the promotion energy $`E_0`$ to be that required to promote an electron from the $`6s`$ level into the $`5d`$ level of an isolated atom b:CRCHandbook ; calc:Simon . The band capacities are $`N_s=2`$, $`N_d=10`$ and the total number of electrons per atom is $`T=1`$. The cutoff radii $`d_s`$ and $`d_d`$ are chosen to be between the second and third nearest neighbours and between the first and second nearest neighbours respectively. The remaining four parameters ($`A_s`$, $`A_d`$, $`C_s`$, $`C_d`$) are available for fitting.

*Ab initio* calculations using pseudopotentials, plane waves and the generalized gradient approximation xie ; SKR:thesis  show that the energy-volume relations for bcc and fcc caesium are almost degenerate, so we have fitted the remaining parameters to the atomic volumes of Cs I (115.9 Å3/atom), Cs II (67.5 Å3/atom) and Cs III (48.7 Å3/atom) b:CRCHandbook  as well as the transition pressure between phases II and III (4.3 GPa) b:CRCHandbook  and the cohesive energy of phase I (-0.704 eV). The cohesive energy is the sum of the heat of formation and the heat of vaporisation b:CRCHandbook .

Initial estimates for the parameters were determined using a symbolic mathematics package. The final parameters were arrived at by an iterative process: Least squares, conjugate gradients and *ab oculo* minimisation techniques were used to determine the best parameters for particular cutoff radii. The cutoff radii were then adjusted by hand to improve the fit. This process was repeated until the optimum fit was achieved. The final parameters are given in table 1.

Figure 2 shows the energy-volume curves for the fcc and bcc structures calculated using the model. The fcc structure is stable everywhere compared to the bcc structure. Experimentally however, Cs I has bcc structure with an equilibrium volume of 116 Å3. The present potential does not have a low pressure bcc phase. However, at ambient pressure the bcc and fcc curves are almost degenerate (0.02 eV/atom difference) in agreement with the *ab initio* calculations. The predicted equilibrium volume of the ambient pressure phase agrees well with the experimental equilibrium volume of phase I of caesium. Although the predicted volumes of phases II and III are larger than the experimental values, the predicted transition pressure and volume collapse are in good agreement with experiment.

### III.2 Elasticity

The elastic moduli are not fitted explicitly, so their behaviour represents a sensitive test of the model. Application of the pressure generalisation of the Born stability criteria is rather confused in the present literature yip ; bbk ; marcus  so we lay this out in appendix A.

Since caesium exhibits an isostructural phase transition the bulk modulus must formally become negative at some volume (where the structure is unstable). It is less clear whether the other Born stability criteria will be violated: ab initio density functional perturbation theory calculations xie  find an instability in the long range acoustic phonons: at zero pressure this is equivalent to an instability in $`C^{}=(C_{11}C_{12})/2`$.

With the present potential the analytic expressions for the elastic constants are complicated: there is no simplification akin to equation 15 for the second derivatives. Consequently, we evaluate the elastic constants numerically from finite strain calculations, see figure 3. We find that the volume collapse is announced by a slight softening of the bulk modulus which then goes negative in the unstable region. Although the shear and tetragonal shear decrease in the unstable region, neither actually goes negative.

### III.3 Defects

One of the major features of many-body and embedded atom type potentials is their ability to describe defects such as surfaces, vacancies and interstitials without the constraints of pair potentials, for which the undercoordination defect energy is simply the sum of broken bond energies, less a small amount from relaxation of atomic positions. Thus the vacancy formation energy is typically the same as the cohesive energy, whereas in real materials it is typically less than half. The present potential is not fitted to any defect configuration, so it is of interest to see what it predicts - moreover the variety of environments associated with defects provides a good check against pathologies.

The relaxed surface energies for four low-energy surfaces are given in table 2. As is usual with many-body potentials the energy is much lower than would be expected from simple bond counting. For the (111) surface the energy for each atom in the plane is increased by about $`1/9`$ of its cohesive energy, while each atom has 1/4 of its bonds broken. There is also transfer of all electrons from $`d`$-like to $`s`$-like states at the free surface. This also leads to an unusual outwards relaxation of the surface atoms.

Likewise the vacancy formation energy and volume (table 3) are rather typical of many-body potential results.

The interstitial formation energy is especially low (table 3), around 1.8 eV depending on the orientation and the amount of relaxation allowed locally, meaning that thermal interstitial formation is possible. This low value is due to the transfer of electrons from $`s`$\- to $`d`$-band on the interstitial atom or dumbbell, similar to the high pressure behaviour. Consistent with recent *ab initio* calculations in various elements madden ; becquart ; SWHan , the calculated interatomic spacing of the dumbbell atoms is much smaller than in the bulk, (around 15%) and much smaller than is typical of standard EAM-type potentials. As a consequence of this, the associated strain fields are considerably smaller.

For a detailed study of point defects in caesium, it would be appropriate to reparametrise the potential with point defects included in the fitting, but the good results obtained here without such fitting suggest that the present model contains the right physics.

## IV Extrapolation to transition metals

In principle, the current formalism should be applicable to $`d`$-band metals. We do not intend to refit the potentials here, but by applying the parameters fitted for caesium with appropriate scaling for ionic charge, and simply varying the total number of electrons we recover the parabolic behaviour of the cohesive energy and bulk modulus which characterises the transition metal series. Considering that there are no fitting parameters, the agreement with experimental results is extraordinarily good. Moreover, the volume collapse phase transition exists only for N=1 and N=2 (consistent with experiment). However, since no information about band shape is included, the sequence of crystal structures cannot be reproduced, so we consider here only the fcc structures, calculating their “experimental” values from the experimental density.

While the extrapolated potentials do not represent the optimal parametrisation for specific transition metals, the recovery of the trends across the group lends weight to the idea that the two-band model correctly reproduces the physics of this series.

## V Conclusions

We have presented a model to describe $`sd`$ transfer within the framework of the empirical second-moment tight binding model. With a very simple parametrisation, our model describes the isostructural phase transition and associated $`sd`$ transfer in Cs and allows study of the elastic instabilities which occur under pressure and uniaxial stress. Throughout, our formalism is guided by the principle of retaining as much electronic detail as possible within the constraint of the computational complexity of short-ranged pair potentials

Our potential form involves evaluation of sums of pair potentials, and minimisation of the energy with respect to band occupation. In the approximation of local density of states evaluated at each atom, it becomes possible to formally write the energy of each atom in terms of the local band occupation, a single parameter $`\eta _i`$. This variational parameter can be explicitly eliminated from the energy and force expressions, meaning that energies and forces can be evaluated by summing pair potentials.

To our knowledge, the present potential represents the most sophisticated energy function which has first derivatives which are analytic pair potentials, and the first exploitation of the Hellman-Feynman theorem in empirical potentials to achieve this. It is thus uniquely well suited for exploitation in standard molecular dynamics codes.

Extension of the formalism to transition metals is straightforward, and the caesium parametrisation gives a surprisingly good description across the group. While providing impressive proof of concept, this may not be of enormous practical use, however, since once there is a significant partial occupation of the $`d`$-band under all conditions the model is similar to the standard many-body potentials which are known to do a good job in describing transition metals.

The precise stable crystal symmetry depends on the shape of the band structure, and is not therefore correctly described in this rectangular band model. Many-body potentials never contain the correct physics to describe phase transitions, except in the martensitic case of freezing in soft phonons Pinsook ; Pinsook2 . However, the crystal structure which has the lowest energy can be determined by judicious choice of functions $`V`$ and $`\varphi `$. Here the simple use of power law repulsion and cubic attraction leads to close packed structures.

The two-band formalism can be easily combined with alternate parameterisations to get effective band shapes, and could be applied to ferromagnetic systems where the two bands represent separate spins, and the atomic $`E_0`$ term is replaced by a term which favours maximum spin.

## Appendix A Elastic instabilities under pressure - energy definition

There is a degree of confusion in the literature regarding the definition of elastic constants under pressure. To some extent this arises because of various definitions of strain (Lagrangian, Eulerian, volume conserving). Here we lay out the definitions in terms of energy. One important point to note is that under pressure the moduli which correspond to long wavelength phonons are volume conserving, while those corresponding to crystal stability are not. At finite pressure, this means they are different.

At pressure, the important quantity for stability in energetic terms is the Gibbs free energy $`G`$.

$$G=UTS+PV$$
(20)
where P is an externally applied hydrostatic pressure. $`UTS`$ is the Helmholtz free energy F. The change in the Gibbs energy due to a distortion of the crystal is

$$\mathrm{\Delta }G=\mathrm{\Delta }F+\mathrm{\Delta }P.V+P\mathrm{\Delta }V.$$
(21)
If the hydrostatic pressure is applied by an external mechanism then $`\mathrm{\Delta }P`$ is zero. Therefore

$$\mathrm{\Delta }G=\mathrm{\Delta }F+P\mathrm{\Delta }V.$$
(22)
Any arbitrary change in the unit cell be expressed in terms of the strain tensor $`\stackrel{}{\epsilon }`$. Using Voigt notation, the matrix representation of the strain tensor is written as

$$\stackrel{}{\epsilon }=\left(\begin{array}{ccc}e_1& e_6/2& e_5/2\\ e_6/2& e_2& e_4/2\\ e_5/2& e_4/2& e_3\end{array}\right).$$
(23)
If we write the equilibrium cell as a matrix comprising the three lattice vectors $`\underset{0}{\overset{}{V}}=(𝐫_1𝐫_2𝐫_3)`$, then any arbitrary strain give a new unit cell:

$$\stackrel{}{V}=\underset{0}{\overset{}{V}}(1+\stackrel{}{\epsilon }).$$
(24)
If the volume of the original cell is $`V_0=|\stackrel{}{V_0}|`$, then the volume of the new cell is

$`V/V_O`$ $`=`$ $`|\stackrel{}{V}|/|\underset{0}{\overset{}{V}}|`$
$`=`$ $`1+e_1+e_2+e_3+e_1e_2+e_2e_3+e_3e_1`$
$`e_4^2/4e_5^2/4e_6^2/4`$
$`+e_1e_2e_3e_1e_4^2/4e_2e_5^2/4e_3e_6^2/4+e_4e_5e_6/4.`$
Which can be expressed concisely in standard strain notation ($`i,j,k,l`$ representing cartesian directions) as:

$$\frac{\mathrm{\Delta }V}{V}=e_{ii}+\frac{1}{4}\left(2\delta _{ij}\delta _{kl}\delta _{ik}\delta _{jl}\delta _{il}\delta _{jk}\right)e_{ij}e_{kl}+𝒪(\epsilon ^3)$$
(26)
where $`\delta _{ij}`$ is the Kroneker delta and the usual implicit sum convention for tensors applies. The change in the Gibbs free energy may then be written

$$\mathrm{\Delta }G=\mathrm{\Delta }F+Pe_{ii}+\frac{PV}{4}\left(2\delta _{ij}\delta _{kl}\delta _{ik}\delta _{jl}\delta _{il}\delta _{jk}\right)e_{ij}e_{kl}.$$
(27)
We define $`C_{ij}`$ as the second order derivatives of the Helmholtz free energy with respect to the Voigt strains. and $`B_{ij}`$ as the second order derivatives of the Gibbs free energy with respect to the Voigt strains. The limiting case of long wavelength phonons is related to $`C_{ij}`$ which can be calculated from the dynamical matrix sjc , the crystal stability criteria are related to the $`B_{ij}`$ bbk . In the standard notation:

$$B_{ijkl}=C_{ijkl}+\frac{P}{2}\left(2\delta _{ij}\delta _{kl}\delta _{ik}\delta _{jl}\delta _{il}\delta _{jk}\right).$$
(28)
In Voigt notation:

$$B_{ij}=\{\begin{array}{cc}C_{ii}P/2\text{if }i=j\text{ and }i>3\hfill & \\ C_{ij}+P\text{if }ij\text{ and }i,j<4\hfill & \\ C_{ij}\text{otherwise}\hfill & \end{array}$$
(29)
The Born elastic stability criteria restrict what values various moduli may have if the crystal is to be stable b:Born . In the case of external pressure the relevant free energy is $`G`$ and the relevant moduli are $`B_{ij}`$. Due to the Voigt symmetry, the elastic constant tensor has in general, only 21 independent components. In cubic crystals however this number is further reduced by symmetry to 3, namely $`B_{11}`$, $`B_{12}`$ and $`B_{44}`$. Consequently, there are three common stability criteria which impose lower bounds on the bulk, shear and tetragonal shear moduli of cubic crystals. These may be written in Voigt notation:

$$\frac{1}{3}(B_{11}+2B_{12})>0,B_{44}>0\mathrm{and}\frac{1}{2}(B_{11}B_{12})>0.$$
(30)
Whence the generalised Born stability criteria for the elastic constants of cubic crystals at pressure become:

$`{\displaystyle \frac{1}{3}}(C_{11}+2C_{12}+2P)`$ $`>`$ $`0,`$
$`C_{44}{\displaystyle \frac{p}{2}}`$ $`>`$ $`0,`$
$`{\displaystyle \frac{1}{2}}(C_{11}C_{12}P)`$ $`>`$ $`0,`$ (31)
where

$$B_{11}=C_{11},B_{12}=C_{12}+P\mathrm{and}B_{44}=C_{44}\frac{p}{2}.$$
(32)
### Calculation of bulk, shear and tetragonal shear moduli from energies

In this paper we evaluate the elastic moduli by applying finite strain and measuring the change in total energy $`U_{tot}`$. The finite strain method has the advantage over analytic differentiation of automatically compensating for non-isotropic movement of the atoms. At zero temperature $`U_{tot}`$ is the same as the classical Helmholtz free energy $`F`$. At finite pressure, however, the system is not at a minimum of $`U_{tot}`$ with respect to strain.

When applying strain to a system under pressure, the first order term in the Gibbs free energy vanishes:

$$\mathrm{\Delta }G=e_{ij}\frac{G}{e_{ij}}+\frac{1}{2}e_{ij}e_{kl}\frac{^2G}{e_{ij}e_{kl}}=\frac{1}{2}e_{ij}e_{kl}\frac{^2G}{e_{ij}e_{kl}}$$
(33)
But the first order term for Helmholtz does not. For a cubic crystal, one can find $`(C_{11}+2C_{12})/3`$, $`(C_{11}C_{12})/2`$ and $`C_{44}`$ directly at any volume by simple distortions and measurement of $`F`$.

For bulk modulus we apply a pair of Voigt strains, with $`e_1=e_2=e_3=e`$, the other elements being zero. Assuming we make a pair of small distortions about a volume V, the difference

$$G(e)+G(e)2G(0)=V\frac{\left(B_{11}+2B_{12}\right)}{3}9e^2+𝒪(e^3).$$
(34)
equivalently, using the Helmholtz free energy

$$\frac{\left(C_{11}+2C_{12}\right)}{3}=\frac{(F(e)+F(e)2F(0))}{9Ve^2}+𝒪(e^3).$$
(35)
For $`(C_{11}C_{12})/2`$, the relevant distortion, volume conserving to first order, is $`e=e_1=e_2`$ and

$$\frac{\left(C_{11}C_{12}\right)}{2}=\frac{(F(e)+F(e)2F(0))}{4Ve^2}+𝒪(e^3).$$
(36)
For $`(C_{44}`$, the relevant distortion is $`e=e_6`$ and

$$C_{44}=\frac{(F(e)+F(e)2F(0))}{Ve^2}+𝒪(e^3).$$
(37)
## Appendix B Derivation of the force on an atom

Since the cohesive energy is variational with respect to the electron distribution $`\eta `$, the derivative of the energy with respect to the electron distribution is zero and so the force may be written as a partial derivative akin to the Hellmann-Feynman theorem p:Feynman . From equation 16

$$𝐟_i=\frac{U_{tot}}{𝐫_i}|_\eta =\underset{j}{}\frac{U_j}{𝐫_i}|_{\eta _j}.$$
(38)
where the sum over $`j`$ includes $`j=i`$. From equations 4, 6 and 8, the cohesive energy of the $`j^{th}`$ atom in terms of the number of electrons in each band ($`n_{j1}`$ and $`n_{j2}`$) is

$$U_j=\left[\frac{n_{j1}}{2}\left(\frac{n_{j1}}{N_1}1\right)W_{j1}+\frac{n_{j2}}{2}\left(\frac{n_{j2}}{N_2}1\right)W_{j2}+n_{j1}\underset{kj}{}V_1(r_{jk})+n_{j2}\underset{kj}{}V_2(r_{jk})+n_{j2}E_0\right],$$
(39)
and the force on atom $`i`$ is thus

$$𝐟_i=\underset{j}{}\left[\frac{n_{j1}}{2}\left(\frac{n_{j1}}{N_1}1\right)\frac{W_{j1}}{𝐫_i}+\frac{n_{j2}}{2}\left(\frac{n_{j2}}{N_2}1\right)\frac{W_{j2}}{𝐫_i}+\left(n_{j1}+n_{i1}\frac{}{𝐫_i}V_1(r_{ij})\right)+\left(n_{j2}+n_{i2}\frac{}{𝐫_i}V_2(r_{ij})\right)\right].$$
(40)
The width of band b on atom $`i`$ is $`W_{jb}=\sqrt{_{kj}\varphi _b(r_{jk})}`$ and its derivative is

$$\frac{W_{jb}}{𝐫_i}=\frac{1}{2W_{jb}}\underset{kj}{}\frac{}{𝐫_i}\varphi _b(r_{jk})=\{\begin{array}{cc}\frac{1}{2W_{jb}}_{kj}\frac{}{𝐫_i}\varphi _b(r_{jk}),\hfill & \mathrm{if}\mathrm{j}=\mathrm{i}\hfill \\ \frac{1}{2W_{jb}}\frac{}{𝐫_i}\varphi _b(r_{jk}),\hfill & \mathrm{if}\mathrm{j}\mathrm{i}.\hfill \end{array}$$
(41)
The pairwise terms are given in Eq.9, the derivatives are straightforward, noting that an atom does not interact with itself so $`V(r_{ii})=0`$

On combining the above equations, the force may then be written:

$`𝐟_i=`$ $`{\displaystyle \underset{ji}{}}[{\displaystyle \frac{n_{j1}}{4W_{j1}}}({\displaystyle \frac{n_{j1}}{N_1}}1){\displaystyle \frac{\varphi _1(r_{ji})}{𝐫_i}}+{\displaystyle \frac{n_{j2}}{4W_{j2}}}({\displaystyle \frac{n_{j2}}{N_2}}1){\displaystyle \frac{\varphi _2(r_{ji})}{𝐫_i}}`$ (42)
$`+n_{j1}{\displaystyle \frac{V_1(r_{ji})}{𝐫_i}}+n_{j2}{\displaystyle \frac{V_2(r_{ji})}{𝐫_i}}]`$
$`[{\displaystyle \frac{n_{i1}}{4W_{i1}}}({\displaystyle \frac{n_{i1}}{N_1}}1){\displaystyle \underset{ki}{}}{\displaystyle \frac{\varphi _1(r_{ik})}{𝐫_i}}+{\displaystyle \frac{n_{j2}}{4W_{j2}}}({\displaystyle \frac{n_{j2}}{N_2}}1){\displaystyle \underset{ki}{}}{\displaystyle \frac{\varphi _2(r_{ik})}{𝐫_i}}`$
$`+n_{j1}{\displaystyle \underset{ki}{}}{\displaystyle \frac{V_1(r_{ik})}{𝐫_i}}+n_{j2}{\displaystyle \underset{ki}{}}{\displaystyle \frac{V_2(r_{ik})}{𝐫_i}}].`$
To simplify the appearance of this equation we define

$$\stackrel{~}{W}_{ib}=\frac{n_{ib}(N_bn_{ib})}{4N_bW_{ib}}$$
(43)
and the force becomes

$`𝐟_i=`$ $`{\displaystyle \underset{ji}{}}\left[\stackrel{~}{W}_{j1}{\displaystyle \frac{\varphi _1(r_{ji})}{𝐫_i}}+\stackrel{~}{W}_{j2}{\displaystyle \frac{\varphi _2(r_{ji})}{𝐫_i}}n_{j1}{\displaystyle \frac{V_1(r_{ji})}{𝐫_i}}n_{j2}{\displaystyle \frac{V_2(r_{ji})}{𝐫_i}}\right]`$ (44)
$`+{\displaystyle \underset{ki}{}}\left[\stackrel{~}{W}_{i1}{\displaystyle \frac{\varphi _1(r_{ik})}{𝐫_i}}+\stackrel{~}{W}_{i2}{\displaystyle \frac{\varphi _2(r_{ik})}{𝐫_i}}n_{j1}{\displaystyle \frac{V_1(r_{ik})}{𝐫_i}}n_{j2}{\displaystyle \frac{V_2(r_{ik})}{𝐫_i}}\right].`$
The summations in the second set of square brackets may be reindexed in terms of $`j`$’s and included in the first set of square brackets to give

$`𝐟_i={\displaystyle \underset{ji}{}}[`$ $`\stackrel{~}{W}_{j1}{\displaystyle \frac{\varphi _1(r_{ji})}{𝐫_i}}+\stackrel{~}{W}_{j2}{\displaystyle \frac{\varphi _2(r_{ji})}{𝐫_i}}n_{j1}{\displaystyle \frac{V_1(r_{ji})}{𝐫_i}}n_{j2}{\displaystyle \frac{V_2(r_{ji})}{𝐫_i}}`$ (45)
$`+\stackrel{~}{W}_{i1}{\displaystyle \frac{\varphi _1(r_{ij})}{𝐫_i}}+\stackrel{~}{W}_{i2}{\displaystyle \frac{\varphi _2(r_{ij})}{𝐫_i}}n_{j1}{\displaystyle \frac{V_1(r_{ij})}{𝐫_i}}n_{j2}{\displaystyle \frac{V_2(r_{ij})}{𝐫_i}}].`$
Since both $`\varphi `$ and $`V`$ are pair potentials $`\varphi _b(r_{ij})=\varphi _b(r_{ji})`$ and $`V_b(r_{ij})=V_b(r_{ji})`$. The derivatives of $`\varphi `$ are

$$\frac{}{𝐫_i}\varphi _b(r_{ij})=\frac{}{r_i}\varphi _b(r_{ij})\widehat{𝐫}_{ij}=\varphi _b^{}(r_{ij})\widehat{𝐫}_{ij}$$

where the $`\varphi _b^{}`$ are scalar quantities and $`\varphi _b^{}(r_{ij})=\varphi _b^{}(r_{ji})`$. Similarly, the derivatives of $`V`$ are

$$\frac{}{𝐫_i}V_b(r_{ij})=\frac{}{r_i}V_b(r_{ij})\widehat{𝐫}_{ij}=V_b^{}(r_{ij})\widehat{𝐫}_{ij}$$

where again the $`V_b^{}`$ are scalar quantities.

Finally, the force on atom i is then
$`𝐟_i={\displaystyle \underset{ji}{}}[`$ $`(\stackrel{~}{W}_{i1}+\stackrel{~}{W}_{j1})\varphi _1^{}(r_{ij})+(\stackrel{~}{W}_{i2}+\stackrel{~}{W}_{j2})\varphi _2^{}(r_{ij})`$  (46)    $`(n_{i1}+n_{j1})V_1^{}(r_{ji})(n_{i2}+n_{j2})V_2^{}(r_{ji})]\widehat{𝐫}_{ij}`$


format: text

text:
# Anisotropy parameters of superconducting MgB2

## 1 Introduction

Physics and applications of new superconducting materials cannot be properly understood and developed without careful characterization of possible anisotropies. High-$`T_c`$ materials are one of the best examples of relevance of the superconducting anisotropy to all aspects of physics and possible applications.

It is a common practice to characterize anisotropic superconductors by a single anisotropy parameter defined as $`\gamma =\xi _a/\xi _c\lambda _c/\lambda _a`$. Here, $`\xi `$ is the coherence length, $`\lambda `$ is the penetration depth, and $`a,c`$ are principal directions of a uniaxial crystal of the interest here. Often $`\xi _a=\xi _b`$ are denoted in literature as $`\xi _{ab}`$; for formal reasons we prefer a single subscript notation. Since the upper critical fields, $`H_{c2,c}=\varphi _0/2\pi \xi _a^2`$ and $`H_{c2,a}=\varphi _0/2\pi \xi _a\xi _c`$, the anisotropy parameter can also be written as $`\gamma =H_{c2,a}/H_{c2,c}`$. Historically, the practice to assign a single anisotropy parameter to each material emerged after the anisotropic Ginzburg-Landau (GL) equations were proposed phenomenologically by Ginzburg  and derived microscopically for arbitrary gap and Fermi surface anisotropies by Gor’kov and Melik-Barkhudarov . Formally, this came out because near the critical temperature $`T_c`$, the same “mass tensor” $`m_{ik}`$ determines the anisotropy of both $`\xi `$ (of the upper critical fields $`H_{c2}`$) and of $`\lambda `$:

$$(\xi ^2)_{ik}m_{ik}^1,(\lambda ^2)_{ik}m_{ik}.$$
(1)
Therefore,

$$\gamma _\lambda ^2=\frac{\lambda _c^2}{\lambda _a^2}=\frac{m_{cc}}{m_{aa}}=\frac{\xi _a^2}{\xi _c^2}=\gamma _\xi ^2.$$
(2)
At arbitrary temperatures, however, the theoretical approach for calculating $`H_{c2}`$ (the position of the second order phase transition in high fields) has little in common with evaluation of $`\lambda `$ (the weak-field relation between the current and the vector potential), so that the anisotropies of these quantities are not necessarily the same. In fact, in materials with anisotropic Fermi surfaces and anisotropic gaps, not only the parameter $`\gamma _H=H_{c2,a}/H_{c2,c}`$ may strongly depend on $`T`$, but this ratio might differ considerably from $`\gamma _\lambda =\lambda _c/\lambda _a`$ at low temperatures .

In this brief review we collect data on anisotropy of MgB2 available to us, discuss briefly the methods used and the results obtained. Then we outline the weak-coupling microscopic approach which can be used to evaluate the anisotropy parameters and their $`T`$ dependence. We conclude with discussion of macroscopic consequences of different $`H_{c2}`$ and $`\lambda `$ anisotropies focussing on the torque problem, the quantity commonly measured to extract the anisotropy parameters.

## 2 Data review

Experimentally, superconducting anisotropy of MgB2 was under scrutiny right away after the discovery of this material. Almost all experimental studies were concerned with the $`H_{c2}`$ anisotropy, $`\gamma _H=H_{c2,ab}/H_{c2,c}`$.

Based on the AC susceptibility and the field dependent magnetization measurements on sorted powders, magnesium diboride was claimed to be an isotropic superconductor . The $`H_{c2}`$ anisotropy for separate particles settled on a flat surface was reported as $`\gamma _H1.73`$ . For a hot-pressed bulk samples values of $`\gamma _H`$ up to 1.2 were obtained from susceptibility and resistance measurements . Magnetotransport measurements on MgB2 thin films with different degrees of $`c`$-axis orientation yielded $`\gamma _H=1.82`$ ($`T_c`$ = 31-37 K) with higher anisotropy for films of higher resistivity . The temperature-independent $`\gamma _H`$ = 1.2 ($`T_c`$ = 39 K) was reported for $`c`$-oriented films, Ref. . Similar results ($`\gamma _H`$ = 1.2 - 1.8) were later obtained for a set of three films pulsed laser deposited (PLD) on different substrates ($`T_c`$ = 31-37 K) . In this case no apparent correlation was observed between the $`H_{c2}`$ anisotropy and residual resistance ratio (RRR) or superconducting transition temperature of the films. Quite different values, $`\gamma _H`$ = 9 - 13, were reported for in-situ grown PLD films with $`T_c`$ suppressed down to 24 - 27 K .

In a different approach, the $`H_{c2}`$ anisotropy of magnesium diboride was evaluated from a number of measurements on randomly oriented powders. Analysis of the conduction electron spin resonance (CESR) data on fine MgB2 powders taken in a wide temperature range yielded an estimate of $`\gamma _H`$ 8 . Similar idea (deconvoluting the signal in a certain temperature range in two components, one corresponding to superconducting state and the other originating from the normal state) was later used for interpretation of the 11B nuclear magnetic resonance (NMR) data in the mixed state of MgB2 . These measurements resulted in $`\gamma _H`$ 6. Subsequent 11B NMR and magnetization studies on Mg1-xAlxB2 ($`x`$ 0.025) reveal a considerable decrease of the anisotropy with Al doping .

In addition to local NMR and CESR probes, global measurements on samples with randomly oriented grains were used for the anisotropy studies. High field reversible DC magnetization data were well described assuming that MgB2 is an anisotropic superconductor with $`\gamma _H`$ = 6 - 9 and using the available experimental estimates of the average penetration depth $`\lambda (0)=(\lambda _{ab}^2\lambda _c)^{1/3}=110140`$ nm . The $`H_{c2}`$ anisotropy was also evaluated from the broadening of the resistive transition in applied magnetic field. Two slightly different theoretical approaches to this problem,  and , based on the anisotropic Ginzburg-Landau and percolation theories for randomly oriented anisotropic superconducting grains, were developed more than a decade ago in early days of the high $`T_c`$ superconductivity. They require somewhat different analysis of the resistivity data but result in similar estimates: $`\gamma _H`$ = 5 - 7 .

The aforementioned determinations of the $`H_{c2}`$ anisotropy utilized traditional techniques, with the exception of those based on CESR and NMR. About a year ago a robust and a simple way to extract this anisotropy from the data on $`T`$ dependence of the magnetization of random powders taken in various fixed fields has been suggested . The method is based on two features in $`(M/T)_H`$: the onset of diamagnetism at $`T_c^{max}`$, which is commonly associated with $`H_{c2}`$ and which is, in fact, $`H_{c2,max}H_{c2,a}`$, and a kink in $`M/T`$ at a lower temperature $`T_c^{min}`$ (see Fig. 1). The origin of these two features can be understood as follows.

Upon cooling a powder sample in a fixed $`H`$, there is a deviation from a roughly $`T`$ independent normal state magnetization (for non-magnetic materials as is the case of MgB2) to an increasingly diamagnetic signal at $`T=T_c^{max}`$. On the $`HT`$ diagram, the point $`(H,T_c^{max})`$ lies at the phase boundary $`H_{c2,max}(T)`$. Repeating this measurement at various fixed field one recovers the full curve $`H_{c2,max}(T)H_{c2,a}(T)`$. This line coincides with the upper critical field determined for polycrystalline samples by, e.g., the standard resistivity measurements. A second sharp feature occurs in $`M/T`$ when the sample temperature passes through $`T_c^{min}`$. This can be readily understood by considering what happens to the sample upon warming. For $`T<T_c^{min}`$ all randomly oriented grains are superconducting, whereas for $`T>T_c^{min}`$ part of them are normal, depending upon their orientation with respect to the applied magnetic field. Therefore, upon warming through $`T_c^{min}`$ we expect some peculiarity, or a kink in $`M/T`$ associated with the onset of normal state properties in an increasing number of appropriately oriented grains. In a similar spirit the field dependent magnetization can be analyzed. This method constitutes a robust procedure independent either of a particular model for describing the anisotropy or of a degree of powder randomness. A kink in $`M/T`$ at $`T=T_c^{min}`$ should be present for any angular distribution of the grains (as long as it is continuous, although not necessarily random).

It should be noted that in superconductors with strong fluctuations (like HTSC) or in the samples with a distribution of $`T_c`$’s or other physical properties due to chemical or/and structural inhomogeneities the region of superconductor-normal transition and the kink in $`M/T`$ may be smeared. In such cases this simple method may not work. Still, given difficulties in growing single crystals of a sufficient size often encountered when a new material is synthesized, the method may prove useful and sometimes the only one available for determination of the $`H_{c2}`$ anisotropy.

The analysis of detailed temperature and field dependent magnetization measurements on MgB2 powders together with the data on polycrystalline $`H_{c2}`$ ($`=H_{c2}^{max}`$) obtained from high field magnetotransport measurements on similar samples yielded a complete anisotropic $`H_{c2}(T)`$ phase diagram , see Fig. 2.

A number of direct measurements of the $`H_{c2}`$ anisotropy on single crystals was reported. Magnetoresistance data on small MgB2 single crystalline platelets  gave $`\gamma _H`$ = 2.6 - 3.0 (in these samples $`T_c`$ = 38.0-38.6 K, $`RRR`$ = 5 - 8). It is argued in Ref.  that electrical transport measurements are not well suited to probe bulk upper critical field of MgB2 due to possible presence of a second phase (related to surface effects) with enhanced superconducting parameters. In this paper, the bulk $`H_{c2}^{ab}(T)`$ and $`H_{c2}^c(T)`$ are evaluated utilizing data on the in-plane thermal conductivity for two orientations of the applied magnetic field. The extrapolated value of $`\gamma _H`$ is $`\gamma _H(T0)`$ = 4.2 .

The torque magnetometry was also employed to study superconducting anisotropy on single crystals . From angular dependent torque measurements taken in different applied fields the anisotropy of the upper critical field was found to be temperature dependent: it decreases from $`\gamma _H`$ 6 at 15 K to $``$ 2.8 at 35 K. Reversible angular dependent torque measurements (with a vortex shaking process ) performed in the temperature range 27 - 36 K and in fields up to 10 kOe were analyzed using a formula derived in  for $`\gamma _H=\gamma _\lambda `$, see Eq. (15) below. As a result, nearly linear field dependence of $`\gamma `$ has been extracted with $`\gamma (H0)`$ 2 and $`\gamma `$(10 kOe) $``$ 3.7 with practically no $`T`$ dependence between 27 K and 36 K. Extensive torque data of Ref.  were also interpreted assuming $`\gamma _H=\gamma _\lambda =\gamma `$ and produced estimates $`\gamma 34`$ which vary somewhat with $`T`$ and $`H`$.

Subsequent analysis  of the anisotropic $`H_{c2}(T)`$ data obtained on random MgB2 powders (Fig. 2) showed $`\gamma _H(T)`$ dependence consistent (Fig. 3) with the one observed in the single crystals . Similar temperature-dependent $`\gamma _H`$ behavior was extracted from magnetotransport, magnetization, ac susceptibility and specific heat in applied magnetic field measurements on single-crystalline platelets  and in thin films .

Summarizing, the reported anisotropy of the upper critical field in magnesium diboride range from $`\gamma _H`$ 1 to $`\gamma _H`$ 13. There are some concerns related to the spread of these results such as purity and homogeneity of materials used, degree of alignment of crystallites/grains, possible sharp angular dependence of $`H_{c2}`$ so that correct evaluation of the anisotropy may suffer from even slight misorientations in measurements on single crystallites, and, last but not least, adequacy of models used to analyze the data on polycrystals and crystals. Even for similar small single crystals different apparently direct measurements give somewhat different estimates of $`\gamma _H`$ (see e.g. ). On the other hand, there is an apparent consistence between recent data on the temperature dependence of $`\gamma _H`$ obtained by at least two groups on single crystals and on high quality powders, similar ”convergence” emerges in data for other physical properties. This indicates that we are coming closer to measuring the intrinsic properties of MgB2. As a side remark it should be mentioned that temperature dependent $`H_{c2}`$ anisotropy is not unique for MgB2, it was observed, for example, in NbSe2 .

For sub-mm size MgB2 single crystal, an anisotropic lower critical field, $`H_{c1}`$, was estimated from the $`M(H)`$ data . The $`H_{c1}^{ab}(T)`$ dependence was reported as near linear with $`H_{c1}^{ab}(0)`$ 384 G, while the $`H_{c1}^c(T)`$ appeared to be non-linear with and $`H_{c1}^c(0)`$ 272 G. Another estimate of the temperature dependent anisotropic $`H_{c1}`$ from magnetization measurements on MgB2 single crystallites  gave close to linear $`H_{c1}(T)`$ for both directions with extrapolated to $`T=`$ 0 values $`H_{c1}^c(0)`$ 2800 G and $`H_{c1}^{ab}(0)`$ 1300 G and near $`T`$ independent ratio $`H_{c1}^c/H_{c1}^{ab}`$ 2.2.

Apparently more work on sample quality, measurement procedures, measurements analysis, and modelling is required to obtain a better physical picture of the superconducting anisotropy of MgB2. Reliable data on the temperature dependence of the anisotropy of the London penetration depth, when available, will serve as an important test of existing models.

## 3 Theory of anisotropies of $`\lambda `$ and $`𝐇_{\mathrm{𝐜𝟐}}`$

To address theoretically the question of anisotropy at arbitrary temperatures one has to turn to a microscopic model. The most complete and sophisticated approach is based on taking into account the actual Fermi surface, the phonon distribution, and the electron-phonon interaction (all of which are anisotropic) in the frame of the Eliashberg theory. This can be and has been done numerically by a few groups, which evaluated the Fermi surfaces  and the gap distribution .

As is shown in Fig. 4, the Fermi surface of MgB2 consists of two nearly cylindrical sheets, ”2D $`\sigma `$-bands”, and two others corresponding to the ”3D $`\pi `$-bands”. The gap on the four Fermi surface sheets of this material has two sharp maxima: $`\mathrm{\Delta }_11.7`$meV at the two $`\pi `$-bands and $`\mathrm{\Delta }_27`$meV at the two $`\sigma `$-bands, see Fig. 5. Within each of these groups, the spread of the gap values is small, and the gaps can be considered as constants, the ratio of which is nearly $`T`$ independent. In this situation, a model with only two Fermi surface sheets and two gaps may prove useful in relating various macroscopic properties of MgB2. Starting with Ref. , the two-band models were studied by many, see, e.g., Ref.  and references therein. We outline below the results of evaluation of the London penetration depth and the upper critical fields within a relatively simple two-gap weak-coupling model which provides consistent physical picture of superconducting anisotropies. We also discuss possible consequences of different anisotropies of $`\lambda `$ and $`H_{c2}`$ upon macroscopic magnetic properties of MgB2.

### 3.1 Two gap model

We begin with the notion that within the weak-coupling BSC theory the macroscopic characteristics of superconductors such as the coherence length and the penetration depth are expressed in terms of quantities averaged over the Fermi surface. This basic consequence of the collective nature of the superconducting condensate is formally expressed by the self-consistency (or gap) equation:

$$\mathrm{\Delta }(𝐫,𝐯,T)=2\pi TN(0)\underset{\omega >0}{\overset{\omega _D}{}}V(𝐯,𝐯^{})f(𝐫,𝐯^{},\omega )_𝐯^{}.$$
(3)
Here, $`N(0)`$ is the total density of states at the Fermi level per one spin; the Fermi velocity $`𝐯`$ relates to a certain position on the Fermi surface, $`\mathrm{\Delta }`$ is the gap function, $`f(𝐫,𝐯,\omega )`$ is Eilenberger Green’s function which describes the condensate, and $`\mathrm{}\omega =\pi T(2n+1)`$ with an integer $`n`$. In the uniform zero-field situation $`f=\mathrm{\Delta }/\sqrt{\mathrm{\Delta }^2+\mathrm{}^2\omega ^2}`$. Substituting this in Eq. (3) one can evaluate how the gap depends on temperature and on position at the Fermi surface provided the pair-pair interaction $`V(𝐯,𝐯^{})`$ is known.

The averages over the Fermi surface weighted with the local density of states $`1/v`$ are defined as

$$X=\frac{d^2𝐤_F}{(2\pi )^3\mathrm{}N(0)v}X(𝐯).$$
(4)
Commonly, the interaction $`V`$ is assumed factorizable , $`V(𝐯,𝐯^{})=V_0\mathrm{\Omega }(𝐯)\mathrm{\Omega }(𝐯^{})`$ (the assumption with far reaching consequences for the two-gap model, see below). Then, one looks for $`\mathrm{\Delta }(𝐫,T;𝐯)=\mathrm{\Psi }(𝐫,T)\mathrm{\Omega }(𝐯)`$. The self-consistency Eq. (3) for the uniform condensate takes the form:

$$1=2\pi TN(0)V_0\underset{\omega >0}{\overset{\omega _D}{}}\frac{\mathrm{\Omega }^2(𝐯)}{\sqrt{\mathrm{\Psi }^2\mathrm{\Omega }^2(𝐯)+\mathrm{}^2\omega ^2}}.$$
(5)
Consider now a model material with the gap anisotropy given by

$$\mathrm{\Omega }(𝐯)=\mathrm{\Omega }_{1,2},𝐯F_{1,2},$$
(6)
where $`F_1,F_2`$ are two sheets of the Fermi surface. Denoting the densities of states on the two parts as $`N_{1,2}`$, and assuming the quantity $`X`$ being constant at each sheet, we obtain for the general averaging:

$$X=\frac{X_1N_1+X_2N_2}{N(0)}=\nu _1X_1+\nu _2X_2,$$
(7)
where we use normalized densities of state $`\nu _{1,2}=N_{1,2}/N(0)`$ for brevity. The Fermi relief of the gap, $`\mathrm{\Omega }`$, can be normalized so that $`\mathrm{\Omega }^2=1`$ which gives for the two-gap model:

$$\mathrm{\Omega }_1^2\nu _1+\mathrm{\Omega }_2^2\nu _2=1,\nu _1+\nu _2=1.$$
(8)
Within this model, the interaction $`V`$ now takes the form of a $`2\times 2`$ symmetric matrix $`V_{ij}/V_0`$ with diagonal components $`\mathrm{\Omega }_1^2`$ and $`\mathrm{\Omega }_2^2`$, and the two equal off-diagonal elements $`\mathrm{\Omega }_1\mathrm{\Omega }_2`$. It is worth noting that if $`V_{12}=0`$, i.e., there is no interband interaction at all, the general gap equation (3) splits in two independent equations for the independent order parameters, $`\mathrm{\Delta }_1`$ on $`F_1`$ and $`\mathrm{\Delta }_2`$ on $`F_2`$. Such a model would imply different, in general, critical temperatures for $`\mathrm{\Delta }`$’s, different phases with many macroscopic consequences . The factorazable interaction $`V_{ij}=V_0\mathrm{\Omega }_i\mathrm{\Omega }_j`$ excludes such possibilities and automatically implies the single $`T_c`$ for both gaps, the feature in fact seen in many tunneling experiments on MgB2 .

### 3.2 $`𝐓`$ dependence of the two gaps

The self consistency equation (5) for the order parameter $`\mathrm{\Psi }(T)`$ can be solved for known $`\mathrm{\Omega }_{1,2}`$ and the densities of states $`\nu _{1,2}`$. Based on the band structure calculations, the relative densities of states $`\nu _1`$ and $`\nu _2`$ of our model are $``$0.56 and 0.44 . The ratio $`\mathrm{\Delta }_2/\mathrm{\Delta }_1=\mathrm{\Omega }_2/\mathrm{\Omega }_14`$, if one takes the averages of $`6.8`$ and $`1.7`$meV for the two groups of distributed gaps as calculated in Ref. . Then, the normalization (8) yields $`\mathrm{\Omega }_1=0.36`$ and $`\mathrm{\Omega }_2=1.45`$. It should be noted that the tunneling data provide somewhat lower ratio $`\mathrm{\Delta }_2/\mathrm{\Delta }_1`$ .

We have now all parameters needed to solve the self-consistency equation (5) for $`\mathrm{\Psi }(T)`$ (the clean case). This is done numerically and the result is shown in Fig. 6 along with two gaps $`\mathrm{\Delta }_{1,2}(T)`$.

The parameter $`\mathrm{\Delta }/T_c`$ commonly used to distinguish strong and weak coupling materials, has a ”strong coupling” value of $`4`$ for $`\mathrm{\Delta }_2/T_c`$, whereas $`\mathrm{\Delta }_1/T_c1`$, which is less than the BCS weak coupling value of 3.5 . Thus, this parameter cannot be used to characterize coupling for materials like MgB2. The $`T`$ dependence of the two gaps similar to shown in Fig. 8 has been reported by a number of groups .

### 3.3 Anisotropy of $`\lambda `$

Having obtained the temperature dependence of both gaps, one can turn to calculation of the weak-field penetration depth $`\lambda (T)`$. To this end, one first solves equations of superconductivity for the zero-field case and then employs the perturbation theory to evaluate the small corrections due to small currents. In this manner one obtains the London relation between the current and the vector potential,

$$\frac{4\pi }{c}j_i=(\lambda ^2)_{ik}^1\left(\frac{2\pi }{\varphi _0}\frac{\theta }{x_k}+A_k\right)$$
(9)
where $`\theta `$ is the phase, $`\varphi _0`$ is the flux quantum, $`𝐀`$ is the vector potential, and

$$(\lambda ^2)_{ik}^1=\frac{16\pi ^2e^2T}{c^2}N(0)\underset{\omega }{}\frac{\mathrm{\Delta }^2v_iv_k}{\beta ^3}$$
(10)
is the tensor of the inverse squared penetration depth. Here, $`i,k=x,y,z`$ and $`\beta ^2=\mathrm{\Delta }^2+\mathrm{}^2\omega ^2`$ (the reader may find details of this calculation in Ref. ). This general expression for the tensor of squared penetration depth is valid for clean superconductors with an arbitrary Fermi surface and gap anisotropies.

The anisotropy parameter now reads:

$$\gamma _\lambda ^2=\frac{\lambda _{cc}^2}{\lambda _{aa}^2}=\frac{v_a^2\mathrm{\Delta }_0^2_\omega \beta ^3}{v_c^2\mathrm{\Delta }_0^2_\omega \beta ^3}.$$
(11)
As $`T0`$, we have $`2\pi T\mathrm{\Delta }_0^2_\omega \beta ^31`$, and

$$\gamma _\lambda ^2(0)=\frac{v_a^2}{v_c^2}.$$
(12)
Note that the gap and its anisotropy do not enter this result. The physical reason for this is in the Galilean invariance of the superfluid flow in the absence of scattering: all charged particles take part in the supercurrent.

Near $`T_c`$, the sum over $`\omega `$ in Eq. (11) is is a $`𝐤_F`$ independent constant because $`\beta \mathrm{}\omega `$, and we obtain the GL result of Ref. :

$$\gamma _\lambda ^2(T_c)=\frac{\lambda _{cc}^2}{\lambda _{aa}^2}=\frac{\mathrm{\Omega }^2v_a^2}{\mathrm{\Omega }^2v_c^2}.$$
(13)
We observe that near $`T_c`$, the gap anisotropy amplifies contribution of the large gap Fermi surface pieces to the anisotropy parameter $`\gamma _\lambda `$. We also see that for isotropic gaps, $`\gamma _\lambda (T_c)=\gamma _\lambda (0)`$ (which is not true for the general case). Thus, the anisotropy parameter depends on $`T`$, the feature absent in superconductors with isotropic gaps.

To apply these results to MgB2, one has to know averages of squared Fermi velocities not only over the whole Fermi surface (as for $`\gamma _\lambda (0)`$ of Eq. (12)), but also the averages over the separate Fermi sheets. These can be taken from the band structure calculations of Ref. :

$`v_a^2_1`$ $`=`$ $`33.2,v_c^2_1=42.2,`$
$`v_a^2_2`$ $`=`$ $`23,v_c^2_2=0.5\times 10^{14}cm^2/s^2.`$ (14)
The numerical results for $`\gamma _\lambda (T)`$ are shown in Fig. 7.

Thus, we expect the London penetration depth of clean MgB2 to be nearly isotropic at low temperatures; $`\gamma _\lambda (T)`$ increases on heating and reaches $`2.6`$ at $`T_c`$. Qualitatively similar behavior of $`\gamma _\lambda (T)`$ is predicted within a more general Eliashberg scheme . As up to date, there is no published data on direct measurements of the $`\lambda `$ anisotropy and its $`T`$ dependence. It is pointed out below that the data on the angular dependence of the torque acting on the anisotropic superconducting crystal in tilted fields cannot provide information about $`\gamma _\lambda (T)`$ without taking into account the $`H_{c2}`$ anisotropy which for MgB2 might be quite different from $`\gamma _\lambda (T)`$. In principle, one can extract $`\gamma _\lambda `$ from the data on the vortex lattice structure in fields tilted relative to the $`c`$ axis of MgB2. These, however, are difficult to acquire and to our knowledge are not yet available.

### 3.4 Anisotropy of $`𝐇_{\mathrm{𝐜𝟐}}`$

Unlike the case of London $`\lambda `$, the evaluation of $`H_{c2}`$ is more complicated a problem because one has to deal with the high-field phase transition from the normal to the superconducting mixed state. Still, major behavior of $`H_{c2}(T)`$ can be obtained within a simplified scheme: Since the actual band structure enters equations for $`H_{c2}`$ via Fermi-surface averages, one can further model the sheets $`F_{1,2}`$ by two spheroids with average Fermi velocities close to the band-structure generated values. As a result one obtains qualitatively (and \- given the spread of existing data - quantitatively) correct behaviors of $`H_{c2}(T)`$ for both principal directions. Details of this evaluation can be found in Ref. . The result for the anisotropy parameter $`\gamma _H(T)=H_{c2,ab}/H_{c2,c}`$ obtained within the scenario of two constant gaps on two Fermi surface sheets, the same used for $`\gamma _\lambda (T)`$ and with the same input parameters, is shown in Fig. 7.

The drop of the $`H_{c2}`$ anisotropy with increasing $`T`$ has been recorded in measurements on single crystals of MgB2 . Nearly the same behavior has been extracted from the data on random polycrystals in the whole temperature range using the $`T`$ dependence of the magnetization of random powders . Recent specific heat measurements  and magnetization data  on single crystals produced similar results. All these data show qualitatively similar behavior to that of the upper curve of Fig. 7.

Physically, the large anisotropy of $`H_{c2}`$ at low temperatures ($`6`$ in our calculation) is related to the large gap value at the Fermi sheet which is nearly two-dimensional. With increasing $`T`$, the thermal mixing with the small-gap states on the three-dimensional Fermi sheet suppresses the anisotropy down to 2.6 at $`T_c`$.

### 3.5  Free energy and torque in materials with different anisotropies of $`\lambda `$ and $`𝐇_{\mathrm{𝐜𝟐}}`$

One of the most sensitive methods used to extract the anisotropy parameter $`\gamma `$ is to measure the torque acting on a superconducting crystal in the mixed state with the applied field tilted relative to the crystal axes. In intermediate field domain, $`H_{c1}HH_{c2}`$, the demagnetization shape effects are weak, and the torque density can be evaluated:

$$\tau =\frac{\varphi _0B(\gamma ^21)\mathrm{sin}2\theta }{64\pi ^2\lambda ^2\gamma ^{1/3}\epsilon (\theta )}\mathrm{ln}\frac{\eta H_{c2,a}}{B\epsilon (\theta )},$$
(15)
where $`\theta `$ is the angle between the induction $`𝐁`$ and the crystal axis $`c`$,

$$\epsilon (\theta )=\sqrt{\mathrm{sin}^2\theta +\gamma ^2\mathrm{cos}^2\theta },$$
(16)
$`\lambda ^3=\lambda _a^2\lambda _c`$, and $`\eta 1`$. This formula can be written as $`\tau =𝐌\times 𝐇`$; since in this field domain the magnetization $`MH`$, one can disregard the distinction between $`𝐁`$ and $`𝐇`$. It has been assumed in the derivation of Eq. (15) that the anisotropies of $`H_{c2}`$ and of the London penetration depth are the same: $`\gamma _H=\gamma _\lambda =\gamma `$. It is worth noting that Eq. (15) describes the system with the stable equilibrium at $`\theta =\pi /2`$, i.e., the uniaxial crystal in the external field positions itself so that the field is parallel to $`ab`$ (if $`\gamma >1`$).

Expression (15) for the torque has been derived within the London approach by employing the cutoff at distances on the order of the coherence length $`\xi `$ where this approach fails. The formula, however, has been confirmed experimentally (for a few high-$`T_c`$ materials close to their $`T_c`$) with a good accuracy as far as the angular dependencies of quantities involved are concerned . Uncertainties of the London approach are incorporated in the parameter $`\eta 1`$; discussion of those can be found, e.g., in Ref. .

We now outline possible effects of different anisotropies of $`\lambda `$ and $`H_{c2}`$ on the torque angular dependence. To this end, we first recall the London expression for the free energy valid for intermediate fields $`H_{c1}HH_{c2}`$ along $`z`$ tilted with respect to the $`c`$ crystal axis over the angle $`\theta `$ toward the crystal direction $`a`$ :

$$F=\frac{B^2}{8\pi }+\frac{B^2m_{zz}}{8\pi \lambda ^2m_a}\underset{𝐆}{}\frac{1}{m_{zz}G_x^2+m_cG_y^2};$$
(17)
here, $`m_{zz}=m_a\mathrm{sin}^2\theta +m_c\mathrm{cos}^2\theta `$, $`m_c/m_a=\gamma _\lambda ^2`$, $`m_a^2m_c=1`$ for uniaxial crystals, and $`𝐆`$ form the reciprocal vortex lattice. The summation is extended over all nonzero $`𝐆`$.

For isotropic case, all reduced masses in Eq.(17) are unity. The sum $`G^2`$ can be approximated by an integral $`(\varphi _0/4\pi ^2B)𝑑G_x𝑑G_y/G^2`$ which diverges logarithmically and should be cut off at a circle in the reciprocal space of a radius $`1/\xi `$ where $`\xi `$ is the radius of the circular vortex core.

In anisotropic material with $`\gamma _\lambda =\gamma _H`$, calculation of this sum is reduced to the isotropic case by rescaling

$$g_x=\sqrt{m_{zz}}G_x,g_y=\sqrt{m_c}G_y.$$
(18)
In this manner one can derive the free energy of a superconductor in an applied field of a certain direction; the torque formula (15) is obtain differentiating the energy with respect to the angle $`\theta `$.

However, for $`\gamma _\lambda \gamma _H`$, the masses in the sum (17) correspond to the anisotropy $`\gamma _\lambda `$ (in other words, the vortex-vortex interaction is determined by the $`\lambda `$ anisotropy) whereas the cut-off reflects the shape of the core (it is determined by the anisotropy of the coherence length $`\xi `$, i.e., by the parameter $`\gamma _H`$). In this situation, the scaling (18) which transforms the summand of Eq. (17) to the isotropic-like form, leaves the shape of the cut-off contour elliptical. Therefore, the energy $`F`$ and its angular dependence acquire extra terms and so does the torque which can be measured.

We refer the reader to the original work  for the details of this calculation and give here only the result. To write explicitly the angular dependence of $`F`$, it is convenient to use the angular functions

$$\mathrm{\Theta }_{\lambda ,H}(\theta )=\epsilon _{\lambda ,H}(\theta )/\gamma _{\lambda ,H}$$
(19)
where $`\epsilon _{\lambda ,H}(\theta )`$ are defined in Eq. (16) with corresponding $`\gamma `$’s. In terms of these functions,

$$F=\frac{B^2}{8\pi }+\frac{\varphi _0B\mathrm{\Theta }_\lambda }{32\pi ^2\lambda _{ab}^2}\mathrm{ln}\frac{2\sqrt{3}\mu _a\varphi _0\mathrm{\Theta }_\lambda }{\xi ^2B(\mathrm{\Theta }_\lambda +\mathrm{\Theta }_H)^2},$$
(20)
where $`\mu _a=\gamma _H^{2/3}`$. The torque density $`\tau =F/\theta `$ follows:

$`\tau `$ $`=`$ $`{\displaystyle \frac{\varphi _0B(\gamma _\lambda ^21)\mathrm{sin}2\theta }{64\pi ^2\lambda ^2\gamma _\lambda ^{4/3}\mathrm{\Theta }_\lambda }}[\mathrm{ln}\left({\displaystyle \frac{\eta H_{c2,c}}{B}}{\displaystyle \frac{4e^2\mathrm{\Theta }_\lambda }{(\mathrm{\Theta }_\lambda +\mathrm{\Theta }_H)^2}}\right)`$ (21)
$``$ $`{\displaystyle \frac{2\mathrm{\Theta }_\lambda }{\mathrm{\Theta }_\lambda +\mathrm{\Theta }_H}}(1+{\displaystyle \frac{\mathrm{\Theta }_H^{}}{\mathrm{\Theta }_\lambda ^{}}})].`$
where $`e=2.718\mathrm{}`$. In the standard case of $`\gamma _H=\gamma _\lambda =\gamma `$, Eq. (21) reduces to the result (15).

Since both $`\mathrm{\Theta }_{\lambda ,H}^{}<0`$, the second contribution to the torque (21) is negative whereas the first one is positive. The positive torque implies that the system energy decreases with increasing $`\theta `$, as in the standard case of $`\gamma _H=\gamma _\lambda `$ for which $`\theta =\pi /2`$ is the stable equilibrium.

On the other hand, it is straightforward to see that for large enough $`\gamma _H`$ and fixed $`\gamma _\lambda `$ (e.g., set $`\gamma _\lambda =1`$ and leave $`\gamma _H>1`$) Eq. (21) yields a negative torque, which means that in the equilibrium the crystal $`c`$ axis is parallel to the applied field. This strange behavior can be tracked down to the angular dependence of the line energy of a single vortex: the energy cost of having an elliptical core for vortices along $`ab`$ for which the current distribution far from the core is near isotropic (as for the example with $`\gamma _\lambda =1`$) is too high.

To illustrate how the angular dependence of the torque varies with anisotropies of $`H_{c2}`$ and $`\lambda `$, we evaluate numerically the torque density (21) for parameters in the range of those for clean MgB2. Fig. 8 shows $`\tau (\theta )`$ for $`\gamma _\lambda =2.2`$ and $`\gamma _H=3`$, the values expected for temperatures somewhat below $`T_c`$. Qualitatively, the dependence is standard; the torque is positive in the whole angular domain, i.e., $`\theta =\pi /2`$ is the stable equilibrium.

With decreasing $`T`$, $`\gamma _H`$ of MgB2 increases whereas $`\gamma _\lambda `$ decreases. In Fig. 9 the torque (21) is plotted for $`\gamma _\lambda =2`$, $`\gamma _H=5`$ (the upper curve) and for $`\gamma _\lambda =1.7`$, $`\gamma _H=5.3`$ (the lower curve). These values correspond roughly to 0.7 and 0.6$`T_c`$ according to Ref. . Clearly, $`\theta =\pi /2`$ as well as $`\theta =0`$ are unstable; the stable equilibrium is shifted to $`0<\theta <\pi /2`$.

Finally, we plot in Fig. 10 the torque density for parameters which correspond to low temperatures, where $`\gamma _\lambda 1.1`$ and $`\gamma _H6`$. The torque is negative for all angles implying the stable equilibrium at $`\theta =0`$.

## 4 Conclusions

Hence, different gaps at different Fermi surface pieces of MgB2 (or generally, anisotropic gaps on anisotropic Fermi surfaces) may lead to profound macroscopic consequences such as those considered above. Moreover, since $`\gamma _H`$ determines the anisotropy of the coherence length and, therefore, of the vortex core, whereas $`\gamma _\lambda `$ describes the ellipticity of the current distribution far from the core, the symmetry of the intervortex interaction should depend on the intervortex spacing, i.e., on the field $`B`$ and its direction. In clean MgB2 at low temperatures in fields along $`ab`$, one expects to have the standard triangular (hexagonal) vortex lattice in low fields ($`\gamma _\lambda 1`$), which should evolve to a distorted triangular (orthorhombic) lattice in increasing fields when the core ellipticity ($`\gamma _H6`$) becomes relevant for the vortex current distribution. The field dependence of the vortex lattice structure should become weaker at elevated temperatures and disappear near $`T_c`$ where $`\gamma _\lambda \gamma _H`$. In intermediate temperature range where the equilibrium vortex lattice orientation shifts from $`\theta =\pi /2`$ to lower angles, one may expect peculiar vortex dynamics (as for fields near parallel to the layers of high-$`T_c`$ materials). This possibility calls for further study.

Thus, the question “what is the anisotropy parameter of MgB2?” does not have a unique answer. To pose the question properly one should specify the quantity of interest. If this is $`H_{c2}`$, the answer is given by the upper curve of Fig. 7 for a clean material; if this is $`\lambda `$, see the dashed line. If this is the magnetization in intermediate fields, $`M(\varphi _0/\lambda ^2)\mathrm{ln}(H_{c2}/H)`$, the main contribution to anisotropy comes from $`\lambda `$; however, the $`H_{c2}`$ anisotropy contributes as well (being smoothed by the logarithm). The last situation should be taken into account while extracting anisotropy from the torque data in tilted fields . It should be noted in closing, that all anisotropies discussed here might be suppressed by impurities.

The authors are grateful to P. Canfield for numerous and helpful discussions. Ames Laboratory is operated for the U. S. Department of Energy by Iowa State University under contract No. W-7405-ENG-82. The work was supported by the Director of Energy Research, Office of Basic Energy Sciences, U. S. Department of Energy.


text:
# Instabilities of off-centered vortices in a Bose-Einstein condensate

## I Introduction

Quantized vorticity which accompanies rotational motion is a characteristic feature of superfluids. Vorticity involves singular lines, vortex lines which exist inside any closed path enclosing finite circulation. The change in the phase of the order-parameter field around the vortex line is quantized in integer units of $`2\pi `$.

In a Bose-Einstein condensate (BEC) of an atomic gas pethick , the first vortex was created by the JILA group matthews . They used a two-component condensate of 87Rb atoms and imprinted a phase winding of $`2\pi `$ onto one component. The position of the vortex core was clearly seen and also the precessional motion of the core was observed matthews ; anderson . In a one-component condensate, the ENS group ens1st  created a vortex state using a rotating trap which consists of an optical spoon and a magnetic trap. Condensates containing up to 4 vortices were observed. Their method can be understood in analogy with fluid motion in a rotating vessel. The idea of a rotating trap has been employed by many groups thereafter.

In addition to these singly quantized vortices, multiply quantized vortices with phase windings $`4\pi `$ and $`8\pi `$ can also be formed topological , using a topological method introduced by Nakahara et al. nakahara . However, the experimentally recorded image of the cloud of atoms does not show any indications for the splitting of multiply quantized vortices. This surprising apparent stability contradicts with what has been widely expected pethick ; pu . Therefore, it is urgent to understand the stability of these off-centered vortices and the multiply quantized vortices.

The axisymmetric and rotationally symmetric (singly quantized) vortex states have been studied and a core-localized excitation with negative excitation energy is found dodd ; tomoyaaxisym1 ; tomoyaaxisym2 ; sami ; sami:tpopov ; svidzinsky . Concerning this excitation, a relation to the instability of a vortex state tomoyaaxisym2 ; svidzinsky  (through the transfer of the population from the condensate mode to the core-localized excitation) and also to the precession motion linn ; feder  of the vortex core have been pointed out. Nevertheless, the axisymmetric and rotationally symmetric systems only mimic the experimental situations because (a) the axisymmetric vortex state can become a vortexfree state only through the off-centered vortex states. (b) a system with a vortex under precessional motion is no longer axisymmetric.

The understanding of the excitation level with negative energy involves a puzzling feature comparison  because the value is connected to both the direction and frequency of the precession motion and the instability of the vortex state. These two cannot be divided as long as we study the axisymmetric state. Therefore, investigations of the off-centered vortex state are necessary in order to understand both the instability and the precessional motion. The off-centered vortex state has been studied through the Gross-Pitaevskii equation butts  and the Thomas-Fermi (TF) approximation svidzinsky ; guilleumas . But no solutions of the Bogoliubov equations pethick  have thus far been presented for a non-axisymmetric vortex state.

The motivation for an analysis of non-axisymmetric doubly quantized vortex may be explained in a similar way: The instability of the doubly quantized vortex state has been pointed out for the rotationally symmetric vortex state pu . It relates the existence of a mode with a complex excitation energy and its exponential increase of population. But the axisymmetric state only lasts for an infinitesimal time if the vortex is unstable. It is not known whether the system keeps having the feature (complex eigenvalue) on which the instability of doubly quantized vortex state depends.

This paper presents a two-dimensional analysis of the excitation spectra supported by the condensate with an off-centered vortex, a doubly quantized vortex, and pairs of split vortices. We discuss the above questions of precession, the instability of singly quantized vortex, and the doubly quantized vortex within a unified numerical framework.

## II Single off-centered vortex

We consider a two-dimensional $`(x,y)`$ system in the rotating frame whose angular velocity is $`\omega `$. The axis of rotation is perpendicular to the $`(x,y)`$ plane. Ideally, this represents not only the static condensate in the rotating vessel but also the precessing vortex. The Hamiltonian is

$`\widehat{H}(𝝎)`$ $`=`$ $`\widehat{H}𝝎\widehat{L}`$ (1)
$`\widehat{H}`$ $`=`$ $`{\displaystyle \widehat{\mathrm{\Psi }}^{}(C^2+V)\widehat{\mathrm{\Psi }}}+{\displaystyle \frac{1}{2}}g\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}^{}\widehat{\mathrm{\Psi }}\widehat{\mathrm{\Psi }}d𝐫`$ (2)
$`\widehat{L}`$ $`=`$ $`{\displaystyle \widehat{\mathrm{\Psi }}^{}(𝐫\times 𝐩)\widehat{\mathrm{\Psi }}𝑑𝐫}`$ (3)
where $`C=\mathrm{}^2/(2m)`$, $`g=4\pi \mathrm{}^2a/m`$, and $`V`$ is the confining potential with $`a`$ denoting the s-wave scattering length and $`m`$ the mass. We use $`𝝎=(0,0,\omega )`$ and $`𝐫=(x,y,0)`$. The time-dependent Gross-Pitaevskii (GP) equation may be written as

$`\{C^2+V(x,y)+g|\varphi (x,y,t)|^2`$
$`𝝎𝐫\times 𝐩\}\varphi (x,y,t)`$ $`=`$ $`i\mathrm{}{\displaystyle \frac{\varphi (x,y,t)}{t}},`$ (4)
where $`\varphi (𝐫)`$ is the condensate wavefunction. The time-independent form is

$`\{C^2+V(x,y)\mu +g|\varphi (x,y)|^2`$
$`𝝎𝐫\times 𝐩\}\varphi (x,y)`$ $`=`$ $`0,`$ (5)
where $`\mu `$ is the chemical potential. We employ an axisymmetric harmonic trapping potential $`V(x,y)=m\omega _{\mathrm{tr}}^2(x^2+y^2)`$ where $`\omega _{\mathrm{tr}}=2\pi \times 200`$. The linear number density of particles $`N=|\varphi (𝐫)|^2𝑑𝐫`$ is fixed to $`1\times 10^{10}m^1`$ in this paper. The mass $`m=38.17\times 10^{27}\text{ kg}`$ and the scattering length $`a=2.75\text{ nm}`$ for Na atoms are employed.

The angular momentum of the condensate per particle

$$L=\frac{\varphi ^{}(𝐫\times 𝐩)\varphi 𝑑𝐫}{\mathrm{}|\varphi (𝐫)|^2𝑑𝐫}$$
(6)
is extensively utilized in the following discussion and plotted in several figures. A condensate with a centered vortex has $`L=1`$. Systems with an off-centered vortex have $`L`$ between 0 and 1. We also use the TF radius $`R_{\mathrm{TF}}=6.79\mu \mathrm{m}`$ as the unit of length.

### II.1 Static Analysis

The time-dependent Bogoliubov equations sami:tpopov

$`i\mathrm{}{\displaystyle \frac{u(x,y,t)}{t}}`$ (7)
$`=`$ $`(C^2+V\mu +2g|\varphi |^2𝝎𝐫\times 𝐩)u(x,y,t)`$
$`g\varphi ^2v(x,y,t),`$
$`i\mathrm{}{\displaystyle \frac{v(x,y,t)}{t}}`$ (8)
$`=`$ $`\left(C^2+V\mu +2g|\varphi |^2+𝝎𝐫\times 𝐩\right)v(x,y,t)`$
$`+g\varphi ^2u(x,y,t),`$
yield the excitation spectra of the condensate which follows the time-dependent GP equation Eq. (4). They reduce to the Bogoliubov equations pethick

$`\left(C^2+V\mu +2g|\varphi |^2𝝎𝐫\times 𝐩\right)u`$
$`g\varphi ^2v`$ $`=`$ $`\epsilon u,`$ (9)
$`\left(C^2+V\mu +2g|\varphi |^2+𝝎𝐫\times 𝐩\right)v`$
$`+g\varphi ^2u`$ $`=`$ $`\epsilon v.`$ (10)
when the system is static.

The criteria for the applicability of the Bogoliubov Eqs. (9) and (10), instead of the time-dependent Bogoliubov Eqs. (7) and (8), depends on how fast the wavefunction of the condensate changes. We measure the amplitude of the transformation using

$$V^{}\frac{\mathrm{max}\left(\left|(C^2+V+g|\varphi |^2𝝎𝐫\times 𝐩\mu )\varphi \right|\right)}{\mathrm{}\mathrm{max}(|\varphi |)}$$
(11)
where the chemical potential $`\mu `$ is fixed here chemical1 ; chemical2 . The numerator is a maximum of the left-hand side of Eqs. (4) and (5). Therefore, in the framework of the time-dependent GP equation

$$V^{}(t)=\mathrm{max}\left(\left|\frac{\varphi (x,y,t)}{t}\right|\right)/\mathrm{max}(|\varphi (x,y,t)|).$$
(12)
For variations of the displacement of the vortex core $`r`$, the $`V^{}`$ has local minima for an axisymmetric vortex state, a vortexfree state, and an off-centered vortex state with certain displacement $`r`$ for a choice of $`\omega `$.

We use off-centered vortex states which satisfy

$$V^{}<0.0005$$
(13)
in the following analysis of the Bogoliubov equations. Equations (12) and (13) imply that the wavefunction $`\varphi `$ varies slowly enough in the time scale of the harmonic trap \[$`V^{}<0.0005\omega _{\mathrm{tr}}/(2\pi )`$\]. Instead of the right-hand side of Eq. (13), for low values of the angular momentum ($`L<0.2`$) the right-hand side is rather relaxed to $`0.03`$.

Each of the relaxation processes, which is numerically equivalent to following the time-dependent GP equation with imaginary time, begins with the TF density profile as the initial configuration on which the phase of an off-centered vortex with various displacements $`r`$ is “imprinted” through multiplication. In the beginning of a numerical simulation, the angular velocity $`\omega `$ is chosen. Then the development with fixed $`\omega `$ is performed without restricting the displacement. The vortex moves freely during the development and ends as an axisymmetric vortex, a vortexfree state, or a vortex state which satisfies Eq. (13). The latter form a set of $`\varphi `$ and $`\omega `$ values which fulfill Eq. (13). Figure 1(a) shows some of the density profiles $`|\varphi |^2`$ of the condensate. The displacement of the vortex core $`r`$ and the angular momentum $`L`$ are calculated through Eq. (6) from each of the resulting wavefunctions $`\varphi `$. Therefore, $`L`$ and $`r`$ are functions of $`\omega `$ through $`\varphi `$. The relation between the angular momentum and the displacement is plotted in Fig. 1(b).

Figure 1(c) displays the normalized additional energy due to the existence of an off-centered vortex, defined through

$$\mathrm{\Delta }E^{}=(EE_0)/(E_1E_0)$$
(14)
where $`E_0`$, $`E_1`$, and $`E`$ are the energies of a vortexfree condensate, a condensate in the presence of a centered vortex, and with an off-centered vortex, respectively. If the condensate is vortexfree, $`\mathrm{\Delta }E^{}=0`$. Off-centered vortex states have $`0<\mathrm{\Delta }E^{}<1`$ and the centered vortex state has $`\mathrm{\Delta }E^{}=1`$.

### II.2 Excitation Spectra

Excitations from a BEC have been interpreted in the framework of the Bogoliubov equations in Ref. vogels . They play a significant role in determining vortex stability. Nevertheless, these modes in off-centered vortices have not been analyzed so far. An excitation whose energy is $`\epsilon `$ is expressed with two components, $`u`$ and $`v`$, of wavefunction each with different momenta. They are determined as the solutions to the eigensystem of Bogoliubov Eqs. (9) and (10). The modes with negative $`(|u|^2|v|^2)𝑑𝐫`$ are ignored. Two of the lowest excitations are the condensate mode (with $`\epsilon =0,u=\varphi ,v=\varphi ^{}`$) and the core-localized mode. We ignore the mode with the smallest $`|\epsilon |`$ as the condensate mode.

The wavefunction expressed with $`u`$ and $`v`$ are obtained through the Bogoliubov equations (9) and (10). We denote their angular momenta using

$`𝒜(u)`$ $``$ $`{\displaystyle \frac{u^{}(𝐫\times 𝐩)u𝑑𝐫}{\mathrm{}|u|^2𝑑𝐫}},`$ (15)
$`𝒜(v)`$ $``$ $`{\displaystyle \frac{v^{}(𝐫\times 𝐩)v𝑑𝐫}{\mathrm{}|v|^2𝑑𝐫}}.`$ (16)
These quantities correspond to the angular momentum of the condensate $`L`$ defined in Eq. (6). Some algebra with the help of Eqs. (5), (9), and (10) shows that the excitation energies $`\epsilon `$ depend linearly on $`\omega `$. We presume that the wavefunctions $`\varphi `$, $`u`$, and $`v`$ do not depend on $`\omega `$ and energies $`\mu `$ and $`\epsilon `$ do not have imaginary part. The coefficient of proportionality is

$$q_\theta \frac{\mathrm{Re}\left[\{𝒜(u)L\}|u|^2𝑑𝐫+\{𝒜(v)+L\}|v|^2𝑑𝐫\right]}{(|u|^2+|v|^2)𝑑𝐫}.$$
(17)
We introduce an excitation energy in the laboratory frame

$$\epsilon _{\mathrm{lab}}\epsilon +\mathrm{}\omega q_\theta .$$
(18)
Equation (17) includes $`𝒜(u)`$, $`𝒜(v)`$, and $`L`$ because an excitation level has two components $`u`$ and $`v`$ of wavefunction and their phases are always affected by the phase of the condensate. For axisymmetric systems, the $`q_\theta `$ definition reduces to an integer quantum number tomoyaaxisym1 ; tomoyaaxisym2 . The excitation energies $`\epsilon `$ and $`\epsilon _{\mathrm{lab}}`$ are normalized by the trap unit $`\mathrm{}\omega _{\mathrm{tr}}`$ from here on. Figure 2(a) displays the computed excitation energies, $`\epsilon `$, and the corresponding angular momenta, $`q_\theta `$. Some of the modes, for example those with $`(\epsilon _{\mathrm{lab}},q_\theta )=(2,0),(1,\pm 1)`$ maintain similar values of $`q_\theta `$ and $`\epsilon _{\mathrm{lab}}`$ upon variations of $`L`$. We call the modes at $`(2,0)`$ and $`(1,\pm 1)`$ the breathing pethick  and the dipole modes, respectively. Some of the core-localized modes have $`(\epsilon _{\mathrm{lab}},q_\theta )=(0,0)`$. These are discussed in Sec. II.3. The modes with $`(\epsilon _{\mathrm{lab}},q_\theta )(1,\pm 1),(1.5,\pm 2)`$ are discussed in Sec. II.4.

### II.3 Core-Localized Mode

For any value of $`L`$, the system has one core-localized mode. The modes always displays a sharp peak in the amplitude $`|u|^2+|v|^2`$ of the wavefunction at the vortex core. The core-localized modes still feature a major difference between the system with $`L=1`$ and those with $`0<L<1`$.

When $`L=1`$, the vortex is in the center of the system. The angular momentum $`q_\theta =1`$ and the excitation energy $`\epsilon _{\mathrm{lab}}`$ is negative \[the lowest bullet at $`q_\theta =1`$ in Fig. 2(a)\]. It has a sharp peak in the vortex core. The sharp localization, the angular momentum, and the negative energy are equivalent with those for axisymmetric calculations tomoyaaxisym1 ; tomoyaaxisym2 .

Once the vortex becomes off-centered ($`0<L<1`$), both of the components $`u`$ and $`v`$ of the wavefunction localize in the vortex core \[Fig. 3(a)\]. Each of them features a sharp peak in the vortex core and its skirt extends towards the surface of the condensate. The angular momenta $`𝒜(u)`$ of the core-localized modes vary between 1 and 4. The corresponding wavefunctions $`v`$ have opposite angular momenta. Their average $`q_\theta `$ almost vanishes ($`|q_\theta |1`$); the points at (0,0) in Fig. 2(a) mean this. The excitation energy plotted in Fig. 2(d) remains below zero ($`0.01<\epsilon <0`$) even in the rotating frame.

The core-localized mode behaves differently between the system with $`L=1`$ and that with $`0<L<1`$. One possible cause is the static approximation of the time-dependent Bogoliubov equations introduced in Sec. II.1. While the system with $`L=1`$ is treated with the Bogoliubov equations, systems with $`0<L<1`$ are described within the time-dependent Bogoliubov equations.

### II.4 Dipole and Quadrupole Modes

Excitations having the lowest positive energy at $`q_\theta =\pm 1`$ and $`\pm 2`$ are classified as the dipole and quadrupole modes. Figure 2(b) shows the energy for these modes. The energies of the dipole modes almost equal one trap unit $`(\mathrm{}\omega _{\mathrm{tr}})`$. The dipole modes are known to be responsible for the center-of-mass motion zambelli . Therefore, they are not affected by various profiles inside the condensate. The angular momenta $`𝒜(u)`$ and $`𝒜(v)`$ are plotted in Fig. 2(c). These excitations show continual transformation as a function of variation of the vortex displacement.

The condensate in a two-dimensional harmonic potential has two quadrupole modes. The splitting between the quadrupole frequencies $`\epsilon _{\mathrm{lab}}(+2)\epsilon _{\mathrm{lab}}(2)`$ is inverted at $`L0.4`$, while it is proportional to $`L`$ in systems with a centered vortex zambelli . The corresponding eigenoscillation of the condensate, called the scissors mode, has been observed experimentally chevy  to measure an angular momentum of the condensate. The inversion may affect the observations of the angular momentum.

Guilleumas et al. guilleumas  estimate the splitting using several methods within the TF limit. Their Eqs. (33) and (42) and our numerical results in Fig. 2(b) agree well when the vortex is close to the axis ($`r/R_{\mathrm{TF}}<0.2`$). Their estimates stay positive and approach zero as the displacement of the vortex approaches $`R_{\mathrm{TF}}`$. But the sign is inverted in our numerical results. This is caused by an increase of the energy $`\epsilon _{\mathrm{lab}}(2)`$ at lower $`L`$, depicted in Fig. 2(b). The wavefunction $`u`$ of the $`2`$ mode displays a sharp peak at the vortex center; for example, in Fig. 3(d). This peak is significant for lower $`L`$ and not seen in the axisymmetric (non)vortex states. We presume that the existence of this peak is pushing up the energy $`\epsilon _{\mathrm{lab}}(2)`$.

For $`\omega 0.35\omega _{\mathrm{tr}}`$ (corresponding to $`L<0.2`$), the excitation energies of the modes with $`q_\theta >2`$ approach each other, and their wavefunctions become numerically indistinguishable from one another. The angular momenta of the excitations in this region ($`L<0.2`$) are omitted in Figs. 2 for the sake of clarity.

## III Two vortices

A vortex with the winding number two has been formed in a condensate of Na atoms topological . Multiply quantized vortices are energetically unfavorable against the array of the singly quantized vortices (e.g. Sec. 9.2.1 in Ref. pethick ), while the experimental data does not indicate any signs of the splitting. In this section, we investigate the excitation spectra of systems with two vortices to inspect their stability in the splitting process.

Each of the GP calculations starts from an initial wavefunction with two slightly displaced vortices under a fixed angular velocity $`\omega `$, which acts as an external adjustable parameter. The initial wavefunction is formed from the TF density profile in the absence of a vortex, but multiplied with the phase factor $`e^{i\theta }`$ around each of the initial vortex positions. The resulting angular momenta are controlled through this $`\omega `$. The GP equation has a static solution for about $`60<\omega /(2\pi )<110`$. The corresponding range of angular momenta is from 1.2 to 1.8. Many vortices enter into the condensate for higher $`\omega /(2\pi )>110`$ and the vortices disappear for $`\omega /(2\pi )<60`$. Figure 4 shows density profiles of the condensate for various angular momenta.

Using the condensate wavefunctions, the excitations from it are calculated with the Bogoliubov equations, Eqs. (9) and (10). We found that one or two conjugate pairs of modes with complex eigenvalues always exists in the calculated range of $`L`$.

Figure 5(d) shows the energy levels for two excitations with the lowest real part in the rotating frame. For $`1.2<L<1.57`$, there occur two conjugate pairs of complex eigenvalues, both of them have a large imaginary part. For larger $`L`$, only one of them has a large imaginary part. Figure 6 displays the wavefunctions of these lowest modes. One of them has an imaginary eigenvalue, while the other one has a real eigenvalue. Both of them display a peak in the vortex core.

A doubly quantized vortex in a rotationally symmetric system has two negative excitations, those for $`q_\theta =1`$ and $`2`$. One of them with $`q_\theta =2`$ sometimes has a complex energy, depending on the particle density in the system, see Ref. pu . This feature is equivalent to those of our results depicted in Fig. 5(d) for higher $`L`$. Therefore, the splitted vortices keep on having unstable nature which exist before the splitting pu .

Figure 5(a) shows that the modes with positive energy stay with similar values of $`\epsilon _{\mathrm{lab}}`$ and $`q_\theta `$ as functions of the separation of the vortices. The lowest positive-energy modes, with $`q_\theta \pm 1`$ and $`q_\theta \pm 2`$, are classified as the dipole and quadrupole modes. Their computed angular momenta and energies in the laboratory frame are presented in Fig. 5(b) and 5(c). The modes near $`(q_\theta ,\epsilon _{\mathrm{lab}})=(0,2)`$ are breathing modes. The consistent behavior of the dipole, quadrupole, and breathing modes throughout the results in systems with an off-centered vortex \[Figs. 2(a)-2(c)\], two vortices \[Figs. 5(a)-5(c)\], and an axisymmetric vortex ($`L=1,2`$) proves the validity of our numerical procedure.

## IV Precession Frequencies

The above results are calculated in the rotating frame. It is not self-evident whether they may be interpreted as results for precessing vortices in the laboratory frame. To estimate the validity of the condensate wavefunction, we carry out time-dependent Gross-Pitaevskii calculations in the laboratory frame. The above results of the static GP equation are allowed to evolve in time according to Eq. (4).

The condensate shows precessional motion of the vortices for both the single-vortex and the two-vortex cases. Figure 7 shows the precessional angular velocity, $`\omega _{\mathrm{pr}}`$. The rotation frequency $`\omega `$ that we obtained in the static calculations in Secs. II and III are also plotted. The two frequencies agree well. It means that our results in the rotating frame can also be identified as those in the laboratory frame. The excitation modes, including the core mode, precess following the vortex core there.

## V Discussion

The excitation spectra of Bose Einstein condensates with one and two off-centered vortices has been computed. The excitation with negative energy and localization at the vortex core always exists not only in the axisymmetric state but also for the off-centered vortex configurations. We depicted the shape of wavefunctions and excitation energy of the core-localized, dipole and quadrupole modes in the off-centered vortex states with both one and two vortices. The invariance of dipole and breathing modes in the nonaxisymmetric cases shows that usage of of the Bogoliubov equation instead of the time-dependent Bogoliubov equation is valid.

The rotation frequencies agree well with the precession frequencies of vortices in the laboratory frame. It enables us to consider the two problems, precession of the vortices and the vortex instability concerning the negative excitation, independently. The precessing vortices pursue to have a core-localized excitation with negative excitation energy in the framework of Bogoliubov equations. Therefore, the condensate with one precessing vortex maintains the instability as long as the vortex exists. In order to pose an answer to the problem comparison  between the precessing mode and the instability, this work is to be extended to treat finite-temperature systems tomoyaaxisym1 ; sami .

Concerning the splitting of the doubly quantized vortex, the existence of two core-localized excitations is confirmed throughout the calculated range of $`L`$. The vortex pair always has an instability pu ; skryabin  because systems with a vortex pair always support core excitations with complex excitation energies, c.f. Fig. 5(d).

Throughout the analyses, the system displays two quadrupole excitations. The excitation energies of the quadrupole modes Figs. 2(b) and 5(b) are useful because these excitation is used to measure the angular momenta of the condensate experimentally.

## Acknowledgements

Authors thank CSC Scientific Computing Ltd (Finland) for computer resources. We are grateful to M. Möttönen, M. Nakahara, C. J. Pethick, T. P. Simula, and S. M. M. Virtanen for stimulating discussions. One of the authors (TI) is supported by the bilateral exchange programme between the Academy of Finland and the Japan Society for the Promotion of Science.