text:
# Focusing a radially polarized light beam to a significantly smaller spot size

## Abstract

We experimentally demonstrate for the first time that a radially polarized field can be focussed to a spot size significantly smaller (0.16(1)$`\lambda ^2`$) than for linear polarization (0.26$`\lambda ^2`$). The effect of the vector properties of light is shown by a comparison of the focal intensity distribution for radially and azimuthally polarized input fields. For strong focusing a radially polarized field leads to a longitudinal electric field component at the focus which is sharp and centered at the optical axis. The relative contribution of this component is enhanced by using an annular aperture.

preprint: APS/123-QED

A large number of optical instruments and devices make use of a sharply focussed light beam. Prominent examples are in lithography, confocal microscopy, optical data storage as well as in particle trapping. A highly concentrated and well matched field is also a necessary requirement for manipulating nanoscopic quantum systems in quantum information processing or cavity quantum electrodynamics KimbleEnk . When reaching for the limits the polarization properties of the electromagnetic field play a dominant role. For a vectorial description of strong focusing using high numerical apertures different methods have been developed RichardsWolf ; MansuripurJOSAA ; KantJModOpt ; SheppardMultipoleJModOpt . It turns out that in general all three mutually orthogonal field components appear at the focal region. The electric energy density patterns of individual polarization components were mapped using single molecules with the absorption dipole axis pointing in the three different directions HechtJMicroscopy ; NovotnyYoungworthPRL ; HechtNovotnyPRL . For linearly polarized light, the energy density distribution of a longitudinally polarized component in the direction of propagation of the beam is not rotationally symmetric. This primarily causes an asymmetric deformation of the focal spot. When using an annular aperture, the relative contribution of the longitudinal component is increased and the asymmetry becomes more pronounced SheppardAnnularApert . These asymmetries were recently confirmed experimentally DornJModOptics . In addition, it was predicted that a special polarization pattern is needed for focusing down to the smallest possible spot QuabisOptComm .

The effects of the vector properties of light on the structure of the focus are best observed if one compares two special input beams with identical intensity distribution but different polarization properties. A good example for two such fields are an azimuthally and a radially polarized field with identical doughnut shaped intensity distributions. When focusing with a high numerical aperture the radially polarized input field leads to a strong longitudinal electric field component in the vicinity of the focus ScullyZubairyPRA ; QuabisOptComm . In contrast, the azimuthally polarized field generates a strong magnetic field on the optical axis NovotnyQuantumDotSpectroscopy  while the electric field is purely transverse and zero at centre YoungworthOpticsExpress . In general, the direction of the polarization may vary largely inside the focal spot. Especially when sub–wavelength–structures are investigated this has to be carefully considered in a quantitative analysis WilsonOptComm . The predicted axial light forces acting on particles in an optical tweezer set up are different in a vectorial treatment as compared to theories based on paraxial and geometric optics NussenzveigOptPinzette . The specific interaction of the different field components at the focus with single molecules can provide information about the orientation of the absorption dipole of the molecules HechtJMicroscopy ; NovotnyYoungworthPRL . Although these field components may also occur in the vicinity of near field apertures BetzigScience , a far field technique can be advantageous as it avoids a modification of the dynamics of the molecule due to an interaction with the probe XieDunnScience ; AmbroseScience . For the spectroscopy of magnetic dipole transitions in quantum dots it was proposed to use an azimuthally polarized field for the excitation and take advantage of the strong magnetic and the vanishing electric field component on axis NovotnyQuantumDotSpectroscopy . To couple efficiently to small quantum systems such as single atoms the incoming field must be well matched in both its amplitude and phase but also in its polarization properties KimbleEnk .

The vector properties not only affect the local field direction but also the intensity distribution at the focus. This is experimentally demonstrated here below by comparing highly resolved measurements for azimuthally and radially polarized input fields. The radially polarized beam can have a narrow central peak due to the appearance of the strong longitudinal field component that is sharply centered around the optical axis. The experimental results reported here confirm for the first time, to the best of our knowledge, that a radially polarized input field leads to the smallest possible spot size in far field focusing observed so far. Three different types of set up have been described with which a radially or azimuthally polarized doughnut mode can be produced 111A radially and azimuthally polarized doughnut mode can be described as a superposition of a TEM01 and a TEM10 Hermite– Gaussian mode with orthogonal polarization.. Most of them either generate the mode inside a laser resonator Pohl ; Niziev  or use a Mach–Zehnder like interferometer Tidwell ; YoungworthSPIE . A third approach involves the use of mode–forming holographic and birefringent elements Tschudi .

Here we use another approach for the generation of this polarization mode. The TEM00–beam of a linearly polarized, single mode helium–neon laser ($`\lambda `$=632.8nm) is expanded and collimated by a telescope with an additional pinhole as a mode cleaner (Fig. 1). The collimated beam is then sent through a polarization converter. This element consists of four half–wave plates, one in each quadrant. The optical axis of each segment is oriented such that the field is rotated to point in the radial direction. However, due to the limited number of segments and some diffraction losses at the boundaries between adjacent segments, the field is not perfectly radially polarized. Nevertheless, this field has a large overlapp with the desired mode and a smaller overlapp with additional higher order transverse modes. These undesired modes are filtered out by sending the beam through a non– confocal Fabry–Perot interferometer that is operated as a mode–cleaner and which is kept resonant only for the radially or azimuthally polarized doughnut–mode. With this method, a purity of $`>`$99% is achieved.

The beam was collimated such that the beam radius at the maximum intensity was 1.2mm. To ensure minimum wave front aberrations we analyzed the beam with a Shack–Hartmann wave front sensor yielding a root mean square wave front error $`<\lambda /35`$. By rotating the polarization of the beam before the polarization converter by 90, it is possible to switch between a radially and an azimuthally polarized beam. For measurements with an annular aperture, a stop is placed directly in front of the microscope objective to block the inner part of the beam. The stop is a high quality glass substrate (surface accuracy of $`<\lambda /20`$) coated with an opaque disc (di=3.0mm or 3.3mm). The diameter of the entrance pupil of the focusing microscope objective (NA=0.9) is 3.6mm. To measure the focal intensity distribution 222In this paper, intensity always refers to electric energy density, which is the part of the field energy that couples to standard photodetectors and photosensitive materials. we used an experimental set up which is based on the knife edge method. The edge is formed by a sharp edged opaque pad that is deposited on the active surface of a photodiode. Therefore no subsequent optical system is needed to collect the transmitted light and effects due to diffraction at the edge are minimized. A more detailed description of the properties of the sensor is given in DornJModOptics . The knife edge method allows one to measure the one dimensional projection i.e. the line integrated intensity of the two dimensional intensity distribution onto the direction along which the edge is moved. By measuring a set of projections onto different directions one acquires the input data needed for a tomographic reconstruction of the two dimensional intensity distribution, using e.g. the Radon back–transformation formalism ImageReconstruction .

For rotationally symmetric input fields the focal spot also shows rotational symmetry and all projections onto any line are identical. Therefore one representative measurement would in principle suffice. But we measured projections onto two orthogonal lines to exclude the possibility that the symmetry of the focussed field is affected by non–spherical aberrations introduced e.g. by the microscope objective. The results were identical within the limits of reproducibility. Line integrals of the intensity distribution were also measured for different positions of the knife edge sensor with respect to the focal plane. This yields a line integrated cross section of the propagating beam containing the optical axis (Fig. 2).

An evidence for the strong longitudinal electric field at the focus of the radially polarized beam is the fact that the intensity projection shows a maximum on the optical axis in the focal plane (Fig. 2a). In contrast, the azimuthally polarized beam has no longitudinal field component and on the optical axis the intensity projection has a minimum. The measured projections in the focal plane are shown in figure 3 for the two modes along with the calculated distributions of the transverse and longitudinal field components. The curves were normalized such that for each case the total intensity in the focal spot equals 1 (arbitrary unit). Cross sections for the tomographically reconstructed intensity distributions (Fig. 3) show that the intensity on the optical axis vanishes for the azimuthally polarized beam. The beam profile at the focus is very similar to the doughnut shaped profile of the input beam just as one would naively expect based on the scalar diffraction theory.

In contrast, the radially polarized mode leads to a narrow peak around the optical axis where the field is basically directed into a longitudinal direction. The transverse field components show a profile very similar to the profile of the azimuthally polarized beam. The power contained in the longitudinal field is calculated to be 49.6% of the total beam power. The relative contribution of the longitudinal component can be enhanced when the numerical aperture of the focusing system is increased.
For a numerical aperture of 0.9 as in this experiment the calculated spot size 333The spot size is defined as the area that is encircled by the contour line at half the maximum value of the intensity. associated with this longitudinal field alone is 0.2$`\lambda ^2`$ which is below the spot size for a linearly polarized input field with homogeneous intensity distribution (0.31$`\lambda ^2`$). The spot size of the total field is still larger (0.47$`\lambda ^2`$) than for linear polarization due to the broad distribution of the additional transverse field component. We note in passing that, in general, a calculation using a scalar field tends to underestimate the spot size for large numerical apertures. The reason is that due to the vector properties of light the plane waves that form the focussed field in image space cannot interfere perfectly at the focus as the electric (and magnetic) field vectors are not all parallel at the focus.

In the case of the radially polarized doughnut mode the rays that propagate under a small angle to the optical axis show only a small longitudinal component. The stronger transverse field components associated with these waves all cancel out on the optical axis and form a broad, doughnut shaped profile which increases the spot size. But when an annular aperture is used for focusing, only waves that propagate under a large angle to the optical axis contribute to the focal field. For a high numerical aperture and a large inner radius of the annular aperture the electric field vectors for a radially polarized

input field are essentially parallel to the optical axis. All rays interfere perfectly and consequently the spot size is close to the spot size which one would calculate using a scalar theory. Cross sections through the tomographically reconstructed intensity distributions are shown for an azimuthally (Fig. 4) and for a radially (Fig. 5) polarized input field when using the annular aperture. Compared to the case without annular aperture the sidelobes are more pronounced in both cases, indicating that the relative contribution of the intensity in the side lobes has increased at the expense of a decreased central maximum.

The spot size for radial polarization with annular aperture is reduced to 0.16(1)$`\lambda ^2`$ (theoretical value 0.17$`\lambda ^2`$) (Fig. 5b). This is well below 0.26$`\lambda ^2`$, the theoretically achievable spot size for linear polarization under the same experimental conditions and also well below 0.22$`\lambda ^2`$, the theoretical value for circularly polarized light. The latter is a little smaller as compared to linear polarization due to the definition of the spot size used here.

Expressed in units of $`\lambda ^2`$ this is to our knowledge the smallest measured spot size that has been reached in far field focusing in air (NA$`<`$1).

It is worth noting that for a radially polarized input field, the magnetic field vectors point in the azimuthal direction just as the electric field vectors do for an azimuthally polarized input beam. The magnetic energy density distribution for a radially polarized beam shows the same distribution as the electric energy density distribution for an azimuthally polarized beam and vice versa.

In summary, we have experimentally verified the influence of the polarization on the shape of the focal spot. For a radially polarized field distribution with annular aperture, the focal spot size was reduced to the smallest focal spot observed up to now when normalizing to the wavelength ($`\lambda ^2`$). The experimentally observed spot size for NA=0.9 is about 35% below the theoretical limit for linearly polarized light. The calculated spot size for the longitudinal field alone is 0.14$`\lambda ^2`$. Therefore, if a surface is covered with a special photosensitive layer which is only sensitive to the narrow longitudinal field distribution QuabisOptComm  even smaller spot sizes can be reached. These ideas may well be combined with other concepts for the increase of transverse resolution such as the solid immersion lens KinoSIL  or coupling to small mesoscopic antennae BouhelierEtAlPRL . A radially polarized field will also provide the input field which is best suited for coupling into the recently proposed all–dielectric waveguide IbanescuScience .

###### Acknowledgements.

We thank G. Döhler, S. Malzer and M. Schardt for the fabrication of the p-i-n diode. We appreciate helpful discussions with H. J. Kimble and C. J. R. Sheppard. We thank J. Pfund for the wavefront measurements with the Shack-Hartmann sensor. This work was supported by the EU grant under QIPC, Project No. IST-1999-13071 (QUICOV).


format: text

text:
# A catalogue of OB associations in IC 1613

## 1 Introduction

A stellar association is ”a single, unbound concentration of early-type luminous stars, embedded in a very young starforming region” (Kontizas et al. 1999). The properties of the associations and the ionized gas clouds in which they are embedded allow tracing of the regions of most recent star formation in the galaxies. IC 1613 is a faint, irregular galaxy within the Local Group, which is resolvable into stars stars from the ground. The young stellar content of IC 1613 was investigated first by Hodge (1978). On the photographic plates taken with the 5 m Hale telescope he identified by eye estimates twenty OB associations. More recently, Freedman (1988) presented color-magnitude diagrams for 11 of them. Georgiev et al. (1999) presented the $`UBV`$ stellar photometry of the northeast sector of the galaxy. The properties of nine Hodge (1978) OB associations in this area are analyzed as well as seven new OB associations determined by cluster analysis. Valdez-Gutierrez et al. (2000) presented $`\mathrm{H}\alpha `$ and $`[\mathrm{SII}]`$ observations using the PUMA scanning Fabry-Perot interferometer. The kinematics of the ionized gas in the complex sample of nebulae was investigated. The ionized gas is distributed in classical HII regions and in a series of superbubbles. They found that almost every superbubble in the NE region has an interior association containing massive stars, suggesting a physical link between them. Rosado et al. (2001) continued the investigation of the influence of the massive stars in the interstellar medium of IC 1613 in the NW and southern region of the galaxy. In the southern region they found that the superbubbles are probably formed by the winds of massive members of associations in spite of the close presence of a WO star. Lozinskaya et al. (2002) obtained spectra of the stars forming the chains and estimated their spectral types and luminosity classes. The stars were found to be at different evolutionary stages and six of them are identified as O stars. Lozinskaya et al. (2003) carried out detailed kinematical studies of the complex of multiple H I and H II shells in region of ongoing star formation in the dwarf irregular galaxy IC 1613.

The purpose of the present paper, the fourth in our IC 1613 series, is to outline the new boundaries of the associations in IC 1613 using observational material of uniform quality and an objective method of identification of stellar associations.

## 2 Observations and data reduction

A set of $`UBV`$ frames of IC 1613 was obtained on the 2-m Richey-Chretien telescope of the Bulgarian National Astronomical Observatory ”Rozhen” and the 2.1 m telescope of the San Pedro Martir observatory in Mexico. The exposure times of the images are between 600 and 1800 sec. The field of view of both telescopes is $`5.6\mathrm{}\times 5.6\mathrm{}`$. The seeing during these observations was $`11.2\mathrm{}`$ with stable and good photometric conditions. Three selected fields in the main body of the galaxy were observed. The Landolt (1992) standards were taken before and after all observations. The IRAF data reduction package was used to carry out the basic image reductions. The stellar photometry of the frames was performed with the point-spread function fitting routine ALLSTAR available in DAOPHOT (Stetson 1993). Complete details of the data reduction and analysis may be found in Georgiev et al. (1999). The zero-point errors of the transformations to the standard $`UBV`$ system are 0.04 mag. We transformed the $`x,y`$ coordinate systems of investigated fields to the one reference local system and then searched for stars in common between datasets. We check the residuals in magnitude and color for the stars in common. They do not show any systematic difference or trend so we can conclude that all datasets are homogeneous in magnitude within the errors. The artificial star technique (Aparicio & Gallart 1995) was used to investigate crowding effects. The completeness can seriously affect our analysis for magnitudes fainter than 21 mag. There are no significant differences in the completeness factors between investigated fields. Since our analysis concerns only bright OB stars, no completeness correction was made. Field stars should not seriously affect the structure of the color-magnitude diagram because IC 1613 is situated far from the Galactic plane and no correction for contamination by field stars was applied. Fortunately IC 1613 has very low reddening - 0.03-0.06 (Freedman 1988, Georgiev et al. 1999), which also dos not affects our results.

## 3 Associations

The automated search for OB associations was carried out using Battinelli’s (1991) algorithm (also called the PLC technique). The method assumes that ”any given two stars in an ensemble of OB stars belong to the same association if and only if it is possible to connect these two stars, by successively linking OB stars located between them, separated from each other by no more than a certain fixed distance parameter, or search radius, called $`d_\mathrm{s}`$”. The distance parameter $`d_\mathrm{s}`$ can be derived from the catalogue of early type stars using the function $`f_\mathrm{p}(d)`$. This function describes the number of groups containing at least $`p`$ stars, for any given value of the distance parameter $`d`$. The function steeply rises until a maximum value is reached, followed by a gradual decrease towards a limiting value of 1, which corresponds to the situation that the value of $`d`$ is so large that all stars can be assigned to one group. The optimum value of the distance parameter $`d_s`$ is defined as the value of $`d`$ corresponding to the maximum of the function $`f_\mathrm{p}(d)`$.

In order to select the catalogue of blue stars we need to adopt brightness and color cut-offs. In IC 1613 we selected as blue stars the stars that have absolute magnitudes $`M_V<2`$ (assuming $`(mM)_0=24.27\pm 0.1`$, Dolphin et al. 2000) and dereddened colors $`(BV)_0<0.0`$ (assuming $`E(BV)=0.03`$, Rienke and Hodge 2001). Recently Pietrzynski et al. (2001) used the same algorithm to investigate associations in NGC 300. They adopted Battinelli’s (1991) method adjusting the two important parameters in the Battinelli algorithm, described in the previous paragraph — distance parameter $`d_\mathrm{s}`$ and the minimum number of stars $`p`$ by numerical analysis. To derive the value of $`d_\mathrm{s}`$ Pietrzynski et al. (2001) applied the PLC technique with several values for the distance parameter and several values for the minimum number of stars, p. Their results show that the position of the maximum of the function $`f_p(d)`$ depends slightly on p and is located in the relatively narrow interval (see their Fig. 2). The second parameter — the minimum number of stars p of the potential association was determined by Pietrzynski et al. (2001) by performing the statistical tests. One hundred random distributions of a number of stars equal to the number of the sample of blue stars, and distributed over the same area, were created and the PLC technique was applied to search for potential groups. This experiment was repeated for different minimum numbers of stars. They found that the spurious detections are less than 10 % if the six stars as a minimum population of OB stars of potential associations are chosen (see their Fig. 3).

Here, we follow Pietrzynski et al. (2001) approach for IC 1613. The value of $`d_s`$ was derived applying the method with different values for the distance parameter. In our case the position of the maximum of the function $`f_p(d)`$ is located in the interval between 17 and 23 pixels. The number of associations obtained with $`d_\mathrm{s}`$ = 17,18…22 and 23 pixels changes by less than $`5\%`$. Thus we confirm the Pietrzynski et al. (2001) result that the number of detected OB associations is insensitive to the adopted value of the search radius in this interval. We adopted $`d_\mathrm{s}=20`$ pixels as a search radius. The next parameter that needs to be specified is the minimum number of stars $`p`$ of the potential association. Again, using the statistical approach of Pietrzynski et al. (2001), described in the previous paragraph we derive four OB stars as a minimum number of stars in some clump in order to have a ”real” association.

With the above described parameters Battinelli’s (1991) algorithm selects 58 associations with sizes between 30 and 130 pc (with accepted distance of 730 kpc, Dolphin et al 2000).

And finally to check our results and to cover the whole galaxy we retrieve the $`UBV`$ images of IC 1613 from the ING Archive, UK Astronomy Data Center, Cambridge. The images were obtained on 11 December 1998 with the 2.5m Isaac Newton Telescope and Wide Field Camera (WFC), which has a field of view of 11.4 x 22.8 arcmin. Unfortunately, they are not of high quality and not as deep as Rozhen images. We reduced the images and transformed the instrumental magnitudes to the standard $`UBV`$ system using about 100 common stars between WFC and Rozhen photometry. Then we ran again Battinelli’s (1991) algorithm in the same way as described above. On the WFC images we found 55 associations. Comparing with the associations obtained from Rohzen images there are 5 omitted associations ( G1, G2, G4, G8 and G26, see Table 1.) and 2 additional (G59 and G60) ones. The WFC images are less deep than Rozhen data by approx. 1 mag in U and there are less than four OB stars in each of the omitted associations. Also, the boundaries of some associations in the fields in common are slightly different. In general the average quantities of the common associations in the two lists however are little affected. The two additional associations are found in the areas not covered by the Rozhen frames.

The map of the association boundaries resulting from the automated search of IC 1613 is shown in Fig. 1 and the properties of the associations are summarized in Table 1. Columns 1, 2 and 3 give the designation and the equatorial coordinates of the detected associations. The number of bright ($`M_V<2`$) OB members in each of the associations, their sizes in parsecs and cross-identification with the catalogue of Hodge (1978) are given in the next columns. The coordinates of the OB associations were determined as the mean coordinates of their separate member stars.

Comparison with the stellar associations outlined by Hodge (1978) shows that most of them coincide with our groups, but in general each Hodge association divides into several smaller groups. The ”new” associations have smaller sizes and look like bright cores within Hodge associations. Hodge (1986) describes the problem of the different sizes of the stellar associations in the Magellanic clouds and other nearby galaxies as M 31 and M 33. He concluded that the mean size of the associations depends on the image scale and the distance to the galaxy. It is known also that young associations in the Galaxy, like the Trapezium, contain small compact subgroups similar to those selected in IC 1613 (Lada & Lada 2003) while the older associations consist of randomly distributed stars and they could not be divided into subgroups. The Battinelli (1991) criterion obviously detects the youngest compact groups and this may explain the smaller mean size of the stellar groups identified in the present paper than those presented by Hodge (1978).

In Fig. 2 (top panel) we show the size distribution of the detected associations. There is a peak located at 60 pc. Comparing the IC 1613 associations with the Magellanic Clouds, M 31, M 33 and NGC 6822 (Bresolin et al. 1998, Ivanov 1996) associations we can see that the distribution of their sizes is similar with a peak between 40 and 80 pc. The mean size of the associations in IC 1613 is 63 pc and is in agreement with the mean sizes of associations in LMC - 60pc, SMC - 70 pc and M33 - 60 pc (Bresolin et al. 1998). All these measurements are performed using the same Battinelli (1991) method. In order to check whether these similar sizes are caused by the identification algorithm we performed a numerical analysis. Using DAOPHOT in IRAF we build the artificial images for our Field II. Field II contains 19 associations with a mean size of 64 pc and an average of 7 OB stars. To conserve the completeness we chose to add 3 artificial associations with 7 ”OB” stars. We created these artificial associations with the PSF obtained for Field II and added them randomly to the image. We then applied Battinelli’s (1991) algorithm. We repeat this procedure many times. As a result we found that the size of an association does not depend on how close the neighboring association is. The same procedure was used to estimate the spurious detections of the associations due to random concentrations of blue stars. Here, we artificially added 21 ”OB” stars ($`10\%`$ of the whole sample in Field II) at random positions and created 20 images. Again, we apply Battinelli’s (1991) algorithm and estimated that between 3-5 OB associations in the whole sample could be randomly concentrated OB stars.

The histogram of the number of OB stars per association is shown in Fig. 2, bottom panel. As can be seen most of the associations have fewer than twelve OB members, with a maximum of the distribution at seven OB members. The most populous association contains 50 OB stars.

## 4 Stellar associations and superbubbles

The associations are not uniformly distributed within the galaxy and they are mostly concentrated in the NE and NW regions (Hodge 1978). On the other hand, the ionized gas in IC 1613 is concentrated in large diameter, expanding superbubbles (Valdez - Gutierrez et al. 2000) in the same regions. Comparison with the velocity map of the superbubbles (Figure 6 in Valdez-Gutierrez et al. 2000) shows that almost every superbubble in IC 1613 has an interior association. These associations are located within or at the periphery of the superbubbles and they are probably physically connected to the gas. The map of the associations with superimposed boundaries of the superbubbles is shown in Fig. 3.

Some of the newly identified stellar associations are found at the boundaries of the superbubbles and/or in the interaction zone between them. This indicates new star formation at the edges of superbubbles. Does this suggest, however, a scenario of self-induced star formation where the shocks that create the larger diameter superbubbles induce the formation of new stars compressing the ambient neutral or molecular gas? And are these more recently formed associations in turn creating new HII regions and bubbles at the periphery of the original superbubble? Lozinskaya et al. (2002) consider that the energy of the stellar wind from the associations is not enough to form the ionized shells and found additional, intense sources of stellar wind at their boundary — early-type supergiants and giants. Lozinskaya et al. (2003) studied in detail the structure and kinematics of the neutral and ionized gas components in the complex of star formation in the NE region. They identified three extended (300350 pc) neutral shells with which the brightest ionized shells in the complex of star formation are associated. They also found evidence of the physical interaction between the H I and H II shells in the region of the chain of stars, early-type giants and supergiants, detected by Lozinskaya et al. (2002). Taking into account the relative positions and ages of the H I and H II shells and OB associations in the complex, they suggested sequential or triggered star formation in the expanding neutral shells. In addition to the three brightest and most prominent H I shells, they found supergiant arches and ring structures in the galaxy whose sizes are comparable to the gaseous-disk thickness and assumed this as a tracer of preceding starbursts in IC 1613. Many OB associations however have no superbubbles. The gas in these associations could be evaporated or perhaps the ionization flux from the stars is not enough to ionize them. In order to clarify these possibilities we are currently in the process of obtaining deeper photometry of the associations and performing spectral identification of the brightest blue stars. The results of such an analysis will be given in a forthcoming paper.

## 5 Summary

We have presented the results of a search for OB associations in the dwarf irregular galaxy galaxy IC 1613. Application of the Battinelli (1991) method resulted in the detection of 60 OB associations with sizes between 30 and 130 pc. Numerical analysis indicates that between 3-5 OB associations in the whole sample could be randomly concentrated OB stars. We detected the expected strong correlation between the spatial distributions of associations and H II regions in IC 1613.

###### Acknowledgements.

J.B is supported by FONDAP Center for Astrophysics grant number 150010003. A part of this work was performed while J.B. was a visiting astronomer in the UNAM, Mexico. The authors gratefully acknowledge the useful comments by an anonymous referee. This research is partially based on data from the ING Archive. L.G. acknowledges the financial support of the CONACyT No. 34422-E. J.B. thanks Marcio Catelan for his help. This research was supported in part by the Bulgarian National Science Foundation grant under contract No. F-812/1998 with the Bulgarian Ministry of Education and Sciences.


text:
# The Different Fates of a Low-Mass X-ray Binary. I: Conservative Mass Transfer

## 1 Introduction

The widely accepted scenario for the formation of millisecond radio pulsars (MSP) is the recycling of an old neutron star (hereafter NS) by a spin-up process. The spin-up is due to the accretion of matter and angular momentum from a Keplerian disc that is fueled via Roche lobe overflow of a binary late-type companion (see Bhattacharya & van den Heuvel 1991 for a review). If enough mass and angular momentum are transferred, the NS spin attains an equilibrium value that is roughly equal to the keplerian angular frequency at the inner rim of the accretion disc (Ghosh & Lamb, 1979). Since the NS has a weak surface magnetic field ($`10^8`$ G), the magnetospheric radius (at which the disc pressure is balanced by the magnetic pressure) truncates the accretion disc close by or at the NS surface, and the equilibrium period is expected to be, in most cases, below one millisecond. Typically $`0.35M_{}`$ are sufficient to reach a spin period of 1 ms (e.g. Burderi et al. 1999). Most donor stars in systems hosting recycled MSPs have certainly lost, during their interacting binary evolution, a mass greater than $`0.35M_{}`$ since they now appear as white dwarfs of mass $`0.150.30M_{}`$ (e.g. Taam, King, & Ritter 2000), whose progenitors are likely to have been stars of $`1.02.0M_{}`$ (Webbink, Rappaport, & Savonije 1983; Burderi, King, & Wynn 1996; Tauris & Savonije 1999). Therefore if the mass transfer is conservative the amount of matter accreted is well sufficient to spin the star up to periods below 1 ms (Cook, Shapiro, & Teukolsky 1994a), or even to produce an accretion induced collapse into a black hole. Once the accretion and spin-up process ends, the magnetospheric radius expands beyond the light cylinder radius (where an object corotating with the NS attains the speed of light). This initiates a phase in which the rotational energy of the NS is emitted via electromagnetic radiation and the star can be observable as a rapidly rotating radio pulsar. According to this model, Low Mass X-ray Binaries (hereafter LMXBs) are the progenitors of MSPs. Indeed the discovery of coherent X-ray pulsations in four transient LMXBs, namely SAX J1808.4–3658 with a spin period $`P=2.5\mathrm{ms}`$ and an orbital period of $`P_{\mathrm{orb}}=2\mathrm{h}`$ (Wijnands & van der Klis 1998), XTE J1751–305 ($`P=2.3\mathrm{ms}`$, $`P_{\mathrm{orb}}=42\mathrm{min}`$, Markwardt et al. 2002), XTE J0929–314 ($`P=5.4\mathrm{ms}`$, $`P_{\mathrm{orb}}=43\mathrm{min}`$, Galloway et al. 2002), and XTE J1807–294 with a spin period $`P=5.2\mathrm{ms}`$ and an orbital period of $`P_{\mathrm{orb}}=40\mathrm{min}`$ (Markwardt et al. 2003) confirmed that NSs in LMXBs can be accelerated to millisecond periods.

Several numerical methods have been developed to solve the Einstein equations for a rotating NS (see Stergioulas 1998 for a review). Stability criteria show that a rapidly rotating NS can support a maximum mass (against gravitational collapse) much higher than the non-rotating mass limit, since the centrifugal force attenuates the effects of the gravitational pull (e.g. Friedman, Ipser, & Sorkin 1988). Conversely, if a rotating NS has a mass that exceeds the non-rotating limit (i.e. a supramassive NS), it will be subject to gravitational collapse if it loses enough rotational energy. Numerical simulations of rotating NSs show that, contrarily to the standard behavior, supramassive NSs spin up just before collapse, even if they lose energy. This effect is known to be stronger for higher mass objects (Cook, Shapiro, & Teukolsky 1992).

The value of the maximum rotating and non-rotating mass depends on the equation of state (EOS) governing the NS matter. On the other hand, the minimum allowed period for a given mass occurs when gravity is balanced by centrifugal forces at the NS equator (mass shedding limit). Thus the spin period can be used to constrain the mass-radius relation for the NS, i.e. its EOS. In the context of the standard (gravitationally-bound) NS models (e.g. Glendenning 2000), several EOSs have been proposed. We usually distinguish different EOSs depending on their stiffness (i.e. the value of $`\mathrm{d}p/\mathrm{d}ϵ`$, were $`p`$ is the fluid pressure and $`ϵ`$ is the energy density). If the EOS is stiff, the matter is less compressible at high densities, resulting in larger NS radii as compared to a soft EOS, and hence in longer minimum rotational periods. Except for few, very stiff cases, most EOS predict minimum rotational periods well below 1 ms. However, no sub-millisecond pulsars have been detected up to date: the shortest observed spin period is $`1.5\mathrm{ms}`$ (Backer et al. 1982), uncomfortably higher than the theoretical predictions.

In an attempt to find an explanation for the apparent clustering of the spin periods of millisecond pulsars around 2 ms, Bildsten (1998) and Andersson (1998) indipendently suggested that LMXBs emit gravitational waves once they reach a critical spin frequency. Burderi and D’Amico (1997) showed that for nonaxisymmetric m-modes, assuming a realistic range of temperatures, the values of the critical spin frequency are remarkably close to the limiting spin frequency determined by the centrifugal limit at the border of the NS. On the other hand, Andersson, Kokkotas and Stergioulas (1999) demonstrated that at a certain spin frequency (much lower than the maximum attainable spin period) an instability to the Rossby waves (r-modes) of the star arises, thus causing emission of gravitational waves. Levin (1999) suggested that the gravitational waves emission causes the onset of a spin-up spin-down cycle, and not of a steady state spin equilibrium: in this scenario the NS undergoes a very rapid spindown (lasting $`1`$ yr) due to the rapid heating during the r-mode excitation, and then starts another cycle of accretion driven spin up. Brown and Ushomirsky (2000) showed that if the NS has a superfluid core the steady state scenario is not viable because the predicted quiescence luminosity in this case is much higher than the observed one in the X-Ray transient Aql X-1. Levin and Ushomirsky (2001) showed that, when keeping into account the presence of the solid crust, the critical spin frequency for the onset of the r-mode varies between $`600`$ Hz and $`200`$ Hz, depending on the core temperature.

In this paper we consider the full evolution of a LMXB, trying to determine how the results of our modeling of the recycling scenario compare to the observations and which effects peculiar to General Relativity are indeed observable. We will also show the differences in the evolution of the system when the r-mode instability is excited during accretion and when it is not excited.

## 2 Evolution equations for the compact object

A rotating NS is unambiguously defined by the boundary conditions for the integration of the Einstein equations or, equivalently, by suitable pairs of resulting integrated quantities, such as the baryonic mass and the angular momentum or the baryonic mass and the total mass–energy of the star. Therefore the evolution of the NS is determined by the evolution of such pairs. We consider the evolution of the NS both during the accretion phase and after the accretion has finished: in the former case it is convenient to solve the evolution equations for baryonic mass and angular momentum, since we know the general formula for the torque exterted by the accreting matter on the NS, while in the latter we solve the evolution equations for baryonic mass and mass–energy, since we can evaluate the luminosity of a magnetodipole rotator using Larmor’s formula.

### 2.1 Evolution equations for the accretion phase

According to accretion theories, matter transferred from the companion star to the NS via Roche-Lobe overflow forms an accretion disc around the compact object. The disc is truncated because of one of the following reasons: i) the interaction with the magnetic field of the NS, which truncates the disc at the magnetospheric radius $`r_m`$; ii) the presence of NS surface itself at $`R_{\mathrm{NS}}`$; iii) the lack of closed Keplerian orbits for radii smaller than the marginally stable orbit radius, $`R_{\mathrm{MSO}}`$ (at few – depending on the mass and spin of the compact object – Schwarzschild radii from the NS centre).

The magnetospheric radius is defined as the radius at which the pressure of the disc equals the pressure of the magnetic field of the NS. The magnetospheric radius can be written as a fraction $`\varphi `$ (see Burderi et al. 1998) of the Alfvén radius $`R_\mathrm{A}`$ (the radius at which the energy density of the, assumed dipolar, NS magnetic field equals the kinetic energy density of the spherically accreting matter) as

$`r_\mathrm{m}`$ $`=`$ $`\varphi R_\mathrm{A}=2.45\times \alpha ^{9/35}m_G^{1/28}\dot{m}_B^{3/70}r_\mathrm{m}^{3/28}R_\mathrm{A}`$ (1)
$`=`$ $`2.244\times \left(\alpha ^{36/5}m_G\dot{m}_B^{6/5}R_\mathrm{A}^{28}\right)^{1/31}\mathrm{cm}`$
where $`\alpha 1`$ is the Shakura-Sunyaev parameterization of the accretion-disc viscosity (for which we will assume a typical value of $`0.1`$), $`m_G`$ is the NS gravitational mass in $`M_{}`$, $`\dot{m}_B`$ is the baryonic accretion rate on to the NS in $`10^8M_{}/\mathrm{yr}`$, and $`R_\mathrm{A}`$ is

$$R_\mathrm{A}=1.24\times 10^6\mu _{26}^{4/7}m_G^{1/7}\dot{m}_B^{2/7}\mathrm{cm}$$
(2)
where $`\mu _{26}`$ is the magnetic dipole moment of the NS in units of $`10^{26}\mathrm{G}\mathrm{cm}^3`$ , defined from $`B_s=\mu /R^3`$, where R is the equatorial radius and $`B_s`$ the surface magnetic field of the NS at its equator.

The NS radius, that is always in the order of $`10^6\mathrm{cm}`$, depends both on the mass of the NS and on its angular momentum. In general smaller radii correspond to larger masses while larger radii correspond to larger angular momenta. Thus a rapidly rotating NS can have a much larger radius than a non-rotating one (the equatorial radius can expand up to $`40\%`$, see Cook, Shapiro, & Teukolsky 1994b).

The marginally stable orbit is the smallest stable orbit possible for a test particle around an axisymmetric, rotating body of gravitational mass $`M_G`$ and angular momentum J. Its radius, following Bardeen (1970), is:

$`R_{\mathrm{MSO}}`$ $`=`$ $`R_\mathrm{g}(3+Z_2\sqrt{(3Z_1)(3+Z_1+2Z_2)});`$
$`Z_2`$ $`=`$ $`\sqrt{3\left({\displaystyle \frac{a}{R_\mathrm{g}}}\right)^2+Z_1^2}`$
$`Z_1`$ $`=`$ $`1+\left(1+\left({\displaystyle \frac{a}{R_\mathrm{g}}}\right)^2\right)^{1/3}+\left(1{\displaystyle \frac{a}{R_\mathrm{g}}}\right)^{1/3}`$
$`a`$ $`=`$ $`{\displaystyle \frac{J}{M_Gc}}`$
$`R_\mathrm{g}`$ $`=`$ $`{\displaystyle \frac{GM_G}{c^2}}.`$
where $`c`$ is the speed of light and $`G`$ is the gravitational constant. Therefore the inner radius of the disc, $`R_D`$, will be:

$$R_D=\mathrm{max}\{r_\mathrm{m},R_{\mathrm{NS}},R_{\mathrm{MSO}}\}.$$
(3)
The position of the inner rim of the disc is crucial for the accretion of matter on to the NS. In particular it is very important in the case $`R_D=r_\mathrm{m}`$: since at $`r_\mathrm{m}`$ the matter is forced by the magnetic field to corotate with the NS, accretion on to a spinning magnetized NS is centrifugally inhibited once $`r_\mathrm{m}`$ lies outside the corotation radius $`r_{\mathrm{co}}`$, the radius at which the Keplerian angular frequency of the orbiting matter is equal to the NS angular frequency: $`r_{\mathrm{co}}=1.50\times 10^6m_G^{1/3}P_3^{2/3}`$ cm, where $`P_3`$ is the spin period in milliseconds. Conversely if the magnetospheric radius is smaller than the corotation radius accretion of matter can proceed undisturbed. However, as $`r_m`$ scales as a negative power of $`\dot{m}_B`$, a decrease in the mass transfer rate, which can occur for instance at the end of the accretion phase, will results in an expansion of the magnetosphere. In this paper we considered NSs with $`\mu _{26}1`$, which is typical of observed MSPs (see e.g. Lorimer 1994), corresponding to surface magnetic fields of $`10^8\mathrm{G}`$. Moreover, in the cases considered here, the accretion rate is always above $`10^{10}M_{}/\mathrm{yr}`$ during the mass transfer phase (see figure 1). In this case, eq. (1) implies a magnetospheric radius of $`10^6`$ cm, that is comparable with the NS radius. Although we considered in our simulations (see below) the possibility of a propeller phase, in these conditions the details of the interaction of the accretion disc matter with the NS magnetosphere during the accretion phase will not affect much our results. In fact the magnetospheric radius lies outside of the corotation radius only for a small time ($`0.1\%`$ of the evolution time) just before of the end of accretion (see figure 2 for details). In other words, the spin evolution of a low magnetized NS is not very different from the spin evolution of a non–magnetized NS, in which no propeller effect is possible since the disc is always truncated near the surface of the NS. We can therefore write the mass evolution equation for $`r_m<r_{co}`$, which is simply given by the condition that the baryonic mass per unit time accreted on to the NS is equal to the mass lost by the companion per unit time due to Roche-Lobe overflow, $`\dot{M}_{RL}`$, i.e.

$$\dot{M}_B=\dot{M}_{RL}.$$
(4)
Lamb, Pethick, & Pines (1973) wrote a general equation describing the flow of angular momentum into the stellar magnetosphere including the material stress at the inner edge of the disc as well as magnetic and viscous stresses in a disc partially threaded by the NS magnetic field. In particular they assumed that the threading occurs in a transition region near the inner rim of the disc where the magnetic field of the NS couples to the accretion disc. However they noticed that, for slow rotators, magnetic and viscous stresses can be neglected with respect to the contribution of the material stress. Ghosh, Lamb, & Pethick (1977) showed that, in the case of a disc rotating in the same sense of the star, if the transition region where the magnetic field couples to the disc is small, these additional contributions are still negligible even for rapidly rotating NS. In the present discussion we assume therefore that the torque exerted on the NS is only due to the contribution of the material stress at the inner edge of the disc (Pringle & Rees 1972). We will study the evolution of high-magnetized NSs with heavily threaded accretion discs in a future paper. Note, however, that Wang (1996) showed that threading effects can modify the maximum achievable period (i.e. the equilibrium period) of the NS only by a few percents with respect to the unthreaded case.

The angular momentum per unit baryonic mass of a particle orbiting around a rotating, axisymmetric object at a distance R is (see Bardeen 1970)

$$j(R)=\frac{(GM_G)^{\frac{1}{2}}(R^22(J/M_Gc)(RGM_G/c^2)^{\frac{1}{2}}+(J/M_Gc)^2)}{R^{\frac{3}{4}}(R^{\frac{3}{2}}3(GM_G/c^2)R^{\frac{1}{2}}+2(J/M_Gc)(GM_G/c^2)^{\frac{1}{2}})^{\frac{1}{2}}}.$$
(5)
Thus we can write the equation for angular momentum evolution of the NS simply as

$$\{\begin{array}{c}\dot{J}=j(R_D)\dot{M}_B,J<J_{max}(M_B)\\ \dot{J}=\frac{\mathrm{dJ}_{\mathrm{max}}}{\mathrm{dM}_\mathrm{B}}\dot{M}_B,J=J_{max}(M_B)\end{array}$$
(6)
The first equation does not apply when the NS is at the mass shedding, i.e. the regime in which the matter at the border of the star has the Keplerian velocity at that radius; for each value of baryonic mass, the mass shedding regime is individuated by the corresponding maximum angular momentum of the star, $`J_{max}(M_B)`$. At mass shedding, the matter of the disc should dissipate angular momentum to accrete on to the star. If the disc matter has an angular momentum that will make the star exceed the mass shedding limit, it will just not be able to accrete since it will not be gravitationally bound to the star, until viscous forces drive the excess of angular momentum to the outer zones of the disc (where it will be given back to the companion through tidal forces) allowing matter to accrete on to the star. In this situation, the second equation holds. Integrating it we get $`J_fJ_i=J_{\mathrm{MAX},f}J_{\mathrm{MAX},i}`$, but since $`J_i=J_{\mathrm{MAX},i}`$ we have $`J_f=J_{\mathrm{MAX},f}`$: the accreted matter will give the NS only the angular momentum that keeps it at mass shedding, and the star will continue to move along the maximum rotation line.

### 2.2 Evolution equations at the end of the accretion phase

The NS is thought to switch on as a radio pulsar when the inner edge of the disc lies outside the light cylinder radius (i.e. the radius at which a particle corotating with the star will have velocity $`c`$), $`r_{lc}=c/\omega _{NS}`$ \- where $`\omega _{NS}`$ is the angular velocity of the NS: this certainly happens when the accretion stops and thus the disc disappears. The emission mechanisms for a radio pulsar are believed to be rotating magnetic dipole radiation and magnetospheric currents associated with the emission of relativistic particles, both depending on the angle i between the NS magnetic moment $`\mu `$ and the spin axis (Goldreich & Julian 1969). These two emission mechanisms compensate in such a way that the total energy emitted is nearly independent of i (Bhattacharya & van den Heuvel 1991). The energy loss per unit time will be $`\dot{E}=\frac{2}{3c^3}\mu ^2\omega _{NS}^4.`$ It is then convenient to describe the NS evolution in terms of baryonic mass and total energy rather than in terms of baryonic mass and angular momentum as we did for the accretion phase. We can write the evolution equations as

$`\dot{M}_B`$ $`=`$ $`0`$
$`\dot{M}_G`$ $`=`$ $`{\displaystyle \frac{2}{3c^5}}\mu ^2\omega _{NS}^4.`$ (7)
Depending on the different conditions, we use either equations (4) and (6) or equations (2.2) to compute the spin evolution of the NS.

### 2.3 Excitation of r-modes

Since our main goal is to compare qualitatively the evolution of systems in which the gravitational waves damping is not present and systems in which it is present, we will not go into the details of the theory of r-modes excitation. We will simply assume that, if the r-mode mechanism is present, any star attaining a period of 2 ms (corresponding to a frequency of 500 Hz) undergoes a rapid spindown to a period of 5 ms due to gravitational waves emission following the scenario proposed by Levin (1999). After this rapid spin-down, the NS starts again its accretion driven spin up process.

### 2.4 Numerical methods of integration

During accretion phases we coupled through equations (4) and (6) a detailed description of the binary evolution of the system obtained with the stellar evolution code with a fully relativistic calculation of NS physical properties given its EOS, its baryonic mass and its angular momentum.

The evolution of the binary system is followed selfconsistently including the full computation of the structure of the secondary star, by means of the ATON1.2 stellar evolution code (D’Antona, Mazzitelli & Ritter, 1989). The equations of stellar structure are numerically solved by a full Newton–Raphson integration from the center up to the basis of the stellar atmosphere. The numerical inputs are described in Mazzitelli (1989). The secondary star is considered to be the mass losing component of the binary system, and its mass loss rate is computed following the formulation by Ritter (1988), as an explicit exponential function of the distance of the stellar radius to the Roche lobe, in units of the pressure scale height. The evolution of the binary parameters can be followed by considering several possible cases for the transfer of mass, and for the loss of mass and associated angular momentum from the system. The orbital evolution also includes losses of orbital angular momentum through magnetic braking, in the Verbunt and Zwaan (1981) formulation, in which the braking parameter is set to $`f=1`$ and through gravitational radiation.

The relativistic computations for the compact object are implemented using a slightly modified version of RNS (Rotating Neutron Stars) public domain code by Stergioulas & Morsink (1999). The RNS routines provide the numerical solution of Einstein equations for a rotating axisymmetric body, integrating the equations via Komatsu-Eriguchi-Hachisu method (Komatsu et al. 1989; see Stergioulas & Friedman 1995 for a comparison with other integration methods). The main problem of this approach is that we have evolution equations for baryonic mass and angular momentum of the compact object, while the boundary conditions for the solution of Einstein equations are the central energy density and the oblateness of the star. To solve this problem, we used a grid of relativistic equilibrium NS models integrating the Einstein equations for a wide range of initial conditions, spanning all allowed values for stable configurations. This was obtained modifying the code to build a complete grid of relativistic equilibrium configurations with the necessary numerical precision. To this we added a stability control routine to exclude unstable configurations.

We obtain the accretion rate $`\dot{M}_B`$ and a corresponding time interval (in which $`\dot{M}_B`$ remains unchanged) from the stellar evolution code. Then we integrate the differential equations (4) and (6) using a finite differences method. For each integration time step the accreted baryonic mass is $`\mathrm{\Delta }M_B=\dot{M}_B(t)\mathrm{\Delta }t`$ and the accreted angular momentum is $`\mathrm{\Delta }J=j(R_D(t))\dot{M}_B(t)\mathrm{\Delta }t`$. For the $`n`$-th time step of the evolution, we search on the grid the single equilibrium configuration with the updated values of $`M_B(t_n)=M_B(t_{n1})+\mathrm{\Delta }M_B`$ and $`J(t_n)=J(t_{n1})+\mathrm{\Delta }J`$; in this way we get the corresponding values of gravitational mass, radius, spin frequency and momentum of inertia of the NS. We fine-tuned the time step so that $`\mathrm{\Delta }M_B`$, $`\mathrm{\Delta }J`$ are always comparable with the distance of two neighbour points in the grid.

Obviously the integration over a grid may introduce numerical uncertainties in our results: in particular, we are implicitly assuming that the spin frequency, the gravitational mass and the radius of the NS remain unchanged between two points on the grid. Therefore we introduce an error in the evaluation of these quantities that is equal to half the distance between two points on the grid, that is always well below $`0.5\%`$ in our simulations. This also yields an uncertainty in evaluating, step by step, the value of the accretion radius given by eq. (3) and the angular momentum of the matter at the inner rim of the disc, that is always in the order of a few percent.

During the pulsar phase, we integrated equations (2.2) in a similar way using the same grid of relativistic equilibrium NS models. To avoid the grid imprecisions we made a cubic spline interpolation on the values of the function $`\omega _{NS}(M_G)`$ for the integration of the equation set (2.2).

## 3 Results

We consider three simple examples that show different possible fates for a LMXB, depending on the characteristics of the companion star. In all cases we start our evolution with a slowly rotating ($`P1\mathrm{s}`$) NS with a gravitational mass of 1.4 $`M_{}`$ – which is a typical value for a newborn NS (Thorsett and Chackrabarty, 1999). The EOS adopted to describe the ultradense matter is that proposed by Friedman, Pandharipande and Smith (FPS, see e.g. Lorenz et al. 1993 for a recent discussion) which has been widely used in the literature and has average values of stiffness and physical parameters when compared to other EOS (see Arnett and Bowers, 1977 for a catalog). For our simulations we consider the maximum mass configuration equal to the maximum rotation configuration (although this is not precisely true, see Stergioulas & Friedman 1995, these configurations differ one from each other for less than $`0.1\%`$ in our case), and we have $`M_{\mathrm{max}}=2.123M_{}`$, $`\omega _{\mathrm{max}}/2\pi =1882\mathrm{Hz}`$. Finally, as already mentioned, we assume that the NS has a magnetic dipole moment of $`10^{26}\mathrm{G}\mathrm{cm}^3`$.

In particular we consider three typical examples of binary evolution in which the accretion rate is low enough to remain conservative, both above and below the bifurcation period – i.e. the period below which the orbital evolution, during the accretion phase, proceeds towards shorter binary periods. The bifurcation period is $`P_{\mathrm{orb},\mathrm{bif}}18`$ hr for a binary composed of a $`1.4M_{}`$ NS plus a $`1.0M_{}`$ secondary (Podsiadlowski et al. 2002). In the cases we considered, the secondary is:

1. A population II donor with an initial mass of $`0.85M_{}`$ and an orbital period above the orbital bifurcation period; it loses $`0.646M_{}`$ during the accretion phase. The initial binary period is $`14.3`$ hours and it steadily increases to about $`124`$ h at the end of the accretion phase, which lasts $`1.7\times 10^9`$ yr (see figure 1a), leaving a helium white dwarf remnant. Note that the orbital bifurcation period for a system composed of a $`1.4M_{}`$ NS plus a $`0.85M_{}`$ is smaller than $`14.3`$ hours, because the bifurcation period depends strongly on the companion mass. In this case, Roche-Lobe overflow is driven by the nuclear evolution of the companion.
2. A population I donor (initially below the bifurcation period) with an initial mass of $`1.15M_{}`$, which loses $`0.91M_{}`$ during the accretion phase. The orbital period evolves from $`10.6`$ to $`3.5`$ hours. In this case the companion star overflows its Roche-Lobe due to angular momentum losses caused magnetic braking. The accretion phase ends when the companion becomes fully convective, according to the classical scenario for the evolution of cataclysmic variables in the period gap (Verbunt & Zwaan 1981). The accretion phase lasts $`1.7\times 10^9\mathrm{yr}`$, see figure 1b.
3. A population I donor (initially above the bifurcation period) with an initial mass of $`1.199M_{}`$; it loses $`0.99M_{}`$ during the accretion phase. The period evolves from $`19`$ to $`92`$ hours, leaving again a helium white dwarf remnant. The accretion phase lasts $`2.2\times 10^9\mathrm{yr}`$ (see figure 1c). In this case, Roche-Lobe overflow is driven by the nuclear evolution of the secondary star.

In these evolutions the mass trasfer rate is smaller than the Eddington limit. Therefore it is reasonable to assume that the mass transfer is conservative. Moreover, the matter flow from the companion is not subject to large fluctuations. The NS mass and spin evolution in more complicated cases, like the ones observed in X-ray transients (see Campana et al. 1998 for a review) or the ones discussed by Burderi et al. (2001), are not considered here and will be discussed in a following paper. We considered both the case in which the accretion process does not excitate the r-mode instability, and the case in which the spin evolution is influenced by the gravitational waves emission.

In figure 3 we plot the evolution of these systems in the gravitational mass – radius plane. The area of this plane that a stable NS equilibrium configuration can span is limited from below by the sequence of stable non-rotating equilibrium configurations, which are the ones with minimum radius for a given gravitational mass, and from above by the mass-shedding sequence – that is the sequence of equilibrium configurations for which the NS angular velocity equals the Keplerian angular velocity at the NS surface radius. These two stability lines are connected by a third curve, the secular stability curve, which limits from the right the region in which equilibrium configurations are stable against gravitational collapse.

The system evolution without gravitational waves emission are indicated by the thick solid and dashed lines in figure 3, where the numbered dots give the corresponding age and spin period. We see that the first phase of accretion (thick solid lines) is similar in all cases: a rapid spin up brings the star to the mass-shedding limit. Thus the star remains on the mass shedding sequence in the M–R diagram; it will detach from it only when the accretion stops and the NS switches on as a radio pulsar (thick dashed lines).

When we consider the spin-down due to r-modes excitation, all evolutionary tracks (dot-dot-dashed lines in figure 3) are even more similar: after a brief spin up phase, the NS attains the critical spin period of 2 ms. After that, the spin-up spin- down cycle begins and it ends either in a direct collapse to a black hole or in a radio pulsar phase. In both hypotheses, the amount of mass accreted before the pulsar lights up will determine the fate of the system.

In figure 4 we show the possible behaviors of NS angular velocity as a function of gravitational mass for given values of baryonic mass. From equation (2.2) we see that during the pulsar phase the spin evolution of the star will only depend on the baryonic mass of the star and on the sign of the derivative $`(\omega /M_G)_{M_B}`$. If there is a non-rotating stable configuration for the given baryonic mass, then $`(\omega /M_G)_{M_B}>0`$ always and the pulsar will spin down until it stops. If such a configuration does not exists for the given baryonic mass, two different evolutionary tracks are possible: either the pulsar spins down until it comes close to the stability limit to gravitational collapse, where it spins up rapidly (in fact $`(\omega /M_G)_{M_B}\mathrm{}`$ at the instability), or it spins up until it collapses, with no spin down phase (i.e. $`(\omega /M_G)_{M_B}<0`$ always). We say that in the former case the pulsar lights up in the supramassive regime, while in the latter it lights up in the extremely supramassive regime. The spin-up phase, that occurs in correspondence to a loss of rotational energy, is caused by the fact that, near the onset of the instability and in the extremely supramassive regime, a little loss of energy (gravitational mass) corresponds to a rapid decrease of the NS radius. Since the momentum of inertia of the NS strongly depends on NS radius R and on its compactness $`GM_G/Rc^2`$ ($`I0.21M_GR^2/(12GM_G/Rc^2)`$, see Ravenhall & Pethick 1994), in these regimes the loss of rotational energy is achieved by means of a reduction of the momentum of inertia rather than via a spin down of the NS, that, instead, will spin up to partially compensate the reduction of the momentum of inertia.

In case a, in the absence of r-mode excitation, the pulsar lights up, with a period of 0.71 ms, when the star has a gravitational mass larger than the maximum non-rotating mass i.e. $`M_{\mathrm{stat}}=1.803M_{}`$ (see figure 4i). As the pulsar loses energy due to dipole radiation according to equation (2.2), it leaves the mass-shedding sequence and returns in the normal gravitational mass range, below $`M_{\mathrm{stat}}`$. Thus the pulsar will slow down (see figure 3a), without collapsing, until it almost stops on a very long time-scale (it reaches $`16\mathrm{ms}`$ in $`8\times 10^{11}`$ yr). If the r-modes are excited during the accretion process, the star lights up as a pulsar with the same baryonic mass (since the same amount of mass has been accreted), but with a different gravitational mass, 1.792 M instead of 1.805 M: the larger rotational period of the newborn pulsar, 3.31 ms, with respect to the period of 0.71 ms of the evolution in absence of gravitational waves, means that there is not enough energy to increase the gravitational mass above the maximum non-rotating mass. In this case, the pulsar slows down on a similar long timescale, reaching $`16\mathrm{ms}`$ in $`7\times 10^{11}`$ yr.

In case 2, if r-modes are not excited the star becomes a pulsar when it is in the extreme supramassive regime (see figure 3b) and thus equations (2.2) imply that the NS is spinning up rather than spinning down. As it loses energy, however, the star begins to shrink and heads towards the secular instability limit which brings the NS to the gravitational collapse. Its radio pulsar phase, characterized by periods well below 1 ms (its initial spin period is 0.581 ms, its period when it crosses the stability line is 0.542 ms) and by the unusual sign of the period derivative, lasts $`5\times 10^7\mathrm{yr}`$. Such a short lifetime, if compared with typical lifetimes for spinning down pulsars, is due to the positive feedback we get in this case for the second of equations (2.2): gravitational mass loss causes a spin up, which in turn will make the term on the right of this equation bigger, causing a further mass loss. If the r-modes are excited no such behavior arises, since spin of the NS is kept above 2 ms, the maximum sustainable mass is smaller than that for a NS rotating at mass shedding (the maximum sutanaible gravitational mass for a NS spinning at a period of 2 ms is $`1.817M_{}`$). Thus the star collapses to a black hole $`5\times 10^8`$ yr before the end of the accretion process.

In case c, the star never lights up as a pulsar, as too much matter is accreted on the star and the maximum mass limit is exceeded. Thus the NS will directly collapse to a black hole (see figure 3c).

All this behaviors are strongly dependent on the EOS adopted to describe the ultradense matter. As a comparison we studied the evolution of a system consisting of the same companion star as in case c and a $`1.4M_{}`$ NS with an ultra-stiff EOS (EOS N by Walecka & Serot, see Arnett and Bowers 1977), which has a maximum non-rotating mass of $`2.634M_{}`$. The result of this evolution different, as we show in figure 3d; in this case (hereafter case d) no spinning up pulsar shows up, nor any accretion induced collapse happens. Instead, we end up with a spinning down submillisecond pulsar with a very long lifetime, comparable with the one we obtained in the first case. If the r-modes are excited (see again figure 3d) the minimum attanaible period is limited to 2 ms, and the pulsar lights up with a period of 2.27 ms. In a following paper, we will discuss in full detail the effects of the EOS of the ultradense matter on the evolution of LMXBs.

## 4 Conclusions

We have shown that, depending on the characteristics of the system, especially on the amount of mass accreted on to the NS, on the EOS adopted to describe the Ns matter, and on the excitability of the r-modes, LMXBs can have quite different fates: they can light up as a spinning down radio pulsars, they can directly collapse to a black hole during the accretion phase or if, at the end of the accretion phase, the NS is left in the extreme supramassive regime, it will light up as an exotic, spinning up submillisecond radio pulsar with a relatively short lifetime.

It is then evident that, in the hypothesis of a conservative mass transfer from the companion onto a low–magnetized NS and in absence of r-modes excitation, the accretion process, if the amount of mass accreted is not enough to collapse into a black hole, will end with a very fast spinning object, as it has been suggested before (Cook, Shapiro, & Teukolsky 1994a). If the NS becomes then detectable as a radio pulsar, it will have a spin period well below one millisecond. In fact our simulations show (cases a and d) that we obtain submillisecond pulsars with long lifetimes (in the former case the pulsar lifetime before reaching a period as long as that one of the fastest millisecond pulsar known to date, PSR B1937+21, is $`3\times 10^9\mathrm{yr}`$, while in the latter it becomes $`5\times 10^9\mathrm{yr}`$).

If r-modes are excited by accretion, pulsars are constrained to spin slower than a critical frequency, and this could explain why no NS spinning at submillisecond periods has been observed to date. However, in this situation any binary system in which enough mass is transferred from the companion to the NS will collapse to a black hole, without lighting up as a pulsar. Thus pulsar formation could be much less favoured than in other cases.

It is likely that millisecond pulsar systems like the ones observed to date (i.e. systems with $`P>1\mathrm{ms}`$) originate from different binary evolution scenarios, in which some critical mechanism has prevented the accretion process to continue until a mass as large as a significant fraction of a solar mass has been transferred. It is probable that the magnetic field of the NS has values much higher than the value we chose (at least at the beginning of the accretion), so that the inner edge of the disc could be outside the corotation radius for at least part of the evolution and magnetic torques could play an important role in the spin evolution of the NS. Moreover, systems in which the mass transfer rate has large fluctuations will light up as pulsars before the end of the accretion process, losing a large amount of mass in a so-called radio ejection phase as proposed by Burderi et al (2001). We will investigate the evolution of such systems in a future paper.

Although there are selection effects that could have prevented the discovery of submillisecond pulsars (Burderi et al 2001), if a self-limiting mechanism like the r-mode instability does not operate, submillisecond radio pulsars should exist and should be detectable in the future. On the other hand we have shown (case 3) that if the NS matter is governed by a moderately stiff EOS like FPS (i.e. with a maximum non–rotating mass of $`1.9M_{}`$), the mass transfer can end in an accretion induced collapse to a black hole if as much as $`1M_{}`$ is accreted. Although not much in known on the range of progenitor masses fot today’s population of LMXB in the Galaxy, a recent study by Pfahl, Rappaport and Podsialowski (2003) argue that a great fraction of observed LMXBs may have descended from intermediate mass X-ray binaries, that is form systems with initial donor mass $`1.5M_{}`$. If this is true, and if we assume that in such a system we have equal probability that the companion transfers any amount of mass between $`0.5M_{}`$ and $`2M_{}`$ on to a $`1.4M_{}`$ NS, for NS governed by EOS N we have a $`100\%`$ probability of obtaining a submillisecond pulsar, while for NS governed by EOS FPS we have only a $`20\%`$ probability of obtaining a spinning down millisecond pulsar, a $`9\%`$ probability of the formation of a spinning up submillisecond pulsar (doomed to gravitational collapse) and a $`71\%`$ probability of direct collapse into a black hole during the accretion phase. Therefore if the EOS governing NS matter is soft, a conservative mass transfer is more likely to end with a direct accretion-induced collapse to a black hole than with the formation of a submillisecond radio pulsar, and thus submillisecond pulsars could be hard to detect because of their low formation probability.

On the other hand, if the r-modes are excited, the spin period will remain well above one millisecond for all of the evolution. However, making the same assumptions as before, the probability of forming a pulsar drops to $`10\%`$ for EOS FPS, while for EOS N we still have a very high probability, $`90\%`$. Therefore this scenario, in which the mass transfer is conservative but the spin frequency is limited by the emission of gravitational waves, implies that the EOS is stiff in order to have an high probability of formation of millisecond pulsars.

Teherefore we should predict that the EOS of NSs is very stiff in order to explain the observational evidence (MSP are formed), if gravitational waves are indeed emitted due to r-modes excitation, while we should predict that the EOS of NSs is soft if they are not emitted, so that submillisecond pulsars are very uncommon, as the observations seem to indicate. If future observations will allow to constrain the stiffness of the NS EOS on an observational basis, this will give an indication on wherther r-modes are indeed excited in LMXBs or not.

We have shown that, if r-modes are not excited in LMXBs, the accretion process can leave us with an extremely supramassive NS, that will spin up during all of its life as a radio pulsar (case b). It is evident that, being the critical baryonic mass for getting to the extremely supramassive regime $`M_{\mathrm{crit}}`$ an EOS-dependent feature, in principle the observation an accelerating (or braking) submillisecond pulsar can allow to exclude several EOS on an observational basis.

To clarify how such an effect can help to constrain the EOS of ultradense matter we need to introduce the new concept of a critical spin period $`P_{\mathrm{crit}}`$ that, together with the minimum period ($`P_{\mathrm{min}}=2\pi /\omega _{\mathrm{max}}`$), is peculiar to each EOS. $`P_{\mathrm{crit}}`$ is the period below which the EOS allows only extremely supramassive stable configurations. In figure 5 we show sequences of equilibrium configurations with constant spin period, together with the critical baryonic mass sequence. It is evident from the figure that $`P_{\mathrm{crit}}`$ is equal to the minimum allowed period to avoid gravitational collapse if the star has $`M_B=M_{\mathrm{crit}}.`$ In fact, any constant period sequence $`P<P_{\mathrm{crit}}`$ will only include stars of baryonic mass greater than the critical one. Thus any NS with $`P<P_{\mathrm{crit}}`$ will accelerate as a consequence of energy loss due to magnetic dipole radiation. Being $`P_{\mathrm{crit}}`$ EOS-dependent, the detection of a submillisecond radio pulsar and the determination of the sign of its period derivative will allow to effectively constrain the equation of state governing ultradense matter.

Thus the detection of a submillisecond radio pulsar can impose two constraints on the EOS of the NS:

1. the spin period must be larger than the minimum allowed period, i.e. the spin period of the maximum rotation configuration, $`P_{\mathrm{min}}`$.
2. if the period is shorter than $`P_{\mathrm{crit}}`$, the radio pulsar must spin up rather than spin down.

Both $`P_{\mathrm{min}}`$ and $`P_{\mathrm{crit}}`$ are EOS dependent and are longer for stiffer EOSs. In fact, the detection of a submillisecond radio pulsar with spin period $`P_{\mathrm{obs}}`$ undergoing a spin up will rule out all the stiff EOSs with $`P_{\mathrm{min}}>P_{\mathrm{obs}}`$. On the other hand the detection of a spinning down submillisecond radio pulsar, with spin period $`P_{\mathrm{obs}}`$, will allow us to rule out all the stiff EOSs with $`P_{\mathrm{crit}}>P_{\mathrm{obs}}`$, because they cannot explain a spinning down radio pulsar with such a short spin period. In this case the limit is more stringent because $`P_{\mathrm{crit}}>P_{\mathrm{min}}`$! As an example, suppose that a spinning down radio pulsar with a period of $`0.713`$ ms (like the one we obtain in case a) will be detected: this will allow us to rule out EOS N, since although the minimum period for this EOS is $`0.69`$ ms, any radio pulsar governed by EOS N with such a low spin period will spin up, being for EOS N $`M_{\mathrm{crit}}=3.63M_{}`$ and $`P_{\mathrm{crit}}=0.74\mathrm{ms}`$.

In summary, in this paper we presented the first results obtained with a new code that allows to study in details the binary system evolution and the spin evolution of the NS, on the basis of fully relativistic calculations. We used this code to study the evolution of systems with conservative mass transfers and confirmed that the large amount of matter that is transferred on to the NS will spin it up to periods well below one millisecond, unless the emission of gravitational waves dissipates the excess of angular momentum. However in this last case the amount of mass acreted onto the NS is easily big enough to cause a direct collapse to a black hole. Therefore we concluded that presumably the recycled systems which give origin to the MSP observed to date should have origin from systems with a highly nonconservative mass transfer. We showed that if the EOS of ultradense matter is not very stiff the direct collapse to a black hole is very likely to happen even if the r-modes are not excited. This could explain the lack of any observation of submillisecond radio pulsars even without invoking gravitational waves emission. As a last remark, since we showed that there is the possibility of obtaining from binary evolution some unusual, accelerating submillisecond radio pulsars, we introduced the new concept of the critical spin period $`P_{crit}`$, peculiar to each EOS, that can allow to effectively constrain the EOS of NS matter if a radio pulsar with a period below one millisecond will be observed in the future.

## Acknowledgements

We thank the anonymous referee for useful discussions that helped us improve this paper. This work was partially supported by the Ministero dell’Istruzione, dell’Università e della Ricerca (MIUR) with PRIN 2001 and by the Università di Palermo.


text:
# Do medium range ensemble forecasts give useful predictions of temporal correlations?

## 1 Introduction

We consider the question of how to make probabilistic forecasts of future temperatures over the 1-10 day timescale. A complete specification of the distribution of future temperatures over this time period would consist of information about the marginal distributions of temperature on each day and about the dependencies between the temperatures on different days. Investigating these marginal distributions and the dependency structure in full generality is rather difficult, and so we will make the approximation that temperature is normally distributed. This simplifies the problem greatly since in this case the distribution of future temperatures is described completely by 10 means, 10 variances, and a 10 by 10 correlation matrix.

Ensemble forecasts can be used to derive predictions of the mean, and such predictions are better than predictions derived from single integrations of forecast models at the same resolution. Ensemble forecasts can also be used to derive predictions of the variance. There are a number of ways this can be done. We have analysed some of these in detail in previous articles, including using linear regression on the ensemble mean, and spread regression on the ensemble mean and spread (see Jewson et al. (2003), Jewson (2003a), Jewson (2003b) and Jewson (2003c)). We have only been able to show that the spread of the ensemble improves the skill of forecasts when the calibrated forecasts are evaluated in-sample (i.e. the calibration and evaluation are performed on the same data that is used to calculate the calibration parameters). For out of sample forecasts it seems that it is very hard to prove that one can beat linear regression on the ensemble mean as a calibration method. We argue that this is because the predictable part of the variability in the variance is small and because only short records of past forecasts are available for training calibration models.

We now address the question of whether ensemble forecasts can be used to predict the correlations between temperatures on different days of the forecast, or whether such predictions should be made using past-forecast error statistics. The potential advantage of deriving the prediction of correlations from the ensemble is that it will then be flow-dependent. It is highly plausible that as the atmospheric state changes the correlations between forecast errors on different days of the forecast should change, and ensembles have the potential to predict that effect whereas past forecast error statistics do not. On the other hand numerical model forecasts are very prone to biases and our experience with trying to predict the variance of forecast errors from the ensemble has taught us that nothing from an ensemble forecast should be taken for granted: everything needs careful analysis and calibration to extract the useful information.

The author became interested in the question of whether ensembles can predict temporal correlations because it arises in a simple weather derivative pricing situation, as described in Jewson and Caballero (2002). Consider a weather option based on December mean temperature. To calculate the fair value of such an option one has to estimate the distribution of the settlement index i.e. the distribution of December mean temperatures. When it is estimated many months prior to the start of the contract, this distribution can be derived entirely from historical data, while when it is estimated immediately prior to or during the contract it should be derived from a combination of historical data and forecasts. For example, imagine that we are estimating this distribution on the 1st of December. The expectation of the settlement index is then given by the sum of two values. The first of these values is the contribution from Dec 1st to Dec 10th (which can be estimated from a forecast) and the second is the contribution from Dec 11th to Dec 31st (which can be estimated from historical data). Thus estimating the mean of the settlement index is rather straightforward. Estimating the variance of the index, however, is more complicated, and involves making estimates of the variances of temperatures on each day of the month (31 values) and the correlations between temperatures on different days of the month (a 31x31 element matrix). One part of estimating this 31x31 matrix is to estimate the correlations between the days of the forecast (a 10x10 element sub-matrix) and this is what motivates the question of whether those correlations can be derived from the ensemble. A more general discussion of how the rest of the correlations and variances in this problem can be calculated is given by Jewson and Caballero (2002).

## 2 Data

We will base our analyses on one year of ensemble forecast data for the weather station at London’s Heathrow airport, WMO number 03772. The forecasts are predictions of the daily average temperature, and the target days of the forecasts run from 1st January 2002 to 31st December 2002. The forecast was produced from the ECMWF model (Molteni et al., 1996) and downscaled to the airport location using a simple interpolation routine prior to our analysis. There are 51 members in the ensemble. We will compare these forecasts to the quality controlled climate values of daily average temperature for the same location as reported by the UKMO.

There is no guarantee that the forecast system was held constant throughout this period, and as a result there is no guarantee that the forecasts are in any sense stationary, quite apart from issues of seasonality. This is clearly far from ideal with respect to our attempts to build statistical interpretation models on past forecast data but is, however, unavoidable: this is the data we have to work with.

Throughout this paper all equations and all values are in terms of double anomalies (have had both the seasonal mean and the seasonal standard deviation removed). Removing the seasonal standard deviation removes most of the seasonality in the forecast error statistics, and partly justifies the use of non-seasonal parameters in the statistical models for temperature that we propose.

## 3 Models

There are potentially many ways that one could address the question of whether or not the ensemble can be used to predict temporal correlations. Since this question has not, apparently, been addressed before, we will take a simple and pragmatic approach, which works as follows.

We will model the mean and variance of the forecast using the spread regression model of Jewson et al. (2003) and will perform in-sample calibration of all the forecasts using this model. We are happy to perform this mean-variance calibration entirely in-sample because we are addressing the question of whether the *correlations* contain useful information, rather than whether the mean and variance of the ensemble do. We note that this calibration does not affect the correlations between days.

We will model the correlation matrix between the days of the forecast as a weighted sum of the two matrices $`C^{\text{past forecast error statistics}}`$ and $`C^{\text{ensemble forecast}}`$ as follows:

$$C_i=\lambda C^{\text{past forecast error statistics}}+(1\lambda )C_i^{\text{ensemble forecast}}$$
(1)
where $`C_i`$ is the modelled correlation matrix on day $`i`$, $`C^{\text{past forecast error statistics}}`$ is a stationary matrix derived from past forecast error statistics, and $`C_i^{\text{ensemble forecast}}`$ is a time-varying matrix derived from the ensemble forecast. We will vary the weighting $`\lambda `$ of these matrices from zero to one and derive the combination that gives the optimum probabilistic forecast, defined as that forecast which maximises the likelihood.

Our first test is an in-sample test that fits the past forecast error based correlation matrix on the whole year of data, and then tests it on the same year of data. One can argue that this test is not very useful, since we are fitting a very large number of parameters (the 56 independent elements of the correlation matrix) on a relatively small amount of data. There is a very large danger of over-fitting.

Our second test avoids this problem by fitting the past forecast error based correlation matrix on the first six months of data and testing it on the second six months of data, and vice versa.

The likelihood score that we will attempt to maximise by combining the two correlation matrices is given by the multivariate normal distribution over all forecast days and all leads. We will assume that forecasts are independent from day to day but not from lead to lead. The likelihood then becomes:

$$L=\underset{i=1}{\overset{i=N}{}}\frac{1}{(2\pi )^{\frac{n}{2}}D_i^{\frac{1}{2}}}\text{exp}\left(\frac{1}{2}(e_i^TC_i^1e_i)\right)$$
(2)
where $`N`$ is the number of days of forecasts, $`n`$ is the number of forecast leads, $`e_i`$ is a vector of forecast errors on day $`i`$ (of length $`n`$), $`C_i`$ is the estimated forecast error covariance matrix on day $`i`$ (of dimension $`n`$ by $`n`$) as given by equation 1, and $`D_i`$ is the determinant of this matrix.

## 4 Results

Our in sample results are shown in figure 1. The horizontal axis shows the weight applied to the past forecast error based correlation matrix (the $`\lambda `$ in equation 1). We have plotted the log-likelihood, given by the log of the likelihood from equation 2. As we vary the correlation matrix, the likelihood changes. We see that the highest values for the likelihood are given when we weight the two correlation matrices roughly in the proportions 90% (for the past forecast error based matrix) to 10% (for the ensemble forecast based matrix). In spite of the caveats we have about the in-sample nature of this test, the results are somewhat interesting in that they certainly imply that the ensemble forecast based correlation matrix contains some useful information. Because it has been performed in sample we would expect this test to strongly favour the past forecast error based matrix, and hence to be biased towards higher values of $`\lambda `$ than out of sample tests.

The out of sample results, in which the data being predicted is different from the data used to calculate the past forecast error correlation matrix, are shown in figures 2 and 3. We see that in both cases the likelihood has a maximum at around 80% (in fact, the exact numbers are 77% and 78%). As we expected, the optimum combination is at a lower level for the weight than for the in-sample tests. The results for the two tests are remarkably consistent, giving us reasonable faith that sampling error is not playing too important a role.

If we were forced to use either the past forecast error based correlation matrix or the ensemble based correlation matrix, then we see clearly that the past forecast error based correlation matrix performs better i.e. gives higher values for the log-likelihood. However there is a wide range of linear combinations of the two matrices that performs better still.

In figure 4 we show an example of the correlation time series derived from the optimum values of lambda. This example is based on data from the first 50 days of the predicted data from the first of the out of sample tests, and shows the correlation between leads 2 and 3. The two solid lines show the correlations from the past forecast errors and the ensemble forecast, while the dotted line shows the weighted combination of these two correlations using the optimum value for lambda. We see that the ensemble based correlation is on average lower than the correlation based on past forecast errors, but that for some values it is higher. The optimum weighting of the two reduces the variability of the predicted correlation very dramatically.

## 5 Summary

We have investigated whether the temporal correlations derived from an ensemble forecast are useful predictors for the correlations between the distributions of possible temperatures on each day of the forecast. We find:

* when we compare probabilistic forecasts derived from a correlation matrix based on past forecast error statistics with forecasts derived from a correlation matrix based on the ensemble, those based on past forecast error statistics are better
* a linear combination of the two correlation matrices performs better than either matrix on its own
* the optimum proportions of the two matrices seem to be fairly robustly given by around 77% for the past forecast error correlation matrix and around 23% for the ensemble derived correlation matrix

We have deliberately used a very simple methodology. In particular, we have used just a single weight for the whole correlation matrix. One could also consider using different weights for the different elements of the matrix. This might give better results, since one could certainly imagine that the ensemble would give better correlations at shorter leads, while the past forecast error based correlations would perform better at longer leads. If having separate weights for each member of the matrix proves unwieldy (there would be 56 different weights) one could consider parametrising the structure of the weights with a smaller number of parameters.

In terms of further work, there are two other questions related to the predictability of correlations that immediately spring to mind. The first is to ask: do ensemble forecasts contain useful correlation information about the correlation *between* stations? One could look at both the instantaneous cross-correlations, and the lag cross-correlations. The second is: do ensembles contain useful information about the correlation between different variables i.e. between temperature and precipitation at the same location and the same lead time, between temperature and precipitation at a different location but the same lead time, or even between temperature and precipitation at a different location and at a different lead time.

## 6 Acknowledgements

Thanks to Ken Mylne for providing the data and Rie Kondo for her patience.