text:
# Molecular Dynamics Study of the Nematic-Isotropic Interface

## I Introduction

Surfaces and interfaces in liquid crystals have been the subject of much interest both from a fundamental point of view and because of their practical importance in the context of liquid crystal devices .

Liquid crystals are formed by anisotropic particles. Depending on temperature and density, they exhibit various liquid phases. Here, we shall be concerned with the isotropic phase (I), which is an ordinary fully symmetric fluid phase, and the nematic phase (N), where the fluid retains the translational symmetry in all directions (no positional order), but exhibits long range orientational order. According to a very general symmetry argument, the transition between the isotropic and the nematic phase is first order. Therefore there exist regions in phase space where nematic and isotropic fluids coexist and are separated by interfaces.

The nematic state is characterized by a so-called nematic director, which specifies the preferred direction of alignment of the particles. Surfaces and interfaces break the isotropy of space and usually favour a certain director orientation. This effect, called surface anchoring, is commonly characterized by the tilt angle $`\theta `$ between the preferred orientation and the surface (interface) normal, and the anchoring strength or anchoring energy. In particular, the cases $`\theta =0^{}`$ (perpendicular alignment) and $`\theta =90^{}`$ (parallel alignment) are often referred to as homeotropic anchoring and planar anchoring, respectively.

Theoretically, particular attention has been given to NI-interfaces in systems of hard particles. Studies of NI-interfaces in hard rod systems based on Onsager’s approach  have indicated that the surface free energy has a minimum when the rods lie parallel to the interface, suggesting that the anchoring in this system is planar. Similar results have been obtained in systems of hard ellipsoids. According to the Onsager theory, the interfaces are flat, i.e., fluctuations of the interface position (capillary waves) are not present.

In general, however, one would expect such fluctuations in fluid-fluid interfaces. The usual simple argument reads as follows: On very large length scales, the free energy of a system containing an interface should be governed by the interfacial tension $`\gamma `$. Let us neglect bubbles and overhangs and parametrize the local position of a fluctuating planar interface by a single-valued function $`h(x,y)`$. Fluctuations enlarge the interfacial area and thus cost the free energy

$$_{\mathrm{CW}}\{h\}=\frac{\gamma }{2}𝑑x𝑑y\left[\left(\frac{h}{x}\right)^2+\left(\frac{h}{y}\right)^2\right],$$
(1)
(assuming small distortions with $`|h/x|`$, $`|h/y|1`$). Since this free energy functional is quadratic, the partition function can be evaluated exactly. $`\{h\}`$ is easily diagonalized by means of a two dimensional Fourier transform. One finds that the squared amplitude of fluctuations with wavevector $`𝐪`$ is on average given by

$$|h(𝐪)|^2=\frac{k_\mathrm{B}T}{\gamma q^2},$$
(2)
i.e., diverges in the long-wavelength limit $`𝐪0`$. The argument thus not only predicts on fairly general grounds the existence of capillary waves at all temperatures, but also that they are actually quite large.

On the other hand, it has been argued that the Onsager approach should be exact in the limit of infinitely elongated particles. The Onsager theory assumes that the structure of the fluid is entirely determined by the second virial coefficient, i.e., by the statistical properties of clusters of two particles. In systems of very long, very thin particles, the probability that one particle interacts with more than one other particle approaches zero, and the Onsager assumption becomes exact. One concludes that capillary waves should vanish in this limit, in contradiction to the argument presented just above.

How can this be? Two explanations are conceivable. First, capillary waves vanish if the interfacial tension diverges. However, this is unlikely here, since the coexistence density of the nematic and the isotropic phase moves towards zero as the particle elongation is increased. Second, capillary waves might be suppressed by long-range interactions, which are disregarded in the functional (1). Indeed, elastic long-range interactions are present in the nematic phase. We recall that the interface orients the director of the nematic liquid in a certain direction. Long-range interface fluctuations therefore lead to long-range elastic distortions of the director field, which are penalized.

In order to shed more light on these issues, we have carried out extensive Molecular Dynamics simulations of an N-I interface in a system of over 100000 particles of elongation 15, and analyzed the capillary wave spectrum in detail.

To our knowledge, there exist only few numerical studies of N-I interfaces . In one of them , a system of molecules with length-to-width ratio $`\kappa =3`$ interacting via a Gay-Berne potential was studied in a simulation box with a temperature inhomogeneity, which served to make and maintain the NI interface. The molecules in the nematic phase were found to align parallel to the interface (planar anchoring). However, due to the use of the temperature inhomogeneity, the system is not at thermodynamic equilibrium. Capillary waves of the interface position fluctuations are obviously suppressed. In the other studies, repulsive ellipsoids of revolution of axis ratio $`\kappa =15`$ were used in a system confined between two hard walls. The interactions of the walls with molecules were systematically varied to study the interplay between surface and interface anchoring. Planar alignment at the N-I interface was again observed. The system sizes considered were too small to allow for an analysis of capillary waves.

The present work builds on Reference . We study NI-interfaces in a very large system of repulsive ellipsoids by means of extensive computer simulations, focusing in particular on the capillary wave fluctuations of the interfaces. The paper is organized as follows: In the next section, we describe the model and some simulation details . Section III explains our way to determine the local position of the interface. The results are described and summarized in section IV.

## II Model and Simulation Details

Our system consists of idealized ellipsoidal particles with elongation $`\kappa =\sigma _l/\sigma _s=15`$ where $`\sigma _l`$ and $`\sigma _s`$ are the length and width of the particles, respectively. The large value of the elongation $`\kappa =15`$ ensures that the order of the transition is fairly strong, which makes it easier to generate and maintain the interface. The interaction between two ellipsoids $`i`$ and $`j`$ is given by

$$V(𝐮_i,𝐮_j,𝐫_{ij})=\{\begin{array}{cc}4ϵ(X^{12}X^6+\frac{1}{4})\hfill & \text{ }X^6>\frac{1}{2}\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array}$$
(3)
with

$$X=\frac{\sigma _s}{r\sigma (𝐮_i,𝐮_j,\widehat{𝐫}_{ij})+\sigma _s}$$
(4)
Here, $`𝐮_i`$ and $`𝐮_j`$ are the orientations of ellipsoids $`i`$ and $`j`$, and $`𝐫_{ij}`$ is the center-center vector between the ellipsoids, of magnitude $`r`$ and direction $`\widehat{𝐫}_{ij}`$. The distance function $`\sigma `$  approximates the contact distance between two ellipsoids and is given by

$`\sigma (𝐮_i,𝐮_j,\widehat{𝐫})`$ $`=`$ $`\sigma _s\{1{\displaystyle \frac{\chi }{2}}[{\displaystyle \frac{(𝐮_i\widehat{𝐫}+𝐮_j\widehat{𝐫})^2}{1+\chi 𝐮_i𝐮_j}}`$ (6)
$`+{\displaystyle \frac{(𝐮_i\widehat{𝐫}𝐮_j\widehat{𝐫})^2}{1\chi 𝐮_i𝐮_j}}]\}^{1/2},`$
where the parameter $`\chi `$ is related to the elongation $`\kappa \sigma _l/\sigma _s`$ through

$$\chi =\frac{\kappa ^21}{\kappa ^2+1}.$$
(7)
Since the attractive tail of the potential eqn. (3) has been cut off, its range is much shorter than that of the Gay-Berne potential. This enables us to simulate systems with a large number of particles, which is essential to study the capillary wave effects of interest to us. For convenience, $`\sigma _s=1`$ defines a unit of length, $`ϵ=1`$ a unit of energy, and we take particle mass $`m=1`$; the particle moment of inertia is set at $`I=50m\sigma _s^2`$.

As a preliminary run, we performed a molecular dynamics (MD) simulation of a system with 7200 molecules at constant volume in a box geometry $`(L_x:L_y:L_z)=(1:1:8)`$ with periodic boundary conditions in all directions. The density was chosen in the coexistence region, $`\rho =0.017/\sigma _s^3`$, and the temperature $`k_\mathrm{B}T/ϵ=1`$. A nematic slab bounded by two interfaces parallel to the $`xy`$ plane assembled as a result. More than $`5.5\times 10^6`$ MD steps were required to equilibrate this system, where each step covers 0.002 time units. The last configuration was then reproduced $`4`$ times in the $`x`$ and $`y`$ direction, which generated a system with $`N=115200`$ molecules and box size $`L_x=L_y=150.1958\sigma _sL`$ and $`L_z=300.3916\sigma _s`$. This was used as the starting point for a MD simulation of $`4.2\times 10^6`$ MD steps.

## III Block Analysis and The Dividing Surface

In order to study the interfacial position fluctuations as a function of the system size, one can either perform simulations of various system sizes or simulate a single (large) system and analyze subsystems of it. In this study, we took the latter approach. We split our system of size $`L\times L\times L_z`$ into columns of block size $`B\times B`$

and height $`L_z`$. The columns are further divided into $`N_{\mathrm{bins}}=200`$ bins in the $`z`$ direction. The local position of the interface $`z_{\mathrm{int}}(x,y)`$ in each column is estimated as follows.

First, we compute the local ordering tensor $`𝐒`$ in each bin (of size $`B\times B\times (L_z/N_{\mathrm{bins}})`$), defined as

$$S_{\alpha \beta }=\frac{1}{n}\underset{i=1}{\overset{n}{}}\frac{1}{2}(3u_{i\alpha }u_{i\beta }\delta _{\alpha \beta }),$$
(8)
where $`u_{i\alpha }`$ is the $`\alpha `$ component of the unit vector which points along the axis of the $`i`$-th molecule ($`\alpha =x,y,z`$), and $`n`$ is the number of molecules in the bin. The nematic order parameter $`S(z)`$ in the bin centered at $`z`$ is given by the largest eigenvalue of the ordering tensor. The further procedure is illustrated in Fig. 1. Once the order parameter profile $`S(z)`$ is obtained (solid line), we compute at least two averaged profiles over $`N_{\mathrm{ave}}`$ bins, where $`N_{\mathrm{ave}}`$ is chosen such that an averaged profile still reflects the existence of the interface, but short range fluctuations are averaged out (dotted and dashed lines in Fig. 1). Finally, the position of the interface is estimated as the intersection of two averaged profiles.

The method is motivated by the following consideration: If the order parameter $`S(z)`$ did not fluctuate at all in the bulk, the intersection would locate a “dividing surface”, where the negative order parameter excess on the nematic side just balances the positive order parameter excess on the isotropic side. It should then be independent of $`N_{\mathrm{ave}}`$, as long as $`LN_{\mathrm{ave}}/N_{\mathrm{bins}}`$ is larger than the interfacial width and smaller than the total width of the nematic slab. This is not exactly true in an actual configuration due to the bulk fluctuations of $`S(z)`$, but the methods still works well even for small block sizes (see Fig. 2).

## IV Results

The fluctuations of the total director (the direction of the eigenvalue corresponding to the largest eigenvalue of the total ordering tensor) turned out to be so slow that the director hardly changed throughout the run. It always pointed in the $`y`$ direction. This made it convenient to resolve wave-vector components along and perpendicular to the director, without the need to apply any kind of director constraint.

We used the procedure sketched in section III to determine the local deviations $`h(x,y)=z_{\mathrm{int}}(x,y)\overline{z}_{\mathrm{int}}`$ of the local interface position from its average for various block sizes $`B`$. The landscape obtained for the block size $`B=L/8`$ was further analyzed and Fourier transformed. Fourier modes $`h(𝐪)`$ are labeled by $`𝐧=(n_x,n_y)`$, where $`n_x`$ and $`n_y`$ are positive integers with $`𝐪=\frac{2\pi }{L}𝐧`$.

First, we inspect the relaxation times and the correlation times of the Fourier modes. The time evolution of the modes $`𝐧=(0,1)`$, $`(1,0)`$, $`(0,4)`$, and $`(4,0)`$, is shown as a function of time (in MD steps) in Fig. 3. From there we estimate the number of MD steps needed to equilibrate the system. For the slowest mode $`𝐧=(0,1)`$, in the direction parallel to the director, the equilibration process seems to require roughly $`1.0\times 10^6`$ MD steps. Based on this information, we have discarded the initial $`1.2\times 10^6`$ MD steps, and collected results over the following $`2.96\times 10^6`$ MD steps only.

One also has to ensure that the total length of our simulation run significantly exceeds the characteristic time

scale of the slowest capillary wave mode. In order to check this, we have computed the autocorrelation functions of the squared amplitude of the Fourier modes $`h(𝐪)`$

$$C_{hh}(t)=\frac{|h(t)|^2|h(0)|^2|h|^2^2}{|h|^4|h|^2^2}e^{t/\tau }$$
(9)
They are shown in Fig. 4 for the slowest Fourier modes in the direction normal to the director $`𝐧=(1,0)`$ and parallel to the director $`𝐧=(0,1)`$. The correlation time is estimated from the slope of the fitted line as $`\tau 1.0\times 10^5`$ MD steps for the $`(1,0)`$-mode and $`\tau 3.0\times 10^5`$ MD steps for the $`(0,1)`$-mode. This is only by a factor of 10 smaller

than the length of the total simulation run. Hence the statistical error of our results is quite large.

We can now proceed to compare our results with the predictions of the capillary wave theory.

First, we consider the distribution of $`h(x,y)`$. From the free energy functional (1), one derives the exact prediction

$$P(h)=\frac{1}{\sqrt{2\pi s^2}}\mathrm{exp}\left(\frac{h^2}{2s^2}\right),$$
(10)
where

$`s^2`$ $`=`$ $`h^2(x,y)`$ (11)
$`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle 𝑑𝐪|h(𝐪)|^2}={\displaystyle \frac{1}{2\pi \gamma }}\mathrm{ln}\left({\displaystyle \frac{q_{\mathrm{max}}}{q_{\mathrm{min}}}}\right).`$ (12)
The lower and upper cutoffs $`q_{\mathrm{max}}=2\pi /B`$ and $`q_{\mathrm{min}}=2\pi /L`$ come into play, because the integral $`𝑑𝐪|h(𝐪)|^2𝑑q/q`$ diverges as $`q0`$ and $`q\mathrm{}`$. The actual local height distribution $`P(h)`$ obtained from our simulations is plotted for several block sizes in Fig. 5. It can be fitted nicely by a Gaussian distribution (10).

Next, we study the width $`\omega `$ of the average order parameter profile as a function of the block size. According to the capillary wave theory, the capillary wave fluctuations broaden it for large blocks $`B`$ according to

$$\omega ^2=\omega _0^2+\frac{\pi }{2}s^2=\omega _0^2+\frac{k_\mathrm{B}T}{4\gamma }\mathrm{ln}\left(\frac{\widehat{q}_{\mathrm{max}}}{\widehat{q}_{\mathrm{min}}}\right).$$
(13)
Here the cutoff wavevectors are given by $`\widehat{q}_{\mathrm{min}}=2\pi /B`$ and $`\widehat{q}_{\mathrm{max}}=2\pi /a_0`$, where $`a_0`$ is a microscopic length which need not be specified here.

The broadening effect is demonstrated for different block sizes $`B`$ in Fig. 6. The interfacial width can be estimated by fitting order parameter profiles such as those shown in Fig. 6 to tanh-profiles,

$$S(z)=\frac{1}{2}\left(S_\mathrm{N}+S_\mathrm{I}\right)+\frac{1}{2}\left(S_\mathrm{N}S_\mathrm{I}\right)\mathrm{tanh}\left(\frac{zz_{\mathrm{int}}}{\omega }\right),$$
(14)
where the parameters $`S_\mathrm{N}`$ and $`S_\mathrm{I}`$ are the values of the order parameter in the bulk nematic and isotopic phases. (Note that $`S_\mathrm{I}`$ is nonzero due to the finite number of particles in a bin). Fig. 7 plots the squared interfacial width $`\omega ^2`$ as a function of block size $`B`$. One clearly observes the logarithmic increase of $`\omega ^2`$ with block size $`B`$ predicted by eqn. (13). From the slope of the line one can estimate the interfacial tension, $`\gamma =0.016\pm 0.002k_\mathrm{B}T/\sigma _s^2`$.

For comparison, we have also determined the interfacial tension from the anisotropy of the pressure tensor, making use of the relation

$$\gamma =_{\mathrm{}}^{\mathrm{}}𝑑z[P_\mathrm{N}(z)P_\mathrm{T}(z)],$$
(15)
where $`P_\mathrm{N}`$ and $`P_\mathrm{T}`$ are the normal and transverse pressure tensor components. Here, we consider particles with pair interactions $`V_{ij}`$, eqn (3) and systems with two interfaces in the $`(xy)`$-plane. Neglecting finite size effects and interactions between the interfaces, eqn. (15) thus reads

$$\gamma =\frac{1}{4L^2}\underset{i<j}{}r_{ij}^x\frac{V_{ij}}{r_{ij}^x}+r_{ij}^y\frac{V_{ij}}{r_{ij}^y}2r_{ij}^z\frac{V_{ij}}{r_{ij}^z}.$$
(16)
The simulation data yield $`\gamma =0.0093\pm 0.003\sigma _s^2/k_\mathrm{B}T`$. This value agrees with earlier results obtained in smaller systems. It is of the same order, but smaller than the value estimated from the interfacial broadening. The quantitative difference possibly stems from attractive interactions between the two interfaces of the nematic slab. More simulations, in which the thickness of the slab is varied systematically, would be needed to elucidate this point.

Finally, we turn to the analysis of the capillary wave spectrum $`|h(𝐪)|^2`$, which is predicted to be inversely proportional to $`q^2`$ according to eqn (2). We study separately the $`𝐪`$-direction parallel ($`𝐪(0,1)`$) and perpendicular ($`𝐪(1,0)`$) to the director. The inverse of $`|h(𝐪)|^2`$ in these two directions is shown in Fig. 8 as a function of the squared wave vectors. In the long wavelength limit $`𝐪0`$, the spectrum appears to be isotropic, and $`1/|h(𝐪)|^2`$ approaches zero in agreement with the capillary wave theory. The initial slope is consistent with eqn. (2) if one uses $`\gamma =0.016k_\mathrm{B}T/\sigma _s^2`$, and still roughly compatible if one uses $`\gamma =0.0093k_\mathrm{B}T/\sigma _s^2`$. As $`𝐪`$ increases, the capillary wave spectrum becomes anisotropic. As one might expect intuitively, the amplitudes $`|h(𝐪)|^2`$ are larger for capillary waves in the direction perpendicular to the director.

Deviations of capillary wave spectra from the straight line predicted by the simple capillary wave theory (2) are often discussed in terms of higher order terms in the capillary wave Hamiltonian (1), i.e., terms proportional to squares of higher order derivatives of $`h(x,y)`$. For example, including the next-to-leading term leads to a prediction of the form $`1/|h(𝐪)|^2\gamma q^2+\delta q^4`$, where $`\gamma `$ is the interfacial tension and $`\delta `$ is a bending rigidity. One might intuitively expect that the NI-interface has bending stiffness, based on the argument that the elastic interactions should penalize interfacial bending.

However, Fig. 8 indicates that the sign of the “bending energy” is negative. The naive argument sketched above

is thus clearly wrong. On the contrary, our unexpected result suggests that the elastic interactions influence the capillary waves on the largest length scale: The interface is rougher on short length scales than one would expect from the long wavelength fluctuations. This observation is consistent with our speculations in the introduction, that the elastic interactions may be responsible for the suppression of capillary wave fluctuations of the NI-interface in the Onsager-limit of infinitely elongated particles.

We note that this is not the first observation of a negative bending rigidity at a fluid-fluid interface. A similar phenomenon has been predicted theoretically for liquid-vapour interfaces  and verified experimentally by C. Fradin et al . The negative bending rigidity is attributed in their case to long range van der Waals interactions. Indications of a negative bending rigidity at a liquid-liquid interface have also been found in Molecular dynamics simulations of simple liquids by Stecki and Toxvaerd.

To summarize, we have studied the NI-interface by large scale molecular dynamics simulation of a system of repulsive ellipsoidal molecules. We find that the standard capillary wave theory explains most of our results very well. Discrepancies are encountered when looking at the amplitudes of capillary wave modes at large wave vectors, i.e., on small length scales. We find that they are smaller than expected and anisotropic.

###### Acknowledgements.

The simulations have been performed on the CRAY T3E of the HLRZ in Jülich. N.A. received financial support from the German Science Foundation (DFG). MPA acknowledges the support of the Alexander von Humboldt Foundation and the Leverhulme Trust. The parallel MD program used in this work was originally developed by the EPSRC Complex Fluids Consortium.


format: text

text:
# Ballooning Instability in Polar Caps of Accreting Neutron Stars

## 1 Introduction

Despite recent progress on the spreading of accreted matter from a disk over an *unmagnetized* neutron star (Inogamov & Sunyaev, 1999; Popham & Sunyaev, 2000), there has been scarcely any effort directed at understanding how the magnetically funneled accretion stream of an accretion-powered X-ray pulsar spreads over the neutron star’s surface. The presence of a strong magnetic field raises the possibility that the accreted matter is confined to the polar regions of the neutron star. This funneling and confinement of the accreted matter has been proposed (Joss & Li, 1980) as the reason for the apparent stability of thermonuclear burning on strongly magnetized neutron stars. The presence of strong magnetic fields, with $`B10^{12}\mathrm{G}`$, is deduced from cyclotron features in their spectrum (for reviews, see White et al., 1995; Heindl et al., 2000) and from regular torque reversals (Bildsten et al., 1997) that suggest the magnetospheric and co-rotation radii are similar. In neutron stars with weak magnetic fields (as inferred by the absence of pulsations or cyclotron features in the persistent emission), unstable burning is observed as type I X-ray bursts (Lewin et al., 1995). The funneling and confinement of the accreted matter by the magnetic field is believed to increase the local accretion rate to near- or super-Eddington magnitudes, for which the burning is stable to temperature perturbations (see Bildsten, 1998, for a review).

Confinement of the accreted matter requires its overpressure be balanced by the magnetic field. Magnetostatic equilibria of this type were discussed by Hameury et al. (1983) and Brown & Bildsten (1998). In these papers, however, the stability of such equilibria was not addressed. This question is considered here.

If the magnetic field has an approximately potential dipole shape (as can be expected as long as the accreted plasma overpressure is low compared to the magnetic pressure), the force of gravity has a nonvanishing component perpendicular to the magnetic field. Such a configuration might be expected to be susceptible to Kruskal-Schwarzschild (“magnetic Rayleigh-Taylor”) interchange instabilities (Kruskal & Schwarzschild, 1954). Indeed, it is easy to show that for characteristic parameters of polar cap plasmas the stabilizing effect of the magnetic curvature is small compared to the gravitational drive. Nevertheless, the line-tying in the neutron star crust eliminates interchange modes. Instead, ballooning modes can be unstable if the ratio of plasma overpressure to magnetic pressure $`\mathrm{\Delta }\beta `$ is sufficiently large. This instability has been a subject of numerous investigations in the context of laboratory (e.g., Ohkawa & Kerst, 1961; Coppi & Rosenbluth, 1966; Connor et al., 1979; Freidberg, 1987), solar (e.g., Hood, 1986; De Bruyne & Hood, 1989; Strauss & Longcope, 1994), magnetospheric (e.g., McNutt et al., 1987; Hameiri et al., 1991) and astrophysical (e.g., Parker, 1966; Shibata et al., 1989) plasmas. In the present context we find that this instability occurs when $`\mathrm{\Delta }\beta `$ exceeds $`a/h`$ where $`a`$ is the polar cap radius and $`h`$ is the density scale height near the neutron star crust. For such a large overpressure the equilibrium magnetic field is strongly distorted away from the potential field.

The structure of this paper is as follows. We begin by reviewing (§ 2) the conditions at the polar cap of an accreting neutron star. We then formulate, in section 3, the mathematical model which serves as the framework of our linear stability analysis. In section 4 we discuss the energy principle and derive a simple intuitive expression for the magnetohydrodynamic (MHD) potential energy, relevant for plasmas in the neutron star polar caps. We then use this expression to derive a criterion for stability of short-wavelength (in the direction transverse to the magnetic field) ballooning modes and discuss the stability of a plasma confined by an approximately potential, dipole magnetic field (§ 5). We then derive an approximate analytic expression for an equilibrium in which a large overpressure distorts the field away from the potential field (§ 6). Using this equilibrium we find the marginal stability criterion (§ 7), from which we compute the amount of mass that must be accreted in order to drive the instability (§ 8). We conclude with the discussion of our results in Section 9.

## 2 The structure of polar caps

We are concerned with the stability of accreted matter on the polar cap of a strongly magnetized neutron star. Figure 2 provides a schematic of the situation. The collimated infalling matter approaches the neutron star at roughly free-fall velocity $`v_{\mathrm{𝑓𝑓}}c`$. At the near- or super-Eddington local accretion rates, some form of radiative shock is likely to form (Basko & Sunyaev, 1976). The region we are interested in is far below where the infalling material decelerates and becomes incorporated into the neutron star atmosphere. The continual deposition of matter into the upper atmosphere pushes the underlying material deeper, and the downward flow velocity is $`v\dot{m}/\rho c_s`$, where $`\dot{m}`$ is the local accretion rate per unit area. The atmosphere is, to a high accuracy, in hydrostatic equilibrium. On time scales for matter to flow through the thickness of the atmosphere, the field lines co-move with the fluid (Brown & Bildsten, 1998). Under these conditions we are justified in describing the structure of the accretion region with a magnetostatic approximation.

The polar cap is characterized by two length scales: its radius and the density scale height. Various estimates of the polar cap radius $`a`$ (Lamb et al., 1973; Elsner & Lamb, 1977; Arons & Lea, 1980) suggest that $`a10^5\mathrm{cm}`$, and we shall adopt this value throughout the remainder of this paper. The density scale height $`h`$ is much smaller than this value of $`a`$ (see below). With these approximations, the mass and radius of the neutron star enter only through the surface gravity $`g`$, for which we use $`g=1.86\times 10^{14}\mathrm{cm}\mathrm{s}^2`$, which is the Newtonian value for a canonical neutron star of mass $`M=1.4M_{}`$ and a radius $`R=10\mathrm{km}`$.

Numerical solutions (Hameury et al., 1983; Brown & Bildsten, 1998) of the magnetostatic equations show that the greatest distortion of the magnetic field occurs near the boundaries where the field lines are assumed to be tied. In the outermost layers of the neutron star, such a condition occurs at the neutron star crust where the ions form a Coulomb lattice. It is the rigidity of the crust that supports the superincumbent mass of the accreted matter in the polar cap and anchors the magnetic field lines. The electrostatic coupling is parameterized by

$`\mathrm{\Gamma }`$ $`=`$ $`{\displaystyle \frac{Z^2e^2}{k_\mathrm{B}T}}\left({\displaystyle \frac{4\pi n_\mathrm{N}}{3}}\right)^{1/3}`$ (1)
$``$ $`1.1\left({\displaystyle \frac{Z^2}{A^{1/3}}}\right)\left({\displaystyle \frac{\rho }{10^8\mathrm{g}\mathrm{cm}^3}}\right)^{1/3}{\displaystyle \frac{10^8\mathrm{K}}{T}}`$
where $`k_B`$ is the Boltzmann constant, $`\rho `$ is the density, and $`T`$, $`n_\mathrm{N}`$, $`Ze`$, and $`A`$ are, respectively, the temperature, the number density, the charge, and the atomic number of nuclei. For $`\mathrm{\Gamma }1`$, the ions form a liquid. For sufficiently large $`\mathrm{\Gamma }`$, the ions form a crystalline lattice. The composition of the deep ocean and crust is unknown (see Schatz et al., 1999); for definiteness we take the ocean and crust to be composed of pure iron, and assume that solidification occurs for $`\mathrm{\Gamma }>173`$, in accordance with recent Monte-Carlo simulations (Farouki & Hamaguchi, 1993). At this location, the crust of the neutron star can support a shear stress and therefore anchor the field lines. In the absence of a full calculation of the crust and magnetic field, we shall treat the crust \[similarly as in previous works (Hameury et al., 1983; Brown & Bildsten, 1998)\] as a rigid, perfectly conducting surface to which the field lines are anchored (see, however, § 9).

As evident in equation (1), the location of the crust is determined by the temperature of the ocean. We assume that the burning of accreted hydrogen and helium is stable (this assumption can be checked *a posteriori* by finding the depth at which matter spreads). Stablity of the H/He burning requires that the accretion rate be locally super-Eddington (see Bildsten, 1998); for the remainder of this work we shall use the local Eddington rate as our fiducial value. In this case, the temperature in the ocean is set by the nuclear flux released from burning of the accreted fuel to iron peak elements (about $`8\mathrm{MeV}`$ per nucleon) and is $`5\times 10^8\mathrm{K}`$ (Brown & Bildsten, 1998; Schatz et al., 1999).

For the density and temperature of interest, the electrons are degenerate, i.e., $`k_\mathrm{B}T<E_Fm_ec^2`$, where $`E_F=m_ec^2[1+(3\pi ^2n_e)^{2/3}(\mathrm{}/m_ec)^2]^{1/2}`$ is the electron Fermi energy, and $`m_e`$ and $`n_e`$ are the electron mass and density, respectively. The principal contribution to pressure $`p`$ is the electron pressure $`p_e`$ and therefore to a good approximation $`pp_en_eE_\mathrm{F}/4\rho ^\gamma `$ where $`\gamma =4/3`$ and the proportionality constant is only a function of the ion charge-to-mass ratio $`Z/A`$. For a plasma of uniform composition this constant is independent of both time and space. We shall use this approximation in the next section.

Using this approximate equation of state, we write the density scale height $`h=\gamma p/\rho g`$ at the top of the crust in terms of $`\mathrm{\Gamma }`$,

$`h`$ $``$ $`\gamma \left({\displaystyle \frac{9\pi }{256}}\right)^{1/3}{\displaystyle \frac{1}{Z^{2/3}A}}{\displaystyle \frac{N_\mathrm{A}k_\mathrm{B}T}{g}}{\displaystyle \frac{\mathrm{\Gamma }}{\alpha _\mathrm{F}}}`$ (2)
$`=`$ $`6.90\times 10^3\mathrm{cm}`$
$`\times \left({\displaystyle \frac{26}{Z}}\right)^{2/3}\left({\displaystyle \frac{56}{A}}\right)\left({\displaystyle \frac{T}{5\times 10^8\mathrm{K}}}\right)\left({\displaystyle \frac{\mathrm{\Gamma }}{173}}\right),`$
where $`\gamma =4/3`$, $`N_\mathrm{A}=6.02\times 10^{23}\mathrm{g}^1`$ is Avogadro’s constant and $`\alpha _\mathrm{F}=e^2/\mathrm{}c`$.

## 3 Model

We model the polar cap as the region of the neutron star atmosphere and ocean that is bounded by the same magnetic flux surface within which accretion occurs. This region is characterized by a somewhat higher pressure than in the surrounding matter, due to the weight of the accreted material. The plasma in the polar cap and in the surrounding region will be treated as a perfectly conducting fluid of uniform composition, bounded by a perfectly conducting, rigid surface (the top of the crust) and immersed in a magnetic field. The magnetic field lines are assumed tied to the crust; the intersection point of the field line with the top of the crust will be referred to as the footpoint. Since the density scale height of the plasma in the polar cap is much smaller than the stellar radius we shall consider the field lines as semi-infinite and having only one footpoint. In the equilibrium, the plasma and the immersing magnetic field are assumed to be axisymmetric with the axis of symmetry coinciding with the magnetic axis; the magnetic field is assumed to be untwisted.

We describe the plasma by the single-fluid MHD theory. In the presence of the gravitational potential $`U`$ and magnetic field $`𝐁`$ the momentum balance equation for a plasma of density $`\rho `$ and pressure $`p`$, moving with velocity $`𝐯`$, is

$$\rho \frac{d𝐯}{dt}=\frac{1}{c}𝐉\times 𝐁p\rho U.$$
(3)
While Ohm’s law is modified by gravity, the induction equation is unaffected (cf. Bernstein et al. (1958)), so that

$$\frac{𝐁}{t}=\times 𝐯\times 𝐁.$$
(4)
The above equations are supplemented by Ampère’s law,

$$𝐉=\frac{c}{4\pi }\times 𝐁,$$
(5)
the continuity equation,

$$\frac{\rho }{dt}+(\rho 𝐯)=0,$$
(6)
and an equation of state which we assume to be adiabatic:

$$\frac{d}{dt}p\rho ^\gamma =0.$$
(7)
As mentioned in §1, the primary contribution to the plasma pressure in the region of interest is due to degenerate, relativistic electrons ($`\gamma =4/3`$). Therefore, in region of uniform composition, $`p\rho ^\gamma `$ is, to a high accuracy, constant in space. For the assumptions in § 2, $`ϵd\mathrm{ln}\rho p^{1/\gamma }/d\mathrm{ln}p3/8p_{\mathrm{ion}}/p_\mathrm{e}=(3/2Z)(k_\mathrm{B}T/E_\mathrm{F})4\times 10^4`$. Thus

$$p\rho ^\gamma =\mathrm{const}.,$$
(8)
to a very good approximation.

The force balance in a static equilibrium is given by

$$\frac{1}{c}𝐉\times 𝐁=p+\rho U,$$
(9)
from which it follows that

$$p+\rho U=\rho F,$$
(10)
where the function $`F`$ is a function constant on the field lines

$$𝐁F=0.$$
(11)
By setting the gravitational potential to zero at the top of the crust, $`F`$ is determined by the sound speed $`c_s=\sqrt{\gamma p/\rho }`$ there, $`F=c_{s0}^2/(\gamma 1)`$.

An untwisted, axisymmetric equilibrium magnetic field can be represented in the form

$$𝐁=\psi \times \theta $$
(12)
where $`\psi `$ is the flux function and $`\theta `$ is the toroidal angle. It is straightforward to show with the aid of Ampère law that in such a geometry the current is perpendicular to the magnetic field: $`𝐉𝐁=0`$. Eq. (11) implies that $`F`$ is a function of the flux function only: $`F=F(\psi )`$.

Upon introducing cylindrical coordinates ($`r,\theta ,z)`$ with the $`z`$-direction along the axis of symmetry, the normal (i.e., parallel to $`\psi `$) component of Eq. (10) yields the Grad-Shafranov-like equation

$$r\frac{}{r}\frac{1}{r}\frac{\psi }{r}+\frac{^2\psi }{z^2}=4\pi r^2\rho (\psi ,z)\frac{dF}{d\psi }.$$
(13)
with $`\rho `$ being determined by equation (10).

In the limit $`aR`$ the curvature of the stellar surface can be neglected on the scale of the accretion region. In this approximation $`U=gz`$, so that

$$\rho (\psi ,z)=\rho _0(\psi )\left[1\alpha \frac{z}{h(\psi )}\right]^{1/\alpha }$$
(14)
where $`z=0`$ corresponds to the top of the crust, $`\alpha \gamma 1`$, $`\rho _0`$ is the plasma density at $`z=0`$, $`h(d\mathrm{ln}\rho /dz)^1|_{z=0}=\alpha F/g`$ is the density scale height at the footpoints and $`g`$ is the gravitational acceleration. Note that $`\rho _0`$ and $`h`$ are functions of $`\psi `$.

## 4 Energy Principle

We now address the question of the linear stability of the polar cap plasma in the approximation described in the previous section. We base our analysis on the MHD energy principle (Bernstein et al. 1958). In the present section we discuss the energy principle as applicable to the polar cap plasmas, generalizing the “intuitive” form of the energy principle (Furth et al., 1965; Greene & Johnson, 1968) to include the effects of a strong gravitational field. In our discussion we follow the approach of Freidberg (1987).

We start by considering the linearized momentum balance equation for a static equilibrium

$$\rho \frac{^2𝝃}{t^2}=𝐅(𝝃)\frac{1}{c}\stackrel{~}{𝐉}\times 𝐁+\frac{1}{c}𝐉\times \stackrel{~}{𝐁}\stackrel{~}{p}\stackrel{~}{\rho }U$$
(15)
where $`𝝃`$ is the plasma displacement and the gravitational potential is unperturbed. Here and in the following the tilde denotes perturbed quantities while symbols without it stand for the equilibrium quantities. Perturbed current $`\stackrel{~}{𝐉}`$, magnetic field $`\stackrel{~}{𝐁}`$, density $`\stackrel{~}{\rho }`$, and pressure $`\stackrel{~}{p}`$ are respectively given by Ampère’s law, the induction equation, the continuity equation, and the equation of state:

$`\stackrel{~}{𝐉}={\displaystyle \frac{c}{4\pi }}\times \stackrel{~}{𝐁}`$ (16)
$`\stackrel{~}{𝐁}=\times (𝝃\times 𝐁)`$ (17)
$`\stackrel{~}{\rho }=(\rho 𝝃)`$ (18)
$`\stackrel{~}{p}=\gamma p𝝃𝝃p.`$ (19)
We shall assume that the plasma boundary (the crust) is rigid and perfectly conducting, so that $`𝝃`$ and $`\stackrel{~}{𝐁}`$ vanish there.

The second-order potential energy is given by (Bernstein et al. 1958)

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle 𝑑\tau 𝝃𝐅(𝝃)}`$
$`={\displaystyle \frac{1}{2}}{\displaystyle }d\tau 𝝃\{{\displaystyle \frac{1}{c}}\stackrel{~}{𝐉}\times 𝐁+{\displaystyle \frac{1}{c}}𝐉\times \stackrel{~}{𝐁}`$
$`+(\gamma p𝝃+𝝃p)+(\rho 𝝃)U\}`$ (20)
where $`d\tau `$ denotes integration over the plasma volume. We now transform the above expression for the potential energy into a form that allows a physical interpretation. This is analogous to what Furth et al. (1965) and Greene & Johnson (1968) have done in the absence of gravitational field. In order to derive it, we follow the standard approach \[e.g., Freidberg (1987)\], with appropriate modifications for the presence of a strong gravitational field.

Integrating by parts and dropping surface terms (which vanish because of the assumed boundary conditions), and exploiting the equilibrium relation given by Eq. (10), one obtains, after some straightforward algebra,

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle }d\tau \{\gamma p(𝝃2𝜿_g𝝃)^2+ϵp(\xi \mathrm{ln}p)^2+{\displaystyle \frac{\stackrel{~}{𝐁}^2}{4\pi }}`$
$`𝝃[{\displaystyle \frac{1}{c}}𝐉\times \stackrel{~}{𝐁}+2𝜿_g(\rho 𝝃F)+(\rho 𝝃F)]\}.`$ (21)
where $`𝜿_g(\rho /\gamma p)U/2`$. The the second term on the RHS, proportional to the Schwarzschild discriminant, can be neglected for highly degenerate plasmas in strong magnetic fields because of the smallness of parameter $`ϵ=d\mathrm{ln}\rho p^{1/\gamma }/d\mathrm{ln}p`$. We shall omit this term in the subsequent analysis (see, however, discussion of finite $`ϵ`$ effects in §9).

The last two terms in the square bracket on the RHS of Eq. (4) can be expanded and recombined to yield an alternative form of the energy integral

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle }d\tau \{{\displaystyle \frac{\stackrel{~}{𝐁}^2}{4\pi }}+\gamma p(𝝃2𝜿_g𝝃)^2`$
$`𝝃[{\displaystyle \frac{1}{c}}𝐉\times \stackrel{~}{𝐁}+{\displaystyle \frac{\rho ^2}{\gamma p}}(𝝃F)\rho F+\rho (𝝃F)]\}.`$ (22)
Since $`𝐁F=0`$ and

$$\frac{1}{c}𝐁𝐉\times \stackrel{~}{𝐁}=\rho 𝐁(𝝃F),$$
(23)
it is clear that the parallel (to $`𝐁`$) component of the expression in the square bracket in Eq. (4), and hence also in Eq. (4), vanishes. Thus the energy integral can be expressed in the form

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle }d\tau \{{\displaystyle \frac{\stackrel{~}{𝐁}^2}{4\pi }}+\gamma p(𝝃2𝜿_g𝝃)^2`$
$`𝝃_{}[{\displaystyle \frac{1}{c}}𝐉\times \stackrel{~}{𝐁}+2(\rho 𝝃_{}F)𝜿_g+(\rho 𝝃_{}F)]\}.`$ (24)
Next, splitting up $`\stackrel{~}{𝐁}`$ into transverse and parallel components, $`\stackrel{~}{𝐁}=\stackrel{~}{𝐁}_{}+\stackrel{~}{B}_{}\widehat{𝐛}`$ ($`\widehat{𝐛}𝐁/B`$), we find

$`𝝃_{}𝐉\times \stackrel{~}{𝐁}=J_{}𝝃_{}\widehat{𝐛}\times \stackrel{~}{𝐁}+\stackrel{~}{B}_{}𝝃_{}𝐉\times \widehat{𝐛}`$ (25)
$`\stackrel{~}{B}_{}=B(𝝃_{}+2𝜿_c𝝃_{})+{\displaystyle \frac{4\pi \rho }{B}}𝝃_{}F,`$ (26)
where $`𝜿_c=\widehat{𝐛}\widehat{𝐛}`$ is the field line curvature. With these expressions the potential energy can be expressed in the following, “intuitive” form

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle }d\tau \{{\displaystyle \frac{\stackrel{~}{𝐁}_{}^2}{4\pi }}+{\displaystyle \frac{𝐁^2}{4\pi }}(𝝃_{}+2𝜿_c𝝃_{})^2`$
$`+\gamma p(𝝃2𝜿_g𝝃)^2{\displaystyle \frac{1}{c}}J_{}𝝃_{}\times \widehat{𝐛}\stackrel{~}{𝐁}`$
$`2𝜿𝝃_{}\rho (𝝃_{}F)\}.`$ (27)
where $`𝜿𝜿_c+𝜿_g`$. The first term is the energy corresponding to field line bending, the second is the field compression energy, the third is the energy of plasma compression and buoyancy, the fourth one is responsible for the kink instability, and the last one drives the curvature and gravitational instabilities. The above expression is a generalization of the result of Furth et al. (1965) and Greene & Johnson (1968) to the case of a degenerate plasma in a strong gravitational field. The above expression reduces to these earlier results if the density scale height is much longer than the plasma length, i.e., if $`Uc_s^2`$.

Thus far we have not taken into account the geometry of polar caps and the derived energy principle (eq. ) applies to an arbitrary configuration of plasmas satisfying the equation of state (7). We shall now perform a reduction of the energy principle as applicable to the configuration of our interest.

First we note that the parallel displacement $`\xi _{}`$ appears only in the nonnegative compressive term. Therefore this term can be minimized, in a standard manner, by considering the perturbation of the compressive energy due to a perturbation of $`\xi _{}`$. The minimum is found for the displacement satisfying the Euler-Lagrange equation which can be expressed in the form

$$𝐁\left[p^{1\frac{1}{\gamma }}\left(𝝃2𝜿_g𝝃\right)\right]=0$$
(28)
Integrating along the field line and imposing the line-tying boundary condition ($`\xi _{}[0]=0`$) yields

$$\xi _{}=\frac{B}{p^{\frac{1}{\gamma }}}_0^l\frac{dl^{}}{B}p^{\frac{1}{\gamma }}\left[p^{\frac{1}{\gamma }1}A+𝝃_{}2𝜿_g𝝃_{}\right]$$
(29)
where $`l`$ is the distance along the field line. The integration constant $`A`$ must be determined from the integrability condition. Note that $`p,\rho 0`$ as $`zz_0h/(\gamma 1)`$ (cf. equations  and ). Therefore, for the kinetic energy of parallel motion to be finite, (i.e., $`𝑑\tau \rho \xi _{}^2<\mathrm{}`$), this requires that

$$A=\frac{_0^{l_0}\frac{dl^{}}{B}p^{\frac{1}{\gamma }}\left(𝝃_{}2𝜿_g𝝃_{}\right)}{_0^{l_0}\frac{dl^{}}{B}p^{\frac{2}{\gamma }1}}$$
(30)
where $`l_0l(z_0)`$. This implies, in particular, that $`\xi _{}(z_0)=0`$. Since

$$𝝃2𝜿_g𝝃=p^{\frac{1}{\gamma }1}A$$
(31)
and $`A0`$, the compressional energy cannot be minimized to zero. This is analogous to that which occurs in flux tubes line-tied at both ends (cf. Hameiri et al 1991) and is, in fact, caused by the effective line-tying in the low density region.

Setting $`J_{}=0`$, as appropriate for an untwisted axisymmetric field, the potential energy becomes

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle }d\tau \{{\displaystyle \frac{\stackrel{~}{𝐁}_{}^2}{4\pi }}+{\displaystyle \frac{𝐁^2}{4\pi }}(𝝃_{}+2\kappa _c𝝃_{})^2`$
$`+\gamma p^{\frac{2}{\gamma }1}\left[{\displaystyle \frac{_0^{l_0}\frac{dl^{}}{B}p^{\frac{1}{\gamma }}\left(𝝃_{}2𝜿_g𝝃_{}\right)}{_0^{l_0}\frac{dl^{}}{B}p^{\frac{2}{\gamma }1}}}\right]^2`$
$`2\rho (𝝃_{}F)𝜿𝝃_{}\}.`$ (32)
## 5 Ballooning instability

We now limit our attention to the short wavelength (in the direction perpendicular to the magnetic field) modes. Following the method of Freidberg (1987), we perform the reduction of the intuitive form of the energy principle (eq. ).

We represent the displacement in the eikonal approximation

$$𝝃_{}=𝜼_{}e^{iS}+c.c.$$
(33)
where $`c.c.`$ stands for complex conjugate; $`𝜼_{}`$ is assumed to be slowly varying while the eikonal $`S`$ is assumed to vary rapidly, i.e.,

$$𝐤_{}=S$$
(34)
is so large that $`k_{}a1`$, where $`a`$ is the equilibrium length scale.

Note that $`𝐁S=0`$. Consequently

$$\stackrel{~}{𝐁}_{}=e^{iS}[\times (𝜼_{}\times 𝐁)]_{}+c.c.$$
(35)
In the limit of $`k_{}\mathrm{}`$, the dominant terms in $`\delta W`$ (Eq.) are those proportional to $`(𝝃_{})^2|𝐤_{}𝜼_{}|^2`$ (second and third term on the RHS of eq. ). Therefore the minimization of the potential energy, to lowest order in $`1/k_{}a`$, requires that $`𝐤_{}𝜼_{}=0`$, which implies that for the zeroth-order displacement

$$𝜼_{}\frac{X}{B}\widehat{𝐛}\times 𝐤_{}.$$
(36)
where $`X`$ is a scalar function. To the next order, the first-order displacement can be chosen so that the sum of two positive terms involving $`𝝃_{}`$ (second and third term on the RHS of eq. ) is minimized. Since we are interested in the limit of $`\beta 8\pi p/B^21`$, we can instead minimize only the second of these terms. A simple choice of a minimizing displacement is

$$𝝃_{}=2𝜿_g𝝃_{}.$$
(37)
For this choice, $`\xi _{}=0`$ as well as the third term on the RHS of eq. (4) vanish so that

$`\delta W={\displaystyle \frac{1}{2}}{\displaystyle }d\tau \{{\displaystyle \frac{\stackrel{~}{𝐁}_{}^2}{4\pi }}+{\displaystyle \frac{𝐁^2}{\pi }}(𝜿𝝃_{})^2`$
$`2\rho (𝝃_{}F)(𝜿𝝃_{})\}.`$ (38)
From equation (36) it follows that

$$\stackrel{~}{𝐁}_{}=e^{iS}(X\times \widehat{𝐛})_{}+c.c.=(\widehat{𝐛}X)\widehat{𝐛}\times 𝐤_{}+c.c.$$
(39)
Consequently, one finds from equation (5), exploiting the hermiticity of the operator $`𝐅`$ (Bernstein et al. 1958), that

$`\delta W={\displaystyle \frac{1}{4\pi }}{\displaystyle }d\tau \{k_{}^2|\widehat{𝐛}X|^2+4(𝜿\times \widehat{𝐛}𝐤_{})^2|X|^2`$
$`{\displaystyle \frac{8\pi \rho }{B^2}}(𝜿\times \widehat{𝐛}𝐤_{})(F\times \widehat{𝐛}𝐤_{})|X|^2\}.`$ (40)
For an axisymmetric system, $`S=m\theta +\widehat{S}(\psi ,l)`$ where $`l`$ is the distance along the field line. In the absence of a toroidal field $`B_\theta `$, $`𝐁S=0`$ implies that $`\widehat{S}=\widehat{S}(\psi )`$. Thus

$$𝐤_{}=k_\theta \widehat{𝜽}+k_n\widehat{𝐧}$$
(41)
where $`k_\theta =m/r`$ and $`\widehat{𝐧}=\psi /|\psi |=\psi /rB`$ is the unit vector normal to the flux surface. Noting that

$`\widehat{𝐛}\times 𝐤_{}=k_\theta \widehat{𝐧}+k_n\widehat{𝜽}`$ (42)
$`\widehat{𝐛}\times 𝐤_{}F=k_\theta rB{\displaystyle \frac{dF}{d\psi }},`$ (43)
the potential energy can be put in form

$`\delta W={\displaystyle \frac{1}{4\pi }}{\displaystyle }d\tau \left\{(k_\theta ^2+k_n^2)\right|{\displaystyle \frac{X}{l}}|^2`$
$`+4k_\theta ^2\kappa _n(\kappa _n{\displaystyle \frac{2\pi r\rho F^{}}{B}})|X|^2\}`$ (44)
($`F^{}=dF/d\psi `$). Minimization with respect to $`k_n`$ is achieved by setting $`k_n=0`$.

The integrand in Eq. (5) involves derivatives of $`X`$ only along the magnetic field but not derivatives of $`\psi `$. Hence $`\psi `$ can be treated as a parameter and the potential energy can be expressed in the form

$`\delta W={\displaystyle 𝑑\psi \delta W(\psi )},`$ (45)
where

$`\delta W(\psi )={\displaystyle \frac{m^2}{4\pi }}{\displaystyle }{\displaystyle \frac{dl}{r^2B}}\left\{\right|{\displaystyle \frac{X}{l}}|^2`$
$`+4\kappa _n(\kappa _n{\displaystyle \frac{2\pi r\rho F^{}}{B}})|X|^2\}.`$ (46)
Equation (45) implies that $`\delta W(\psi )<0`$ for some $`\psi `$ is both a necessary and sufficient condition for instability (Bernstein et al., 1958).

Let us first address the question of stability of a plasma confined by a stellar dipole magnetic field. We approximate $`F^{}`$ by $`\mathrm{\Delta }F/\psi _a`$ where $`\psi _a=a^2B_0/2`$ and $`\mathrm{\Delta }F=\mathrm{\Delta }p_0/\rho _0`$, with $`\mathrm{\Delta }p_0`$ and $`\rho _0`$ being, respectively, the overpressure and density at the top of the crust. Since $`\kappa _{gn}=(U/n)/2c_{s0}^2\widehat{b}_r/2h`$, $`\kappa _{cn}\widehat{b}_r/R`$, $`\widehat{b}_rr/R`$ (for a dipole field) and the density scale height $`h=c_{s}^{2}{}_{0}{}^{}/g`$ is much smaller than the neutron star radius $`R`$, the second term on the RHS of Eq. (5) is destabilizing, if $`\mathrm{\Delta }\beta a^2/hR`$.

For the interchange modes $`X/l=0`$. Because of line tying, however, $`X(0)=0`$ with $`l=0`$ being chosen to correspond to the top of the crust. Expanding into Taylor series, we can approximate $`Xl`$ for $`lh`$. From the structure of the integral in Eq. (5) we expect that the most unstable displacement vanishes for $`lh`$.

We shall now estimate the integral in Eq. (5). To lowest order in $`a/R`$, where $`a`$ is the radius of the polar cap, $`U=c_s^2l/h`$. We then find that $`\delta W(\psi _a)`$ becomes negative only if

$`\mathrm{\Delta }\beta {\displaystyle \frac{R}{h}}`$ (47)
where $`\mathrm{\Delta }\beta =8\pi \mathrm{\Delta }p_0/B_0^2`$. The above equation implies that for a low overpressure ($`\mathrm{\Delta }\beta 1`$), the plasma confined by the dipole magnetic field of a neutron star is stable with respect to ballooning modes. Equation (47) indicates that the plasma may become unstable at high $`\mathrm{\Delta }\beta 1`$. For such a high plasma overpressure, however, the equilibrium magnetic field is distorted away from the dipole field and has to be computed self-consistently. We consider this case in the next section.

## 6 A simple high-$`\mathrm{\Delta }\beta `$ equilibrium

In the limit $`aR`$ the dipole magnetic field in the polar cap is well approximated by a uniform field. Let us therefore assume that prior to the loading of the accretion column with plasma, the magnetic field $`𝐁=B_0\widehat{z}`$. We shall assume that at $`z=0`$ the flux at the is unaffected by the accretion, i.e., $`\psi (r,0)=r^2B_0/2`$. The loading is assumed to occur on flux surfaces with $`\psi \psi _aa^2B_0/2`$.

Until now, the overpressure profile at the base of the accretion region, which determines the function $`F`$, has not been specified. As a simple example we shall consider the following form of $`F`$,

$$F(\psi )=F_{ext}+\mathrm{\Delta }F\left[1\left(\frac{\psi }{\psi _a}\right)^2\right]^{1+\epsilon }$$
(48)
for $`\psi \psi _aa^2B_0/2`$, with $`\epsilon 1`$; $`F`$ is a constant $`F_{ext}`$ for $`\psi >\psi _a`$. The gradient of $`F`$ is approximately linear for $`\psi <\psi _a(1\epsilon )`$ and goes smoothly to zero at $`\psi =\psi _a`$. In the following we shall restrict our attention to the case of an overpressure $`\mathrm{\Delta }p_0`$ that is small compared to the total pressure $`p_0`$, so that $`\mathrm{\Delta }F=\mathrm{\Delta }p_0/\rho _0`$. In the same approximation, the scale height $`h`$ is, to lowest order, independent of $`\psi `$.

For such a pressure profile, the Grad-Shafranov equation (Eq. ) inside the accretion column \[$`\psi <\psi _a(1\epsilon )`$\] becomes approximately linear,

$$r\frac{}{r}\left(\frac{1}{r}\frac{\psi }{r}\right)+\frac{^2\psi }{z^2}=\frac{4\mathrm{\Delta }\beta }{a^2}\left(\frac{r}{a}\right)^2\left(1\alpha \frac{z}{h}\right)^{\frac{1}{\alpha }}\psi .$$
(49)
We anticipate that because $`\mathrm{\Delta }\beta 1`$ and $`ha`$, the radial excursion of the flux surfaces $`\mathrm{\Delta }r=𝑑l\widehat{b}_r`$ ($`h`$) will be small compared to $`r`$. We can therefore neglect the $`r`$-dependence on the RHS of Eq. (49) and treat $`\chi r/a`$ as a parameter. With these approximations, we can then solve Eq. (49) by separation of variables,

$$\psi (r,z)=\psi _0\chi ^2f(\zeta )|_{\chi =r/a,\zeta =z/h},$$
(50)
where $`f`$ satisfies the equation

$$\frac{d^2f}{d\zeta ^2}=\lambda ^2\chi ^2(1\alpha \zeta )^{\frac{1}{\alpha }}f$$
(51)
($`\lambda 2h\sqrt{\mathrm{\Delta }\beta }/a`$). We solve this equation subject to the boundary conditions $`f(0)=1`$ and $`f^{}(\alpha ^1)=0`$. The first of these conditions expresses normalization; the second follows from the requirement that the field be vertical in the region of vanishing overpressure. By changing variables, we can transform the above equation into a modified Bessel equation which has the solution

$$f(\zeta )=u^\nu (\zeta )\left\{𝒜I_\nu [u(\zeta )]+I_\nu [u(\zeta )]\right\},$$
(52)
where $`u(\zeta )=2\lambda \chi (1\alpha \zeta )^{\frac{1}{2\nu }}/(1+2\alpha )`$, $`\nu =\alpha /(1+2\alpha )=(\gamma 1)/(2\gamma 1)`$, $`I_{\pm \nu }`$ are the modified Bessel functions, and $`𝒜`$ and $``$ are the integration constants. From the condition $`f^{}(\alpha ^1)=0`$ it follows that $`=0`$; the boundary condition $`f(0)=1`$ then implies that

$$f(\zeta )=\left(1\alpha \zeta \right)^{\frac{1}{2}}\frac{I_\nu \left[\frac{2\lambda \chi }{1+2\alpha }\left(1\alpha \zeta \right)^{\frac{1}{2\nu }}\right]}{I_\nu \left(\frac{2\lambda \chi }{1+2\alpha }\right)}.$$
(53)
An example of flux surfaces given by equations (50) and (53) is shown in Figure 2. We see that even for relatively large $`\lambda `$, the radial excursion $`\mathrm{\Delta }r`$ of flux surfaces is small compared to $`r`$, in accord with our assumptions.

Let us now consider the case of $`a/h\mathrm{\Delta }\beta a^2/h^2`$ which, as we shall see, is valid at the onset of instability. In this case, $`\lambda 1`$ and we can expand the Bessel functions into a power series. In this approximation we find that

$$\psi (r,z)=\psi _0\frac{r^2}{a^2}\left[1\lambda ^2\frac{r^2}{a^2}\frac{1\left(1\alpha z/h\right)^{2+\frac{1}{\alpha }}}{\left(1+\alpha \right)\left(1+2\alpha \right)}\right].$$
(54)
With the aid of the above expression we can now demonstrate the self-consistency of our earlier assumptions. We find that the first term on the LHS of Eq. (49) is smaller than the second one as long as $`\rho /\rho _0>(2h/a)^2`$. Thus for $`ha`$ our solution is indeed self-consistent in the region of plasma where the most unstable displacement is expected to be localized (see §7).

From (54) the magnetic field components are

$`B_r`$ $``$ $`\mu \left({\displaystyle \frac{\rho }{\rho _0}}\right)^\gamma B_0`$ (55)
$`B_z`$ $``$ $`B_0`$ (56)
where $`\mu =2\mathrm{\Delta }\beta hx^3/\gamma a`$. From the above it follows that the angle $`\varphi `$ between the magnetic and gravitational fields is given by

$`\mathrm{tan}\varphi `$ $``$ $`\mu \left(1\alpha {\displaystyle \frac{z}{h}}\right)^{\gamma /\alpha }.`$ (57)
From equation (57) one sees explicitly that the field line distortion becomes large ($`\mathrm{tan}\varphi 1`$) when $`\mu 1`$, i.e., $`\mathrm{\Delta }\beta a/h`$, as originally observed by Hameury et al. (1983). This is, in fact, a simple consequence of balancing overpressure $`\mathrm{\Delta }p/a`$ with magnetic tension $`\kappa B^2/4\pi `$ since $`\kappa h1`$ when the field distortion becomes large.

## 7 Criterion for the onset of instability

We now use the equilibrium found in the previous section to evaluate the energy integral. Let us first change the integration variable in Eq. (5) from $`l`$ to $`z`$, such that $`dl=dz/\widehat{b}_z`$. For the nearly parabolic pressure profile considered earlier,

$`\delta W(\psi )={\displaystyle \frac{m^2}{8\pi }}{\displaystyle }{\displaystyle \frac{dz}{r^2B_z}}\left\{\widehat{b}_z^2\right|{\displaystyle \frac{X}{z}}|^2`$
$`+4\kappa _n(\kappa _n+{\displaystyle \frac{\gamma \mu }{2h}}{\displaystyle \frac{B_0}{B}}{\displaystyle \frac{\rho }{\rho _0}})|X|^2\}.`$ (58)
In the approximation of Eq. (54), we find that to lowest order in $`\lambda `$

$$\kappa _n(1+2\gamma \widehat{b}_z^2)\left(\frac{\rho _0}{\rho }\right)^\alpha \frac{\widehat{b}_r}{2h}.$$
(59)
Changing the integration variable to $`s=\mathrm{tan}\varphi `$, as given by equation (57), we then find, with the aid of (55), (56), (57) and (59),

$`\delta W(\psi )={\displaystyle \frac{\mu ^{11/\gamma }m^2}{8\pi r^2\gamma hB_0}}{\displaystyle _0^\mu }ds{\displaystyle \frac{s^{1/\gamma }}{1+s^2}}\left\{\gamma ^2\right|{\displaystyle \frac{X}{s}}|^2`$
$`(\gamma 1){\displaystyle \frac{(s^2\mu _0^2)(s^2+2\gamma +1)}{(s^2+1)^2}}\left|X|^2\right\},`$ (60)
where $`\mu _0=\sqrt{(\gamma +1)/(\gamma 1)}`$. We observe that the integrand in equation (7) is positive for $`s\mu _0`$. Consequently, the most unstable displacement is expected to vanish for $`s\mu _0`$ which implies the instability is localized to the region where

$`\rho \left({\displaystyle \frac{\mu _0}{\mu }}\right)^{\frac{1}{\gamma }}\rho _0.`$ (61)
Since the line-tying in the crust requires that $`X(s=\mu )=0`$, we therefore choose as a test function

$$X(s)=X_0\mathrm{sin}\left(\pi \frac{s\mu _0}{\mu \mu _0}\right)$$

For this test function and $`\gamma =4/3`$, as appropriate for a plasma in the vicinity of the crust, we find that $`\delta W(\psi )<0`$ for $`\mu >11.67`$ for which it follows that

$`\mathrm{\Delta }\beta >7.8{\displaystyle \frac{a}{h}}.`$ (62)
is a sufficient condition for instability. For $`\mu =11.67`$, equation (61) implies that the instability is localized to within 0.93 scale height from the crust.

The instability criterion (62) is obtained for the equilibrium given by equation (54), derived in the approximation of $`\mathrm{\Delta }\beta (a/h)^2`$. Thus self-consistency implies that the validity of (62) is limited to $`h/a0.13`$. This condition is satisfied, albeit marginally, for the parameters discussed in §2. The instability threshold for more general equilibria, with higher values of $`h/a`$, is at present unknown.

## 8 Mass of confined accreted matter

Assuming that the crust is perfectly rigid, the amount of mass needed to drive the instability is, from inequality (62),

$`\mathrm{\Delta }M`$ $`=`$ $`\pi a^2{\displaystyle \frac{B^2}{8\pi }}{\displaystyle \frac{\mathrm{\Delta }\beta }{g}}`$ (63)
$`=`$ $`0.73B^2a^3{\displaystyle \frac{\rho }{p}}.`$
Substituting Eq. (2) in Eq. (63), we have

$`\mathrm{\Delta }M`$ $`=`$ $`3.8\times 10^{13}M_{}\left({\displaystyle \frac{B}{10^{12}\mathrm{G}}}\right)^2\left({\displaystyle \frac{a}{10^5\mathrm{cm}}}\right)^3\times `$ (64)
$`\left({\displaystyle \frac{Z}{26}}\right)^{2/3}\left({\displaystyle \frac{A}{56}}\right)\left({\displaystyle \frac{5\times 10^8\mathrm{K}}{T}}\right)\left({\displaystyle \frac{173}{\mathrm{\Gamma }}}\right).`$
At a given accretion rate $`\dot{M}`$, the time necessary to accumulate this mass is $`\mathrm{\Delta }M/\dot{M}=33\mathrm{hr}(10^{10}M_{}\mathrm{yr}^1/\dot{M})`$. At present we cannot say what, if any, are the observable effects of the instability, but the timescale to accumulate a critical mass $`\mathrm{\Delta }M`$ is well-suited for monitoring.

## 9 Discussion

Using the MHD energy principle, we have shown that the accreted material on the polar cap of a strongly magnetized neutron star is subject to the ballooning instability when the overpressure in the vicinity of the crust exceeds the magnetic pressure by the factor $`7.8a/h`$, where $`a`$ is the horizontal pressure gradient length scale (identified with the polar cap radius) and $`h`$ is the density scale height. For typical polar cap parameters, the critical overpressure is $`10^4`$ of the total plasma pressure. The instability is localized to within a density scale height of the field line footpoint (at the top of the crust) which implies that plasma is perturbed at densities $`(\rho /\rho _0)>0.33`$; for $`\rho _0`$ found from equation (1) for a model atmosphere (cf. Brown & Bildsten 1999), with a temperature of $`7.6\times 10^8\mathrm{K}`$, assuming a locally Eddington accretion rate. This implies $`\rho >10^{10}\mathrm{g}\mathrm{cm}^3`$. This density is higher than where the accreted hydrogen and helium ignite which suggests that the ballooning instability near the threshold for onset will not affect the nuclear burning. It is not clear, however, whether the instability can degrade the confinement enough to prevent a further accumulation of mass well above the instability threshold.

In our stability analysis we have omitted the contribution of the Schwarzschild term to the MHD potential energy arguing that it is small in highly degenerate plasmas of uniform composition immersed in strong magnetic fields. To test the validity of this approximation we have performed a calculation along the same lines as in §4 and §5, including the buoyancy contribution, for $`ϵ=4\times 10^4`$, for a realistic atmosphere (cf. Brown & Bildsten 1999). We found that this inclusion increases the numerical factor on the RHS of inequality (62) by 17% for $`B=10^{12}`$ G and by less than 1% for $`B5\times 10^{12}`$. Thus the effect of buoyancy is indeed weak for magnetic fields of interest, as long as $`ϵ`$ is as small as assumed.

The latter assumption, as well as that of composition uniformity, is violated in the electron-capture layers. The location of these layers is not precisely known. Haensel & Zdunik (1990) found that electron capture layers closest to the top of the crust occur at densities $`10^{10}`$ g cm<sup>-3</sup> and $`8\times 10^{10}`$ g cm<sup>-3</sup>. Since the instability found in the present paper is localized within a density scale height from the crust which for the model atmosphere implies localization between densities $`10^{10}`$ and $`3\times 10^{10}`$ g cm<sup>-3</sup>, we expect that the electron capture effects would not significantly modify our conclusions.

In our analysis the crust was modelled as a rigid boundary in which the magnetic field is anchored. For this approximation to be justified the lateral distortion of the crust in response to the overpressure must be much smaller than the magnetic field line distortion in the overlying fluid. This is true when $`\mathrm{\Delta }\beta /\beta `$ exceeds $`\alpha _s\mathrm{\Delta }a/a`$ where $`\mathrm{\Delta }a(h^2/a)\mathrm{\Delta }\beta `$ is the displacement of flux surfaces in the ocean (cf. §6), i.e., when $`\beta >(a/h)^2/\alpha _s`$; here $`\alpha _sp`$ is the crust shear modulus. Estimates (Pandharipande et al., 1976; Strohmayer et al., 1991) find $`\alpha _s10^310^2`$, which for the characteristic density scale height and polar cap radius (cf. §1) implies that $`\beta 10^410^5`$. For the model atmosphere, this condition is satisfied for $`B5\times 10^{12}`$ G.

The imperfect rigidity of the crust may affect the estimate of the accreted mass necessary to destabilize the ballooning mode (eq. ). Ushomirsky et al. (2000) found that the crust tended to sink in response to a lateral density discontinuity. While their findings are not directly applicable to our case, it is nevertheless conceivable that the crust could deform sufficiently to significantly reduce the applied overpressure. The mass required to drive the instability would then be higher than our estimate.

The effect of the instability is an enhanced cross-field transport of matter and a resulting decrease in the plasma confinement time in the polar cap. To find the magnitude of this anomalous transport requires analysis of nonlinear effects which are beyond the scope of this paper. The analysis of these effects and, in particular, the question of the existence of a steady state, in which the transport is sufficiently rapid to balance the influx of the accreted matter, are relegated to future publications.

The authors thank Ellen Zweibel for pointing out the potential importance of the stabilizing effect of the Schwarzschild term in equation (4). Helpful discussions with Thierry Emonet, Timur Linde, and Dmitri Uzdensky are also acknowledged. This research was supported by the ASCI Center for Astrophysical Thermonuclear Flashes at the University of Chicago under Department of Energy contract B341495. EFB acknowledges support from an Enrico Fermi Fellowship at the University of Chicago.


text:
# 1 Basic Definitions

## 1 Basic Definitions

In this section we recall some basic definitions and statements from :
Let $`D`$ be a closed disk on $`^2`$, $`KD`$ finite set, $`uD`$. Any diffeomorphism of $`D`$ which fixes $`K`$ and is the identity on $`D`$ acts naturally on $`\mathrm{\Pi }_1=\mathrm{\Pi }_1(DK,u)`$. We say that two such diffeomorphisms of $`D`$ (which fix $`K`$ and equal identity on $`D`$) are equivalent if they define the same automorphism on $`\mathrm{\Pi }_1(DK,u)`$. This equivalence relation is compatible with composition of diffeomorphism and thus the equivalence classes form a group.

###### Definition 1.1

Braid Group $`B_n[D,K]`$.

Let $`D,K`$ be as above, $`n=\mathrm{\#}K`$, and let $``$ be the group of all diffeomorphisms $`\beta `$ of $`D`$ such that $`\beta (K)=K`$, $`\beta |_D=Id_D`$. For $`\beta _1,\beta _2`$ we say that $`\beta _1`$ is equivalent to $`\beta _2`$ if $`\beta _1`$ and $`\beta _2`$ define the same automorphism of $`\mathrm{\Pi }_1(DK,u)`$. The quotient of $``$ by this equivalence relation is called the Braid group $`B_n[D,K]`$.
Equivalently, if we take the canonical homomorphism $`\psi :Aut(\mathrm{\Pi }_1(DK,u))`$ then $`B_n[D,K]=Im(\psi )`$. The elements of $`B_n[D,K]`$ are called braids.

###### Lemma 1.2

If $`K^{}D^{}`$ in another pair as above with $`\mathrm{\#}K^{}=\mathrm{\#}K=n`$ then $`B_n[D^{},K^{}]`$ is isomorphic to $`B_n[D,K]`$.

This gives rise to the definition of $`B_n`$:

###### Definition 1.3

$`B_n=B_n[D,K]`$ for some $`D,K`$ with $`\mathrm{\#}K=n`$.

###### Definition 1.4

$`H(\sigma )`$, half-twist defined by $`\sigma `$.

Let $`D`$ and $`K`$ be defined as above. Let $`a,bK`$, $`K_{a,b}=K\backslash \{a,b\}`$ and $`\sigma `$ be a simple path (without a self intersection) in $`D\backslash D`$ connecting $`a`$ with $`b`$ such that $`\sigma K=\{a,b\}`$. Choose a small regular neighborhood $`U`$ of $`\sigma `$ such that $`K_{a,b}U=\varphi `$, and an orientation preserving diffeomorphism $`\psi :^2^1`$ such that $`\psi (\sigma )=[1,1]=\{z^1|Re(z)[1,1],Im(z)=0\}`$ and $`\psi (U)=\{z^1||z|<2\}`$. Let $`\alpha (r),r0`$, be a real smooth monotone function such that $`\alpha (r)=1`$ for $`r[0,3/2]`$ and $`\alpha (r)=0`$ for $`r2`$. Define a diffeomorphism $`h:^1^1`$ as follows: for $`z^1,z=re^{i\varphi }\text{let}h(z)=re^{i(\varphi +\alpha (r)\pi )}`$. It is clear that the restriction of $`h`$ to $`\{z^1||z|3/2\}`$ coincides with the $`180^{}`$ positive rotation, and that the restriction to $`\{z^1||z|2\}`$ is the identity map. The diffeomorphism $`\psi ^1h\psi `$ is inducing a braid called half-twist and denoted by $`H(\sigma )`$

###### Definition 1.5

Frame of $`B_n[D,K]`$

Let $`K=\{a_1,\mathrm{},a_n\}`$ and $`\sigma _1,\mathrm{},\sigma _{n1}`$ be a system of simple smooth paths in $`DD`$ such that $`\sigma _i`$ connects $`a_i\text{ with }a_{i+1}`$ and $`L=\sigma _i`$ is a simple smooth path. The ordered system of half-twists $`(H_1,\mathrm{},H_{n1})`$ defined by $`\{\sigma _i\}_{i=1}^{n1}`$ is called a frame of $`B_n[D,K]`$

###### Theorem 1.6

Let $`(H_1,\mathrm{},H_{n1})`$ be a frame of $`B_n[D,K]`$ then: $`B_n[D,K]`$ is generated by $`\{H_i\}_{i=1}^{n1}`$ and the following is a complete list of relations:
$`H_iH_j=H_jH_i\text{ if }|ij|>1`$ (Commutative relation)
$`H_iH_jH_i=H_jH_iH_j\text{ if }|ij|=1`$ (Triple relation)

Proof can be found for example in .
This theorem provides us with Artin’s algebraic definition of Braid group.

## 2 Properties of Hurwitz Moves

In this section we define Hurwitz Move and Hurwitz Equivalence and prove a certain theorem on the Hurwitz Equivalence of words in the Braid group which is correct for any Braid group of arbitrary order.

###### Definition 2.1

Hurwitz move on $`G^m`$ ($`R_k,R_k^1`$)

Let $`G`$ be a group, $`\stackrel{}{t}=(t_1,\mathrm{},t_m)G^m`$. We say that $`\stackrel{}{s}=(s_1,\mathrm{},s_m)G^m`$ is obtained from $`\stackrel{}{t}`$ by the Hurwitz move $`R_k`$ (or $`\stackrel{}{t}`$ is obtained from $`\stackrel{}{s}`$ by the Hurwitz move $`R_k^1`$) if

$$s_i=t_i\text{ for }ik,k+1,$$

$$s_k=t_kt_{k+1}t_k^1,s_{k+1}=t_k.$$

###### Definition 2.2

Hurwitz move on factorization

Let $`G`$ be a group and $`tG`$. Let $`t=t_1\mathrm{}t_m=s_1\mathrm{}s_m`$ be two factorized expressions of $`t`$. We say that $`s_1\mathrm{}s_m`$ is obtained from $`t_1\mathrm{}t_m`$ by the Hurwitz move $`R_k`$ if $`(s_1,\mathrm{},s_m)`$ is obtained from $`(t_1,\mathrm{},t_m)`$ by Hurwitz the move $`R_k`$.

###### Definition 2.3

Hurwitz equivalence of factorization

The factorizations $`s_1\mathrm{}s_m`$, $`t_1\mathrm{}t_m`$ are Hurwitz equivalent if they are obtained from each other by a finite sequence of Hurwitz moves. The notation is $`t_1\mathrm{}t_m\stackrel{HE}{}s_1\mathrm{}s_m`$.

###### Definition 2.4

Word in $`B_n`$

A word in $`B_n`$ is a representation of braid as a sequence of the frame elements and their inverses.

We will need a certain result of Garside:

###### Claim 2.5

(Garside): Every two positive words (all generators with positive powers) which are equal are transformable into each other through a finite sequence of positive words, such that each word of the sequence is positive and obtained from the proceeding one by a direct application of the commutative relation or the triple relation.

Proof: .

The following is obvious, never the less we give a short proof for clarification:

###### Claim 2.6

$`G`$ is a group $`g_1,g_2G`$:
1. If $`g_1g_2=g_2g_1`$ then $`g_1g_2\stackrel{HE}{}g_2g_1`$
2. If $`g_1g_2g_1=g_2g_1g_2`$ then $`g_1g_2g_1\stackrel{HE}{}g_2g_1g_2`$

Proof:
1. $`g_1g_2\stackrel{R_1}{}g_1g_2g_1^1g_1=g_2g_1`$.
2. $`g_1g_2g_1\stackrel{R_2}{}g_1g_2g_1g_2^1g_2\stackrel{R_1}{}g_1g_2g_1g_2^1g_1^1g_1g_2`$ but $`g_1g_2g_1=g_2g_1g_2`$ so $`g_1g_2g_1g_2^1g_1^1=g_2`$ and we get $`g_1g_2g_1\stackrel{HE}{}g_2g_1g_2`$

From Garside and the above we get the following proposition:

###### Proposition 2.7

Let $`H_1,\mathrm{},H_{n1}`$ be a set of generators of $`B_n`$ and $`H_{i_1}\mathrm{}H_{i_p}=H_{j_1}\mathrm{}H_{j_p}`$ two positive words (with $`i_k,j_k\{1,\mathrm{}n1\}`$) then $`H_{i_1}\mathrm{}H_{i_p}\stackrel{HE}{}H_{j_1}\mathrm{}H_{j_p}`$.

Proof:
Apply 2.5 on $`H_{i_1}\mathrm{}H_{i_p}=H_{j_1}\mathrm{}H_{j_p}`$ we get a finite sequence of positive words $`\{W_r\}_{r=0}^q`$ s.t. $`W_0=H_{i_1}\mathrm{}H_{i_p}`$, $`W_q=H_{j_1}\mathrm{}H_{j_p}`$ and $`W_{r+1}`$ is obtained from $`W_r`$ by a single application of the commutative relation or the triple relation.
An application of the triple relation (as in 2.6) is equal to an application of 2 Hurwitz moves $`R_{t+1},R_t`$. An application of the commutative relation is equal to an application of one Hurwitz move (on the commuting elements). Thus, $`W_0\stackrel{HE}{}W_q`$

###### Definition 2.8

$`\mathrm{\Delta }_n^2B_n[D,K]`$

$`\mathrm{\Delta }_n^2=(H_1\mathrm{}H_{n1})^n`$ where $`\{H_i\}_{i=1}^{n1}`$ is a frame.

We apply 2.7 on $`\mathrm{\Delta }_n^2`$:

###### Corollary 2.9

All $`\mathrm{\Delta }_n^2`$ factorizations $`H_{i_1}\mathrm{}H_{i_{n(n1)}}i_k\{1,\mathrm{},n1\}`$ are Hurwitz Equivalent.

## 3 $`B_3`$ Properties

This section concentrates on $`B_3`$ and we prove the following new result:

###### Theorem 3.1

Let $`\{H_1,H_2\},\{F_1,F_2\}`$ two frames of $`B_3`$, then for every two words $`\mathrm{\Delta }_3^2=H_{i_1}\mathrm{}H_{i_6}=F_{j_1}\mathrm{}F_{j_6}`$ (with $`i_k,j_k\{1,2\}`$) the factorizations $`H_{i_1}\mathrm{}H_{i_6}`$, $`F_{j_1}\mathrm{}F_{j_6}`$ are Hurwitz equivalent.

The proof of the main theorem is at the end of the section following from propositions 3.7 and 3.9.

Let $`D`$ be a closed disk on $`^2`$ and $`KD`$, $`K=\{a,b,c\}`$. We fix two parallel lines $`L_1,L_2`$ in $`D`$ which separate the points in $`K`$ (see figure 1).

With the choice of $`L_1`$ and $`L_2`$ we define the length and the minimal paths related to a half-twist:

###### Definition 3.2

Let $`\sigma `$ be a simple path in $`D\backslash D`$ that defines the half-twist $`H(\sigma )`$. Then $`m(\sigma )=\mathrm{\#}\{\sigma (L_1L_2)\}`$ as demonstrated in figure 1.

###### Definition 3.3

Let $`H`$ be a half-twist in $`B_n[D,K]`$. We define the length of $`H`$, $`m(H)=min\{m(\sigma )|H=H(\sigma )\}`$ and we call such a path $`\sigma `$ which gives the minimal value, the minimal path of $`H`$.

###### Notation 3.4

Any frame $`\{H_1,H_2\}`$ can be defined by the minimal paths $`\sigma _1,\sigma _2`$ which also satisfy the relations in 1.5.

###### Definition 3.5

We say that a frame $`\{H_1,H_2\}`$ is economic if there exist $`\sigma _1,\sigma _2`$ minimal paths of $`H_1,H_2`$ respectively s.t. for every small neighborhood $`U`$ of $`(\sigma _1K)\backslash \sigma _2`$, $`\sigma _1\sigma _2`$ outside $`U`$ (see 2)

###### Lemma 3.6

Let $`\{H_1,H_2\}`$ be economic frame then $`m(H_1)m(H_2)`$

Proof: Trivial.

###### Proposition 3.7

Any frame $`\{H_1,H_2\}`$ which generates $`B_3[D,K]`$ s.t $`max(m(H_1),m(H_2))>1`$ is economic.

Proof:
Suppose that there exists a frame $`\{H_1,H_2\}`$ as above which is not economic. From notation 3.4, there exist $`\sigma _1,\sigma _2`$ as defined in 1.5 s.t. $`\sigma _1,\sigma _2`$ are minimal paths of $`H_1,H_2`$ respectively and the two paths are ’separated’ at some point.
We will show that such a ’separation’ is not possible. Starting from $`s=\sigma _1\sigma _2`$ we will examine the first ’separation’ where the two paths are separated by one of the points in $`K`$.
The cases where $`s=a`$ and $`s=c`$ are symmetric. Proving the case where $`s=b`$ is similar.
We will assume that both paths start at $`a`$

Case 1: The two paths entering a neighborhood of the same point but with different orientation (one clockwise and the other counter-clockwise).

(1.a) Entering and turning back as shown in figure 3. The two paths create a cycle with one point inside, the internal path has no where to turn, and thus cannot end. This case is also relevant when one path is turning back and the other continues to the next point.

(1.b) Entering and continue to the next point:
(1.b.i) The two paths continue to the next point with the same orientation (figure 4). The two paths create a cycle that starts at $`a`$ and contains $`b`$. One path has to end at $`b`$ and will have to change it’s orientation through $`c`$. This will close a cycle with all the points inside, and the external path has no where to end.
(1.b.ii) The two paths continue to the next point with different orientation (figure 5). The two paths create a cycle with all points inside so the external path has no where to end.

Case 2: Entering to the same point with the same orientation and leaving to different points.

This can only happen in point $`b`$. The internal path must go back to $`a`$ and the external will go to $`c`$.
(2.a) The external path goes to $`c`$ by changing orientation (figure 6). This enforces the external path to go around $`a`$ and create a cycle with the internal path. Thus, the two paths create a cycle with all points inside and the external path has no where to end.

(2.b) The external path goes to $`c`$ with the same orientation (figure 7). The external path cannot go to the common part (this will create an empty cycle and the path has no where to end) of the paths so both paths continue to $`a`$. In this case we get an equivalent cycle as in 1.b.i.
This conclude the proof of proposition 3.7.

###### Lemma 3.8

Let $`\{H_1,H_2\}`$ be a frame generates $`B_3[D,K]`$ and $`H_\zeta `$ half-twist, than $`\{H{}_{\zeta }{}^{}{}_{}{}^{1}H_{1}^{}H_\zeta ,H{}_{\zeta }{}^{}{}_{}{}^{1}H_{2}^{}H_\zeta \}`$ is also a frame.

Proof:
By III.1.0 of ,

$$H_1^1H_2H_1=H(H_1(\sigma _2))\text{ where }H_2=H(\sigma _2).$$

$`\{H_1=H(\sigma _1),H_2=H(\sigma _2)\}`$ is a frame, therefore $`\sigma _1\sigma _2`$ is a simple smooth path. $`H_\zeta `$ is a diffeomorphism so $`H_\zeta (\sigma _1)H_\zeta (\sigma _2)=H_\zeta (\sigma _1\sigma _2)`$ is also a simple smooth path.
$`K\sigma _1\sigma _2,H_\zeta (K)=K`$ and therefore $`H_\zeta (\sigma _1),H_\zeta (\sigma _2)`$ connecting all points in $`K`$ appropriately.

###### Proposition 3.9

Let $`\{H_1,H_2\}`$ be a frame of $`B_3[D,K]`$ then there exist a frame $`\{F_1,F_2\}`$ s.t.

$$H_1H_2H_1H_2H_1H_2\stackrel{HE}{}F_1F_2F_1F_2F_1F_2$$

$$\text{and }max(m(H_1),m(H_2))>max(m(F_1),m(F_2)).$$

Proof:
From proposition 3.7 we get that $`\{H_1,H_2\}`$ is economic and from lemma 3.6 $`m(H_1)m(H_2)`$ we will assume that $`m(H_1)<m(H_2)`$.

Since $`\{H_1,H_2\}`$ is economic and $`m(H_1)<m(H_2)`$, operating $`H_1`$ on $`\sigma _2`$ (conjugation) will cause $`\sigma _2`$ to start coming back from the starting point ($`\sigma _1\sigma _2K`$) along the common part ($`\sigma _1`$) as shown in figure 8. Therefore, we get that either

$$m(H_1^1H_2H_1)=m(H(H_1(\sigma _2)))=m(H_2)m(H_1)\text{ or }$$

$$m(H_1H_2H_1^1)=m(H(H_1^1(\sigma _2)))=m(H_2)m(H_1)$$

In the case where $`m(H(H_1(\sigma _2)))=m(H_2)m(H_1)`$:
By performing the Hurwitz moves $`R_1R_0R_4R_3`$ we get:

$$H_1H_2H_1H_2H_1H_2\stackrel{HE}{}H_2H_1H_2H_1H_2H_1$$

and by performing $`R_0^1R_2^1R_4^1`$:

$$H_2H_1H_2H_1H_2H_1\stackrel{HE}{}H_1H_1^1H_2H_1H_1H_1^1H_2H_1H_1H_1^1H_2H_1$$

Since $`m(H_1^1H_2H_1)=m(H(H_1(\sigma _2)))`$ and $`m(H_2)>m(H_1)`$,

$$max(m(H_1),m(H_2))>max(m(H_1),m(H_1^1H_2H_1))$$

From Lemma 3.8, the set $`\{H_1,H_1^1H_2H_1\}`$ is a frame, we get an equivalent factorization of the same form but with a reduced length.
In the case where $`m(H(H_1^1(\sigma _2)))=m(H_2)m(H_1)`$:

$$H_1H_2H_1H_2H_1H_2\stackrel{HE}{}H_1H_2H_1^1H_1H_1H_2H_1^1H_1H_1H_2H_1^1H_1$$

By performing the Hurwitz moves $`R_0R_2R_4`$. Once again, we get a frame where the maximum length is reduced.

Proof of the main theorem:
Let $`H_{i_1}H_{i_2}H_{i_3}H_{i_4}H_{i_5}H_{i_6}`$ and $`F_{j_1}F_{j_2}F_{j_3}F_{j_4}F_{j_5}F_{j_6}`$ be two factorizations as above.
From 2.7 we get that:

$$H_{i_1}H_{i_2}H_{i_3}H_{i_4}H_{i_5}H_{i_6}\stackrel{HE}{}H_1H_2H_1H_2H_1H_2\text{and}$$

$$F_{j_1}F_{j_2}F_{j_3}F_{j_4}F_{j_5}F_{j_6}\stackrel{HE}{}F_1F_2F_1F_2F_1F_2$$

By 3.9,both $`H_1H_2H_1H_2H_1H_2`$ and $`F_1F_2F_1F_2F_1F_2`$ are Hurwitz equivalent to $`X_1X_2X_1X_2X_1X_2`$ where $`\{X_1,X_2\}`$ is a frame with $`m(X_1)=m(X_2)=1`$. Since such a frame is unique we get that

$$H_{i_1}H_{i_2}H_{i_3}H_{i_4}H_{i_5}H_{i_6}\stackrel{HE}{}F_{j_1}F_{j_2}F_{j_3}F_{j_4}F_{j_5}F_{j_6}$$


text:
# PERIODIC-ORBIT BIFURCATIONS AND SUPERDEFORMED SHELL STRUCTURE

## Abstract

We have derived a semiclassical trace formula for the level density of the three-dimensional spheroidal cavity. To overcome the divergences occurring at bifurcations and in the spherical limit, the trace integrals over the action-angle variables were performed using an improved stationary phase method. The resulting semiclassical level density oscillations and shell-correction energies are in good agreement with quantum-mechanical results. We find that the bifurcations of some dominant short periodic orbits lead to an enhancement of the shell structure for “superdeformed” shapes related to those known from atomic nuclei.

Keywords: Single-particle level density, periodic orbit theory, Gutzwiller’s trace formula, bifurcations, superdeformations.

PACS numbers: 03.65.Ge, 03.65.Sq, 05.45.Mt

Introduction — The periodic orbit theory (POT)  is a nice tool for studying the correspondence between classical and quantum mechanics and, in particular, the interplay of deterministic chaos and quantum-mechanical behavior. But also for systems with integrable or mixed classical dynamics, the POT leads to a deeper understanding of the origin of shell structure in finite fermion systems from such different areas as nuclear , metallic cluster , or mesoscopic semiconductor physics . Bifurcations of periodic orbits may have significant effects, e.g., in connection with the so-called “superdeformations” of atomic nuclei , and were recently shown to affect the quantum oscillations observed in the magneto-conductance of a mesoscopic device .

In the semiclassical trace formulae that connect the quantum-mechanical density of states with a sum over the periodic orbits of the classical system , divergences arise at critical points where bifurcations of periodic orbits occur or where symmetry breaking (or restoring) transitions take place. At these points the stationary-phase approximation, used in the semiclassical evaluation of the trace integrals, breaks down. Various ways of avoiding these divergences have been studied , some of them employing uniform approximations . Here we employ an improved stationary-phase method (ISPM) for the evaluation of the trace integrals in the phase-space representation, based on the studies in Refs. , which we have derived for the elliptic billiard . It yields a semiclassical level density that is regular at all bifurcation points of the short diameter orbit (and its repetitions) and in the circular (disk) limit. Away from the critical points, our result reduces to the extended Gutzwiller trace formula  and is identical to that of Berry and Tabor  for the leading-order families of periodic orbits.

The main purpose of the present note is to report on the extension of our semiclassical approach to the three-dimensional (3D) spheroidal cavity , which may be taken as a simple model for a large deformed nucleus  or a (highly idealized) deformed metal cluster , and to specify the role of orbit bifurcations in the shell structure responsible for the superdeformation. Although the spheroidal cavity is integrable (see, e.g., Ref. ), it exhibits all the difficulties mentioned above (i.e., bifurcations and symmetry breaking) and therefore gives rise to an exemplary case study of a non-trivial 3D system. We apply the ISPM for the bifurcating orbits and succeed in reproducing the superdeformed shell structure by the POT, hereby observing a considerable enhancement of the shell-structure amplitude near the bifurcation points.

Theory — The level density $`g(E)`$ is obtained from the semiclassical Green function  by taking the imaginary part of its trace in $`(𝐈,𝚯)`$ action-angle variables :

$`g(E)={\displaystyle \underset{i}{}}\delta \left(E\epsilon _i\right)\mathrm{Re}{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{\mathrm{d}𝐈^{}\mathrm{d}𝚯^{\prime \prime }}{\left(2\pi \mathrm{}\right)^3}\delta \left(EH\right)}`$ (1)
$`\times \mathrm{exp}\left\{{\displaystyle \frac{i}{\mathrm{}}}\left[S_\alpha (𝐈^{},𝐈^{\prime \prime },t_\alpha )+\left(𝐈^{\prime \prime }𝐈^{}\right)𝚯^{\prime \prime }\right]i{\displaystyle \frac{\pi }{2}}\mu _\alpha \right\}.`$ (2)
Here $`\left\{\epsilon _i\right\}`$ is the single-particle energy spectrum and $`H=H(𝐈)`$ is the classical Hamiltonian. The sum is taken over all classical trajectories $`\alpha `$ specified by the initial actions $`𝐈^{}`$ and final angles $`𝚯^{\prime \prime }`$. $`S_\alpha (𝐈^{},𝐈^{\prime \prime },t_\alpha )=_𝐈^{}^{𝐈^{\prime \prime }}𝐈d𝚯`$ is the action integral and $`t_\alpha `$ the time for the motion along the trajectory $`\alpha `$, and $`\mu _\alpha `$ is the Maslov index related to the caustic and the turning points . In the spheroidal variables $`\{u,v,\phi \}`$, the action $`𝐈`$ has the components

$`I_u`$ $`=`$ $`{\displaystyle \frac{pc}{\pi }}{\displaystyle \underset{u_c}{\overset{u_c}{}}}du\sqrt{\sigma _1\mathrm{sin}^2u\sigma _2/\mathrm{cos}^2u},`$ (3)
$`I_v`$ $`=`$ $`{\displaystyle \frac{pc}{\pi }}{\displaystyle \underset{v_c}{\overset{v_t}{}}}dv\sqrt{\mathrm{cosh}^2v\sigma _1\sigma _2/\mathrm{sinh}^2u},`$ (4)
$`I_\phi `$ $`=`$ $`pc\sqrt{\sigma _2}.`$ (5)
Hereby $`p=(2mE)^{1/2}`$ is the particle’s momentum and $`c=(b^2a^2)^{1/2}`$ half the distance between the foci; $`b`$ and $`a`$ are the semiaxes (with $`b>a`$) of the spheroid with its volume fixed by $`a^2b=R^3`$ and the axis ratio $`\eta =b/a`$ as deformation parameter; and $`\pm u_c`$ (or $`v_c`$) and $`v_t`$ are the caustic and turning points, respectively. In Eq. (5) we use the dimensionless “action” variables $`\sigma _1`$, $`\sigma _2`$  in which the torus of the classical motion is given by

$`\sigma _2^{}=0`$ $`\sigma _2`$ $`{\displaystyle \frac{1}{\eta ^21}}=\sigma _2^+,`$ (6)
$`\sigma _1^{}=\sigma _2`$ $`\sigma _1`$ $`{\displaystyle \frac{\eta ^2}{\eta ^21}}\sigma _2\left(\eta ^21\right)=\sigma _1^+.`$ (7)
In the ISPM, we expand the action $`S_\alpha `$ in (2) up to second-order terms around its stationary points and keep the pre-exponential factor at zero order, taking the integrations over the torus within the finite limits given by Eq. (7). For the oscillating (“shell-correction”) part of the level density $`\delta g(E)=g(E)\stackrel{~}{g}(E)`$, where $`\stackrel{~}{g}(E)`$ is its smooth part , we obtain

$$\delta g(E)\frac{1}{E_0}\mathrm{Re}\underset{\beta }{}A_\beta (E)\mathrm{exp}\left(ikL_\beta i\frac{\pi }{2}\mu _\beta \right)w_\beta ^\gamma ,$$
(8)
where $`k=p/\mathrm{}`$ is the wave number and $`E_0=\mathrm{}^2/2mR^2`$ our energy unit. The amplitudes $`A_\beta `$ will be specified below. The sum over $`\beta `$ includes all 2-parameter families of three-dimensional (3D) periodic orbits and elliptic and hyperbolic 2D orbits lying in a plane containing the symmetry axis (all with degeneracy parameter $`𝒦`$=2), the 1-parameter families of (2D) equatorial orbits lying in the central plane perpendicular to the symmetry axis (with $`𝒦`$=1), and the (1D) isolated long diameter (with $`𝒦`$=0). In Eq. (8), $`L_\beta `$ is the length of the orbit $`\beta `$ at the stationary point ($`\sigma _1^{},\sigma _2^{}`$) which for the 3D orbits lies inside the physical region of the torus (7), and is analytically continued outside this region. The $`\sigma _2=0`$ boundary of (7) is occupied by the 2D orbits with $`𝒦`$=2. The stationary points are determined by the roots of the periodicity conditions $`\omega _u/\omega _v=n_u/n_v`$ and $`\omega _u/\omega _\phi =n_u/n_\phi `$; hereby $`\omega _\kappa =H/I_\kappa `$ are the frequencies and $`n_\kappa `$ are co-prime integers which specify the periodic orbits $`\beta =M(n_v,n_\phi ,n_u)`$, where $`M`$ is the repetition number. The factor $`w_\beta ^\gamma =\mathrm{exp}(\gamma ^2L_\beta ^2/4R^2)`$ in Eq. (8) is the result of a convolution of the level density with a Gaussian function over a range $`\gamma `$ in the variable $`kR`$. This ensures the convergence of the POT sum (8) by suppressing the longer orbits which are not relevant for the coarse-grained gross-shell structure .

For Strutinsky’s shell-correction energy $`\delta U`$ , we obtain (with time reversal symmetry and a spin factor 2)

$`\delta U`$ $`=`$ $`\mathrm{\hspace{0.33em}2}{\displaystyle \underset{i=1}{\overset{N/2}{}}}\epsilon _i2{\displaystyle _0^{\stackrel{~}{E}_F}}E\stackrel{~}{g}(E)dE`$ (9)
$``$ $`\mathrm{\hspace{0.33em}8}R^2E_F\mathrm{Re}{\displaystyle \underset{\beta }{}}{\displaystyle \frac{A_\beta (E_F)}{L_\beta ^2}}\mathrm{exp}(ik_FL_\beta i{\displaystyle \frac{\pi }{2}}\mu _\beta ).`$ (10)
The Fermi energies $`E_F`$ (and with it $`k_F`$) and $`\stackrel{~}{E}_F`$ are determined by the particle number conservation $`N=2_0^{E_F}g(E)dE=2_0^{\stackrel{~}{E}_F}\stackrel{~}{g}(E)dE`$. Due to the factor $`L_\beta ^2`$, the PO sum in (10) may converge faster for the shortest orbits than the level density (8) for small $`\gamma `$. Any enhancement of the amplitudes $`A_\beta `$ of the most degenerate short periodic orbits – e.g., due to bifurcations or to symmetry restoring, as discussed below – therefore leads to an enhancement of the shell structure and hence to an increased stability of the system.

We present here only the amplitudes of the leading contributions to (8) and (10). For further details (including, e.g., explicit expressions for the Maslov indices $`\mu _\alpha `$), we refer to a forthcoming, more extensive publication .

For the amplitudes $`A_\beta `$ of the most degenerate ($`𝒦`$=2) families of periodic 3D and 2D orbits, we obtain

$$A_\beta ^{(𝒦=2)}=\frac{icL_\beta [I_u/\sigma _1]_{\sigma _n^{}}}{\pi (4MRn_v)^2\sqrt{K_\beta \sigma _2^{}}}\underset{n=1}{\overset{2}{}}\text{erf}(x_n^{},x_n^+).$$
(11)
The quantity $`K_\beta =K_\beta ^{(1)}K_\beta ^{(2)}`$ is related to the main curvatures $`K_\beta ^{(n)}`$ of the energy surface $`E=H(\sigma _1,\sigma _2)`$ in the “action” plane $`(\sigma _1,\sigma _2)`$, given by

$$K_\beta ^{(n)}=[\frac{^2I_v}{\sigma _n^2}+\frac{\omega _u}{\omega _v}\frac{^2I_u}{\sigma _n^2}+\frac{\omega _\phi }{\omega _v}\frac{^2I_\phi }{\sigma _n^2}]_{\sigma _n^{}}.(n=1,2)$$
(12)
In Eq. (11), the arguments of the two-dimensional error function erf$`(x,y)=2_x^y𝑑ze^{z^2}/\sqrt{\pi }`$ are given by the turning points in the action plane

$$x_n^\pm =\sqrt{i\pi Mn_vK_\beta ^{(n)}/\mathrm{}}\left(\sigma _n^\pm \sigma _n^{}\right);(n=1,2)$$
(13)
see Eq. (7) for the boundaries $`\sigma _n^\pm `$. All quantities in (11) can be expressed analytically in terms of elliptic integrals. For the 3D orbits, our result (11) is in agreement with that obtained by exact Poisson summation over the EBK spectrum (cf. Refs. ).

For the contribution of the $`𝒦=1`$ families of equatorial orbits to (8), we obtain the amplitudes

$$A_\beta ^{(𝒦=1)}=\sqrt{\frac{i\mathrm{sin}^3\varphi _\beta }{\pi Mn_vkR\eta F_\beta }}\underset{n=1}{\overset{3}{}}\text{erf}(x_n^{},x_n^+),$$
(14)
where $`\varphi _\beta =\pi n_\phi /n_v`$, $`F_\beta `$ is their stability factor , $`\sigma _1^{}=\sigma _2^{}=\mathrm{cos}^2\varphi _\beta /(\eta ^21)`$, and

$$x_3^+=kc\sqrt{\frac{i\pi \mathrm{}F_\beta }{64Mn_v(\sigma _2^{}+1)K_\beta ^{(1)}}},x_3^{}=0.$$
(15)
The contribution of the isolated long diameter orbit, which may be expressed in terms of incomplete Airy integrals , is not important for deformations of the order $`\eta 1.22`$.

Discussion of results — In Fig. 1 we show $`|A_\beta |`$ versus deformation $`\eta `$ (at $`kR=40`$) for a pair of orbits involved in a typical bifurcation scenario. At the critical point $`\eta =1.618\mathrm{}`$ the equatorial “star” orbit (5,2) undergoes a bifurcation at which the 3D orbit (5,2,1) is born; the latter does not exist below $`\eta =1.618\mathrm{}`$

In the standard stationary-phase approach (SSPM; dash-dotted lines), the amplitude of the (5,2) orbit diverges at $`\eta =1.618\mathrm{}`$, whereas that of the bifurcated orbit (5,2,1) is finite but discontinous. As seen in Fig. 1, the ISPM (solid lines) leads to a finite amplitude $`A_{(5,2)}^{(𝒦=1)}`$ for the (5,2) orbit. This is because the factor $`F_\beta `$ in the denominator of (14), which goes to zero at the bifurcation, is cancelled by the same factor in the numerator of $`x_3^+`$ (15) via the third error function in (14). A similar result was found for the short diameter orbit 2(2,1) in the elliptic billiard . Furthermore, the ISPM softens the discontinuity for the (5,2,1) orbit, leading to a maximum amplitude slightly above the critical deformation.

The relative enhancement of these amplitudes $`A_\beta `$ near the bifurcation point can also be understood qualitatively from the following argument. At the bifurcation of the equatorial (5,2) orbit, its degeneracy parameter $`𝒦=1`$ locally increases to 2, because it is there degenerate with the orbit family (5,2,1) that has $`𝒦=2`$ at all deformations $`\eta 1.618\mathrm{}`$ This is similar to a symmetry restoring transition. An increase of the symmetry parameter $`𝒦`$ by one unit leads to one more exact integration compared to the SSPM, and thus the amplitudes (11) and (14) acquire an enhancement factor $`\sqrt{kL_\beta }\sqrt{pR/\mathrm{}}`$ (cf. Refs. ).

A similar enhancement of the double triangle 2(3,1) and the 3D orbit (6,2,1) is found near their bifurcation point $`\eta =\sqrt{3}=1.732\mathrm{}`$ However, the curvature $`K_\beta ^{(1)}`$ (12) for orbits like $`M(3t,t,1)`$ ($`t=2,3`$,…) is identically zero and hence the SSPM is divergent for all deformations $`\eta 1`$, in contrast to the situation with orbits like (5,2,1) with finite $`K_\beta ^{(1)}`$. Here we have to take into account the next nonzero 3rd-order terms in the expansion of $`S_\alpha `$, although the $`(3t,t,1)`$ ISPM amplitude (11) is finite and continuous everywhere. The amplitude can then be expressed in terms of incomplete Airy and Gairy integrals with finite limits . For the equatorial orbits $`t`$(3,1), like for the double triangles 2(3,1), one has a zero curvature $`K_\beta ^{(1)}`$ only at the bifurcation point $`\eta =\sqrt{3}`$. Here $`F_\beta /K_\beta ^{(1)}0`$, and a similar mechanism of cancellation of singularities for other orbits takes place through Eqs. (13-15). But the relative enhancement of the ISPM amplitudes (11,14) of such orbits at the bifurcations is of order $`kL_\beta `$ because of a change of the degeneracy parameter $`𝒦`$ by two units (see  for details). In this sense we avoid here a double singularity related to a double restoring of symmetry.

In Figs. 2 and 3, we present semiclassical level densities $`\delta g(E)`$ (8) versus $`kR`$ and shell-correction energies $`\delta U`$ (10) versus $`N^{1/3}`$ for various critical deformations (heavy dotted lines), and compare them to the corresponding quantum-mechanical results (thin solid lines). We observe a very good agreement of the gross-shell structure at all deformations. The most significant contributions to these results near the critical deformations are coming from bifurcating orbits with lengths smaller than about 10$`R`$, in line with the convergence arguments for the POT sums (8) and (10) mentioned above. For the bifurcation at $`\eta =1.618\mathrm{}`$, the orbits (5,2,1) and (5,2) give contributions comparable with other 2D orbits. For $`\eta =\sqrt{3}`$, the bifurcating orbits (6,2,1) and (6,2) are also important.

The role of the bifurcating orbits increases for larger deformations and is dominating at the superdeformation $`\eta =2`$. For this deformation, the most important orbits in the present spheroidal cavity model are the 3D orbits (5,2,1), (6,2,1), (7,2,1), and (8,2,1).

These results are in agreement with both heights and positions of the peaks in the length spectra obtained in Ref.  from the Fourier transforms of the quantum level densities $`g(kR)`$ at the same deformations.

Summary and conclusions — We have obtained an analytical trace formula for the 3D spheroidal cavity model, which is continuous through all critical deformations where bifurcations of periodic orbits occur. We find an enhancement of the amplitudes $`|A_\beta |`$ at deformations $`\eta 1.62.0`$ due to bifurcations of 3D orbits from the shortest 2D orbits. We believe that this is an important mechanism which contributes to the stability of superdeformed systems. Our semiclassical analysis may therefore lead to a deeper understanding of shell structure effects in superdeformed fermionic systems – not only in nuclei or metal clusters but also, e.g., in deformed semiconductor quantum dots whose conductance and magnetic susceptibilities are significantly modified by shell effects.

Acknowledgements — A.G.M. gratefully acknowledges the financial support provided under the COE Professorship Program by the Ministry of Education, Science, Sports and Culture of Japan (Monbu-sho), giving him the opportunity to work at the RCNP, and thanks Prof. H. Toki for his warm hospitality and fruitful discussions. Two of us (A.G.M. and S.N.F.) acknowledge financial support by the Regensburger Universitätsstiftung Hans Vielberth.