text:
# Skyrmions on the Two-Sphere.

## 1 Introduction

Topological solitons are usually studied on flat space, but there are several reasons why the curved-space setting is interesting (apart from the most obvious one of possible cosmological applications). First, it often reveals useful mathematical features; there are now two length scales, namely the size of the solitons and the radius of curvature of the underlying geometry, and there is an interplay between these scales. For monopole-like objects such as as vortices and Yang-Mills-Higgs monopoles, one can work on hyperbolic space: in the case of vortices, the equations for static vortices become integrable for a particular value of the curvature , ; and for YMH monopoles, the system is simpler on hyperbolic space than on flat space . For textures, such as the Skyrmions discussed in this paper, one may study the problem on either hyperbolic spaces or spheres (ie on spaces of constant negative or positive curvature), and this similarly throws up interesting geometrical features , , , , , , .

There are also more physical reasons for studying textures on compact spaces, namely that this models a finite density of solitons, and the transistion between the high-density and the low-density phases. One may do this without introducing curvature, by taking space to be a flat torus; but it is in some ways simpler to use a sphere. This was first done for the Skyrme model in three spatial dimensions, which is a model of nuclear matter. But there has also been interest in the two-dimensional version, both as a ‘toy’ model and also because it has several potential applications in condensed-matter physics. In this paper, we study the Skyrme model on the two-sphere of radius $`L`$, and in particular the dependence on the ratio $`\beta =L/(\mathrm{size}\mathrm{of}\mathrm{Skyrmion})`$.

In the three-dimensional Skyrme model (without a pion mass term), the Skyrmion for small $`\beta `$ is homogeneously spread out over space; when $`\beta `$ reaches a certain critical value, the Skyrmion localizes around a point , . If the winding number is greater than one, or if a zeroth-order potential term (such as a mass term) is present, then this transistion is less clear-cut. In the two-dimensional case, such a potential term is obligatory, and so this smearing effect is inevitable (cf ). The quantity $`\mathrm{\Phi }`$ which monitors the deviation from homogeneity tends to zero as $`\beta `$ does, but $`\mathrm{\Phi }`$ is never quite zero, and so is not an order parameter in the strict sense. There is, however, another quantity $`\mathrm{\Psi }`$ which does serve as an order parameter, namely one which monitors the breakdown of a certain reflection symmetry , ; and the associated phase transition is sharp (as long as the potential has this symmetry, which was not the case for the system studied in ). We shall investigate these phenomena, for several different potentials, by finding rotationally-symmetric solutions numerically.

## 2 Skyrme models on the two-sphere

In this section, we review the Skyrme system on the two-sphere, the topological Bogomol’nyi bound on its energy, and the imposition of rotational symmetry. Let $`M`$ denote the standard 2-sphere of radius $`L`$, with local coordinates $`(x^1,x^2)`$ and metric $`ds^2=g_{jk}dx^jdx^k`$. The area element on $`M`$ is $`\omega =hdx^1dx^2`$, where $`h=[det(g_{jk})]^{1/2}`$. The Skyrme field $`\varphi `$ is a map from $`M`$ to the unit sphere $`S^2`$, and may be thought of as a unit vector field $`\stackrel{}{\varphi }(x^j)`$ (so $`\stackrel{}{\varphi }\stackrel{}{\varphi }=1`$). We are considering only the static problem, so $`\stackrel{}{\varphi }`$ depends only on the spatial variables $`x^j`$, and not on time. Let $`_j\stackrel{}{\varphi }`$ denote the partial derivative $`\stackrel{}{\varphi }/x^j`$. The winding number (degree) of $`\stackrel{}{\varphi }`$ is denoted by the integer $`k`$. If we define the function $`\rho [\stackrel{}{\varphi }]`$ as the triple scalar product $`\rho [\stackrel{}{\varphi }]=2\stackrel{}{\varphi }(_1\stackrel{}{\varphi })\times (_2\stackrel{}{\varphi })`$, then (assuming that the $`x^j`$ cover almost all of $`M`$) we have

$$k=\frac{1}{8\pi }_M\rho [\stackrel{}{\varphi }]𝑑x^1𝑑x^2.$$
(1)
The (normalized) energy $`E_k`$ of $`\stackrel{}{\varphi }`$ is taken to be a sum of three terms $`E_k=E^{(2)}+E^{(4)}+E^{(0)}`$, with each being the integral over $`S^2`$ of the corresponding density function:

$$E^{(n)}=\frac{1}{4\pi |k|}_M^{(n)}[\stackrel{}{\varphi }]\omega $$
(2)
for $`n=2,4,0`$. The densities $`^{(n)}(\stackrel{}{\varphi })`$ are given by

$`^{(2)}[\stackrel{}{\varphi }]`$ $`=`$ $`{\displaystyle \frac{1}{2}}g^{jk}(_j\stackrel{}{\varphi })(_k\stackrel{}{\varphi }),`$ (3)
$`^{(4)}[\stackrel{}{\varphi }]`$ $`=`$ $`{\displaystyle \frac{1}{8}}\gamma h^2\rho ^2,`$ (4)
$`^{(0)}[\stackrel{}{\varphi }]`$ $`=`$ $`\alpha V(\stackrel{}{\varphi }),`$ (5)
where $`\gamma `$ and $`\alpha `$ are positive constants, and where $`V(\stackrel{}{\varphi })`$ is some suitable function of $`\stackrel{}{\varphi }`$ (not involving its derivatives). We regard the normalized energy as dimensionless. The constants $`\gamma `$ and $`\alpha `$ have dimensions of length2 and length-2 respectively; the dimensionless combination $`\alpha \gamma `$ determines the energy of a soliton, while its size is of order $`(\gamma /\alpha )^{1/4}`$. In flat space, this latter quantity simply sets the length scale; but on the sphere $`M`$, there already is a length scale $`L`$. So one has a dimensionless parameter $`\beta =L(\alpha /\gamma )^{1/4}`$, which is the ratio between the size of space and the size of the soliton. From now on, we shall set $`\gamma =1`$.

The energy of any configuration of nonzero degree satisfies a generalized Bogomol’nyi bound . First, $`E^{(2)}`$ satisfies the standard O(3) sigma-model bound $`E^{(2)}1`$. To get a bound on $`E^{(4)}+E^{(0)}`$, write $`^{(0)}=G^2/2`$ and observe that

$`E^{(4)}+E^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{8\pi |k|}}{\displaystyle _M}\left[(h^1\rho /2G)^2+h^1\rho G\right]\omega `$
$``$ $`{\displaystyle \frac{1}{8\pi |k|}}{\displaystyle _M}\rho G𝑑x^1𝑑x^2`$
$`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _{S^2}}G(\stackrel{}{\varphi })\stackrel{}{\varphi }𝑑\stackrel{}{\varphi }\times 𝑑\stackrel{}{\varphi };`$
the last equality follows from the degree theorem $`\varphi ^{}(G\mathrm{\Omega })=kG\mathrm{\Omega }`$, with $`\mathrm{\Omega }=\stackrel{}{\varphi }d\stackrel{}{\varphi }\times d\stackrel{}{\varphi }`$ being the area element on the target sphere $`S^2`$. If $`G`$ depends only on the third component $`\varphi ^3`$ of $`\stackrel{}{\varphi }`$, which is the case for all the examples that we consider, then the bound on the total energy takes the form

$`E_k`$ $``$ $`1+{\displaystyle \frac{1}{2}}{\displaystyle _1^1}G(\varphi ^3)𝑑\varphi ^3`$ (6)
$`=`$ $`1+\sqrt{\alpha /2}{\displaystyle _1^1}\sqrt{V(\varphi ^3)}𝑑\varphi ^3.`$
In this paper, we study rotationally-symmetric configurations. Let us use polar coordinates $`(\theta ,\phi )`$ on $`M`$, so that $`ds^2=L^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2)`$ and $`h=L^2\mathrm{sin}\theta `$. We impose rotational symmetry by taking $`\stackrel{}{\varphi }`$ to have the form

$$\stackrel{}{\varphi }=(\mathrm{sin}(f)\mathrm{cos}(k\phi ),\mathrm{sin}(f)\mathrm{sin}(k\phi ),\mathrm{cos}(f)),$$
(7)
where $`k`$ is an integer and $`f=f(\theta )`$. Assuming that this profile function $`f(\theta )`$ satisfies the boundary conditions $`f(0)=\pi `$ and $`f(\pi )=0`$, the configuration (7) has winding number $`k`$. Its normalized energy is

$$E_k=\frac{1}{4|k|}_0^\pi \left[(f^{})^2+k^2\left(\frac{\mathrm{sin}f}{\mathrm{sin}\theta }\right)^2+\left(\frac{kf^{}}{L}\frac{\mathrm{sin}f}{\mathrm{sin}\theta }\right)^2+2L^2^{(0)}(\mathrm{cos}f)\right]\mathrm{sin}\theta d\theta ,$$
(8)
where $`f^{}=df/d\theta `$ and where we are (as mentioned above) taking $`^{(0)}`$ to be a function of $`\varphi ^3=\mathrm{cos}f`$ only.

We are interested in the stationary points of the functional (8), and we find these numerically (with an accuracy better than $`0.1\%`$). Two independent methods were used: first, a direct minimization of (8) by a conjugate-gradient method; and secondly, solving the corresponding Euler-Lagrange equations.

One particular ‘trial’ profile which is of interest is the one corresponding to a rotationally-symmetric conformal map (cf , , ). Such a profile can be used for any value of $`k`$, but we shall restrict here to the case $`k=1`$. The profile function $`f(\theta )`$ is then given by

$$\mathrm{cot}(f/2)=\lambda \mathrm{tan}(\theta /2),$$
(9)
where $`\lambda `$ is a positive constant. Notice that if $`\lambda =1`$, then $`f(\theta )=\pi \theta `$, and this corresponds to the identity map from $`MS^2`$ to $`S^2`$. This identity map is only a solution of the equations of motion if $`^{(0)}`$ is constant (ie not depending on $`\stackrel{}{\varphi }`$). Another way of putting this is that if (and only if) $`^{(0)}`$ is constant, then Skyrmions spread out to fill the whole of space homogeneously, in the sense that the energy density $``$ is constant. The configuration (9) has $`E^{(2)}(\lambda )=1`$ and $`E^{(4)}(\lambda )=(\lambda ^4+\lambda ^2+1)/(6L^2\lambda ^2)`$. Notice that $`E^{(4)}(\lambda )`$ has its minimum at $`\lambda =1`$; but the location (and existence) of a minimum in the total energy $`E(\lambda )`$ depends on the choice of the potential term $`E^{(0)}`$. The question is: for a given $`L`$ and choice of $`E^{(0)}`$, what are the minima of $`E(\lambda )`$, and in particular are they at $`\lambda =1`$ or at $`\lambda 1`$?

We now wish to define the two quantities $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ mentioned in the introduction. The first of these is

$$\mathrm{\Phi }=(\varphi ^3)^2\frac{1}{3}=\frac{1}{4\pi L^2}_M(\varphi ^3)^2\omega \frac{1}{3}.$$
(10)
This is zero for the homogeneous configuration $`f(\theta )=\pi \theta `$, and so gives an indication of the deviation from homogeneity. The second quantity is

$$\mathrm{\Psi }=\varphi ^3^2=\left[\frac{1}{4\pi L^2}_M\varphi ^3\omega \right]^2.$$
(11)
This is zero for configurations which have the reflection symmetry $`\varphi ^3\varphi ^3`$, but nonzero if the Skyrmion is localized around the point $`\theta =0`$ or $`\theta =\pi `$.

## 3 Two asymmetric examples

In this section, we give a brief discussion of two systems, arising from two possible choices of $`^{(0)}`$. Each is ‘asymmetric’, in the sense that neither posesses the symmetry $`\varphi ^3\varphi ^3`$. So a Skyrmion in these systems will never have $`\mathrm{\Psi }=0`$; the field prefers to be near $`\stackrel{}{\varphi }=(0,0,1)`$ rather than $`\stackrel{}{\varphi }=(0,0,1)`$. The first example (for corresponding flat-space studies see , ) is motivated by the question: can one saturate the Bogomol’nyi bound (6)? We shall not give a complete analysis of this question here; but in the $`k=1`$ case, it is easy to derive the following fact: the configuration (9) saturates the Bogomol’nyi bound if and only if $`^{(0)}`$ is given by

$$^{(0)}=\frac{1}{2}\left[\frac{(\lambda ^2+1)(\lambda ^21)\varphi ^3}{2L\lambda }\right]^4.$$
(12)
(Notice that the $`\lambda =1`$ case is degenerate: $`^{(0)}`$ is then constant.) For this system, the bound is

$$E1+(\lambda ^4+\lambda ^2+1)/(3L^2\lambda ^2).$$
(13)
If, for example, we choose $`\lambda =2L`$, then the configuration (9) is a solution with energy

$$E=\frac{28L^4+4L^2+1}{12L^4}=\frac{7}{3}+O(L^2);$$
(14)
in the limit $`L\mathrm{}`$, this corresponds to ,  with their parameters $`\theta _1`$ and $`\theta _2`$ set equal to $`0.5`$ and $`1`$ respectively (note also that their expression for energy $`E`$ is $`2\pi `$ times ours). The energy (14) is a monotonic-decreasing function of $`L`$; this remains true for variations of (12) such as $`^{(0)}=(1\varphi ^3)^4`$, and for $`k>1`$.

In this system, there is a repulsive force between solitons , , . Consequently, one expects that the minimum-energy configurations for $`k>1`$ will not be rotationally-symmetric (or, equivalently, that all rotationally-symmetric $`k>1`$ configurations are unstable). So investigation of the $`k>1`$ sectors requires a full two-dimensional study.

The comments of the previous paragraph apply equally well to our other example in this section, which was the one studied in , ,  (on flat space) and  (on the sphere). For this system, we have $`V(\varphi ^3)=1\varphi ^3`$. The Bogomol’nyi bound (6) becomes $`E_k1+4\sqrt{\alpha }/3`$ (note that this is stronger than the bound given in ). Here we give a brief description of the $`k=1`$ case, taking $`\alpha =0.1`$. For the $`\lambda `$-approximation (9), we have

$$E^{(0)}(\lambda )=\frac{2L^2\alpha }{1\lambda ^2}\left[1+\frac{\lambda ^2\mathrm{log}(\lambda ^2)}{1\lambda ^2}\right],$$
(15)
which is a decreasing function of $`\lambda `$. The total energy $`E(\lambda )`$ has a minimum at $`\lambda =\lambda _{\mathrm{min}}`$, with $`\lambda _{\mathrm{min}}>1`$ depending on the parameter $`\beta =L\alpha ^{1/4}`$. In figure 1, we plot the energy $`E(\lambda _{\mathrm{min}})`$, the actual Skyrmion energy $`E=E_1`$ (the minimum of (8)), and the quantities

$`\mathrm{\Psi }`$, $`\mathrm{\Phi }`$ and $`\lambda _{\mathrm{min}}`$, all as functions of $`L`$. Notice that the lowest value of $`E`$ (when $`L=1.65`$ and $`\beta =0.93`$) is $`1.424`$, which less than $`0.2\%`$ greater than the Bogomol’nyi bound $`1.4216`$. At this value of $`L`$, we have $`\lambda _{\mathrm{min}}=1.39`$ and $`E(\lambda _{\mathrm{min}})=1.425`$. As $`L\mathrm{}`$, the energy of the Skyrmion tends to its flat-space value of $`1.549`$ .

As was noted in , there is a regime ($`\beta <1`$) where the Skyrmion is approximately spread out over $`M`$ and well-approximated by the configuration (9), and one ($`\beta >1`$) where the Skyrmion is localized and not well-approximated by (9). But the transition between these two regimes is not sharply-defined, and is not a phase transistion in the usual sense. This is partly due to the asymmetry of the system. In our next examples, we shall see a sharper transition.

## 4 Two symmetric examples

In this section, we study the systems which are defined, respectively, by $`V(\varphi ^3)=1(\varphi ^3)^2`$ and $`V(\varphi ^3)=(\varphi ^3)^2[1(\varphi ^3)^2]`$. In these systems, Skyrmions attract one another , ; and consequently the solutions with $`k>1`$ have rotational symmetry. Since these potentials have the symmetry $`\varphi ^3\varphi ^3`$, one expects (in the rotationally-symmetric case, and for large $`L`$) that there will be a stable Skyrmion solution localized at a point (either $`\theta =0`$ or $`\theta =\pi `$); and also that there exists a solution which is symmetric under $`\theta \pi \theta `$ (and which may or may not be stable). For small $`L`$, however, only the symmetric solution exists.

First take $`V(\varphi ^3)=1(\varphi ^3)^2`$, with the value $`\alpha =0.08873`$ (this corresponds to the value used in the flat-space study of , and allows a quantitative comparison with its results). The Bogomol’nyi bound (6) is $`E1+\pi \sqrt{\alpha /8}=1.3309`$. Let us look first at the $`k=1`$ sector. For the $`\lambda `$-approximation (9), we have

$$E^{(0)}(\lambda )=\frac{4\alpha L^2\lambda ^2}{(1\lambda ^2)^2}\left[2\frac{1+\lambda ^2}{1\lambda ^2}\mathrm{log}(\lambda ^2)\right],$$
(16)
which is a positive function with a maximum at $`\lambda =1`$ (recall that $`E^{(4)}(\lambda )`$ has a minimum at $`\lambda =1`$). Notice that $`E(\lambda ^1)=E(\lambda )`$; this corresponds to the $`\varphi ^3\varphi ^3`$ symmetry of $`E[\stackrel{}{\varphi }]`$ in this case. So $`\lambda =1`$ gives either a minimum or a local maximum of $`E(\lambda )`$, depending on the value of $`L`$. A simple calculation reveals that if $`\beta ^45/2`$, then

$$E(1)=1+\frac{1}{2L^2}+\frac{2L^2\alpha }{3}$$
(17)
is a minimum of $`E(\lambda )`$; while if $`\beta ^4>5/2`$, then there is a local maximum at $`\lambda =1`$, and the minima are at $`\lambda _{\mathrm{min}}>1`$ and $`1/\lambda _{\mathrm{min}}`$. Consequently, if we restrict to the special profiles (9), then there is a sharp phase transition at $`\beta ^4=5/2`$ between a homogeneous phase ($`\lambda =1`$) and one where the Skyrmion localizes at a point (in this case, one of the poles) of the sphere $`M`$. Notice also that the lowest value of $`E(1)`$, namely $`E=1+\sqrt{4\alpha /3}=1.344`$, is attained when $`\beta ^4=3/4L=1.705`$.

This is just an approximation; the true Skyrmion energy is strictly less than $`E(\lambda )`$. However, there is indeed a minimum at $`L=1.705`$ (with $`E=1.341`$, only $`0.2\%`$ less than the $`\lambda `$-approximation); and there is a transition at $`\beta ^4=5/2L=2.304`$, as one can see from figure 2.

Figure 2(a) depicts the energies of two static $`k=1`$ solutions, as functions of $`L`$. The solid curve is the energy of the ‘symmetric’ solution, ie the one satisfying $`f(\pi \theta )=\pi f(\theta )`$. This is obtained by starting with the configuration $`f(\theta )=\pi \theta `$ and iterating towards a stationary point; our procedures are able to converge to a saddle-point (which this is for large $`L`$). The dashed curve is the energy of an asymmetric solution, obtained by starting with a highly-asymmetric configuration. We see that for $`\beta ^4<5/2`$, the two solutions coincide; but for $`\beta ^4>5/2`$, the former becomes unstable (numerical experiment indicates that it is a saddle-point rather than a local minimum). The corresponding parameters $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ for the minimum-energy solution are plotted in figures 2(b) and 2(c); note that when $`\beta ^4<5/2`$, we have $`\mathrm{\Psi }=0`$ (the solution is symmetric), but $`\mathrm{\Phi }0`$ (it is not homogeneous). There is a second-order phase transition at $`\beta ^4=5/2`$. In figure 2(d), the Skyrmion energies (minima of $`E`$) are plotted for $`k=1,2,3`$. It might be noted that $`E_k`$ gets very close to the Bogomol’nyi bound, at $`L=1.7`$, $`L=2.9`$ and $`L=4.3`$ for $`k=1,2,3`$ respectively (only $`0.1\%`$ above for $`k=2`$). As $`L\mathrm{}`$, we have

$$E_11.564,E_21.405,E_31.371.$$
(18)
Now let us turn to the system with $`V(\varphi ^3)=(\varphi ^3)^2[1(\varphi ^3)^2]`$, taking $`\alpha =0.04`$ (which is consistent with the choice made in ). The Bogomol’nyi bound is $`E1+\sqrt{2\alpha }/3=1.0943`$. In the $`k=1`$ sector, the $`\lambda `$-approximation (9) has

$$E^{(0)}(\lambda )=\frac{4L^2\alpha \lambda ^2\left[(10+18\lambda ^218\lambda ^410\lambda ^6)+(3+21\lambda ^2+21\lambda ^4+3\lambda ^6)\mathrm{log}(\lambda ^2)\right]}{3(\lambda ^21)^5},$$
(19)
which has a local minimum at $`\lambda =1`$, with $`E^{(0)}(1)=2L^2\alpha /15`$. Consequently, $`E(\lambda )`$ always has a local minimum at $`\lambda =1`$. Notice that for $`L=[15/(4\alpha )]^{1/4}=3.11`$ (where the true energy has its lowest value $`1.1026`$), the approximation has energy $`E(1)=1+\sqrt{4\alpha /15}=1.1033`$.

There are two values $`L_1`$ and $`L_2`$ which are critical in the following sense. For $`LL_1`$, $`E(\lambda )`$ has only one stationary point (the minimum at $`\lambda =1`$). For $`L>L_1`$, $`E(\lambda )`$ has three local minima: a symmetric one ($`\lambda =1`$) and a degenerate asymmetric one ($`\lambda 1`$). For $`L_1<L<L_2`$, the asymmetric minimum has higher energy than the symmetric one: see figure 3(a), where $`E`$ is plotted as a function of $`\mu =(\lambda ^21)/(\lambda ^2+1)`$, for $`L=7.4`$. Finally, for $`L>L_2`$, the asymmetric minimum has the lower energy (cf. figure 3(b), for $`L=8.4`$); and the energy of the symmetric minimum tends to infinity as $`L\mathrm{}`$.

This approximate picture suggests that for $`L`$ large enough, there are two stable Skyrmion solutions (symmetric and asymmetric); but that only the latter survives in the limit as $`L\mathrm{}`$. Our numerical investigation shows that this is indeed the case. Results from minimizing (8) are presented in figure 4.

In figure 4(a), the solid curve is the energy of the symmetric $`k=1`$ solution; while the dashed curve is the energy of the asymmetric solution, which exists only for $`LL_16`$. For $`L_1<L<L_26.9`$, the energy of the asymmetric solution is greater than that of the symmetric solution. Numerical evidence indicates that each of these solutions is stable (ie is a local minimum of the energy), at least in the rotationally-symmetric class. The stability of the symmetric solution, even for large $`L`$, is a consequence of the form of the potential $`V=(\varphi ^3)^2[1(\varphi ^3)^2]`$, which favours $`\varphi ^3=0f(\theta )=\pi /2`$. It also favours $`\varphi ^3=\pm 1`$, which stablizes the asymmetric solution. Between these local minima there will be a saddle-point solution, but we have not investigated this.

Figures 4(b) and 4(c) show the quantities $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ for the minimal-energy $`k=1`$ solution; and in figure 4(d), the minimal energy is plotted for $`k=1,2,3`$. As $`L\mathrm{}`$, we have

$$E_11.281,E_21.166,E_31.143.$$
(20)
## 5 Concluding Remarks

We have studied rotationally-symmetric static solutions of the Skyrme model on the two-sphere of radius $`L`$. Even in the sector with winding number $`k`$ equal to unity, the properties of the solutions, and of the transition between the high-density (small $`L`$) and low-density (large $`L`$) phases, depends crucially on the choice of potential term.

It should be instructive to investigate the $`k>1`$ solution spaces in more detail (cf , , ). For symmetric potentials (such as those of section 4), this can be done within the rotationally-symmetric class; this is the analogue of the rational-map ansatz  in three dimensions. But for asymmetric potentials such as those of section 3, one expects a different picture: the Skyrmions should separate and form some pattern on $`S^2`$ as their optimal configuration. Flat-space studies  indicate that this pattern will, in general, be rather complicated.


format: text

text:
# Intensity limits for stationary and interacting multi-soliton complexes

## Abstract

We obtain an accurate estimate for the peak intensities of multi-soliton complexes for a Kerr-type nonlinearity in the (1+1) – dimension problem. Using exact analytical solutions of the integrable set of nonlinear Schrödinger equations, we establish a rigorous relationship between the eigenvalues of incoherently–coupled fundamental solitons and the range of admissible intensities. A clear geometrical interpretation of this effect is given.

Pulse or beam evolution under the action of competing linear dispersion and nonlinear self-action is a fundamental phenomenon which can be observed in optical fibers, nonlinear crystals or gases, atomic Bose-Einstein condensates, etc. The corresponding wave properties can often be described by the nonlinear Schrödinger equation (NLSE). In the (1+1)-dimensional case, if the homogeneous medium has a Kerr-type nonlinear response, this equation is integrable , and its solutions can be found analytically. They describe a nonlinear superposition of localized waves, called solitons, and radiation modes. Such analysis allows us to reveal some general features of nonlinear localized waves.

Localized solutions of the NLSE are among the main objects of interest for telecommunication and information processing applications. All possible localized solutions of integrable systems are composed of solitons. Hence it is important to understand their interactions and collective properties. For a single NLSE describing the evolution of a coherent field, soliton interactions are phase-sensitive, so that constructive or destructive nonlinear interference can be observed. Since the profiles of interacting solitons are distorted according to the nonlinear superposition principle, large variations of the peak intensities can occur. Indeed, it has been demonstrated that the peak intensity can vary by a factor of $`N^2`$ for $`N`$ interacting solitons . In some sense, phase sensitivity can be “amplified” due to nonlinearity, so that it is greater than in the case of linear interference between mutually coherent sources.

Interactions of fundamental solitons can also be phase-insensitive if they correspond to incoherently coupled components. Then, fundamental solitons can be composed into a stationary multi-soliton complex (MSC). To have a rough idea of an MSC, we can use the following analogy. If we think of solitons as elementary particles, then we can imagine an MSC as a more complicated structure, such as an atom. An optical incoherent soliton is an example of an MSC. These objects have been studied extensively both experimentally and theoretically during recent years (see Refs.  and references therein); some other examples of MSCs are listed in Ref. . The question of the maximal intensity of stationary and interacting MSCs is of primary interest for applications. However, to the best of our knowledge, this problem has not been addressed yet. In this work, we derive a rigorous relationship between the parameters of incoherently–coupled solitons and the range of admissible intensities for MSCs.

The evolution of spatial incoherent solitons in a medium with a Kerr-type self-focusing nonlinearity can be described by a set of NLSEs for mutually–incoherent components :

$$i\frac{u_m}{z}+\frac{1}{2}\frac{^2u_m}{x^2}+Iu_m=0,$$
(1)
where $`u_m`$ is the normalized amplitude of the $`m`$-th component ($`1mM`$), $`M`$ is the number of components, $`I=_{j=1}^M|u_j|^2`$ is the total intensity, $`z`$ is the coordinate along the direction of propagation, and $`x`$ is the transverse coordinate (we are considering a (1+1)-dimensional case).

Eqs. (1) are completely integrable by means of the inverse scattering technique (IST)  and, therefore, any localized solution can be represented as a nonlinear superposition of solitary waves. Every fundamental soliton (labeled $`j`$) is characterized by a complex eigenvalue $`k_j=r_j+i\mu _j`$, and a shift in the coordinate plane $`(x_j,z_j)`$. The simplest single–soliton solution can be written as: $`u_j(x,z)=r_j\mathrm{sech}(\beta _j)e^{i\gamma _j}`$, where $`\beta _j=r_j(\overline{x}_j\mu _j\overline{z}_j)`$, $`\gamma _j=\mu _j\overline{x}_j+(r_j^2\mu _j^2)\overline{z}_j/2`$, and $`(\overline{x}_j,\overline{z}_j)=(xx_j,zz_j)`$ are the shifted coordinates. The solution for bright MSCs consisting of $`M`$ incoherently–interacting fundamental solitons can be found by solving a set of linear equations :

$$\underset{m=1}{\overset{M}{}}\frac{e_je_m^{}u_m}{k_j+k_m^{}}+\frac{1}{2r_j}u_j=e_j,$$
(2)
where $`e_j=\chi _j\mathrm{exp}\left(\beta _j+i\gamma _j\right)`$ and the coefficients $`\chi _j`$ define the relative coordinate centers for the fundamental solitons. The seemingly straightforward process of inverting the set of linear equations results in a complicated expression.

Using the technique described in Ref. , the solution of Eqs. (2) can be obtained in a compact explicit form. However, it is not possible to find a general analytical expression for the maximum intensity, and it can only be determined by numerically solving transcendental equations. Therefore, a different approach is needed to make an analytical estimate for the peak intensity levels. In order to do this, we first multiply Eqs. (2) by $`u_j^{}/(k_ji\mu )`$, add the complex conjugate, and sum over the fundamental soliton numbers $`1jM`$ (we note that such a transformation was suggested earlier  as an intermediate step for deriving bright and dark soliton solutions on a background). Finally, we obtain the following inequality

$$\underset{j=1}{\overset{M}{}}\left|\frac{u_j}{k_j^{}+i\mu }\right|^21,$$
(3)
which is valid for any (real) $`\mu `$. This key result makes it possible to estimate the limitations on the total intensity. In particular, we conclude that the peak intensity is limited from above, i.e. $`I_{max}(z)\mathrm{max}_xI(x,z)I_{inf}`$, where the upper boundary is

$$I_{inf}=\underset{j}{\mathrm{max}}[r_j^2+(\mu _j\mu )^2],$$
(4)
and $`\mu `$ is chosen to minimize the value of $`I_{inf}`$. We will now consider several examples illustrating the physical meaning of Eqs. (3) and (4).

Stationary MSC. Let us first analyse the properties of a single MSC, composed of several solitons with identical velocities, $`\mu _n`$. In such a case, the optimal choice of the free parameter in Eq. (4) is $`\mu =\mu _n`$, and we have $`I_{inf}=\mathrm{max}_j(r_j^2)`$. Note that this value can be interpreted as the minimal squared radius of a circle which has its center at the point $`(0,\mu )`$ in the parameter space $`(r,\mu )`$ and which contains all the soliton eigenvalues within it \[see Fig. 1(a)\]. Thus $`I_{inf}`$ is proportional to the area ($`\pi R^2/2`$) of this semi-circle. On the other hand, we note that for a stationary MSC the amplitude profiles satisfy a self-consistent eigenvalue problem, $`d^2U_m/dx^2r_m^2U_m+2IU_m=0`$, where the functions $`U_m=u_m\mathrm{exp}(i\gamma _m)`$ are real. Then, since solutions should be localized in the transverse direction ($`x`$), we conclude that $`I_{max}>I_{sup}=\mathrm{max}_j(r_j^2/2)`$. This means that the variations of the peak intensity are strictly limited by the largest eigenvalue in an MSC, as illustrated in Fig. 1. Note that $`I_{max}=I_{inf}`$ if there is only one fundamental soliton (Fig. 1, left). The peak intensity decreases if several solitons compose an MSC, but always remains above the lower limit (Fig. 1, right). This happens despite the fact that the individual fundamental soliton intensities are always superimposed, i.e. destructive interference is not possible, in contrast to the case of coherent interactions. The observed decrease of total intensity underlines the complicated nature of the nonlinear superposition phenomenon and occurs because the profiles of individual solitons are strongly distorted due to the nonlinear self-action effect.

Interacting MSCs. Our general results can also be applied to the case when solitons have different velocities, $`\mu _n`$. In other words, we can estimate the peak intensity changes during the collision of several MSCs. As follows from the form of Eq. (4), the limiting value $`I_{inf}`$ depends only on the maximum eigenvalue in each of the colliding MSCs. Again, Eq. (4) has a clear geometrical interpretation: the optimal value of $`\mu `$ in Eq. (4) must be chosen to minimize the area of a half-circle which has its center at the point $`(0,\mu )`$ and which contains all the soliton eigenvalues. Two examples, corresponding to a collision between a single fundamental soliton and an MSC with different velocities ($`\mu _n`$), are shown in Fig. 2. When the relative velocity is small (Fig. 2, left), the intensity at the impact area of the collision decreases. We have already observed this effect for a stationary MSC when the relative velocity is zero. However, for larger relative velocities, the peak intensity increases (Fig. 2, right).

In order to understand the differences in the interaction pattern, we study the dependence of the peak intensities on the relative velocities of the colliding MSCs. To illustrate the key features, we consider interactions of two identical MSCs with a relative velocity $`2\mu _1`$ \[see Fig. 3(a)\]. Although Eq. (4) can be used to obtain an estimate for the peak intensities, as in the previous examples, we find that the results are not optimal for large $`\mu _1`$. To obtain a more accurate estimate, we recall that Eq. (3) is satisfied for all (real) $`\mu `$ simultaneously. We now choose $`\mu =\pm \mu _1`$, add up the corresponding inequalities together, and obtain the following upper limit, $`I_{inf}=\mathrm{max}_j[r_j^2(r_j^2+4\mu _1^2)/(r_j^2+2\mu _1^2)]`$. We thus see that $`I_{inf}(\mu _1+\mathrm{})2\mathrm{max}_j(r_j^2)`$, which is a simple sum of the upper bounds for individual MSCs. This result means that interaction of solitons with large relative velocities is weak, and their intensities are added together similarly to the linear case. This conclusion is confirmed by numerical simulations, the results of which are presented in Fig. 3(b). We also illustrate the evolution of the intensity profiles at small and large relative velocities in Figs. 3(c) and (d), respectively

Finally, we note that these results can be applied to bright solitons existing on top of a background , and also to dark solitons which can exist in media with a self-defocusing Kerr-type nonlinearity . Indeed, it has recently been demonstrated  that the total intensity profile is uniquely determined by the eigenvalues of the fundamental solitons, and does not depend on the angular distribution of radiation modes which form the background. Therefore, in these cases, our estimates are valid for the maximum changes of intensity relative to the background level $`I_b`$, i.e. $`\mathrm{max}_x|I(x,z)I_b|I_{inf}`$.

In conclusion, we have obtained an accurate estimate for the peak intensities of multi-soliton complexes for a homogeneous medium with a Kerr-type nonlinearity in the (1+1)–dimensional situation. We have demonstrated that incoherent coupling can result in a decrease of peak intensities when the relative soliton velocities are small (or zero) and nonlinear interaction is strong. On the other hand, solitons with large relative velocities, where the interaction is weak, can roughly be superimposed as linear waves.

The authors are members of the Australian Photonics Co-operative Research Centre (APCRC). We are grateful to Adrian Ankiewicz for a critical reading of this manuscript.


text:
# New results on the temporal structure of GRBs

## Introduction

The variability of GRBs provided the main evidence for the internal-external shocks scenario. External shocks cannot produce efficiently such variability 1 . Internal shocks can produce such temporal structure provided that there are two time scale within the “inner engine” - a short time scale that produces the variability and a long time scale that determines the duration of the burst. So far variability was shown only for long bursts. It is an open question whether short bursts arise in internal shocks as well. Using a new algorithm 2  we study their variability. We also present some new results on the temporal structure of long bursts. Our results provide further support for the internal shocks scenario and show that three different time scales operate within the “inner engine”.

## Variability of short bursts

We analyze the distribution of $`\delta t_{min}/T`$ (where $`\delta t_{min}`$ is the duration of the shortest pulse in a burst, and $`T`$ is the total duration of the burst) in a sample of the brightest 33 short bursts (peak flux in 64ms$`>`$4.37$`ph/(seccm^2)`$) with a good TTE data coverage (for BATSE data types review 3 ). The TTE data is binned into 2msec time bins. We compare the results to a sample of 34 long bursts with the same peak flux, using the 64ms concatenate data to which we have added a Poisson noise so that the signal to noise ratio of both samples would be similar. We call this later sample the *’noisy long’* set.

Fig 1 depicts $`\delta t_{min}/T`$ in both data sets: *’short’* and *’noisy long’*. In the *’short’* set the median of $`\delta t_{min}/T`$ is 0.25. 35% of bursts have $`\delta t_{min}/T<0.1`$ and 35% of the bursts show a smooth structure ($`\delta t_{min}/T=1`$). This result could mislead us to the conclusion that a significant fraction of the short bursts have a smooth time profile. But a comparison with the *’noisy long’* results show that also in this group more than 20% of the bursts are single pulsed, while there were no such bursts in the original (without the added noise) *’long’* set.

We conclude that short bursts are variable and hence are most likely produced in internal shocks. While the observed variability is not as large as seen in long bursts one has to remember that when studying variability of short bursts we are approaching the instrumental limitations both in terms of the time scales and of the signal to noise ratio. It is possible that 10%-20% of the short bursts are produced by external shocks.

## The pulses’ width and the intervals between pulses

According to the internal shocks model 4  the source ejects relativistic shells with different velocities and shocks arise when faster shells catch slower ones. We show in 5  that both the pulses‘ width $`\delta t`$ and the intervals between pulses $`\mathrm{\Delta }t`$ are proportional to the same parameter - the separation between two following shells, namely the variability time scale of the “inner engine”. Therefore both distributions should be similar. Moreover, any interval should be correlated to the width of its neighboring pulses.

We have applied our algorithm to a sample of the 68 brightest long bursts in BATSE 4B catalog (peak flux in 64ms$`>`$10.19$`ph/(seccm^2)`$) . This resulted in 1330 pulses (1262 intervals). Our null hypothesis was that both $`\delta t`$ and $`\mathrm{\Delta }t`$, have lognormal distributions. The $`\chi ^2`$ test gives a probability of $`0.52`$ that the pulses width were taken from a lognormal distribution with $`\mu =0.065\pm 0.04`$ ($`\overline{\delta t}1sec`$) and $`\sigma =0.77\pm 0.03`$ ($`1\sigma `$ corresponds to $`\delta t`$ between 0.5 and 2.3sec).

The $`\mathrm{\Delta }t`$ distribution shows, however, an excess of long intervals relative to a lognormal distribution. The $`\chi ^2`$ probability for a lognormal distribution is $`1.210^{10}`$. McBreen6 ; 7  and Li & Fenimore8  suggest that this deviation is due to the limited resolution (64ms). However, fitting the intervals above the median with a half Gaussian fails. The inconsistency is not due to the resolution.

Many of the long intervals are dominated by a quiescent time: periods within the burst with no observable counts above the background noise. When excluding *all* the intervals that contained a quiescent time the $`\chi ^2`$ probability that the data is lognormal is $`0.27`$, with $`\mu =0.257\pm 0.051`$ ($`\overline{\mathrm{\Delta }t}1.3sec`$) and $`\sigma =0.90\pm 0.04`$ ($`1\sigma `$ corresponds to $`\mathrm{\Delta }t`$ between 0.53 and 3.1sec).

The similarity between the parameters of both distributions is remarkable. Moreover, we find, as predicted by the internal shocks model, a linear correlation, $`r`$, between intervals and the following pulses. The average $`r`$ is 0.48, showing a strong correlation.

## Conclusions

For most short bursts $`\delta t_{min}/T1`$. This suggests that these bursts are produced by internal shocks. If, later, the ejecta encounters a surrounding ISM then we expect it to produce an external shock and emit an afterglow. For some (30% of our sample) short bursts $`\delta t_{min}T`$. However, a comparison with the *’noisy long’* set, shows that this feature could very well be due to the noise. We cannot rule out the possibility that 10%-20% of the short bursts are produced by external shocks or by a single internal collision.

The distribution of interval between pulses shows an excess of long intervals relative to a lognormal distribution. After removing intervals that include quiescent times the distribution is consistent with a lognormal distribution with comparable parameters to the pulse width distribution. This result suggests that the $`\mathrm{\Delta }t`$ distribution is made from the sum of two different distributions: A lognormal distribution that is also compatible with the $`\delta t`$ distribution and the quiescent times distribution. As $`\mathrm{\Delta }t`$ reflects the central engine behavior, this suggests that there are two different mechanisms operating within the source. A short time scale mechanism, with a lognormal distribution and a longer time scale mechanism that turns the central engine on and off and is responsible for the quiescent times. Results that support this suggestion were obtained by Ramirez & Merloni 9  . The correlation between the interval and the following pulse confirms this suggestion and is in an excellent agreement with the internal shocks model.


text:
# Transition temperature of a dilute homogeneous imperfect Bose gas

\[

## Abstract

The leading-order effect of interactions on a homogeneous Bose gas is theoretically predicted to shift the critical temperature by an amount $`\mathrm{\Delta }T_\mathrm{c}\mathrm{\#}a_{\mathrm{sc}}n^{1/3}T_0`$ from the ideal gas result $`T_0`$, where $`a_{\mathrm{sc}}`$ is the scattering length and $`n`$ is the density. There have been several different theoretical estimates for the numerical coefficient $`\mathrm{\#}`$. We claim to settle the issue by measuring the numerical coefficient in a lattice simulation of O(2) $`\varphi ^4`$ field theory in three dimensions—an effective theory which, as observed previously in the literature, can be systematically matched to the dilute Bose gas problem to reproduce non-universal quantities such as the critical temperature. We find $`\mathrm{\#}=1.32\pm 0.02`$.

preprint: UW/PT-01-06
\]

The phase transition temperature $`T_0`$ of an ideal three-dimensional Bose-Einstein gas at fixed density is something that every physicist learns to calculate in graduate school, if not before. It is amusing that the first correction to that result, from arbitrarily weak interactions, is sufficiently challenging that there has not yet been theoretical agreement on its magnitude. It is understood that, in the weak interaction limit, the correction $`\mathrm{\Delta }T_\mathrm{c}T_\mathrm{c}T_0`$ behaves parametrically as

$$\frac{\mathrm{\Delta }T_\mathrm{c}}{T_0}ca_{\mathrm{sc}}n^{1/3},$$
(1)
where $`c`$ is a numerical coefficient. (A clean argument may be found in Ref. .) Here $`a_{\mathrm{sc}}`$ is the scattering length, and the weak interaction (or dilute) limit is $`a_{\mathrm{sc}}n^{1/3}1`$: that is, $`a_{\mathrm{sc}}`$ is small compared to the typical separation between particles. (We will assume that the interaction is repulsive.) However, there has been little agreement in estimates of the coefficient $`c`$, a variety of which  are shown in Fig. 1. Some of these estimates are advertised as rough, but others are not. The difficulty arises because the phase transition is second order, and perturbation theory typically breaks down at second-order phase transitions: the physics that determines $`\mathrm{\Delta }T_\mathrm{c}`$ is non-perturbative. As we shall briefly review, the problem of finding $`\mathrm{\Delta }T_\mathrm{c}`$ in the weak interaction limit can be related to solving static three-dimensional O(2) scalar $`\varphi ^4`$ field theory . In this paper, we present results from using standard, numerical, lattice Monte Carlo methods to solve that theory. In principle this provides an exact method for computing $`c`$ to any desired precision; in practice, one is limited by computer time and memory. Working on desktop computers, we find $`c=1.32\pm 0.02`$, which is the grey bar in Fig. 1.

It is long distance physics that determines $`\mathrm{\Delta }T_\mathrm{c}`$ in the weak interaction limit. It is well known  that, at distance scales large compared to the scattering length $`a_{\mathrm{sc}}`$, an appropriate effective theory for a dilute Bose gas is the second-quantized Schrödinger equation, together with a chemical potential $`\mu `$ that couples to particle number density $`\psi ^{}\psi `$, and a $`|\psi |^4`$ contact interaction that reproduces low-energy scattering. The corresponding Lagrangian is

$$=\psi ^{}\left(i\mathrm{}_t+\frac{\mathrm{}^2}{2m}^2+\mu \right)\psi \frac{2\pi \mathrm{}^2a_{\mathrm{sc}}}{m}(\psi ^{}\psi )^2.$$
(2)
Corrections to this effective theory (due, for instance, to energy dependence of the cross-section or 3-body interactions) may be ignored for the propose of computing the leading-order result for $`\mathrm{\Delta }T_\mathrm{c}`$. To study (2) at finite temperature, apply the imaginary time formalism, so that $`t`$ becomes $`i\tau `$ and imaginary time $`\tau `$ is periodic with period $`\mathrm{}\beta =\mathrm{}/k_\mathrm{B}T`$. The field $`\psi `$ can then be decomposed into frequency modes with Matsubara frequencies $`\omega _n=2\pi n/\mathrm{}\beta `$. For distances large compared to the thermal wavelength $`\lambda =\mathrm{}\sqrt{2\pi \beta /m}`$, and sufficiently near the transition so that $`|\mu |T`$, the non-zero Matsubara frequencies decouple from the dynamics, leaving behind an effective theory of only the zero-frequency modes $`\psi _0`$, with the action $`S=\mathrm{}^1_0^\mathrm{}\beta 𝑑\tau d^3x`$ becoming

$$\beta d^3x\left[\psi _0^{}\left(\frac{\mathrm{}^2}{2m}^2\mu _{\mathrm{eff}}\right)\psi _0+\frac{2\pi \mathrm{}^2a_{\mathrm{sc}}}{m}(\psi _0^{}\psi _0)^2\right],$$
(3)
up to corrections that again do not affect the leading-order result for $`\mathrm{\Delta }T_\mathrm{c}`$. This action can be thought of as the $`\beta H`$ of a classical three-dimensional field theory. Finally, it is convenient to rewrite $`\psi _0=(\varphi _1+i\varphi _2)\sqrt{2\pi }/\lambda `$ so that the effective action becomes a conventionally normalized O(2) field theory:

$$S=d^3x\left[\frac{1}{2}|\mathbf{}\varphi |^2+\frac{1}{2}r\varphi ^2+\frac{u}{4!}(\varphi ^2)^2\right],$$
(4)
where $`\varphi `$ is understood to be a 2-component real vector $`(\varphi _1,\varphi _2)`$ and $`u=96\pi ^2a_{\mathrm{sc}}/\lambda ^2`$.

As noted by Baym et al., it is technically somewhat more convenient, in this formalism, to calculate the shift $`\mathrm{\Delta }n_\mathrm{c}(T)`$ in the critical density at fixed temperature instead of the shift $`\mathrm{\Delta }T_\mathrm{c}(n)`$ in the critical temperature at fixed density. The two are trivially related at first order in $`\mathrm{\Delta }T_\mathrm{c}`$ by $`\mathrm{\Delta }T_\mathrm{c}/T=\frac{2}{3}\mathrm{\Delta }n_\mathrm{c}/n`$, where the factor of $`\frac{2}{3}`$ arises from the ideal gas relation $`T_0n^{2/3}`$. In the field theory, $`n`$ is given by $`\psi ^{}\psi `$, which is proportional to $`\varphi ^2`$. Putting everything together ,

$$\frac{\mathrm{\Delta }T_\mathrm{c}}{T_0}=\frac{2mk_\mathrm{B}T_0}{3\mathrm{}^2n}\mathrm{\Delta }\varphi ^2_\mathrm{c},$$
(5)
where

$$\mathrm{\Delta }\varphi ^2_\mathrm{c}\left[\varphi ^2_\mathrm{c}\right]_u\left[\varphi ^2_\mathrm{c}\right]_0$$
(6)
is the difference between the effective theory value of $`\varphi ^2`$, at the critical point, for the cases of (i) $`u`$ small and (ii) the ideal gas $`u=0`$. Unlike $`\varphi ^2_\mathrm{c}`$, the difference $`\mathrm{\Delta }\varphi ^2_\mathrm{c}`$ is an infrared quantity, independent of how the effective theory (4) is regularized in the ultraviolet (UV).

Finding $`\left[\varphi ^2_\mathrm{c}\right]_u`$ in the effective theory (4) corresponds to fixing $`u`$, varying $`r`$ to reach the critical point, and then measuring $`\varphi ^2`$. The only parameter of this problem is $`u`$. By dimensional analysis, physics is therefore non-perturbative at the infrared length scale $`1/u`$, and, again by dimensional analysis, $`\mathrm{\Delta }\varphi ^2_\mathrm{c}`$ is proportional to $`u`$. Putting this together with the ideal gas formula $`T_0=(2\pi \mathrm{}^2/k_\mathrm{B}m)[n/\zeta \left(\frac{3}{2}\right)]^{2/3}`$, one may summarize the relationship between the weak interaction limit for $`\mathrm{\Delta }T_\mathrm{c}`$ and the O(2) effective theory as $`\mathrm{\Delta }T_\mathrm{c}/T_0ca_{\mathrm{sc}}n^{1/3}`$ with

$$c=\frac{128\pi ^3}{\left[\zeta (\frac{3}{2})\right]^{4/3}}\frac{\mathrm{\Delta }\varphi ^2_\mathrm{c}}{u}.$$
(7)
To compute $`\mathrm{\Delta }\varphi ^2_\mathrm{c}/u`$, we put the O(2) theory (4) on a lattice. For example, the most straightforward discretization would use the action

$`S_\mathrm{U}=a_{\mathrm{lat}}^3{\displaystyle \underset{\mathrm{𝐱𝐲}}{}}{\displaystyle \frac{1}{2}}\left[\left({\displaystyle \frac{\varphi _{1𝐱}\varphi _{1𝐲}}{a_{\mathrm{lat}}}}\right)^2+\left({\displaystyle \frac{\varphi _{2𝐱}\varphi _{2𝐲}}{a_{\mathrm{lat}}}}\right)^2\right]`$ (8)
$`+a_{\mathrm{lat}}^3{\displaystyle \underset{𝐱}{}}\left[{\displaystyle \frac{r_{\mathrm{lat}}}{2}}(\varphi _{1𝐱}^2+\varphi _{2𝐱}^2)+{\displaystyle \frac{u}{4!}}(\varphi _{1𝐱}^2+\varphi _{2𝐱}^2)^2\right],`$ (9)
on a simple cubic lattice, where $`\mathrm{𝐱𝐲}`$ represents all nearest-neighbor pairs and $`a_{\mathrm{lat}}`$ is the lattice spacing (unrelated to $`a_{\mathrm{sc}}`$). The dimensionless coupling of the lattice theory is $`ua_{\mathrm{lat}}`$, and the continuum limit is $`ua_{\mathrm{lat}}0`$.

Our simulations use an improved action to reduce lattice spacing errors. (The subscript U on $`S`$ in Eq. (9) stands for “unimproved.”) Details concerning the action and our simulations are given in Ref. . Our simulations use a combination of heat bath and multi-grid updates . At finite volume, we use the method of Binder cumulants  to determine a nominal critical value for $`r`$.

Because three-dimensional scalar theory requires UV regularization of its $`\varphi ^2`$ interactions, the $`r_{\mathrm{lat}}`$ in the lattice action (9) is not simply the $`r_{\mathrm{eff}}`$ in the effective theory action (4), which in turn is not simply related to the chemical potential $`\mu `$ in the original action (2). However, for our purposes here, we are only interested in adjusting $`r_{\mathrm{lat}}`$ to find the critical point, for a given $`u`$, and measuring $`\mathrm{\Delta }\varphi ^2`$ there. The actual relation between $`r_{\mathrm{lat}}`$, $`r_{\mathrm{eff}}`$, and $`\mu `$ is unnecessary . The $`\varphi ^4`$ interactions, in contrast, do not require UV regularization: in the limit $`ua_{\mathrm{lat}}0`$, the coupling $`u`$ of the lattice theory may be identified with the continuum coupling $`u=96\pi ^2a_{\mathrm{sc}}/\lambda ^2`$ introduced earlier.

For a given $`u`$, we compute $`\varphi ^2_\mathrm{c}`$ by Monte Carlo numerical simulations of the lattice theory. The $`u=0`$ piece of the difference $`\mathrm{\Delta }\varphi ^2_\mathrm{c}`$ can be easily computed without simulations:

$$\mathrm{\Delta }\varphi ^2_\mathrm{c}=\underset{a_{\mathrm{lat}}0}{lim}\left[\varphi ^2_{\mathrm{lat}}_{𝐩\mathrm{BZ}}G_{\mathrm{lat}}^{ii}(𝐩)\right],$$
(10)
where $`G_{\mathrm{lat}}^{ij}(𝐩)`$ is the free lattice propagator, and the momentum $`𝐩`$ is integrated over the Brillouin zone (BZ). Such integrals are reviewed in Ref. . For the unimproved lattice theory (9), for example ,

$$\mathrm{\Delta }\varphi ^2_\mathrm{c}=\underset{a_{\mathrm{lat}}0}{lim}\left[\varphi ^2_{\mathrm{lat}}\frac{\mathrm{\Sigma }_\mathrm{U}}{2\pi a_{\mathrm{lat}}}\right],$$
(11)
where $`\mathrm{\Sigma }_\mathrm{U}=3.175911535625\mathrm{}`$.

There are two limits that must be taken of lattice Monte Carlo data: the continuum limit $`ua_{\mathrm{lat}}0`$ and the infinite volume limit $`Lu\mathrm{}`$. Fig. 2 shows the dependence of our data on system size ($`Lu`$) at $`ua_{\mathrm{lat}}=6`$. Fig. 3 shows the dependence on lattice spacing ($`ua_{\mathrm{lat}}`$) at $`Lu=576`$. From these two figures, it is reasonably clear that our raw data includes reasonably large volumes and reasonably small lattice spacings. We will discuss extrapolations of the infinite-volume continuum limit using finite volume scaling, the known critical exponents of this model, and an analysis of finite lattice-spacing errors.

First, however, we wish to show that one can obtain a quick estimate from the raw data without relying on anything fancy. Table I shows values associated with a selected subset of the largest $`Lu`$ and smallest $`ua_{\mathrm{lat}}`$ data points from the figures. We take the $`(Lu,ua_{\mathrm{lat}})=(576,6)`$ result as a starting point for our estimate. The finite-volume correction is at least as big as the difference with the value at $`(1152,6)`$ \[which corresponds to doubling $`Lu`$\] but is unlikely to be double this difference. This difference is roughly $`0.00015`$ (ignoring the small statistical errors), so we might estimate the finite-volume correction to the $`Lu=576`$ value to be somewhere between $`0.00015`$ and $`0.00030`$. From a similar comparison of the (576,6) and (576,3) data, we might estimate the finite lattice spacing correction to (576,6) to be between roughly $`+0.00005`$ and $`+0.00011`$. Adding our corrections and the original $`(576,6)`$ data point, the final rough estimate of the continuum infinite-volume value would be

$$\frac{\mathrm{\Delta }\varphi ^2_\mathrm{c}}{u}0.00119\pm 0.00011,$$
(12)
which, by (7) would translate to $`c=1.31\pm 0.12`$.

We now summarize a more careful analysis of corrections, detailed in ref. . Our strategy is to start again with $`Lu=576`$ data, extrapolate the continuum limit, and estimate the finite volume correction. To improve the approach to the continuum limit, we have analytically calculated the $`O(ua_{\mathrm{lat}})`$ corrections to $`\mathrm{\Delta }\varphi ^2_\mathrm{c}`$ and the relation between lattice and continuum values of $`u`$, using lattice perturbation theory . Fig. 4 shows $`Lu=144`$ data which clearly demonstrates our control of lattice spacing errors. The uncorrected data clearly has a linear dependence on $`ua_{\mathrm{lat}}`$. But the corrected data fits, to high confidence level, the assertion that the remaining error scales as $`(ua_{\mathrm{lat}})^2`$. Based on a similar fit to the data of Fig. 3, we estimate the $`ua_{\mathrm{lat}}0`$ result at $`Lu=576`$ as

$$\left[\frac{\mathrm{\Delta }\varphi ^2_\mathrm{c}}{u}\right]_{Lu=576}=0.000957\pm 0.000015.$$
(13)
Finite scaling arguments predict that the large $`Lu`$ corrections to $`\mathrm{\Delta }\varphi ^2_\mathrm{c}`$ should scale as

$$\mathrm{\Delta }\varphi ^2_\mathrm{c}A+BL^{(1\alpha )/\nu }$$
(14)
if the method of Binder cumulants is used to determine the transition point in finite volume. Here $`\alpha 0.01`$ and $`\nu =(2\alpha )/30.67`$ are the specific heat and correlation length critical exponents of the O(2) model . Further discussions of fits, and an analysis of corrections to scaling, may be found in ref. . On the basis of these fits, we estimate the finite size correction at $`Lu=576`$ to be $`0.000241\pm 0.000007`$ . Putting this together with the $`Lu=576`$ continuum extrapolation (13), we obtain the final infinite-volume continuum result

$$\frac{\mathrm{\Delta }\varphi ^2_\mathrm{c}}{u}=0.001198\pm 0.000017.$$
(15)
Using Eq. (7) for the weak-interaction limit of $`\mathrm{\Delta }T_\mathrm{c}`$, the result for $`\mathrm{\Delta }\varphi ^2_\mathrm{c}`$ translates to

$$\frac{\mathrm{\Delta }T_\mathrm{c}}{T_0}(1.32\pm 0.02)a_{\mathrm{sc}}n^{1/3}.$$
(16)
It is interesting to compare our numerical results with results from the large $`N`$ expansion, depicted by triangles in Fig. 1. In large $`N`$, one generalizes the O(2) effective field theory to an O($`N`$) theory of $`N`$ real fields, calculates results in powers of $`1/N`$, and hopes the expansion will be useful for the case of interest, $`N=2`$. The leading-order (LO) result $`c2.33`$  is off by 77%, but the next-to-leading-order (NLO) result $`c1.71`$  moves in the right direction and is off by only 30%. This is surprisingly good for an expansion that treats $`N=2`$ as large.

We should comment on the discrepancy of our results with previous numerical calculations in the literature , shown in Fig. 1, which used a radically different starting point. Rather than using field theory methods and the grand canonical ensemble, they start with the path integral for a large, fixed number of particles in a box. The recent work of Holzmann and Krauth , however, makes a flawed assumption at the very beginning: they assume that the integrand of the path integral can be expanded perturbatively in the interaction, and keep only the leading term. This is wrong because interactions cannot, generically, be treated perturbatively at a second-order phase transition.

We believe one likely problem with the older simulations of Grüter et al. is inadequate system size . Reppy et al. have reported $`c=5.1\pm 0.9`$ from experimental data on He-Vycor systems, but cautions about the data’s interpretation may be found in Ref. .

As we completed this work, another paper appeared  which uses techniques very similar to ours and obtains the statistically compatible result $`c=1.29\pm 0.05`$.

This work was supported by the U.S. Department of Energy under Grant Nos. DE-FG03-96ER40956 and DE-FG02-97ER41027.


text:
# The binary fraction of extreme horizontal branch stars.

## 1 Introduction

Surveys for blue stars brighter than $`B16`$ are dominated by subdwarf-B (sdB) stars (Green, Schmidt & Liebert 1986). The effective temperatures (Teff) and surface gravities ($`\mathrm{log}g`$) of the majority of these stars place them on the extreme horizontal branch (EHB), i.e., they appear in the same region of the Teff – $`\mathrm{log}g`$ plane as evolutionary tracks for core helium burning stars with core masses of about 0.5$`\mathrm{M}_{}`$ and extremely thin ($`0.02\mathrm{M}_{}`$) hydrogen envelopes (Heber 1986; Saffer et al. 1994). We make a distinction in this paper between the nomenclature “sdB star”, which is a spectral classification, and “EHB star” which is an interpretation of the evolutionary state of a star.

The observed dispersion of core masses for EHB stars is very low ($`<0.04\mathrm{M}_{}`$, Saffer et al. 1994). It is thought that the eventual fate of an EHB star is to cool to form a white dwarf with a mass of about 0.5$`\mathrm{M}_{}`$, which is low compared to the typical mass for white dwarfs (Bergeron, Saffer & Liebert 1992). The formation of low mass white dwarfs is, in general, thought to involve interactions with a binary companion star, e.g., a common envelope phase, in which a companion to a red giant star is engulfed by the expanding outer layers. The resulting friction causes the companion to spiral in towards the core of the red giant, ejecting the envelope at the expense of orbital binding energy (Iben & Livio 1993). If this process occurs while the red giant is within $``$0.4 magnitudes of the tip of the red giant branch, the core can go on to ignite helium, despite the dramatic mass loss, and may then appear as an EHB star (D’Cruz et al. 1996; Mengel, Norris & Gross 1976).

The binary fraction of sdB and EHB stars is expected to be high given the scenario outlined above. Allard et al. (1994) found that 31 of their sample of 100 sdB stars show flat spectral energy distributions which indicate the presence of companions with spectral types in the range late-G to early-M. They infer a binary fraction for main-sequence companions of 54 – 66 percent, although the companions in their survey appear to be over-luminous compared to normal main-sequence stars. A similar conclusion was reached by Ferguson, Green & Liebert (1984) using a similar argument and by Jeffery & Pollacco (1998) based on the detection of spectral features due to cool companions. What is not clear from these observations is whether the cool companion is sufficiently close to the EHB star to be implicated in the mass loss process that is supposed to form the EHB star. These techniques are also insensitive to white dwarf companions and faint M-dwarf companions. Companions to EHB stars can also be detected in eclipsing systems such as the short-period EHB – M-dwarf binaries HW Vir (Wood & Saffer 1999) and PG 1336$``$018 (Kilkenny et al. 1998) and in EHB – white dwarf binaries which show ellipsoidal variability, e.g., KPD 0422+5421 (Koen et al. 1998) and KPD 1930+2752 (Maxted, Marsh & North 2000). These binaries are extremely useful for studying the properties of EHB stars, but they do not offer a useful method for finding binary EHB stars in general because the probability of such a binary showing eclipses or a measureable ellipsoidal effect decreases rapidly for increasing orbital periods.

Radial velocity surveys are an excellent method for identifying binary stars in general, particularly since the effeciency of this technique can be accurately quantified and the selection effects are well understood. This technique can be applied to many types of star, from main-sequence stars (Duquennoy & Mayor 1991) to white dwarfs (Maxted & Marsh 1999). For EHB stars in particular, this is the method of choice because the short-period binaries which are expected to result from a common-envelope phase are the easiest to identify using this method. If the radial velocities are measured to an precision of a few km s-1 and the observations are obtained over a baseline of weeks or months the technique has the potential to identify binaries with much longer periods ($``$ 100 d) even if the companion is a low mass M-dwarf. Saffer, Livio & Yungelson (1998) have shown the potential for this technique with their observations of 46 sdB stars. The precision of their radial velocity measurements was modest (20–30 km s-1) and the three spectra they obtained for each star over a baseline of 1–2 days were compared by-eye, yet they found that at least 7 of their sample of 46 sdB stars show radial velocity variations. Several of these binaries have subsequently had their orbital periods determined (Moran et al. 1999, Maxted et al. 2000), although further observations are required to determine the nature of the companions in these binaries.

In this paper we present the results of a radial velocity survey of binary EHB stars. We have used observations of the H$`\alpha `$ line to measure 205 precise radial velocities for 36 EHB stars from the Palomar-Green survey (Green, Schmidt & Liebert 1986). We positively identify 22 short-period binary EHB stars, 20 of which are new discoveries. We conclude that at least 60$`\pm `$8 percent of of EHB stars are short-period binary stars. If we allow for the detection efficiency of our survey, we find that at least 69$`\pm `$9 percent of EHB stars are binaries. This is strong evidence that binary star evolution is fundamental to the formation of the majority of EHB stars. We also observed 5 stars which we identify as post-EHB stars and found that one of these stars is a binary.

## 2 Observations and reductions

### 2.1 H$`\alpha `$ spectra

Targets were selected from objects in the Palomar-Green catalog (Green, Schmidt & Liebert 1986) classified as sdB stars. We avoided stars where follow-up observations have shown the classification was in error or that the star is not an EHB star. Observations were obtained with the 2.5m Isaac Newton Telescope on the Island of La Palma. Spectra were obtained with the intermediate dispersion spectrograph using the 500mm camera, a 1200 line/mm grating and a TEK charge coupled device (CCD) as a detector. The spectra cover 400Å around the H$`\alpha `$ line at a dispersion of 0.39Å per pixel. The slit width used was 0.97 arcsec which gave a resolution of about 0.9Å. Spectra of the targets were generally obtained in pairs bracketed by observations of a copper-neon arc. We obtained a total of 243 spectra for 43 stars over a total of about 7 nights during the interval 2000 April 10 – 21. The seeing was good ($``$ 1 arcsec) on most of these nights.

Extraction of the spectra from the images was performed automatically using optimal extraction to maximize the signal-to-noise of the resulting spectra (Marsh 1989). The arcs associated with each stellar spectrum were extracted using the same weighting determined for the stellar image to avoid possible systematic errors due to the tilt of the spectra on the detector. The wavelength scale was determined from a fourth-order polynomial fit to measured arc line positions. The standard deviation of the fit to the 8 arc lines was typically 0.09Å. The wavelength scale for an individual spectrum was determined by interpolation to the time of mid-exposure from the fits to arcs taken before and after the spectrum to account for the small amount of drift in the wavelength scale ($`<0.1`$Å) due to flexure of the instrument. Statistical errors on every data point calculated from photon statistics are rigorously propagated through every stage of the data reduction.

### 2.2 Blue spectra

In order to measure the effective temperature and surface gravity of some our targets we also obtained blue spectra of our targets with the same telescope and instrument. We did not attempt to measure radial velocities from these spectra.

Spectra of PG 1032+406, PG 1043+760, PG 1051+501, PG 1039+219, PG 1043+760 and PG 1110+294 were obtained over the wavelength range 3810–5020Å using a 400 line/mm grating. The observations were obtained while the stars were at low airmass in good seeing with a vertical 1.5 arcsec slit. The resolution in pixels was determined from the width of the spatial profile, which was typically 3 – 4 pixels which corresponds to a resolution of about 4Å.

We used an EEV CCD on the 235mm camera and a 900 line/mm grating to obtain spectra of PG 1505+074, PG 1512+244 and PG 1553+273 on the night 2000 July 16 and of PG 1616+144, PG 1627+017, PG 1632+088, PG 1647+056, and PG 1653+131 on the night 2000 August 15. The useful region of the spectra cover the wavelength range 3590 – 5365Å at a dispersion of 0.63Å per pixel. We used a vertical, 1 arcsec wide slit which gave a resolution of 1.6Å. We also obtained spectra of PG 0907+123 and PG 1116+301 on the night of 2001 February 3 with the same instrument covering the wavelength range 3850 – 5200Å. One other spectrum of PG 0907+123 was also obtained and reduced for us by Martin Altman using the Calar Alto 2.2m telescope and the CAFOS spectrograph with a B100 grism at a dispersion of 100Å/MM at lower spectral resolution than our INT spectra, but covering the Balmer lines from H$`\beta `$ to H10.

## 3 Analysis

### 3.1  Effective temperatures, surface gravities and helium abundances.

For those stars for which we have blue spectra we measured the effective temperature, Teff, the surface gravity $`\mathrm{log}g`$ and the helium abundance by number, $`y`$. The simultaneous fitting of Balmer line profiles by a grid of synthetic spectra has become the standard technique to determine the atmospheric parameters of hot high gravity stars (Bergeron et al. 1992). The procedure has been extended to include helium line profiles and applied successfully to sdB stars by Saffer et al. (1994). We have applied Saffer’s procedure to the Balmer lines (H$`\beta `$ to H 9), and the He I (4026Å, 4388Å, 4471Å, 4713Å, 4922Å) and He II 4686Å lines.

A grid of synthetic spectra derived from H and He line blanketed NLTE model atmospheres (Napiwotzki 1997) was matched to the data to simultaneously determine the effective temperature, surface gravity and helium abundance. For stars cooler than 27 000 K we used the metal line-blanketed LTE model atmospheres of Heber, Reid, & Werner (2000). The synthetic spectra were convolved beforehand with a Gaussian profile of the appropriate width to account for the instrumental profile. The adopted values of Teff, $`\mathrm{log}g`$ and $`y`$ for these stars and all other stars where values could be found are given in Table 1. Examples of the observed spectral lines and synthetic spectrum fits for six sdB stars are shown in Fig. 1. The values given for PG 1040+234 are only approximate because of the contamination by the companion star, particularly in the H$`ϵ`$ line, which we excluded from the fit. The values given for PG 1701+359, PG 1722+286 and PG 1743+477 are based on updated fits to the spectra described Theissen et al. (1993). The values for PG 0907+123 are an average of the values derived from our INT spectrum and the spectrum taken by Martin Altman, which agree very well.

The measured values of Teff and $`\mathrm{log}g`$ are compared to the evolutionary tracks for extreme horizontal branch stars of Dorman et al. (1993) in Fig. 2. It can be seen that most the targets lie in or near the band defined by the zero-age extreme horizontal branch (ZAEHB), the terminal-age extreme horizontal branch (TAEHB) and the helium main-sequence (HeMS) and are therefore EHB stars. The errors in the atmospheric parameters are estimated to be 3% for Teff and 0.1 dex for $`\mathrm{log}g`$ and the helium abundance $`y`$. We consider a star to be a post-EHB star if it lies above the TAEHB by more than these error bars, otherwise we regard it as an EHB star. The spectra of PG 1040+234 and PG 1701+359 are contaminated by the cool companion. The principal effect of this contamination is to bias the value of $`\mathrm{log}g`$ to lower values, so we treat the value obtained as an upper limit in Table 1 and in Fig. 2. Both these stars are included as EHB stars in our discussion of the binary fraction of EHB stars. PG 1632+088 is too cool to appear in Fig. 2. It is probably a normal horizontal branch star so we exclude it from our discussion of the binary fraction of EHB stars. We also exclude PG 0909+164, PG 1000+408, PG 1051+501, PG 1505+074 and PG 1553+273 because they lie too far from the EHB in the Teff – $`\mathrm{log}g`$ plane. They appear to be more evolved than EHB stars so we classify them as post-EHB stars and discuss them separately from the EHB stars. We also exclude PG 1631+267 from our discussion of the binary fraction of EHB stars because no Teff – $`\mathrm{log}g`$ measurement is available for this star.

### 3.2 Radial velocity measurements.

To measure the radial velocities we used least-squares fitting of a model line profile. This model line profile is the summation of three Gaussian profiles with different widths and depths but with a common central position which varies between spectra. Only data within 2000 km s-1 of the H$`\alpha `$ line is included in the fitting process and the spectra are normalized using a linear fit to the continuum either side of the H$`\alpha `$ line.

To measure the radial velocities we use a fitting process with four steps to determine an optimum set of radial velocities. We use a least-squares fit to one of the spectra to determine an initial shape of the model line profile. A least squares fit of this profile to each spectrum in which the position of the line is the only free parameter gives an initial set of radial velocities. We use these initial radial velocities to fix the position of the H$`\alpha `$ line in a simultaneous fit to all the spectra to obtain an improved model line profile. A least squares fit of this profile to each spectrum yields the radial velocities given in Table 4. The uncertainties quoted are calculated by propagating the uncertainties on every data point in the spectra right through the data reduction and analysis. These uncertainties are reliable in most cases, but some caution must be exercised for quoted uncertainties $`2`$ km s-1. This corresponds to about $`1/10`$ of a pixel in the original data, so systematic errors such as telluric absorption features, uncertainties in the wavelength calibration and motion of the star within the slit during good seeing are certain to be a significant source of uncertainty for these measurements. An example of the observed spectra and multiple Gaussian fits for one star is shown in Fig. 3.

We rebinned all the spectra onto a common wavelength scale allowing for the measured radial velocity shifts and then formed the average spectrum of each star shown in Fig. 5.

### 3.3 Criterion for variability.

For each star we calculate a weighted mean radial velocity. This mean is the best estimate of the radial velocity of the star assuming this quantity is constant. We then calculate the $`\chi ^2`$ statistic for this “model”, i.e., the goodness-of-fit of a constant to the observed radial velocities. We can then compare the observed value of $`\chi ^2`$ with the distribution of $`\chi ^2`$ for the appropriate number of degrees of freedom. We then calculate the probability of obtaining the observed value of $`\chi ^2`$ or higher from random fluctuations of constant value, $`p`$. To allow for the systematic errors described above, we have added a 2 km s-1 external error in quadrature to all the radial velocity uncertainties prior to calculating these statistics. If we find $`\mathrm{log}_{10}(p)<4`$ we consider this to be a detection of a binary. In our sample of 36 EHB stars, this results in a less than 0.4 percent chance of random fluctuations producing one or more false detections.

### 3.4 Results

The observed values of $`\chi ^2`$ and $`\mathrm{log}_{10}(p)`$ and the number of measured radial velocities, N, are given for all the targets in our sample in Table 2. Stars which were observed but are not EHB stars are shown in parentheses. Stars which we consider to be binaries are denoted by displaying $`\mathrm{log}_{10}(p)`$ in bold type. In column 3 we give the maximum difference between the observed radial velocities, $`\mathrm{\Delta }`$. In column 6 (SLY98) we note whether Saffer, Livio & Yungelson (1998) saw a marginal detection (2) or a positive detection (1) of a radial velocity shift or failed to detect any radial velocity shift ($`\times `$). In column 7 (JP98) we note whether Jeffery & Pollacco (1998) saw spectral features due to a cool companion (Y) or failed to detect a companion ($`\times `$). In column 8 (AWFBL98) we note whether the BVRI photometry of Allard et al. (1994) failed to detect a companion ($`\times `$) or note the spectral type of the companion if it was detected. In column 8 (UT98) we note stars for which Ulla & Thejll did not detect any infrared excess due to a companion from their JHK photometry ($`\times `$). Stars for which comments can be found in section 4 are noted in column 9. There are 36 EHB stars in our sample, 21 of which are binaries. With the exception of PG 1716+426 and, perhaps, PG 0839+399, these are all new detections.

## 4  Notes on individual objects.

In this section we note previous results for our targets and any other remarkable or peculiar characteristics.

Late-type spectral features can be seen in the average spectrum of this star shown in Fig. 5.

There are weak spectral features due to a late-type companion visible in our spectra.

This star in listed in Jeffery & Pollacco as Ton 1273.

Spectral features due to the companion are seen in our blue spectra of this star, notably the G-band and Ca II H&K spectral lines, and some weak features can also be seen around H$`\alpha `$ (Fig. 5).

This is a pulsating sdB variable star (Billères et al. 1997).

Saffer, Livio & Yungelson list this star as PG 1114+072.

This star has a bright G-type companion which dominates the spectrum around the H$`\alpha `$ line so the radial velocities quoted here refer to the G-star companion to the sdB star, which we denote PG 1631+267 B.

Theissen et al. noted spectral features from a cool star in their spectra. Spectral features from a cool star are also visible in our spectra around the H$`\alpha `$ line.

Geffert (1998) considers the galactic orbit of this star based on the HIPPARCOS astrometry and a radial velocity measurement of $`10.6\pm 30`$ km s-1. Clearly, this calculation needs to be revised.

## 5 Discussion

### 5.1 Estimating the binary fraction.

The probability of detecting $`N_B`$ binaries in a sample of $`N`$ stars which have a binary fraction of $`f`$ is

$$\frac{N!}{(NN_B)!N_B!}(\overline{d}f)^{N_B}(1\overline{d}f)^{NN_B}$$

where $`\overline{d}`$ is the fraction of all binaries detected by the survey averaged over all orbital periods. For our survey, $`N_B=21`$ and $`N=36`$. We can set a lower limit to $`f`$ by assuming $`\overline{d}=1`$, i.e., the lower limit to $`f`$ is set by assuming we have detected all the binaries in our sample. In this case we expect that the lower limit to the binary fraction will be about 21/33=58.3 percent. In fact, the distribution of $`f`$ calculated with the expression above with $`\overline{d}=1`$ is approximately Gaussian with a maximum at $`f=0.60`$ and a standard deviation of 0.08, i.e., the absolute lower limit to $`f`$ from our survey is 60$`\pm `$8 percent.

We calculated the fraction of all binaries of a given orbital period, $`P`$, detected by our survey, $`d`$, as follows. We assume that the EHB star and its companion both have a mass of 0.5$`\mathrm{M}_{}`$. We can then calculate the orbital speed of the EHB star, $`V_{\mathrm{orb}}`$, assuming a circular orbit. We assume the that orbits are circular because a common envelope phase will quickly reduce the eccentricity of an orbit and no post-common envelope systems are observed to have any appreciable eccentricity. For a given star for which we have $`N_{\mathrm{obs}}`$ radial velocity measurements we can then use the actual dates of observation, $`T_j,j=1\mathrm{}N_{\mathrm{obs}}`$ to calculate radial velocities for a hypothetical binary with an edge-on orbit from $`V_{\mathrm{orb}}\mathrm{sin}(\varphi _i)`$, where $`\varphi =(T_jT_0)/P`$. These values are used to calculate the value of chi-squared for this hypothetical binary, $`\chi _{\mathrm{max}}^2`$, using the actual radial velocity uncertainties for the $`N_{\mathrm{obs}}`$ observations given in Table 4 including 2 km s-1 additional systematic uncertainty. The calculation is repeated for 50 values of the $`T_0`$ and the average value of $`\chi _{\mathrm{max}}^2`$ is taken. We can then compare the value of $`\chi _{\mathrm{max}}^2`$ for this hypothetical binary to the value of chi-squared required to exactly satisfy our detection criterion, $`\chi _{\mathrm{crit}}^2`$. If $`\chi _{\mathrm{max}}^2<\chi _{\mathrm{crit}}^2`$ then no binaries with that orbital period, mass and eccentricity will be detected by our observations. 111This is not strictly true. Adding random fluctuations to the $`V_{\mathrm{orb}}\mathrm{sin}(\varphi _i)`$ values can result in detections in cases where $`\chi _{\mathrm{max}}^2`$ is only slightly less then $`\chi _{\mathrm{crit}}^2`$. Similarly, noise can prevent some detections when $`\chi _{\mathrm{max}}^2`$ is slightly greater than $`\chi _{\mathrm{crit}}^2`$. The overall effect is negligible when the detection effeciency is averaged over a wide range of orbital periods as we have done. Otherwise, we can calculate the projected orbital velocity for which $`\chi _{\mathrm{max}}^2=\chi _{\mathrm{crit}}^2`$, $`K_{\mathrm{crit}}=V_{\mathrm{orb}}\mathrm{sin}i`$ for some orbital inclination $`i`$. For randomly oriented orbits, $`i`$ is distributed as $`\mathrm{cos}i`$ so the fraction of binaries detected for this combination of observations, period, mass etc. is simply $`d=\sqrt{1(K_{\mathrm{crit}}/V_{\mathrm{orb}})^2}`$. We have calculated this detection efficiency for 20 000 orbital periods distributed uniformly in $`\mathrm{log}_{10}(P)`$ over the range $`1.5\mathrm{log}_{10}(P/\mathrm{d})2`$ for every EHB star we observed and used these values to calculate the average detection efficiency for stars in our sample, $`d`$. The results are shown in Fig. 4, where the value of $`d`$ has been binned into 400 groups of 50 periods.

To calculate the binary fraction of EHB stars, we need to know $`\overline{d}`$, the weighted mean of $`d`$ over the period distribution of EHB binaries. Unfortunately, the period distribution of EHB binaries is very poorly known. The existing observational data for EHB stars with measured orbital periods is rather scarce, but is summarised in Table 3 and is also shown in Fig. 4. We are not aware of any reliable predictions for the orbital period distribution based on models of the evolution of EHB stars. From the size of the radial velocity shifts given in Table 2 we can set an upper limit to the orbital period of the binaries we have found. These are typically tens-of-days, so the actual orbital periods are likely to be $`10`$ d.222Observations to determine the actual orbital periods are being undertaken at the time of writing.

In the absence of any good determination of the orbital period distribution of binary EHB stars for longer orbital periods, we consider the binary fraction for short orbital period binaries only ($`P10`$ d). We can see from Fig. 4 that any reasonable period distribution will give a mean detection efficiency over $`1.5\mathrm{log}_{10}(P/\mathrm{d})1`$ of about 85 percent. For example, the unweighted average over this range of $`\mathrm{log}_{10}(P)`$ is 86.6 percent, which implies a binary fraction of 69$`\pm `$9 percent. If, for the sake of argument, we assume a distribution for $`\mathrm{log}_{10}(P)`$ which is a Gaussian function with a mean of $`\mathrm{log}_{10}(P/\mathrm{d})=0`$ and a full-width at half-maximum of $`\mathrm{log}_{10}(P/\mathrm{d})=2`$, the mean detection efficiency over $`1.5\mathrm{log}_{10}(P/\mathrm{d})1`$ is 84.7 precent so we obtain a binary fraction of 70$`\pm `$9 percent. If we change the mean of the Gaussian function to $`\mathrm{log}_{10}(P/\mathrm{d})=1`$, the mean detection efficiency over the same range of $`\mathrm{log}_{10}(P)`$ is 91.1 percent and the binary fraction is 65$`\pm `$9 percent. These are all ad-hoc assumptions for the period distribution of binary EHB stars, but they do show that the fraction of EHB stars which are short-period binaries is about 2/3 for any reasonable period distribution.

This calculation is also insensitive to the assumed mass ratios of the binaries. If we assume the companions have a mass of 1$`\mathrm{M}_{}`$, the lower limits to the binary fraction we derive are reduced by 1–2 percent. The mean companion mass is unlikely to be larger than 1$`\mathrm{M}_{}`$ because a main-sequence or sub-giant star of this mass would be easily visible in the spectrum and the upper limit to the white dwarf companion mass is, of course, the Chandrasekhar mass of 1.4$`\mathrm{M}_{}`$. We do not expect there to be large numbers of neutron star or black hole companions to EHB stars. If the companions have a lower mean mass, our detection efficiency would be lower than the value calculated, so the the minimum binary fraction we would derive would be higher. In summary, the minimum binary fraction implied by our observations is about 69$`\pm `$9 percent and this result is not strongly dependent on the assumed distributions of period or mass ratios for short period EHB binaries.

Of course, if there are also binary EHB stars with longer periods, these would not be detected as frequently by our survey as the shorter period binaries, so the binary fraction may be much higher than 2/3. At some point the orbital period is too long for the binary to be relevant to this discussion. The binaries of interest are those for which the orbital separation now is less than the size of a red giant star near the tip of the red giant branch (RGB). In these cases, we can say that the companion has influenced the formation process of the EHB star. The radii of red giants near the tip of the RGB are $`100\mathrm{R}_{}`$, which corresponds to orbital periods of a hundred days or more. We can see from Fig. 4 that most binary EHB stars with orbital periods of tens-of-days would be missed by our survey. Therefore, if about 1/3 of EHB stars are binaries with orbital periods 10 d$`P`$100 d, then our results are consistent with all EHB stars being formed through interactions with a companion star.

### 5.2 Other surveys.

We note that there is no significant radial velocity shift in any of the stars for which there is evidence of a cool companion (PG 0749+658, PG 1018$``$047, PG 1040+234, PG 1647+056 and PG 1701+359). There is a bias in our sample in the sense that the orbital separation of a binary with two 0.5$`\mathrm{M}_{}`$ stars may be too small to contain a main-sequence or sub-giant K-star for the shorter orbital periods where our survey has the greatest sensitivity. However, our observations are still quite sensitive to periods of a day or more, at least for PG 1040+234 and PG 1701+359. The orbital separation for an orbital period of a few days is several solar radii. This may suggest that there is a real trend for EHB stars with K-type companions to have long orbital periods.

Of the six stars observed by SLY which also appear in our sample, we find four to be binaries. Two of these binaries were not detected by SLY (PG 0918+029 and PG 0919+273), one was noted as a marginal detection (PG 0839+399) and one as positive detection (PG 1716+426). The shifts seen by us for PG 0918+029 suggest that it is at the limit of detection for the method employed by SLY. Curiously, the shift we observed for PG 0839+399, for which SLY note a marginal detection, is smaller than that of PG 0918+029, which suggests that our observations have not yet sampled the full range of radial velocity shifts for this star. The shift of only 25.8 km s-1 we measured for PG 0919+273 would not have been seen by SLY.

### 5.3 Post-EHB stars.

Of the stars observed which we have excluded from this discussion, 5 appear to be stars which have evolved away from the extreme horizontal branch. One of these post-EHB stars a good candidate to be a binary from our data (PG 1000+408). Observations by Green (2000) have shown that PG 1000+408 is indeed a binary with an orbital period near 1 day. There are too few stars in this sub-sample to derive useful limits on the binary fraction of these stars, though this would obviously be an interesting number because we would expect it to be similar to the binary fraction for normal EHB stars if the two groups of stars are related as we have suggested.

### 5.4 Selection effects.

One advantage of choosing objects from the PG survey is that we were able to choose brighter stars based on their photographic magnitudes without introducing a bias in our sample towards short period binaries. This is because the majority of the short period binaries have white dwarf or K/M dwarf companions, both of which contribute a negligible amount of light in the blue region of the spectrum on which the PG survey is based. One type of binary we are biased against are those containing brighter F/G-type companions. Most of these binaries were excluded from the PG survey because stars showing the Ca II K line were assumed to be main-sequence subdwarfs with normal colours which appeared bluer due to the substantial uncertainty in the photographic photometry on which survey is based (Green, Schmidt & Liebert 1986). In fact, a substantial fraction of these stars may be sdB stars with F/G-type companions (Kilkenny et a., 1997). Some sdB stars with F/G-type companion were are included in the PG survey, e.g., PG 1631+267. Although we cannot measure the radial velocities of the sdB star from the H$`\alpha `$ line in these cases, we can measure the radial velocities of the F/G star from its H$`\alpha `$ line, if it dominates at the wavelength, or from the many other absorption lines should the H$`\alpha `$ line be a blend of the sdB star and the F/G star H$`\alpha `$ line. In fact, these radial velocities are preferable to the sdB velocities in some ways as they are more accurate and we can estimate the mass of the companion star from its spectral type.

### 5.5 Triple stars

The effect of the selection criterion against F/G companion stars applied to the PG survey is difficult to judge without knowing the actual fraction of sdB stars with F/G-type companions that have been “lost” from the survey. We note that this effect also biases the results of Allard et al. (1994), Ferguson, Green & Liebert (1984) and Jeffery & Pollacco (1998). Nevertheless, these authors still find 1/2 – 2/3 of their sample have cool companion stars. We have already noted that the EHB stars we have observed to be binaries do not have cool companions. The total binary fraction is then 2/3 short period sdB stars without cool companions plus 1/2 – 2/3 with cool companions in the PG survey plus the “lost binaries” with cool companions excluded from the PG survey minus a small fraction of short period sdB stars with cool companions. This number is clearly greater than 1, an apparent paradox which is easily explained by some of the sdB stars being triple stars, i.e., short period sdB stars with M-dwarf or white dwarf companions and a distant F/G-type cool companion which was not involved with the evolution of the inner pair.

## 6 Conclusion

We have measured 205 precise radial velocities for 36 extreme horizontal branch stars to look for variability due to a close binary companion. We found that 21 of our stars are positively identified as short-period binaries. All but one or two of these are new identifications. We conclude that at least 2/3 of all EHB stars are short-period binaries. The orbital separations of these binaries is much less than the size of the star during the red giant phase which almost certainly preceded its emergence as an EHB star. We conclude that some kind of interaction with a binary companion, perhaps in a common envelope phase, is fundamental to the formation process for the majority of EHB stars.

## Acknowledgements

PFLM was supported by a PPARC post-doctoral grant. The Isaac Newton Telescope is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. PFLM would like to thank E.M. Green and R.A. Saffer for many informative discussions on the subject of sdB stars. We would like to thank Rex Saffer for providing us with his results for sdB stars. We would like to thank Martin Altmann for his spectrum of PG 0907+123.