text:
# Phase Transitions in Isolated Vortex Chains

## I Introduction

The vortex system in layered superconductors,Blatterreview  which includes the high-$`T_c`$ cuprates, displays a rich set of physical phenomena, such as a thermodynamic melting transition,Zeldov  a pinning-induced disordering transition,vanderBeek  and various structural transitions at different angles.Mirkovic ; Ooi2000  The average density and orientation of the vortices are controlled by the magnetic field, because each vortex carries one quantum of flux, $`\mathrm{\Phi }_0=hc/2e`$. In this paper we are concerned with the vortex chains that appear in a certain regime of tilted magnetic fields. These chains, which consist of a high density of flux lines perpendicular to the layers, were first observed with the Bitter decoration technique by Bolle et al.Bolle  A qualitative explanation followed shortlyHuse  in terms of the proposed crossing-lattice state in tilted fields.Theodorakis ; Bulaevskii92  This state consists of a lattice of flux lines perpendicular to the layers crossed by a lattice of flux lines along the layers. The in-plane lattice is strongly distorted due to the anisotropy,Campbell  and has a large spacing along the layers. Huse surmised that a possible attractive interaction between the two species of flux line would lead to a higher density of out-of-plane flux lines along chains, with an inter-chain separation equal to the in-plane flux-line spacing. This picture seems to be consistent with the experimental observations.Grigorieva  The energy of the crossing-lattice state was considered by Benkraouda and Ledvij,Benkraouda  who found a transition from a single lattice of tilted flux lines to the crossing-lattice state as the tilt angle is increased for sufficiently large anisotropy. Their work, however, neglected interactions between the two crossing lattices.

Interest in the crossing-lattice state has been revitalized in the last couple of years. KoshelevKoshelev99  has shown that the regime for crossing lattices is larger than previously expected. This is because a correct treatment must include the interactions between the perpendicular flux lines, and a distortion away from ideal crossing lattices has a lower energy. Experimental evidence comes in the form of the unusual dependence of the melting temperature of the vortex lattice as a function of magnetic field angle,Schmidt ; Ooi1999  which is consistent with the crossing-lattice state rather than a single tilted flux-line lattice.Koshelev99  In addition, the work of Koshelev explained quantitatively the attraction of perpendicular flux lines that leads to the vortex-chain state in a certain field range. Apart from more detailed work on the melting transition,Mirkovic  and on magnetization curves,Ooi2000  these concepts have also inspired recent scanning Hall-probe measurementsGrigorenko  on the vortex-chain state. These experiments are similar to the Bitter decoration technique in that they only probe the field distribution emanating from the surface of the superconductor. However, they have the advantage of speed and resolution, and are not a “one-off” measurement, so that the system can be finely tuned to observe different effects. These new experiments have shown unexpected properties within the chain state. In particular, it is possible to tune to a field range where all of the out-of-plane flux lines are arranged in chains. It is this “isolated-chain state” that we will study in this paper.

We will use the following geometry. The $`z`$-axis is perpendicular to the layers, and the component of magnetic induction $`B_z`$ determines the density of “pancake” vorticesPancakes  in every superconducting layer. The in-plane field is along the $`x`$-axis, and $`B_x`$ gives the density of the flux lines in this direction. These in-plane flux lines have their centers in the spacing between layers, and so are commonly called Josephson vortices.Jvortex  A $`z`$-directed flux line is made from a stack of pancake vortices, and contains circulating currents in the layers up to a distance of the penetration depth, $`\lambda _{\mathrm{ab}}\lambda `$, from the vortex center. A Josephson vortex in the $`x`$-direction has an ellipsoidal current pattern, with the flux and currents decaying over the much larger distance $`\lambda _\mathrm{c}=\gamma \lambda `$ in the $`y`$-direction. Here, $`\gamma `$ is the anisotropy ratio, which is large for weakly coupled layers. The phase singularity of the Josephson vortex is confined between two neighboring layers, with separation $`d`$, where the phase difference across the layers changes by $`2\pi `$ over a distance of the Josephson length $`\gamma d`$.

The most simple structure of the crossing-lattice state is when an ideal triangular lattice (with spacing $`a_\mathrm{P}=[(2/\sqrt{3})\mathrm{\Phi }_0/B_z]^{1/2}`$) of pancake-vortex stacks crosses a stretched triangular lattice of Josephson vortices, with separation in the $`z`$-direction of $`a_{\mathrm{J}z}=[(2/\sqrt{3})\mathrm{\Phi }_0/B_x]^{1/2}/\sqrt{\gamma }`$. Koshelev has shown how the interactions between pancake stacks and Josephson vortices lead to an effective attraction.Koshelev99  At the low pancake densities we are interested in here ($`B_z<\mathrm{\Phi }_0/\lambda ^2`$) this means that there will be large distortions from the ideal triangular pancake-vortex lattice, with a higher density of pancake stacks located over the centers of the Josephson vortices. This leads to the high-density chains of pancake stacks that are observed by Bitter decoration,Bolle ; Grigorieva  with an inter-chain spacing of $`a_{\mathrm{J}y}=\sqrt{\gamma }[(\sqrt{3}/2)\mathrm{\Phi }_0/B_x]^{1/2}`$. Eventually, for small $`B_z\mathrm{\Phi }_0/\lambda ^2`$, all of the pancakes lie over the Josephson vortices. This isolated-chain state has a pancake-stack separation along the chains of $`a_{\mathrm{ch}}=(\mathrm{\Phi }_0/B_z)/a_{\mathrm{J}y}`$, see Fig. 1.

The crossing lattices must compete energetically with the more conventional tilted lattice of vortices. This is quite similar to the vortex lattice in a continuous superconductor, only there is a kinked structure along each vortex with a periodicity $`L=d(B_x/B_z)`$. The anisotropy of the currents mean that there is a stretched aspect ratio of this lattice. In addition, this distortion may be different from that expected from simple rescalingBlatterGL  due to the attraction of tilted vortices along the direction of tilt. This may lead to a “tilted chain state”,TiltedChains  which is distinct from the crossing chains mentioned above. Koshelev showed that for intermediate pancake densities $`B_z>\mathrm{\Phi }_0/\lambda ^2`$ there should be a transition at a small value of $`B_x=B_x^{}0.01B_z`$ between the tilted lattice ($`B_x<B_x^{}`$) and the crossing lattices ($`B_x^{}<B_x`$).

In this work we will concentrate on very small out-of-plane fields and $`B_zB_x<\mathrm{\Phi }_0/\lambda ^2`$ when all of the pancake stacks become attached to the centers of the Josephson vortices, giving the isolated-chain state observed in recent scanning Hall probe measurements.Grigorenko  We find that the stability of the isolated-chain state is quite delicate. We first note in Section II that the existence of a stable crossing configuration in the isolated chains is only possible for large enough anisotropy, layer spacing, and pancake separation. We also find that the energy of the isolated chain as a function of pancake-stack separation has a minimum, i.e., there is an optimum density of pancake stacks along a chain, and there will be an instability towards clustering when the total density is low. As the density of pancake stacks along the chain increases (at a fixed density of Josephson vortices), the chain will buckle (Section III), and then eject vortices at a critical minimum separation (Section IV). More surprisingly, at very low pancake density the chain state may be replaced by a chain of tilted stacks (Section V). These results are summarized in Figs. 2 and 3, where we plot the calculated phase boundaries for three first-order transitions: clustering, buckling and tilting as well as a continuous ejection transition. Finally, in Section VI we discuss the extent to which these transitions have been observed experimentally, and consider the effect of fluctuations due to finite temperature or quenched pinning disorder.

We conclude this introduction with a note on parameters. The experiments we have referred to,Bolle ; Grigorieva ; Grigorenko  were performed on the extremely anisotropic cuprate Bi2Sr2CaCu2Ox (BiSCCO). Where we make explicit calculations we will take the following appropriate parameters $`\gamma =500`$, $`\lambda =2000`$ Å, and $`d=15`$ Å. This gives a Josephson length (the size of the non-linear core of a Josephson vortex) of $`\gamma d=7500`$ Å $`\lambda `$.

## II Structure of the Isolated Vortex Chain

Koshelev has estimated the structure of a pancake stack that crosses the center of a single Josephson vortex.Koshelev99  In the crossing configuration there is a Lorentz force on the pancakes in a direction parallel to the Josephson vortex, which must balance the attraction of each pancake to its stack. This distorts the pancake stack, with the largest displacements by the two pancakes immediately above and below the Josephson core. The distortions give an energy gain of $`\mathrm{\Delta }E_\times 8.4\epsilon _0d(\lambda /\gamma d)^2/\mathrm{ln}(3.5\gamma d/\lambda )`$. Here $`\epsilon _0d=(\mathrm{\Phi }_0/4\pi \lambda )^2d`$ is the typical energy scale for pancake interactions. This causes the pancake stack to be attracted to a stack of Josephson vortices within the crossing-lattice state. Koshelev’s result uses a quadratic approximation for the pancake-to-stack attraction. Recent workDodgson\_Crete  using the correct potential has shown that, while Koshelev’s estimate for the crossing energy is close to the full result, a stable crossing configuration only exists for extreme enough anisotropy, $`\gamma >\gamma _{\mathrm{min}}=2.86\lambda /d`$. In this section we include the interaction between pancakes in different stacks separated by $`a_{\mathrm{ch}}`$ along an isolated chain. For finite $`a_{\mathrm{ch}}`$ we find a further reduction of the stability limit of crossing found in Ref. Dodgson\_Crete, . We also find a minimum in the energy of the chain as a function of stack separation, which leads to a clustering instability when $`a_{\mathrm{ch}}\lambda `$.

To calculate the crossing configuration of the chain, we use the following assumptions: First, we neglect the effect of induced Josephson currents due to displacing pancakes from their stacks (reasonable for large values of $`\gamma d/\lambda `$ and small enough displacements $`u<\gamma d`$). Second, we utilize the long range of the remaining electromagnetic pancake interactions,Pancakes  of order $`\lambda `$ in the $`z`$-direction. There are many ($`\lambda /d10^2`$) pancakes that contribute to the current distribution in one layer, determining the potential felt by a given pancake. We can therefore reduce this many-body optimization problem to a one-dimensional problem, considering the displacement of a row of pancakes in one layer under the potential due to the ideal chain of pancake stacks (the corrections due to the displacements in other layers is small when the number of pancakes with large displacements is much less than $`\lambda /d`$).

Within this scheme, the energy profile (per stack) for displacing by $`u_n`$ a row of pancakes in the $`n`$th layer when there is a Josephson vortex between layers $`0`$ and $`1`$ crossing the chain is,

$$\mathrm{\Delta }E_n(u_n)=\frac{\mathrm{\Phi }_0d}{c}J_n^yu_n+V_{\mathrm{em}}^{\mathrm{row}}(u_n).$$
(1)
The first term here comes from the Lorentz force on the pancakes due to the in-plane current density from the Josephson vortex, $`J_n^y`$, and tends to pull the pancakes away from their stacks. The form of this current will be discussed more in Section III, but the numerical value is taken from Ref. Koshelev\_preprint, , e.g., the current in the $`n=1`$ layer (immediately above the Josephson vortex) is $`J_1^y=2.28\epsilon _0c/\mathrm{\Phi }_0\gamma d`$. The second term is the attractive magnetic interaction of the pancake row with the remainder of all the stacks in the chain. Clem showed for a single pancake,Clem  that this has the same interaction energy as the sum of a pancake with a full stack, $`V_{\mathrm{em}}^{\mathrm{stack}}(R)=2\epsilon _0dK_0(R/\lambda )`$ plus a pancake with its anti-image $`V_{\mathrm{em}}^{\mathrm{pc}\mathrm{pair}}=2\epsilon _0d\mathrm{ln}(R/L)`$, so that for the entire row we find

$`V_{\mathrm{em}}^{\mathrm{row}}(u)`$ $`=`$ $`2\epsilon _0d\left[K_0\left(u/\lambda \right)+\mathrm{ln}(u/2\lambda )+\gamma _\mathrm{E}\right]`$
$`+2\epsilon _0d{\displaystyle \underset{j0}{}}K_0\left[{\displaystyle \frac{ja_{\mathrm{ch}}+u}{\lambda }}\right]K_0\left({\displaystyle \frac{ja_{\mathrm{ch}}}{\lambda }}\right)+\mathrm{ln}\left[{\displaystyle \frac{ja_{\mathrm{ch}}+u}{ja_{\mathrm{ch}}}}\right].`$
\[Euler’s constant $`\gamma _\mathrm{E}`$ is needed to fix the zero of energy to that of fully aligned stacks, c.f. the small $`x`$ expansion of the modified Bessel function $`K_0(x)`$\].

Such energy profiles for the $`n=1`$ pancake row are shown in Fig. 4 for different values of $`a_{\mathrm{ch}}`$ with the choice of the ratio $`\gamma d/\lambda =3.75`$ (a reasonable choice for BiSCCO). For dilute chains there is a stable minimum (as expected from the results of Ref. Dodgson\_Crete, ), which determines the optimal displacement. However, this stability disappears once the separation is below the critical value of $`a_{\mathrm{ch}}^{\mathrm{min}}=5.3\lambda `$. There can be no stable isolated-chain state for densities higher than this critical value.relax\_caveat  Fig. 4 also suggests that the crossing configuration is only metastable, with a lower energy at $`u_1=0.5a_{\mathrm{ch}}`$, but this is not reliable as the simple Lorentz force argument does not hold for pancake displacements of such a large fraction of the stack-separation. For such displacements the original Josephson vortex becomes completely fragmented, and a different approach is needed. In fact, the competing state here is a “tilted chain”, and in Section V we will calculate the energy of the tilted chain, to compare to the isolated chain state. (A third possibility is a soliton-like structureKoshelev\_private  that we will not consider here.)

The displacements of the entire stacks are shown in Fig. 1 for $`a_{\mathrm{ch}}=10\lambda `$, and a Josephson-vortex separation of $`a_{\mathrm{J}z}=20d`$ (this corresponds to $`B_x=53`$ G for our parameters). All of the displacements contribute to the total crossing energy gain. We find that the crossing energy per pancake stack only lowers from $`\mathrm{\Delta }E_\times (\mathrm{})0.21\epsilon _0d`$ for an isolated stack to $`\mathrm{\Delta }E_\times (a_{\mathrm{ch}}^{\mathrm{min}})0.26\epsilon _0d`$ as the stacks reach the critical separation $`a_{\mathrm{ch}}^{\mathrm{min}}`$. Therefore within the region of stable crossing, there is a pinning energy to the Josephson vortex centers per unit length of pancake stack,

$$\epsilon _{\mathrm{Jv}\mathrm{pin}}=\mathrm{\Delta }E_\times /a_{\mathrm{J}z}0.2\epsilon _0(d/a_{\mathrm{J}z}),$$
(3)
where the last result is with our parameters for BiSCCO.

In Fig. 5a we plot the total energy per pancake of the isolated chain $`E_{\mathrm{ch}}=\frac{1}{2}_{j0}V_{\mathrm{em}}^{\mathrm{stack}}(ja_{\mathrm{ch}})+(d/a_{\mathrm{J}z})\mathrm{\Delta }E_\times (a_{\mathrm{ch}})`$. Note that there is a minimum at $`a_{\mathrm{ch}}=8.2\lambda `$, reflecting the fact that at large separations the pancake stacks attract each other. This unusual feature has been explained by Buzdin and BaladieBuzdin\_Baladie  by considering the distorted pancake stacks of Fig. 1 as the superposition of a straight stack plus a series of pancake–anti-pancake dipoles. The straight stacks give a repulsive term, but this is exponentially small for $`a_{\mathrm{ch}}\lambda `$, while the dipoles give a weak attractive term that only falls off like $`1/a_{\mathrm{ch}}^2`$ (note there is some similarity to the attraction between tilted flux linesTiltedChains ). Therefore there is a net attraction at large distances, which will destabilize a homogeneous chain at very low densities.

To formally describe this instability we should consider the free energy density as a function of out-of-plane field, $`f(B_z)=f_{\mathrm{ch}}(B_z)+\epsilon _{\mathrm{stack}}B_z/\mathrm{\Phi }_0`$. Here $`\epsilon _{\mathrm{stack}}`$ is the line energy of a pancake stack, $`f_{\mathrm{ch}}(B_z)=E_{\mathrm{ch}}(a_{\mathrm{ch}})/a_{\mathrm{ch}}a_{\mathrm{J}y}d`$ and $`B_z=\mathrm{\Phi }_0/a_{\mathrm{ch}}a_{\mathrm{J}y}`$. In Fig. 5b we plot $`f_{\mathrm{ch}}(B_z)`$ using the same data in Fig. 5a. We also note that the thermodynamically stable phase is determined by the minimum Gibbs free energy $`g(H_z)=f[B_z(H_z)]B_zH_z/4\pi `$ with $`H_z=4\pi f/B_z`$. By a geometric construction we see that the point $`B_z=B_z^1`$ on $`f_{\mathrm{ch}}(B_z)`$ where a straight line from the origin connects with the same gradient must have $`g=0`$, and there is a first-order phase transition between a Meissner phase with no pancake stacks and a finite density of stacks with $`a_{\mathrm{ch}}=a_{\mathrm{ch}}^1`$. All the points at lower density $`B_z<B_z^1`$ have $`g>0`$ and so are thermodynamically unstable. In fact, similar arguments determine that $`f_{\mathrm{ch}}(B_z)`$ must always be a convex function. The dotted line represents the region where the curvature of $`f_{\mathrm{ch}}(B_z)`$ is positive, meaning that the solution here is always unstable.

The straight line that determines $`B_z^1`$ is given by $`f=B\delta H_{\mathrm{c1}}^z/4\pi `$ where $`\delta H_{\mathrm{c1}}^z`$ is the change in the lower critical field due to the attraction of pancake stacks to the Josephson vortices. In Fig. 5b it is given by $`\delta H_{\mathrm{c1}}^z=2.510^3H_{\mathrm{c1}}^z`$ which is hardly measurable. In real experiments however, the geometry of the samples often makes demagnetization effects important such that the average $`B_z`$ becomes fixed to the external field $`H_{\mathrm{ext}}`$, rather than to $`H`$. This means that small values of $`B_z`$ are accessible, and that for $`B_z<B_z^1`$ we may expect a coexistence of $`B_z=0`$ and $`B_z=B_z^1`$ phases (c.f. the intermediate state in type-I superconductorsLandauLifshitz ). Alternatively we could describe this mixed regime in terms of “clusters” of pancake stacks with separation $`a_{\mathrm{ch}}^1`$. The relative proportion of space taken up by these clusters is determined by the value of $`B_z`$, but the size of individual clusters depends on the energy of the “domain wall” between the cluster of stacks and the region of no stacks, compared to the magnetic energy cost of large clusters with the wrong flux density. Also, in the experimental situation one can have different pancake densities on the different chains, so the inhomogenous state may have coexisting empty and filled chains. Our results for the critical field separating the clustered phase from the homogeneous isolated straight chains are shown in Figs. 2 and 3.

## III Buckling Instability within Vortex Chain

It is important to realize that, while the isolated-chain state gains energy due to the Lorentz force of the Josephson vortex currents, there is an energy penalty for the pancake stacks to be close to each other. One consequence of this is a maximum density above which pancake stacks will be ejected from the chain. This critical density will be derived in the next section. Below this density, however, the chain can already react to the stack repulsion by buckling. Note that we will assume that the Josephson vortices remain straight, as in experiments they are held in place by the interactions with Josephson vortices in neighboring chains, and the chain separation is much smaller than the interaction range for Josephson vortices, $`a_{\mathrm{J}y}\gamma \lambda `$. In contrast, the interactions between pancake stacks on different chains are in the dilute limit, $`a_{\mathrm{J}y}\lambda `$, and so are easily displaced.

To calculate this buckling, we need to know the energy gain from a crossing event when the pancake stack is displaced away from the center of a Josephson vortex. We therefore need the full current profile of a Josephson vortex, which is a solution of the non-linear equations that arise from the London-Lawrence-Doniach model.Jvortex ; Bulaevskii92  An accurate numerical solution for a Josephson vortex is described in Appendix B of Ref. Koshelev93, . For a vortex directed along $`\widehat{𝐱}`$ centered at $`y=0`$ and between the $`n=0,1`$ layers, the current in the $`y`$ direction can be writtenKoshelev\_preprint  (ignoring screening, i.e. $`n<\lambda /d`$ and $`y<\lambda _c`$),

$$J_n(y)=\frac{2c\epsilon _0}{\mathrm{\Phi }_0\gamma d}p_n(y/\gamma d),$$
(4)
with $`p_n(\stackrel{~}{y})=\varphi _n^{}(\stackrel{~}{y})`$ the reduced superfluid momentum, where the phase has the form,Koshelev2000

$`\varphi _n(\stackrel{~}{y})`$ $`=`$ $`\mathrm{tan}^1\left[{\displaystyle \frac{(n\frac{1}{2})}{\stackrel{~}{y}}}\right]+{\displaystyle \frac{0.35(n\frac{1}{2})\stackrel{~}{y}}{[(n\frac{1}{2})^2+\stackrel{~}{y}^2+0.38]^2}}`$ (5)
$`+{\displaystyle \frac{8.81(n\frac{1}{2})\stackrel{~}{y}[\stackrel{~}{y}^2(n\frac{1}{2})^2+2.77]}{[(n\frac{1}{2})^2+\stackrel{~}{y}^2+2.02]^4}}.`$
When $`y=0`$, this gives the results used in Section II, with $`J_n(0)=(2c\epsilon _0/\mathrm{\Phi }_0\gamma d)C_n/|n\frac{1}{2}|`$, and $`C_1=0.57`$, $`C_2=0.86`$, $`C_3=0.99`$, and $`C_n1`$ for $`n>3`$.

We can now recalculate the crossing energy for an arbitrary distance between a pancake stack and the center of a Josephson vortex. Within the quadratic approximation used in Refs.Koshelev99,  and Koshelev\_preprint,  one can easily solve the crossing configuration for a single stack. Writing $`\mathrm{\Delta }E_n(u_n)=(\mathrm{\Phi }_0d/c)J_n(y)u_n+\frac{1}{2}\alpha u_n^2`$, where we take $`\alpha =(\epsilon _0d/\lambda ^2)\mathrm{ln}(\lambda /r_w)`$ with $`r_w`$ a short distance cut-off, we find the relaxation energy,

$$\mathrm{\Delta }E_\times (y)=2\left(\frac{\lambda }{\gamma d}\right)^2\frac{1}{\mathrm{ln}(\lambda /r_w)}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\left[p_n(y/\gamma d)\right]^2.$$
(6)
Koshelev takes $`r_w=u_1`$ ($`0.29\lambda ^2/\gamma d`$ for $`y=0`$), and in Ref. Dodgson\_Crete,  we have shown this to be a good approximation to the result when the full form of the pancake-stack interaction is used.

This result for the crossing energy as a function of stack displacement $`u^y`$ is shown as the dashed line in Fig. 6, where there is one Josephson vortex for every 25 pancakes along a stack. As might be expected, the energy increases quadratically at small $`u^y`$. More interesting is the fact that the crossing energy vanishes only as $`\mathrm{\Delta }E_\times \lambda ^2/u_{}^{y}{}_{}{}^{2}`$ for $`u^y\gamma d`$ (note that an exponential suppression will occur at extremely long distances $`u^y>\lambda _c`$). In contrast, also shown in Fig. 6 (the dotted line) for $`a_{\mathrm{ch}}=5.5\lambda `$ is the repulsive energy cost for displacing one pancake stack away from a chain, $`E_{\mathrm{rep}}=_{j0}V_{\mathrm{em}}^{\mathrm{stack}}\left(\sqrt{(ja_{\mathrm{ch}})^2+u_{}^{y}{}_{}{}^{2}}/\lambda \right).`$ This contribution vanishes exponentially when $`u^y>\sqrt{a_{\mathrm{ch}}\lambda }`$, and so the total energy is always negative for large enough $`u^y`$. This gives the possibility of two minima in the total energy (see the full line in Fig. 6), and a first-order transition where the displacement jumps to a finite value as $`a_{\mathrm{ch}}`$ decreases. While we have started by considering displacing a single stack from the chain, it is this instability that drives a buckling of the whole chain.

To estimate the buckling transition we use a variational approach for the buckled configuration (See the inset of Fig. 7a). This configuration is only characterized by a single displacement parameter $`u^y`$, and so we can look at the energy profiles in $`u^y`$ for different $`a_{\mathrm{ch}}/\lambda `$. First, the energy per pancake is calculated before crossing relaxations,

$$E_0^{\mathrm{b}\mathrm{c}}(u^y)=\frac{\epsilon _0d}{N}\underset{ij}{}K_0\left(\sqrt{(ij)^2a_{\mathrm{ch}}^2+(y_iy_j)^2}/\lambda \right)$$
(7)
with $`y_i=0`$ for $`i`$ even, and $`y_i=(1)^{(i+1)/2}u^y`$ for $`i`$ odd. Next we checked for the stability of crossing, as in Section II, and found that the crossing configuration becomes unstable if $`a_{\mathrm{ch}}<2.65\lambda `$. Finally, to calculate the energy gain from crossing the Josephson vortices, we use the result (6) giving the total energy of the buckled chain,

$$E^{\mathrm{b}\mathrm{c}}(u^y)=E_0^{\mathrm{b}\mathrm{c}}(u^y)+\frac{d}{2a_{\mathrm{J}z}}[\mathrm{\Delta }E_\times (0)+\mathrm{\Delta }E_\times (u^y)].$$
(8)
In Fig. 7a we plot this energy for $`a_{\mathrm{ch}}/\lambda `$ from $`5.2`$ to $`6.0`$, with $`a_{\mathrm{J}z}=25d`$, showing how the minimum energy of a buckled chain crosses the energy of a straight chain as $`a_{\mathrm{ch}}`$ decreases. Note that the minimum for a finite $`u^y`$ is quite shallow, and so we might expect large fluctuations in the extent of buckling along the chain due to random vortex pinning, which is discussed in Section VI. Also shown in Fig. 7b is the optimal displacement $`u_{\mathrm{opt}}^y`$ as a function of $`a_{\mathrm{ch}}`$, which has a smallest value of $`u_{\mathrm{opt}}^y=2.3\lambda `$ at the first appearance of (meta)stable buckling at $`a_{\mathrm{ch}}=5.8\lambda `$ and then increases at higher densities.

It has been shown in Section II that one should plot the energy density of the chain as a function of the out-of-plane flux density $`B_z`$, in order to determine the stable thermodynamic phases. Such a plot is shown in Fig. 8 comparing the straight and buckled chains. It shows that the isolated straight chains at $`B_z=B_z^21.3`$ G has the same Gibbs free energy as the buckled-chain state at $`B_z=B_z^{2+}2.5`$ G. Therefore at fields with $`B_z^2<B_z<B_z^{2+}`$ the thermodynamic phase will be a coexistence of straight and buckled regions, with the relative proportion linearly dependent on $`B_z`$. Note that $`B_z^2`$ is only slightly higher than $`B_z^1`$ so that there is only a small regime where pure isolated straight chains are the stable phase. The same procedure has been followed for a range of $`a_{\mathrm{J}z}/d`$, and the resulting boundaries to the mixed buckled chains are shown in Figs. 2 and 3.

## IV Ejection of Pancake Stacks from Isolated Chain

We now consider a simple model for the energy of the composite state where a fraction of the pancake stacks are not located on the chains. This model will show that there is a continuous phase transition between this composite state and the isolated-chain state as a function of $`B_z`$. The order parameter of this transition is the density of the dilute lattice of pancake stacks not trapped on chains. The simple model assumes that the chain spacing $`a_{\mathrm{J}y}`$, and the spacing between pancake stacks is much greater than the penetration depth $`\lambda `$. In this case we can separate the total interaction energy density of pancake stacks in the form $`f_{\mathrm{tot}}=f_{\mathrm{ch}}(\nu _{\mathrm{ch}})+f_{\mathrm{dil}}(\nu _{\mathrm{dil}})`$, where $`\nu _{\mathrm{ch}}`$ and $`\nu _{\mathrm{dil}}`$ are the (2D) densities of pancake stacks on and off the chains respectively, and the total density is $`\nu _{\mathrm{ch}}+\nu _{\mathrm{dil}}=\nu _{\mathrm{tot}}=B_z/\mathrm{\Phi }_0`$. The separation of pancakes along the chain is $`a_{\mathrm{ch}}=1/a_{\mathrm{J}y}\nu _{\mathrm{ch}}`$. The chain energy $`f_{\mathrm{ch}}(\nu _{\mathrm{ch}})`$ is to be calculated as in Sections II and III for straight and buckled chains respectively. We will see that the phase transition to the composite state occurs when there is a minimum in $`f_{\mathrm{ch}}(\nu _{\mathrm{ch}})`$, and so near this transition we will expand,

$$f_{\mathrm{ch}}(\nu _{\mathrm{ch}})=ϵ^{}\lambda ^2(\nu _{\mathrm{ch}}\nu _{\mathrm{ch}}^\mathrm{c})^2.$$
(9)
The density of pancake stacks not on a chain is small near the transition, so that the interaction energy density is simply (using the limit $`K_0(x)\sqrt{\pi /2x}e^x`$ and only including nearest neighbors),

$$f_{\mathrm{dil}}(\nu _{\mathrm{dil}})=3\epsilon _0\sqrt{2\pi \lambda }\nu _{\mathrm{dil}}^{5/4}e^{1/\lambda \nu _{\mathrm{dil}}^{1/2}}.$$
(10)
Note that in the small $`\nu _{\mathrm{dil}}`$ limit, all terms in a power series expansion of $`f_{\mathrm{dil}}(\nu _{\mathrm{dil}})`$ are zero, i.e. the function is extremely flat. For this reason the critical total density at the transition is only determined by the minimum in $`f_{\mathrm{ch}}(\nu _{\mathrm{ch}})`$ at $`\nu _{\mathrm{ch}}^\mathrm{c}`$. For $`\nu _{\mathrm{tot}}<\nu _{\mathrm{ch}}^\mathrm{c}`$, all of the pancake stacks are on chains. At densities just above $`\nu _{\mathrm{ch}}^\mathrm{c}`$ the energy is minimized with a dilute off-chain density,

$$\nu _{\mathrm{dil}}=\delta \nu _{\mathrm{tot}}\frac{3\sqrt{2\pi }}{4}\frac{\epsilon _0}{\epsilon ^{}}\frac{1}{\lambda ^{5/2}\delta \nu _{\mathrm{tot}}^{1/4}}e^{1/\lambda \delta \nu _{\mathrm{tot}}^{1/2}}.$$
(11)
The second term is extremely small for the fields we are interested in $`B_z\mathrm{\Phi }_0/\lambda ^2`$, so that we can say in the composite state at $`\nu _{\mathrm{tot}}>\nu _{\mathrm{ch}}^\mathrm{c}`$ there is a fixed density on the chains $`\nu _{\mathrm{ch}}=\nu _{\mathrm{ch}}^\mathrm{c}`$ and the remaining density is $`\nu _{\mathrm{dil}}=\nu _{\mathrm{tot}}\nu _{\mathrm{ch}}^\mathrm{c}`$. The phase boundary $`\nu _{\mathrm{ch}}^\mathrm{c}`$ between the composite state and the isolated-chain state is calculated from the minimum of $`f_{\mathrm{ch}}`$. In Figs. 2 and 3 this boundary is plotted, but is not distinguishable by eye from the transition between mixed buckled and pure buckled chain. Therefore the pure buckled chain has a very narrow range of existence. It is worth stating that an ejection transition has been observed experimentallyGrigorenko  with results that seem to be consistent with the above, although a quantitative analysis of the transition has not yet been published.

## V Crossing-to-Tilted transition within isolated chains

In this section we calculate the energy of a chain of tilted vortices, and compare to the chain of crossing vortices. The tilted chain is an alternative configuration to the crossed chain with the same density of pancake and Josephson vortices (see Fig. 9) where the tilt angle is determined by $`\mathrm{tan}\theta =L/d=a_{\mathrm{ch}}/a_{\mathrm{J}z}`$. These tilted stacks have an increased energy of pancake interactions but there may be an energy gain from replacing fully developed Josephson vortices with many short Josephson segments of length $`L<\gamma d`$. To calculate the energy of the tilted chain, we use the linear approximation of the London-Lawrence-Doniach model,Bulaevskii92 ; Koshelev93  where we ignore the non-linear effects from regions with large phase differences across neighboring layers. This approximation is justified as long as the length $`L`$ is much smaller than the Josephson length $`\gamma d`$ and the pancake separation $`a_{\mathrm{ch}}`$. Note that the condition $`L<\gamma d`$ corresponds to $`B_x/B_z<\gamma `$, which is the case we are interested in.

The energy due to pancake interactions in the linear approximation is,

$$E_{\mathrm{pc}}^{\mathrm{tilt}}=\frac{1}{8\pi }\frac{d^2k}{(2\pi )^2}\frac{dq}{2\pi }\frac{\left[1+(\gamma \lambda )^2(k^2+\stackrel{~}{q}^2)\right]|S_z|^2}{\left[1+\lambda ^2(k^2+\stackrel{~}{q}^2)\right]\left[1+(\gamma \lambda )^2k^2+\lambda ^2\stackrel{~}{q}^2\right]}$$
(12)
with $`S_z(𝐤,q)=\mathrm{\Phi }_0d_{j,n}e^{iqnd}e^{ik_x(ja_{\mathrm{ch}}+nL)}`$ and $`\stackrel{~}{q}^2=(2/d)^2\mathrm{sin}^2(qd/2)`$. The integral should be cut off due to the pancake cores at $`k>\pi /\xi _{\mathrm{ab}}`$ and due to the layered structure at $`q>\pi /d`$. Similarly, the Josephson segments interact with an energy contribution,

$$E_{\mathrm{Jv}}^{\mathrm{tilt}}=\frac{1}{8\pi }\frac{d^2k}{(2\pi )^2}\frac{dq}{2\pi }\frac{|S_x|^2}{\left[1+(\gamma \lambda )^2k^2+\lambda ^2\stackrel{~}{q}^2\right]}$$
(13)
with $`S_x=\frac{2\mathrm{\Phi }_0}{k_x}\mathrm{sin}\left(\frac{k_xL}{2}\right)_{j,n}e^{iq(n+\frac{1}{2})d}e^{ik_x[ja_{\mathrm{ch}}+(n+\frac{1}{2})L]}`$.

The integrals can all be done exactly, leaving a sum over the values of $`k_x=2\pi m/a_{\mathrm{ch}}`$. We have evaluated these sums numerically to find the energy of the tilted chain. This is shown for $`a_{\mathrm{J}z}=25d`$ as the dotted line in Fig. 10. Also shown as the full line is the energy per pancake of the crossing chain with the same density of pancakes and Josephson vortices, $`E_{\mathrm{cross}}=E_{\mathrm{ch}}+E_{\mathrm{Jv}}^{\mathrm{cross}}`$, where $`E_{\mathrm{ch}}`$ is defined in Section II and

$$E_{\mathrm{Jv}}^{\mathrm{cross}}=\frac{\epsilon _0d}{\gamma }\frac{a_{\mathrm{ch}}}{a_{\mathrm{J}z}}\left[1.55+\mathrm{ln}\left(\frac{\lambda }{d}\right)+\underset{j0}{}K_0\left(\frac{ja_{\mathrm{J}z}}{\lambda }\right)\right].$$
(14)
The tilted chain has a lower energy from Josephson vortex contributions (dominant at low density) but a higher energy from pancake vortex contributions (dominant at high density), and there is a first-order transition between the isolated chain of straight stacks and the chain of tilted stacks. In fig. 10 it is shown that for $`a_{\mathrm{ch}}<45\lambda `$ (higher density) the preferred state is the crossing chain, while for $`a_{\mathrm{ch}}>45\lambda `$ (lower density), the lowest energy is for the tilted chain. Note that Savelev et al.,Savelev  also find a “reentrant” transition back to a tilted lattice for fields close to the $`ab`$-plane, when $`\lambda >\gamma d`$.

A strict treatment should follow Sections II and III and consider the points of the same Gibbs free energy density. This is not possible within our level of treatment, however, as we would need accurate results at very low $`B_z<B_x/\gamma `$ where non-linear effects in the tilted chain are important. Here, we can only show that there must be a low-density transition to the tilted chain, but we cannot plot the coexistence boundaries. Even so, we plot the line where the energies cross in Figs. 2 and 3 separating the clustered chain from a tilted chain phase. In real experiments there is a finite 2D density of Josephson vortices, which cannot be considered isolated in the sense of the pancake vortex chains (the interaction range $`\lambda _c`$ is greater than the chain spacing $`a_{\mathrm{J}y}`$). Therefore a full calculation of the 3D tilted lattice requires a minimization over the aspect ratio, which we leave for a future work.

## VI Discussion

We now discuss the extent to which these transitions have been observed experimentally. In the recent scanning Hall probe experiments of Grigorenko et al.,Grigorenko  there does seem to be a sharp transition as a function of $`B_z`$ between a composite state of chains with a dilute lattice and the isolated chain state, as derived in Section IV. At lower fields there is some evidence for buckling and clustering (Sections II and III) although the influence of pinning disorder may be contributing to this. Finally, Grigorenko et al. report a strange transition at very low $`B_z`$ where the chains are replaced by faint, homogeneous “stripes” of flux. These stripes would be consistent with the flux distribution from a tilted chain, and it is likely that a transition from the isolated-chain state to a tilted-chain state (as found in Section V) takes place. Experimentally, a large jump in $`B_z`$ (i.e. in the pancake density) is seen at the transition to tilted chains; unfortunately we were not able to calculate the jump in this paper, as it requires an accurate treatment of non-linear effects in the small $`B_z`$ limit. Other experimental features are that the tilted chains are much straighter than the crossing chains, and less pinned. This is consistent with the fact that the tilted chains are really part of a 3D line lattice with stronger interactions that dominate over pinning disorder, whereas in the crossing case the pancake stacks are in the extreme dilute limit, and therefore easily pinned in random low-energy sites.

Finally, we briefly discuss the effect of fluctuations. The calculations in this paper have found the lowest-energy vortex configurations of the chains. In reality there will be distortions of these states due to thermal fluctuations and random pinning to inhomogeneities in the underlying crystal. It is well established that thermal fluctuations can lead to a melting of the vortex crystal.Blatterreview  A melting transition to completely decoupled layers is studied for $`B_x=0`$ in Ref. Dodgson\_Substrate, , involving short-wavelength fluctuations of pancake stacks with $`k_z>\pi /\lambda `$. While there will be some modification in the presence of Josephson vortices in tilted fields, in the low out-of-plane fields studied here, $`B_z\mathrm{\Phi }_0/\lambda ^2`$, this melting will take place at a high temperature, close to the Berezinskii-Kosterlitz-Thouless transitionBKT  at $`T_{\mathrm{BKT}}=\epsilon _0d/2`$. If we consider long wavelength fluctuations, $`k_z0`$, then we may expect some kind of melting of the chain at very low fields due to the exponentially small interaction of pancake stacks. However, the mechanism of such a transition must be different from the “entanglement” proposed in the original derivation of low-field vortex lattice melting by Nelson.Nelson  Considering the nature of the crystalline order of the isolated vortex chain, we should recognize the two-dimensional nature of this state. A simple consequence is that thermal fluctuations at long wavelength will lead to a quasi-long-range ordered state. It may then be possible to have a 2D continuous melting transition via the unbinding of dislocations. However, a dislocation for this system of stacks corresponds to a stack that terminates in the middle of the chain, and this should cost an energy linear in the system size, rather than the usual logarithmic energy of a dislocation. This will suppress a dislocation unbinding transition.

In the experiments of Grigorenko et al. the pinning disorder seems to be more important in disturbing the chain states than thermal fluctuations. The fact that the pancake stacks are observed in fixed positions is due to pinning (otherwise thermal fluctuations would smear out the average density in the chain state), and this also tends to disorder the chains. On general grounds we expect pinning-induced wandering of the pancake stacks within the chains so that there is only a short range order. There will also be significant transverse displacements, and the buckling effect should be enhanced by disorder.

It is a pleasure to thank Simon Bending and Sasha Grigorenko for enthusiastically introducing me to the world of chains and crossing lattices, as well as Dima Geshkenbein for useful discussions and Alex Koshelev for demonstrating how to calculate the currents near a Josephson vortex. The author is supported by an EPSRC Advanced Fellowship AF/99/0725. Work at the ETH-Zürich was financially supported by the Swiss National Foundation.


format: text

text:
# A microscopic model for thin film spreading

\[

## Abstract

A microscopic, driven lattice gas model is proposed for the dynamics and spatio-temporal fluctuations of the precursor film observed in spreading experiments. Matter is transported both by holes and particles, and the distribution of each can be described by driven diffusion with a moving boundary. This picture leads to a stochastic partial differential equation for the shape of the boundary, which agrees with the simulations of the lattice gas. Preliminary results for flow in a thermal gradient are discussed.

\]

Spreading of involatile liquid drops on surfaces and fibers plays a significant role in many technologies, the efficient deployment of which requires detailed understanding of the underlying process . Behavior at a macroscopic scale is accurately described by hydrodynamics , but experiments performed over the last 15 years or so have produced a surprise: in the case of complete wetting, the spreading drop, when examined on an atomic length scale in the direction normal to the substrate, is found to be preceded by a precursor film about one molecule thick; this can be followed laterally out to an extension of the order of $`10^7`$ molecules diameter. At the horizontal resolving power of the technique used (ellipsometry), the precursor film is found to be flat and homogeneous. Its radius advances with time as $`\sqrt{t}`$. Typical systems showing this class of behavior which can be investigated by ellipsometry are various silanes spreading over atomically flat Si(111) wafers with highly pure oxydised surfaces. With selected silanes, such systems even show dynamical layering with up to four superposed precursor films advancing as $`\sqrt{t}`$, the layer directly in contact with the substrate being much faster than the others. As Ball observed over a decade ago , at that time there was hardly even a formative theory of such phenomena. Such theory must capture the experimentally determined diffusive behavior and, at the same time, explain how an extremely viscous, involatile material is transported from the reservoir to the precursor edge.

A partial solution was first achieved by applying Molecular Dynamics to a droplet composed of spherical molecules with Lennard-Jones interactions . The strength of these is adjusted to achieve sufficient involatility. Experimentally, the droplet is placed in contact with a substrate composed of much smaller units; this is therefore treated as a continuum for calculating the interaction of the spreading molecules with the substrate. Such a model shows a precursor film with the correct, diffusive behavior . Monte Carlo (MC) lattice gas simulations in 3d  confirm this, and suggest a dual mechanism of matter transport involving on the one hand particles located in a supernatant layer directly placed above the precursor film, and on the other, holes in the precursor film itself. At the same time, essentially all the configurational change occurs at the boundary of the drop. However, so far these ideas just constitute a picture without quantitative implications.

In this Letter, we introduce an associated 2d driven Ising lattice gas model which can be simulated far more effectively. Extensive continuous time MC simulations are consistent with the $`\sqrt{t}`$ behavior found in all the experiments to which we refer . They also allow us to formulate a probabilistic continuum model of independent hole and particle diffusion in the supernatant and precursor films, which accounts for the matter transfer between the reservoir and the precursor edge (PE). The motion of the latter in the driven lattice gas is then equivalent to a moving boundary in the continuum model. Starting from the constitutive equations of the moving boundary problem, we derive the effective dynamics of the PE as a stochastic partial differential equation. We thus obtain analytical predictions which are consistent both with experiments and with the original 2d Ising model, and which can be extended to the asymptotic regime thus far unexplored both theoretically and experimentally.

We now describe our MC lattice gas simulations for this problem. At each cell on a square lattice $`𝒓=(x,y,z)`$ with unit lattice spacing, we define the ocuppation number $`n(𝒓,t)`$ to be $`0`$ or $`1`$ depending whether the cell is occupied or not. Moreover, we restrict $`z`$ to the values $`z=1`$ (precursor film) and $`z=2`$ (supernatant layer). The interaction energy for a given configuration is

$$=J\underset{|𝒓𝒔|=1}{}n(𝒓,t)n(𝒔,t)A\underset{z=1,2}{}\frac{1}{z^3}n(𝒓,t),$$
(1)
where the sums are over $`𝒓,𝒔`$ on the lattice. The last term is the van der Waals interaction with the continuum substrate characterized by a Hamaker constant $`A>0`$. We are interested in $`J/k_BT`$ large enough to achieve a high degree of involatility and $`A/k_BT`$ large enough to be in the complete wetting regime . Continuous time MC Kawasaki dynamics are used  with initial conditions: (i) $`n(𝒓,0)=1`$ for $`y=1`$, but with any $`x`$ and for both $`z=1`$ and $`z=2`$; (ii) $`n(𝒓,0)=0`$ for $`y>1`$. If, as a result of an allowed MC move at any later stage, $`n(x,1,z,t)=0`$, then the implied vacancy is filled immediately from the putative reservoir droplet, $`n(x,1,z,t)=1`$ being maintained in this sense for all $`t>0`$. Thus, the detailed evolution of the reservoir droplet is ignored. In reality, it would equilibrate relatively rapidly in profile by motion of particles near its surface and deflate by emission of supernatant particles and absorption of precursor holes . The main advantage of the boundary condition we use, and a very significant one, is in computability. We believe our model to be new — it is two superposed Ising lattice gases driven from one edge.

A typical snapshot of the simulations is shown in Fig. 1. There is a compact film advancing in both the $`z=1`$ and $`z=2`$ layers, but the former (shown black) advances much faster that the latter (shown gray). It is easy to understand why this should be so: if a hole in the precursor film diffuses to lie beneath the compact film in $`z=2`$, at an appropriately selected update it will be filled by creating a hole in film above, which diffuses and may be trapped at the edge $`y=1`$, or at the supernatant film boundary, causing this film to shrink. The precursor advances because of emission of holes into itself, but also from the arrival of supernatant particles at its edge; such particles drop into the $`z=1`$ plane with concomitant Hamaker stabilization .

To study the dynamics of the PE, definitions of spin-percolative type are needed. Thus, the precursor film is the particle cluster $`𝒫(t)`$ in $`z=1`$ which is connected along nearest neighbor bonds to the boundary line $`y=1`$. The PE is at $`y=h(x,t)`$ for $`x=1,\mathrm{},L_x`$, where $`h(x,t)`$ is given by the maximum value of $`y`$ among those cells $`(x,y)`$ that belong to $`𝒫(t)`$. In Fig. 2 we present results for the dynamics of the PE: the velocity $`V(t)=d\overline{s}/dt`$ \[panel (a)\], where $`\overline{s}(t)`$ is the the mean PE displacement $`\overline{s}(t)=(1/L_x)_{x=1}^{L_x}h(x,t)`$, and the roughness $`w^2(L_x,t)=(1/L_x)_{x=1}^{L_x}\zeta ^2(x,t)`$ \[panel (b)\], where $`\zeta (x,t)h(x,t)\overline{s}(t)`$ and $`\mathrm{}`$ stands for average over different MC runs.

After a short initial transient $`t=\tau _010^3`$ MC units (independent of $`L_x`$), we recover the law $`\overline{𝒔}(t)=\gamma t^{1/2}`$, with $`\gamma `$ a time independent constant. This is in clear qualitative agreement with the experimental results, and the confirmation of the universal $`\sqrt{t}`$ law is considerably more precise than in the 3-d simulations . On the other hand, for times longer than the initial transient $`t>\tau _0`$ —before which the roughness oscillates due to the discreteness of the model and the flat initial condition on the PE—, the PE displays kinetic roughening as reflected in the behavior $`w(L_x,t)t^\beta `$ . The growth exponent $`\beta `$ seems to cross over from the value $`\beta =1/6`$ at intermediate times to $`\beta =1/8`$ at long times. Due to the restricted range of computationally accessible values of $`L_x`$, the roughness exponent $`\alpha `$ in the stationary state relation $`w(L_x,\mathrm{})L_x^\alpha `$  is more accurately extracted from the behavior of the structure factor $`S(k,t)=\zeta _k(t)\zeta _k(t)`$, where $`\zeta _k(t)`$ is the Fourier transform of $`\zeta (x,t)`$. For a rough interface, $`S(k,t)`$ is expected  to behave as $`S(k,t)k^{(2\alpha +1)}`$ for long enough times. In Fig. 3(a) we show the structure factor of the PE measured for various times. As is clear from the figure, $`S(k,t)`$ tends to a scale invariant behavior at long times, characterized by a roughness exponent $`\alpha =1/2`$.

Inspection of typical runs shows a rather sparse distribution of particles \[concentration $`c_\mathrm{p}(𝒓,t)`$\] and holes \[concentration $`c_\mathrm{h}(𝒓,t)`$\], which attenuates very rapidly away from their respective sources. Since the concentrations are typically low, the dynamics suggest to introduce an independent and non-interacting particle-hole diffusion model, with identical equations for particles and holes. We will treat the $`𝒓`$ coordinate as a continuous variable. Thus, on $`𝒫(t)`$, we have

$$_tc_\mathrm{h}=D\mathrm{\Delta }c_\mathrm{h}𝒒,$$
(2)
where $`𝒒(𝒓,t)`$ is a noise term representing fluctuations in the hole diffusion current within $`𝒫(t)`$ . The interface velocity along the outward normal $`𝒏`$ is given by

$$V_n=2(Dc_\mathrm{h}𝒒)𝒏B_s^2\kappa _s𝒑.$$
(3)
The first term on the right in (3) comes from diffusion of holes away from the PE. The prefactor two reflects the dual rôle of particles and holes and that both have the same diffusion constant. The second term represents surface diffusion due to the tendency for particles/holes to move along the interface from regions of negative to positive curvature, the constant $`B`$ being proportional to the interface concentration of mobile species. Here $`\kappa `$ is the mean curvature. Finally, $`𝒑(𝒓,t)`$ is a noise term describing fluctuations in the interface diffusion current along the precursor edge . The concentration obeys also the boundary condition $`c_\mathrm{h}(x,0,t)=0`$ and the Gibbs-Thompson condition

$$c_\mathrm{h}(𝒓,t)|_{y=h(x,t)}=c_0+\mathrm{\Gamma }\kappa ,$$
(4)
where $`c_0`$ is the average concentration of holes at the precursor edge and $`\mathrm{\Gamma }=c_0\sigma /k_BT`$, where $`\sigma `$ is the surface tension, which to leading order is taken to be isotropic.

This moving boundary problem presents formidable obstacles to direct attack. Neglecting the noise terms $`𝒒`$ and $`𝒑`$, assuming a flat interface, and averaging over $`x`$, a classical Stefan problem  is recovered for the mean hole concentration $`\overline{c}(y,t)`$. On $`𝒫(t)`$, $`\overline{c}(y,t)`$ satisfies the diffusion equation $`_t\overline{c}=D_y^2\overline{c}`$, but now with far simpler boundary conditions $`\overline{c}(y,t)|_{y=\overline{s}(t)}=c_0`$, $`\overline{c}(0,t)=0`$ and

$$\frac{d\overline{s}}{dt}=2D_y\overline{c}|_{y=\overline{s}(t)}.$$
(5)
This problem admits a similarity solution  in which $`\overline{s}(t)=\gamma t^{1/2}`$, and

$$\overline{c}(y,t)=\frac{\gamma }{4D}_0^{y/t^{1/2}}e^{(\gamma ^2u^2)/4D}𝑑u,$$
(6)
where $`\gamma `$ is the unique solution of

$$c_0=\frac{\gamma }{4D}_0^\gamma e^{(\gamma ^2u^2)/4D}𝑑u.$$
(7)
Equation (7) then gives a relation between $`\gamma `$ and $`c_0`$. Now, both $`\gamma `$ and $`c_0`$ can be read off from the simulations, with different choices of temperature giving different $`(c_0,\gamma )`$ points which are compared in Fig. 3(b) with the predictions of (7). The way the simulations of the lattice gas are set up implies that, for comparison with the diffusive system, we should have $`D=1/3`$. The excellent agreement in Fig. 3(b) furnishes very significant support for our continuum model.

The spatio-temporal fluctuations of the PE, expressed by $`w(L_x,t)`$ and $`S(q,t)`$ in Figs. 2(b) and 3(a), respectively, require full implementation of the boundary conditions (3) and (4), a formidable mathematical problem. Progress can be made when $`\gamma `$ is small \[see Fig. 3(b)\] by (i) setting $`_tc_\mathrm{h}=0`$ for $`y<h(x,t)`$, the quasistationary approximation ; (ii) perturbing around the spatially-averaged solution given by (6) and (7) : we expand Eq. (3) to quadratic order in the fluctuation variable $`\zeta (x,t)=h(x,t)\overline{s}(t)`$ describing the precursor edge —taking the noise terms to be of the same order as $`\zeta (x,t)`$ —, to obtain

$`{\displaystyle \frac{\zeta _k}{t}}`$ $`=`$ $`\left[D|k|\left(\mathrm{\Gamma }k^2+{\displaystyle \frac{c_0}{\gamma t^{1/2}}}\right)\mathrm{coth}(|k|\gamma t^{1/2})+Bk^4\right]\zeta _k`$ (9)
$`+{\displaystyle \frac{V(t)}{2}}_k[(_x\zeta )^2]+\eta _k(t),`$
where $`_k[f(x)]`$ is the Fourier transform of $`f(x)`$, and $`\eta _k(t)`$ is a white noise with correlations

$`\eta _k(t)\eta _k^{}(t^{})`$ $`=`$ $`2c_0\delta _{k,k^{}}\delta (tt^{})`$ (10)
$`\times `$ $`\left\{D\left[1+{\displaystyle \frac{|k|^2\gamma t^{1/2}}{\mathrm{sinh}^2(|k|\gamma t^{1/2})}}\right]+{\displaystyle \frac{B}{\mathrm{\Gamma }}}k^2\right\}.`$ (11)
In Eq. (9) the linear terms stabilize PE fluctuations, while the nonlinear term is of the Kardar-Parisi-Zhang (KPZ) type , with a coefficient $`V(t)=\gamma /t^{1/2}`$ that decays in time. In the asymptotic regime $`t\mathrm{}`$, and neglecting terms proportional to $`c_0`$, $`\mathrm{\Gamma }𝒪(\gamma ^2)`$, Eqs. (9) and (11) reduce to the conserved linear equation

$$\frac{\zeta _k}{t}=Bk^4\zeta _k+\eta _k,$$
(12)
with noise correlations $`\eta _k(t)\eta _k^{}(t^{})=2(c_0B/\mathrm{\Gamma })k^2\delta _{k,k^{}}\delta (tt^{})`$. It is easy to solve this equation exactly and obtain $`\beta =1/8`$ and $`\alpha =1/2`$, which provide accurate fits at long times for the behavior of the roughness \[Fig. 2(b)\] and the structure factor \[Fig. 3(a)\], respectively. Moreover, the behavior of $`S(k,t)`$ at long but finite times \[$`e.g.`$, $`t=10^6`$ MC units in Fig. 3(a)\] is also consistent with that predicted by Eq. (12). For sufficiently long times, the rate of arrival of particles becomes essentially zero. In this case, emission of particles to the exterior of the PE and holes to the interior are in balance, so the mean interface is stationary. In such a situation, up to order $`𝒪(\gamma ^2)`$, the only fluctuations on the PE are due to the conserved surface diffusion current, which is precisely described by Eq. (12). However, for intermediate times ($`e.g.`$, $`t10^5`$ MC units), the behavior of the fluctuations deviates from that predicted by Eq. (12), see Figs. 2(b) and 3(a). A plausible explanation is that the KPZ term in Eq. (9) is still non-negligible, in which case it should dominate the large scale properties of the PE. Using results for the KPZ equation , the roughness should then be given in this time regime by $`w[V(t)t]^{1/3}t^{1/6}`$, again in rather good agreement with the numerical results for the lattice gas, see Fig. 2(b).

In summary, using simulations of our new driven lattice gas model we have (i) recaptured the $`\sqrt{t}`$ universal mean precursor displacement law; (ii) shown that the spatio-temporal fluctuations of the PE demonstrate dynamical scaling; (iii) shown that the mechanisms of matter transport from the reservoir droplet to the PE involve low density of particles and holes, described by diffusion equations with moving boundaries; (iv) derived a stochastic partial differential equation for the PE which brings in the fluctuation regimes found previously in a natural way, and which tends asymptotically to a Langevin description for a free interface with conserved dynamics as in (12), possibly observable experimentally, even though it presents a KPZ behavior at intermediate times. In preliminary work, we have investigated our lattice gas model for a single layer ($`z=1`$) in a thermal gradient, finding a mean spatial displacement growing linearly in time, generated by biased hole diffusion. This has relevance for microfluidics .

We thank A. Sánchez for participation in the early stages of this work, and M. Castro for discussions. E.M. acknowledges the EU fellowship No. HPMF-CT-2000-0487. This research has been supported by EPSRC (UK) Grant No. GR/M04426 and by DGES (Spain) grants Nos. HB1999-0018 and BFM2000-0006.


text:
# KUV 01584-0939: A Helium-transferring Cataclysmic Variable with an Orbital Period of 10 MinutesThis paper uses observations made from the South African Astronomical Observatory (SAAO).

## 1 INTRODUCTION

The star KUV 01584-0939 (hereafter referred to as Cet3, as in the Downes, Webbink & Shara (1997) catalogue of Cataclysmic Variables (CVs)) was discovered during the Kiso survey (Kondo, Noguchi & Maehara 1984) for ultraviolet rich objects. A spectrum of Cet3 covering the wavelength range 4000 – 7100 Å was obtained by Wegner, McMahan & Boley (1987) who drew attention to the great strength of He II emission and the weakness of the Balmer emission lines. We have observed Cet3 as part of a high time resolution photometric survey of faint CVs (Woudt & Warner 2001).

## 2 OBSERVATIONS AND RESULTS

Our high speed photometric observations were made with the University of Cape Town CCD photometer (O’Donoghue 1995) mounted on the 74-in Radcliffe reflector at the Sutherland site of the South African Astronomical Observatory on the nights of 18 – 21 October 2001. Cet3 is quite faint ($`m_V16.9`$) so we used ‘white light’, counting all photons detected by the CCD. The light curves for the three long observing runs, obtained on consecutive nights, are shown in Figure 1.

From our observations we find that Cet3 is a variable star with a period of 620.26 $`\pm `$ 0.02 secs and a peak-to-peak range of $``$0.15 mag. The Fourier transform (FT) of our brightness measurements shows a fundamental and its harmonics. The presence of harmonics implies that the modulation is not sinusoidal. In Figure 2 we present a compact view of the FTs and the mean pulse profiles for each of the 4 observing runs. The mean pulse profiles are phased according to the ephemeris given in Eqn. 1. There are only small changes in average profile from night to night, but individual pulse profiles vary as though there is some rapid flickering present. In Figure 1 it can also be seen that slow changes in mean brightness can also occur. We have found no persistent features in the FTs other than the 620 s fundamental and its harmonics. The ephemeris for maximum light of Cet3, derived from a non-linear fit of sinusoids to the fundamental and the first two harmonics, is

$$\mathrm{HJD}_{\mathrm{max}}=2452201.391401+0\stackrel{\mathrm{d}}{\mathrm{.}}0071788(\pm 2)\mathrm{E}.$$
$`(1)`$
Photometric periods of 102 – 103 s in stellar sources can originate from a number of different structures: including, e.g., orbital periods of low mass X-ray binaries (Bradt & McClintock 1983) containing a neutron star, or of the AM CVn stars mentioned above, or the rotation periods of magnetic white dwarfs in CVs (known as intermediate polars (Warner 1995a)). We can exclude a number of other possibilities (e.g., pulsations) because of the exceptional nature of the spectrum of Cet3. The latter also helps us to select from among the orbital or rotational options.

The most striking feature of the optical spectrum (Wegner et al. 1987) of Cet3 is the strength of the He II lines at 4026 and 4686 Å, and their presence at 4100, 4199 and 5412 Å. The C IV emission line at 5802 Å and the very blue continuum emphasize that the source is very hot. These identifications raise the possibility that the weak apparent Balmer emission is not from H I: it is more probably the Pickering series of He II, in which case Cet3 is hydrogen deficient, probably belongs to the AM CVn class of helium transferring interacting binaries, and the very blue continuum is that of an accretion process. No low mass X-ray binary or intermediate polar has a spectrum resembling that of Cet3. The slow changes in mean brightness and the presence of flickering are characteristics of a mass-transferring close binary system.

## 3 DISCUSSION

AM CVn stars evolve from normal composition binaries which have experienced two phases of mass exchange, exposing the helium cores of the original stars. During their orbital evolution such systems pass through a minimum period near 4 mins and thereafter evolve to longer periods (Tutukov & Yungelson 1996). The driving mechanism for orbital evolution in such short period binaries is loss of angular momentum through emission of gravitational radiation (Paczynski 1967). This determines the rate of mass transfer, $`\dot{M}`$, and predicts (Warner 1995b) values $`1\times 10^8`$ M y-1 at $`P_{orb}`$ = 10 min, decreasing steeply to $`4\times 10^{11}`$ M y-1 at $`P_{orb}`$ = 40 min. These rates are compatible with the range of observed properties of the seven optically observed (Warner 1995b; Ruiz et al. 2001) AM CVn stars and the X-ray system (Cropper et al. 1998). The latter (V407 Vul = RX J1914+24) has $`P_{orb}`$ = 569 s, is thought to be a polar (i.e., the primary is strongly magnetic, so its rotation is synchronised with the revolution of the secondary), but is very faint optically because it is obscured by an interstellar cloud. In contrast, Cet3 has a similar period and is relatively accessible at $`m_V=16.9`$.

If Cet3 has $`\dot{M}1\times 10^8`$ M y-1, its accretion luminosity will be $``$ 5 L for a 0.7 M white dwarf primary. Most of this will be radiated at the inner boundary layer of the accretion disc at a temperature of $`3\times 10^5`$ K, so Cet3 should be a soft X-ray source of high intrinsic luminosity. This also accounts for the strength of the helium and carbon ionic spectra in Cet3. We calculate that d$`P_{orb}`$/d$`t`$ $`6\times 10^{12}`$, which implies that d$`P_{orb}`$/d$`t`$ should be measurable in only a few years of observations if the modulation indeed has an orbital origin. We have no means at present to estimate the distance of Cet3, but we expect it to be a strong source of gravitational radiation (Warner 1995b; Hils & Bender 2000).

We can estimate parameters for Cet3 in the following way. As the secondary must fill its Roche lobe, a $`P_{orb}`$ of 620 s implies a mean density $`\overline{\rho }3.6\times 10^3`$ gm cm-3 (see equation 2.3b of Warner 1995a). From the mass-radius relation for low mass white dwarfs (Tutukov & Yungelson 1996), this gives a mass for the secondary $`M(2)`$ = 0.066 M. For an assumed primary mass $`M(1)`$ = 0.7 M, we have $`q=M(2)/M(1)`$ = 0.094 and a binary separation of $`9.9\times 10^9`$ cm. Cet 3 is a very compact system: from standard formulae (Warner 1995a) we find an accretion disc radius $`5.4\times 10^9`$ cm, the radius $`R(2)`$ of the secondary $`2.1\times 10^9`$ cm, and the radius of the primary $`R(1)0.85\times 10^9`$ cm. The accretion disc therefore extends only $`5.4R(1)`$ above the surface of the primary and will be strongly irradiated by the $`3\times 10^5`$ K source at the inner edge of the disc. In addition, for $`\dot{M}1\times 10^8`$ M y-1 the maximum temperature in the disc itself is $`6.2\times 10^4`$ K. The surface of the secondary intercepts $``$ 1% of the soft X-ray emission and, with allowance for shielding by the accretion disc, the side facing the primary will as a result be heated to $`1.8\times 10^4`$ K.

These parameters are of use in attempting to understand the modulation profile in Cet3. For a very small mass-ratio ($`q`$ = 0.094) the accretion disc nearly fills the Roche lobe of the primary (Warner 1995a) and partial eclipses of the disc are quite probable. The shallow but wide minima in the Cet3 light curve could therefore be eclipse features. The profile could alternatively be a combination of reflection effects from the secondary and aspect effects of the disc and bright spot. However, a small mass ratio $`q`$ and high $`\dot{M}`$ usually result in perturbation of the accretion disc into an elliptical form with consequent precession (Warner 1995a). The AM CVn stars with high $`\dot{M}`$ show ‘superhumps’ with periods a few percent longer than $`P_{orb}`$, which result from the tidal stresses in a precessing disc (Patterson et al. 2001). It is probable, therefore, that the modulation in Cet3 is a permanent superhump, in which case $`P_{orb}`$ will be slightly shorter than the period we have measured. Alternatively, Cet3 may be a polar like V407 Vul (the absence of two distinct periods disfavours an intermediate polar interpretation). It should be possible to distinguish between these various interpretations when high time resolution spectroscopy is available.

## ACKNOWLEDGEMENTS

We thank Dr. D. Kilkenny for helpful comments. This research is supported by grants from the University of Cape Town.


text:
# An attempt to model globular cluster red giant abundance anomalies with a simulated hydrogen shell instability

## 1 Introduction

As long ago as 1947, Popper  had noted that the globular cluster red giant L199 in M13 was CN-strong. Further work by Harding , Osborn  and Hesser, Hartwick and McClure  confirmed that this anomaly is common. So for over fifty years it has been known that the surface abundances of elements such as C and N can vary from one red giant star to another within an individual globular cluster. It has become clear in the intervening years, however, that these anomalies are not restricted to C and N: heavier elements such as O, Na, Mg and Al also show star to star variations. Most other elements studied show little variation, and the degree of variation of C, N, O etc. differs from cluster to cluster and is much less apparent in field stars .

Observations by many groups, including Kraft et al. , Smith and Kraft , Langer , Shetrone \[1996a\] \[1996b\] , have identified this abundance anomaly problem and the most commonly suggested explanation for such high Na and Al are the deep mixing  and primordial enhancement  hypotheses. The deep mixing hypothesis suggests that the rotation of the star can give rise to meridional circulation currents that in turn create a mixing zone in the radiative region separating the top of the hydrogen burning shell (HBS) and the base of the convective envelope (BCE). This would allow elements processed in the HBS to be mixed to the surface as the star proceeds up the RGB. The primordial enhancement hypothesis suggests that the stars are born with these anomalies already present. The stars are formed from material ejected by intermediate mass ($`310M_{}`$) AGB stars from an earlier epoch. Without going into the full mechanisms of each process, either could in principle explain the observed anomalies, but modelling and observations have indicated that probably both processes occur to some extent.

Langer, Hoffman and Zaidins  suggested that the Na and Al abundance anomaly problem may be solved if the H shell reached higher temperatures for at least part of its lifetime whilst on the RGB. The idea stems from work by Von Rudloff and Vandenburg  and even earlier work by Bolton and Eggleton . Von Rudloff and Vandenburg  found that the H shell of a non-rotating small mass red giant model was stable but only marginally. A rotating model could be unstable. The link with rotation is particularly encouraging because this mechanism could explain both the deep mixing and the thermally unstable shell. The work performed by Langer et al  shows that a thermally unstable H shell would allow for the greater production of 27Al at the expense of 24Mg and that there would also be extra 23Na available for mixing to the surface.

The conclusion was that if the H shell reached temperatures of approximately 70-75 million K then the material mixed to the surface would produce a composition in close agreement with observation. In fact temperatures of this order and no lower would be required to match the Na observations because at the lower temperatures too much 23Na would be produced. Also it is only at these temperatures that significant 27Al is produced and 24Mg is only destroyed at temperatures above 70 million K. Therefore, to have any surface depletions in 24Mg the temperature of the H shell has to reach at least 70 million K unless primordial changes are invoked or the accepted reaction rates are altered. The preferred temperature range was 72-73 million K (see section 2 of Langer et. al ). This is the temperature at which they believed all the observed surface abundance anomalies could be matched.

As the structure of a model is heavily dependent upon the temperature profile any change to the temperature must be included in the evolution modelling. The following sections discuss models where temperatures in the H shell are artificially increased by decreasing the energy generation rate. The next section looks at oscillating temperature rises in the HBS; i.e. a pulsating model where each pulse immediately follows the previous pulse. In the following section pulses in the temperature of the HBS are created with a time gap between each pulse.

## 2 Computer Codes

The calculations reported here make use of two separate computer codes: the Monash/Mt Stromlo stellar evolution code, with which we model the evolution of stars from contraction to the main sequence up to the RGB tip, and a post-processing nucleosynthesis code , with which we model the abundances within the star for the same epoch. The nucleosynthesis code uses the independent variables from the evolution code for its structural basis, with typically every sixth evolution model being used. The full details of how these two codes work and interface can be found in Messenger .

The postulated deep mixing has been included in the nucleosynthesis code, and is described by three parameters - the depth, the speed and the time of onset of the deep mixing. We measure the depth in the same way as Wasserburg, Boothroyd and Sackmann , by mixing down to a temperature $`\mathrm{}\mathrm{log}T`$ above the base of the HBS, which is defined as the first point (moving outward from the center) to have a non-zero hydrogen abundance. This mixing commences once the HBS burns through the discontinuity left from the first dredge-up episode earlier in the evolution. Prior to this time it is assumed that this molecular weight discontinuity will inhibit mixing: once it is removed, the mixing is assumed to begin. Finally, we found (in agreement with previous researchers) that provided the speed of mixing was above a small critical value, it made almost no difference to the abundance patterns. We used a mixing speed of $`10^4M_{}`$/year in all calculations reported here.

In the work of Messenger  it was found that varying the mixing depth as the star ascends the RGB gives better agreement with observations of Na. Therefore this type of modelling is adopted in this paper, and the mixing depth was changed linearly from an initial setting $`\mathrm{}\mathrm{log}T=0.18`$ to $`\mathrm{}\mathrm{log}T=0.03`$ as the star ascended the RGB.

## 3 Increased H shell Temperatures

### 3.1 Algorithm

To increase the temperature of the H shell it is not simply a matter of arbitrarily adding temperature to the mesh points within the shell. The temperature is one of the dependent structure variables and affects other variables such as the density and pressure through the equation of state. To self-consistently increase the temperature we chose to change the energy generation of the shell by increasing or decreasing the energy production rate for each point if it is within the boundaries of the H shell.

An increase of the energy generation rate causes the star to expand and the temperature to decrease. A reduction in the energy generation rate causes the opposite; the star contracts and the temperature increases. Therefore it was necessary to decrease the energy rate from H burning in order to artificially create a temperature increase.

A factor for the energy generation rate reduction, designated as $`g`$, was introduced. This reduction is only initiated after the H shell passes through the $`\mu `$-discontinuity remaining from the first dredge-up episode. The desired nuclear energy generation rate at a point in the H shell (when this feature is switched on) is

$$ϵ_{used}=ϵ_{rate}(1g)$$
(1)
However, it is not simply a matter of reducing the energy generation rate over the full shell. This would create a large jump in energy generaqtion from the point just below the shell and also from the point at the top of the shell.111We define the top of the shell to be the point that first has a hydrogen abundance 0.98 of the surface abundance moving from the centre of the model outwards. It is necessary to linearly interpolate from 0 to the value of $`g`$ over a number of mesh points at the bottom of the shell and from $`g`$ to 0 at the top of the shell.

Let $`h_m`$ be the mesh-based linear interpolation factor. This varies from zero at the bottom of the HBS to a maximum of $`g`$ over $`n`$ mesh-points, remains at this maximum over $`2n`$ mesh points, and then decreases linearly over the next $`n`$ mesh points to a value of zero again. This is shown in Figure 1.

Once $`h_m`$ has been determined the energy generation rate remaining at any point within the H shell (as long as the H shell instability feature is switched on) is

$$ϵ_{used}=ϵ_{rate}(1h_m)$$
(2)
This alone caused a temporal discontinuity in the models and subsequently created convergence problems. Remember that this arbitrary temperature rise is non-physical. To avoid these convergence problems it was necessary to smoothly introduce the reduced energy generation rate in time as well as in mass (i.e. to gradually turn on the energy reduction). We introduce the reduction rate gradually over time with a parameter $`h_t`$ which is linearly increased from 0 to 1 over 20 timesteps. Thus the energy generation rate within the H shell is

$$ϵ_{used}=ϵ_{rate}(1h_mh_t)$$
(3)
Note that parametrisations in terms of mesh point (rather than mass) and model number (rather than time) are purely numerical conveniences. In reality, of course, these variations should be dependent on mass and time, but without any identified mechanism for the instability, and in view of the exploratory nature of the calculations reported here, we believe this simplification is justified and in no way affects our conclusions, which depend only on the shell temperature.

### 3.2 Results

All calculations reported here are carried out for a star of $`0.8M_{}`$ with a metallicity of $`Z=0.0005`$. Parameters relevant to the deep mixing were given above.

The temperature in the H shell is determined by the requirement that the energy production rate be high enough to support the outer layers of the star. With a decreased rate of energy generation at a given temperature, the star’s response is to increase the temperature in the H shell. As the star evolves up the RGB, the H shell progresses outwards and thins. The higher temperatures mimic this process and the final result is an acceleration of the star’s evolution. This is shown in figure 2, where plots of luminosity against time for differing values of $`g`$ are shown. The greater the energy rate reduction, the quicker the star evolves to the RGB tip.

This result will not assist in the abundance anomaly problem because any enhancements/depletions will still occur to about the same degree as before but just at an earlier epoch. The higher temperatures within the H shell mimic normal evolution but reach it at an earlier stage. All this manages to do is force the model to the He flash more quickly. No extra Na, Mg or Al would be produced in this situation. in fact, with a reduced timescale for the burning, the effect may be the opposite of what is required, but in any event it will be quantitatively negligible.

## 4 Intermittent Temperature Rises

Langer et. al.  suggested that the H shell is a thermally unstable region and that under certain conditions the shell may give rise to higher temperatures than standard physics suggests. The modelling performed in the previous section is not that of a thermally unstable region, but just a consistently higher temperature. To model a thermal instability the energy generation rate modification must not be steady but should be periodic or at least variable. A small modification to the algorithm of section 3 was therefore made to simulate a thermally unstable H shell.

The simplest way to do this is to continue to use the time based rate reduction algorithm but just make it periodic. Here the energy generation rate reduction was introduced over five timesteps, held at the peak for 20 timesteps and then decreased over another five timesteps. This was then repeated until the RGB tip was reached. For an initial test and to determine how pulsation may affect convergence there was no gap between each pulse. Apart from the increase and decrease of each pulse the affect would be identical to the continuous algorithm described in the previous section. We would therefore expect little change from the model for the previous algorithm and that is exactly the case. Figure 3 which has $`g=0.8`$ shows that the result of the above algorithm is to accelerate the evolution just as it did for the case of a continuous H shell temperature rise. The problem with the above approach is that the decrease in the energy generation rate causes the shell to contract and temperatures to rise. This shifts the star to a different evolutionary track. Once the energy generation rate returns to normal values the shell begins to expand again but there is not enough time for the star to settle back to its original (or close to original) evolutionary track before the next pulse begins. The end result is that the star shifts from its original evolutionary track and diverges from it henceforth.

If we do not want to significantly alter the evolution of the star, then it is necessary to introduce a time delay between pulses which is long enough for the star to resettle to its original track.

We use a simple method, but it is efficient and can test the concept of a H shell instability without the introduction of large amounts of complex code which is hard to justify since there is no confirmed mechanism for the instability. The purpose of our algorithm is to test the principle of a periodic H shell instability, not the details. We use a technique similar in principle to that shown in figure 1 but we introduce periodicity and time delays between pulses. We firstly select the desired number of pulses between FDU and the tip of the giant branch. We thus determine the inter-pulse duration and distribute the pulses evenly over the giant branch life-time. After some experimentation we decided to use 20 pulses because this allows the star to settle back to its standard evolutionary track and it should illustrate any differences likely to occur due to the pulses. We then take as a characteristic time-scale the time-step used to construct the current model, which we denote as $`dt`$. Next, the duration of a pulse is taken as $`30dt`$, comprising $`5dt`$ over which we linearly increase $`h_t`$ from zero to unity, then $`20dt`$ using $`h_t=1`$, and a further $`5dt`$ when $`h_t`$ is linearly decreased back to zero. Of course, the time-steps used by the code vary, and are not all equal to $`dt`$ during the pulse, but $`dt`$ gives a convenient timescale for the test. Also, because the time-steps vary with evolution along the giant branch, the $`dt`$ used for each of the 20 pulses will be slightly different, but this is not important for our exploratory calculations.

The result of the introduction of this algorithm is shown in figure 4. The points that should be noted are the decrease in the luminosity shown during a pulse accompanied by a temperature increase and the fact that the curve and hence the model settles back to its original track. It is only at the latter stages of the RGB evolution that the curves begin to diverge.

The conclusion drawn from this model is that the H shell temperature can be increased intermittently without significantly affecting the evolution of the star. But how will it affect the nucleosynthesis of the model? Is the temperature increase high enough and long enough to produce significant amounts of 27Al and reduce 24Mg as originally hoped? It is necessary to look at the temperatures within the H shell to see whether temperatures of approximately 70 million K are reached, the temperature at which 27Al enhancements are reasonable and 24Mg depletions begin. Figure 5 shows the temperature at the base of the H shell as the star ascends the RGB for an energy generation rate reduction of 50% ($`g=0.5`$). The peak temperature reached by any of the pulses is only increased very slightly, well short of the required temperatures for 24Mg destruction. Even a 100% reduction of the nuclear burning energy generation rate ($`g=1.0`$) during a pulse, as shown in figure 6, cannot produce high enough temperatures within the shell for the processes discussed here. For a test case, some of the gravitational and neutrino energy generation rates were also reduced ($`g=1.5`$). This created a minor increase in the temperature (see figure 7), but the peak temperature still fell short of the required temperatures.

Also if there are to be enhancements from this temperature increase, they must reach such temperatures early enough for significant composition differences to be seen in the envelope. The peak temperatures from these models are only seen during the very latest stages of RGB evolution, probably not early enough to match the data from M13 for example.

It seems that the required temperature increase cannot be produced via this method. The suggested instability has some desirable features but to self consistently produce the required temperatures proves difficult: our simulated instability does not succeed, and the question of whether a genuine instability would succeed remains open until such a mechanism can be found, and consistently modelled. We are therefore reluctant to dismiss this possibility, based on the technique attempted in this paper. Also a genuine instability may be different to what we have modelled here and may create higher temperatures in the HBS.

Despite this, temperature rises are seen whilst the pulse is active, and these could alter the nucleosynthesis. It would not be expected that any 24Mg would be destroyed at these temperatures, but 27Al production starts at lower temperatures and peak 23Na production at even lower temperatures. Even though the temperature pulses are not high enough to test the hypothesis espoused by Langer et al , there may still be some effect on the surface abundances of Na, O and Al as well as of the other lighter elements including such isotopes of Mg as 25Mg and 26Mg.

Table 1 lists all the models run where the energy generation rate was reduced. Case 6 was run through the nucleosynthesis code with deep mixing included. Figures 8, 9 and 10 show the surface abundances for Na, Mg and Al respectively with the curves for the standard model (Case 0) superimposed for comparison. It is quite clear that the temperature pulses have had a negligible effect on the surface abundances.

The temperature rises are either too small or too short lived to have any effect. Longer lived or more pulses only accelerate the evolution of the star and does not assist in maintaining higher temperatures and thus providing more 27Al for example. To test this, one of the models with constant energy generation rate reduction, Case 2 from table 1 was also run through the nucleosynthesis code, with its results and that of a standard model shown in figures 11, 12 and 13. Once again there is little effect on the surface abundances, as expected, although the Na values fit a little better. The process accelerates the evolution but does not significantly alter the surface abundances.

## 5 Conclusion

It appears that a thermally unstable H shell, modelled with the methods of this paper, cannot produce the high temperatures of about 70 million K required to produce the observed aluminium abundances. We succeeded only in accelerating the evolution with negligible effect on the interior temperatures and the surface abundances. Although a genuine instability, may manifest itself in a way which is quite different to these exploratory calculations, we have tried to artificially reproduce the most favourable conditions without success.

It appears that primordial abundances variations among the stars combined with deep mixing currently offer greater potential in making progress on this problem.