text:
# The First Substellar Subdwarf? Discovery of a Metal-poor L Dwarf with Halo Kinematics

## 1 Introduction

Subdwarfs are metal-deficient stars, classically defined as lying below the stellar main-sequence in optical color-magnitude diagrams (Kuiper, 1939). These objects are in fact not subluminous but rather hotter (i.e., bluer in optical colors) than equivalent mass main-sequence dwarfs, a consequence of their reduced metal opacity (Chamberlin & Aller, 1951; Sandage & Eggen, 1959). Cool subdwarfs (spectral types sdK and sdM) are typically found to have halo kinematics ($`V=202`$ km s-1; Gizis 1997), and these objects are presumably relics of the early Galaxy, with ages of 10–15 Gyr. Because low mass subdwarfs have lifetimes far in excess of the age of the Galaxy, they are important tracers of Galactic chemical history, and are representatives of the first generations of star formation.

All of the coolest subdwarfs (\[Fe/H\] $`1.2\pm `$0.3) and extreme subdwarfs (\[Fe/H\] $`2.0\pm `$0.5; Gizis 1997) currently known have been identified in optical proper motion surveys, most notably Luyten’s Half-Second Catalog (Luyten, 1979a, LHS) and Two-Tenths Catalog (Luyten, 1979b), the APM Proper Motion Survey (Scholtz et al., 2000), and the Galactic plane survey of Lepine, Shara, & Rich (2002). The substantial space velocities of halo subdwarfs allow them to stand out in these surveys amongst the overwhelming multitude of similarly faint, but slowly moving, background stars. However, only a handful of very cool subdwarfs with spectral types sdM6/esdM6 and later111Based on the Gizis (1997) classification scheme. have been found (Gizis, 1997; Gizis & Reid, 1997; Schweitzer et al., 1999; Lepine, Rich, & Shara, 2003; Lepine, Shara, & Rich, 2003). These stars exhibit characteristic spectral signatures of strong metal hydrides (CaH, MgH, AlH, FeH), weak and/or absent metal oxides (TiO, VO, CO), and enhanced collision-induced (CIA) H2 absorption (Mould & Hyland, 1976; Liebert & Probst, 1987; Saumon et al., 1994; Leggett et al., 2000). With T$`{}_{eff}{}^{}`$ 2900 K (Leggett et al., 2000), the coolest subdwarfs known have masses just above the Hydrogen Burning Minimum Mass (HBMM)222The HBMM is defined here as the minimum mass for which core Hydrogen fusion balances luminosity at ages later than 10 Gyr., which ranges from 0.072 M for Solar composition ($`Z=Z_{\mathrm{}}`$) to 0.092 M for $`Z=0`$ (Chabrier & Baraffe, 1997; Burrows et al., 2001).

In contrast, hundreds of significantly cooler Population I, or disk, dwarfs have been identified333A listing of known M, L, and T dwarfs is maintained by J. D. Kirkpatrick at http://spider.ipac.caltech.edu/staff/davy/ARCHIVE/index.html., predominately in the wide-field optical and near-infrared (NIR) surveys DENIS (Epchtein et al., 1997), 2MASS (Skrutskie et al., 1997), and SDSS (York et al., 2000). These discoveries include many dozens of substellar objects, and extend well beyond the M spectral class into the L (Kirkpatrick et al., 1999; Martín et al., 1999) and the T (Burgasser et al., 2002; Geballe et al., 2002) classes. The success of recent surveys in identifying very cool dwarf star and brown dwarfs derives from their use of red-sensitive CCDs and infrared array detectors. As we proceed to cooler temperatures, the majority of emitted flux from both disk dwarfs and subdwarfs moves to NIR wavelengths. Unfortunately, proper motion surveys of sufficient temporal breadth to identify halo stars rely primarily on optical, typically photographic, imaging. Factoring in their intrinsic rarity (0.3% number density compared to disk dwarfs; Reid & Hawley 2000) and substantial evolution toward lower luminosities (Burrows et al., 1993; Chabrier & Baraffe, 1997, see Figure 3), it is clear why few very cool subdwarfs have been found.

In this article, we present the discovery of a new cool subdwarf, identified in the 2MASS database. This object, 2MASS J05325346+8246465444Source designations for the 2MASS Point Source Catalog are given as “2MASS Jhhmmss\[.\]ss$`\pm `$ddmmss\[.\]s”. The suffix conforms to IAU nomenclature convention and is the sexigesimal R. A. and decl. at J2000 equinox. We adopt a shorthand notation of “2MASS hhmm$`\pm `$ddhh” throughout the remainder of this Letter., is a high motion source with a spectral morphology consistent with a metal deficient L dwarf. In $`\mathrm{\S }`$ 2 we discuss the identification of 2MASS 0532+8246 and subsequent NIR and optical observations obtained at the Palomar 60” Telescope. In $`\mathrm{\S }`$ 3 we describe spectroscopic observations of 2MASS 0532+8246 obtained with the Near-IR Spectrometer (McLean et al., 1998, 2000, hereafter NIRSPEC) and Low Resolution Imaging Spectrograph (Oke et al., 1995, hereafter LRIS) at Keck Observatory. In $`\mathrm{\S }`$ 4 we derive the space motion of this object from the imaging and spectroscopic data. In $`\mathrm{\S }`$ 5 we discuss the substellarity of 2MASS 0532+8246 based on subsolar metallicity models from Burrows et al. (2001). We discuss our results in $`\mathrm{\S }`$ 6.

## 2 Identification of 2MASS 0532+8246

2MASS 0532+8246 was identified in a search for T dwarfs using the 2MASS Working database (Burgasser et al., 2002). It was selected for its blue NIR colors ($`JK_s=0.26\pm 0.16`$) and lack of an optical counterpart at the 2MASS position in both POSS-I and POSS-II photographic plates. A finder chart is given in Figure 1, and astrometric and photometric data from 2MASS are listed in Table 1. Follow-up NIR imaging using the Palomar 60” Infrared Camera (Murphy et al., 1995, hereafter IRCam) on 25 September 1999 (UT) determined that this object is not an uncatalogued minor planet. Gunn $`r`$\- and $`z`$-band images obtained with the Palomar 60” Facility CCD Camera on 18 November 2002 and 25 January 2003 (UT) further confirmed the red optical/NIR colors of this object. Aperture photometry from 18 November 2002, calibrated with observations of the SDSS photometric standard BD +33$`\mathrm{°}`$ 4737 (Smith et al., 2002, correcting ($`r^{},z^{}`$) $``$ (gunn $`r`$, gunn $`z`$) using photometry from Kent 1985) yield $`r>20.9`$ (5$`\sigma `$ upper limit) and $`z=18.30\pm 0.16`$. As the NIR and optical/NIR colors of this object are consistent with those of other identified late-type L and T dwarfs (Kirkpatrick et al., 1999; Fan et al., 2000; Dahn et al., 2002), 2MASS 0532+8246 was initially believed to be a firm T dwarf candidate.

## 3 Spectroscopic Observations

### 3.1 NIR Spectrum

We obtained NIR spectra of 2MASS 0532+8246 on 24 December 2002 (UT) using NIRSPEC on the Keck-II 10m Telescope. Conditions were clear with seeing of 0$`\stackrel{}{\mathrm{.}}`$6. We observed the object in three spectral configurations, N3 (1.14 – 1.38 $`\mathrm{\mu m}`$), N5 (1.51 – 1.79 $`\mathrm{\mu m}`$), and N6b (1.94 – 2.32 $`\mathrm{\mu m}`$), corresponding roughly to the J, H, and K NIR bands. Spectral resolution was R $``$ 2000 for our 0$`\stackrel{}{\mathrm{.}}`$38 slit. For each observation, two exposures of 300s each were obtained, dithering 20$`\mathrm{}`$ along the slit. All observations were made in the airmass range 2.6–3.0. Immediately after each sequence of exposures, internal flat field and NeAr arc lamps were observed for calibration, and the A0 star HD 34360 was observed just before or after the target source in each spectral order at similar airmasses. Data were reduced using the REDSPEC package555See http://www2.keck.hawaii.edu/inst/nirspec/redspec/index.html.; a detailed description of our reduction procedures is given in McLean et al. (2003). Final flux calibration of the individual spectral orders was done using 2MASS photometry; however, because the orders do not span the entire photometric bandpasses (Cutri et al., 2003), our final flux calibration may be somewhat underestimated at J-band, but reasonably correct for the H- and K-bands.

The resulting NIR spectrum is shown in Figure 2a. The 1.5–2.4 $`\mathrm{\mu m}`$ spectral energy distribution is fairly smooth and blue, reminiscent of an early-type main-sequence star or white dwarf, and quite unlike the much redder NIR spectra of late-type M and L dwarfs (cf. overlay of the L7 DENIS 0205$``$1159AB in Figure 2a; Delfosse et al. 1997; Kirkpatrick et al. 1999; McLean et al. 2003). However, the red optical/NIR colors and optical molecular bands (see below) of 2MASS 0532+8246 are inconsistent with a hot stellar atmosphere. We instead attribute the shape of the NIR spectrum to enhanced absorption by the CIA H2 1-0 quadrupole band centered near 2.5 $`\mathrm{\mu m}`$ (Saumon et al., 1994; Borysow, Jørgensen, & Zheng, 1997). This pressure-sensitive band, present in late-type L and T dwarfs (Tokunaga & Kobayashi, 1999; Burgasser et al., 2002) and cool subdwarfs (Leggett et al., 2000), is generally broad and featureless, consistent with the observed spectrum. Strong H2 absorption may explain the absence of CH4 absorption bands at 1.6 and 2.2 $`\mathrm{\mu m}`$, the defining features of T dwarfs; and the 2.3 $`\mathrm{\mu m}`$ CO band, a key spectral feature in M, L, and early-type T dwarf spectra. On the other hand, H2O absorption appears to be present at 1.3–1.5 and 1.7–2.0 $`\mathrm{\mu m}`$, and minor features throughout the H-band are likely attributable to weaker H2O lines. Two weak features at 1.573 and 1.626 $`\mathrm{\mu m}`$ are present, possibly attributable to FeH absorption (Cushing et al., 2003) or poor correction of telluric OH lines; these features are not coincident with the 1.6 $`\mathrm{\mu m}`$ CH4 band. Overall, this portion of the spectrum of 2MASS 0532+8246 is quite unlike any known late-type M, L, or T dwarf, based on its blue slope and absence of CH4 and CO bands.

The J-band spectrum, shown in detail in Figure 2b, is far more rich, with atomic lines of K I (1.1690, 1.1773, 1.2543, 1.2522 $`\mathrm{\mu m}`$) and Fe I (1.1528, 1.1886, 1.1976, 1.2128, 1.3210 $`\mathrm{\mu m}`$), and strong molecular bands of FeH (1.194, 1.239 $`\mathrm{\mu m}`$) and H2O (1.3 $`\mathrm{\mu m}`$) observed. Above the spectrum we plot an opacity spectrum of FeH (Dulick et al., 2003), which shows that much of the fine structure observed is caused by this molecule. There are no lines of Mg I, Ca I, or Al I present in this spectral region, and the $`\varphi `$ TiO system around 1.25 $`\mathrm{\mu m}`$ (Jørgensen, 1994) is either absent or obscured by FeH absorption. The overall spectral morphology at J-band is quite similar to the L7 DENIS 0205$``$1159AB, with the notable exception of stronger FeH and H2O absorption.

### 3.2 Optical Spectrum

We obtained a red optical spectrum of 2MASS 0532+8246 on 3 January 2003 (UT) using LRIS on the Keck-I 10m Telescope. Conditions were clear with sub-arcsecond seeing. Two exposures of 1200 sec each were obtained through the red channel, dithered by 2$`\mathrm{}`$ between exposures along the 1$`\mathrm{}`$ (4.7 pixels) slit. We employed the 400 lines mm-1 grating blazed at 8500 Å, yielding 6300–10100 Å spectra with 7 Å resolution (R $``$ 1200). Dispersion on the chip was 1.9 Å pixel-1. The OG570 order-blocking filter was used to suppress higher-order light. Observations of the B1 V flux standard Hiltner 600 (Hamuy et al., 1994) were obtained for flux calibration, and the G0 V star HD 38847 was observed after the target at similar airmass (2.4) for telluric calibration. HeNeAr arc lamp exposures were taken immediately after the target observations for wavelength calibration, and quartz lamp flat-field exposures (reflected off of the interior dome) were observed at the start of the night to correct for detector response. Data reduction procedures, using standard IRAF666Image Reduction and Analysis Facility (IRAF) is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. routines, are discussed in detail in Kirkpatrick et al. (2003).

The final flux-calibrated and telluric-corrected optical spectrum for 2MASS 0532+8246 is included in Figure 2a and shown in detail in Figure 2c. Strong metal hydride bands of FeH (8692, 9896 Å), CrH (8611, 9969 Å), and CaH (6750 Å) are clearly discernable. The 9896 Å FeH Wing-Ford band is the strongest seen in any cool spectrum to date, and absorption from the shorter wavelength FeH and CrH bands extends to roughly 9100 Å. Weaker bands of TiO (7053, 8432 Å) and H2O (9250, 9400 Å) are also seen. We point out an absorption feature between 9570–9700 Å which is too strong to be attributable to intrinsic H2O absorption, based on the weakness of the 9250 Å band; it is also not properly placed for FeH, CrH, TiH, or CH4 absorption. This band is weakly present in late-type M and early-type L dwarf spectra, and may be a hitherto unrecognized metal hydride band. It is also possible, although unlikely, that this feature is attributable to residual telluric absorption, as telluric corrections were made only to 9650 Å (the extent of the strong terrestrial H2O bands; Stevenson 1994). Individual alkali lines of Na I (8183, 8195 Å), Rb I (7800, 7948 Å), and Cs I (8521, 8943 Å) are present, while the resonance K I doublet (7665, 7699 Å) is exceptionally strong and pressure-broadened, as seen in late-type L and T dwarfs (Kirkpatrick et al., 1999; Liebert et al., 2000). The red optical/NIR colors of 2MASS 0532+8246 are clearly a consequence of the K I red wing; integrating the spectrum over the $`I_C`$ bandpass (Bessell, 1990) yields $`I_C`$ $``$ 19.2 and $`I_CJ`$ $``$ 4.0, similar to mid- and late-type L dwarfs (Dahn et al., 2002). The blue end of the spectrum appears to be similarly suppressed by the pressure-broadened Na I doublet (5890, 5896 Å). Both Li I (6708 Å) absorption and H$`\alpha `$ (6563 Å) emission lines are absent. The overall spectral morphology, excluding hydride bands and TiO, is again quite similar to the L7 DENIS 0205$``$1159AB, suggesting a rather cool atmosphere.

### 3.3 Characterization and Spectral Classification

2MASS 0532+8246 is not a T dwarf, based on the absence of the 1.6 and 2.2 $`\mathrm{\mu m}`$ CH4 bands (Burgasser et al., 2002; Geballe et al., 2002), and the presence of strong FeH and CrH and weak TiO and Na I features in the optical (Burgasser, Kirkpatrick, & Liebert, 2003). It appears that this object is instead a metal-deficient, late-type L dwarf. Enhanced CIA H2 and strong hydride bands in the optical and NIR spectra of 2MASS 0532+8246 are consistent with a metal-poor atmosphere. Both are consequences of the overall reduced metal opacity, resulting in larger column abundances of the remaining chemical species and increased relative abundance of metal hydrides over double-metal species. The presence of the 7053 and 8432 Å TiO bands is somewhat at odds with this interpretation, particularly since TiO is generally weak or absent in disk dwarfs later than L5 (Kirkpatrick et al., 1999). One possibility is that the mechanism of Ti depletion in these objects – incorporation into multiple-metal condensates of CaTiO3, Ca3Ti2O7, Ca4Ti3O10, and Ti2O3 (Lodders, 2002) – may be inhibited in the low-metallicity atmosphere of 2MASS 0532+8246. Chemical equilibrium calculations incorporating subsolar abundances are required to address this hypothesis. The TiO bands may also indicate a somewhat warmer atmosphere than the L7 comparison object DENIS 0205$``$1159AB.

Classifying 2MASS 0532+8246 is therefore not straightforward, particularly as there are no other metal-deficient L dwarfs currently known for direct comparison. The M subdwarf and extreme subdwarf sequences of Gizis (1997) terminate at sdM7 and esdM5.5, respectively, although recent late-type discoveries are classified via extrapolation and limited modifications of this scheme (Lepine, Rich, & Shara, 2003; Lepine, Shara, & Rich, 2003). Metal classification (d, sd, esd) in the Gizis scheme is based on the strength of the 7053 Å TiO band, which is exceedingly weak in the spectrum of 2MASS 0523+8246. Numerical classification is based on the 6750 Å CaH band, which is suppressed in this spectrum by Na I and K I absorption. Thus, a simple extrapolation of the Gizis scheme fails in the L dwarf regime. While there is general agreement between the spectral morphologies of 2MASS 0532+8246 and the DENIS 0205$``$1159AB, the optical/J-band hydride features and TiO bands are more similar to mid-type L dwarfs (Kirkpatrick et al., 1999; McLean et al., 2003). Additionally, there are no cool dwarf standards that match the 1.3–2.4 $`\mathrm{\mu m}`$ spectral region. We therefore tentatively classify this object a late-type sdL. Identification of additional metal-poor L subdwarfs may eventually allow the definition of a robust classification scheme.

## 4 Space Motion

Examination of the IRCam and CCD followup images show that 2MASS 0532+8246 has moved significantly since it was first imaged by 2MASS. We therefore measured its proper motion by deriving its coordinates on each of the images, using the 2MASS coordinates of background stars for astrometric calibration. Results are listed in Table 2. A linear fit yields $`\mu `$ = 2$`\stackrel{}{\mathrm{.}}`$60$`\pm `$0$`\stackrel{}{\mathrm{.}}`$15 yr-1 at position angle 130.0$`\pm `$1.8$`\mathrm{°}`$, making 2MASS 0532+8246 one of the highest proper motion stars known (Bakos, Sahu, & Németh, 2002). This object is not present in the LHS because of its optical faintness, although a faint counterpart may be present in the POSS-II R-band image at the extrapolated position (Figure 1). Residuals in the proper motion fit are indicative of parallactic motion, but positional uncertainties are too large for a reliable measure.

The radial velocity of 2MASS 0532+8246 was derived by cross-correlating its J-band spectrum with those of five early- and mid-type L dwarfs (chosen for their strong NIR FeH absorption) with measured radial velocities: 2MASS 0746+2000AB, 2MASS 1439+1929, G 196-3B, Kelu 1, and DENIS 1228$``$1547AB (Basri et al., 2000; Reid et al., 2002). All five of these objects have been observed as part of the NIRSPEC Brown Dwarf Spectroscopic Survey (McLean et al., 2003), and are therefore wavelength calibrated in the same manner as the data presented here. We derive a heliocentric $`v_{rad}=195\pm 11`$ km s-1, where the uncertainty is derived from the scatter of values amongst the comparison objects. This value is consistent with heliocentric velocity offsets in the 1.1690/1.1773 $`\mathrm{\mu m}`$ K I lines ($`v_{rad}=175\pm 17`$ km s-1) and the optical Cs I and Rb I lines ($`v_{rad}=192\pm 24`$ km s-1).

The substantial radial velocity and proper motion of 2MASS 0532+8246 identifies it as a high velocity source. Determining its total space motion, however, requires an estimate of its distance. If we assume 2MASS 0532+8246 is bound to the Galaxy (Galactocentric $`v_{tot}500`$ km s-1; Carney, Latham, & Laird 1988), we derive an upper limit distance of 54 pc assuming ($`U,V,W`$) = (10.00, 5.25, 7.17) km s-1 (Dehnen & Binney, 1998) and Local Standard of Rest (LSR) velocity $`V_{LSR}=220`$ km s-1 (Kerr & Lynden-Bell, 1986). A spectrophotometric parallax estimate is highly uncertain, given the substantial redistribution of flux in the NIR and uncertainty in the numeric subtype for this object. However, assuming that the $`M_{I_C}`$ magnitude of 2MASS 0532+8246 is similar to that of mid- to late-type L dwarfs (Dahn et al., 2002, $`17M_{I_C}19`$), we derive a conservative estimate of 10–30 pc. The corresponding heliocentric space velocities at the median distance are ($`U,V,W`$) = ($``$13, $``$301, 28) km s-1. The substantial V velocity, between $``$390 and $``$212 km s-1 over the adopted distance range, indicates zero or retrograde motion with respect to the LSR, and lends additional support that this object is a member of the halo population.

## 5 Is 2MASS 0532+8246 Substellar?

Disk L dwarfs with spectral types later than $``$L4 are predominately substellar (Gizis et al., 2000); hence, the late-type spectral morphology of 2MASS 0532+8246 suggests that it too may be a brown dwarf. To test this hypothesis, we compared the estimated Teff of this object to subsolar metallicity evolutionary models from Burrows et al. (2001). In Figure 3, we plot Teff versus mass for ages 10 and 15 Gyr and metallicities $`Z`$ = 0.1, 0.3, and 0.01 $`Z_{\mathrm{}}`$. At these late ages, the transition between stellar and substellar masses spans a substantial range in temperature, and is more dramatic for lower metallicities due to the reduced atmospheric opacities and hence enhanced luminosities (Burrows et al., 1993).

Deriving an accurate Teff for 2MASS 0532+8246 is difficult without adequate distance or bolometric flux measurements. Solar-metallicity mid- to late-type L dwarfs have temperatures $`1400T_{eff}1800`$ K (Kirkpatrick et al., 1999; Basri et al., 2000; Pavlenko, Zapatero Osorio, & Rebolo, 2000; Leggett et al., 2001; Stephens et al., 2001; Burgasser et al., 2002; Dahn et al., 2002). In general there is little difference between the temperatures of M-type subdwarfs and dwarf stars of the same numeric type (Gizis, 1997; Leggett et al., 2000), although late-type M subdwarfs (e.g., the sdM7 LHS 377) may be 100-300 K hotter (Leggett et al., 2000). We therefore adopt a conservative temperature range of $`1400T_{eff}2000`$ K. We further assume that the metallicity of 2MASS 0532+8246 lies in the range $`0.01Z0.1`$ $`Z_{\mathrm{}}`$, typical for sdM and esdM stars (Gizis, 1997; Gizis & Reid, 1997). For these parameters, we derive a mass range of 0.077 $``$ M $``$ 0.085 M. Larger masses correspond to lower metallicities, so that in each case the derived mass is below the HBMM. Therefore, we conclude the 2MASS 0532+8246 is probably substellar, although more refined estimates of its atmospheric properties through further observations and spectral modelling are required to validate this result.

## 6 Discussion

The existence of a substellar subdwarf is not unexpected, as mass functions for the nearby halo population (Gizis & Reid, 1999) and globular clusters (Piotto & Zoccali, 1999) are shallow but rising: dN/dM $``$ M-(0.5-1.3) at M $``$ 0.15 M, similar to the disk population (Reid et al., 1999; Chabrier, 2002). However, the steep drop in Teff, and hence luminosity, across the metal-poor substellar boundary implies that the vast majority of these halo brown dwarfs will have T$`{}_{eff}{}^{}1300`$ K; i.e., in the T dwarf regime and cooler (Burgasser et al., 2002). Indeed, one possible T-type subdwarf has already been identified, 2MASS 0937+2931 (Burgasser et al., 2002), an object which also exhibits unusually strong CIA H2 and FeH absorption (Burgasser, Kirkpatrick, & Liebert, 2003) and a substantial proper motion (Burgasser et al., 2003). Relatively few L-type subdwarfs should exist, however, as they span a much narrower range in mass at ages later than 10 Gyr. The discovery of 2MASS 0532+8246 may therefore be quite fortuitous.

Given its optical faintness, we expect that few counterparts to this late-type L subdwarf will be identified in current proper motion surveys. NIR proper motion surveys, on the other hand, could potentially reveal more of these objects by detecting them at the peak of their spectral energy distribution. A sizeable fraction of the sky ($``$30%) has been scanned more than once by 2MASS during its three years of observations, and analysis of these data may identify new, very cool high-motion stars. A second generation 2MASS survey could greatly improve on this work, adding substantially to the short list of very cool subdwarfs currently known.

A. J. B. acknowledges useful discussions with K. Cruz and J. Wei-Chun during the preparation of the manuscript, and thanks our anonymous referee for her/his helpful criticisms. We are grateful for the assistance of our Keck Observing Assistants Joel Aycock and Gary Puniwai, and Instrument Specialist Paola Amico; and our Palomar Telescope Operators, Karl Dunscombe and Barrett “Skip” Staples. We also acknowledge the 2MASS staff for their laudable efforts on the 2MASS database, without which none of this research could have been possible. A. J. B. acknowledges support by NASA through Hubble Fellowship grant HST-HF-01137.01 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. A. S. B. acknowledges funding through NASA grants NAG5-10760 and NAG5-10629. Portions of the data presented herein were obtained at the W. M. Keck Observatory which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. This publication makes use of data from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center, funded by the National Aeronautics and Space Administration and the National Science Foundation. 2MASS data were obtain through the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. POSS-I and POSS-II images were obtained from the DSS image server maintained by the Canadian Astronomy Data Centre, which is operated by the Herzberg Institute of Astrophysics, National Research Council of Canada. The authors wish to extend special thanks to those of Hawaiian ancestry on whose sacred mountain we are privileged to be guests.


format: text

text:
# Atmospheric neutrino oscillations from upward throughgoing muon multiple scattering in MACRO

## 1 Introduction

00footnotetext:  $`a`$ Also INFN Milano, 20133 Milano, Italy
$`b`$ Also IASF/CNR, Sez. di Bologna, Italy
$`c`$ Also Department of Physics, Pennsylvania State University, University Park, PA 16801, USA
$`d`$ Also Università di Trieste and INFN, 34100 Trieste, Italy
$`e`$ Also U. Minn. Duluth Physics Dept., Duluth, MN 55812
$`f`$ Also Dept. of Physics, MIT, Cambridge, MA 02139
$`g`$ Also Intervideo Inc., Torrance CA 90505 USA
$`h`$ Also Department of Physics, SLIET,Longowal,India
$`i`$Also Dipartimento di Fisica dell’Università della Calabria, Rende (Cosenza), Italy
$`l`$ Also RPD, PINSTECH, P.O. Nilore, Islamabad, Pakistan
$`m`$ Also Institute for Space Sciences, 76900 Bucharest, Romania
$`n`$ Also Università della Basilicata, 85100 Potenza, Italy
$`o`$ Also Dip. di Energetica, Universit‘a di Roma, 00185 Roma, Italy.
$`p`$ Also Dipartimento di Scienze e Tecnologie Avanzate, Università del Piemonte Orientale, Alessandria, Italy
$`q`$ Also Resonance Photonics, Markham, Ontario, Canada
$``$ Corresponding authors: SCAPPARONE@BO.INFN.IT,
SIOLI@BO.INFN.IT

The results obtained by experiments looking for neutrinos coming from natural or artificial far sources gave the first indication of physics beyond the standard model. Neutrino flavor changing is in fact the most straightforward explanation for electron and muon neutrino disappearance and for the evidence of non electron neutrinos from the sun. In this scenario, the study of atmospheric neutrinos plays an important role to improve our understanding of the neutrino oscillation mechanism.

The MACRO detector studied  three categories of events as shown in Fig. 1: (1) upward throughgoing muons, (2) upward semicontained muons, and (3+4) upward stopping muons plus downward semicontained muons. The (3+4) categories cannot be separated experimentally, because the events have only one time measurement. However, upward stopping muons and downward semicontained muons have similar parent neutrino energy.

An atmospheric muon neutrino deficit was found in all of the three categories and an angular distribution distortion was found in category (1). Such results are explained by the neutrino oscillation hypothesis with parameters $`\mathrm{\Delta }\mathrm{m}^2=2.5\times 10^3\mathrm{eV}^2`$ and $`\mathrm{sin}^22\theta `$=1, in good agreement with the Super-Kamiokande results . A detailed study of the upward throughgoing muon angular distribution by MACRO  and by Super-Kamiokande  allowed the exclusion at the 99$`\%`$C.L. of the muon neutrino into sterile neutrino oscillation mechanism, compared to the muon neutrino into the tau neutrino.

Considering a two flavor neutrino oscillation, the probability for a $`\nu _\mu `$ to oscillate to $`\nu _\tau `$ is given by:

$$P_{\nu _\mu \nu _\tau }=\mathrm{sin}^22\theta \mathrm{sin}^2\left(\frac{1.27L_\nu \mathrm{\Delta }\mathrm{m}^2}{E_\nu }\right)$$
(1)
where $`L_\nu `$(km) is the distance between the neutrino production and interaction points and $`E_\nu `$($`\mathrm{GeV}`$) is the neutrino energy. The different categories of events quoted above have median energy from $`4\text{GeV}`$ to $`50\text{GeV}`$, which provides evidence of the dependence on neutrino energies, as required by the oscillation hypothesis. In this paper we address for the first time the estimate of the upward throughgoing muon energy by using multiple scattering. The implementation of this method to the MACRO data is discussed in Section 2. Two different analysis were performed. The first one, using the tracking system in digital mode, is described in Section 3. Section 4 shows the results obtained with the streamer tube in drift mode, where electronic readout of the hit time was used to improve the position resolution and thereby the range of muon energies that could be estimated.

## 2 The muon energy estimate

The MACRO detector was extensively described in . The detector consisted of a lower half (the lower detector) and an upper half (the “Attico”). The lower detector had 10 horizontal planes of streamer tubes separated by either crushed rock absorbers ($`60\text{g/}\mathrm{cm}^2`$) or (at its top and bottom) a plane of scintillator counters. The $`Attico`$ had 4 horizontal planes of streamer tubes separated into two groups with a plane of scintillator counters between them.

Upward throughgoing muons are mainly produced in neutrino deep inelastic scattering (DIS) in the rock below the detector. The recoil hadrons are lost and the muon energy is degraded in the propagation to the detector. Nevertheless, Monte Carlo simulations show a linear relation between the parent neutrino energy $`E_\nu `$ and the muon residual energy $`E_\mu `$ at detector level.

The momentum resolution obtainable from multiple coulomb scattering (MCS) measurements is the result of two different contributions: the number of sampling planes $`N`$ and the space resolution of the detector $`\sigma _0`$. In MACRO the number of tracking planes interleaved with rock absorbers gives $`N=8`$. The resolution of the streamer tube system used in digital mode is $`\sigma _0`$$`1\text{cm}`$ and was improved to $`\sigma _0`$$`0.3\text{cm}`$ by using the drift time. The other six horizontal tracking planes are separated by a negligible amount of material, and do not contribute to the MCS measurements.

In a tracking detector with equispaced tracking planes, separated by a given absorber and neglecting the energy losses, the characteristic momentum scale can be written as :

$$p_{MCS}(\mathrm{GeV})=\frac{0.015\mathrm{\Delta }\sqrt{\mathrm{\Delta }/\mathrm{X}_0}}{\sigma _0}$$
(2)
where $`\mathrm{\Delta }`$ is the distance between tracking planes and $`X_0`$ is the material radiation length. For p$`<`$$`p_{MCS}`$, the main limitation to the momentum reconstruction comes from the number of sampling planes while for p$`>`$$`p_{MCS}`$ the space resolution dominates. Under the conditions quoted above the relative momentum error  can be written as:

$$\frac{\sigma _\mathrm{p}}{\mathrm{p}}\frac{1}{\sqrt{2\mathrm{N}}}$$
(3)
for p$`<<`$$`p_{MCS}`$, and

$$\frac{\sigma _\mathrm{p}}{\mathrm{p}}\frac{1}{\sqrt{2\mathrm{N}}}\left(\frac{\mathrm{p}}{\mathrm{p}_{\mathrm{MCS}}}\right)^2,$$
(4)
for p$`>>`$$`p_{MCS}`$. As shown in  the 1/p distribution approaches a Gaussian only for N$``$30; therefore the momentum error for MACRO is not expected to have a Gaussian behaviour. The eight horizontal streamer tube planes interleaved with rock absorbers (planes 2 to 9, from the bottom) are equispaced and have the same space resolution, hence $`p_{MCS}^{MACRO}`$=$`2.2\text{GeV/c}`$. Since about 90$`\%`$ of upward throughgoing muons have p $`>`$ $`2.2\text{GeV/c}`$, the intrinsic chamber resolution dominates. The MCS-based momentum estimates scale therefore as the square of the space resolution $`\sigma _o`$.

The r.m.s. of the lateral displacement of a relativistic muon crossing a layer of material with depth $`X`$ is proportional to the inverse of muon momentum $`p_\mu `$:

$$\sigma _{proj}^{MCS}\frac{X}{\sqrt{3}}\frac{13.6\text{MeV}}{p_\mu \beta c}\sqrt{\frac{X}{X_o}}\left(1+0.038\mathrm{ln}\frac{X}{X_o}\right).$$
(5)
In MACRO, X$``$25$`X_o`$/$`\mathrm{cos}\mathrm{\Theta }`$, giving on the vertical $`\sigma _{proj}^{MCS}`$$`10\text{cm}`$/$`E_\mu (\mathrm{GeV})`$. Therefore, a saturation point above $`10\text{GeV}`$ requires a space resolution better than $`1\text{cm}`$. An improvement in the space resolution increases the maximum energy value (saturation point) above which the MCS method is not effective.

Consequently we decided to improve the streamer tube resolution: in the second analysis, described in Section 4, we took advantage of the QTP-TDC electronics, developed for magnetic monopole searches . In this case, we obtained an improved space resolution, $`\sigma _x`$$``$$`0.3\text{cm}`$, which allowed us to estimate the muon energy up to $`E_\mu `$$``$$`40\text{GeV}`$.

For both analyses, we used the whole sample of upward throughgoing muon events collected with the complete apparatus in a period of data taking equivalent to 5.5 years of live time. We studied upward muon events selected by the time-of-flight measured by planes of scintillators combined with the standard MACRO tracking algorithm. The original upward throughgoing muon data set has been described previously . To make a comparison between real data and expectation we performed a Monte Carlo simulation using the Bartol neutrino flux  and the GRV94 DIS parton distributions  for deep inelastic scattering. For low energy neutrino interactions we used the cross sections given in . The energy loss for muons propagating through rock is taken from Lohmann et al.  while the muon simulation inside the detector was performed with GMACRO (the GEANT 3.21 based detector simulation). For the second analysis, where the streamer system is used in drift mode, another simulation chain has been implemented. In this case neutrino interactions are randomly distributed in a rock semi-sphere below the detector. Muons are then transported to the detector using the FLUKA99 package. A Monte Carlo statistics corresponding to a live time of $``$2700 years was produced (500 equivalent experiments). This simulation was compared with that used in, obtaining a satisfactory result.

## 3 The first analysis: streamer tubes in digital mode

Because of MCS in the lower detector, with this analysis we expect a measurable deflection between the incoming and the outgoing track directions for muons with energies smaller than $``$$`10\text{GeV}`$. This corresponds to an effective arm-lever of $``$$`4\text{m}`$ between the lower detector and the $`Attico`$.

The eight lowest streamer tube planes were used, through a track refit, to estimate the direction of the incoming muon. The five upper streamer tube planes in the lower detector and the four $`Attico`$ planes were used to estimate the direction of the outgoing muon. The distance $`r_w`$ between the intercepts of the two tracks in the $`z=0`$ plane, and the difference $`\mathrm{\Delta }\mathrm{\Phi }`$ between the two slopes depend on the muon energy $`E_\mu `$. We divided the upward throughgoing muons in three sub samples, according to the values of $`r_w`$ and $`\mathrm{\Delta }\mathrm{\Phi }`$: sample L= Low energy, if $`r_w>`$$`3\text{cm}`$ and $`\mathrm{\Delta }\mathrm{\Phi }>0.3^{}`$; H= High energy if $`r_w`$$`3\text{cm}`$ and $`\mathrm{\Delta }\mathrm{\Phi }0.3^{}`$. The remaining events were classified as M=Medium Energy. These cuts were optimized using two large samples of real atmospheric muons: $`i)`$ downward throughgoing (average energy of $``$ $`300\text{GeV}`$ ), and $`ii)`$ downward going stopping muons (average energy of $``$$`1\text{GeV}`$).

The aforementioned cuts were applied on the upward throughgoing muons (crossing the whole apparatus, lower detector \+ Attico), both real and simulated. They were selected from the same analysis chain  and had the same data format. From the simulation we expect 430 events: 178 of Low energy, with $`<`$$`E_\nu `$$`>=`$$`11\text{GeV}`$, 59 of Medium energy, with $`<`$$`E_\nu `$$`>=`$$`33\text{GeV}`$ and 193 of High energy, with $`<`$$`E_\nu `$$`>=`$$`72\text{GeV}`$. From the real data, 316 events were selected: 101 of Low, 51 of Medium and 164 of High energy.

The $`L`$ and $`H`$ events were further divided according to their zenith angle $`\mathrm{\Theta }`$: events from the vertical direction $`1.<\mathrm{cos}\mathrm{\Theta }<0.8`$, and events with $`\mathrm{cos}\mathrm{\Theta }>0.8`$. For each event topology the number of detected and expected events were determined and ordered with decreasing value of the $`<`$$`L_\nu `$$`>`$/ $`<`$$`E_\nu `$$`>`$ ($`<`$$`L_\nu `$$`>`$= average value of neutrino path length). The relative systematic uncertainty on each of the five $`<`$$`L_\nu `$$`>`$/$`<`$$`E_\nu `$$`>`$ values is 12.5$`\%`$, coming from the uncertainties on the neutrino spectrum and angular shape (discussed in Sect. 4.3), from the detector related effects and analysis cuts. This includes the number of planes used for the refit, the definition of the cuts, the fluctuations in the streamer tube and scintillator efficiencies, and detector acceptance uncertainties.

We evaluated $`\chi ^2`$ for the hypothesis of no neutrino oscillations using the five $`<`$$`L_\nu `$$`>`$/$`<`$$`E_\nu `$$`>`$ bins, plus the additional point from the analysis of semi contained upward going events (IU)  (upgoing neutrinos with $`<`$$`E_\nu `$$`>`$$``$$`4\text{GeV}`$). We found that the distribution agrees with the no oscillation hypothesis with a probability lower than 2%). For two-flavor oscillations with parameters $`\mathrm{\Delta }\mathrm{m}^2=2.5\times 10^3`$ e$`V^2`$ and $`\mathrm{sin}^22\theta =1`$ we get a probability of 45%.

## 4 The second analysis: streamer tubes in drift mode

The performance of the streamer tube system, read out in drift mode, is described in . Here an absolute muon energy calibration was performed at the PS-T9 and SPS-X7 beams, where a slice of the MACRO detector was reproduced. The MCS information was handled by a neural network, based on JETNET 3.0, calibrated with the muon beams quoted above. The neural network(NN) output obtained using test beam data, was compared with that expected from the Monte Carlo simulation, obtaining a satisfactory agreement. The application of such analysis to the MACRO data resulted in an improvement of the space resolution from $`\sigma _x`$$``$$`1\text{cm}`$ to $`\sigma _x`$$``$$`0.3\text{cm}`$.

Figure 2 shows the Monte Carlo prediction for NN downward throughgoing muons and downward going stopping muons compared with experimental data. A nice agreement between simulation and real data is found for both categories, representing a large energy interval.

The result of the neural network output energy calibration shows that NN output is almost linear with $`\mathrm{log}_{10}(E_\mu )`$, increasing with the muon energy up to $`E_\mu `$$``$$`40\text{GeV}`$, where a saturation effect occurs. With respect to the approach used in few details of the neural network were optimized: the neural network training and the energy calibration was performed separately for events with hits in the upper part of the detector(Attico).

Although smeared by the energy carried away by hadrons and by energy loss in the rock, the detected neutrino induced muons still carry on average $``$40$`\%`$ of the original neutrino energy. By using the full Monte Carlo simulation quoted in Sect. 2, we calibrated the NN output as a function of $`\mathrm{log}_{10}(E_\nu `$): the calibrated NN output is linear up to $`\mathrm{log}_{10}(E_\nu (\mathrm{GeV}))`$$``$2.15. We may write:

$$\delta \mathrm{log}_{10}(E_\nu )=\mathrm{log}_{10}(e)\delta \mathrm{ln}E_\nu =\mathrm{log}_{10}(e)\frac{\delta E_\nu }{E_\nu },$$
(6)
where $`\delta \mathrm{log}_{10}(E_\nu )`$ is the difference between the $`\mathrm{log}_{10}`$ of the real and of the reconstructed neutrino energy. Taking into account that the NN output was calibrated as a function of $`\mathrm{log}_{10}(E_\nu )`$, the energy resolution $`\delta E_\nu `$/$`E_\nu `$ was obtained by plotting the quantity $`\delta \mathrm{log}_{10}(E_\nu )`$/$`\mathrm{log}_{10}(e)`$.

In Fig. 3 we show the resolution that could be obtained with an ideal muon energy resolution(dotted line) and that obtained with the present analysis(continuous line).

The precision of the neutrino energy estimate obtained with an ideal muon energy resolution detector is $`\delta E_\nu `$$`/`$$`E_\nu `$$``$70$`\%`$, while with the present method a resolution of $`\delta E_\nu `$/$`E_\nu `$$``$150$`\%`$ is obtained. The asymmetry present in both curves comes from neutrino interactions occurring far from the detector, for which a large fraction of the muon energy is lost during the transport.

### 4.1 Data Selection

We used the 783 upward throughgoing muon data set collected in the full detector run started in 1994. Table 1 describes further event selection to arrive at the sample for MCS analysis. We required a single track in both the wire and the strip view. We selected hits belonging to the track and made of a single fired tube, to associate unambiguously its QTP-TDC time information. This cut is effective for muon tracks with large zenith angles ($`\mathrm{\Theta }`$ $`>`$ 30), while it is quite loose around the vertical: we restricted the present analysis to events with $`\mathrm{\Theta }`$$``$$`60^{}`$. Comparisons were performed between simulated and experimental downward going muon data, to ensure that the selection efficiency, in the used angular window, is the same in the Monte Carlo and in the real data. Spurious background hits have been avoided by requiring a time window of 2 $`\mu `$s around the trigger time. Finally, we selected events with at least four streamer tube planes with valid QTP-TDC data. We fitted the drift circles using the same tracking developed to analyse test beam muons. A minimum path length of $`200\text{cm}`$ in the lower detector is required for tracks hitting the $`Attico`$ and $`400\text{cm}`$ for tracks not hitting it. These geometric requirement ensure a minimum depth of material where muons may experience a measurable amount of multiple scattering and a lever arm long enough for a comfortable tracking. After the selection cuts 300 events survived, giving an overall efficiency of 38.3$`\%`$.

### 4.2 Qualitative tests

We used the information provided by the neural network to separate the neutrino events into four energy subsamples, as shown in Table 2. The same selection was applied to simulated events.

Figure 4 shows the zenith angle distributions of the upward throughgoing muons in the four energy subsamples, compared with the expectations of the Monte Carlo simulation, assuming no-oscillations(dotted line) and oscillations with parameters $`\mathrm{\Delta }\mathrm{m}^2`$= 2.5$`\times `$$`10^3`$$`eV^2`$ and $`\mathrm{sin}^22\theta `$1 (wall boxes). The Monte Carlo including the oscillation hypothesis, with the parameters quoted above, reasonably reproduces the real data in each subsample. We point out the strong difference between the oscillations/no-oscillations hypotheses at low energies, while such difference is reduced by increasing the reconstructed neutrino energy.

Finally, we used information on the ratio $`L_\nu /E_\nu `$. The distance $`L_\nu `$, travelled by the neutrinos from production to the interaction points, was measured by MACRO relying on the muon zenith angle determination, with a precision $`\mathrm{\Delta }L_\nu `$/$`L_\nu 3\%`$. The resolution on the ratio $`L_\nu /E_\nu `$ is therefore fully dominated by the uncertainty in the neutrino energy estimate, giving a relative error of $`150\%`$. The ratio DATA/Monte Carlo as a function of $`\mathrm{log}_{10}(L_\nu /E_\nu )`$, is plotted in Fig. 5. The black circles are the real data over the Monte Carlo predictions, assuming no oscillations, the shaded regions represent the Monte Carlo predictions with oscillations, $`\mathrm{\Delta }\mathrm{m}^2`$ = 2.5$`\times `$$`10^3`$eV2 and $`\mathrm{sin}^22\theta `$=1, divided by the Monte Carlo predictions with no oscillation.

The $`\mathrm{log}_{10}`$($`\frac{L_\nu }{E_\nu }`$) distribution of the neutrinos detected by MACRO spans from 0.8 to 3.5. The left-pointing arrow at low $`\mathrm{log}_{10}`$($`\frac{L_\nu }{E_\nu }`$) represents the effect of the neural network saturation. Since the maximum energy reconstructed by the NN is $`E_\nu `$$``$$`140\text{GeV}`$ and since the minimum $`L_\nu `$ for the present analysis is $`L_\nu ^{min}`$$``$6500 km, the minimum estimated value is $`\mathrm{log}_{10}`$($`\frac{L_\nu }{E_\nu })^{min}`$$``$1.7. As far as the right-pointing arrow at high $`\frac{L_\nu }{E_\nu }`$ is concerned, the reconstruction of neutrino energy based on the residual muon energy, does not allow one to reconstruct the ratio L/E beyond a given limit. This is due to far neutrino interactions, originated by high energy muons, which lost a large fraction of its energy. An ideal detector with ideal muon residual energy resolution, requiring $`E_\mu <`$$`2\text{GeV}`$, would select an average $`\mathrm{log}_{10}`$($`\frac{L_\nu }{E_\nu }`$)=3.2. Finally, the neural network resolution in muon energy estimate results in a maximum value $`\mathrm{log}_{10}`$($`\frac{L_\nu }{E_\nu }`$)=3. The dashed line at 1 is the expectation without oscillations. The last point(empty circle) is obtained from semicontained upward going muon rate. It is not used for the evaluation of the oscillation probabilities (Sect. 4.4) nor in the allowed region plot(Fig. 7). Good agreement is found with the oscillations expected with the parameters quoted above.

### 4.3 Experimental results

To quantify the contribution of the multiple scattering measurement in stand-alone mode, we performed a blind analysis using the Monte Carlo data, looking for the variable, based on the energy estimate, showing the maximum sensitivity to separate the oscillation from the no-oscillation hypothesis. We found that the best performance is given by the ratio:

$$R=N_{low}/N_{high}$$
(7)
where $`N_{low}`$ and $`N_{high}`$ are the number of events with $`E_\nu ^{rec}`$$`<`$$`30\text{GeV}`$ and $`E_\nu ^{rec}`$$`>`$$`130\text{GeV}`$ respectively.

We considered systematic uncertainties due to the neutrino flux calculation and neutrino cross section. Due to the large uncertainty in the absolute flux , we use in this paper only the angular distribution and the ratio between different event categories selected according to the reconstructed energy. Nevertheless, the theoretical uncertainty of the existing neutrino flux, based on the old CR spectrum, has to be accounted for. We varied the input primary cosmic ray spectral index $`\gamma `$ in our simulation by $`\mathrm{\Delta }\gamma `$=$`\pm `$0.05, obtaining a theoretical error on R, $`\left(\frac{\mathrm{\Delta }R}{R}\right)_{flux}`$=$`\pm `$13$`\%`$.

Another source of systematics comes from the neutrino cross section. We looked at R varying the Monte Carlo input cross sections. The most important contribution to the error comes from the “Low” category, since the cross section at low neutrino energy is more uncertain. By comparing the cross section computed under several hypotheses, varying for instance the structure function (,) in the deep inelastic scattering, including or neglecting the contribution of resonant scattering, we found $`\left(\frac{\mathrm{\Delta }R}{R}\right)_\sigma `$=9 $`\%`$. We estimated a total theoretical error $`\left(\frac{\mathrm{\Delta }R}{R}\right)_{theo}`$= $`\sqrt{\left(\frac{\mathrm{\Delta }R}{R}\right)_{flux}^2+\left(\frac{\mathrm{\Delta }R}{R}\right)_\sigma ^2}`$=16$`\%`$. The systematic error on R reconstruction, evaluated in 6$`\%`$, includes the uncertainties on absorber density, drift velocity, streamer tube efficiency and detector acceptance.

Figure 6 shows the ratio R as a function of $`\mathrm{\Delta }\mathrm{m}^2`$, assuming maximal mixing, for the Monte Carlo simulation: the area between the two solid curves represents the 17$`\%`$ systematic error. The Monte Carlo prediction in case of no oscillations is $`R^{}`$=1.5$`\pm `$0.25(theo.+sys.), while for $`\mathrm{\Delta }\mathrm{m}^2`$=2.5$`\times 10^3eV^2`$ and $`\mathrm{sin}^22\theta `$ =1, R=1.00$`\pm `$0.17(theo.+sys.). As pointed out in, the ratio does not have a Gaussian distribution. The errors on the ratio are therefore also not Gaussian: they are reported just to give a crude estimates of the significance. The experimental ratio is $`R^{exp}`$=0.85$`\pm `$0.16(stat.). The one sided probability of measuring a value smaller than the measured one, was computed by using two different methods. In the first one we let the Monte Carlo simulation (no oscillations) to fluctuate according to statistical and systematic errors of the considered ratio. We then evaluated the fraction of events giving a value smaller than the measured one. In a more pessimistic view, rejecting the hypothesis that the lower number of events detected by MACRO with respect to the Monte Carlo expectation is due to oscillation effects, the total number of Monte Carlo events (sum of the numerator and the denominator in the ratios), were normalized to the experimentally measured sum and then they were allowed to fluctuate.

The corresponding one sided probability of measuring a value smaller than $`R^{exp}`$ according to method 1 (2) is 0.75$`\%`$(1.9$`\%`$) corresponding to 2.4$`\sigma `$(2.1 $`\sigma `$).

Finally, we combined the information coming from the energy estimate and from the angular distribution. We considered the ratio:

$$A=N_{vert}/N_{hor},$$
(8)
where $`N_{vert}`$ is the number of upward throughgoing muon events with $`\mathrm{cos}\theta `$$``$-0.7 and $`N_{hor}`$ is the number of events with $`\mathrm{cos}\theta `$$``$-0.4, as discussed in . The probability of measuring a value $`A^{}`$ lower than $`A^{exp}`$ according to the method 1 (2) is P(A’$`<`$$`A^{exp}`$)=0.001$`\%`$(0.01$`\%`$), corresponding to 4.3 $`\sigma `$(3.7$`\sigma `$). It is worth noting that, combining the two independent probabilities on the ratios R and A, the probability that a fluctuation of the expected values(no oscillations) generates the experimental results is P($`A^{}`$$``$$`R^{}`$$`<`$$`A^{exp}`$$``$$`R^{exp}`$)=1.3$`\times `$$`10^6`$(3.2$`\times `$$`10^5`$), corresponding to 4.7 $`\sigma `$(4 $`\sigma `$).

The 90% confidence level allowed regions in the oscillation parameter space have been computed according to the prescription given in . Figure 7 shows the 90$`\%`$ C.L. for the ratio R, the angular distribution  and for their combination.

## 5 Conclusions

Muon multiple coulomb scattering has been used to estimate the energy of the neutrino induced upward throughgoing muons detected by MACRO. The different tests performed on the data, using the estimated energy (separation in

subsamples with different energies, $`L_\nu /E_\nu `$ estimates, etc.) give a consistent picture, all of them supporting the neutrino oscillation hypothesis with parameters $`\mathrm{\Delta }\mathrm{m}^2=`$ 2.5$`\times `$$`10^3eV^2`$ and $`\mathrm{sin}^22\theta `$$``$1. To quantify such effect, the ratio $`R=N_{low}/N_{high}`$ was used, in stand-alone mode and in combination with the angular distribution, $`A=N_{vert}/N_{hor}`$. Both of them are independent of the neutrino absolute flux. The significance of the MACRO observation of the neutrino oscillations is above 4$`\sigma `$.

Acknowledgements

We gratefully acknowledge the support of the Director and of the staff of the Laboratori Nazionali del Gran Sasso and the invaluable assistance of the technical staff of the Institutions participating in the experiment. We thank the Istituto Nazionale di Fisica Nucleare (INFN), the U.S. Department of Energy and the U.S. National Science Foundation for their generous support of the MACRO experiment. We thank INFN, ICTP (Trieste), WorldLab and NATO for providing fellowships and grants (FAI) for non Italian citizens.


text:
# An existence of simple choreographies for N-body problem - a computer assisted proof

## 1. Introduction

In this paper we consider the problem of finding periodic solution to the $`N`$-body problem in which all $`N`$ masses travel along a fixed curve in the plane. The $`N`$-body problem with $`N`$ equal unit masses is given by a differential equation

(1.1)
$$\ddot{q_i}=\underset{ji}{}\frac{q_jq_i}{r_{ij}^3}$$

where $`q_i^n`$, $`i=1,\mathrm{},N`$, $`r_{ij}=q_iq_j`$. The gravitational constant is taken equal to 1.

We consider planar case (n=2) only, we set $`q_i=(x_i,y_i)`$, $`\dot{q_i}=(v_i,u_i)`$. Using this we can express (1.1) by:

(1.2)
$$\{\begin{array}{c}\dot{v_i}=_{ji}\frac{x_jx_i}{r_{ij}^3}\hfill \\ \dot{u_i}=_{ji}\frac{y_jy_i}{r_{ij}^3}\hfill \\ \dot{x_i}=v_i\hfill \\ \dot{y_i}=u_i\hfill \end{array}$$

Recently, this problem received a lot attention in literature see \[M, CM, CGMS, S1, S2, S3, MR\] and papers cited there.

By a simple choreography \[S1, S2\] we mean a collision-free solution of the N-body problem in which all masses move on the same curve with a constant phase shift. This means that there exists $`q:^2`$ a $`T`$-periodic function of time, such that the position of $`k`$-th body ($`k=0,\mathrm{},N1`$) is given by $`q_k(t)=q(t+k\frac{T}{N})`$ and $`(q_0,q_1,\mathrm{},q_{N1})`$ is solution of the $`N`$-body problem. The simplest choreographies are Langrange solutions in which the bodies are located at the vertices of a regular $`N`$-gon and move with constant angular velocity. Another simple choreography, a figure eight curve (see Figure 1), was found numerically by C.Moore \[M\]. A.Chenciner and R.Montgomery \[CM\] gave a rigorous existence proof of the Eight in 2000. In December 1999, J.Gerver found orbit for $`N=4`$ called ’Super-Eight’ (Figure 6). After that C.Simó found a lot of simple choreographies with different shape and number of bodies ranging form 4 to several hundreds (see \[S, S1, CGMS\] for pictures, animations and more details).

Up to now the only choreographies whose existence has been rigorously established \[MR\] are Lagrange solutions and the Eight solution. While the Lagrange solution is given analytically, the existence of Eight was proven in \[CM\] using variational arguments and there is still a lot of open questions about it \[Ch, MR\]. For example uniqueness (up to obvious symmetries and rescaling ) and convexity of the lobes in Eight. In Section 3 we give a computer assisted proof of existence of Eight, its local uniqueness and convexity of the lobes.

In Sections 4, 5 and 6 we concentrate on choreographies called doubly symmetric linear chains. They are symmetric with respect to the $`x`$ and $`y`$ coordinate axes and all points of the curve self-intersection are on the $`x`$ axis. We give a computer assisted proof of an existence (and local uniqueness) of doubly symmetric linear chains for four (Gerver SuperEight) and six bodies.

Proofs given in this paper are computer assisted. By this we mean that we use computer to provide rigorous bounds for solutions of (1.1). The problem of proving an existence of a choreography is reduced to finding a zero for a suitable function. For this purpose we use an interval Newton method and a Krawczyk method \[A, K, Mo, N\], (see Section 2) which are apparently rather unknown outside the interval arithmetic community. To integrate equations (1.1) we use a $`C^1`$-Lohner algorithm \[ZLo\]. All computations were performed on AMD Athlon 1700XP  with 256 MB DDRAM memory, with Windows 98SE operating system. We used CAPD package \[Capd\] and Borland C++ 5.02 compiler. A total computation time for 6-bodies was under 90 seconds and considerably smaller for the Eight and SuperEight solutions - see Section 7 for more details.

## 2. Two zero finding methods

The main technical tool used in this paper in order to establish an existence of solutions of equations of the form $`f(x)=0`$ is an interval Newton method\[A, Mo, N\] and a Krawczyk method\[A, K, N\]. The interval Newton method was used to prove an existence of the Eight (Figure 1) and Gerver orbit (Figure 6). The Krawczyk method for Gerver orbit and an orbit with 6 bodies in a linear chain (Figure 8).

### 2.1. Notation

In the application of interval arithmetics to rigorous verification of theorems single valued objects, like numbers, vectors, matrices etc are in the formulas replaced by sets containing sure bounds for them. In the sequel, we will not use any special notation for single valued object and sets. For a set $`S`$ by $`[S]`$ we denote the interval hull of $`S`$, i.e. the smallest product of intervals containing $`S`$. For a set which is an interval set (i.e. can be represented as a product of intervals) we will also use square brackets to stress its interval nature. For any interval set $`[S]`$ by $`mid([S])`$ we denote a center point of $`[S]`$. For any interval $`[a,b]`$ we define a diameter by $`diam[a,b]:=ba`$. For an interval vector (matrix) $`S=[S]`$ by $`diamS`$ we denote a vector (matrix) of diameters of each components.

### 2.2. Interval Newton method

###### Theorem 2.1.

\[A, N\] Let $`F:^n^n`$ be a $`C^1`$ function. Let $`[X]=\mathrm{\Pi }_{i=1}^n[a_i,b_i]`$, $`a_i<b_i`$. Assume the interval hull of $`DF([X])`$, denoted here by $`[DF([X])]`$, is invertible. Let $`\overline{x}X`$ and we define

(2.1)
$$N(\overline{x},[X])=[DF([X])]^1F(\overline{x})+\overline{x}$$

Then

1. if $`x_1,x_2[X]`$ and $`F(x_1)=F(x_2)`$, then $`x_1=x_2`$
2. if $`N(\overline{x},[X])[X]`$, then $`!x^{}[X]`$ such that $`F(x^{})=0`$
3. if $`x_1[X]`$ and $`F(x_1)=0`$, then $`x_1N(\overline{x},[X])`$
4. if $`N(\overline{x},[X])[X]=\mathrm{}`$, then $`F(x)0`$ for all $`x[X]`$

### 2.3. Krawczyk method

We assume that:

* $`F:^n^n`$ is a $`C^1`$ function,
* $`[X]^n`$ is an interval set,
* $`\overline{x}[X]`$
* $`C^{n\times n}`$ is a linear isomorphism.

The Krawczyk \[A, K, N\] operator is given by

(2.2)
$$K(\overline{x},[X],F):=\overline{x}CF(\overline{x})+(IdC[DF([X])])([X]\overline{x}).$$

###### Theorem 2.2.

1. If $`x^{}[X]`$ and $`F(x^{})=0`$, then $`x^{}K(\overline{x},[X],F)`$.

2. If $`K(\overline{x},[X],F)int[X]`$, then there exists in \[X\] exactly one solution of equation $`F(x)=0`$.

3. If $`K(\overline{x},[X],F)[X]=\mathrm{}`$, then $`F(x)0`$ for all $`x[X]`$

### 2.4. Algorithm for Newton and Krawczyk method

Theorems 2.1 and 2.2 can be used as a basis for an algorithm for rigorous enclosing for solution of equation $`F(x)=0`$. Let $`T(\overline{x},[X])=N(\overline{x},[X])`$ if we are using the interval Newton method and $`T(\overline{x},[X])=K(\overline{x},[X],F)`$ in case of Krawczyk method.

First, we need to have a good guess for $`x^{}^n`$. For this purpose we use a nonrigorous Newton method to obtain $`\overline{x}=(\overline{x}_1,\overline{x}_2,\mathrm{},\overline{x}_n)`$. Then we choose interval set $`[X]`$ which contains $`\overline{x}`$ and perform the following algorithm:

* Compute $`T(\overline{x},[X])`$.
* If $`T(\overline{x},[X])[X]`$, then return success.
* If $`XT(\overline{x},[X])=\mathrm{}`$, then return fail. There are no zeroes of $`F`$ in $`[X]`$.
* If $`[X]T(\overline{x},[X])`$, then modify computation parameters (for example: a time step, the order of Taylor method, size of $`[X]`$). Go to Step 1.
* Define a new $`[X]`$ by $`[X]:=[X]T(\overline{x},[X])`$ and a new $`\overline{x}`$ by $`\overline{x}:=mid([X])`$, then go to Step 1.

In practical computation it is convenient to define a maximum number of iteration allowed and return fail if the actual iteration count is larger.

Observe that the third assertion in both Theorems 2.1 and 2.2 can be used to exclude an existence of zero of $`F`$. This has been used by Galias \[G1, G2\] to find all periodic orbits up to a given period for Hénon map and Ikeda map. It seems possible to obtain similar results for N-body problem in the future.

## 3. The Eight - existence, local uniqueness and convexity

Existence of the Eight has been shown in \[CM\] using mixture of symmetry and variational arguments. Here we give an another existence proof and in addition we obtain local uniqueness and the convexity of each lobe of the Eight. We follow \[CM\] in the use of the symmetry, but other component of the proof is different - we use an interval Newton method discussed in Section 2.

In notation, but only for the Eight, we follow \[CM\].

Let $`T`$ be any positive real number. We define action of the Klein group $`_2\times _2`$ on $`/T`$ and on $`^2`$ as follows: if $`\sigma `$ and $`\tau `$ are generators,

(3.1)
$$\sigma (t)=t+\frac{T}{2},\tau (t)=t+\frac{T}{2},\sigma (x,y)=(x,y),\tau (x,y)=(x,y).$$

For loop $`q:(/T)^2`$ and $`i=1,2,3`$ we define the position for $`i`$-th body by

(3.2)
$$q_i(t)=q\left(t+(3i)\frac{T}{3}\right),$$

Here $`q_i`$ is just position of $`i`$-th body.

The following theorem without the uniqueness part was proved in \[CM\].

###### Theorem 3.1.

There exists an ”eight”-shaped planar loop $`q:(/T)^2`$ with the following properties:

1. for each t,

$$q_1(t)+q_2(t)+q_3(t)=0;$$
2. q is invariant with respect to the action of $`_2\times _2`$ on $`/T`$ and on $`^2`$:

$$q\sigma (t)=\sigma q(t)\text{ and }q\tau (t)=\tau q(t);$$
3. the loop $`x:/T^6`$ defined by

$$x(t)=(q_1(t),q_2(t),q_3(t))$$

is $`T`$-periodic solution of the planar three-body problem with equal masses.

Moreover, $`q`$ is locally unique (up to obvious symmetries and rescaling).

###### Remark 3.2.

If conditions (1),(2),(3) are satisfied then

1. $`\dot{q}_1(t)+\dot{q}_2(t)+\dot{q}_3(t)=0`$ for each t,
2. $`\dot{q}\sigma (t)=\sigma \dot{q}(t)\text{ and }\dot{q}\tau (t)=\sigma \dot{q}(t)`$ for each t,
3. $`q_3(0)=(0,0)`$ and $`\dot{q_3}(0)=2\dot{q_1}(0)`$
4. $`q_1(0)=q_2(0)`$ and $`\dot{q}_1(0)=\dot{q}_2(0)`$,
5. $`q_1(T/12)`$ is on the X axis and $`\dot{q}_1(T/12)`$ is orthogonal to the X axis,
6. $`q_3(T/12)=\tau q_2(T/12)`$ and $`\dot{q}_3(T/12)=\sigma \dot{q}_2(T/12)`$,

The following lemma describes the symmetry reduction for the Eight

###### Lemma 3.3.

Assume that $`\stackrel{~}{q}:[0,\stackrel{~}{T}]^6`$ is a solution of the three body problem, such that $`\stackrel{~}{q}=(\stackrel{~}{q}_1,\stackrel{~}{q}_2,\stackrel{~}{q}_3)`$ satisfies conditions (c),(d),(e),(f) in Remark 3.2, then exists $`q:(/T)^2`$ satisfying condition (1),(2),(3) in Theorem 3.1 with $`T=12\stackrel{~}{T}`$.

Proof: We define

(3.3)
$$\widehat{q}(t)=\{\begin{array}{cc}\stackrel{~}{q}_3(t)\hfill & \text{for }t[0,\stackrel{~}{T}]\hfill \\ \tau \stackrel{~}{q}_2(2\stackrel{~}{T}t)\hfill & \text{for }t[\stackrel{~}{T},2\stackrel{~}{T}]\hfill \\ \sigma \stackrel{~}{q}_1(t2\stackrel{~}{T})\hfill & \text{for }t[2\stackrel{~}{T},3\stackrel{~}{T}]\hfill \end{array}$$

and

(3.4)
$$q(t)=\{\begin{array}{cc}\widehat{q}(t)\hfill & \text{for }t[0,3\stackrel{~}{T}]\hfill \\ \tau \widehat{q}(\tau ^1(t))\hfill & \text{for }t[3\stackrel{~}{T},6\stackrel{~}{T}]\hfill \\ \sigma \widehat{q}(\sigma ^1(t))\hfill & \text{for }t[6\stackrel{~}{T},9\stackrel{~}{T}]\hfill \\ \sigma \tau \widehat{q}(\tau ^1\sigma ^1(t))\hfill & \text{for }t[9\stackrel{~}{T},12\stackrel{~}{T}]\hfill \end{array}$$

Let $`f_1`$ and $`f_2`$ be two solutions of 1.1 on intervals $`[t_1,t_2]`$ and $`[t_2,t_3]`$ respectively. If $`f_1(t_2)=f_2(t_2)`$ and $`\dot{f}_1^{}(t_2)=\dot{f}_2^+(t_2)`$ then $`f=\{f_1,f_2\}`$ is solution on interval $`[t_1,t_2]`$. To show that $`q(t)`$ is a solution it is enough to show that ”pieces fit together smoothly”.

For $`t=\stackrel{~}{T}`$ from 3.2.f we have

(3.5) $`\stackrel{~}{q}_3(\stackrel{~}{T})`$ $`=`$ $`\stackrel{~}{q}_3({\displaystyle \frac{T}{12}})=\tau \stackrel{~}{q}_2({\displaystyle \frac{T}{12}})=\tau \stackrel{~}{q}_2(\stackrel{~}{T})=\tau \stackrel{~}{q}_2(2\stackrel{~}{T}\stackrel{~}{T})`$
(3.6) $`\dot{\stackrel{~}{q}}_3^{}(\stackrel{~}{T})`$ $`=`$ $`\sigma \dot{\stackrel{~}{q}}_2^+(\stackrel{~}{T})`$
For $`t=2\stackrel{~}{T}`$ from 3.2.d we have

(3.7) $`\stackrel{~}{q}(2\stackrel{~}{T})`$ $`=`$ $`\tau \stackrel{~}{q}_2(0)=\sigma \stackrel{~}{q}_1(0)`$
(3.8) $`\dot{\stackrel{~}{q}}^{}(2\stackrel{~}{T})`$ $`=`$ $`\sigma \dot{\stackrel{~}{q}}_2^{}(\stackrel{~}{T})=\sigma \dot{\stackrel{~}{q}}_1^+(2\stackrel{~}{T})=\dot{\stackrel{~}{q}}^+(2\stackrel{~}{T})`$
From (3.5-3.8) it follows that $`\widehat{q}(t)`$ is smooth curve for $`t[0,3\stackrel{~}{T}]`$. We construct $`q(t)`$ from that curve.

For $`t=3\stackrel{~}{T}`$ from 3.2.e we obtain

$`q(3\stackrel{~}{T})`$ $`=`$ $`\sigma \stackrel{~}{q}_1(\stackrel{~}{T})=\sigma \tau \stackrel{~}{q}_1(\stackrel{~}{T})=\sigma \tau \widehat{q}(\tau ^1(3\stackrel{~}{T}))`$
$`\dot{q}^{}(3\stackrel{~}{T})`$ $`=`$ $`\sigma \dot{\stackrel{~}{q}}_1(\stackrel{~}{T})=\dot{\stackrel{~}{q}}_1(\stackrel{~}{T})=\dot{q}^+(3\stackrel{~}{T})`$
Using 3.2.c for $`t=6\stackrel{~}{T}`$ we infer that

$`q(6\stackrel{~}{T})`$ $`=`$ $`\tau \widehat{q}(\tau ^1(6\stackrel{~}{T}))=\tau \stackrel{~}{q}_3(0)=\sigma \stackrel{~}{q}_3(0)=\sigma \widehat{q}(\sigma ^1(6\stackrel{~}{T}))`$
$`\dot{q}^{}(6\stackrel{~}{T})`$ $`=`$ $`\sigma \dot{\stackrel{~}{q}}_1(\stackrel{~}{T})=\dot{q}^+(6\stackrel{~}{T})`$
Now we show that $`q(t)`$ is a closed curve.

$`q(12\stackrel{~}{T})=\sigma \tau \widehat{q}(\tau ^1\sigma ^1(12\stackrel{~}{T}))=\sigma \tau \widehat{q}(0)=\widehat{q}(0)=q(0)`$
$`\dot{q}^{}(12\stackrel{~}{T})={\displaystyle \frac{}{t}}[\sigma \tau \widehat{q}(\tau ^1\sigma ^1(12\stackrel{~}{T}))]=\dot{\widehat{q}}(0)=\dot{q}^+(0)`$
Hence it is easy to see that $`q(t)`$ can be extended to a $`T`$-periodic curve, such that $`x(t)=(q_1(t),q_2(t),q_3(t))`$ is a $`T`$-periodic solution of the three body problem.

Condition (2) in Theorem 3.1 follows easily from definition of $`\widehat{q}(t)`$ and properties of $`\sigma `$ and $`\tau `$ (for T-periodic orbit $`\sigma ^1(t)=\sigma (t)`$ and $`\tau ^1(t)=\tau (t)`$). For $`t[0,3\stackrel{~}{T}]`$ we have

(3.9) $`q(\sigma (t))`$ $`=`$ $`\sigma \widehat{q}(\sigma ^1\sigma (t))=\sigma \widehat{q}(t)=\sigma q(t)`$
(3.10) $`q(\tau (t))`$ $`=`$ $`\tau \widehat{q}(\tau ^1\tau (t))=\tau \widehat{q}(t)=\tau q(t)`$
For $`t[3\stackrel{~}{T},6\stackrel{~}{T}]`$ using (3.10) we obtain

(3.11) $`q(\sigma (t))=\sigma \tau \widehat{q}(\tau ^1\sigma ^1\sigma (t))=\sigma \tau \widehat{q}(\tau ^1(t))=\sigma q(t)`$
(3.12) $`q(\tau (t))=\widehat{q}(\tau (t))=\tau \widehat{q}(\tau \tau (t))=\tau q(t)`$
We omit other two cases, because the proof is very similar.

Because there aren’t any external force so the center of mass can only move with a constant velocity. But from 3.2.c and 3.2.d we have $`q_1(0)+q_2(0)+q_3(0)=0`$ and $`\dot{q}_1(0)+\dot{q}_2(0)+\dot{q}_3(0)=0`$. Hence we obtain condition (1) in Theorem 3.1. ∎

Hence to prove Theorem 3.1 it is enough to show that there exists a locally unique (up to obvious degeneracies) function satisfying assumptions of Lemma 3.3. For this end we rewrite these assumptions as a zero finding problem to which we apply an interval Newton method in the reduced space.

Our original phase space is 12 dimensional, the state of bodies is given by $`(q_1,q_2,q_3,\dot{q_1},\dot{q_2},\dot{q_3})`$. The center of mass is fixed at the origin. Hence one body’s position and velocity is determined by other two bodies. We start from a collinear position with the third body at the origin and with equal velocities of the first and the second body (see 3). Hence it is enough to know the position and the velocity of the first body to reconstruct initial condition of other bodies. Moreover, if we have one solution we could get another solution by a suitable rotation (both have the same shape). To remove this degeneracy we place the first body on the X axis (we will make computation for rotated Eight, see Figure 3). In addition we fix the size of trajectory. This fixes the period of the solution. But if we know one periodic solution then by Kepler law we may obtain a solution of any period just by rescaling. Finally, we set $`q_1(0)=(1,0)`$, hence initial conditions are fixed by the velocity of the first body.

The reduced space for Eight is two dimensional and is parameterized by velocity of second body. We define map from the reduced space to full phase space $`E:^2^{12}`$, which expands velocity of first body, given by $`(v,u)`$, to the initial conditions $`(x_1,y_1,x_2,y_2,x_3,y_3,v_1,u_1,v_2,u_2,v_3,u_3)`$ for equation (1.1)

$$E(v,u)=(1,0,1,0,0,0,v,u,v,u,2v,2u).$$

For each such initial condition exists a solution of the tree body problem defined on some interval. To each initial configuration, following that solution, we associate, if it exists, a configuration in which for the first body the position vector is orthogonal to its velocity vector for the first time. This defines the Poincaré map $`P:^{12}\mathrm{\Omega }^{12}`$.

Now, we define map $`R:^{12}^2`$, by

$$R(q_1,q_2,q_3,\dot{q}_1,\dot{q}_2,\dot{q}_3)=(q_2q_1^2q_3q_1^2,(\dot{q}_2\dot{q}_3)\times q_1),$$

where by $`\times `$ we denote vector product, and map $`\mathrm{\Phi }:^2^2`$ by

$$\mathrm{\Phi }=RPE.$$

###### Remark 3.4.

If $`R(q_1,q_2,q_3,\dot{q}_1,\dot{q}_2,\dot{q}_3)=0`$, then $`q_2q_1=q_3q_1`$ and if in addition $`y_1=0`$ but $`x_10`$, then $`\dot{y}_2=\dot{y}_3`$. The bodies are then located as in Fig. 4 and conditions (e) and (f) in Rem. 3.2 are satisfied in a suitably rotated coordinate frame.

The following lemma, which is crucial for the proof of Theorem 3.1 is obtained with computer assistance.

###### Lemma 3.5.

There exists a locally unique $`(v,u)^2`$ that $`\mathrm{\Phi }(v,u)=(0,0)`$.

Proof: We use an interval Newton method (Theorem 2.1). First we come close to zero, starting with some rough initial condition (for example from \[S1\]) using a non-rigorous Newton method. Once we have a good candidate $`x_0=(v_0,u_0)`$, we set $`[X]=[v_0\delta ,v_0+\delta ]\times [u_0\delta ,u_0+\delta ]`$ and compute rigorously $`\mathrm{\Phi }(x_0)`$ and $`\frac{\mathrm{\Phi }([X])}{x}`$. For this purpose we use $`C^1`$-Lohner algorithm described in \[ZLo\]. In this computation we used a time step $`h=0.01`$ and the order $`r=7`$.

It turns out that the assumption of assertion 1 in Theorem 2.1 holds for
$`x_0=(0.347116768716,0.532724944657)`$ and $`\delta =10^6`$. Numerical data from this computation are listed in Table 1.

Moreover, from Theorem 2.1 we know that this zero is unique in the set X.∎

Proof of Theorem 3.1: From Lemma 3.5 there exists an initial condition $`(v,u)`$ in reduced space that $`\mathrm{\Phi }(v,u)=(0,0)`$. Hence there exists solution $`\overline{q}(t)`$ of the three body problem in some interval $`[0,\stackrel{~}{T}]`$ that in $`t=0`$ bodies are in collinear position with third body in the origin and in $`t=\stackrel{~}{T}`$ they are in isosceles configuration ($`r_{12}=r_{13}`$). By rotation we could get solution $`\stackrel{~}{q}(t)`$ that in $`t=\stackrel{~}{T}`$ first body is on X axis. Obviously we still start in collinear position and end in isosceles. Now from Remark 3.4 $`\stackrel{~}{q}(t)`$ satisfy all condition in Lemma 3.3 and hence there exist $`q(t)`$ satisfying condition (1),(2),(3).

Local uniqueness follows from local uniqueness in Lemma 3.5.∎

### 3.1. Convexity of the Eight

###### Theorem 3.6.

Each lobe of the Eight is convex.

Proof: For the proof it is enough to show that the only inflection point on the curve $`q(t)`$ is the origin.

From Lemma 3.5 we get set $`[X]`$ which includes initial condition for the Eight in reduced space. We expand $`[X]`$ to the full space, we set $`\overline{X}=E([X])`$. In the coordinate frame in which the Eight looks like in Fig. 1 and symmetries $`\sigma `$ and $`\tau `$ are reflections with respect to coordinate axes (see 3.1) we see immediately that the symmetry properties of the Eight imply that at the origin $`\frac{^2y_i}{x_i^2}=0`$ and $`\frac{^2x_i}{y_i^2}=0`$. To prove convexity of the Eight we follow rigorously the trajectory of set $`\overline{X}`$ and show that the only point in which $`\frac{^2y_i}{x_i^2}=0`$ and $`\frac{^2x_i}{y_i^2}=0`$ is the origin. The same is true if we rotate the Eight, as in Figure 3, to the coordinate system in which we performed the actual computations and in which we will work for the remainder of this proof.

Let $`q_i(t)=(x_i(t),y_i(t))`$ be position of i-th body . If $`0\frac{x_i}{t}[t_{k1},t_k]`$ then, from the implicit function theorem, we can write $`y_i`$ as a function of $`x_i`$ for interval between $`x_i(t_{k1})`$ and $`x_i(t_k)`$. Otherwise we try to represent $`x_i`$ as a function of $`y_i`$. For small enough time steps at least one of this representations is always possible for the Eight (this is really verified during the rigorous computations) . Now we derive the formulas for the derivatives of $`y_i(t(x_i))`$ with respect to $`x_i`$ (for derivatives of $`x_i(t(y_i))`$ we just exchange $`x`$ and $`y`$ variables in formulas).

$`{\displaystyle \frac{y_i(t(x_i))}{x_i}}`$ $`=`$ $`{\displaystyle \frac{y_i}{t}}\left({\displaystyle \frac{x_i}{t}}\right)^1`$
$`{\displaystyle \frac{^2y_i(t(x_i))}{x_i^2}}`$ $`=`$ $`\left({\displaystyle \frac{^2y_i}{t^2}}{\displaystyle \frac{^2x_i}{t^2}}{\displaystyle \frac{y_i}{x_i}}\right)\left({\displaystyle \frac{x_i}{t}}\right)^2`$
$`{\displaystyle \frac{^3y_i(t(x_i))}{x_i^3}}`$ $`=`$ $`{\displaystyle \frac{\frac{^3y_i}{t^3}\frac{x_i}{t}\frac{^3x_i}{t^3}\frac{y_i}{t}+2\left(\frac{^2x_i}{t^2}\right)^2\frac{y_i}{x_i}2\frac{^2x_i}{t^2}\frac{^2y_i}{t^2}}{\left(\frac{x_i}{t}\right)^4}}{\displaystyle \frac{\left(\frac{^2x_i}{t^2}\frac{^2y_i}{x_i^2}\right)}{\left(\frac{x_i}{t}\right)^2}}`$
From this equations we need to know the derivatives $`x_i`$ and $`y_i`$ with respect to time. We obtain them easily by differentiation of (1.1) with respect to time. In fact this is done during each step of $`C^1`$-Lohner algorithm.

By $`\phi (t,x_0)`$ we denote the state of bodies (the position and the velocity) at time $`t`$ with the initial condition $`x_0`$ at $`t=0`$. Let $`h_k`$ be length of $`k`$-th time step, $`t_k=h_1+\mathrm{}+h_k`$ \- the time after $`k`$ steps, $`[q^k]=[q_1^k]\times [q_2^k]\times [q_3^k]\times [\dot{q}_1^k]\times [\dot{q}_2^k]\times [\dot{q}_3^k]^{12}`$ be product of intervals, such that $`\phi (t_k,\overline{X})[q^k]`$ and $`[Q^k]=[Q_1^k]\times [Q_2^k]\times [Q_3^k]\times [\dot{Q}_1^k]\times [\dot{Q}_2^k]\times [\dot{Q}_3^k]^{12}`$, where $`[Q_i^k]=[X_i^k]\times [Y_i^k]`$, be product of intervals, such that $`\phi ([t_{k1},t_k],\overline{X})[Q^k]`$ Note that $`[Q^k]`$ can be seen as an interval enclosure for whole trajectory between $`[q^{k1}]`$ and $`[q^k]`$, which is computed during each step of $`C^1`$-Lohner algorithm.

For $`k`$-th step and $`i`$-th body (except first step for third body) we check if at least one of the following conditions is true

(3.13) $`0{\displaystyle \frac{x_i}{t}}[X_i^k]\text{ and }0{\displaystyle \frac{^2x_i}{y_i^2}}[X_i^k]`$
(3.14) $`0{\displaystyle \frac{y_i}{t}}[Y_i^k]\text{ and }0{\displaystyle \frac{^2y_i}{x_i^2}}[Y_i^k].`$
For first step for third body (starting in the origin) we check if one of following conditions is satisfied

(3.15) $`0{\displaystyle \frac{x_3}{t}}[X_3^1]\text{ and }0{\displaystyle \frac{^2x_3}{y_3^2}}[X_3^1]\text{ and }0{\displaystyle \frac{^3x_3}{y_3^3}}[X_3^1]`$
(3.16) $`0{\displaystyle \frac{y_3}{t}}[Y_3^1]\text{ and }0{\displaystyle \frac{^2y_3}{x_3^2}}[Y_3^1]\text{ and }0{\displaystyle \frac{^3y_3}{x_3^3}}[Y_3^1].`$
To verify above conditions we follow rigorously using $`C^1`$-Lohner algorithm the trajectory of set $`\overline{X}`$ until reaching the section described in proof of Lemma 3.5 and for each time step and each body we check suitable condition. It turns out that for each time step and each body these conditions were satisfied. This finishes the proof. ∎

### 3.2. Some numerical data from the convexity of Eight computation

The parameters of the methods were: a time step $`h=0.01`$, order $`r=7`$. We needed 53 steps to cross the section, below we show data for characteristic cases in order to show that it was really quite easy to verify, with algorithms we used.

Step 1. We start in collinear position with third body in the origin. So we are in an inflection point.

Numerical data for this case are given in Table 2. We see that for third body we have $`0\frac{^2x_i}{y_i^2}[X_i^1]`$, but this derivative is monotonic ($`\frac{^3x_i}{y_i^3}[X_i^1]<0`$), hence there can be only one zero of it in interval $`[X_3^1]`$. But we know that one zero is at $`(0,0)[X_3^1]`$.

Steps 2-36 and 38-53 In this case all second derivatives do not contain 0. Table 3 contains data obtained in second step.

Step 37. In this case we cannot represent $`y_1`$ as function of $`x_1`$, hence we interchange variables and represent $`x_1`$ as function of $`y_1`$ and check condition 3.14 instead of 3.13. The first body is in the rightmost position. Table 4 contains the derivatives for this case.

## 4. Doubly symmetric choreographies with even number of bodies

### 4.1. Symmetries

Many of the choreographies found by Simo \[S1\] have at least one symmetry. Just as in case of the Eight this is not only a symmetry of trajectory image, but also of how the curve describing the trajectory is parameterized with time.

In this section we introduce a notation for symmetries which will be used till the end of this paper. By $`S_x`$, $`S_y`$, $`S_0`$ we want to denote a symmetry with respect to the x axis, the y axis and to the origin. These spatial symmetries act also on time variable parameterizing curves.
Let $`T`$ be any positive real number. We define actions of $`S_x`$, $`S_y`$ and $`S_0`$ on $`/T`$ and on $`^2`$ as follows:

$$\begin{array}{cc}S_x(t)=t+\frac{T}{2},\hfill & S_x(x,y)=(x,y),\hfill \\ S_y(t)=t,\hfill & S_y(x,y)=(x,y),\hfill \\ S_0(t)=t+\frac{T}{2},\hfill & S_0(x,y)=(x,y).\hfill \end{array}$$

It follows from this definition that $`S_0=S_xS_y=S_yS_x`$.

Let $`q(t):/T^2`$. We say that $`q(t)`$ is invariant (equivariant) with respect to the action of S if $`S(q(t))=q(S(t))`$ for all $`t`$. Moreover we assume that $`q(t)`$ is $`C^2`$ function. If $`q(t)`$ is invariant with respect to $`S_x`$ (resp. $`S_y`$) then $`\dot{q}(S_x(t))=S_y(\dot{q}(t))`$ (resp. $`\dot{q}(S_y(t))=S_x(\dot{q}(t))`$). Hence $`S_0`$ invariance implies that $`\dot{q}(S_0(t))=S_0(\dot{q}(t))`$.

From now on we will enumerate bodies starting from 0. We set $`q_i(t)=q\left(t+\frac{T}{N}i\right)`$ for $`i=0,\mathrm{},N1`$.

### 4.2. Doubly symmetric choreographies with even number of bodies

We will consider only cases with even number of bodies.

Let $`T=N\overline{T}>0`$. We search for a function $`q(t):/T^2`$ which has following properties:

1. for each $`t`$ the origin is center of mass,

(4.1)
$$\underset{i=0}{\overset{N1}{}}q\left(t+\overline{T}i\right)=0,$$
2. $`q(t)`$ is invariant with respect to

1. $`S_x`$ i.e. $`q(S_x(t))=S_x(q(t))`$,
    2. $`S_0`$ i.e. $`q(S_0(t))=S_0(q(t))`$,
3. function $`x=(q_0(t),q_2(t),\mathrm{},q_{N1})`$ where $`q_i(t)=q\left(t+\overline{T}i\right)`$ for $`i=0,\mathrm{},N1`$, is a T-periodic solution of the Newtonian N-body problem (1.1).

Observe that condition P3 says that bodies trace each other with a constant phase shift.

The following Lemma, which is analogous to Lemma 3.3, gives necessary and sufficient conditions for an existence of a choreography satisfying P1,P2 and P3.

###### Lemma 4.1.

Let $`N=2n`$ be a number of bodies. There exists a function q(t) with properties P1, P2, P3 if and only if there are functions $`q_i:[0,\overline{T}/2]R^2`$ for $`i=0,1,\mathrm{},N1`$ such as:

1. $`q_0(0)=(0,y_0)`$ for some $`y_00`$ ($`q_0(0)`$ is on Y axis),
2. At time $`t_0=0`$ for $`i=0,1,\mathrm{},N/21`$ we have

1. $`q_i(t_0)=S_x(q_{\frac{N}{2}i}(t_0))`$
    2. $`\dot{q}_i(t_0)=S_y(\dot{q}_{\frac{N}{2}i}(t_0))`$
    3. $`q_i(t_0)=S_0(q_{\frac{N}{2}+i}(t_0))`$
    4. $`\dot{q}_i(t_0)=S_0(\dot{q}_{\frac{N}{2}+i}(t_0))`$
3. At time $`t_1=\overline{T}/2`$ for $`i=0,1,\mathrm{},N/21`$ we have

1. $`q_i(t_1)=S_x(q_{\frac{N}{2}i1}(t_1))`$
    2. $`\dot{q}_i(t_1)=S_y(\dot{q}_{\frac{N}{2}i1}(t_1))`$
    3. $`q_i(t_1)=S_0(q_{\frac{N}{2}+i}(t_1))`$
    4. $`\dot{q}_i(t_1)=S_0(\dot{q}_{\frac{N}{2}+i}(t_1))`$
4. $`x(t)=(q_0(t),q_1(t),\mathrm{},q_{N1}(t))`$ is a solution of the Newton N-body problem for $`t[0,\overline{T}/2]`$.

Proof: If we have a trajectory with properties P1, P2, P3 then it is easy to show that functions $`q_i(t)`$ defined in P3 satisfy conditions (1)-(3). On the other hand if we have a functions $`q_i(t)`$ which have properties (1)-(3) then we define $`q(t)`$ by

(4.2)
$$q(t)=\{\begin{array}{cc}q_i(ti\overline{T})\hfill & \text{for }t[i\overline{T},(i+\frac{1}{2})\overline{T}]\text{ and }i=0,1,\mathrm{},N1\hfill \\ S_x(q_{\frac{N}{2}i}(i\overline{T}t))\hfill & \text{for }t[(i\frac{1}{2})\overline{T},i\overline{T}]\text{ and }i=1,2,\mathrm{},\frac{N}{2}\hfill \\ S_x(q_{\frac{3N}{2}i}(i\overline{T}t))\hfill & \text{for }t[(i\frac{1}{2})\overline{T},i\overline{T}]\text{ and }i=\frac{N}{2}+1,\mathrm{},N\hfill \end{array}$$

From Lemma 4.1 it follows that the proof of an existence of a doubly symmetric choreography is equivalent to some boundary value problem for the N-body equation. We will now formulate this problem as a zero finding problem for a suitable map.

Original phase space for the planar N-body problem has $`4N`$ dimensions. It turns out that, if initial conditions satisfy all conditions in point 2 of Lemma 4.1 it is enough to know values of only $`N`$ variables to recover rest of them. We still may obtain solutions of any period, to determine this we can fix the size of curve by fixing one variable. Hence our reduced space be $`(N1)`$-dimensional. In the next paragraph we will be more specific.

We define map $`E:^{N1}^{4N}`$, which expands, using symmetries from (2), initial conditions from reduced space to the full phase space:
$`(x_0,y_0,\dot{x}_0,\dot{y}_0,x_1,y_1,\dot{x}_1,\dot{y}_1,\mathrm{},x_{N1},y_{N1},\dot{x}_{N1},\dot{y}_{N1})`$. We consider two cases: $`N=4k`$ and $`N=4k+2`$.

For $`N=4k`$ we set

$`E(\dot{x}_0\times {\displaystyle \underset{i=1}{\overset{k1}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i)\times (x_k,\dot{y}_k))=(0,a,\dot{x}_0,0)\times {\displaystyle \underset{i=1}{\overset{k1}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i)`$
$`\times (x_k,0,0,\dot{y}_k)\times {\displaystyle \underset{i=1}{\overset{k1}{}}}(x_{ki},y_{ki},\dot{x}_{ki},\dot{y}_{ki})`$
$`\times (0,a,\dot{x}_0,0)\times {\displaystyle \underset{i=1}{\overset{k1}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i)`$
$`\times (x_k,0,0,\dot{y}_k)\times {\displaystyle \underset{i=1}{\overset{k1}{}}}(x_{ki},y_{ki},\dot{x}_{ki},\dot{y}_{ki})`$
For $`N=4k+2`$ we set

$`E(\dot{x}_0\times {\displaystyle \underset{i=1}{\overset{k}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i))=(0,a,\dot{x}_0,0)\times {\displaystyle \underset{i=1}{\overset{k}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i)`$
$`\times {\displaystyle \underset{i=0}{\overset{k1}{}}}(x_{ki},y_{ki},\dot{x}_{ki},\dot{y}_{ki})\times (0,a,\dot{x}_0,0)`$
$`\times {\displaystyle \underset{i=1}{\overset{k}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i)\times {\displaystyle \underset{i=0}{\overset{k1}{}}}(x_{ki},y_{ki},\dot{x}_{ki},\dot{y}_{ki})`$
In both cases $`a`$ is a parameter fixing a size of an orbit.

Again, like for the Eight, for every such initial condition exists a solution of the $`N`$-body problem (1.1) defined on some interval. To every initial configuration, following that solution, we associate, if it exists, a configuration in which for the first time

* for $`N=4k`$: bodies $`k`$ and $`k1`$ have equal $`x`$ coordinate ($`x_k=x_{k1}`$),
* for $`N=4k+2`$: $`k`$-th body is on the X axis ($`y_k=0`$).

This defines the Poincaré map $`P:^{4N}\mathrm{\Omega }^{4N}`$.

Now, we define the reduction map $`R:^{4N}^{N1}`$ in such way that $`R`$ has zeroes in points satisfying conditions (3a) and (3b) in Lemma 4.1 and only in such points. Observe that we don’t have to worry about conditions (3c) and (3d) in Lemma 4.1, because from the properties of (1.1) if follows that if (2c) and (2d) holds, then (3c) and (3d) are satisfied for any $`t_1`$.

For $`N=4k`$ we set

$`R\left({\displaystyle \underset{i=0}{\overset{N1}{}}}(x_i,y_i,\dot{x}_i,\dot{y}_i)\right)=(y_k+y_{k1},\dot{x}_k+\dot{x}_{k1},\dot{y}_k\dot{y}_{k1})`$
$`\times {\displaystyle \underset{i=0}{\overset{k2}{}}}(x_ix_{2ki1},y_i+y_{2ki1},\dot{x}_i+\dot{x}_{2ki1},\dot{y}_i\dot{y}_{2ki1})`$
For $`N=4k+2`$ we set

$$R\left(\underset{i=0}{\overset{N1}{}}(x_i,y_i,\dot{x}_i,\dot{y}_i)\right)=\{v_k\}\times \underset{i=0}{\overset{k1}{}}(x_ix_{2ki},y_i+y_{2ki},\dot{x}_i+\dot{x}_{2ki},\dot{y}_i\dot{y}_{2ki})$$

We define map $`\mathrm{\Phi }:^{N1}E^1(\mathrm{\Omega })^{N1}`$ by

$$\mathrm{\Phi }=RPE.$$

###### Theorem 4.2.

If for some $`x^{N1}`$ $`\mathrm{\Phi }(x)=0`$, then there exists a trajectory with properties P1, P2, P3.

Proof: Let $`x`$ be point in which we have $`\mathrm{\Phi }(x)=0`$. Then $`E(x)`$ satisfies condition (1),(2) in Lemma 4.1 and from construction of $`\mathrm{\Phi }`$ we obtain that there exist trajectories ($`q_i(t)`$) starting with initial condition $`E(x)`$ and ending with $`\widehat{x}`$ such that $`R(\widehat{x})=0`$, hence condition (3) in Lemma 4.1 is satisfied. Now the assertion follows from Lemma 4.1. ∎

## 5. Existence of the SuperEight - Gerver orbit

The Gerver orbit (Figure 6, 7) is a choreography with 4 bodies forming a linear chain. It’s the simplest trajectory after the Eight.

With a computer assistance we proved the following

###### Theorem 5.1.

The Gerver SuperEight exists and is locally unique (up to obvious symmetries and rescaling).

We show the existence of SuperEight using an approach described in Section 4, i.e. we verify assumptions of Theorem 4.2. Below we give some details.

We set

(5.1)
$$E(x_1,\dot{x}_0,\dot{y}_1)=(0,a,\dot{x}_0,0,x_1,0,0,\dot{y}_1,0,a,\dot{x}_0,0,x_1,0,0,\dot{y}_1),$$

where $`a`$ is a parameter fixing the size of the orbit.

(5.2)
$$R(x_0,y_0,\dot{x}_0,\dot{y}_0,x_1,y_1,\dot{x}_1,\dot{y}_1,x_2,y_2,\dot{x}_2,\dot{y}_2,x_3,y_3,\dot{x}_3,\dot{y}_3)=(y_1+y_0,\dot{x}_1+\dot{x}_0,\dot{y}_1\dot{y}_0),$$

Poincáre section $`S`$ is defined by

(5.3)
$$S(x_0,y_0,\dot{x}_0,\dot{y}_0,x_1,y_1,\dot{x}_1,\dot{y}_1,x_2,y_2,\dot{x}_2,\dot{y}_2,x_3,y_3,\dot{x}_3,\dot{y}_3)=x_1x_0=0.$$

At first we have proved an existence of a zero of $`\mathrm{\Phi }(x)`$ using an interval Newton method, but in the paper we present data from proof based on Krawczyk method. In case of Gerver orbit the choice of the method isn’t so important, because if we take a time step small enough or smaller set $`[X]`$, then the computed matrix $`\left[\frac{\mathrm{\Phi }}{x}([X])\right]`$ becomes invertible and the proof usually goes through. But in general the main problem of application of an interval Newton method is that of an invertibility of $`\left[\frac{\mathrm{\Phi }}{x}([X])\right]`$. We avoid this using the Krawczyk method.

To check the assertion (2) in Theorem 2.2 we used the $`C^1`$-Lohner algorithm with order $`r=6`$. The time step was set to $`h=0.002`$. As matrix $`C`$ we used an inverse of the monodromy matrix computed nonrigorously in a point $`\overline{x}`$, i.e. $`C=\frac{\mathrm{\Phi }(\overline{x})}{x}^1`$.

Tables 5 and 6 contain numerical data from this proof.

## 6. Existence of the ’Linear chain’ orbit for the 6 bodies

Figure 8 displays a linear chain choreography with 6 bodies.

In this section we report about the computer assisted proof of the following

###### Theorem 6.1.

The linear chain for $`6`$ bodies exists and is locally unique (up to obvious symmetries and rescaling).

Again, we show an existence of this orbit using an approach described in Section 4 with some minor changes. To speed up calculation and to increase accuracy we take into account in computation that all time $`q_3(t)=q_0(t)`$, $`q_4(t)=q_1(t)`$ and $`q_5(t)=q_2(t)`$ (we change equation (1.1) doing a suitable substitution). So, our full space for ODE has now 12 dimensions. We use also a different time parameterization (we do time shift of $`\frac{1}{4}`$ of the period). To use exactly approach described in Section 4 we interchange axes (see Fig. 8 and 9). From Lemma 4.1 we will get, in this coordinate frame, a doubly symmetric periodic solution $`q(t)`$. It is easy to see that $`\overline{q}=q(t\frac{T}{4})`$ is then a solution sharing needed symmetries in the original coordinate frame. All data are given in frame with interchanged axes.

We set

(6.1)
$$E(\dot{x}_0,x_1,y_1,\dot{x}_1,\dot{y}_1)=(0,a,\dot{x}_0,0,x_1,y_1,\dot{x}_1,\dot{y}_1,x_1,y_1,\dot{x}_1,\dot{y}_1),$$

where $`a`$ is a parameter fixing the size of the orbit.

(6.2)
$$R(x_0,y_0,\dot{x}_0,\dot{y}_0,x_1,y_1,\dot{x}_1,\dot{y}_1,x_2,y_2,\dot{x}_2,\dot{y}_2)=(\dot{x}_1,x_0x_2,y_0+y_2,\dot{x}_0+\dot{x}_2,\dot{y}_0\dot{y}_2).$$

Poincáre section $`S`$ is defined by

(6.3)
$$S(x_0,y_0,\dot{x}_0,\dot{y}_0,x_1,y_1,\dot{x}_1,\dot{y}_1,x_2,y_2,\dot{x}_2,\dot{y}_2)=y_1=0.$$

To find a zero of $`\mathrm{\Phi }(x)`$ we use Krawczyk method. We check the assertion (2) in Theorem 2.2 using $`C^1`$-Lohner algorithm \[ZLo\] of order $`r=9`$ and a time step $`h=0.0025`$ for a computation in point $`\overline{x}`$ and $`h=0.001`$ for a computation on set $`[X]`$. As matrix $`C`$ we used an inverse of the monodromy matrix computed nonrigorously in a point $`\overline{x}`$, i.e. $`C=\frac{\mathrm{\Phi }(\overline{x})}{x}^1`$.

In Tables 7 and 8 we present data from computation of Krawczyk method.

## 7. Some technical data

All computations were performed on AMD Athlon 1700XP  with 256 MB DDRAM memory, with Windows 98SE operating system. We used CAPD package\[Capd\] and Borland C++ 5.02 compiler.

In the listing below $`r`$ is an order and $`h`$ is a time step used in the $`C^1`$-Lohner algorithm \[ZLo\].

The computation times for The Eight, $`h=0.01`$, $`r=7`$

* in point $`\overline{x}`$ : 1.417 sec
* for set $`[X]`$ : 2.66 sec
* convexity : 1.15501 sec

For the proof of an existence of Gerver solution in 4-body problem we used $`r=6`$, $`h=0.002`$. The computation times for both $`\overline{x}`$ and set $`[X]`$ were approximately equal to $`30.5`$ seconds.

For the proof of linear chain for 6-body problem

* computation of Poincaré Map for set $`[X]`$ took $`57.5`$ seconds with $`h=0.001`$ and $`r=9`$
* computation for $`\overline{x}`$ took $`23.8`$ seconds with $`h=0.0025`$ and $`r=9`$

The programm performing the proofs is available on http://www.ap.krakow.pl/~tkapela