text:
# Atomic–Molecular Condensates with Large Positive Scattering Length

## Abstract

We show that in the limit of large and positive atom–atom scattering length the properties of an atomic–molecular Bose–Einstein Condensate (amBEC) are determined by an universal energy density functional (EDF). We find that the optimal conditions for the formation of a stable amBEC are in the regime where there are no shallow trimers and the atom–dimer scattering length is negative and comparable in magnitude with the atom–atom scattering length. At temperatures lower than $`T_c`$ the chemical potentials for the atoms and molecules can be specified independently. Besides three–body recombinations processes into dimers of large size, inelastic processes involving the formation of deeply bound small size molecular states are possible. These inelastic processes do not lead to an efficient heating of the amBEC and can be used for its mostly non–destructive monitoring.

The theoretical prediction , along with the experimental confirmation  of the fact that one can manipulate relatively easily the interaction of slow atoms by means of external fields, opened the doors to a long series of possible experiments. Perhaps among the most spectacular experimental results obtained lately are the creation of the so–called Bosenova  and the evidence that something related to a mixed atomic–molecular Bose–Einstein Condensate (BEC) was created , an object envisioned by theory several years before , see also Refs.  for a couple of theoretical analyses of this experiment. In the present work we identify a new universal regime, which appears in the case when the two–body scattering length $`a`$ is positive ($`a>0`$) and large, $`ar_0`$, where $`r_0`$ is the effective range, appearing the the low energy parameterization of the $`s`$–wave phase shift $`k\mathrm{cot}\delta (k)=1/a+r_0k^2/2+\mathrm{}`$ Typically $`r_0`$ is of the order of the so called van der Waals length $`(C_6m/\mathrm{}^2)^{1/4}`$ . The universality of this new regime is encoded in only two parameters, which define the effective Hamiltonian (except for the mass), one of them the two–body scattering length $`a`$ and the other a three–body characteristic, which can be chosen in various ways. As we will show below, in this regime, when $`na^31`$, the properties of a cold mixture of atoms and shallow dimers are described by a very simple theoretical model. In the other limit, when $`a<0`$ and $`|a|r_0`$, there is another universal regime and a new class of dilute quantum liquids exists . The basic ingredients of the model are derived from the well known Efimov effect , which pertains to the universal properties of the three–body systems with large two–body scattering lengths.

In order to derive the energy density of an atomic–molecular gas let us first consider the dilute regime where $`nr_0^31`$, with $`n`$ a typical atomic or molecular density and $`a=𝒪(r_0)`$. The details of the potential are not important and the atom–atom interactions can be described by a contact term

$$H_a=\psi _a^{}\frac{\mathrm{}^2\mathbf{}^2}{2m}\psi _a+\frac{1}{2}\lambda _2\psi _a^{\mathrm{\hspace{0.33em}2}}\psi _a^2+\frac{1}{3}\lambda _3\psi _a^{\mathrm{\hspace{0.33em}3}}\psi _a^3,$$
(1)
where $`\psi _a^{},\psi _a`$ are creation and annihilation operators for atoms. The value of $`\lambda _2`$ is determined (after fixing a short distance regulator) by the scattering length $`a`$. $`\lambda _3`$ depends also on a genuinely three–body length scale, denoted below by $`a_3^{}`$, that can be determined only through the value of a three–body observable. Terms with more derivatives or higher body forces can be included, but their effect on a dilute gas is small as long as the relevant momenta satisfy the condition $`kr_01`$ . The Hamiltonian in Eq. (1) contains, in principle, all the information necessary to describe the system in the dilute regime, including the molecular states/dimers with binding energy $`\mathrm{}^2/ma^2`$ and radius $`a`$ when $`a>0`$. In the case $`ar_0`$ we are interested in, the presence of bound states within the regime of validity of Eq. (1) evidences the fact that perturbation theory breaks down. Since particle–particle diagrams should be summed up to all orders, a simple mean–field approximation is not legitimate as well and two–particle correlations have to accounted for explicitly. Thus, even though $`H_a`$ contains all the information needed to describe molecular states, it is inconvenient to use it directly .

Let us consider now a regime even more dilute, where $`na^31`$. The relevant typical momenta are of the order $`p=\mathrm{}k\sqrt{na^3}\mathrm{}/a\mathrm{}/a\mathrm{}/r_0`$ . We will argue later that transitions altering the number of molecular states are “very slow” in this regime and, consequently, for time scales shorter than the transition rate, the number of atoms and molecules are separately conserved. An appropriate effective Hamiltonian describing the system in this regime is

$`H_{am}`$ $`=`$ $`\psi _a^{}{\displaystyle \frac{\mathrm{}^2\mathbf{}^2}{2m}}\psi _a\psi _m^{}({\displaystyle \frac{\mathrm{}^2\mathbf{}_a^2}{4m}}\epsilon _2)\psi _m`$ (2)
$`+`$ $`{\displaystyle \frac{1}{2}}\lambda _{aa}\psi _a^{\mathrm{\hspace{0.33em}2}}\psi _a^2+\lambda _{am}\psi _a^{}\psi _m^{}\psi _a\psi _m+{\displaystyle \frac{1}{2}}\lambda _{mm}\psi _m^{\mathrm{\hspace{0.33em}2}}\psi _m^2.`$ (3)
Both $`H_a`$ and $`H_{am}`$ are applicable in the regime $`kakr_01`$. $`H_{am}`$ has a great advantage over $`H_a`$ however. Perturbation or mean–field theory is not valid for $`H_a`$ when $`p\mathrm{}/a`$, but it is valid for $`H_{am}`$ since all non–perturbative physics occurring at the scale $`a`$ and leading to the formation of the bound state (dimers of size $`a`$) and to the Efimov effect (trimers of size $`a`$) is already neatly encapsulated in the constants $`\lambda _{aa},\lambda _{am}`$ and $`\lambda _{mm}`$. To make that more explicit, let us determine $`\lambda _{aa},\lambda _{am}`$ and $`\lambda _{mm}`$ in terms of $`\lambda _2`$ and $`\lambda _3`$ by considering zero momentum scattering of two atoms, atom–molecule and two molecules. Using $`H_a`$ these processes are described by an infinite number of diagrams shown in Fig. 1, corresponding to the solution of the Lippmann–Schwinger, Fadeev and Yakubovski equations respectively . The explicit solutions of the Lippmann–Schwinger and Fadeev equations are known. If we use dimensional regularization , where divergences proportional to powers of the momentum cutoff vanish, the same processes are described using $`H_{am}`$ simply by the Born amplitude, see the Fig. 1. This leads to the identification

$`\lambda _{aa}={\displaystyle \frac{4\pi \mathrm{}^2a}{m}},\epsilon _2={\displaystyle \frac{\mathrm{}^2}{ma^2}},`$ (4)
$`\lambda _{am}={\displaystyle \frac{3\pi \mathrm{}^2a_{am}}{m}}={\displaystyle \frac{3\pi \mathrm{}^2a}{m}}\left[c_1+c_2\mathrm{cot}\left(s_0\mathrm{ln}{\displaystyle \frac{a}{a_3^{}}}\right)\right],`$ (5)
$`\lambda _{mm}={\displaystyle \frac{2\pi \mathrm{}^2a_{mm}}{m}}={\displaystyle \frac{4\pi \mathrm{}^2a}{m}}c_3,`$ (6)
where the values of the universal constants $`s_01.00624`$, $`c_11.46`$ and $`c_22.15`$ are known  and $`|\epsilon _2|`$ is the true dimer binding energy. $`a_3^{}>0`$ is the value of the two–body scattering length for which a trimer energy is exactly at the atom–dimer threshold and this constant is system dependent, see Fig. 2. The specific value of $`a_3^{}>0`$ is determined by short–range three–particle physics, when all three atoms are in a volume of order $`r_0^3`$. The functional form of the atom–dimer scattering length was established by Efimov . It is natural to expect also that the dimer–dimer scattering length is of order $`a`$, as there are no known four–body shallow states near the $`42+2`$ threshold and $`aa_3^{}`$. The value of $`c_3`$ has not yet been calculated for bosons, but for fermions Randeria has shown at mean–field level that $`c_3=1`$ . This result, in connection with Gribakin and Flambaum’s reasoning , leads us to expect that in most cases the dimer–dimer scattering length is positive. Notice that the matching described above is performed in vacuum (length scales of $`𝒪(a)`$) and it is not modified by the many–body physics (length scales of $`𝒪(a/\sqrt{na^3})𝒪(a)`$). The fact that the dimensionally regulated $`H_{am}`$ is perturbative in the $`ka1`$ regime, unlike $`H_a`$, is firmly established .

The leading order contributions to the energy density functional are given by the mean field approximation of $`H_{am}`$, see Fig. 1 ( NB, a mean–field approximation for $`H_a`$ would be highly inaccurate even at these low densities!),

$$=\frac{1}{2}\lambda _{aa}n_a^2+\lambda _{am}n_an_m+\frac{1}{2}\lambda _{mm}n_m^2+\epsilon _2n_m,$$
(7)
where $`n_a,n_m`$ are the atom and dimer densities (see below the discussion on trimers contribution). In terms of many–body diagrams built using Eq. (1) this expression corresponds to the exact summation of the linear and quadratic terms in densities of all the two–particle, three–particle and four–particle diagrams, computed with the initial Hamiltonian $`H_a`$. Using well–known arguments , it is straightforward to show that higher order corrections, due to mean–field fluctuations, are controlled by the diluteness parameter $`\sqrt{na^3}1`$. This EDF expression is remarkable, since, in the Efimov regime when $`ar_0`$, it is universal. Only the exact value of the parameter $`c_3`$ for bosons is not explicitly known. Notice, however, that the condition $`na^31`$ precludes its use arbitrarily close to a Feshbach resonance. By changing the units of energy to $`\mathrm{}^2/ma^2`$ and length to $`a`$ this EDF acquires the following dimensionless and thus universal form:

$$=2\pi n_a^2+\frac{3\pi \xi }{2}n_an_m+2\pi c_3n_m^2n_m,$$
(8)
with $`\xi =c_1+c_2\mathrm{cot}[s_0\mathrm{ln}(a/a_3^{})]`$ tunable.

In the $`a>0`$ regime only the atom–dimer scattering length can change dramatically from negative to positive values and even become infinite. If $`\lambda _{aa}>0`$, $`\lambda _{mm}>0`$ and $`\lambda _{aa}\lambda _{mm}\lambda _{am}^2>0`$, then a mixed dilute atom–molecule system is a metastable gas. This can be seen as well by computing the pressure

$$P=\frac{1}{2}\lambda _{aa}n_a^2+\lambda _{am}n_an_m+\frac{1}{2}\lambda _{mm}n_m^2.$$
(9)
Phase separation occurs when $`\lambda _{am}>\sqrt{\lambda _{aa}\lambda _{mm}}`$ and domain wall structures are expected to appear . A (meta)stable trimer phase might also form if the (yet unknown) trimer–trimer interaction is repulsive. Another type of domain walls can appear in the opposite limit . It is easy to see that there is no region where $`P=0`$ and thus no stable liquid phase exists, unlike Refs. . An amBEC system would be in chemical equilibrium if

$$(2\lambda _{aa}\lambda _{am})n_a+(2\lambda _{am}\lambda _{mm})n_m=\epsilon _2.$$
(10)
This type of chemical equilibrium can be established through the reaction $`A+A+AA+A_2`$, which has a low rate for dilute systems $`\mathrm{}a^4n_a^2/m`$, see Refs.  and below. At temperatures $`T<T_c(na^3)^{2/3}\mathrm{}^2/ma^2|\epsilon _2|=\mathrm{}^2/ma^2`$, the rhs of the above chemical equilibrium condition is parametrically much larger than the lhs when $`na^31`$. The reverse reaction $`A_2+AA+A+A`$ is therefore exponentially suppressed, as it requires an activation energy larger then the thermal energy. Consequently, these type of collisions thus lead to a very small heating (and very slow chemical equilibration of the system) and its ultimate depletion . Since these two reactions (and $`A_2+A_2A_2+A+A`$, etc. as well) are very slow in an amBEC, the chemical potentials for atoms and dimers could be specified independently

$$\mu _a=\lambda _{aa}n_a+\lambda _{am}n_m,\mu _m=\lambda _{am}n_a+\lambda _{mm}n_m+\epsilon _2.$$

In a trap one has to add the trap potentials for atoms and dimers in the rhs of these two relations .

A number of other processes can occur in such a system. It is desirable to consider three–body inelastic collisions involving one or more dimers in the initial state. Particularly interesting are however the inelastic atom–dimer $`A+A_2A^{}+A_2^{}`$ (see Ref.  for a rate calculation) and dimer–dimer $`A_2+A_2A^{}+A^{}+A_2^{}`$ collisions, where $`A^{}`$ and $`A_2^{}`$ stand for a fast atom (with its kinetic energy significantly exceeding the thermal energy $`3T/2`$) and $`A_2^{}`$ for a fast dimer respectively (with a binding energy of the order $`\mathrm{}^2/mr_0^2\mathrm{}^2/ma^2`$). If the temperature of the system is significantly lower then the dimer binding energy, i.e. $`T<T_c\mathrm{}^2/ma^2`$, these inelastic collisions are not leading to a significant heating of the system, but mostly to its depletion, in a manner similar to the depletion in the BEC of metastable helium . If an atom and a dimer collide and this leads to the formation of a dimer with a binding energy of order $`\mathrm{}^2/mr_0^2\mathrm{}^2/ma^2`$ a relative large amount of energy is released (NB $`T<T_c\mathrm{}^2/ma^2\mathrm{}^2/mr_0^2`$). It is conceivable that transitions to various final bound states could be resolved and thus the spectrum of fast ejectiles would be discrete. The momenta of the outgoing atom and of the deeply bound dimer are large as well, of the order $`\mathrm{}/r_0\mathrm{}/a`$, and such fast atoms and dimers will interact relatively weakly with the rest of the atoms and dimers. While the slow/thermal atoms and shallow dimers interact with each other with cross sections of the order of $`a^2`$, a fast atom and a deeply bound and fast dimer will interact with the amBEC with a cross section of order $`r_0^2a^2`$. One can expect a similar behavior in a dimer–dimer inelastic reaction. The fact that the products of these inelastic processes are fast and weakly interacting objects can be used for a non–destructive monitoring of the amBEC. The rates of these inelastic processes are obviously controlled by the parameter $`r_0^3/a^31`$, since in order to form a deeply bound dimer two atoms in the shallow dimer have to change their separation from distances of order $`a`$ to distances of order $`r_0`$. One cannot help here but to make a parallel with the similar role played by neutrinos (which are also fast and weakly interacting products of inelastic processes) in studying the interior of the Sun.

A possible final channel of a dimer–dimer collision could be a trimer–atom final state if $`a_{am}>0`$. If the trimer energy is close to the atom–dimer threshold, then the atom–dimer cross section could be significantly larger than the atom–atom and dimer–dimer cross sections ($`|a_{am}|a,a_{mm}`$). If this is the case, one can expect that the atom–dimer elastic collisions would be chiefly responsible for the thermalization of the system. However, besides the elastic channel there is also an inelastic channel present if $`a_{am}>0`$: a dimer–dimer collision leading to a shallow trimer and an atom in the final state. This process will lead to a relatively small release of energy since the trimer energy is very close to the threshold, see Fig. 2. With good accuracy the ground state energy of a trimer is given in this limit by $`\epsilon _3\epsilon _23\mathrm{}^2/4ma_{am}^2.`$  Moreover, the trimer under these conditions is basically a very large two–body object of size $`a_{am}a`$, made of a dimer and a relatively loosely bound atom to it. In this limit the atoms and dimers will phase separate, since $`\lambda _{am}\sqrt{\lambda _{aa}\lambda _{mm}}`$. A pure trimer phase could appear instead, iff the trimer–trimer interaction is repulsive. When $`a_{am}>0`$ and $`\lambda _{am}\sqrt{\lambda _{aa}\lambda _{mm}}`$ dimer–dimer collisions with a trimer in the final state will lead to an energy release of the order of $`\mathrm{}^2/ma^2`$, and thus to the heating of the system, since the momenta of the emerging products are of the order $`\mathrm{}/a`$ in this case. Consequently, the best conditions to create a stable amBEC are expected when $`a_{am}<0`$ (no shallow bound trimers and thus no $`A_2+A_2A_3+A`$ reactions) and $`|a_{am}|<a`$. As a last note, there is essentially no overlap between present results and those of Refs. .

Discussions with David B. Kaplan and Lan Yin are greatly appreciated.


format: text

text:
# Edge critical currents of dense Josephson vortex lattice in layered superconductors

## I Introduction

Atomically layered superconductors, such as Bi2Sr2CaCu2Ox (BSCCO), behave as stacks of Josephson junctions. This so-called intrinsic Josephson effect has been a subject of intense research in the past decade (see reviews IJJReviews ; YurgensReview ). In spite of atomic-size separation between the neighboring junctions, in zero magnetic field the system approximately behaves like an array of independent junctions. A weak dynamic interaction between the neighboring junctions appears due to nonequilibrium effects nonequil .

Situation is very different in the magnetic field applied along the layer direction, which generates the Josephson vortex lattice. Recently, dynamic properties of this lattice in BSCCO have been subject of extensive experimental research Lee ; Hechtfischer ; Latyshev ; Ooi01 . At high fields (above 0.5 tesla for BSCCO) the Josephson vortices homogeneously fill all layers BulClemPRB91 ; KorshLarkPRB92  (dense Josephson vortex lattice, see Fig. 1). Josephson vortex arrays in neighboring layers have strong inductive interaction Inductive , which strongly influences static and dynamic properties of the lattice. Due to this interaction the magnetic properties of the stack are very different from the magnetic properties of an isolated junction.

A transport current flowing across the layers exerts a Lorentz force on the Josephson vortices along the layer direction. An important parameter is the *critical current* above which the lattice starts to move producing a finite voltage. This current is determined either by bulk pinning or by interaction with the boundaries. In this paper we consider the case of a homogeneous junctions and neglect bulk pinning. The simplest and most known case is a single small Josephson junction without inhomogeneities, where the field dependence of the critical current is given by the Fraunhofer dependence, $`I_c(\mathrm{\Phi })=I_{c0}|\mathrm{sin}(\pi \mathrm{\Phi }/\mathrm{\Phi }_0)|/(\pi \mathrm{\Phi }/\mathrm{\Phi }_0)`$, with $`\mathrm{\Phi }`$ being the magnetic flux through the junction. Observation of this dependence has been considered as an important confirmation of the dc Josephson effect. At present, the Fraunhofer field dependence is routinely used as an indicator of the junction homogeneity. The same dependence is also expected for the junction stack with the lateral size smaller than the Josephson length BulPRB92 . In a single long junction the field dependence of the critical current has rather complicated structure due to multiple coexistent states of the lattice LongJosJunct .

In this paper we consider the critical current and lattice structure near boundaries at high magnetic fields in the regime of dense lattice in the stack of Josephson junctions. In this regime the Josephson coupling can be treated as a small perturbation, which allows to develop a full analytical description of this system. Another important simplification is that at high fields one can neglect the current self-field and neglect the difference between the magnetic field inside the junction and the external applied magnetic field.

We find that the boundary induces an alternating deformation of the lattice. Averaging out the rapid phase oscillations, we derive that the lattice deformation obeys the Sine-Gordon equation. The deformation decays inside superconductor at the typical length $`l_B`$, which is larger than the Josephson length and increases proportional to the magnetic field, $`l_BB`$. Structure of the surface deformation and the total current flowing along the surface is uniquely determined by the lattice position far away from the boundary. The total current flowing through the stack is given by the sum of two independent surface currents flowing at the edges of the sample. Due to the commensurability effects, the maximum current through the stack has oscillating field dependence, which resembles the Fraunhofer dependence: it has strong oscillations and overall $`1/B`$ dependence. However, the period of these oscillations corresponds to adding *one flux quantum per two junctions*.

We also calculate the flux-flow voltage for the Josephson vortex lattice slowly moving through a finite stack. Due to the interaction with the boundaries this voltage also has an oscillating contribution with period one flux quantum per two junctions. Recently, such oscillations have been observed in BSCCO mesas Ooi01 , and also were seen in numerical simulations Machida02 . We find that the oscillating voltage normalized to the bare flux-flow voltage has universal dependence on the magnetic flux per one junction and the current density normalized to the Josephson current density. At high currents the relative amplitude of voltage oscillations decreases with current and field as $`1/J^2H^2`$. We also show that in the case of dominating in-plane dissipation, typical for BSCCO, the *absolute* amplitude of voltage oscillations weakly depends on magnetic field in a wide field range, in agreement with experiment Ooi01 .

The new length scale $`l_B`$ increasing with the magnetic field also sets the new field scale $`B_L`$, at which this length becomes of the order of the junction length $`L`$. Above this field the system crosses over to the regime of small junction and the field dependence of the critical current crosses over to the usual Fraunhofer dependence. Regimes of different behavior for a finite stack in magnetic field are summarized in the phase diagram shown in the right part of Fig. 1.

The paper is organized as follows. In Section II we write general relations describing the phase distribution for layered superconductor in the high magnetic field applied along the layers. In Section III we consider lattice structure near the surface, calculate the surface energy and surface current. In Section IV we calculate the oscillating field dependence of the critical current in a finite stack. In Section V we consider the field oscillations of the flux-flow voltage due to the interaction with the boundaries for slow lattice motion.

## II  General relations

Consider a stack of the Josephson junctions with length $`L\lambda _J\gamma s`$, where $`\gamma `$ is the anisotropy parameter, in high magnetic field $`BB_{cr}=\mathrm{\Phi }_0/2\pi \gamma s^2`$ applied along the direction of the layers (see Fig. 1). At high field the screening effects can be neglected and the phase distribution $`\phi _n(y)`$ is described by the energy (per junction)

$$_J=\frac{w}{N}\underset{n}{}_0^L𝑑y\left[\frac{J}{2}\left(\frac{d\phi _n}{dy}\right)^2E_J\mathrm{cos}\left(\phi _{n+1}\phi _n\frac{2\pi sB}{\mathrm{\Phi }_0}y\right)\right].$$
(1)
where $`J`$ is the in-plane phase stiffness, $`E_J`$ is the Josephson energy per unit area, and $`w`$ is the stack size in the field direction. To facilitate analysis, we introduce dimensionless parameters

$$u=y/\lambda _J,\stackrel{~}{L}=L/\lambda _J,e_J=\frac{_J}{E_J\lambda _Jw},h=\frac{2\pi \gamma s^2B}{\mathrm{\Phi }_0}.$$

The reduced energy $`e_J`$ is given by

$$e_J=\frac{1}{N}\underset{n}{}\underset{0}{\overset{\stackrel{~}{L}}{}}𝑑u[\frac{1}{2}\left(\frac{d\phi _n}{du}\right)^2\mathrm{cos}\left(\phi _{n+1}\phi _nhu\right)].$$
(2)
In a stable state the phase obeys the following equation, which expresses the current conservation,

$$\frac{d^2\phi _n}{du^2}+\mathrm{sin}\left(\phi _{n+1}\phi _nhu\right)\mathrm{sin}\left(\phi _n\phi _{n1}hu\right)=0$$

with the boundary conditions (vanishing of the in-plane current at the boundaries)

$$\frac{d\phi _n}{du}|_{u=0,L}=0.$$
(3)
In principle, a stable state can support a finite Josephson current flowing through the stack

$$j=j_J\lambda _Jw_0^{\stackrel{~}{L}}𝑑u\mathrm{sin}\left(\phi _{n+1}\phi _nhu\right).$$
(4)
In a large-size sample this current is concentrated near the edges. In the following section we will derive a simple closed expression for the edge current.

## III Lattice structure near boundaries. Surface energy and surface current.

Far from the boundaries (exact criterion will be established below) the ground state configuration corresponds to the rhombic lattice (Fig. 1) and the phase distribution can be easily calculated using expansion with respect to the Josephson current BulClemPRB91 ; KorshLarkPRB92

$$\phi _nn\alpha +\frac{n(n+1)}{2}\pi \frac{2}{h^2}\mathrm{sin}\left(hu\alpha \pi n\right)$$
(5)
The bulk part of the energy is degenerate with respect to the phase shift $`\alpha `$. Change of $`\alpha `$ corresponds to translational displacement of the lattice. This degeneracy is eliminated by the interaction with the boundaries. To study this interaction, we introduce the new phase variable $`\varphi _n(u)`$

$$\phi _n=n\alpha +\pi \frac{n(n+1)}{2}+\varphi _n$$

and impose the condition that $`\varphi _n`$ contains only oscillating contribution far away from the boundaries, $`\varphi _n(u)(1)^n(2/h^2)\mathrm{sin}\left(hu\alpha \right)`$. This new phase $`\varphi _n(u)`$ obeys the following equation

$$\frac{d^2\varphi _n}{du^2}+\mathrm{sin}\left(\varphi _{n+1}\varphi _nhu+\alpha +\pi n\right)+\mathrm{sin}\left(\varphi _n\varphi _{n1}hu+\alpha +\pi n\right)=0.$$
(6)
We assume that the boundary does not change the alternating nature of the vortex lattice and will look for solution of this equation in the form

$$\varphi _n(u)=(1)^n\varphi (u).$$

For the phase $`\varphi (u)`$ we obtain the following equation KorshLarkPRB92

$$\frac{d^2\varphi }{du^2}+2\mathrm{cos}2\varphi \mathrm{sin}\left(hu+\alpha \right)=0$$
(7)
with the boundary condition $`d\varphi /du|_{u=0,L}=0`$. The energy and total Josephson current in terms of $`\varphi (u)`$ are given by

$`e_J(\alpha )`$ $`={\displaystyle _0^{\stackrel{~}{L}}}𝑑u[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{d\varphi }{du}}\right)^2\mathrm{sin}\left(2\varphi \right)\mathrm{sin}(hu+\alpha )],`$ (8)
$`j(\alpha )`$ $`=j_J\lambda _J{\displaystyle _0^{\stackrel{~}{L}}}𝑑u\mathrm{sin}(2\varphi )\mathrm{cos}\left(hu+\alpha \right)`$
$`=j_J\lambda _J{\displaystyle \frac{e_J}{\alpha }}`$ (9)
Note that the energy and current have the symmetry $`\alpha \alpha +\pi `$, $`\varphi (u)\varphi (u)`$ (change of $`\alpha `$ by $`\pi `$ is equivalent to vertical displacement of the lattice by one junction). Therefore, the energy is $`\pi `$-periodic function of $`\alpha `$, $`e_J(\alpha +\pi )=e_J(\alpha )`$.

We now focus on the lattice structure near the boundary $`u=0`$. In the limit of high field $`h1`$ one can separate and average out the rapidly oscillating part of the phase. This technique has been used in Ref. KorshLarkPRB92,  to study melting of the Josephson vortex lattice. We split $`\varphi `$ into the smooth ($`v`$) and rapidly changing ($`\stackrel{~}{\varphi }`$) components,

$$\varphi =v+\stackrel{~}{\varphi },$$

where the smooth component describes deformation of the lattice induced by the boundary and satisfies conditions $`dv/duv`$ and $`v0`$, at large $`u`$. The local lattice compression is given by $`(v/u)/h`$. The oscillating part by definition obeys the following equation

$$\frac{d^2\stackrel{~}{\varphi }}{du^2}+2\mathrm{cos}\left(2v\right)\mathrm{sin}\left(hu+\alpha \right)=0.$$
(10)
Neglecting a weak coordinate dependence of $`v`$ in comparison with the rapidly oscillating sine function, we obtain an approximate solution of the last equation

$$\stackrel{~}{\varphi }\frac{2}{h^2}\mathrm{cos}\left(2v\right)\mathrm{sin}\left(hu+\alpha \right)$$
(11)
In the limit $`h1`$ the oscillating phase $`\stackrel{~}{\varphi }`$ can be treated as a small perturbation. Subtracting Eq. (10) from Eq. (7) we obtain equation for $`v`$ in the first order with respect to $`\stackrel{~}{\varphi }`$

$$\frac{d^2v}{du^2}4\stackrel{~}{\varphi }\mathrm{sin}\left(2v\right)\mathrm{sin}\left(hu+\alpha \right)=0$$

Substituting expression for $`\stackrel{~}{\varphi }`$ from Eq. (11) and averaging out rapid oscillations, we finally obtain Sine-Gordon equation for the smooth lattice deformation $`v(u)`$,

$$\frac{h^2}{2}\frac{d^2v}{du^2}\mathrm{sin}\left(4v\right)=0.$$
(12)
Using Eq. (11), we also obtain the boundary condition for $`v(u)`$,

$$\frac{dv}{du}|_{u=0}=\frac{d\stackrel{~}{\varphi }}{du}|_{u=0}=\frac{2\mathrm{cos}2v_0\mathrm{cos}\alpha }{h}$$
(13)
with $`v_0v(0)`$. As follows from Eq. (12) the boundary deformation decays at the typical length $`l_h=h/2\sqrt{2}`$ or, in the real units,

$$l_B=\gamma s\frac{\pi \gamma s^2B}{\sqrt{2}\mathrm{\Phi }_0}.$$
(14)
This important length scale appears as an interplay between the shear and compression stiffnesses of the lattice. A finite system is in the large-size limit if $`L>2l_B`$. This condition always brakes at sufficiently large field, $`B>B_L`$,

$$B_L=\frac{L}{\gamma s}\frac{\mathrm{\Phi }_0}{2\pi \gamma s^2}.$$
(15)
Averaging out the rapid oscillations, we also derive the surface energy $`e_s(\alpha )`$ and the surface current $`j_s(\alpha )`$ (in units of $`j_J\lambda _Jw`$) in terms of $`v(u)`$

$`e_s(\alpha )`$ $`{\displaystyle \frac{1}{h}}\mathrm{sin}\left(2v_0\right)\mathrm{cos}\left(\alpha \right)`$
$`+{\displaystyle _0^{\mathrm{}}}𝑑u\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{dv}{du}}\right)^2+{\displaystyle \frac{1\mathrm{cos}4v}{2h^2}}\right]`$
$`j_s(\alpha )`$ $`={\displaystyle \frac{1}{h}}\mathrm{sin}\left(2v_0\right)\mathrm{sin}\left(\alpha \right)`$
Equation (12) has the well-known soliton solution

$$\mathrm{tan}v=\mathrm{tan}v_0\mathrm{exp}\left(2\sqrt{2}u/h\right)$$
(16)
Employing the boundary condition (13), we derive relation between the boundary deformation $`v_0`$ and the bulk phase shift $`\alpha `$

$$\mathrm{tan}\left(2v_0\right)=\sqrt{2}\mathrm{cos}\alpha $$
(17)
From this equation we can see that the maximum deformation $`v_0=\mathrm{arctan}(\sqrt{2})/20.478`$ occurs at $`\alpha =0,\pi `$ and the deformation vanishes at $`\alpha =\pi /2`$.

Equations (16) and (17) determine the surface deformation $`v(u)`$ for the given bulk phase shift $`\alpha `$ (i.e., lattice position). From the surface deformation we can restore all other surface properties of the lattice. Fig. 2 shows distribution of the Josephson current near the boundary in the two neighboring junctions for two magnetic fields, $`h=3`$ and $`5`$, and different values of the $`\alpha `$. To illustrates the surface deformation, we also show the vortex center positions extrapolated from the bulk region assuming undeformed lattice.

From the obtained solution for $`v(u)`$ we derive the closed analytical expressions for the surface energy and the surface Josephson current

$`e_s(\alpha )`$ $`={\displaystyle \frac{1}{\sqrt{2}h}}\left(1\sqrt{2+\mathrm{cos}2\alpha }\right)`$ (18)
$`j_s(\alpha )`$ $`={\displaystyle \frac{1}{\sqrt{2}h}}{\displaystyle \frac{\mathrm{sin}2\alpha }{\sqrt{2+\mathrm{cos}2\alpha }}}`$ (19)
The $`\alpha `$-dependence of the surface current is plotted in the left panel of Fig. 3. The maximum surface current, $`j_{s,\mathrm{max}}`$, is achieved at $`\mathrm{cos}2\alpha =2+\sqrt{3}`$ and in real units is given by

$$j_{s,\mathrm{max}}=\sqrt{2\sqrt{3}}j_J\frac{\mathrm{\Phi }_0}{2\pi sB}.$$
(20)
Two possible directions of this current correspond to the Lorentz force on Josephson vortices directed from the boundary and towards the boundary, promoting either entrance or exit of vortices. In contrast to the regime of dilute lattice at $`B<B_{cr}`$, in which the maximum values of the entrance and exit current are very different (see, e.g., Ref. BurPRB96, ), in the dense lattice they are approximately the same. It is interesting to note that, up to numerical factor, the expression for the entrance current in the dilute lattice regime coincides with Eq. (20).

## IV Critical current of a finite stack

Consider now the total current flowing through a finite stack. For the given phase shift $`\alpha `$ (lattice position) the total current, $`J_L(\alpha )J_L(\alpha ,hL,h)`$, is the sum of the two surface contributions coming from the boundaries $`u=0`$ and $`u=L`$,

$$J_L(\alpha )=\frac{1}{\sqrt{2}h}\left(\frac{\mathrm{sin}2\alpha }{\sqrt{2+\mathrm{cos}2\alpha }}\frac{\mathrm{sin}2(hL\alpha )}{\sqrt{2+\mathrm{cos}2(hL\alpha )}}\right)$$
(21)
Introducing notation $`\chi =hL=2\pi \mathrm{\Phi }/\mathrm{\Phi }_0`$, where $`\mathrm{\Phi }=BsL`$ is the magnetic flux through one junction, we obtain that the maximum current $`J_c(\chi ,h)=f(\chi )/h`$ is given by

$$J_c(\chi ,h)=\underset{\alpha }{\mathrm{max}}\left[J_L(\alpha ,\chi ,h)\right]$$
(22)
Numerically obtained dependence of $`f=hJ_c`$ vs $`\chi `$ is plotted in the right panel of Fig. 3. In the real units, the field dependence of the maximum Josephson current in the field range $`\mathrm{\Phi }_0/2\pi \gamma s^2<B<L\mathrm{\Phi }_0/\pi \gamma ^2s^3`$ is given by

$$I_c(B)=I_J\frac{\mathrm{\Phi }_0}{2\pi sLB}f\left(\frac{2\pi sLB}{\mathrm{\Phi }_0}\right),$$
(23)
where $`I_J=j_JLw`$ is the maximum Josephson current through the stack at zero field, and oscillating function $`f(\chi )`$ is plotted in the right panel of Fig. 3. To facilitate comparison with experiment we also obtain an interpolation formula for $`f(\chi )`$ for $`0<\chi <\pi /2`$, $`f(\chi )0.128+0.888\mathrm{cos}(\chi )+0.021\mathrm{cos}(3\chi )`$. The derived field dependence of the critical current somewhat resembles the well-known Fraunhofer dependence in a single small-size Josephson junction: it has overall $`1/B`$ dependence with large oscillations. However, it has a very different physical origin and also has several qualitative differences. The most important difference is that the period of oscillations corresponds to adding *one flux quantum per two junctions*, i.e., it is two times smaller than for the Fraunhofer oscillations. This kind of oscillations have been recently observed in the flux-flow voltage of BSCCO mesas Ooi01 . Also, contrary to the Fraunhofer dependence, the points where the magnetic flux inside the junction equals to integer flux quanta, correspond to the local maxima of the critical current. Dependence of the critical current on the magnetic flux per junction is plotted in Fig. 4. For comparison, the Fraunhofer dependence is also shown. The dependence (23) holds until two edges give independent contributions to the total current. This condition breaks when the magnetic field exceeds $`B_L=L\mathrm{\Phi }_0/(2\pi \gamma ^2s^3)`$. At higher fields the field dependence crosses over to the usual Fraunhofer dependence. Recently, this crossover has been studied by numerical simulations Machida02 .

## V Oscillating flux-flow voltage

When the external current exceeds the critical current, the lattice starts to move. At slow motion the surface deformation has time to adjust to the current lattice position. In this case the surface energy produces a periodic potential for the moving lattice and one can use static results to predict the I-V dependencies. Time variation of the lattice phase shift obeys equation

$$\nu _{ff}\frac{d\alpha }{dt}+J_L(\alpha )=J,$$
(24)
where the total surface current $`J_L(\alpha )J_L(\alpha ,hL,h)`$ is given by Eq. (21) (for brevity we again skip in equations its dependence on the magnetic field and size) and the viscosity coefficient, $`\nu _{ff}`$, is related to the flux-flow resistance of the stack, $`R_{ff}`$,

$$\nu _{ff}=\frac{N\mathrm{\Phi }_0}{2\pi cR_{ff}}.$$

where $`N`$ is the number of junctions in the stack. Voltage drop per one junction $`U`$ is related to $`d\alpha /dt`$ by the Josephson relation

$$U=\frac{\mathrm{\Phi }_0}{2\pi c}\frac{d\alpha }{dt}.$$

Solution of Eq. (24) is given by the implicit relation

$$_0^\alpha \frac{\nu _{ff}d\alpha ^{}}{JJ_L(\alpha ^{})}=t,$$

from which we obtain the average phase change rate

$$\overline{\frac{d\alpha }{dt}}=\left[\frac{1}{\pi }_0^\pi \frac{\nu _{ff}d\alpha }{JJ_L(\alpha )}\right]^1$$

and the flux-flow voltage

$$\frac{U}{U_{ff}}=\left[\frac{J}{\pi }_0^\pi \frac{d\alpha }{JJ_L(\alpha )}\right]^1$$
(25)
with $`U_{ff}=R_{ff}J`$ being the bare flux-flow voltage without periodic potential. Because the surface current $`J_L(\alpha )J_L(\alpha ,hL,h)`$ oscillates with the magnetic field, this flux-flow voltage will also experience similar field oscillations. Such oscillations have been recently observed by Ooi *et. al.* Ooi01  Using scaling property of the current $`J_L(\alpha ,hL,h)=LF(\alpha ,hL)`$ we can see from Eq. (25) that the reduced flux-flow voltage $`U/U_{ff}`$ depends only on two parameters: the current density normalized to the Josephson current density, $`iJ/Lj/j_J`$, and the magnetic flux through one junction, $`\mathrm{\Phi }=BsL`$. This allows for a unified description of the voltage oscillations in junction stacks with different sizes. Oscillating dependencies of $`U/U_{ff}`$ on the magnetic flux $`\mathrm{\Phi }`$ for different current densities are shown in Fig. 5.

At large currents one can expand this result with respect to surface current and obtain a small correction to the flux-flow voltage

$$\frac{U}{U_{ff}}1\frac{g(0)+g(\chi )}{\pi J^2h^2}$$
(26)
where the oscillating part is described by dimensionless function

$`g(\chi )`$ $`={\displaystyle _0^\pi }{\displaystyle \frac{\mathrm{sin}2\alpha \mathrm{sin}2(\alpha \chi )}{\sqrt{2+\mathrm{cos}2\alpha }\sqrt{2+\mathrm{cos}2(\chi \alpha )}}}{\displaystyle \frac{d\alpha }{\pi }},`$
$`0.263\mathrm{cos}2\chi +0.0046\mathrm{cos}4\chi .`$
Therefore, the absolute amplitude of voltage oscillations is given by

$$\delta U=0.527\frac{R_{ff}}{J}\left(j_J\lambda _J\right)^2\left(\frac{\mathrm{\Phi }_0}{2\pi s\lambda _JB}\right)^2.$$
(27)
Independently from the dissipation mechanism, the relative amplitude of oscillations $`\delta U/R_{ff}J`$ decreases with field as $`1/B^2`$. In the regime of dominating in-plane dissipation, typical for BSCCO, the flux-flow resistance is given by JosFluxFlow

$$R_{ff}=R_c\frac{B^2}{B^2+B_\sigma ^2};B_\sigma =\sqrt{\frac{\sigma _{ab}}{\sigma _c}}\frac{\mathrm{\Phi }_0}{\sqrt{2}\pi \gamma ^2s^2}$$

where $`\sigma _{ab}`$ and $`\sigma _c`$ are the components of quasiparticle conductivity. Using this relation one can rewrite Eq. (27) as

$$\delta U=\frac{R_{ab}J_M^2}{2J}\frac{1}{1+B^2/B_\sigma ^2}$$
(28)
with $`J_M=c\mathrm{\Phi }_0/8\pi ^2\lambda _{ab}^2`$, $`R_{ab}=Ns/S\sigma _{ab}`$, and $`S=wL`$ is the junction area. For BSCCO, typically, $`B_\sigma 10`$T. In the region $`B<B_\sigma `$ the oscillation amplitude weakly depends on the magnetic field, in agreement with experiment (see Fig. 2 in Ref. Ooi01, ). It interesting to note that in this regime it is mainly determined by the in-plane parameters of superconductor.

## VI Summary and acknowledgements

In conclusions, we found that the edge current of the dense Josephson lattice is uniquely determined by the magnetic field and position of the lattice in the bulk. Near the surface the lattice has alternating deformation, which decays inside superconductor at the typical length which is proportional to the magnetic field. Due to the rhombic lattice structure, both the critical current and the flux-flow voltage at small velocities have oscillating field dependencies with the period of one flux quantum per two junctions up to the size-dependent magnetic field.

The author thanks K. Hirata and M. Machida for useful discussions. This work was supported by the U. S. DOE, Office of Science, under contract # W-31-109-ENG-38.


text:
# Correction to scaling analysis of diffusion-limited aggregation

## 1 Introduction

Since its introduction by Witten and Sander in 1981 , diffusion-limited aggregation (DLA) has been the fundamental stochastic model of quasistatic growth processes where the growth is limited by a diffusion process. The model can be described in simple terms: a rigid aggregate grows by the capture of a low density of Brownian particles, which attach to it on first contact. A highly ramified branching structure is produced, which—at least on first sight—appears to be fractal.

One of the most basic questions asked about DLA is whether the growing clusters obey simple scaling, i.e. are they indeed simple fractals? Based on numerical simulations, it has been suggested that the scaling is more complex: multiple divergent length scales might be present , the ensemble variance of cluster radii might have anomalous scaling , there could be more than one fractal dimension , or the clusters might obey multiscaling where the fractal dimension continuously depends on the position .

We claim that these anomalous scaling claims are *wrong*, they are misled by finite size transients. In particular, we will show in detail that subdominant terms in the scaling account for the apparent “multiscaling” observed in small to medium size simulations. For clarity we should mention that here *multiscaling* refers to the space dependent fractal dimension (anomalous scaling). This should not be confused with the well established *multifractality* of the harmonic (growth) measure, which is consistent with the simple asymptotic scaling of the clusters.

## 2 Simple scaling

In this section we will look at the scaling of various characteristic lengths: the *deposition radius* $`R_{\mathrm{dep}}=r`$ (the average distance of newly arriving particles from the center), the cluster’s *gyration radius* $`R_{\mathrm{gyr}}=\sqrt{\frac{1}{N}_{N^{}=1}^Nr^2_N^{}}`$, the *root-mean-square radius* $`R_2=\sqrt{r^2}`$, and the *penetration depth* $`\xi `$ (the width of the active zone ring, where newly deposited particles land). According to our results, DLA obeys simple scaling; all length scales scale with the same fractal dimension. To illustrate this, we look at one of the anomalous scaling claims mentioned in the introduction: that the cluster radius $`R_{\mathrm{dep}}`$ does not scale in the same way as the penetration depth $`\xi `$  — although it is worth mentioning that this claim has been questioned very soon . The ratio of the two, often called *relative penetration depth*, $`\mathrm{\Xi }=\xi /R_{\mathrm{dep}}`$, in our measurements obeys the asymptotic form  for large $`N`$:

$$\mathrm{\Xi }(N)\mathrm{\Xi }_{\mathrm{}}(1+CN^\nu ).$$
(1)
On Figure 1 we plot $`\mathrm{\Xi }`$ against $`N^\nu `$ with an appropriately chosen $`\nu `$; the linear behavior at $`N^\nu 0`$ clearly indicates the validity of Eq. (1). It is not easy to obtain numerically the exponent $`\nu `$: systematic errors (fitting data far from the asymptotic point) have to be balanced with large statistical errors (fitting close to the asymptotic point). Nevertheless, all data presented in this paper is consistent with a single “universal” exponent $`\nu =0.33\pm 0.06`$.

When the ratios of various lengths defined on DLA obey Eq. (1), and the lengths have an asymptotic power-law dependence on $`N`$, then they can be written in a scaling form with a leading subdominant term :

$$R(N)\widehat{R}N^{1/D}(1+\stackrel{~}{R}N^\nu )$$
(2)
A numerical proof of this finite size scaling is plotted on Figure 2. For completeness, we collected the corresponding coefficients $`\widehat{R}`$ and $`\stackrel{~}{R}`$ for many characteristic lengths in Table 1.

The coefficients are not independent, it is easy to derive some relations between them by neglecting higher order corrections: $`\widehat{R}_2=\sqrt{\widehat{R}_{\mathrm{dep}}^2+\widehat{\xi }_0^2}`$ and $`\stackrel{~}{R}_2=(\stackrel{~}{R}_{\mathrm{dep}}\widehat{R}_{\mathrm{dep}}^2+\stackrel{~}{\xi }_0\widehat{\xi }_0^2)/(\widehat{R}_{\mathrm{dep}}^2+\widehat{\xi }_0^2)`$, or for the gyration radius $`\widehat{R}_{\mathrm{gyr}}=\widehat{R}_2/\sqrt{1+2/D}`$ and $`\stackrel{~}{R}_{\mathrm{gyr}}=\stackrel{~}{R}_2(1+2/D)/(1\nu +2/D)`$. The measured coefficients satisfy these relations within error.

## 3 Correction to scaling analysis of multiscaling

We start with the multiscaling assumption: it has been suggested  that the particle density of an $`N`$-particle cluster at distance $`xR_{\mathrm{gyr}}`$ away from the center scales with $`R_{\mathrm{gyr}}`$ with an $`x`$-dependent co-dimension:

$$g_N(xR_{\mathrm{gyr}})=C(x)R_{\mathrm{gyr}}^{d+D(x)}$$
(3)
From this scaling law $`D(x)`$ can be obtained as

$$d+D(x)=\frac{\mathrm{ln}g_{N(R_{\mathrm{gyr}})}(xR_{\mathrm{gyr}})}{\mathrm{ln}R_{\mathrm{gyr}}}|_x.$$
(4)
In direct numerical measurements  a non-trivial $`D(x)`$ was obtained, forming the basis of multiscaling claims.

Now we consider the distribution of $`r`$, the distance of the $`N`$-th particle from the seed. From the definitions in Table 1, the mean of $`r`$ is $`R_{\mathrm{dep}}`$ and its standard deviation is $`\xi _0`$. The probability density of $`r`$ can be written as

$$\frac{1}{\xi _0(N)}h\left(\frac{rR_{\mathrm{dep}}(N)}{\xi _0(N)}\right),$$
(5)
if we assume that the shape $`h`$ of the distribution is independent of $`N`$. After replacing the sum over the particles with an integral, for the particle density we obtain

$$2\pi rg_N(r)=_0^N𝑑N^{}\frac{1}{\xi _0(N^{})}h\left(\frac{rR_{\mathrm{dep}}(N^{})}{\xi _0(N^{})}\right)$$
(6)
A formula similar to this has been suggested earlier . At this point we can calculate  the function $`D(x)`$ from Eqs. (4) and (6), because we already know the correction to scaling approximation (2) of $`R_{\mathrm{gyr}}(N)`$, $`R_{\mathrm{dep}}(N)`$ and $`\xi _0(N)`$. The only extra ingredient needed is the functional form of $`h`$, which we measured directly. It is a normalized probability density of zero mean and unit variance, and as shown on Figure 3, it turns out to be well approximated by the standard normal distribution.

Now we compare $`D(x)`$ calculated with correction to scaling forms with that of earlier direct measurements on Figure 4a: the agreement is rather good. However, our method indicates (Figure 4b) that for larger size clusters $`D(x)`$ collapses to a constant: for $`N\mathrm{}`$ the radii approach pure scaling. From this we conclude that the observed “multiscaling” is a small size transient caused by the strong correction to scaling of the radii (mostly of $`\xi _0`$).

## 4 Summary

We have seen that current numerical measurements are consistent with the simple scaling picture of DLA. Subdominant terms, however, are strongly present: earlier anomalous scaling claims—including divergent length scales and multiscaling—were misled by them. The correction to scaling analysis, calculating the effect of the dominant correction, explains these earlier observations, even for a complicated quantities like $`D(x)`$. It remains a challenge for the future to predict theoretically the—so far only empirical—correction to scaling parameters.

This research has been supported by a Marie Curie Fellowship of the EC programme “Improving Human Potential” under contract number HPMF-CT-2000-00800.


text:
# Observations of r- and s-process elements in Population II stars

## 1 Introduction

Elements beyond the iron peak are formed in neutron-capture reactions. Traditionally these are split into two categories, the slow- (s-) process when beta-decay rates exceed neutron-capture rates, and the rapid- (r-) process when neutron-capture rates are high enough to produce neutron-rich unstable isotopes far from the valley of beta-stability.

The main site for the r-process is currently thought to be supernovae, though neutron-star collisions may have played a role in later stages of Galactic enrichment . The s-process has contributions from two sites: thermally-pulsing (TP-) AGB stars are responsible for the so-called ‘main’ s-process fraction, while an additional contribution from high-mass He-core-burning stars is required to explain the first neutron-capture peak (nuclei with closed neutron shells having a neutron number $`N`$ = 50; Sr, Y, Zr). A third s-process, the so-called ‘strong’ s-process, that was invoked to explain the third neutron-capture peak (closed-neutron-shell nuclei with $`N`$ = 126; Pb and Bi), appears now to be associated with metal-poor AGB stars .

In metal-poor environments, such as existed during the formation of the oldest stars in the Galaxy, the s-process has long been regarded as unimportant. This interpretation was established some 20-25 years ago. The observational work of Spite and Spite showed that the abundance of Ba dropped faster than the abundance of Eu as progressively more-metal-poor stars were observed. As Eu is an almost pure r-process element, whereas Ba can have contributions from both processes, this was interpreted as evidence for the disappearance of the s-process contribution to Ba in metal-poor stars. The theoretical basis for the framework was provided by Truran who pointed out that the lack of seed nuclei would block the s-process in the metal-poor counterparts of the objects which, at solar-metallicity, are primarily responsible for the s-process. The r-process, on the other hand, would be capable of generating its own seed nuclei even in metal-poor environments.

This basic framework was supported by various observations and modelling efforts in subsequent years. Gilroy et al.  showed that in metal-poor giants, neutron-capture elements have the same relative abundances as in the inferred r-process contribution to the solar composition. Furthermore, detailed calculations of Galactic chemical evolution  show that AGB stars would not contribute significantly to Galactic enrichment of s-process elements until \[Fe/H\] had climbed to $``$ $`1`$.

Observations of \[Ba/Eu\] ratios were extended to \[Fe/H\] $``$ $`3`$ and showed that the solar r-process ratio persisted . The star CS 22892-052, with \[Fe/H\] $``$ $`2.5`$ and having abundance ratios of r-process to iron-peak elements forty times higher than normal, likewise had a composition in which r-process abundances were seen . The relative weakness of iron-peak elements meant that a vast array of neutron-capture elements could be studied, rather than just the more easily detected ones such as Ba and Eu, and hence the existence of an r-process pattern was established far more thoroughly for this star than for any other metal-poor object. This included filling in the gaps in the coverage of elements between the neutron capture peaks .

A view developed that the r-process might be “universal”, i.e. that the element yields in the r-process were independent of the metallicity of the gas from which the star providing the site of nucleosynthesis formed. This view was supported by findings that the element ratios near the second neutron-capture peak (closed-neutron-shell nuclei with $`N`$ = 126; e.g. Ba, La) resembled the solar r-process even at low \[Fe/H\] . As we shall see below, the concept of universality is still entertained, but observations that Sr exhibits a very different behaviour  have seen the idea of universality qualified with the restriction that is does not apply to elements in the first neutron-capture peak.

CS 22892-052 soon became, and probably rightly so, one of the most talked about stars in the field of Pop. II nucleosynthesis. Less justifiably, it was often spoken of as an archetypal metal-poor star for the purposes of investigating neutron-capture nucleosynthesis, an epithet which overlooked two key observational results: (1) its fame stemmed from it being abnormally rich in neutron-capture elements, very definitely not representative of the rest of Pop. II, and (2) it was one of a growing class of carbon-rich, metal-poor stars.

Carbon-rich, metal-poor stars were found  in follow-up observations of stars originally identified as metal-poor on objective-prism spectra of the Ca H and K and H$`ϵ`$ lines. The objective prism spectra covered none of the CH or CN molecular features, so contained no information on carbon abundances. However, blue, $``$1 Å-resolution spectra obtained subsequently showed that a significant fraction of such objects contained unusually strong CH bands; up to 25% of very metal-poor stars were suspected of having significant carbon excesses. The excesses were eventually shown  to extend up to two orders of magnitude, giving \[C/Fe\] = 2.0.

In work aiming to clarify the relationship between the C and r-process enhancements in CS 22892-052, and with the hope of identifying more stars whose r-process elements could be studied in detail, Norris, Beers and Ryan set out to observe a number of C-rich stars identified from 1 Å-resolution spectra of objective-prism-selected metal-poor objects. Their early observations , and extended studies involving many of the other targets , showed not a single additional C-rich star with an r-process enhancement, but revealed instead large numbers of metal-poor, C-rich stars with clear s-process signatures (see also ), as well as others with no neutron-capture anomalies .

In the remainder of this paper, we discuss two challenges to the standard framework for understanding neutron-capture nucleosynthesis in Pop. II stars: the existence of a significant number of Pop. II stars with enhanced carbon and s-process elements, and the apparent independence of Sr from the behaviour of Ba.

## 2 C-rich, metal-poor stars

The abundant s-process elements found in many of the C-rich stars were unexpected given the strong arguments (presented above) that favoured an r-process origin for the neutron-capture species in metal-poor stars. High abundances of C and s-process elements initially suggested that AGB stars may have been active, since AGB stars are the major sites of both C and s-process element production in the Galaxy. However, two of the first few C-rich, metal-poor stars examined in detail were subgiants rather than giants, having being found originally in a study of subdwarfs selected on the basis of high-proper motion , so it was shown quickly that the stars being observed with high C and s-process abundances could not themselves be responsible for the enrichment. Since the observed stars could not be AGB stars, and the contribution of AGB stars to Galactic chemical enrichment is believed be insignificant at \[Fe/H\] $``$ $`2.5`$, the search turned for means of directly contaminating them with the products of AGB nucleosynthesis after their formation. Suspected mechanisms therefore centred on the accretion of mass lost from an AGB companion, either via Roche-lobe overflow or perhaps more commonly via wind accretion. This mechanism matches the binary mechanism responsible for the origin of Pop. II CH stars and Pop. I Ba stars .

If this mechanism were to be able to explain the origins of these stars, we would require binary periods long enough (i.e. separations large enough) to avoid mass transfer on the RGB instead of the AGB, and yet short enough for wind accretion to be efficient. These factors imply periods in the range 1-10 yr . We would expect to see the chemical signature of a low-metallicity AGB star, and radial velocity variations on the stated timescale due to the presence of a white dwarf companion, being the compact remnant of the AGB star.

A mechanism proposed for the CH stars, before their binarity and its role was discovered, was that the He flash that terminates RGB evolution, or a He shell flash on the TP-AGB, triggered a deep mixing event that refuelled the core and re-ignited core H burning, hence providing it with a second life . Although this possibility fell out of favour, the idea has been revived in a slightly different context, and mixing events triggered by the He flash in zero-metallicity (Pop. III) stars have been investigated more recently . They may provide an alternative means of producing C-rich, metal-poor objects.

Observations of s-process enhanced stars have provided mixed support for the binary-star, mass-transfer hypothesis. LP 625-44, which has \[Fe/H\] = $`2.7`$, \[C/Fe\] = 2.0 and \[Ba/Fe\] = 2.6, has been found to have a low-luminosity companion with an orbital period $`P`$ $`{}_{}{}^{}{}_{}{}^{>}`$ 13 yr (; see also Figure 1 below). Furthermore, the s-process abundance ratios in this object from Ba to Pb are consistent with production in a low-metallicity AGB star (where the efficiency of the s-process is treated as a free parameter) . The observations suggest a low mass for the AGB progenitor.

However, LP 706-7, which has almost identical abundance ratios to LP 625-44 (\[Fe/H\] = $`2.7`$, \[C/Fe\] = 2.0 and \[Ba/Fe\] = 2.0) shows no significant radial-velocity variations over 11 years of observation; see Figure 1. LP 706-7 has an effective temperature of 6250 K and lies on the subgiant branch between the main-sequence turnoff and the base of the giant branch. Like most other Pop. II stars in this phase of evolution, it retains the primordial Li abundance. This last point is very difficult to reconcile with either a mass-transfer or mixing scenario for the origin of its C and s-process excesses. The reason is that Li is extremely fragile, and survives only in the outermost few percent (by mass) of the star. Any process which involves mixing of the surface with interior layers, or requires that the surface layers originated in a more evolved star, would fail to account for the existence or primordial abundance of this element. Although Li can be produced in evolved stars, producing it at just the right level to match the primordial value would be a severe challenge, and hence is unlikely to explain this case.

Observations of other C-rich, metal-poor stars close to the main-sequence turnoff  likewise show a significant lack of radial-velocity variations. This, combined with the problem of Li survival in LP 706-7, emphasises that while the AGB mass-transfer explanation for the C-rich, s-process-rich stars may be attractive in that it fits with established phenomena in the evolution of stars, clearly it is not the answer for all of the s-process-rich, C-rich Pop. II stars. Even if AGB stars can produce the element ratios seen in these objects, the problem of how it is incorporated into the atmospheres of these objects is not yet resolved.

## 3 Neutron-capture elements in C-normal, very metal-poor stars

As discussed in the introduction, wide variations (up to a factor of 100) in the \[Sr/Fe\] and \[Sr/Ba\] ratios of stars with \[Fe/H\] $`<`$ $`3`$ have been observed . Star-to-star variations in \[$`r`$/Fe\] have also been seen for a range of other r-process elements, e.g. , but these seem to only rarely reach large extremes, whereas Sr regularly exhibits large variations. It was the decoupling of the behaviour of Sr from that of Ba which first indicated that a universal r-process could not exist. Current evidence is consistent with a universal r-process for stable elements at the second neutron-capture peak and above, but clearly not for the first neutron-capture peak. At present there are too few observations of elements between the first and second peaks in stars with \[Fe/H\] $`<`$ $`3`$ to know at which atomic number universality is established. That is, we do not know whether silver correlates with Sr, with Ba, or with neither in the most metal-poor stars.

The lack of data on Ag may be expected since observations of the lines of this element, which are in the UV, are extremely rare. Europium, on the other hand, is frequently measured in the blue spectra of Pop. II stars, but the lines become very weak at \[Fe/H\] $`<`$ $`3`$. Unfortunately, that is precisely where the Sr anomalies become most evident. The result of this lack of data on Eu in very metal poor stars is that we do not usually have measurements much above the second r-process peak in the most metal-poor stars. As Eu is almost a pure r-process element, it avoids any ambiguity associated with the s- and r-process fractions of Ba. The origin of Ba was called into question when analysis of its hyperfine structure suggested an odd-to-even isotope ratio more in keeping with the s-process than the r-process , though this result has been challenged by new observations . More importantly, Eu probes a different region of the chart of nuclides than Ba.

The lack of observations of Eu has until very recently  meant that the region with \[Fe/H\] $`<`$ $`3`$ was devoid of Eu detections, constrained only by quite high upper limits on \[Eu/Fe\]. In addition to the new data cited, observations with the new HDS spectrograph on the Subaru telescope have permitted an extension of the observations down towards \[Fe/H\] = $`3.4`$. These observations  confirm that \[Eu/Fe\] is only slightly sub-solar, $``$ $`0.2`$. This, combined with the previously known \[Ba/Fe\] decline at low metallicity, indicates that \[Ba/Eu\] is close to the solar r-process fraction in these objects. In contrast, \[Sr/Eu\] remains above the solar r-process fraction.

This finding, combined with the high and unpredictable \[Sr/Fe\] values in such metal-poor stars, adds weight to the suggestion that the r-process sets a lower limit on the production of Sr, which one may infer from the \[Sr/Ba\] ratio that would be found if the r-process contribution were truly universal even for elements at the first neutron-capture peak. On top of this, an additional Sr source is required that appears to be highly variable from one nucleosynthesis site to the next .

Possible mechanisms  are split between the known and the unknown. The former category includes the weak s-process, which is known to produce only first neutron-capture-peak elements, but which requires pre-existing metals both as the 22Ne neutron source and for the neutron-capture seed nuclei, and hence is expected to be ineffective in very metal-poor stars. The alpha-process (alpha-rich freeze-out) may be an alternative if nucleosynthesis can extend to Sr in significant numbers. Unknown processes that may mimic the weak s-process yields include a weak-r-process , and possibly a neutron-capture process that resembles neither of the traditional extremes (s- and r-) of neutron-capture timescales.

An element that may help resolve the origin of the Sr excesses is zinc. At atomic number $`Z_{\mathrm{Zn}}`$ = 30, this lies between the iron peak ($`Z_{\mathrm{Fe}}`$ = 26) and Sr ($`Z_{\mathrm{Sr}}`$ = 38). Contributions to its nucleosynthesis come from both the alpha-process and the weak-s-process, and hence it may be able to indicate whether one or other of these processes has been active.

The possible link between Zn and Sr was first tested with the star CS 22987-008 , which was previously discovered to be very Sr-rich . The star was indeed found to have \[Zn/Fe\] above solar, but around the same time evidence emerged that so may large numbers of very metal-poor stars . The comparison in Figure 2 shows data for CS 22897-008  and the normal star HD 140283  against metal-poor giants investigated by Johnson and Bolte . Although  could not calculate Sr abundances, they found high \[Y/Ba\] values in several of their objects. One might expect these same objects to have high \[Sr/Ba\] values; this could be tested. The high-Y stars have been highlighted in the figure. Although the genuinely Sr-rich giant, CS 22897-008, has a \[Zn/Fe\] value even higher than the rising trend, the high-\[Y/Ba\] stars do not. If high-Sr and high-Y stars do belong together, then the lack of distinction between the Y-rich and Y-normal stars in the \[Zn/Fe\] vs \[Fe/H\] diagram may indicate that the high Zn abundance of CS 22897-008 is not related to its high Sr. Analysis of these objects is continuing, but it appears that the interpretation of the observations will not be simple! The nucleosynthesis of elements above the iron peak retains at least some of its secrets 20 years after the framework for understanding metal-poor stars was established.

## 4 Acknowledgments

This paper reflects collaborations and discussions with numerous scientists, especially W. Aoki, T. C. Beers, L. A. J. Blake, R. Gallino, Y. Ishimaru, J. E. Norris, S. Tsangarides, and S. Wanajo, to whom I am grateful for their contributions. I am also very gratefully to the organising committee for financial assistance that made my attendance possible.


text:
# 1 Introduction

## 1 Introduction

This paper is intended to be an introduction to the use of KtJet. For a detailed explanation of the $`K_{}`$ clustering algorithm, its properties, and the physics motivation behind it, the reader is referred to  and references therein. The KtJet package implements all the features of the Fortran implementation of the $`K_{}`$algorithm . The philosophy throughout has been to design an interface that users of the Fortran code will recognise, but which exploits the advantages of object-oriented design. Therefore the names and functions of the input parameters are retained wherever possible. The KtJet library, examples of its use and detailed documentation are all available on the KtJet website, http://www.ktjet.org/, which also provides a link to a CVS repository from which the latest version of the source code may be obtained.

## 2 The $`𝑲_{\mathbf{}}`$ algorithm

The $`K_{}`$algorithm can function in several distinct ways depending on the nature of the colliding beams and the physics to be studied. We deal first with hadron-hadron and lepton-hadron collisions. There are two modes of operation, the ‘inclusive’ and ‘exclusive’ modes. The difference between the two cases, described in sections 2.1 and 2.2, is in the definition of the hard final state jets, and the separation of these jets from the beam remnants.

For the analysis of $`e^+e^{}`$data, in contrast, there is of course no concept of beam remnants. The algorithm proceeds in a very similar way to that employed in a jet substructure analysis. We deal with these two cases together in section 2.3.

### 2.1 The inclusive mode

The algorithm proceeds as follows:

1. For every final state object111Final state objects could be, for example, partons, hadrons, calorimeter cells, tracks etc. $`h_k`$ and for every pair $`h_k`$ and $`h_l`$, compute the resolution variables $`d_{kB}`$ and $`d_{kl}`$. The precise definition of these variables can be chosen by the user (using the parameter angle), and will be described in detail in section 2.4. They always have the property that in the small angle limit they reduce respectively to the squared relative transverse momentum of the object with respect to the beam direction, and the squared relative transverse momentum of one object with respect to the other.

$`d_{kB}`$ $``$ $`E_k^2\theta _{kB}^2`$
$``$ $`k_{kB}^2,\mathrm{for}\theta 0`$
$`d_{kl}`$ $``$ $`\mathrm{min}(E_k^2,E_l^2)\theta _{kl}^2`$
$``$ $`k_{kl}^2,\mathrm{for}\theta _{kl}0`$
At this stage, a dimensionless parameter $`R`$ is introduced, which plays a radius-like role in defining the extent of the jets. This is usually set to 1.0.
2. Scale the $`d_{kB}`$ by $`R^2`$

$$d_k=d_{kB}R^2$$
(3)
3. Find the smallest value among the $`d_k`$ and $`d_{kl}`$. If a $`d_{kl}`$ is the smallest, $`h_k`$ and $`h_l`$ are combined into a single object with momentum $`p_{(kl)}`$ according to a user specified recombination scheme (parameter recom) which will be described in section 2.5. As an example, recom = 1 would correspond to 4-vector addition. If a $`d_k`$ is the smallest, object $`k`$ is defined to be a jet and is removed from the list of objects to be merged.
4. Repeat until all objects have been included in jets.

### 2.2 The exclusive mode

In this mode, the algorithm separates the ‘hard final state’ from the soft ‘beam remnants’ explicitly. Jets are defined in the hard final state by a stopping parameter $`d_{\mathrm{𝑐𝑢𝑡}}`$ with dimensions of energy squared.

1. $`d_{kB}`$ and $`d_{kl}`$ are defined as in section 2.1.
2. Find $`d_{min}`$, the smallest value among the $`d_{kB}`$ and $`d_{kl}`$. If $`d_{min}>d_{\mathrm{𝑐𝑢𝑡}}`$, all remaining objects in the event will be classified as jets, and the algorithm is complete.
3. If a $`d_{kl}`$ is the smallest, $`h_k`$ and $`h_l`$ are combined into a single object with momentum $`p_{(kl)}`$ according to a user specified recombination scheme (parameter recom). If a $`d_{kB}`$ is the smallest, object $`k`$ is included in a ‘beam jet’ and removed from the list. Go back to (2).

Note that the stopping parameter $`d_{\mathrm{𝑐𝑢𝑡}}`$ defines the hard scale of the process ($`\mathrm{\Lambda }_{QCD}^2d_{\mathrm{𝑐𝑢𝑡}}s`$, where $`\sqrt{s}`$ is the centre-of-mass energy). One can therefore think of the hard subprocess being factorised from the low-$`p_{}`$ scattering fragments, which are thrown away into the beam jets.

Instead of specifying the stopping scale $`d_{\mathrm{𝑐𝑢𝑡}}`$, one can choose to stop merging when a given number of jets is reached.

### 2.3 Subjet analysis and $`𝒆^\mathbf{+}𝒆^{\mathbf{}}`$ mode

Once the hadronic final state has been decomposed into jets, the structure of the jets themselves can be investigated. This is commonly known as a subjet analysis. The procedure is physically (and practically) identical to that employed in the analysis of an $`e^+e^{}`$event. In the $`e^+e^{}`$case, all final state objects are used as input to the algorithm. In a subjet analysis on a particular jet, only the particles within that jet are used as input. The steps are as follows;

1. Define a resolution parameter

$$y_{\mathrm{𝑐𝑢𝑡}}=Q_0^2/E_{\mathrm{𝑐𝑢𝑡}}^2.$$
(4)
The parameter $`E_{\mathrm{𝑐𝑢𝑡}}^2`$ may be chosen by the user, but is conventionally taken to be the square of the total energy of the $`e^+e^{}`$event (or the $`p_t`$ of the jet) in the frame in which the algorithm is run. These are the default settings in KtJet.
2. For each pair of objects $`h_k`$ and $`h_l`$, construct the rescaled resolution variable $`y_{kl}`$

$$y_{kl}=d_{kl}/E_{\mathrm{𝑐𝑢𝑡}}^2$$
(5)
where $`d_{kl}`$ is defined as in section 2.1.
3. Find $`y_{min}`$, the smallest of the $`y_{kl}`$. If $`y_{\mathrm{𝑚𝑖𝑛}}<y_{\mathrm{𝑐𝑢𝑡}}`$, $`h_k`$ and $`h_l`$ are combined into a single object with momentum $`p_{(kl)}`$ according to a user specified recombination scheme (parameter recom).
4. Repeat the procedure until all pairs of objects have $`y_{kl}>y_{\mathrm{𝑐𝑢𝑡}}`$. The remaining objects are called subjets.

As described in section 2.2, one can choose to stop merging when a given number of jets is reached, instead of specifying the stopping scale $`y_{\mathrm{𝑐𝑢𝑡}}`$.

### 2.4 Jet resolution variables

The small-angle behaviour of the resolution variables $`d_{kB}`$ and $`d_{kl}`$ given in equations 1 and 1 controls the behaviour of the $`K_{}`$algorithm in the soft and collinear limits. Away from the collinear limit, the user has a choice of angular definitions, controlled by the parameter angle. Note that the resolution variables are defined such that they remain monotonic functions of angle at large angles. Here, and throughout this paper, we use $`\eta `$ to denote true rapidity, i.e. $`\eta =\frac{1}{2}\mathrm{ln}\frac{E+p_z}{Ep_z}`$.

#### 2.4.1 The angular scheme, angle = 1

$`d_{kB}`$ $`=`$ $`2E_k^2(1cos\theta _{kB}),`$
$`d_{kl}`$ $`=`$ $`2\mathrm{min}(E_k^2,E_l^2)(1cos\theta _{kl}).`$ (6)
This definition was originally proposed for jets in DIS  (where the $`K_{}`$algorithm should be implemented in the Breit frame).

#### 2.4.2 The $`\mathrm{\Delta }R`$ scheme, angle = 2

$`d_{kB}`$ $`=`$ $`p_{tk}^2,`$
$`d_{kl}`$ $`=`$ $`\mathrm{min}(p_{tk}^2,p_{tl}^2)R_{kl}^2,`$ (7)
where

$$R_{kl}^2=(\eta _k\eta _l)^2+(\varphi _k\varphi _l)^2.$$
(8)
This definition corresponds to that used in cone algorithms, and is the most common choice for hadron-hadron collisions.

#### 2.4.3 The QCD emission scheme, angle = 3

An alternative definition of $`R_{kl}^2`$, motivated by the form of the QCD matrix elements for multi-parton emissions:

$$R_{kl}^2=2[cosh(\eta _k\eta _l)cos(\varphi _k\varphi _l)].$$
(9)
### 2.5 Recombination schemes

The recombination scheme defines how two objects $`h_k`$ and $`h_l`$ are merged into a single object with 4-momentum $`p_{(kl)}`$. Five schemes are available, selected by the parameter recom.

#### 2.5.1 The $`𝑬`$ scheme, recom=1

Simple 4-vector addition:

$$p_{kl}=p_k+p_l.$$
(10)
This procedure results in massive final state jets.

#### 2.5.2 The $`𝒑_𝒕`$ scheme, recom=2

$`p_{t(kl)}`$ $`=`$ $`p_{tk}+p_{tl},`$
$`\eta _{kl}`$ $`=`$ $`{\displaystyle \frac{p_{tk}\eta _k+p_{tl}\eta _l}{p_{t(kl)}}},`$
$`\varphi _{kl}`$ $`=`$ $`{\displaystyle \frac{p_{tk}\varphi _k+p_{tl}\varphi _l}{p_{t(kl)}}}.`$ (11)
This definition constrains only the 3 spatial components of the object’s 4-vector. The combined object is made massless by setting its energy equal to the magnitude of its 3-momentum. If massive objects are input (for instance particle 4-vectors from a simulated event) they are made massless in the same way before the $`d_{kB}`$ and $`d_{kl}`$ are calculated. Compare the $`E_t`$ scheme (section 2.5.4).

#### 2.5.3 The $`𝒑_𝒕^\mathrm{𝟐}`$ scheme, recom=3

$`p_{t(kl)}`$ $`=`$ $`p_{tk}+p_{tl},`$
$`\eta _{kl}`$ $`=`$ $`{\displaystyle \frac{p_{tk}^2\eta _k+p_{tl}^2\eta _l}{p_{tk}^2+p_{tl}^2}},`$
$`\varphi _{kl}`$ $`=`$ $`{\displaystyle \frac{p_{tk}^2\varphi _k+p_{tl}^2\varphi _l}{p_{tk}^2+p_{tl}^2}}.`$ (12)
This definition constrains only the 3 spatial components of the object’s 4-vector. The energy is made equal to to the magnitude of its 3-momentum, thus making the combined object massless. Note also that this definition of the $`p_t^2`$ scheme is that used in the Fortran implementation of the algorithm . It is not equivalent to the monotonic $`p_t^2`$ scheme defined in equations 16 and 17 of reference . We discuss the issue of monotonicity further in section 3.1.1.

#### 2.5.4 The $`𝑬_𝒕`$ scheme, recom=4

$`E_{t(kl)}`$ $`=`$ $`E_{tk}+E_{tl},`$
$`\eta _{kl}`$ $`=`$ $`{\displaystyle \frac{E_{tk}\eta _k+E_{tl}\eta _l}{E_{t(kl)}}},`$
$`\varphi _{kl}`$ $`=`$ $`{\displaystyle \frac{E_{tk}\varphi _k+E_{tl}\varphi _l}{E_{t(kl)}}}.`$ (13)
For massless input objects this definition is identical to the $`p_t`$ scheme (section 2.5.2). It differs solely in the way it deals with massive input objects. If massive objects are input, the $`p_t`$ scheme uses their transverse momentum, whereas the $`E_t`$ scheme uses the transverse energy $`E\mathrm{sin}\theta `$. This can have a significant effect for low $`E_t`$ jets, particularly for steeply falling distributions close to a cut-off value. All combined objects are massless in both cases.

#### 2.5.5 The $`𝑬_𝒕^\mathrm{𝟐}`$ scheme, recom=5

$`E_{t(kl)}`$ $`=`$ $`E_{tk}+E_{tl},`$
$`\eta _{kl}`$ $`=`$ $`{\displaystyle \frac{E_{tk}^2\eta _k+E_{tl}^2\eta _l}{E_{tk}^2+E_{tl}^2}},`$
$`\varphi _{kl}`$ $`=`$ $`{\displaystyle \frac{E_{tk}^2\varphi _k+E_{tl}^2\varphi _l}{E_{tk}^2+E_{tl}^2}}.`$ (14)
This bears the same relationship to the $`p_t^2`$ scheme that the $`E_t`$ scheme bears to the $`p_t`$. All combined objects are massless.

## 3 Practical implementation of the algorithm

KtJet uses the HepLorentzVector class of the CLHEP package . The input 4-vectors (the objects on which the algorithm will be run) may be HepLorentzVectors or KtLorentzVectors. The KtLorentzVector class, which inherits from HepLorentzVector, carries an internal index which allows the user to determine, for example, to which final state jet a particular input particle belongs. Examples of the use of KtLorentzVector are given in section 3.2. The output 4-vectors (the jets) are instantiated as KtLorentzVectors. The KtLorentzVector class has the constructors

```
KtLorentzVector(const HepLorentzVector &);
KtLorentzVector(const KtLorentzVector &);
KtLorentzVector(float Px, float Py, float Pz, float E);

```

The first step in running the algorithm is the creation of a KtEvent object. There are separate constructors to run the algorithm in inclusive and exclusive modes respectively:
KtEvent(const std::vector\<KtLorentzVector\> particles &, int type, int angle, int recom, float rparameter);
KtEvent(const std::vector\<KtLorentzVector\> particles &, int type, int angle, int recom);
There are also versions of these constructors taking vectors of HepLorentzVectors instead of KtLorentzVectors. The parameter type should be set to the colliding beam type, as defined in table 1. The parameters angle and recom are described in sections 2.4 and 2.5. rparameter is the parameter $`R`$ defined in section 2.1 and should almost always be set to $`1.0`$, although values smaller than $`1.0`$ have been suggested for certain applications .

### 3.1 KtEvent methods

For the inclusive case, the final state jets are now fully defined, and the jets can be recovered using the methods described below. In the exclusive case, the final state jets are defined by setting either the stopping parameter $`d_{\mathrm{𝑐𝑢𝑡}}`$:
void findJetsD(float dcut);
or by forcing the final state to decompose into $`N`$ jets:
void findJetsN(int N);
The $`d_{min}`$ value at which the final state changes from $`N+1`$ to $`N`$ jets can be recovered using the method
float getDMerge(int N) const;

In the case of an $`e^+e^{}`$analysis, the variable $`y_{\mathrm{𝑐𝑢𝑡}}`$ may be used. The final state jets are defined by the method
void findJetsY(float yCut);
and the $`y_{\mathrm{𝑚𝑖𝑛}}`$ value at which the final state changes from $`N+1`$ to $`N`$ jets can be recovered using the method
float getYMerge(int N) const;
By default, $`y_{\mathrm{𝑐𝑢𝑡}}`$ is defined as in equation 4, where $`E_{\mathrm{𝑐𝑢𝑡}}`$ is taken to be the total energy in the event. The $`E_{\mathrm{𝑐𝑢𝑡}}`$ value can be set by the user if required:
void setECut(float ECut);

The KtEvent object has the following methods for recovering its final state jets:
std::vector\<KtLorentzVector\> getJets() const;
Returns jets without sorting.
std::vector\<KtLorentzVector\> getJetsE() const;
Returns jets in order of decreasing energy.
std::vector\<KtLorentzVector\> getJetsEt() const;
Returns jets in order of decreasing transverse energy.
std::vector\<KtLorentzVector\> getJetsPt() const;
Returns jets in order of decreasing transverse momentum.
std::vector\<KtLorentzVector\> getJetsRapidity() const;
Returns jets in order of decreasing rapidity.
std::vector\<KtLorentzVector\> getJetsEta() const;
Returns jets in order of decreasing pseudorapidity.
KtLorentzVector getJet(const KtLorentzVector &) const;
Returns the final state jet which contained the input particle KtLorentzVector. Note that this method is only available if KtLorentzVectors were used as input in the KtEvent constructor.
std::vector\<const KtLorentzVector\*\> & getConstituents() const;
Returns a vector of pointers to all the input objects in the KtEvent.
std::vector\<KtLorentzVector\> copyConstituents() const;
Returns a vector of copies of the input objects in the KtEvent.
int getNConstituents() const;
Returns the number of objects in the KtEvent.
The parameters of a particular KtEvent can be retrieved using the methods
float getETot() const;
int getType() const;
int getAngle() const;
int getRecom() const;
float getECut() const;
bool isInclusive() const;

#### 3.1.1 Monotonicity

Although in practice most users will not have to face the issue, it is worth bearing in mind that recombination schemes 1 to 5 are not guaranteed to lead to monotonic resolution variables, i.e. it is NOT necessarily true that $`\mathrm{min}\{d_{kB},d_{kl}\}\mathrm{min}\{d_{kB}^{},d_{kl}^{}\}`$, where $`d_k`$ and $`d_k^{}`$ are the resolution variables before and after recombination respectively. This means that, physically, the question “How many jets are there at a particular scale $`d`$?” may not have a unique answer. In KtJet, as should be clear from section 2.2, for a particular value of $`d_{\mathrm{𝑐𝑢𝑡}}`$ set by findJetsD, the largest value of $`N`$ (i.e. the first place at which $`d_{\mathrm{𝑚𝑖𝑛}}>d_{\mathrm{𝑐𝑢𝑡}}`$) will be returned by getNJets. Similarly, there may be no value of $`d_{\mathrm{𝑐𝑢𝑡}}`$ for which a particular event has $`N`$ jets. There will however always be a $`d_{min}`$ value at which the event changed from $`N+1`$ to $`N`$ jets, which is the value returned by getDMerge. The above discussion also applies to the $`y`$ variables and associated methods.

### 3.2 KtLorentzVector Methods

These methods would normally be used to investigate the structure and constituents of final state jets defined in a particular KtEvent.
const std::vector\<const KtLorentzVector\*\> & getConstituents() const;
Returns a reference to the vector of pointers to all the objects in the KtLorentzVector.
std::vector\<KtLorentzVector\> copyConstituents() const;
Returns a vector of copies of the objects in the KtLorentzVector.
int getNConstituents() const;
Returns the number of objects in the KtLorentzVector.
bool contains (const KtLorentzVector &) const;
Check if a jet contains a particular object. For example,
if (JET.contains(PARTICLE)) { /\* do something \*/ }
where JET is the KtLorentzVector of a final state jet and PARTICLE is the KtLorentzVector of an input object in the KtEvent.
KtLorentzVector & operator+= (const KtLorentzVector &);
Adds a KtLorentzVector constituent to a jet using the $`E`$ scheme (4-vector addition) and maintains an internal record of constituents (so for example the getConstituents method will work on the resulting KtLorentzVector).
void add (const KtLorentzVector &, int recom);
Adds a KtLorentzVector constituent to a jet using any of the available recombination schemes and maintains an internal record of constituents.

Note that the momentum of a jet (KtLorentzVector) will be given according to the recombination scheme used in its construction. If the user wishes to reconstruct the momenta according to a different recombination scheme (for example, to recover the mass for jets which were found using a massless recombination scheme) the getConstituents() or copyConstituents() methods may be used. The constituents can then be recombined in a new scheme using the add method.

## 4 Subjet analysis

A subjet analysis is performed on a particular final state jet by constructing a new KtEvent object:
KtEvent(const KtLorentzVector jet &, int type, int angle, int recom);
The subsequent analysis and methods are identical to those described above for the $`e^+e^{}`$ case: for example, subjets can be defined using the findJetsY or findJetsN methods. $`E_{\mathrm{𝑐𝑢𝑡}}`$ is taken to be the transverse momentum of the jet, although this can be changed by the user using the setECut method.

## 5 Advanced features of the KtLorentzVector class

The KtLorentzVector class has several features which are used internally by KtJet, but which may be useful for the user.

A KtLorentzVector can be constructed in two different ways. When instantiated using the constructor KtLorentzVector(const HepLorentzVector &) (or KtLorentzVector(float Px,float Py,float Pz,float E)), the KtLorentzVector is simply a copy of a single HepLorentzVector with an internal index added, i.e. it has no constituents. This index is then used, for example, in the getJet method, to ascertain to which final state jet a particular input KtLorentzVector belongs. There are two methods which compare the index of two KtLorentzVectors:

```
bool operator== (const KtLorentzVector &) const;
bool operator!= (const KtLorentzVector &) const;

```

These would be useful if, for example, one wanted to ascertain whether the highest $`p_t`$ particle in a particular jet was a given input KtLorentzVector. A KtLorentzVector may also be instantiated using the constructor
KtLorentzVector();
The KtLorentzVector can then be built up by adding other KtLorentzVectors to it using either the assignment operator += or the add method described in section 3.2. The KtLorentzVector will then also carry a list of pointers to the constituent KtLorentzVectors. To find out whether a KtLorentzVector has constituents, use the method
bool isJet() const;
which will be true for a KtLorentzVector with a constituent list.

## 6 Using KtJet

In this section we demonstrate a simple example analysis using the KtJet library. We run the inclusive $`K_{}`$algorithm on a proton-proton collision and do a subjet analysis on the highest $`E_t`$ jet in the event.

```
// This header file must be included.
#include "KtJet.h"

std::vector<KtLorentzVector> jetvec;

// Loop over input particles.
for (int i=0; i<npart; i++) {
  KtLorentzVector  r = KtLorentzVector(p[i][0],p[i][1],p[i][2],p[i][3]);
  jetvec.push_back(r);
}
// We now have a vector of KtLorentzVector, jetvec,
// containing all the particle 4-vectors, p, in the event.

// Run the inclusive KT algorithm in PP mode using the covariant E-scheme
// (type=4, angle=2, recom=1, rparameter=1).
KtEvent ev(jetvec,4,2,1,1.0);

// Get jets, sorted in Et.
std::vector<KtLorentzVector> jetsEt = ev.getJetsEt();

// Perform subject analysis on highest Et jet (angle=2, recom=1).
KtEvent jet1(jetsEt[0],2,1);

// Decompose jet1 into 2 subjets.
jet1.findJetsN(2);

// Get the KtLorentzVectors of the 2 subjets
std::vector<KtLorentzVector> subjetsEt = jet1.getJetsEt();

// Write out the pseudorapidity of the highest Et subjet.
cout << subjetsEt[0].eta() << endl;

```

## 7 Optimisation issues

The time taken to process events grows with the multiplicity $`n`$ of the input objects as $`n^3`$. As the multiplicity grows with energy, and current and future colliders have ever higher event rates, timing issues are therefore of increasing practical importance.

A few compilation options are worth considering if speed is critical. A compile option to switch from single to double precision arithmetic is available, but single precision is used by default as this increases the speed by around 10%. In gcc, the option -O2 is essential; higher or lower optimisation levels seriously degrade performance. Choice of architecture is also important. For example using the gcc compiler flag -march=athlon on an Athlon processor gives a performance gain of around 5%. More details on performance issues are given at http://www.ktjet.org/.

## 8 Adding new functionality

KtJet interfaces with its recombination and jet resolution schemes via purely abstract base classes called KtRecom and KtDistance respectively. The packaged code comes supplied with 5 Recombination and 3 jet resolution schemes, shown in table 2 along with their angle and recom flags.

Users can define their own recombination and jet resolution scheme classes which inherit from the base classes KtRecom and KtDistance and use these instead of the supplied schemes.

To use their own scheme the user needs to pass pointers to the base classes, instantiated as their own scheme objects, in the KtEvent constructor instead of the integer flags.

As an example, if the user defines a new jet resolution scheme in the class KtDistanceNew, and wishes to use the $`E`$ recombination scheme in an $`ep`$ collision, the inclusive KtEvent constructor should be called as follows:

```
KtDistance* distance_scheme = new KtDistanceNew(2);
KtEvent ev(particles, 2, distance_scheme, 1, rparameter);

```

If the user defines a new recombination scheme in the class KtRecomNew, and wishes to use the angle=2 scheme in a $`pp`$ collision, the inclusive KtEvent constructor should be called as follows:

```
KtRecom* recom_scheme = new KtRecomNew();
KtEvent ev(particles, 4, 2, recom_scheme, rparameter);

```

Users are invited to submit any such additions to the authors for inclusion in future releases.

## Acknowledgements

We would like to thank Mike Seymour for useful discussions, suggestions and encouragement throughout the project. This work was funded in the UK by PPARC and HEFCE.