text:
# The ELODIE survey for northern extra–solar planetsBased on observations made at the Haute–Provence Observatory (operated by French CNRS), the 1.2–m Euler swiss telescope at ESO–La Silla Observatory (Chile) and the 1.52–m ESO telescope also at La Silla Observatory. The ELODIE measurements discussed in this paper are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.125.5) or via http://cdsweb.u-strasbg.fr/Abstract.html

## 1 Introduction

The ELODIE Planet Search Survey, an extensive radial–velocity northern survey of dwarf stars, has been underway for several years at the Haute–Provence Observatory (OHP, CNRS, France) using the ELODIE high-precision fiber-fed echelle spectrograph (Baranne et al. 1996) mounted on the Cassegrain focus of the 1.93–m telescope. This survey was initiated in 1994 by M. Mayor and D. Queloz in order to detect very low-mass stellar companions and allowed them to discover the first extra–solar planet (51 Peg b) orbiting a solar-type star (Mayor & Queloz 1995). This discovery was followed by a number of subsequent detections of extra–solar planets with ELODIE. Some of these results have previously been published.

This survey is part of a large effort aiming at extra–solar planet search through radial–velocity measurements since several other planet surveys with a sensitivity $`K_1`$$`>`$ 10 m s-1 are underway elsewhere. A list of the most productive ones includes the Geneva Southern Planet Search Programme, carried out with the CORALIE spectrograph mounted on the 1.2–m Euler swiss telescope (Queloz et al. 2000b; Udry et al. 2000a) at ESO–La Silla observatory, the California & Carnegie Planet Search with the HAMILTON spectrograph at the Lick Observatory (Marcy & Butler 1992) and the HIRES spectrograph mounted on the 10–m Keck–1 telescope (Vogt et al. 1994) at the W. M. Keck Observatory, the G-Dwarf Planet Search also performed with the HIRES spectrograph at Keck Observatory (Latham 2000), the Anglo-Australian Planet Search with the UCLES echelle spectrograph mounted on the 3.92–m Anglo-Australian Telescope (Tinney et al. 2001), the AFOE Planet Search (Korzennik et al. 1998) using the AFOE spectrograph mounted on the 1.5–m telescope at the Whipple Observatory and the McDonald Planetary Search (Cochran & Hatzes 1994) using the McDonald Observatory 2.7–m telescope coude spectrograph.

115 planet candidates, including 11 planetary systems, with minimum masses lower than 18 MJup and down to sub-Saturnian planets (Marcy et al. 2000; Pepe et al. 2002), have been detected so far, with periods from slightly less than 3 days to about 15 years. This should permit us to obtain more robust conclusions about the mass function estimation, and notably its fast rise towards lower masses, the relative metal content of the parent stars and constraints on the so-called brown dwarf desert at separations of a few AU (Halbwachs et al. 2000) and on the emerging period desert between 10 and 70 days (Udry et al. 2002), crucial ingredients to the understanding of planet formation mechanisms.

This paper is the first of a series devoted to the discovery of extra–solar planet candidates in the framework of the ELODIE Planet Search Survey. It reports on the detection of six new candidates, with minimum planetary masses ($`m_2\mathrm{sin}i`$) in the range 2 to 11 Jovian masses. Five of these are orbiting F8 or G0 dwarf stars and one is found around a G5IV star. They are all but one long-period objects on little to moderately elongated orbits: three periods are found to be greater than 1500 days, the shortest one is 227.55 days and the orbital eccentricities are of the order of 0.5 for four objects and close to 0.3 for the two other ones. These discoveries, together with 3 others announced during the Scientific Frontiers in Research on extra–solar Planets conference in Washington (June 2002), bring the number of planets discovered with ELODIE up to 18.

The next section describes the sample of the ELODIE Planet Search Survey and summarizes the planet candidates it has allowed to detect up to now.

The basic properties of the stars harboring the six new companions are given in Sect. 3. Sect. 4 presents the radial–velocity data and the orbital solutions for the companions. The case of HD 13507, suspected for some time to harbor a planet but that actually has a low mass star companion, is exposed in Appendix A.

## 2 The ELODIE Planet Search Survey

The original sample of Mayor & Queloz (1995) contained 142 stars, out of which 51 Peg. The main selection criterion of these G and K dwarfs was their radial–velocity non–variability according to the anterior CORAVEL survey (Duquennoy et al. 1991; Duquennoy & Mayor 1991). The sample was largely modified in 1996 and later, so that the to–date survey sample size amounts up to 372 stars. Among the present sample, only 71 stars remain from the original list (active and/or too faint stars of this list were excluded) and 259 additional stars were included with the following criteria: $`m_\mathrm{V}`$$``$ 7.65, $`\delta `$$`>`$ 0 $`\mathrm{°}`$, $`v\mathrm{sin}i`$ (from CORAVEL) $`<`$ 4 km s-1, CORAVEL non–variability, spectral type from F8 to M0. Both lists are now forming the main programme with 330 objects.

In addition to these main programme stars, for follow–up or phase–coverage reasons based on the large access one can have to a 2 m–class telescope observing time, the sample was more recently complemented with 42 stars common with other programmes, among which 32 stars, with V magnitudes up to 9, were identified as variable by the Keck G–Dwarf Planet Search (Latham 2000). One star, HD 209458, was in our main programme and now also belongs to the collaborative effort with the G–Dwarf Planet Search.

The Geneva Southern Planet Search programme and the ELODIE Planet Search Survey have distinct stars samples, based on the declination, although the former may occasionally support the monitoring of a usually faint, low declination ELODIE programme star.

The planetary companion detection is based on the stellar radial–velocity reflex motion produced by gravitational interaction between companion and star. The ELODIE spectrograph, operating between 3850 and 6800 Å, has a spectral resolution of about 42 000. The planet survey makes use of its dual–fiber mode where a star and a reference thorium–argon calibration lamp are simultaneously observed to monitor and correct the instrumental drift and to achieve the required high precision. The signal–to–noise ratio per pixel for $`m_\mathrm{V}`$ = 7.65 stars is about 100 for a 15-minute observation under good atmospheric conditions at OHP. In order to achieve the most homogeneous observations, some exposures are duplicated in order to reach the target signal–to–noise. None, except for some special cases, are longer than 15 minutes in order that no differential instrumental drift between the two beams dominate the resultant accuracy. The total exposure time per source is dictated by the goal that the photon–noise be lower than the instrumental limit. The radial–velocity measurement is produced on–line by the ELODIE automatic reduction software which cross-correlates the observed stellar spectrum with a binary mask template (Baranne et al. 1996). Some improvements have been implemented in 2001 in those regions of the mask particularly rich in telluric absorption features (see for details Pepe et al. 2002) leading to a current instrument–limited precision better than 6 m s-1.

The improved precision is demonstrated by the distribution of the weighted radial-velocity rms measured for 111 constant stars in our sample (Fig. 3). We corrected the measured dispersions for the contribution of the photon noise. The instrumental error can be well estimated by the rms values of the low-dispersion tail of the histogram. The median of this distribution is 7.8 m s-1 and we find that 20% of these constant stars exhibit a velocity dispersion lower than 6 m s-1. We also show in Fig. 4 our data for 4 constant stars with dispersions lower than 6.5 m s-1 (photon noise included). Finally, the residuals to the orbital solution obtained with ELODIE for 70 Vir (Naef et al. 2003a) are 6.1 m s-1 (photon noise included). The photon noise error of the 70 Vir measurements is $`ϵ_{\mathrm{phot}}`$ = 2.8 m s-1. The instrumental error inferred from this long series of values is: $`ϵ_{\mathrm{instr}}`$ = (6.12$``$ 2.82)1/2 = 5.4 m s-1.

More recent improvements in the telluric absorption features treatment have allowed us to reduce systematic yearly velocity effects that were still visible for some stars in our sample. This is for example the case for HD 33636 where the residuals to the fitted orbit could be reduced by a factor of almost two.

A total of more than 9000 measurements, spread over observing periods that are programmed about once every month, have been performed since the beginning of the programme at OHP. All stars but a few ones have been observed more than 10 times and in average almost 20 times but those presenting indications of variations may have been observed more than 50 times.

Table 1 lists all the planet candidates discovered so far with ELODIE having $`m_2\mathrm{sin}i`$$``$ 11 MJup including the six presented in this paper and four others also recently announced. More massive companions, with minimum masses in the range 20–75 MJup, among which HD 127506, HD 174457 and HD 185414, have been discovered with ELODIE (Naef et al. 2000b, Naef et al. in prep). A peculiar case, HD 166435, has also emerged from the ELODIE Planet Search Survey observations. This star presents a chromospheric activity stable enough to mimic a short period Keplerian signature. Therefore, any radial–velocity variation potentially due to a short period companion is now first checked against such activity using line bisector variations, that have been shown to be a robust diagnostic by Queloz et al. (2001).

The list in Table 1 includes HD 209458 b (Mazeh et al. 2000), the only extra–solar planet for which transits were observed. Photometric transits were obtained (Charbonneau et al. 2000) with timing predictions computed from ELODIE monitoring and communicated to D. Charbonneau and collaborators. The first spectro-photometric transit was observed with ELODIE (Queloz et al. 2000a). The photometric and spectroscopic transits yielded notably the planet radius, the rotation direction with regard to that of its primary and the mean density (0.3 g cm-3), essential to assess the gaseous nature of the planets. Two planet candidates in Table 1 were found around stars not belonging to the ELODIE main programme: HD 178911 Bb, a planetary companion in a stellar triple system (Zucker et al. 2002), and HD 80606 b, the extra–solar planet having the highest orbital eccentricity known so far (Naef et al. 2001a). They are both belonging to the G-Dwarf Planet Search programme carried out at the Keck telescope. Also, Gl 876 b, which is presently in the closest extra–solar planetary system known, was discovered as a by-product of an M–dwarf multiplicity statistical study carried out with the OHPELODIE and ESOFEROS spectrographs and with the Canada–France–Hawaii Telescope (CFHT) PUE’O adaptive optics imaging system by a Grenoble–Geneva team (Delfosse et al. 1998), and independently by Marcy et al. (1998). Later, the presence of two planets on resonant orbits with periods in the ratio 2:1 that were masqueraded by a unique, more eccentric one, was recognized by Marcy et al. (2001).

Fig. 1 shows the Hertzsprung–Russell diagram of the 317 stars in the ELODIE Planet Search Survey main programme sample for which an HIPPARCOS entry exists. The six stars harboring planet candidates discovered with ELODIE, presented in this paper, are shown with distinctive symbols. $`T_{\mathrm{eff}}`$ is derived from the HIPPARCOS $`BV`$ by using the calibration of Flower (1996)111Note that the coefficients of the calibration printed in Flower (1996) are incorrect; we rather used the coefficients directly obtained from the author., except for the six planet-candidate hosts for which the effective temperatures listed in Table 2 have been used, and $`L`$ from HIPPARCOS apparent magnitude and parallax and the bolometric correction given by Flower’s calibration.

In Fig. 2, the histograms of $`m_\mathrm{V}`$ and $`BV`$ for the same stars as in Fig. 1, show a clear cut-off in magnitude dictated by the ELODIE performances and a range in spectral type with a pronounced bias toward G0–G5 dwarfs and a strict cut-off at the F8 spectral type.

## 3 Stellar properties

The six stars found to harbour planetary companions are HD 8574, HD 23596, HD 33636, HD 50554, HD 106252 and HD 190228. Their basic characteristics are summarized in Table 2.

Spectral types, apparent magnitudes, color indexes, parallaxes and proper motions were all taken from the HIPPARCOS catalogue (ESA 1997). The bolometric corrections of all stars were derived from $`BV`$ by using the calibration of Flower (1996).

The atmospheric parameters $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$, $`\xi _\mathrm{t}`$ and $`[`$Fe/H$`]`$ and the stellar masses $`M_{}`$ are from Santos et al. (2003). The projected stellar rotational velocities, $`v\mathrm{sin}i`$, were obtained using the mean ELODIE cross-correlation functions (CCF) dip width (Baranne et al. 1996), except for HD 106252, and the calibration by Queloz et al. (1998). For HD 106252, $`v\mathrm{sin}i`$ has been derived from CCFs obtained with CORALIE and the calibration in Santos et al. (2002).

Table 2 also gives for each star the equivalent width $`W_{\lambda ,\mathrm{Li}}`$ of the $`\lambda 6707.8`$ ÅLi i line, as measured on the co-addition of all the available ELODIE spectra, corrected for the spectrograph straight light. It can be seen that all stars exhibit a lithium absorption feature which is rather strong except for two of them, HD 106252 and HD 190228. Lithium abundances, scaled to $`\mathrm{log}n(\mathrm{H})`$ = 12, were derived following the method described in Naef et al. (2001b) by using the curves of growth published by Soderblom et al. (1993) and our tabulated $`T_{\mathrm{eff}}`$. The errors were estimated by assuming temperature and equivalent width variations of $`\pm `$ 1$`\sigma `$. The maximum abundance values are obtained by using $`T_{\mathrm{eff}}`$ \+ 1$`\sigma `$ and $`W_{\lambda ,\mathrm{Li}}`$ \+ 1$`\sigma `$ and the minimum values correspond to $`T_{\mathrm{eff}}`$$``$ 1$`\sigma `$ and $`W_{\lambda ,\mathrm{Li}}`$$``$ 1$`\sigma `$. The resulting errors are not symmetrical because the error on $`W_{\lambda ,\mathrm{Li}}`$ is supposed to be symmetrical and equal to $`\pm `$ 5 mÅ.

The Ca ii H absorption feature was extracted for the six stars as a chromospheric activity diagnostic. Fig. 5 displays the spectra of the Ca ii H line region obtained by the co-addition of all ELODIE spectra available for each star. None of the stars shows an appreciable emission. However, in this region, the ELODIE spectra do not usually have a high signal–to–noise ratio preventing us from the extraction of meaningful values of the $`\mathrm{log}(R_{HK}^{^{}})`$ indicator. The apparent absence of emission does therefore not preclude a mild activity.

## 4 Orbital solutions

The achieved precision (photon noise + instrumental error) of the ELODIE radial–velocity measurements, although improved since the beginning of the ELODIE Planet Search Survey (see Sect. 2), is of the order of 10 m s-1 for these targets over the time span of the reported observations. The 254 individual high-precision ELODIE radial velocities (+5 CORALIE measurements) obtained for the six planet-candidate host stars are available in electronic form at the CDS in Strasbourg. All these measurements are displayed in Fig. 6 and 7 which also show the orbital solutions, resulting from the best Keplerian fit and described below for each star.

We checked all solutions against false-alarms with Fourier transform analysis and found negligible (below the 1% level) false-alarm probabilities. We also checked their fit errors with Monte–Carlo simulations provided within the ORBIT programme described in Forveille et al. (1999). To this end, for each orbital solution, we ran 10 000 trials to derive Monte–Carlo estimates of the 1 and 3 $`\sigma `$ confidence intervals. None of the solution displays a noticeable discrepancy with these simulations except in three cases for a few parameters, as mentioned in the text below.

### 4.1 HD 8574

We obtained a total of 41 ELODIE high-precision radial–velocity measurements for HD 8574 which was observed since Jan. 11, 1998 (HJD = 2 450 825). We list in Table 3 the fitted orbital elements to these measurements together with the computed minimum mass $`m_2\mathrm{sin}i`$ = 2.11 MJup of the planetary companion, assuming a primary star mass of 1.17 M. Our discovery of a companion of $`m_2\mathrm{sin}i`$ = 1.95 MJup had been announced by the April 4th ESO Press Release222www.eso.org/outreach/press-rel/pr-2001/pr-07-01.html. The orbital solution presented here results from additional measurements obtained since then which further improved it. The errors listed in Table 3 are derived from Monte–Carlo simulations. They do not differ appreciably from the fit errors assessing the quality of this solution despite the rather high value of $`\chi _{\mathrm{red}}^2`$.

The companion of HD 8574 has the lowest period (227.55 d), well constrained (0.3%) and the best covered in terms of number of orbital cycles of the six planet candidates. Accordingly, it has also the smallest semi-major axis ($`a`$ = 0.77 AU, i.e. 17 mas). The slightly abnormal residuals from the fitted orbit yield a moderate ”jitter” of 8.2 m s-1, although the primary has no known activity and does not exhibit any clear one from the Ca ii H line core (Fig. 5).

### 4.2 HD 23596

We obtained a total of 39 ELODIE high-precision radial–velocity measurements for HD 23596 which was observed since Jan. 10, 1998 (HJD = 2 450 824). Table 3 gives the obtained orbital elements together with the minimum mass of the companion ($`m_2\mathrm{sin}i`$ = 8.10 MJup assuming $`M_{}`$ = 1.30 M). The computed semi-major axis of the planetary orbit is $`a`$ = 2.88 AU, i.e. 55 mas. This solution seems to be well constrained as judged from the absence of abnormal residuals and the agreement between fit errors and Monte–Carlo simulations.

The planet of HD 23596 has the largest semi-major axis of the six candidate planets if we put aside the peculiar case of HD 33636. With a period of 4.28 yr, the candidate around HD 23596 belongs to the category of long period planets which is naturally growing with the time span of the on-going surveys. It is also a massive planet. We emphasize the metal-rich nature of the parent star, a characteristics not shared by any of the six other stars.

### 4.3 HD 33636

The fitted orbital elements for HD 33636 to its 42 ELODIE radial–velocity measurements obtained since Feb. 18, 1998 (HJD = 2 450 863), plus 5 recent CORALIE measurements, are presented in Table 3. The velocity offset between the two instruments is an additional free parameter. Its fitted value is $`\mathrm{\Delta }RV_{\mathrm{E}\mathrm{C}}`$ = 0.038 $`\pm `$ 0.004 km s-1. The variability of the HD 33636 radial velocities has been evident for several years since the first measurements were done while the planet was at its periastron. With a primary mass of 1.12 M given by Santos et al. (2003), we compute a minimum mass of $`m_2\mathrm{sin}i`$ = 10.58 MJup and a semi-major axis of $`a`$ = 4.08 AU, i.e. 142 mas, for the planet. It must be emphasized that the period is certainly not yet covered (see Fig. 7), a fact corroborated by the change of the fitted period with the addition of the measurements of our last observing campaign. At this stage, most of the constraints come from the impact of the eccentricity on the radial velocity curve shape. Therefore the period of 7.74 yr is presently badly constrained but other parameters should be less affected. Interestingly, Monte–Carlo simulations show that the fit errors (listed in Table 3) are somewhat overestimated. We adopt the asymmetric Monte-Carlo errors for the period and keep the fit errors for the other parameters.

According to our preliminary orbital solution, HD 33636 b has the largest semi-major axis (4.08 AU) of the six planets. Despite its uncertainty, the period is long and therefore this new planet also belongs to the growing category of long period planets. It is also be the most massive of our six candidates.

Vogt et al. (2002) have announced the independent detection of a companion of $`m_2\mathrm{sin}i`$ = 7.7 MJup. We have superimposed their orbital solution to ours in Fig. 7. The discrepancy between the two solutions is clearly due to their largely underestimated period. Our last campaign provides several measurements of good quality with which their orbital curve departs by 6 to 8 standard deviations of our most recent measurements. Moreover, the period in Vogt et al. (2002) is shorter than the time span of our measurements.

### 4.4 HD 50554

A total of 41 ELODIE precise radial–velocity measurements were obtained for HD 50554 since Nov. 16, 1997 (HJD = 2 450 769). We list in Table 3 the fitted orbital elements together with the computed minimum mass $`m_2\mathrm{sin}i`$ = 5.16 MJup, assuming a primary star mass of 1.11 M. The confidence interval resulting from Monte–Carlo simulations is a bit larger than the fit error for the period. This is probably due to the peculiar weight of the first season measurement in the period estimation (see Fig. 6). We therefore took the confidence interval for the period instead of its fit error. The semi-major axis of 2.41 AU, i.e. 77 mas, is derived for this new planetary companion. The mean radial–velocity measurement error $`ϵ_{\mathrm{RV}}`$ and the weighted rms radial–velocity deviation $`\sigma (\mathrm{O}\mathrm{C})`$ differ notably, yielding an appreciable ”jitter” of 6.3 m s-1.

Our discovery of a companion of $`m_2\mathrm{sin}i`$ = 4.9 MJup had been announced by the ESO Press Release nb 07/01 (April 4, 2001). The orbital solution presented here results from additional measurements obtained since then which further improved it.

Fischer et al. (2002) have confirmed the detection of a companion with $`m_2\mathrm{sin}i`$ = 3.7 MJup. Except for their $`K_1`$ and $`\omega `$ values that differ significantly from ours, the two orbital solutions are compatible. The better coverage of the radial–velocity maximum by our measurements probably explains the difference. The $`m_2\mathrm{sin}i`$ discrepancy results from the rather large $`K_1`$ difference between the two solutions.

The Fischer et al. (2002) value of the chromospheric activity tracer $`\mathrm{log}(R_{HK}^{^{}})`$ ($``$4.94), indicating a low activity for HD 50554, is not useful for explaining the observed ”jitter”, specially since the level of ”jitter” is not well correlated to the $`\mathrm{log}(R_{HK}^{^{}})`$ activity indicator for low-activity stars (Saar et al. 1998; Santos et al. 2000b).

### 4.5 HD 106252

We gathered a total of 40 ELODIE high-precision radial–velocity measurements for HD 106252 which was observed since March 1, 1997 (HJD = 2 450 509). In Table 3 are presented the fitted orbital elements to our measurements. With a primary mass of $`M_{}`$ =1.02 M, we compute a minimum mass of $`m_2\mathrm{sin}i`$ = 7.56 MJup and a semi-major axis of 2.70 AU, i.e. 72 mas, for HD 106252 b. This solution is well constrained as confirmed by the absence of abnormal residuals. The Monte-Carlo confidence intervals are in good agreement with the fit errors.

The discovery of a companion with $`m_2\mathrm{sin}i`$ = 6.8 MJup had also been announced by the ESO Press Release nb 07/01 (April 4, 2001). The orbital solution presented here results from additional measurements obtained since then which further improved it. With a period of 4.38 yr, this new planet also belongs to the growing category of long period planets. It is also a rather massive candidate.

Fischer et al. (2002) have confirmed the detection of a companion with $`m_2\mathrm{sin}i`$ = 6.96 MJup, a value in rather good agreement with ours. Their period is at more than 3$`\sigma `$ from our value. This difference probably results in their poorer temporal coverage.

The Fischer et al. (2002) value of the chromospheric tracer $`\mathrm{log}(R_{HK}^{^{}})`$ ($``$4.97), indicating a low activity for HD 106252, agrees with the absence of abnormal residuals.

### 4.6 HD 190228

We obtained a total of 51 ELODIE high-precision radial–velocity measurements for HD 190228 which was observed since Jul. 22, 1997 (HJD = 2 450 652). We list in Table 3 the fitted orbital elements to these measurements together with the computed minimum mass $`m_2\mathrm{sin}i`$ = 3.58 MJup assuming the primary star mass of 0.83 M given by Santos et al. (2003). The semi-major axis $`a`$ = 2.02 AU, i.e. 33 mas, is given as well. The confidence interval resulting from Monte–Carlo simulations is somewhat larger (indeed asymmetrically) than the fit error for the eccentricity and $`K_1`$. This could be caused by the bad coverage of the periastron (see Fig. 6). We therefore took the confidence interval for those two parameters instead of their fit error.

The discovery of a companion with $`m_2\mathrm{sin}i`$ = 3.4 MJup had been announced in Sivan et al. (2000). The orbital solution presented here results from additional measurements obtained since then which slightly improved it.

Among the previous discoveries, more than a dozen planet candidates orbit sub-giant stars i.e. almost 20% of the total. Not yet statistically meaningful, this ratio however appears to be much larger than the proportion of sub-giant stars in the survey sample (19 stars among the ELODIE Planet Search Survey sample i.e. 5% of the total). The planet of HD 190228 is therefore not peculiar in this respect. On the contrary, HD 190228 is peculiar in terms of metallicity since it is the unique case in this paper sample with a metallicity value ($``$0.24 $`\pm `$ 0.06) quite below the solar abundance. It is also peculiar with respect to its spectral sub-type, since it is the coolest star of this sample.

## 5 Discussion

### 5.1 Existence of the companions

The orbital periods derived for the six planetary companions, spanning from 0.62 to 4.38 yr for five of them and even more for HD 33636 whose period is fairly imprecise, are much longer than the possible rotation periods of the primary stars, even though we have no individual estimation for most of them. This makes it very unlikely that any of the radial–velocity curves might be produced by activity spot phenomena as it has been observed for the short period variable HD 166435 (Queloz et al. 2001). Nevertheless we checked this for each star by computing the line bisectors of the observed CCFs and the linear correlation between the radial–velocity measurements and the bisector span. With maximum absolute values of the correlation coefficient below 0.23, we found no significant correlation, down to the 1$`\sigma `$ level, for any of the six stars. The interpretation of the Doppler signature by the presence of a companion thus casts no doubt for the six planet candidates, even for the mildly active stars, a conclusion in agreement with that of Naef et al. (2000a) about long-period planet detectability around active stars.

### 5.2 $`m_2\mathrm{sin}i`$ distribution and mass function

With six new planet candidates, we can usefully update the distribution of $`m_2\mathrm{sin}i`$. Fig. 8 displays the histogram of $`m_2\mathrm{sin}i`$ of the 115 objects known so far, with $`m_2\mathrm{sin}i`$$`<`$ 18 MJup, including our six new candidates and about thirty other new candidates announced very recently. These additional objects, half of which below $`m_2\mathrm{sin}i`$ = 2 MJup globally tend to smooth the cumulative frequency distribution and to fill some gaps present about $`m_2\mathrm{sin}i`$ = 8 MJup (see Fig. 1 of Jorissen et al. (2001)). Because the techniques used to estimate the deconvolved $`M_2`$ distribution (e.g. Jorissen et al. 2001) is sensitive to the sampling quality, one can infer that such an approach will be even more robust with the new, larger input frequency distribution. More specifically, the latter is likely to result in a smoother $`M_2`$ distribution, specially around and above 10 MJup that could make it less obvious that 10 MJup would be a quasi upper $`M_2`$ limit. Since this is only qualitative, it will be very interesting to know whether this study, applied to the new extra–solar planet sample, would rather lead to a vanishing distribution in the 15 MJup region.

### 5.3 Eccentricity versus period diagram

As seen above, the periods of all the companions are long, exceeding three years except for HD 8574. Their orbital eccentricities are close to 0.5 for four of the objects and about 0.3 for the two other ones. We have placed these six companions on a $`e`$ vs $`P`$ diagram together with all the other known companions with $`m_2\mathrm{sin}i`$$`<`$ 18 MJup (Fig. 9). Their location is dictated by the large values of the periods which is the natural effect of the bias towards long periods while the time span of the observing programmes increases. Apart from this, the diagram shows no particular trend for the six companions which are spread in a region of intermediate eccentricities, already well populated. The locations of our six new companions are in agreement with, and reinforce, the general behavior that tends to reproduce the spectroscopic binaries $`e`$ vs $`P`$ distribution for long periods (Mayor et al. 2001; Mazeh & Zucker 2001), a rather paradoxical constraint on planet formation mechanisms.

### 5.4 Metallicity

It has been shown that the metallicity distribution of the planet-bearing stars is statistically 0.25 dex higher than for solar neighborhood stars (e.g. Santos et al. 2001). Since four of the six stars with new companions discussed here have solar metallicity, they will not help to statistically improve this conclusion. Moreover, while HD 23596 is over-metallic by 0.32 dex and therefore clearly falls in the upper part of the planet-bearing star distribution of Santos et al. (2001), HD 190228 is under-metallic by a similar amount ($``$0.24 dex). Both stars therefore almost compensates each other’s effect on the distribution. The six stars will rather serve as additional material for further studies of metallicity distribution based on larger statistical samples.

### 5.5 Direct imaging feasibility

With separations at the orbital elongation larger than 60 mas, four out of the six new planet candidates could in principle be targets for diffraction–limited, i.e. adaptive optics, K band near infrared imaging on 8 m–class telescopes. The feasibility of such measurements is of course strongly dependent on the achievable contrast and therefore very questionable for long periods, cool Jupiter–like planets. Such imaging would anyway place useful constraints on the maximum mass and the orbital inclination of the companions. Altogether, the objective to image such companions clearly calls for very high contrast dedicated instruments (see e.g. Mouillet et al. 2001).

###### Acknowledgements.

We acknowledge support from the Swiss National Research Found (FNRS), the Geneva University and the French CNRS. We are grateful to the Observatoire de Haute-Provence, the CFGT french time allocation committee and the Programme National de Planétologie of CNRS for the generous time allocation and their continuous support to this long–term project. We thank T. Forveille and D. Ségransan for the adaptive optics observation and data treatment. T. Forveille also helped for the Monte–Carlo simulations with his ORBIT programme. We are very much indebted to the Observatoire de Haute-Provence technical and telescope support team for the years-long dedication to the ELODIE operation, specially to A. Vin for his critical software support, and to the various night assistants who have taken an essential part in the very demanding telescope operation of this observing programme. We are grateful to the referee for his comments to the original version of the paper. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

## Appendix A HD 13507

According to the HIPPARCOS catalogue, HD 13507 is a $`m_\mathrm{V}`$ = 7.12 G0 dwarf star with a parallax of 38.12 $`\pm `$ 0.89 mas but its absolute magnitude $`M_\mathrm{V}`$=5.10 $`BV`$ (0.672$`\pm `$ 0.007) rather suggest a G4 spectral type.

The atmospheric parameters and the mass of this star have been derived from a high signal–to–noise ELODIE spectrum by N.C. Santos (priv. comm.) in the way described in Santos et al. (2000a): $`[`$Fe/H$`]`$ = 0.01$`\pm `$ 0.08, $`T_{\mathrm{eff}}`$ = 5775$`\pm `$ 80 K, $`\mathrm{log}g`$ = 4.71 $`\pm `$ 0.17 (cgs), $`\xi _\mathrm{t}`$ = 1.33$`\pm `$ 0.11 km s-1 and $`M_{}`$ = 1.09 M.

A total of 19 ELODIE precise radial–velocity measurements were gathered for HD 13507 since Jan. 8, 1998 (HJD = 2 450 822). A preliminary but apparently good solution, based on a linear drift of 341 $`\pm `$ 8 m s-1 yr-1 and a Keplerian variation compatible with a planetary companion, that was announced as such, could be derived until mid 2002.

The later measurements obtained between july and december 2002 (see Fig. 10) invalidated this interpretation and revealed instead a classical spectroscopic binary velocity curve, caused by a low mass star companion. This possibly explains the noticed discrepancy between the spectral type from HIPPARCOS and its photometry. The computed orbital parameters obtained from a still preliminary fit are: period $`P`$$``$ 3000 days, $`e`$ = 0.14, $`K_1`$ = 0.694 m s-1, minimum mass of the companion $`m_2\mathrm{sin}i`$ = 52 MJup, assuming a primary star mass of 1.09 M, and semi-major axis $`a`$ = 4.3 AU, i.e. 164 mas. This solution is superimposed to the measurements, displayed in Fig. 10.

We have searched for a companion to HD 13507 by adaptive optics imaging. An image was obtained in Bracket$`\gamma `$ with PUE’O at the CFHT on August, $`7^{th}`$ 2001 (HJD = 2 452 128.61) but no companion shows up at a 3$`\sigma `$ level. A unique negative measurement leaves place to a possible alignment of the components but monitoring the binary with high dynamics imaging over a few years should resolve it and permit to derive the companion mass, a useful information for studies of sub–stellar objects statistics.


format: text

text:
# A brief (blazar oriented) overview on topics for multi-wavelength observations with TeV photons

## 1 Introduction

With the upcoming of the new generation of Čerenkov Telescopes (CT, e.g. H.E.S.S., VERITAS, CANGAROO III, MAGIC), TeV astronomy is about to greatly improve its scientific capabilities. Multi-wavelength observations are an essential tool to investigate how TeV sources work. Below, I will briefly present some recent topics and perspectives for multi-wavelength observations with TeV photons, with main emphasis on blazars since they represent nearly all the extragalactic sources so far detected. The scheme is the following: at first I will very briefly recall different mechanisms to produce TeV photons, to identify the most interesting frequencies to accompany TeV observations. Then I will discuss the objects’ main variability characteristics, which determine the observation strategies, and will focus on some issues related to the simultaneous coverage.

## 2 TeV emission ingredients

The first obvious requirement to produce TeV photons is to have particles with TeV energies, for energy conservation. According to which particle is the main carrier of this energy, the recipes for TeV emission can be divided in two classes. In the so called “hadronic” models, protons play the main role, producing TeV photons through hadronic interactions (as $`pp\pi ^02\gamma `$ pohl , or $`p\gamma \pi \gamma e^+e^{}`$ with a subsequent electromagnetic cascade mann ), or directly through synchrotron radiation in strong magnetic fields (B$``$100G), with emission at lower frequencies contributed also by secondary pair-produced electrons ahap . For these models the overall spectral energy distribution (SED) is determined by the details of the cascade or the secondary electrons production, as well as the initial particle spectrum and source conditions, so that basically all lower frequencies are equally interesting to distinguish among them and to constrain the physical parameters.

In the “leptonic” models, instead, electrons are the main particles, producing TeV photons through the inverse Compton (IC) process (more economic than the synchrotron one). Allowing for beaming effects, to give origin to an observed TeV photon (i.e. a TeV/$`\delta `$ photon in the comoving frame), an electron needs at least an energy $`\gamma m_ec^2\mathrm{TeV}/\delta `$ (i.e. $`\gamma 10^5`$ for $`\delta 10`$). The seed photons most effective to be up-scattered by these electrons are of frequency $`h\nu m_ec^2/\gamma `$ (in the co-moving frame), since above they scatter in the Klein-Nishina regime. In the presence of magnetic field, these same electrons will also produce synchrotron radiation of observed energy $`h\nu =1.5B(\mathrm{Gauss})(\gamma /10^5)^2(\delta /10)`$ keV. In this scenario, therefore, two other frequency bands stand out in a natural way: the X-ray one, which can provide number density and energy distribution of the electrons also emitting in the TeV band, and the “optical” one (in a broad sense, i.e. from UV to IR), where we can have information about the most effective seed photons. Multi-wavelength observations in these three bands represent therefore a powerful diagnostic tool for the acceleration mechanism and the details of the high energy emission (as well as a test for the leptonic scenario itself).

Note anyway that, whatever scenario is adopted, the conditions must be compatible with the TeV photons survival (since we detect them), i.e. the source must be transparent with respect to the $`\gamma `$-$`\gamma `$ absorption.

## 3 Variability properties

Multi-wavelength observations of variable sources require an understanding of their variability characteristics to decide the best observing strategy. The most variable and potentially strong extragalactic sources of TeV photons are Gamma-Ray Bursts (GRBs) and blazars. For GRBs (whose flux decays as $`t^1`$) the choice is rather straightforward (albeit difficult to carry out): to observe as quick as possible and until it fades away.

Blazars instead are characterized by very different variability timescales, from years to minutes, and not only among different objects but also in the single sources. In high energy peaked BL Lacs (HBLs, which comprise all the detected TeV sources) we are now aware of at least 3 main timescales, well exemplified by the Mkn 501 behaviour in X-rays (Fig. 1): a long-term one, from months to years (e.g. the 1997 huge flare); a “middle-term” one, with a typical timescale of $``$1 day (as revealed by the ASCA long observations, also in Mkn 421 and PKS 2155-304 chiharu ); and a short-term one, from hours down to even minutes, as observed by RXTE in May 1998, with a hard X-ray doubling timescale of $``$6 minutes catanese . Extremely short variability is also a property of the TeV emission (e.g. Mkn 421, gaidos ; aha2 ), with doubling timescales as short as 20 minutes.

Connected with the flux variability, HBLs show also strong spectral variations, which often change the position of the SED peaks. Although there seems to be a trend (in the X-ray and TeV bands) for “harder when brighter” spectra, there is not a unique relation between flares and (synchrotron) peak changes. In fact, we can outline 3 types of behaviour: a “Mkn421–like” one (Fig. 2 left), where the peak frequency changes very little despite large flux variations gfoss ; a “Mkn501–like” one (Fig. 2 right), where during flares the peak shifts towards higher energies by orders of magnitude (e.g from $`<`$1 to $``$100 keV); and a “1ES 1426+428–like” one (Fig. 3), which seems characterized by large peak variations (of the same magnitude of Mkn501) but during roughly constant flux states (at least in the X-ray band). A “Mkn501–like” behaviour has been recently followed by 1ES 1959+650 during the May 2002 strong flare henric , with striking similarities also in the ’pivoting’ of the SED around the optical–soft-X range horn19 .

All these differences and similarities suggest that there are common but intrinsically different parameters/mechanisms driving the variability, whose understanding requires to follow the overall SED changes on all typical timescales.

## 4 The X-ray–TeV connection

The interest of X-ray and TeV observations is that they sample electrons of roughly equal energy, emitting through two different processes. In HBLs, moreover, the bulk of the luminosity is emitted in those bands and it appears to be produced by the same electron population, as suggested by the often observed tight correlation between the two emissions. A strong confirmation of this last hypothesis has come from a simultaneous observation of Mkn 421 in 1998 (Fig. 4), where the peaks of the emissions were correlated within one hour (although with different halving times; maraschi ).

However, new results are now showing that this relation may be more complex than expected, since the two bands does not always follow the same behaviour. An example is given again by Mkn 421, in 2001 (Fig. 4 right): although there is a general correlation, there are periods where we see a strong decrease (MJD $``$51989) or a rapid flare (MJD $``$51990) of the TeV flux not accompanied by relevant changes in the X-ray band. This behaviour seems also confirmed by 1ES 1959+650 in its 2002 active state, where evidence was found for a so called “orphan” $`\gamma `$-ray flare (i.e. not accompanied by an X-ray flare) henric . More observations, possibly over several days, are necessary, to see if there are also X-ray flares not accompanied by TeV flares, and to unveil possible lags between them and with the seed photons fluxes.

For both diagnostics and model testing, the intra-day activity of TeV blazars is perhaps the most challenging, since the shortest TeV timescales pose the tightest constraints to the size of the emitting region and to the $`\gamma `$-$`\gamma `$ transparency celotti . Moreover, recent HEGRA data on Mkn 421 have shown an entire “zoo” of flares with quite different rise and decay times (e.g. Fig. 5), hint of a complex interplay between acceleration and cooling/escape times. The improvement in sensitivity ($``$1 order of magnitude) of the new generation CT, coupled with the continuous coverage provided by the large area X-ray telescopes (CHANDRA, XMM), are now making possible to follow the spectral evolution in both bands on hour timescales and without gaps, allowing to study the inner jet conditions with unprecedented detail.

Another fundamental advantage of the new generation CT is their lower energy thresholds ($``$150 GeV), which will provide data also at energies not (or less) affected by absorption due to $`\gamma `$-$`\gamma `$ collisions with the cosmic Extragalactic Background Light (EBL, steck ; primack ). Together with the study of the correlated X-ray–TeV variability, from simultaneous observations we hope to be able to constrain the whole SED model parameters such to “measure” the EBL by comparing the predicted intrinsic TeV spectrum with the observed one coppi ; kraw .

## 5 Observation strategies

It is clear at this point that we need simultaneous observations in all the three main bands with spectral information. Spectra are precious also in the “optical” band, to outline the overall SED shape (and peak position) and to follow its evolution during flares. The “optical” band is important also because changes in the seed photons flux or/and spectrum may give origin to different behaviours between the TeV and X-ray emissions.

The ’long-term’ variability requires of course a long term monitoring program of a number of objects (as now done by RXTE ASM), which should continue to be done preferably at X-rays, since an only optical monitor would have missed the “Mkn501-like” flares, due to the pivoting of the SED at low frequencies. Simultaneous campaigns over several days, maximizing the continuous coverage through the coordination of ground CT and optical telescopes, will allow instead to sample both the middle and short-term variabilities, and to reveal possible lags. Long observations will allow also to investigate the relations between these two variabilities, e.g. if the shortest one appears only on top of the other or not, and if with different characteristics. For the same reason, simultaneous campaigns should be performed not only during the flares detected by the long term monitoring, but also in quiescent periods.

## 6 A warning about X-ray coverage

Together with CHANDRA, XMM is in principle the best X-ray instrument for this type of studies, thanks to its large collecting area (providing spectra on shorter timescales) and its long orbital period, which allows for continuous observations without gaps up to $``$37 hrs. Regrettably, the strong constraints on the solar panels orientation (the sun avoidance angle has to be $`70^{}\theta 110^{}`$) severely affects the satellite performance in two for us important areas: the ability to perform observations triggered by flaring states and simultaneously with ground telescopes. As shown in Fig. 6 (left), the visibility of an object during the year is limited for most of the sky to $``$2-3 months, with correspondingly low probability that a flare will occur in the right epoch. Even more important, the constraints involve also the anti-solar direction (contrary to CHANDRA), thus making impossible to observe a source when it is visible under the best conditions from ground telescopes (i.e. for the full night).

An example is given by PKS 2155-304 with H.E.S.S. (Fig. 6, right): the source can be pointed by XMM only in May (28/4–27/5) and November (31/10–26/11), i.e. when it can be observed by H.E.S.S. for only $`2`$ hours at the end or beginning of the night. The same problems affect also INTEGRAL, although with slightly less constraints on the sun angle and thanks to the wide field of view of its instruments (JEM-X f.o.v. is 4.8, fully coded). Since SWIFT will be devoted mainly to GRBs studies (but see gehrels  for its possible use also with flaring blazars), this leaves the TeV community for the next years (i.e. when the new generation CTs are becoming operative) with only CHANDRA and RXTE as the main X-ray satellites for simultaneous campaigns.


text:
# Scattering of periodic solitons

## 1 Introduction

Physics in (2+1) dimensions is an area of much active research, covering topics that include Heisenberg ferromagnets, the quantum Hall effect, superconductivity, nematic crystals, topological fluids, vortices and solitary waves . Most of these systems are non-linear. In their mathematical description the well-known family of sigma models plays a starring role. The simplest Lorentz-covariant soliton model in (2+1) dimensions is the $`CP^1`$ or non-linear $`O(3)`$ sigma model. Its solutions, sometimes called ‘lumps’, are realisations of harmonic maps, a long-established area of research in pure mathematics. However, analytical $`O(3)`$ solitons solutions have only been found for the static case; the full time-dependent model must be studied using numerical methods and/or other approximation procedures .

Sigma models are also useful as low dimensional analogues of field theories in higher dimensions. In effect, the $`O(3)`$ model in two dimensions exhibits conformal invariance, spontaneous symmetry breaking, asymptotic freedom and topological solitons, properties similar to those present in a number of important field theories in (3+1) dimensions –like the Skyrme model of nuclear physics .

We are concerned with the planar $`CP^1`$ model (both in its original and Skyrme-like versions) with periodic boundary conditions where the solitons are harmonic maps $`T_2S_2`$. A rich diversity of phenomena has been found in this model , going beyond the two-soliton and annular structures one might expect by analogy with the model with standard boundary conditions, where the soliton fields are harmonic maps $`\mathrm{}_2\{\mathrm{}\}S_2S_2`$.

For planar systems, $`N`$ identical lumps initially placed at the vertices of a rectangular $`N`$-gon, fired with equal speed to collide head-on at the centre of the polygon, scatter and emerge on the vertices of the dual polygon. Such $`\pi /N`$ scattering, studied for the $`CP^1`$ model on $`\mathrm{}_2`$ in reference , may be understood on symmetry grounds: the initial data has $`D_N`$ (dihedral group) symmetry and the time evolution respects that. But imposing periodic boundary conditions breaks the foresaid symmetry, and the interesting question whether dual-polygon scattering still holds for this case must be investigated.

Through numerical simulations on $`T_2`$ we confirm the $`\pi /N`$ scattering for 3 and 4 identical lumps. The case $`N=2`$ has being considered elsewhere  and it conforms to the well-documented scattering at 90. We also look into non-symmetrical configurations (they do not scatter at $`\pi /N`$) and explain the results using the geodesic approximation.

In the following section we introduce our periodic $`CP^1`$ model. The numerical procedure is explained in section 3. Section 4 analyses collisions between three solitons, and the case $`N=4`$ is considered in section 5. In section 6 we present our version of the geodesic approximation which includes the case of initial configurations that are not symmetrical. We find that our predictions (based on this approximation) are in agreement with what is observed in numerical simulations. We close with some concluding remarks in section 7.

## 2 The $`CP^1`$ model on the torus

The non-linear $`O(3)`$ model involves three real scalar fields $`\stackrel{}{\varphi }(x^\mu )(\varphi _1,\varphi _2,\varphi _3)`$ which satisfy the constraint that $``$ $`x^\mu (x^0,x^1,x^2)=(t,x,y)T_2`$ the fields lie on the unit sphere $`S_2`$:

$$\stackrel{}{\varphi }.\stackrel{}{\varphi }=1.$$
(1)
Subject to this constrain, the Lagrangian density of the system is

$$=\frac{1}{4}(_\mu \stackrel{}{\varphi }).(^\mu \stackrel{}{\varphi }).$$
(2)
The invariance of the model under global $`O(3)`$ rotations in field-space is apparent.

The model is conveniently recast in terms of one independent complex field $`W`$ (the $`CP^1`$ formulation) related to $`\stackrel{}{\varphi }`$ via

$$\stackrel{}{\varphi }=(\underset{\varphi _1}{\underset{}{\frac{W+\overline{W}}{|W|^2+1}}},\underset{\varphi _2}{\underset{}{i\frac{W+\overline{W}}{|W|^2+1}}},\underset{\varphi _3}{\underset{}{\frac{|W|^21}{|W|^2+1}}}).$$
(3)
The Lagrangian (2) now reads

$$=\frac{|_tW|^22|_zW|^22|_{\overline{z}}W|^2}{(1+|W|^2)^2},$$
(4)
where $`z=x+iy`$ and $`\overline{z}=xiy`$. For all $`t`$, the fields $`W`$ are mappings from the torus $`T_2`$ to the sphere $`S_2`$, ie, they satisfy the periodic boundary conditions

$$W(z+mL+inL)=W(z),$$
(5)
where $`m,n=0,1,2,\mathrm{}`$ and the period $`L`$ denotes the size of the square torus.

The *static* soliton solutions of the model are thus doubly periodic functions of $`z`$, that is, elliptic functions that may be expressed through Weierstrass $`\sigma `$-function as

$$W=\lambda \underset{j=1}{\overset{\kappa }{}}\frac{\sigma (za_j)}{\sigma (zb_j)},\underset{j=1}{\overset{\kappa }{}}a_j=\underset{j=1}{\overset{\kappa }{}}b_j.$$
(6)
The complex number $`\lambda `$ is related to the overall size of the solitons, the zeros $`a_j`$ and poles $`b_j`$ determine their sizes and positions on $`T_2`$, and the positive integer $`\kappa `$ is the order of the elliptic function $`W`$.

The static energy density (or potential energy density) associated with (6) can be read-off from (4):

$$E=2\frac{|_zW|^2+|_{\overline{z}}W|^2}{(1+|W|^2)^2}.$$
(7)
Pictures of this energy distribution reveal the familiar $`CP^1`$ lumps, as those of figure 1 for example.

The energy is related to the topological charge by the Bogomolnyi bound

$$E2\pi |N|.$$
(8)
The instanton solutions correspond to the equality in (8): solutions carrying $`N>0`$ ($`N<0`$) imply $`_{\overline{z}}W=0`$ ($`_zW=0`$), the Cauchy-Riemann conditions for $`W`$ being an analytic function of $`z`$ ($`\overline{z}`$). Note that the simplest non-trivial elliptic function is of order two, hence there are no single-soliton solutions on $`T_2`$.

We utilise the pure model (4) for $`N=3`$, but for $`N=4`$ the energy involved is so large that the well-known instability of the planar model (recall its conformal invariance means that the lumps can have any width) breaks the numerical procedure fast. A stabilising Skyrme term must be introduced for this case.

## 3 Numerical procedure

In this paper we want to discuss time dependent solutions of the model (4) and its Skyrme version \[see equation (5)\], where the time dependence describes the movement of solitons. As we do not have analytical time-dependent solutions, we resort to numerical simulations. For this we take fields of the form (6) as the initial conditions. Since during the simulations the field $`W`$ may become arbitrarily large, we have preferred to run our simulations in the $`\varphi `$-formulation of the model. One returns to the $`CP^1`$ formulation inverting (3): $`W=(1\varphi _3)/(\varphi _1+i\varphi _2)`$.

Strictly speaking, truly independent solitons can only be obtained in the asymptotic regime of large soliton separation, which really never happens on a compact manifold such as $`T_2`$. However, each factor $`W_j\sigma (za_j)/\sigma (zb_j)`$ in (6) (when $`a_jb_j`$) roughly represents one soliton, providing a setting to studying more or less independent structures. The present work is limited to systems in the topological classes 3 and 4. These systems move, collide and scatter off upon being set into motion by boosting each $`W_j`$ separately. The initial-value problem is then completely specified by giving both $`W(t)`$ and $`_tW(t)`$ at the initial time $`t=0`$.

For a square torus we have the so-called lemniscatic case  where $`\sigma `$ posseses the simple Laurent expansion

$$\sigma (u)=\underset{j=0}{\overset{\mathrm{}}{}}G_ju^{4j+1},G_jG_j(L)\mathrm{}.$$
(9)
For the accuracy of our calculations it is sufficient to compute the series (9) up to $`G_5`$ (our coefficients for $`j>5`$ are negligible):

$$\begin{array}{ccc}G_0\hfill & =\hfill & 1\hfill \\ G_1\hfill & =\hfill & 0.7878030\hfill \\ G_2\hfill & =\hfill & 0.221654845\hfill \\ G_3\hfill & =\hfill & 9.36193\times 10^3\hfill \\ G_4\hfill & =\hfill & 7.20830\times 10^5\hfill \\ G_5\hfill & =\hfill & 2.37710\times 10^5\hfill \end{array}\}.$$

We have employed the fourth-order Runge-Kutta method and approximated the spatial derivatives by finite differences. The Laplacian has been evaluated using the standard nine-point formula and, to further check our results, a 13-point recipe has also been used. Respectively, the Laplacians are:

$$^2=\frac{\left[\begin{array}{cccc}1& 4& 1& \\ 4& 20& 4& \\ 1& 4& 1& \end{array}\right]}{6\times (\delta x)^2},^2=\frac{\left[\begin{array}{ccccc}& & 1& & \\ & 1& 12& 1& \\ 1& 12& 48& 12& 1\\ & 1& 12& 1& \\ & & 1& & \end{array}\right]}{10\times (\delta x)^2}.$$
(10)
The discrete model has been evolved on a $`n_x\times n_y=200\times 200`$ square periodic lattice with spatial and time steps $`\delta x`$=$`\delta y`$=0.02 and $`\delta t`$=0.005, respectively. The vertices of the fundamental period cell we have used for our simulations were at

$$(0,0),(0,L),(L,L),(L,0),L=n_x\times \delta x=4.$$
(11)
Unavoidable round-off errors have gradually shifted the fields away from the constraint $`\stackrel{}{\varphi }.\stackrel{}{\varphi }=1`$. So, like in the planar case , to correct for this we have rescaled $`\stackrel{}{\varphi }\stackrel{}{\varphi }/\sqrt{\stackrel{}{\varphi }.\stackrel{}{\varphi }}`$ every few iterations. Each time, just before the rescaling operation, we have evaluated the quantity $`\mu \stackrel{}{\varphi }.\stackrel{}{\varphi }1`$ at each lattice point. Treating the maximum of the absolute value of $`\mu `$ as a measure of the numerical errors, we have found that max$`|\mu |`$ $``$ 10-8. This magnitude is useful as a guide to determine how reliable a given numerical outcome is. Usage of unsound numerics in the Runge-Kutta evolution shows itself as a rapid growth of max$`|\mu |`$; this also occurs, for instance, when the solitons pinch-off.

The parameter $`\lambda `$ in (6) has been set to $`\lambda =(1,0)`$ all throughout.

## 4 Three solitons

First we have considered states with three solitons. Our initial configuration is given by taking $`\kappa =3`$ in (6), the elliptic function of order 3

$$W=\frac{\sigma (za_1)}{\sigma (zb_1)}\frac{\sigma (za_2)}{\sigma (zb_2)}\frac{\sigma (za_3)}{\sigma (zb_3)},\underset{j=1}{\overset{3}{}}a_j=\underset{j=1}{\overset{3}{}}b_j.$$
(12)
The values of $`a_j,b_j`$ have been selected in such a way that the solitons lie symmetrically around a circle in the network (11). This is easily achieved by fixing $`a_1,b_1`$ to reasonable values and setting

$$\begin{array}{ccc}a_1^{}\hfill & =\hfill & a_1c\hfill \\ b_1^{}\hfill & =\hfill & b_1c\hfill \end{array}\},$$
(13)
where $`c=(2,2)`$ stands for the centre of the period cell. Then

$$\begin{array}{ccc}a_j^{}\hfill & =\hfill & a_1^{}\mathrm{exp}(i\beta _j)\hfill \\ b_j^{}\hfill & =\hfill & b_1^{}\mathrm{exp}(i\beta _j)\hfill \end{array}\};\beta _j=(j1)\frac{2\pi }{3},j=1,2,3.$$
(14)
This symmetrical arrangement gives solitons of the same size, and satisfies the selection rule in (12) for any choice of the complex numbers $`a_1^{},b_1^{}`$:

$`{\displaystyle \underset{\kappa =1}{\overset{3}{}}}a_j^{}`$ $`=`$ $`(0,0)`$ (15)
$`=`$ $`{\displaystyle \underset{\kappa =1}{\overset{3}{}}}b_j^{}.`$
Next we go back to $`a_j`$ and $`b_j`$ through

$$\begin{array}{ccc}a_j\hfill & =\hfill & a_j^{}+c\hfill \\ b_j\hfill & =\hfill & b_j^{}+c\hfill \end{array}$$
(16)
and supply the system with an initial speed $`v_0`$ by boosting

$$\begin{array}{ccc}a_j\hfill & \hfill & a_j+v_jt\hfill \\ b_j\hfill & \hfill & b_j+v_jt\hfill \end{array}\};v_j=v_0\mathrm{exp}(i\beta _j).$$
(17)
It is now possible to evaluate the time derivative of $`W`$ at the initial time. Inserting (17) into (12) we get

$$_tW(t)_{t=0}=\underset{j=1}{\overset{3}{}}v_j[\frac{_z\sigma (za_j)}{\sigma (za_j)}\frac{_z\sigma (zb_j)}{\sigma (zb_j)}]W_j,$$
(18)
where $`W_j\sigma (za_j)/\sigma (zb_j)`$ denotes the $`j`$-th soliton. Our initial-value problem is defined by (12) and (18).

Choosing

$$\begin{array}{ccc}a_1=(3,2)\hfill & & \\ b_1=(1.45,1.95)\hfill & & \\ v_0=0.35\hfill & & \end{array}\},$$
(19)
gives an initial configuration whose energy density is exhibited in figure 1. The solitons are placed on the vertices of an equilateral triangle; the first lump is situated along the line $`y`$=2 at $`\beta _1=0^{}`$, with respect to which solitons 2 and 3 are rotated $`\beta _2=120^{}`$ and $`\beta _3=240^{}`$, respectively. These angles are readily checked from the picture with the help of a protractor.

The results of simulations for this case are depicted in figure 2 for various times. The three solitons are fired towards the centre of the triangle and in so doing they expand: the peak $`E_{max}`$ of the energy density decreases. The solitonic trio collides head-on and coalesces in a ringish structure, then emerging towards the vertices of the dual triangle, that is, the initial line of approach of a given incoming soliton forms an angle of $`\pi /3`$ with the line along which an outgoing, emerging soliton progresses. This $`\pi /3`$ scattering can be best appreciated in figure 3, where both the initial state ($`t`$=0) and the final state ($`t=2`$) are displayed together.

Returning to figure 2 we see that $`E_{max}`$ becomes narrower with time, particularly after the lumps scatter off and start drawing away from each other. At $`t=2`$, for instance, the maximum value of the energy density goes up to $`E_{max}=579.37`$. Soon after $`t>2`$ this peak gets so spiky that the numerics breaks down: the instability of the planar $`O(3)`$ model takes over and leads to singularity formation.

Our result is noteworthy: the initial configuration, although positioned at the vertices of an equilateral triangle in the period cell, does not produce $`D_3`$ symmetry because the torus itself, being homogeneous but not isotropic, has no such symmetry (the fundamental grid has directed sides). One could reasonably expect the lumps to scatter along directions that need not respect $`D_3`$ symmetry. Therefore, the rationale applied to explain $`\pi /N`$ scattering for the model on $`\mathrm{}_2`$ is no longer valid for the model with periodic boundary conditions.

However, we observe that as the solitons are well localised the boundary conditions may not be very important. Note that our numerical results are consistent with this expectation. To test this further we could place our solitons nearer the edges of the grid and see whether we still observe the 600 scattering. We hope to investigate this issue in the future.

## 5 Four solitons

Next we have looked at the $`N=4`$ configurations. The initial field is given by the elliptic function of order 4

$$W=\frac{\sigma (za_1)}{\sigma (zb_1)}\frac{\sigma (za_2)}{\sigma (zb_2)}\frac{\sigma (za_3)}{\sigma (zb_3)}\frac{\sigma (za_4)}{\sigma (zb_4)},\underset{j=1}{\overset{4}{}}a_j=\underset{j=1}{\overset{4}{}}b_j,$$
(20)
where $`a_j,b_j`$ are chosen so that the solitons sit symmetrically at the vertices of a square in the basic cell. A treatment parallel to that of the previous section leads to

$$\begin{array}{ccc}a_j^{}\hfill & =\hfill & a_1^{}\mathrm{exp}(i\phi _j)\hfill \\ b_j^{}\hfill & =\hfill & b_1^{}\mathrm{exp}(i\phi _j)\hfill \end{array}\};\phi _j=(j1)\frac{\pi }{2},j=1,2,3,4,$$
(21)
The condition between the zeros and poles is again verified $`a_1^{},b_1^{}`$ :

$`{\displaystyle \underset{j=1}{\overset{4}{}}}a_j^{}`$ $`=`$ $`(0,0)`$
$`=`$ $`{\displaystyle \underset{j=1}{\overset{4}{}}}b_j^{}.`$
Both the initial velocity and the time dependence are introduced by boosting:

$$\begin{array}{ccc}a_j\hfill & \hfill & a_j+v_jt\hfill \\ b_j\hfill & \hfill & b_j+v_jt\hfill \end{array}\};v_j=v_0\mathrm{exp}(i\phi _j),$$
(22)
where equation (16) should be kept in mind.

As pointed out in section 2, four lumps involve a large energy and while evolving in time they become too spiky before we can learn anything much about the scattering process. We have therefore studied this system in the Skyrme version of the theory, where the solitons are stable and may be examined for as long as required.

Instead of the Lagrangian density (4) we have thus taken:

$$=\frac{|W_t|^22|W_z|^22|W_{\overline{z}}|^2}{(1+|W|^2)^2}8\theta _1\frac{|W_z|^2|W_{\overline{z}}|^2}{(1+|W|^2)^4}(|W_t|^2+|W_z|^2|W_{\overline{z}}|^2).$$

The configuration (6) is no longer an exact solution of the field equation derived from (5), albeit it is a very good approximation to it. We should also stress that the presence of a Skyrme term does not affect the trajectory of the lumps before the collision or their scattering angle.

Taking

$$\begin{array}{ccc}a_1=(2.05,2.55)\hfill & & \\ b_1=(2.95,2.95)\hfill & & \\ \theta _1=10^3\hfill & & \\ v_0=0.35\mathrm{exp}(i\pi /4)\hfill & & \end{array}\},$$
(23)
entails the state of identical ‘baby skyrmions’ illustrated in figure 4. The solitons $`W_2`$, $`W_3`$ and $`W_4`$ are rotated $`90^{},180^{}`$ and $`270^{}`$ with respect to $`W_1`$, which we have conveniently placed in the first quadrant, roughly on the central diagonal joining the grid points (0,0)-(4,4).

The system (20) gets moving via (22) and zeroes in on the middle of the mesh, where the four skyrmions bump head-on into each other. This dynamics makes the solitons scatter off and emerge towards the vertices of the dual square, emerging at $`45^{}`$ with respect to the initial direction of motion, as depicted in diagram 5. Figure 6 superimposes both the initial state ($`t`$=0) and the final state ($`t=3`$), allowing greater clarity in the appreciation of the scattering angle $`\pi /4`$.

Unlike the outcome of the previous section, the dual square scattering on $`T_2`$ is not so surprising. For although the initial data doesn’t have $`D_4`$ symmetry, it does have 4-fold rotational symmetry and zero angular momentum. Note that the boundary conditions (5) break the $`SO(2)`$ rotational symmetry of the plane into a 4-fold rotational symmetry.

## 6 Geodesic approximation

Note that, in analogy to the $`S_2`$ case, when $`v`$ in (17) and (22) is set to zero our initial configuration is a static solution of the equations of motion. If we change $`a_j`$, $`b_j`$ to a new value given by (17) or (22) for a particular value of $`t`$ the new configuration is again a static solution of the equations of motion. Hence as $`v`$ changes the changes (17) and (22) connect configurations which correspond to static solutions of the equations of motion. Thus it is reasonable to expect that a system set off with a small $`v`$ will follow such a change. This expectation goes under the name of geodesic approximation (the system evolves by changing its zero-mode parameters). Its validity is not expected to depend too much on whether the model is defined on $`S_2`$ or $`T_2`$.

So far we have been concerned with solitons of equal size. Let us now look into the more general and interesting situation of energy lumps of different sizes, illustrating the proceedings by studying the case $`N=3`$.

First we set up the initial configuration and evolve it through our standard numerical simulation. Then we choose a set of collective coordinates to reproduce the results of our simulations (trajectory, scattering), offering an explanation in the framework of the geodesic approximation.

### 6.1 Numerical simulation

Our 3-soliton system is still given by a function of the form (12), but with a layout not as symmetrical as before. Instead of (14) we put

$$\begin{array}{ccc}a_j^{}\hfill & =\hfill & a_1^{}\mathrm{exp}(i\alpha _j)\hfill \\ b_j^{}\hfill & =\hfill & b_1^{}\mathrm{exp}(i\alpha _j)\hfill \end{array}\};\alpha _j=(j1)\frac{2\pi }{3}+(1)^{(j1)}(1\delta _{1j})\xi ,$$
(24)
with $`j`$=1,2,3 and the numbers $`a_1^{},b_1^{}`$ being fixed as customary. The initial three soliton configuration can be written as

$`W`$ $`=`$ $`{\displaystyle \frac{\sigma (za_1)}{\sigma (zb)}}{\displaystyle \frac{\sigma (za_2)}{\sigma (zb_2)}}{\displaystyle \frac{\sigma (za_3)}{\sigma (zb_3)}}`$ (25)
$`=`$ $`{\displaystyle \frac{\sigma (z^{}a_1^{})}{\sigma (z^{}b^{})}}{\displaystyle \frac{\sigma (z^{}a_1^{}e^{i\alpha _2})}{\sigma (z^{}b_1^{}e^{i\alpha _2})}}{\displaystyle \frac{\sigma (z^{}a_1^{}e^{i\alpha _2})}{\sigma (z^{}b_1^{}e^{i\alpha _2})}},`$
where

$$b^{}=a_1^{}+2(a_1^{}b_1^{})\mathrm{cos}\alpha _2[\alpha _2=\frac{2\pi }{3}\xi ]$$
(26)
ensures that the zeros and poles in (25) comply with the constraint in (6). As usual we switch between primed and unprimed numbers through formula (16), plus $`b^{}=bc`$ and $`z^{}=zc`$. We have also used the fact that $`\mathrm{exp}(i\alpha _3)=\mathrm{exp}(i\alpha _2)`$.

Note that equating $`\xi `$ to zero simplifies $`\alpha _j`$ to $`(j1)2\pi /3=\beta _j=`$ of (14), whereupon the elliptic function (25) reduces to the field (12) simply because $`\mathrm{cos}\alpha _2=\mathrm{cos}\beta _2=1/2b^{}=b_1^{}`$ according to (26).

An angle $`\xi 0`$ generates solitons of different sizes which no longer enjoy the positional symmetry boasted by the arrangement (14). In the set-up (14), where $`\xi `$=0 and the solitons have the same size, the required relationship $`a_j^{}=b_j^{}`$ looks after itself as shown in (15). For $`\xi 0`$ a demand of the kind (26) is needed.

Now we have evolved the solitons (25) in the pure version of the $`CP^1`$ model, with the solitons being sent into collision in the regular fashion:

$$\begin{array}{ccc}a_j\hfill & \hfill & (a_1^{}v_0t)\mathrm{exp}(i\alpha _j)+c\hfill \\ b\hfill & \hfill & (b^{}v_0t)+c\hfill \\ b_s\hfill & \hfill & (b_1^{}v_0t)\mathrm{exp}(i\alpha _s)+c\hfill \end{array}\};j=1,2,3;s=2,3.$$
(27)
The associated energy distribution is shown for a typical case by the contour plot of figure 7, which corresponds to a choice of parameters (19) in addition to $`\xi =10^{}`$.

The starting configuration is represented by the three wider structures whereas the final, scattered state is the narrower trio plotted therein. The sense of motion is clearly indicated by the arrows. We note that for both configurations the solitons are rotated an extra angle $`\xi `$ with respect to $`W_1`$, as compared to their partners of figure 3, And, unlike the latter, it is not the case that the three initial lumps have the same size; nor they are situated at the vertices of an equilateral triangle. Apparent as well from our simulations is that the scattering angle differs from $`\pi /3`$. This is a consequence of considering nonsymmetrical solitons, whose collisions are not elastic and thus involve energy transfer.

Clearly the distance between the solitons differ for different pairs of them and so their interactions are not the same resulting in a different scattering pattern. Can we explain this difference? In the next subsection we show that an explanation can be provided in terms of, appropriately chosen, collective coordinates and the motion following appropriate geodesics.

### 6.2 Collective coordinates

In order to make a wise selection of collective coordinates let us consider closely the positions and sizes of the solitons defined in the previos subsection.

The location of the solitons (25) are clearly given by

$$\frac{a_1^{}+b^{}}{2},\frac{a_1^{}+b_1^{}}{2}e^{i\alpha _2},\frac{a_1^{}+b_1^{}}{2}e^{i\alpha _2},$$
(28)
and their sizes are

$$|\frac{a_1^{}b^{}}{2}|,|\frac{a_1^{}b_1^{}}{2}|,|\frac{a_1^{}b_1^{}}{2}|.$$
(29)
Put

$$a_1^{}=k\chi ,b_1^{}=k+\chi $$
(30)
so the positions (28) adopt the form

$$k\chi [1+2\mathrm{cos}(\alpha _2)],ke^{i\alpha _2},ke^{i\alpha _2},$$
(31)
while the sizes (29) read

$$|2\chi \mathrm{cos}(\alpha _2)|,|\chi |,|\chi |.$$
(32)
Consistently, for $`\xi =0`$ the description (31)-(32) corresponds to objects of the same size $`|\chi |`$ which are symmetrically situated at $`k`$ and $`k\mathrm{exp}(\pm i2\pi /3)`$ (with respect to the centre). For simplicity, we have taken two lumps with the same size and a third lump with a different size.

A possible collective coordinate description involves treating $`k`$, $`\chi `$ and $`\alpha _2`$ as collective coordinates. Thus, in the simulation we expect $`\alpha _2`$ and $`\chi `$ to remain approximately constant and $`k`$ to vary. We are suggesting that the scattering can be understood as proceeding (on average) with only $`k`$ depending on time and varying from (take $`k`$ real for simplicity) $`k>0`$ to $`k<0`$. It is easy to show, although a bit tedious in practice (this involves estimating various elliptic integrals which can be done partially analytically, partially numerically), that both on $`\mathrm{}_2`$ or on our torus the kinetic energy of the motion just described is finite. So such a behaviour is possible.

We have carried out our collective coordinates “motion” using the specific values

$$\chi =(0.55,0),k=(u,0),u[1,1],$$
(33)
where $`u`$ varies across the interval in steps of 0.2. We see from (31) and (33) that the configurations for $`u>0`$ ($`u<0`$) correspond to incoming (outgoing/scattered) lumps. The value $`u0`$ represents the situation where the solitons have collided and are on top of each other, coalescing in the centre of the mesh. It is largely the behaviour of $`W`$ at $`u0`$ what determines the scattering angle.

The state for $`k`$=(1,0) (incoming lumps) and $`k`$=(-0.2,0) (“scattered” lumps) are depicted in figure 8 for the case $`\xi =10^{}`$. As long as the path followed by the lumps and the scattering angle is concerned, the similarity with the numerically evolved situation of figure 7 is clear. We stress that since we are mostly interested in the relative position of the solitons and their scattering angles, it is immaterial that the breadths of the solitons in figure 7 differ from the solitons of figure 8.

A more detailed comparison can be made by viewing a trajectory plot, one showing consecutive snapshots \[corresponding to $`k`$=(1,0), (0.8,0),…,(-0,2),(-0.4,0)…(-1,0)\] of lump positions as they approach each other and scatter off.

This is shown in figure 9, where the full collective coordinate motion has been ticked with small circles and the motion according to the time evolution of subsection 6.1 has been sketched with solid lines. Note that after scattering the path of the numerically evolved lumps (continuos lines) cannot be followed much farther; this is because the solitons get very spiky and the simulation breaks down.

Our simulations have been carried out for several values of $`\xi `$ and have been compared with the corresponding collective coordinate motion. We have always found the trajectories from both approaches to be in agreement, thus supporting our choice of collective coordinates. Another illustration is provided by figure 10 for the case $`\xi =20^{}`$.

## 7 Conclusions

In this article we have studied head-on collisions between solitons in the (2+1)-D $`CP^1`$ model with periodic boundary conditions, that is, with the model defined on a flat torus $`T_2`$.

Through numerical simulations we have found that solitons of equal sizes (at the initial time) scatter at an angle $`\pi /N`$, or dual-polygon scattering, where $`N`$ is the soliton number or topological charge of the system. In this paper we have focused on $`N=3,4`$ (the case $`N=2`$ has been considered previously  and it has been found to comply with the $`\pi /2`$ scattering).

Unlike the usual model on $`\mathrm{}_2`$, our model with periodic boundary conditions breaks the $`SO(2)`$ rotational invariance of the plane, leaving us with a numerical mesh with directed sides where the initial soliton configuration has no symmetry under the dihedral group $`D_N`$. So it is remarkable to still find dual-triangle scattering in this format. On the other hand, since the $`SO(2)`$ symmetry breaks into a four-fold rotational symmetry, the dual-square scattering is less unexpected.

We have also considered solitons of different sizes (at the initial time) and have observed that the scattering angle in no longer $`\pi /N`$, outcome arising from the fact that there is energy transfer in collisions between unsymmetrical solitons.

By reparametrising the quantities describing the positions of the solitons using a juditious set of collective coordinates, we have been able to reproduce the above numerical results, thus offering an explanation of the scattering process. We have illustrated this approach using a 3-soliton field.

These results raise important questions pertaining to the the interplay between the symmetry of the initial configuration and the lack of symmetry of the torus itself. As pointed out at the end of section 4, the non-isotropy of the torus might affect the evolution of the lumps if they are initially placed near the boundary of the mesh. What about systems with $`N`$=5,6,…? Numerical experiments on a periodic, rectangular grid would also be worth performing. We hope to report on these matters in the near future.

Acknowledgements

Work for this paper was carried out during the authors’ visits to Durham and Maracaibo, supported by the ROYAL SOCIETY, CONICIT and LUZ-CONDES. The paper was completed at IISC, Bangalore-India, where $`\mathrm{}`$J Cova stayed thanks to a grant TWAS-UNESCO. He thanks Prof Uberoi and Prof Rangarajan for their hospitality at IISc.