text:
# Precise Measurements of Atmospheric Muon Fluxes with the BESS Spectrometer

## 1 Introduction

The evidence for atmospheric neutrino oscillation has been reported from the Super-Kamiokande collaboration by using high-statistics samples of muon neutrino events . There are two major sources of systematic errors in evaluating the neutrino flux; the flux of primary cosmic-rays and the production cross sections of secondary mesons; pions and kaons . Recently, the fluxes of primary cosmic-ray particles, mainly consisting of protons and helium nuclei, have been measured precisely by two independent and consistent observations . Although the details of interaction model itself is hard to be determined, the measurement of atmospheric muons plays crucial role in evaluating the flux of atmospheric neutrinos because muons and muon neutrinos are produced always in pairs as decay products of mesons and the kinematics of meson and muon decay is well known.

The muon flux at sea level has been measured by many groups. However, there are large discrepancies among those measurements much larger than the statistical error quoted in each publication. Therefore it is conceivable that the difference comes from systematic effects such as uncertainties in momentum determination, geometrical factor, exposure time, particle identification, trigger efficiency and normalization procedure.

We report here precise measurements of the absolute flux of atmospheric muons at sea level at Tsukuba ($`36.2^{}N`$,$`140.1^{}E`$), Japan and Lynn Lake ($`56.5^{}N`$,$`101.0^{}W`$), Canada by using the BESS spectrometer . The data were collected in ’95 (at Tsukuba) and in ’97, ’98 and ’99 (at Lynn Lake). The cutoff rigidities are 11.4 GV (at Tsukuba) and 0.4 GV (at Lynn Lake).

## 2 Spectrometer Setup

The BESS spectrometer was designed as a high resolution spectrometer with a large geometrical acceptance to perform precise measurements of primary and secondary cosmic-rays as well as a sensitive search for rare exotic particles of primary origin. Cross sectional views of ’95 and ’99 configurations are shown in Fig. 1. The spectrometer configuration was updated in ’97 as described below, and was kept similar in ’98 and 99 except for shower counters installed in ’99.

The thin superconducting coil  (4.70 g/cm2 thick including the cryostat) produces a uniform axial magnetic field of 1 Tesla. A jet-type drift chamber (JET), inner drift chambers (IDCs) and outer drift chambers (ODCs) are located inside and outside the coil. These chambers are operated with a slow gas (CO2 90 %, Ar 10 %). Tracking signal from the drift chambers are read out by flash ADCs. The $`r\varphi `$-tracking is performed by fitting up to 28 hit-points, each with a spatial resolution of 200 $`\mu `$m. Tracking in the $`z`$-coordinate is made by fitting points in IDC measured with vernier pads with an accuracy of 470 $`\mu `$m and points in the JET chamber measured using charge-division with a spatial resolution of 20 mm. By using these data, we performed the continuous and redundant 3-dimensional track information. In order to get momentum of particle, we used 28 hit-points of the JET chamber and IDCs in the magnetic field. The ODCs provide extra hit positions outside the magnet and is used to calibrate the JET chamber and IDCs. In addition, all drift chambers have capabilities to distinguish the multi-hit. This feature enables us to recognize multi-track events, thus we could see the tracks having interactions and scatterings. The time-of-flight (TOF) scintillator hodoscopes measured the velocity of particles with a time resolution of 110 ps in ’95. The acrylic Čerenkov shower counter consists of acrylic and lead plate (12 mm). These counters are placed outside the lower TOF counter. The acrylic Čerenkov shower counters, used to separate electron and muons, were installed only for the ground observation. The total material thickness from outside the pressure vessel, passing through superconducting magnet coil, inside the JET chamber was 9.03 g/cm2.

Since ’97 experiment, we installed a newly developed threshold-type Čerenkov counter with silica-aerogel radiator, after removing the ODCs . The resolution of TOF was improved to 75 ps by using new photomultipliers (PMTs) with a larger diameter for better light collection . In ’99 experiment, we installed a part of the shower counter just below the superconducting magnet.

## 3 Data Samples

The ’95 “ground” experiment was carried out at KEK, Tsukuba ($`36.2^{}N`$,$`140.1^{}E`$), Japan, from December 23 to 28. KEK is located at 30 m above sea level. The vertical cutoff rigidity is 11.4 GV ($`\lambda =26.6^{}N`$ at geomagnetic latitude ). The mean atmospheric pressure in this experiment was 1010 hPa (1030 g/cm2). The scientific data were taken for a live time period of 291,430 sec and 9,148,104 events were recorded on magnetic tapes. The ’97, ’98 and ’99 ground experiments were carried out in Lynn Lake ($`56.5^{}N`$, $`101.0^{}W`$), Canada, on July 22, August 16 and July 26, respectively. The experimental site in Lynn Lake is located at 360 m above sea level. The vertical cutoff rigidity is 0.4 GV ($`\lambda =65.5^{}N`$ at geomagnetic latitude ). The mean atmospheric pressures in Lynn Lake experiments in ’97, ’98 and ’99 were 980.6 hPa (1000 g/cm2), 990.5 hPa (1010 g/cm2) and 964.9 hPa (983.9 g/cm2), respectively. The total scientific data were obtained for a period of 21,304 sec (7,011 sec, 3,949 sec and 10,344 sec) of live time and 242,934, 137,629 and 354,869 events were recorded on the magnetic tapes, respectively.

The trigger was provided by a coincidence between the top and the bottom scintillators of TOF counters. All triggered events were gathered in the magnetic tapes. The core information (momentum, T.O.F., etc.) was composed and extracted from the original data. There were two kinds of efficiencies so as to gather atmospheric cosmic-ray data ; trigger efficiency ($`\epsilon _{trigger}`$), track reconstruction efficiency ($`\epsilon _{reconstruction}`$).

## 4 Data Analysis

At first, the following off-line selections were applied for the recorded events.

(i) One or two counters are hit in each layer of the TOF hodoscope and only one track should be found in the JET chamber.

(ii) Track should be fully contained in the fiducial region, namely the number of hits in the JET chamber expected from the trajectory should be 24 and the extrapolated track should cross the fiducial region of TOF scintillators ($`|z|<43.0`$ cm).

We call an efficiency that pass through these selection by the name of single track efficiency ($`\epsilon _{single}`$). We used Monte Carlo calculation in order to obtain efficiency which depend on the momentum.

Next, we selected muon tracks from tracks that pass through the above selection. In order to select the muon tracks, we used the time-of-flight and rigidity information obtained by the TOF scintillation counters and drift chambers, respectively as shown in Fig. 2. We selected the muon tracks using ”muon $`\beta ^1`$-band cut” which are defined by :

$$\frac{1}{\beta }=\sqrt{\left(\frac{m}{R}\right)^2+1}\pm \mathrm{\hspace{0.33em}3.89}\sigma .(0.01\%)$$

Here, $`\beta `$ is velocity of particle, $`m`$ is muon mass and rigidity($`R`$) is momentum per charge. We selected particles which pass through this requirement ($`<|3.89\sigma |`$), and we call this selection efficiency muon selection efficiency ($`\epsilon _{\mu select}`$). From the plots, protons, electrons and positrons are major sources of background events that contaminates the muon bands.

The rejection of electrons and positrons was performed by utilizing the Acrylic Čerenkov shower counter. Proton events could be eliminated in a rigidity range below 1.4 GV by muon $`\beta ^1`$-band cut. Above this rigidity, we used other experimental data of proton flux  to reduce the contamination of protons into muon $`\beta ^1`$-band. A contamination of protons in the muon $`\beta ^1`$-band was estimated to be 2.0 % at 1.4 GV and decreased rapidly with rigidity. According to the work of R. L. Golden et al. , the proton flux at sea level follows the power spectrum with an index of about $`3.0`$ from 2 to 20 GV, steeper than the index of muon flux on the ground level. The protons were subtracted from observed muon $`\beta ^1`$-band using the result of R. L. Golden et al. normalized to number of protons below 1.4 GV obtained in this experiment. A contamination of electrons and positrons in the muon $`\beta ^1`$-band cut was estimated by using Acrylic Čerenkov shower counter. The electrons interacts with lead plates (12 mm) and it generates shower in acrylic plates. Therefore Čerenkov light yielded by electrons are distinguished from that of muons . A contamination was estimated to be about 2 % at 0.5 GV and decreased drastically with rigidity, because the electron flux had steeper index than that of the the muon index. We subtracted electrons and positrons from muon $`\beta ^1`$-band cut by using the estimation obtained by analysing the Acrylic Čerenkov shower counter. The systematic errors of these subtraction was 1 % for protons at 1.4 GeV/$`c`$ and 0.5 % for electrons and positrons at 0.6 GeV/$`c`$. The systematic errors decreased drastically as rigidity increases and it was negligible at 20 GeV/$`c`$. We used muon events that pass through only these selections to obtain muon energy spectrum. As we shall see later, we had about 98.9 % efficiency to take the muon events.

Based on these muon events, we obtained the muon rigidity spectrum at the top of instrument (TOI) in the following way: The TOI energy of each event was calculated by tracing back the particle through the spectrometer material and correcting energy loss by using GEANT 3.21. The corrections were usually small, about 10 MeV for a 1 GeV event.

Among the factors necessary to obtain the flux, the geometrical acceptance can be calculated reliably by Monte Carlo (M.C.) methods due to the simple geometry and the uniform magnetic field of the BESS spectrometer. The geometrical acceptance for the vertical muons ($`\mathrm{cos}\theta 0.98`$) taken in Tsukuba (’95) was about 0.03 m2sr above 2 GeV/$`c`$ and decreased gradually at lower momentum. Because east and west effect is not important at high latitude, we analyze the data taken at Lynn Lake (’97, ’98 and ’99), in a range of $`\mathrm{cos}\theta 0.90`$ (0.09 m2Sr). The systematic error caused by the east-west effect in Tsukuba was estimated to be 1.0 % by comparison with the experimental data and the isotropic M.C. calculation, and to be negligible in Lynn Lake. The mean value of zenith angle distribution of muon flux was $`\mathrm{cos}\theta =0.990`$ for Tsukuba data ($`\mathrm{cos}\theta 0.98`$), and $`\mathrm{cos}\theta =0.955`$ for Lynn Lake data ($`\mathrm{cos}\theta 0.90`$). We estimated that the total systematic error of the geometrical acceptance was 0.4 %. The geometrical acceptance can be calculated reliably both by an analytical method and by Monte Carlo methods at simple geometries, such as circle, quadrangle, etc. The difference of the results obtained by both methods was negligible (less than 0.2%). A systematic error of geometrical acceptance due to imperfect alignment was dominant. The livetime fraction of exposure time was 94.9 % (’95) and 99.3 % (’97, ’98 and ’99); the error due to this factor was negligibly small.

In summary, the efficiencies used in deriving the muon flux were trigger efficiency ($`\epsilon _{trigger}`$), track reconstruction efficiency ($`\epsilon _{reconstruction}`$), single track efficiency ($`\epsilon _{single}`$) and muon selection efficiency ($`\epsilon _{\mu select}`$). The trigger was provided by a coincidence between the top and the bottom scintillators, with the threshold set at 1/3 of the pulse height from vertically incident minimum ionizing particles. $`\epsilon _{trigger}`$ was obtained from pulse height distribution of the TOF counter. The efficiency for the trigger ($`\epsilon _{trigger}`$) was estimated to be 99.95 %. All triggered events were recorded in magnetic tape, thereafter data summary tape (DST) was constructed by using the calibration data base. The DST contains information of the track (momentum, track length, etc.), therefore only reconstructed events were filled in DST and we analyzed the muon flux by using DST. In order to estimate $`\epsilon _{reconstruction}`$, we made off-line scanning (eye scanning) for about 1000 tracks by using data made before DST and $`\epsilon _{reconstruction}`$ was found to be 99.5 %. The single track selection efficiency ($`\epsilon _{single}`$) was obtained from the M.C. simulation and the systematic error was estimated by examining agreements between observed and simulated distributions of the values used in the single track selection. $`\epsilon _{single}`$ was found to be 99.5 %. The M.C. data agreed with the real data within 1.5 % in total. In order to select the muon events, we utilized the $`\beta `$ band cut that has a width of 3.89 $`\sigma `$, thus the muon selection efficiency ($`\epsilon _{\mu select}`$) was 99.99 %. From the efficiencies mentioned above, the total efficiency was found to be 98.9 %.

In order to eliminate possible influence of the momentum resolution to the muon flux, we used momentum up to 20 GeV/$`c`$. The momentum resolution of BESS spectrometer was $`\mathrm{\Delta }P/P=0.005P`$ (M.D.M.$`=200`$ GeV/$`c`$). Therefore the errors of the muon flux was 1 % at 20 GeV/$`c`$ by M.C. calculation if we assumed the spectral index of muon flux is $`2.7`$. Our previous paper  discussed about this spectrum deformation effect of the BESS spectrometer. As the momentum decreases, this error decreases. Then the errors caused by this effect was negligible at 0.6 GeV/$`c`$.

Summation of all the estimated systematic errors were 2.4 % for positive muons and 2.2 % for negative muons in Tsukuba and 2.2 % for positive muons and 1.9 % for negative muons in Lynn Lake.

## 5 Atmospheric Effect

Variations in cosmic-ray flux by the change of the atmospheric conditions is called ”atmospheric effect”. It has been known that there are two main sources of this effect  due to variations of the atmospheric pressure and temperature. Denoting integral flux of the muons at depth $`x_0`$ (g/cm2) as $`I(\overline{E_0},x_0,\theta )`$, and the changes of atmospheric pressure and temperature as $`\delta P`$ (mb) and $`\delta T(x)`$ at $`x`$ (g/cm2) $`(x<x_0)`$, we have a relation of

$$\delta I(\overline{E_0},x_0,\theta )/I(\overline{E_0},x_0,\theta )=\beta (\overline{E_0},x_0,\theta )\delta P+_0^{x_0}\alpha (x,\overline{E_0},x_0,\theta )\delta T(x)𝑑x,$$

Here, $`E_0`$ is the total energy of muons at $`x_0`$, $`\overline{E_0}`$ is the threshold energy, and $`\beta (\overline{E_0},x_0,\theta )`$ and $`\alpha (x,\overline{E_0},x_0,\theta )`$ are so called “barometric coefficient” and “partial temperature coefficient”, respectively.

In order to get the barometric coefficient, we used two sets of ’95 experimental data taken at different atmospheric pressures with a deviation of 25 hPa. Fig. 3 shows the barometric coefficient for the integral muon flux. The barometic effect has a negative correlation, and then flux decreases if the atmospheric pressure increases. A specific negative correlation due to the increases of the $`\mu e`$ decay is expected dominant below 2 GeV/$`c`$ and another specific negative effect due to the absorption by the ionization loss becomes dominant above 2 GeV/$`c`$. The observed coefficient seemed to be consistent with calculated values as shown in Fig. 3. The effect at the 25 hPa pressure-difference on the muon flux amounts to be 2.5 % below 1 GeV/$`c`$ and less than 1 % above 5 GeV/$`c`$.

The temperature effect was calculated by using a temperature coefficient reported by S. Sagisaka , and observed variations in ’95 experimental data are shown in Fig. 4. We used high altitude temperature data observed by using a radio sonde data taken at Tateno Meteorological Observatory ($`36.1^{}N`$, $`140.1^{}E`$, 10 km south of KEK) . In order to analyse the temperature effect, we used two data sets which were taken at different temperature at the ’95 experiment. The observed variation seemed to be consistent with the calculated variation. The variation of muon flux due to the temperature effect in the period of this experiment was less than 1 %.

## 6 Solar Modulation

Fig. 5 shows annual variation of the muon flux in Lynn Lake. In order to distinguish small difference of each variation, the ’97 and ’98 muon fluxes are divided by the ’99 muon flux, and the flux of each year were combined in a wide momentum region to reduce the statistical errors. It is clearly shown that the ’99 flux is lower than other fluxes. The ’99 experiment was performed at the lowest ambient pressure among other three measurements. Since the barometric effect has the negative correlation, the reduction of the flux in ’99 can not be explained by the barometric effect. We need to takes into account an effect of the solar modulation. The solar activity varies globally with the 11 year solar cycle and the solar minimum was ’96 - ’97 and the solar maximum happened between ’00 and ’01 according to observations of sunspot numbers . However, this effect appears about one year later in neutron monitor data . Not only the muon flux, but also muon charge ratios ($`\mu ^+/\mu ^{}`$) decrease below 3.5 GV if the low energy primary proton flux decreases by the solar modulation. These charge ratios of ’97, ’98 and ’99 experiments in this energy region (0.58 - 3.44 GeV/$`c`$) were $`1.258\pm 0.017`$, $`1.235\pm 0.022`$ and $`1.218\pm 0.014`$. These decreases were consistent to decreases of the muon flux. The decrease of charge ratios and muon fluxes are caused by decreasing of primary proton flux, therefore the effect of solar modulation were observed. The BESS spectrometer observed variation of the primary proton flux due to the solar modulation effect at an altitude of 37 km (launched from Lynn Lake) from ’97 to ’00 . These variation shows about 20 % decrease at 2 GeV/$`c`$ from ’97 to ’99 and the difference of these fluxes becomes much smaller at higher momentum. The mean energy of the primary proton responsible to muons of 0.5 GeV/$`c`$ at sea level is about 20 - 30 GeV, and then the degree of muon flux is much smaller than that of primary proton flux at the same energy.

These difference of muon flux were 3 % around 1 GeV/$`c`$. These differences were within statistic and systematic errors with small bins and it was important to obtain spectral shape with small statistic errors. Therefore ’97, ’98 and ’99 experimental data sets were combined to obtain the muon flux in Lynn Lake.

## 7 Results

Fig. 6 shows the resultant positive and negative muon fluxes, and Table. 1 and Table. 2 summarize those data with systematic and statistic errors. We have observed the vertical fluxes of the positive and negative muons in a momentum range from 0.6 to 20 GeV/$`c`$ with an estimated systematic error of 2.4 % for the positive muons and 2.2 % for the negative muons in Tsukuba, and 2.2 % for the positive muons and 1.9 % for the negative muons in Lynn Lake. The cutoff rigidity at Tsukuba is much higher than Lynn Lake. By comparing the data collected at two different geomagnetic latitudes, we have observed an effect of cutoff rigidity.

Fig. 7 shows the total (positive and negative) differential muon spectra at Tsukuba and Lynn Lake, together with previous measurements . Our data on the muon fluxes those which were multiplied by $`p^2`$ at sea level are shown in Fig. 8. From these figures, it is clearly seen that the muon flux measured in Tsukuba and in Lynn Lake were different in lower momentum ranged below 3.5 GeV/$`c`$, but were in good agreement in higher momentum beyond 3.5 GeV/$`c`$. This is because the cutoff rigidity for primary cosmic-rays does not affect in higher momentum.

Fig. 9 shows ratios of positive and negative muons together with the previous measurements , and the results are summarized in Table. 3 and in Table. 4. It was seen that the charge ratio obtained in Tsukuba decreased below 3.5 GeV/$`c`$ while the charge ratio obtained in Lynn Lake remained almost constant value even in this energy range. This difference comes from the influence of the geomagnetic cutoff rigidity. Because the low momentum muons must be generated by the higher momentum protons at Tsukuba. The muon charge ratio observed in Tsukuba had to include systematic errors of proton subtraction and east-west effect. On the other hand, the muon charge ratio observed in Lynn Lake had systematic errors due to only proton subtraction. The systematic error due to proton contamination was less than 2 %, therefore the muon charge ratio observed in Lynn Lake had very small systematic errors.

## 8 Discussion

The obtained momentum spectrum appeared to be good agreement with recent CAPRICE 94  data using the instruments of magnetic spectrometer. These data agreed well within the systematic and statistic errors. But the results of previous experiments were about 20 % larger than these recent experimental data. Most of these previous experiments needed normalization point in order to determine absolute muon fluxes. Since the spectrum shape is similar enough among those experiments, systematic errors of the absolute fluxes are supposed to be the main cause of these difference of absolute fluxes. Therefore we performed our observations with a great care to evaluate the efficiencies to detect the muon tracks. We then have muon fluxes with much smaller systematic errors.

The atmospheric effect was clearly observed. Our result agreed well with the expectation of the analytical calculation. The barometric coefficient had about $`0.1`$ %/hPa at 1 GeV/$`c`$. A temperature effect had less influence on the muon flux in comparison with the barometric effect. Our observation of the temperature effect can also interpreted quantitatively with the expectation of an analytical calculation. Therefore for a precise calculation of the atmospheric neutrino flux precisely, we could include these effect in an analytical way.

The solar modulation effects to the muon flux at the ground level was clearly observed in our experiments. Not only decreasing of the total flux, but also decreasing of the muon charge ratio has been observed. Decreasing of the muon flux should be due to decreasing of primary proton flux according to a temporal variation of the solar modulation. In the atmospheric neutrino flux calculation, it may be important to consider even small changes of the muon fluxes caused by the solar modulation.

## 9 Conclusion

The vertical absolute fluxes of atmospheric muons have been precisely measured with systematic errors of 2.4 % or smaller. We observed the geomagnetic effect by comparing the muon fluxes observed at Tsukuba, Japan and Lynn Lake, Canada. Muon charge ratios obtained at these two sites also showed the geomagnetic effect. The precise measurement of the muon flux at sea level is very important to understand cosmic-ray interactions inside the atmosphere and to decide fundamental parameters to study atmospheric neutrino oscillation.

BESS experiment has been supported by Grants-in-Aid from Ministry of Education, Culture, Sports,Science and Technology, (MEXT) and Heiwa Nakajima Foundation in Japan and by NASA in the U.S.A. The analysis was performed with the computing facilities at ICEPP, University of Tokyo. The data of the CLIMAX neutron monitor were provided by Space Physical Data System of University of Chicago, supported by National Science Foundation Grant ATM-9912341.


format: text

text:
# A very massive spectroscopic binary in the LH 54 OB association in the LMC

## 1 Introduction

Though we can trust that the present development of stellar interiors modelling reproduces qualitatively the consecutive stages along the evolution of stars, there still remain appreciable quantitative discrepancies between theoretical predictions and empirical evidence, especiallly in that concern to the higher mass stars. That is not surprising since, on one hand, their modelling is intrinsically more difficult due to their more complex evolution, and on the other hand, very massive stars are rare because of its shorter lifetimes, so that there are scarce observational constraints.

During the last years, we have performed a campaign of observations of early Magellanic stars in order to enrich the empirical knowledge of these objects, particularly in low metallicity regions.

One of the new binaries discovered for which its light curve was completed is a $`V13`$ mag star in the LH 54 OB association (Lucke & Hodge 1970). This star ($`\alpha =5^\mathrm{h}26^\mathrm{m}24^\mathrm{s}`$, $`\delta =67\mathrm{°}30\mathrm{}13\mathrm{}`$) is bright enough and has sufficiently big radial velocity excursions to be at scope of the 2.15-m telescope at CASLEO111 Complejo Astronómico El Leoncito, operated under agreement between the Consejo Nacional de Investigaciones Científicas y Técnicas de la República Argentina and the National Universities of La Plata, Córdoba and San Juan. (San Juan, Argentina) for radial velocity measurements.

Hill et al. (1994) identified this object as LH 54-425 and performed $`UBV`$ photometry, obtaining $`V=13.19\pm 0.01`$, $`BV=0.31\pm 0.01`$ and $`UB=1.02\pm 0.01`$. Later on, Oey (1996a) obtained $`V=13.13\pm 0.007`$, $`BV=0.215\pm 0.014`$ and $`UB=1.010\pm 0.015`$, listing this star as L54S-4. CCD spectra were also obtained by Oey (1996b) with the ARGUS multifiber spectrograph at the CTIO 4-m telescope. The spectral type assigned by Oey to this star was O4 III (f). She did not detect double lines at the phase of her observation.

LH 54-425 is the earliest star of the OB association LH 54, located at the East-side of the superbubble DEM 192 (Davies, Elliot & Meaburn 1976) (or N51D, Henize 1956) that is probably related with its growth (Oey & Smedley 1998).

## 2 Observations

LH 54 had been included in one of our selected fields for searching new eclipsing binaries. It was observed during four observing runs between 1998 and 2001.

LH 54-425 showed no eclipses, but exhibited small periodic luminosity changes of the order of $`0.1`$ mag, suggesting an ellipsoidal nature of the variability. For this reason we included it as subject for further spectroscopic study as soon as its light curve was sufficiently complete to compute a reasonably reliable ephemeris.

### 2.1 Photometry

CCD images of the region of LH 54 were obtained at CASLEO between 1998 and 2001, with the same instrumentation described in Ostrov et al. (2000). We derived a $`V`$ light curve from aperture photometry performed with DAOPHOT (Stetson 1987, 1991), and tied all the measurements to a unique instrumental system using several stars as local standards to reduce errors, as described in Ostrov et al. (2000). We performed absolute photometry during a photometric night (Dec 3, 1998), obtaining $`V=13.05\pm 0.015`$, and $`BV=0.121\pm 0.011`$. Table 1 displays our photometric measurements of LH 54-425, corrected to the standard system, together with the internal errors derived from the local standards, the $`FWHM`$ and the airmass for each frame.

### 2.2 Spectroscopy

Spectroscopic observations were obtained by means of a REOSC spectrograph in its single dispersion mode, attached to the 2.15-m telescope at CASLEO. The CCD detector was the same Tektronik used for photometry. We used a 600 mm-1 grating giving a reciprocal dispersion of 1.63 Å/pixel on the range from 3900 to 5500 Å. Our resolution, measured from the Cu-Ne-Ar comparison lamp lines, was 2 pixels.

We obtained spectra near the quadratures in series of three 1200 s exposures. The spectrograms were processed and extracted using IRAF222 IRAF software is distributed by NOAO, operated by AURA for NSF. routines. Each series of three observations was combined in one spectrogram, before performing the concomitant measurements.

## 3 Analysis

### 3.1 Ephemeris

We applied to the photometric data two period search methods: one derived from the Lafler & Kinman’s (Lafler & Kinman 1965) and the one described in Shwarzenberg-Czerny (1997). Because of our short time base and the intrinsically small light variations, we can not derive a precise ephemeris. We found that the most probable period is $`P=2.2475`$ days, with aliases each 0.007 days. We used this period and $`E_0=2452216.597`$ to plan the spectroscopic observations.

### 3.2 Spectroscopic classification

We classified the spectra following the criteria of Walborn & Fitzpatrick (1990). The earliest component shows no lines of He i, corresponding thus an O3 spectral type. The other component shows He i 4471 hardly visible due our poor signal to noise ratio (in fact, it is only detectable in our best spectrograms, obtained during dark nights) and more intense He ii lines which is in agreement with an approximately O5 type. With regard to the luminosity classes, we combined all the spectrograms acquired in each quadrature to increment the signal to noise ratio. In the averaged spectrograms, N iv 4058 emission from the primary is visible, together with Si iv $`\lambda \lambda `$ 4089 or 4116 depending on the phase (this is alternatively blended with the H$`\delta `$ absorption of the secondary). We are not able to assign a luminosity class to the secondary. So, we classify the system as O3 III(f)+O5:. Fig. 1 shows the CCD spectrograms of LH 54-425 obtained during dark nights, corresponding to both quadratures.

### 3.3 Radial velocity measurements

The He ii 4686 line appears clearly double in all our spectrograms. He ii 4200 and He ii 4542 appear double in some spectrograms and severely blended in others. The He ii absorption pairs were measured by means of a double gaussian fit using the “deblend” function of the “splot” IRAF task. Only the easily separable He ii pairs were used for the radial velocity determinations. Table 2 lists the heliocentric radial velocities derived for both binary components, together with that corresponding to the nebular \[O iii\]$`\lambda `$ 5007 emission, used to check the stability of the system.

### 3.4 Geometry system’s exploration

Since the system shows no true eclipses, it is not possible to derive reliable physical parameters just from the light curve and radial velocities. For this reason, we assumed an absolute magnitude and investigate which system configurations were compatible with the observations.

We computed the models with the Wilson-Devinney codes (Wilson 1990, Wilson & Devinney 1971). For the bolometric albedos and gravity darkening coefficients we assumed values of $`A=1.0`$ and $`g=1.0`$, respectively, which are adequate for radiative envelopes (Rucinski 1969, Lucy 1976). For the limb darkening a square root law (Díaz-Cordovés & Giménez, 1992) was used, taking the corresponding coefficients from tables by Díaz-Cordovés et al. (1995).

#### 3.4.1 Detached configurations

We first explored the possible detached configurations, so the mode of operation 2 of the Wilson-Devinney program was chosen. In order to reduce the degrees of freedom of the problem, we fixed the total bolometric magnitude of the system. We considered three cases, with $`M_{bol}=10,10.3`$ and $`10.6`$. For the models with $`M_{bol}=10`$ we assigned temperatures of 52500 K and 44500 K, corresponding to spectral types O3 V and O5 V according with Schmidt-Kaler (1982), while for those models with $`M_{bol}=10.6`$ we adopted the temperature scale of Chlebowski & Garmany (1991) for O3 III and O5 III stars, 46500 K and 42300, respectively. For the models with intermediate bolometric magnitude we used the mean of the two temperature scales. With this procedure, the extreme cases in distance moduli and temperature scales were considered.

For each adopted orbital inclination, temperatures and potential $`\mathrm{\Omega }_1`$ of the primary, the secondary’s potential $`\mathrm{\Omega }_2`$ was chosen according with the assumed total bolometric magnitude of the system.

Fig. 2 displays the results of the experiments with detached configurations. Points of bigger size represent smaller O-C light curve residuals. The lower limit corresponds to the primary being 2.5 magnitudes brighter than the secondary component, which can be rejected according to the relative intensity of the corresponding spectral lines. The upper limit corresponds to semi-detached configurations, with the secondary star filling its Roche-lobe. Models above the solid curve have the secondary brighter than the primary, therefore can be also disregarded.

The upper panel, that corresponds to highest temperatures and lower luminosity, shows that the best fitting models are near the semi-detached configuration, with the secondary filling its Roche-lobe. These models are not realistic from an evolutionary perspective, since they have secondary components of bigger size than the primary. If the less massive star is the more evolved one, we have to conclude that mass inversion has occurred, and the system is now semi-detached.

In the middle panel we considered models with intermediate luminosities and temperatures. The best solutions are aligned at the upper limit of the models considered, indicating a semi-detached configuration. There are also good models with $`i58\mathrm{deg}`$ and $`\mathrm{\Omega }_13`$. It means that the primary star is almost filling its Roche-lobe. However, it is not very plausible since the magnitude difference between both stars is near 2.5 mag, which does not agree with the spectroscopic evidence 333 It does not allows to disregard completely such solutions. For example, Niemelä & Bassino (1994) called “primary” at the Roche-lobe filling component of HV 2543, based on the intensity of their spectral lines, while from the W-D modelling that star is approximately 0.9 $`V`$ magnitudes fainter than its companion..

The residuals from models with lowest temperature scale and highest total luminosity (which allows the biggest sizes of the stars) are presented in the bottom panel. In this case, besides of the semi-detached configurations, there exist plausible detached solutions, with $`i47\mathrm{deg}`$. The lowest value considered for the potential $`\mathrm{\Omega }_1`$ (2.85) implies that the primary component is almost filling its Roche-lobe, although (as mentioned above) in these models the luminosity difference between both stars is excessive.

#### 3.4.2 Semi-detached configurations

Considering that some results of the experiments suggested that the system is semi-detached, with the cooler and less massive star filling its Roche-lobe, we explored these configurations, using therefore the mode of operation 5 in the Wilson-Devinney programs. Figure 3 shows the residuals of these models, corresponding to the high and low temperature scales above mentioned, respectively. As the margins of plausible models are wide, it is necessary to consider limits to the total luminosity to restrict the amount of solutions. The upper and bottom curves drawn on figure 3 correspond to models with $`M_{bol}=10`$ and $`M_{bol}=10.6`$, respectively, while models above the horizontal line have the secondary more luminous than the primary component. Configurations with $`\mathrm{\Omega }_1`$ smaller than those plotted in the graphics are in overcontact, and can be also rejected according to the assumed limits in the total luminosity of the system.

## 4 Results

In spite of the wide range of models that can fit well the observations, we can state with reasonable confidence that the inclination is roughly $`4550\mathrm{°}`$, since bigger values would cause eclipses. This implies that the masses are $`100`$ and $`50M_{\mathrm{}}`$ for the O3 and the O5 component, respectively. Table 3 summarises the parameters that best fit suitable detached and semi-detached solutions. It must be noted that the errors quoted in Table 3 correspond to fixed values of $`i`$, $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, i.e., the errors of the masses reflect only the errors in the radial velocity measurements.

The modelled light and radial velocity curves are displayed in Fig. 4 and Fig. 5, together with the observations and O-C residuals. The models correspond to the best fitting detached configuration, although there are not visible difference with the best fitting semi-detached solution.

In the more modest case, the inclination is $`i58\mathrm{°}`$ and the corresponding masses are $`64`$ and $`31M_{\mathrm{}}`$, although this model have a primary 2.3 magnitudes brighter than the secondary component, that is not supported by the spectroscopic data. (The other extreme case, in that both stars have similar brightness, corresponds to $`i39\mathrm{°}`$ and gives masses of $`160`$ and $`80M_{\mathrm{}}`$). All possible solutions indicate that this system contains the more massive star measured at the present.

In a previous work (Ostrov 2001) we found from the analysis of another early Magellanic binary with O3If and O6V components, masses of 41 and 27 $`M_{\mathrm{}}`$. That was an overcontact system, with a still shorter period, of only 1.4 days. These results suggests that extremely early stars have a wide range of masses, at least when they belong to very close binaries, where mass exchange and mass loss play a fundamental role in the evolution.

Indeed, observations of these systems with higher signal to noise and dispersion than that reachable at CASLEO would provide valuable data to investigate the similitudes and differences between these objects.

## Acknowledgments

I wish especially to thank Dr. Claudia Giordano for help with English, and Dr. Nidia Morrell for helpful comments. The author acknowledge use at CASLEO of the CCD and data acquisition system supported through U.S. National Science Foundation grant AST-90-15827 to R. M. Rich. The focal reducer in use at CASLEO was kindly provided by Dr. M. Shara. This research has made use of the Astronomical Data Center catalogs. I would also like to thank the referee, Dr Otmar Stahl, for suggestions that allowed to clarify some issues in the paper.


text:
# Signatures of polaronic excitations in quasi-one-dimensional LaTiO3.41

## ACKNOWLEDGEMENTS

We thank P. Haas for providing additional data. This work was supported by the BMBF (project No. 13N6918/1), and the Deutsche Forschungsgemeinschaft.


text:
# Solar Seismic Model as a New Constraint on Supersymmetric Dark Matter

## 1 Introduction

The dark matter problem has been an important unresolved problem in astrophysics for several decades. It has motivated much activity in cosmology, in both theoretical and observational areas, leading to significant progress in our understanding of the evolution of the universe. To agree with the measured abundances of helium, deuterium and lithium, the baryonic content of the Universe $`\mathrm{\Omega }_b`$ must satisfy $`\mathrm{\Omega }_b=0.04\pm 0.01`$. The dark matter is dominated by non-baryonic relic particles created in the early stages of the Big Bang, amounting to $`\mathrm{\Omega }_m=0.3\pm 0.1`$.

The best-motivated non-baryonic dark matter candidates are Weakly Interacting Massive Particles (WIMPs), which are stable particles, neutral and weakly interacting with ordinary matter. WIMPs were copiously produced in the early universe through their weak interactions with other forms of matter and radiation. As the universe expanded and cooled, their number density became too low for the annihilation processes to keep up with the Hubble expansion rate. A relic population of WIMPs should exist. Lee & Weinberg (1977) showed that if such a stable particle exists, its relic abundance is $`\mathrm{\Omega }_xh^23\times 10^{27}cm^3s^1/\sigma _av`$, where $`\sigma _av`$ is the thermally-averaged product of annihilation cross-section and relative velocity. A particle which interacts with baryonic matter within the range of the weak interactions will lead to $`\mathrm{\Omega }_x`$ of the order of unity, in agreement with cosmological determinations. Specifically, a recent analysis (Melchiorri and Silk 2002) points to $`\mathrm{\Omega }_xh^2=0.12\pm 0.04`$, thereby allowing us to restrict ourselves to study the case of annihilating WIMPs with annihilation cross-section, $`\sigma _av10^{27}cm^3s^1`$.

In this paper we choose as our favorite dark matter candidate the so-called neutralino, arising in supersymmetric (SUSY) extensions of the Standard Model of electroweak interactions. More precisely, we work in the framework of the Minimal Supersymmetric Standard Model (MSSM) as implemented by Gondolo et al. (2001). The neutralino is expected to be the Lightest Supersymmetric Particle (LSP) in the MSSM (Ellis 2001), almost independently of the further simplifying assumptions one has to make to reduce the huge number of free parameters (63) contained in the MSSM.

Neutralinos or other WIMP candidates have not as yet been detected in accelerators, although LEP measurements set a lower bound of about 50 GeV on the WIMP mass (Baltz and Gondolo 2001). If WIMPs populate the halo of the Milky Way, then they can be detected either directly in low-background laboratory detectors or indirectly via observation of anomalous cosmic-ray antiprotons, positrons, and gamma rays from WIMPs that have annihilated in the Galactic halo (Silk & Srednicki 1984). Furthermore, WIMPs may be captured in the Sun (Press and Spergel 1985; Silk, Olive and Srednicki 1985) or in the Earth (Freese 1986; Gaisser, Steigman and S. Tilav 1986; Krauss, M. Srednicki and F. Wilczek 1986; Gould, Frieman and Freese 1989) and annihilate, thereby producing high-energy neutrinos. Annihilations balance capture, leading to an an equilibrium concentration in the solar core ( Gilliland, Faulkner, Press and Spergel 1986; Krauss, Freese , Spergel, Press 1985; Jungman, Kamionkowski and Griest 1996; Lopes, Silk and Hansen 2002). We henceforth only consider solar captures since the resultant neutrinos dominate the terrestrial muon detector signal.

The accretion by the Sun of WIMPs from the galactic background population leads to efficient mechanisms of transport of energy in the solar interior: such WIMPs provide an additional process for transferring energy of the deeper layers to the external layers, depending upon the scattering cross-section of the WIMPs on the solar nuclei. By virtue of their long mean free paths, WIMPs transport energy radially, and tend to produce an isothermal core. WIMP-accreting solar models that have a current relatively low concentration number of WIMPs reduce the temperature at the star’s core and produce a reduction in the solar neutrinos generated by some nuclear reactions of the pp chain and CNO cycle.

This leads to a reduction of solar neutrino counting rates measured on Earth. If such WIMPs exist, a decrease of the solar neutrino fluxes arising from the existence of a WIMP isothermal core could be successfully measured in the coming years by future solar neutrino experiments (Lopes & Silk 2002). Furthermore, if heat transport by WIMPs is significantly changing the structure of the solar core, a particular signature of this peculiar Sun’s core should be expected in solar seismic data (Lopes, Silk and Hansen 2002).

The seismic diagnostics of the Sun’s interior have been for many years the most important constraint on the internal thermodynamic structure of the Sun. Indeed, such research has led to significant improvements in the microphysics such as an update of the equation of state and the opacity calculations, and to a better determination of specific cross-sections of the $`pp`$ chain (Turck-Chièze & Lopes 1993; Christensen-Dalsgaard et al. 1996; Turck-Chièze, Nghiem, Couvidat & Turcotte 2001; Bahcall, Pinssonneault, Basu 2001; Provost, Berthomieu & Morel 2000). This research, motivated by helioseismology, has led to a significant improvement of the so-called Standard Solar Model. In addition, through helioseismic inversion diagnostics, it is possible to build solar models that evolve in an halo of WIMPs, which are consistent with oscillation data and hence predict the annihilating neutrino flux. Furthermore, there exist WIMP-accreting solar evolution models with different masses, scattering cross-sections, and annihilation cross-sections, for which the present solar structure is inconsistent with solar seismic data. The present article discusses the validity of such models and the independent method of diagnostic that seismology can provide for probing the SUSY parameter space.

## 2 The evolution model of the Sun

Application of Newton’s laws to our Galaxy tell us that the luminous disk and bulge must be immersed in a dark halo with a local density of 0.3 GeV cm<sup>-3</sup> and that dark matter particles move with velocities comparable to the local circular speed. Additional theoretical arguments suggest that the velocity distribution of these particles is locally nearly isotropic and nearly a Maxwell-Boltzmann distribution. When a WIMP enters a star, it may interact with nuclei and lose enough kinetic energy to be trapped by the gravitational potential well. The evolution of the Sun is performed within a sea of WIMPs. The WIMP gas tends towards thermalization with baryonic matter with a time-scale much shorter than the time scale of stellar evolution. The WIMP spatial distribution is then simply the barometric equilibrium density, i.e., near a Gaussian with a typical length scale of $`r_x0.13\sqrt{\left(1GeV/m_x\right)}R_{}`$. The more massive are the WIMPs, the more they are concentrated in a very small region within the core of the star. The trapped WIMPs supply another means of transferring energy from the energy producing core to the outside, and thus supplement the usual transfer by photon diffusion (Lopes, Silk and Hansen 2002). WIMPs might have spin-independent (scalar) interactions in which case they would interact with all chemical elements in the Sun, or they might have only spin-dependent (axial) interactions in which case they would interact essentially only with hydrogen.

In fact, at present, direct laboratory searches are restricted to SUSY particles with scalar interactions. For this reason, the research reported in this work will focus only on this type of interactions. Current detectors (DAMA, CDMS, UKDMC) are sensitive to scattering cross-sections $`\sigma _s\genfrac{}{}{0pt}{}{_>}{^{}}10^{42}\mathrm{cm}^2`$, with one to two orders of magnitude improvement possible in the near future. Future detectors (e.g. GENIUS, Cryoarray) plan to be sensitive down to $`\sigma _s\genfrac{}{}{0pt}{}{_>}{^{}}10^{45}\mathrm{cm}^2`$. The WIMP-accreting solar models are computed in a way similar to the standard solar model, the only difference being the existence of an alternative mechanism of transport of energy supplemented by the presence of WIMPs.

As usual, in stellar evolution through the main sequence, we assume that the star is in hydrostatic equilibrium and is spherically symmetric, and that effects of rotation and magnetic fields are negligible. The evolution of the star starts on the pre-main sequence, 0.05 Gyr from the ZAMS. The solar structure and evolution are calculated starting from an initially homogenous star with a given composition. A Henyey method is used to solve the system of nonlinear differential equations describing the stellar structure (Morel 1997). Starting with a standard primordial chemical composition the present solar luminosity and radius is reached at its present age 4.6 Gyr, by readjusting the initial helium abundance and the mixing length parameter (Lopes, Silk and Hansen 2002). It follows that the final WIMP-accreting solar model is very similar to the solar standard model, the difference between the two models possibly being identifiable via seismic diagnostics. The transport of energy by WIMPs is strongly dependent on the mass and scattering cross-sections.

In summary, the heat transport is optimized for $`\sigma _s\sigma _c`$ when the WIMP scale height is roughly equal to its mean free path. $`\sigma _c`$ is a natural geometrical scattering cross-section, depending on the proton mass $`m_p`$ and the radius and mass of the star, $`\sigma _c=m_p/MR^2=8\times 10^{36}\mathrm{cm}^2`$. In order to be effective in heat transport, the WIMPs must have mean scattering cross-section per baryon in the range of $`10^{43}\mathrm{cm}^2\genfrac{}{}{0pt}{}{_<}{^{}}\sigma _\mathrm{s}\genfrac{}{}{0pt}{}{_<}{^{}}10^{33}\mathrm{cm}^2`$, depending upon the annihilation cross-section and mass of the WIMP. The transport of energy by WIMPs falls rapidly outside this range, and it becomes very difficult to test this effect on the solar structure against the solar seismic data. At higher cross-sections the energy is transported locally and the conductivity falls as $`\sigma _c/\sigma _s`$. At lower cross-sections the conductivity falls as $`\sigma _s/\sigma _c`$ and in addition only a fraction are captured by the Sun (Lopes, Silk and Hansen 2002).

## 3 Seismic diagnostic of the solar interior

The WIMPs are thermalized within the solar core and are on Keplerian orbits around the solar center, interacting through elastic scattering with solar nuclei, such as hydrogen and helium, thereby providing an alternative mechanism of energy transport other than radiation. The result is a nearly flat temperature distribution, leading to an isothermal core. Consequently, the central temperature is reduced. This reduction of temperature has two main consequences: since central pressure support must be maintained, the central density is increased in the WIMP-accreting models, and since less hydrogen is burnt at the centre of the Sun, the central helium abundance and the central molecular weight are smaller than in standard solar models. The increase of the central density and hydrogen partially offset the effect of lowering the central temperature in the central production of energy. In fact, this is the reason why minor changes are required to the initial helium abundance and the mixing-length parameter in order to produce a solar model of the Sun with the observed luminosity and solar radius (Lopes, Silk and Hansen 2002). This readily leads to a balance between the temperature, $`T`$, and the molecular weight, $`\mu `$, in the core, leading to the peculiar profile of the square of the sound speed, $`c_s^2T/\mu `$, and the density, $`\rho \mu /T`$, within the solar interior. This seems to be the case for most of the WIMP-accreting solar models.

Figure 1 illustrates the differences between the radial profile of the square of the sound speed and the radial profile of density for WIMP-accreting solar models and the Sun, as predicted by the solar standard model theory. In the same figure, we illustrate the square of the sound speed and the density as inferred for the present Sun by using the data of Global Oscillations at Low Frequency (GOLF) and Michelson Doppler Imager (MDI) experiments. The average mean hydrostatic structure of the Sun has been obtained by an optimally localized averaging inversion method by Kosovichev (1999). The method for independently determining the radial dependencies of the sound speed and density also yields the radial dependence of the first adiabatic index or the chemical composition. All the different methods of inversion are extremely sensitive to the quality of the frequency measurements, and very accurate seismic data for the low-degree acoustic modes is necessary for a precise inversion. The nuclear region is probed by as many as 120 acoustic modes that are significantly influenced by the turbulence, non-adiabatic effects and magnetic field perturbations at the surface layers (Lopes & Gough 2001). Nevertheless, the long duration of continuous measurements has reduced the uncertainty related with dynamics of the outer layers.

In particular the global modes, dipole modes and quadruple modes measured by the GOLF experiment (Bertello et al. 2000), coupled with high-degree modes obtained by MDI experiment data (Rhodes et al. 1997), enable us to significantly constrain the central region of the Sun, where the presence of WIMPs could be detected. However, the inversion of the sound speed seems to be in better agreement with the SSM than the inversion of density. This is due to the fact that this inversion uses only acoustic modes. Naturally, in such a case it is adequate to infer the sound speed profile rather than the density profile in the nuclear region. To successfully obtain the same level of accuracy for the density profile as in the case of the square of the sound speed profile, it is necessary to use gravity modes, which have not yet been unequivocally detected (Turk-Chièze et al. 2002).

Obviously, another reason for the density difference could come from the physics of the solar standard model, and not only from the density profile inversion. However, this is more unlikely because the sound speed difference is very small, reenforcing the view that the physics in the solar nuclear region is already described with the necessary accuracy.

Even if the answers to some questions about the inversion are still unclear, it seems very likely that there exists a class of annihilating WIMPs with annihilation cross-section $`\sigma _av`$ of the order of $`10^{27}\mathrm{cm}^3/\mathrm{sec}`$ and relatively small masses that are excluded by the present seismic results. WIMPs with masses smaller than 60 GeV and scattering cross-sections between $`10^{38}\mathrm{cm}^2`$ and $`10^{40}\mathrm{cm}^2`$ seem to significantly modify the profile of the sound speed near the core.

This result is also reenforced by the density profile inversion. If we accept these results, both inverted quantities reject the existence of WIMPs with masses smaller than $`30\mathrm{GeV}`$ in the proposed scattering cross-section range. It follows that the acoustic spectrum of the previous WIMP-accreting solar models is incompatible with the observed solar spectrum measured by the SOHO seismic experiments. Usually, the solar model of the present Sun which best reproduces the observed acoustic spectrum is referred as the solar seismic model. In particular the solar standard model leads to the best representation of the solar seismic model but it is not the only solution.

## 4 The neutrino flux of annihilating WIMPs

Model-independent predictions can be made for neutrinos from the centre of the Sun, where neutralinos may have been gravitationally trapped and therefore their density enhanced. Today, the rate of change of the number of neutralinos between capture and annihilation in the Sun is in equilibrium. As they annihilate, many of the possible final states produce, after decays and hadronization, energetic neutrinos which propagate out from the interior of the Sun. In particular, the muon neutrinos are useful for indirect detection of neutralino annihilation processes, since muons have quite a long range in a suitable detector medium like ice or water. They can be detected through their Cherenkov radiation after having been produced at or near the detector. Detection of neutralino annihilation into neutrinos is one of the most promising indirect detection methods, and it will be subject to extensive experimental investigations in view of the new neutrino telescopes such as AMANDA, ICeCube, Baikal, BAKSAN, MACRO, NESTOR and ANTARES planned or under construction (Halzen 1997). The neutrino-induced muon flux may be detected in a neutrino telescope by measuring the muons that come from the direction of the centre of the Sun or Earth. The energy of these muons will typically be between 1/2 and 1/3 of the neutralino mass, so they will be much more energetic than ordinary solar neutrinos. These neutrinos have energies of the order of a GeV, well above the energy of solar neutrinos which is of the order of MeV.

To investigate this question, we concentrate our attention on supersymmetric dark matter and evaluate the expected flux of neutrinos (and consequent muon fluxes in detectors) from annihilating neutralinos in the solar core. The numerical results discussed in this work are obtained in the framework of the Minimal Supersymmetric Standard Model as implemented in the DarkSUSY code (DMSSM; Gondolo et al. 2000), which takes into account the most recent particle physics constraints, such as the LEP lower bounds on the lightest Higgs and chargino masses. We extended our analysis to some benchmark points of a different supersymmetric scenario, the Constrained Minimal Supersummetric Standard Model (CSSM; Ellis et al. 2000). In this work we are mainly concerned with WIMP candidates, such as the neutralinos capable of producing changes in the structure of the solar core which can be tested against the solar seismic model. Indeed, even after all the particle accelerator and relic abundance constraints are taken into account, there are large numbers of SUSY models which can produce neutralino dark matter. It is in this large parameter space of candidates proposed by the different extensions of the standard model of particle physics that the Sun can provide another independent diagnostic of the SUSY parameter space.

One should be aware that the choice of nuclear form factors needed to compute neutralino-nucleus elastic scattering, as well as other specific quantities related to hadronic physics which relate quarks/gluons with nucleons and also the quantities related with the step from nucleons to nuclei, are at best approximate. This is the reason why detailed processes such as scattering cross-sections are difficult to obtain with any generality. A more sophisticated treatment would, however, change the values by much less than the spread due to the unknown super-symmetric parameters. Indeed, our understanding of SUSY models is still developing, so predictions of annihilation rates in the early Universe, and thus relic neutralino densities may require modification. More significantly, another source of incertitude in modeling the incoming flux of annihilating neutrinos is related with the poor description of the solar core usually assumed in the computations. However, the present and future capability of solar seismology will provide the means to reduce this source of error.

Future solar seismic experiments will be able to detect deviations of order $`10^5`$ from the luminosity predicted in the solar standard model. If the microphysics of the solar standard model is understood with the necessary precision, then at this level of accuracy a large portion of the supersymmetric parameter space would be ruled out by means of the seismic diagnostics. We show in Fig. 2 how this analysis would affect the expected muon flux, induced by high energy neutrinos from neutralino annihilation in the solar core, in the two different supersymmetic scenarios discussed above: DMSSM and CMSSM. A wide portion of the models obtained in DMSSM, namely the models leading to the highest fluxes, would be probed at $`10^5`$ accuracy using solar seismic data. The benchmark points for CMSSM lie outside most of the parameter space where solar seismic data are sensitive enough to the neutralino annihilation luminosity. The SUSY parameter space that is constrained by helioseismology will be discussed in more detail in a future paper (Bertone, Lopes, Sigl & Silk 2002, in preparation).

## 5 Discussion and Conclusions

Although important questions still remain that need to be addressed regarding the structure of the Sun, such as the asymmetric macroscopic motions in the core, the dynamical effects in the nuclear reaction rates and the chemical abundances in the nuclear region (Turck-Chéze et al. 2001), seismic analysis has been a powerful diagnostic tool for the Sun’s interior. This has led us to constrain the WIMP parameters based on the proposition that the present standard solar model is an accurate approximation of the observed Sun.

The WIMP-accreting solar models of this article are not certainly an unequivocal constraint on the SUSY parameters. If we had modified some of the solar parameters within the error bars, such as for example the age of the Sun, we would probably have found other solutions. However, the general results, as far as the WIMPs are concerned, would be qualitatively the same. The results discussed in the article highlight the powerful tool that solar seismology represents in order to constrain the WIMP parameters within the expected experimental values. Such results strongly favor the non-existence of annihilating WIMPs with $`\sigma _av10^{27}\mathrm{cm}^3/\mathrm{sec}`$, masses lower than 60 GeV, and with scalar cross-sections from $`10^{40}\mathrm{cm}^2`$ to $`10^{38}\mathrm{cm}^2`$.

It should be noticed that the resolution in the inner core is relatively poor; about $`0.05R_{}`$, which is insufficient to detect small discontinuities, probably due to the WIMP isothermal core. This certainly justifies further searching for gravity modes. The modelling of WIMP-accreting solar models is central to the prediction of the muon flux that is expected to be measured by the neutrino telescopes such as AMANDA, ICeCube, Baikal, BAKSAN, MACRO, NESTOR and ANTARES (Halzen 1997 and references therein). Assuming that the WIMP is a neutralino and the density of neutralinos is in agreement with CMB and large-scale structure observations, the prediction of the muon flux is within the range of $`1`$ to $`10^4`$ $`\mathrm{km}^2\mathrm{yr}^1`$, for neutralinos with masses between 60 GeV and 400 GeV. Theoretically, it is the low-degree gravity modes that are the most sensitive to the conditions of the nuclear region, at present the only region where substantial deviations from the standard models occur. Indeed, accreting-WIMP solar models have gravity modes period spacing that are markedly different from that of other solar models, at least for the case of smaller WIMP masses. The asymptotic period spacing of high-order gravity modes of low degree is proportional to an integral of the ratio of the buoyancy frequency by the radius, and its present value for the solar standard model is of the order of $`36min`$ (Lopes & Turk-Chièze 2002, in preparation). Following the period separation between modes of the same degree and consecutive radial-order, the period separation is expected to be relatively smaller than in the case of solar standard models, by as much as a few percent for the WIMP-accreting solar models discussed in this work. It is the sensitivity to the buoyancy that can be used to invert the density in the very deep layers of the solar core. Therefore, gravity mode observations hold the promise of a sensitivity test, although their current interpretation is difficult between the possible candidates detected (Turck-Chièze et al. 2002). Furthermore, the gravity mode data reduction is done assuming that the Sun is well represented by the solar standard model, which could be a major difficulty in terms of detection, if dark matter changes the structure of the solar core.

In any case, the physical processes presented here are meant to be indicative of what one might expect for realistic WIMP-accreting models. Furthermore, since the physics of solar models, the seismic data reduction, as well as the SUSY model predictions are themselves evolving, the detailed model results quoted here should be taken as indicative of the general order of magnitude of expectations for muon fluxes.

## Acknowledgments

The authors wish to thank P. Morel for using the CESAM code, and P. Gondolo, J. Edsjo, L. Bergstrom, P. Ullio, E. A. Baltz, for using the darkSUSY code. Thanks also to S. Hansen, G. Sigl and S. Turck-Chièze for stimulating discussions on the physics of super-symmetric models of particle physics and solar models. IPL is grateful for support by a grant from Fundação para a Ciência e Tecnologia.


text:
# 1 Preference lists are constructed with a step-like structure. In the first positions there is a fraction 𝑓(𝑛₊) of individuals with non-preferred traits. The parameter n controls how efficiently individuals with a biased arrangement of alleles (𝑛₊ away from 1/2) discriminate mates with preferred traits. The function 𝑓_𝜈 is shown for 𝜈={0.1,0.3(∙),1(◆)} and 4. Dots and diamonds show the discrete values of 𝑓_𝜈 for the results of figure 3 (preferences with 30 alleles).

Mating Instabilities Lead to Sympatric Speciation

Catarina R. Almeida\* & Fernão Vistulo de Abreu

\*Depart. de Biologia, Universidade de Aveiro, 3800 Aveiro, Portugal

Depart. de Física, Universidade de Aveiro, 3800 Aveiro, Portugal

First Version: February 2002; Revised: May 2002

Corresponding author: Fernão Vistulo de Abreu abreu@fis.ua.pt

ABSTRACT

One of the most challenging issues of evolutionary biology concerns speciation, the emergence of new species from an initial one. The huge amount of species found in nature demands a simple and robust mechanism. Yet, no consensus has been reached concerning a reasonable disruptive selection mechanism that prevents mixing genes among the emerging species, especially when they live in sympatry. Usually it is assumed that females select males according to their displaying traits, but males perform no selection on female traits. However, recent experimental evidence accumulates towards the existence of male choice. Here we propose a robust mechanism for sympatric speciation, based on the assumption that sexual selection operates in two directions: selection of males by females and of females by males. Complex mating instabilities emerge, creating differential fitness depending on the individuals displaying traits and preferences. When a secondary sexual trait is introduced in a population, due to mutations, the activation of previously neutral genes or due to a different perception of already existent displaying traits, sympatric speciation may result (together with a species recognition system) from a competitive exclusion principle. We suggest that potential candidates to test our theory could be yeasts.

Key index words: sympatric speciation, sexual competition, mating systems

TEXT

Traditionally, two main classes of models try to explain speciation. Allopatric speciation models are the most consensual and assume that the initial population is suddenly divided in two geographically isolated subpopulations, that then diverge until they become reproductively isolated, even after a secondary contact (Mayr 1963). There is however, mounting evidence (Albertson et al. 1999; Panhuis et al. 2001; Schliewn et al. 1994) that speciation might have emerged in sympatry (without geographical isolation). Examples like the cichlids in lake Victoria or many migratory birds, do not seem to fit the basic requirement of long periods of geographical isolation needed for allopatric speciation. Laboratory experiments have also shown that in principle, sympatric speciation is possible (Rice & Hostert 1993). Understanding how sympatric speciation can be driven, has thus attracted much theoretical effort. However, some early works showed that finding a biologically reasonable and robust model seems not to be an easy subject (Felsenstein 1981; Maynard Smith 1966).

Recently, several models have been proposed to explain sympatric speciation. Nevertheless, no consensus has been reached concerning the main driving mechanism, nor the conditions required in practice. One class of models concentrates on sexual competition for mates. So far, these models have been the least convincing models (Turelli et al. 2001) as they have been incapable of providing a robust mechanism, valid under reasonable initial conditions (Takimoto et al. 2000; Turelli et al. 2001). The other class of models assumes that competition for resources is essential to create disruptive selection (Doebeli & Dieckmann 2000; Kondrashov & Kondrashov 1999; see a discussion in van Doorn & Weissing). However it also requires a non-trivial association between traits relevant for ecological adaptation and traits used for sexual discrimination. Further, and in particular in the Doebeli and Dieckmann theory, phenotypes using scarce resources may substitute phenotypes using more abundant resources, as a result of a strong competition for resources in the later case. This is a somehow counter-intuitive result. Again it remains uncertain which species should be more speciation prone, in practice.

Here we propose an alternative approach that leads to a robust theory for sympatric speciation. We give a special emphasis to pair formation. Pair formation is an essential issue in any sexual population, as it is at the basis of the species reproductive success. As we will show, modelling pair formation is essential to create linkage disequilibrium due to the emergence of complex mating instabilities.

In our model we assume that any individual possesses a preference list, where all individuals of the opposite sex are ranked, in descending priority order. The goal of any individual is to optimise his satisfaction by mating with convenient partners. How the preference list is established depends on a response to a sensory stimulus induced by varied secondary sexual traits displayed by the opposite sex individuals. In principle, it could also have an ecological component, to model, for instance, the preference for certain habitats (Bush 1969).

As the number of displaying traits can be quite large and perceived differently by different individuals, a simple assumption will be to consider that each individual has his own preference list and that this is approximately random. In this way we assume that sexual selection has already eliminated non-favoured traits, so that only neutral traits, concerning sexual selection, remain. Note that sexual selection does not necessarily mean the selection of traits relevant for adaptation. Sexual selection may only improve the fitness of a population, by optimising the consequences of the mating conflicts created.

Sexual conflicts lead to intricate dynamics, as both males and females pursue selfish goals, with conflicting wishes. Consider the following preference lists, for a population with three males and three females.

|  | M1 | M2 | M3 |
| --- | --- | --- | --- |
| 1 | F1 | F2 | F2 |
| 2 | F2 | F3 | F1 |
| 3 | F3 | F1 | F3 |

|  | F1 | F2 | F3 |
| --- | --- | --- | --- |
| 1 | M2 | M1 | M2 |
| 2 | M1 | M3 | M3 |
| 3 | M3 | M2 | M1 |

If all individuals search continuously for better partners, then we can imagine the following algorithm. First, males propose sequentially to females ranking higher in their preference lists, and females accept the new partners only if they improve their satisfaction. Then females and males switch roles. In this case, both males and females adopt active strategies towards mating (they both search for better mates). We could sketch that M1 mates F1, M2 mates F2, but then F2 divorces M2 and mates with M3 because he is higher in her preference list, and so on. It is easy to conclude that no stable arrangement is reached. That is, the system has no Nash equilibrium (Gale & Shapley 1962; Omero et al. 1997). While these instabilities are not the rule (for instance, a couple can be stable if they prefer each other), they produce differential costs for reproduction, as fitness depends on the stability of the couples. There is an algorithm, due to Gale and Shapley (Gale & Shapley 1962; Omero et al. 1997), that leads to stable arrangements. This happens if males propose to females, and females dispose, accepting to mate new pretenders only if they improve their satisfaction. Nevertheless, even if a stable solution exists, the time needed to reach it, tGS, tends to be quite large, due to the large number of degrees of freedom and the complex dynamics. In practice, the results we will present do not depend on the strategy adopted by females towards mating.

Methods

We consider discrete generations of a population with an equal number, $`N_{tot}`$, of males and females. The genetic information on each individual controls several characters. A sexual locus may control the strategy towards mating if one would want to study the selection of these strategies, but this, we checked, does not change the results concerning speciation. Henceforth males and females are assumed to take active and passive strategies towards mating, respectively.

One locus with two alleles (we assume dominance) accounts for a secondary sexual trait (A or B). Several independent loci with two alleles ($`+`$ or $``$) define a quantitative preference for one trait. Each individual ranks in a list of preferences individuals of the opposite sex. Preference lists are defined with a simple step-like structure. Calling $`n_+`$ the fraction of positive alleles, then a positively (negatively) biased arrangement corresponds to a preference for mates with trait A (B), such that there is a fraction $`f_\nu (n_+)=[4n_+(1n_+)]^\nu /2`$ of non-preferred mates ranked within the first $`N_A`$ ($`N_B`$) positions ($`N_A`$ and $`N_B`$ being the number of individuals respectively with phenotype A and B). The positive exponent n controls how efficiently individuals with a biased arrangement of alleles discriminate mates with preferred traits. For small values of $`\nu `$, deviations from an unbiased arrangement produces poor discrimination, as can be seen in figure 1.

On each iteration, randomly picked individuals (with an active strategy towards mating) are given one opportunity to find a better partner. New matings happen only if they improve the satisfaction of both elements of the couple. The process repeats until $`2N_{tot}`$ offspring are born. One offspring of each sex is born if a couple stays together during a courtship time tc, plus a reproduction time tr (results concerning the emergence of speciation do not change if offsprings are only allowed to born after a stabilization time). During reproduction time, females cannot engage in new matings, but males may find a better partner (in relation to speciation, the relevant parameter is simply $`t_c+t_r`$). Offspring genotypes are found using Mendel genetic rules.

By introducing a probability m of changing each allele controlling the trait of the offspring, we can model phenomenologically a many genes dependence on the trait or the effect of mutations. A similar procedure can be applied to the loci controlling the preference.

Finally, ornamentation costs can be introduced within our framework by assigning a birth probability, related to the parent’s trait. We checked that speciation is not prevented if the ornamentation costs are not too different, and thus this point will not be discussed further here.

Results and Discussion

Consider first a population with only one trait. In figure 2, the distribution function of mating lifetimes is shown. It exhibits clear power law decay. This scale invariant behaviour means that there is a large amount of stable matings, but also a large amount of very unstable short-lived matings, relatively to what one would expect if matings were random. Thus pair formation will certainly play a role in a sexual selection theory, as it distinguishes individual fitnesses over several scales. In what concerns speciation, another important consequence is that couples stability is context dependent: the same couple may have very different lifetimes depending on the preferences present in the whole population.

Note that scale invariant laws are often associated to self-organized criticality (Jensen 1998). Here, however, it results from the existence of preference lists. It is nevertheless conceivable that other alternative rules exist leading to the emergence of a similar scale invariant behaviour in a self-organized way.

For speciation to emerge as a result of sexual selection, a new sexual secondary trait must be introduced in a population (due to mutations, or to the activation of certain genes (Rutherford & Lindquist 1998), induced by new ecological conditions), or already existent displaying traits should be perceived differently under new ecological conditions (for instance, water quality or light intensity in a lake (Boughmann 2001)), differentiating the individuals. Thus, a population subject to new ecological conditions may be driven into a sexual selection process, from which, we will show, sympatric speciation can be the outcome.

Consider the limit case in which an individual would rank all the mates with the preferred trait first. In general, if a new trait (A or B) is introduced in a population, some individuals may prefer a mate with trait A, while others, one with trait B. We can classify males and females in four groups, AA, AB, BA and BB, the first letter standing for the trait, the second for the preference. If only males adopt active strategies towards mating, then both, AA and BA males, compete for the same females, AA or AB. A potential opportunity cost exists if, for instance, a BA male mates with an AA female, as there is a high probability that this mating will not last enough time to accomplish reproduction. We have checked that the mating instabilities reduce strongly the lifetime of couples where traits and preferences do not match. From the female’s point of view, AA and BA females prefer AA or AB males, which produces an indirect competition between AA and AB males. Neglecting the effects of recombination and assuming that the average number of males $`ij`$, $`m_{ij}`$, is equal to the average number of females with matched preferences, $`f_{ji}`$ , then $`m_{ij}`$ grows roughly with the probability of forming stable couples, $`dm_{ij}\alpha m_{ij}^2/m_j,`$ where $`m_j`$ is the total number of males with preference for $`j`$. Sexual selection will only keep those individuals reproducing at the highest rate, leading to a *competitive exclusion principle* with a symmetry breaking, where only AA$`\times `$AA and BB$`\times `$BB couples or BA$`\times `$AB and AB$`\times `$BA, are selected. Introducing recombination with one dominant locus for each, trait and preference, this competitive scenario can be tackled analytically as shown in the appendix. Then, only AA$`\times `$AA and BB$`\times `$BB couples are selected, which corresponds to the emergence of homotypic preferences.

Our general model considers a quantitative genetic preference, determined by several additive independent loci. The number of these loci does not need to be only one, as considered in the simple case in the appendix. If we consider a large number of loci controlling the preferences, then random configurations produce only small preference tendencies. The more biased arrangements of preference alleles are, the larger will be the number of individuals with the preferred trait ranked in the first positions. As discussed above, the function $`f_\nu (n_+)`$, establishes the relation between bias in preference loci and the intensity of discrimination. In figure 3a we show how sympatric speciation emerges, starting from a population with random genotypes. Even for species with a considerable number of independent genes controlling the preferences, speciation is remarkably achieved within just a few generations. Due to trait dominance, only AA or BB individuals last in the end, all heterozygous being eliminated. This corresponds to the emergence of two new species and a species recognition system. Prezygotic isolation is thus achieved giving way to species divergence. In this example, trait discrimination is the same for males and females (same $`\nu `$). However, for this example speciation could also be achieved if we had chosen higher trait discrimination for females (higher $`\nu `$, $`\nu _{females}=1.5`$ and $`\nu _{males}=0.5`$) as it is likely to happen in nature. In fact, trait discrimination needs to operate on both sexes, but its intensity may vary.

If $`f_\nu (n_+)`$ departs slowly from the unbiased arrangement of alleles corresponding to a difficult discrimination, disruptive sexual selection can be prevented (figure 3b). It is also possible to prevent speciation if traits are not straightfully inherited, due to complex dependencies on a large number of other genes (a trait may result from intricate interactions of a large number of proteins, leading to a complex inheritance) or as a result of mutations. These effects are taken into account through the m parameter. Large values of $`\mu `$ (for the example of figure 3a, $`\mu 6\times 10^3`$ mutations per generation and per allele) reduce the effects of the non-random pair formation created by the complex dynamics. Consequently above a certain threshold value speciation is prevented and a population with genotype frequencies in the Hardy-Weinberg proportions remains.

Our theory assumes that in a population, each individual presents many phenotypic traits, each contributing in a complicated way to establish the other individuals mate choice. Many of these traits may not have a strong impact in mate choice, in such a way that a phenotypic polymorphism remains in a population (they would be approximately neutral in what concerns sexual selection). In figure 4a we show a population in which discrimination relatively to a given trait is poor in the first 200 generations. The population consists of individuals with both traits and small preferences relatively to both traits. There is no relevant correlation between traits and preferences in the individuals. At generation 200, discrimination for the trait is increased and a quick speciation event emerges. An alternative picture is presented in figure 4b. Starting from a monomorphic population and sufficiently strong trait discrimination, mutations and a release of natural selection pressure (as a result of change of habitat, for instance) allows genetic drift. Preferences acquire variability, until individuals with preferences for both traits exist in the population. Then sympatric speciation emerges. In practice this later scenario needs more generations to reach the speciation disruptive configuration.

The theory presented here explains, under minimal assumptions, how sympatric speciation emerges or not, together with a species recognition system. It is a robust mechanism, not requiring special assumptions on the initial species (Doebeli & Dieckmann 2000; Kondrashov & Kondrashov 1999; Lande 1982; Turner & Burrows 1995; van Doorn & Weissing), on male-female incompatibilities (Gavrilets 2000; Parker & Partridge 1998), on the ecological conditions (Turner & Burrows 1995), or on non obvious connections between sexual and ecological traits (Doebeli & Dieckmann 2000; Felsenstein 1981; Maynard Smith 1966; Turelli et al. 2001; van Doorn & Weissing). Our theory assumes that males and females have preferences over the same traits. While female preferences have been extensively studied (Andersson 1994; Futuyma 1998), male preferences are usually neglected. Nevertheless, there is now mounting experimental evidence for the existence of male choice (Andersson 1994; Burley et al. 1982; Rutowski 1982; Roulin et al. 2000; Katvala & Kaitala 2001; Amundsen & Forsgren 2001). This is actually a natural outcome of natural selection as males would otherwise incur in potential opportunity costs: those having the ability to elicit healthier females, or predators mimicry, should certainly be favoured. However, in practice, they should be more difficult to detect, as males usually adopt active strategies towards mating, and thus are seen performing courtship to many females. This is certainly an issue deserving more research.

A unified understanding of sexual selection, starting from simple organisms like yeasts to complex organisms like mammals, may be attempted with our model. Indeed, in yeasts it is now clear that both mates in a couple choose their partners (Jackson & Hartwell 1990). Further, many genes can be involved in the several complex steps of this process (White & Rose 2001). Hence, analogues to the preference lists used in our model may already be built in. For this reason we propose yeasts as potential candidates to test our theory in the laboratory. And if, as predicted, robust sympatric speciation is observed, then many implications can result not only on a conceptual level. Indeed we believe that sympatric speciation could be used as a general framework for studying the interactions between receptors and pheromones, and to construct a better understanding on how these correlate with the intracellular state. This issue will be addressed in a forthcoming publication.

ACKNOWLEDGMENTS

We acknowledge Francisco Dionísio and Benoit Douçot for critical reading of the manuscript. We thank R.M.Conde for earlier collaboration, S. van Doorn, A. Luis, M. A. Santos and I. Pimentel for discussions.

Appendix: Mathematical approach of the sexual competitive scenario for a population with two locus, for traits and preferences, with a dominant allele each.

Here we outline a mathematical approach to the emergence of an exclusion principle from our results. The dynamical correlations in the model reduce strongly the fitness of non-optimal couples. Assume then that fitness is proportional to the probability of forming couples where traits and preferences match. In a two-locus, two-allele model, we have 9 genotypes, denoted by $`ij`$, respectively associated with displaying traits and preferences, with $`i,j=1,2,3`$. Genotype $`2`$ is heterozygous, and $`1`$ is a dominant homozygous. Denote the number of males (females) with genotype $`ij`$ by $`m_{ij}(f_{ij})`$. The evolution of the number of genotypes $`ij`$ is proportional to: $`dm_{ij}_{klmn}\alpha _{ij}^{klmn}m_{kl}f_{mn}/f_M`$, where $`f_M`$ is the number of females with trait corresponding to genotype mn $`(M=A,B)`$ and $`\alpha _{ij}^{klmn}`$ is a matrix giving the reproduction frequencies following Mendel rules. This equation is quite general but simplifies considerably in the present case, as most entries in the matrix are zero (see below). For instance, the equation associated with phenotype BB (genotype $`m_{33}`$) is: $`dm_{33}`$ $`m_{33}f_{33}/f_B+m_{23}f_{32}/(4f_B)+m_{32}f_{23}/(4f_B)+m_{22}f_{22}/(16f_A).`$

An exclusion principle emerges as genotypes try to grow the fastest possible (following similar growth equations) while interacting and as we concentrate on populations with a fixed number of individuals. Only those growing at the highest rate will remain.

If we sum all genotypes contributing to the same phenotype, we arrive to similar equations to the evolution of the phenotypes. It is possible to check that all terms contributing to $`dm_{\mathrm{𝐴𝐵}}=dm_{13}+dm_{23}`$ and $`dm_{\mathrm{𝐵𝐴}}=dm_{31}+dm_{32}`$ appear in $`dm_{\mathrm{𝐴𝐴}}=dm_{11}+dm_{12}+dm_{21}+dm_{22}`$ but with larger coefficients in the later. Thus $`m_{\mathrm{𝐴𝐵}}`$ and $`m_{\mathrm{𝐵𝐴}}`$ are eliminated (unless the starting population had only genotypes $`13`$ and $`31`$). Phenotypes $`m_{\mathrm{𝐴𝐴}}`$ and $`m_{\mathrm{𝐵𝐵}}`$ can coexist as both have a non-mixed contribution corresponding to the union of genotypes $`m_{11}+f_{11}`$ and $`m_{33}+f_{33}`$.

The non-zero matrix entries are: ($`\alpha _{ij}^{klmn}=`$ $`\alpha _{ij}^{mnkl}`$)

$`1/2=\alpha _{111121}=\alpha _{111211}=\alpha _{121112}=\alpha _{121212}=\alpha _{211121}=\alpha _{212121}`$

$`1/4=\alpha _{111122}=\alpha _{111212}=\alpha _{111221}=\alpha _{112121}=\alpha _{121122}=\alpha _{121221}=\alpha _{121222}=`$

$`=\alpha _{131212}=\alpha _{211122}=\alpha _{211221}=\alpha _{212122}=\alpha _{221122}=\alpha _{221221}=\alpha _{221222}=`$

$`=\alpha _{222122}=\alpha _{222222}=\alpha _{222332}=\alpha _{232332}=\alpha _{312121}=\alpha _{322332}=\alpha _{332332}`$

$`1/8=\alpha _{111222}=\alpha _{112122}=\alpha _{122122}=\alpha _{122222}=\alpha _{131222}=\alpha _{211222}=\alpha _{212222}=\alpha _{231222}=`$

$`=\alpha _{232222}=\alpha _{312122}=\alpha _{322122}=\alpha _{322222}`$

$`1/16=\alpha _{112222}=\alpha _{132222}=\alpha _{312222}=\alpha _{332222}`$

$`1=\alpha _{111111}=\alpha _{221331}=\alpha _{333333}`$

REFERENCES

Albertson, R.C., Markert, J.A., Danley, P.D. & Kocher, T.D. 1999 Phylogeny of a rapidly evolving clade: the cichlid fishes of Lake Malawi, East Africa. Proc. Natl. Acad. Sci. USA 96, 5107-5110.

Amundsen, T. & Forsgren, E. 2001 Male mate choice selects for female coloration in a fish. Proc. Natl. Acad. Sci. USA 98, 13155-13160.

Andersson, M. 1994 Sexual Selection. Princeton University Press.

Boughmann, J. W. 2001 Divergent sexual selection enhances reproductive isolation in sticklebacks. Nature 411, 944-948.

Burley, N., Krantzberg, G. & Raduran, P. 1982 Influence of colour-banding on the conspecific preferences of zebra finches. Anim. Behav. 27, 686-698.

Bush, G. L. 1969 Sympatric host race formation and speciation in frugivorous flies of the genus Rhagoletis (Diptera, Tephritidae). Evolution 23, 237-251.

Doebeli, M. & Dieckmann, U. 2000 Evolutionary branching and sympatric speciation caused by different types of ecological interactions. Am. Nat. 156, S77-S101.

Felsenstein, J. 1981 Skepticism towards Santa Rosalia, or why are there so few kinds of animals? Evolution 35, 124-138.

Futuyma, D.J. 1998 Evolutionary Biology. Sunderland: Sinauer.

Gale, D. & Shapley, L. S. 1962 College admission and the stability of marriage. Am. Math. Monthly 69, 9-14.

Gavrilets, S. 2000 Rapid evolution of reproductive barriers driven by sexual conflict. Nature 403, 886-889.

Jackson, C. L. & Hartwell, L. H. 1990 Courtship in S. cerevisae: both cell types choose mating partners by responding to the strongest pheromone signal. Cell 63, 1039-1051.

Jensen, H. J. 1998 Self-Organized Criticality. Cambridge Univ. Press.

Katvala, M. & Kaitala, A. 2001 Male choice for current female fecundity in a polyandrous egg-carrying bug. Anim. Behav. 61, 133-137.

Kondrashov, A. S. & Kondrashov, F. A. 1999 Interactions among quantitative traits in the course of sympatric speciation. Nature 400, 351-354.

Lande, R. 1982 Rapid origin of sexual isolation and character divergence in a cline. Evolution 36, 213-223.

Maynard Smith, J. 1966 Sympatric speciation. Am. Nat. 916, 637-650.

Mayr, E. 1963 Animal Species and Evolution. Harvard University Press.

Omero, M. J., Dzierzawa, M., Marsili, M. & Zhang, Y. - C. 1997 Scaling behavior in the stable marriage problem. J. Phys.I France 7, 1723-1729.

Panhuis, T. M., Butlin, R., Zuk, M. & Tregenza, T. 2001 Sexual selection and speciation. Trends Ecol. Evol. 16, 364-371.

Parker, G. A. & Partridge, L. 1998 Sexual conflict and speciation. Phil. Trans. R. Soc. London B 353, 261-274.

Rice, W. R. & Hostert, E. E. 1993 Perspective: Laboratory experiments on speciation: what have we learned in forty years?. Evolution 47, 1637-1653.

Roulin A, Jungi TW, Pfister H, Dijkstra C. 2000 Female barn owls (Tyto alba) advertise good genes. Proc R Soc Lond B 267, 937-941

Rutherford, S. L. & Lindquist, S. 1998 Hsp90 as a capacitor for morphological evolution. Nature 396, 336-342.

Rutowski, D.L. 1982 Epigamic selection by males as evidenced by courtship partner preferences in the checkered white butterfly (Pieris protodice). Anim. Behav. 30, 108-112.

Schliewen, U. K., Tautz, D. & Paabo, S. 1994 Sympatric speciation suggested by monophyly of crater lake cichlids. Nature 368, 629-632.

Takimoto, G., Higashi, M. & Yamamura, N. 2000 A deterministic genetic model for sympatric speciation by sexual selection. Evolution 54, 1870-1881.

Turelli, M., Barton, N. H. & Coyne, J. A. 2001 Theory and speciation. Trends Ecol. Evol. 16, 330-343.

Turner, G. F. & Burrows, M. T. 1995 A model of sympatric speciation by sexual selection. Proc. R. Soc. Lond. B 260, 287-292.

van Doorn, G. S. & Weissing, F. J. Ecological versus sexual selection models of sympatric speciation. Selection in press.

White, J.M. & Rose, M.D. 2001 Yeast mating: getting close to membrane merger. Curr. Biol. 11, R16-R20.