text:
# The Discovery of an Outburst and Pulsed X-ray Flux from SMC X-2 from RXTE Observations

## 1 Introduction

The first three X-ray point sources to be discovered in the Small Magellanic Cloud (SMC) were SMC X-1, SMC X-2 and SMC X-3, all of which are thought to be X-ray binaries. While SMC X-1 is a persistent pulsing source (e.g. Wojdowski et al. 1998 and references therein), SMC X-2 and X-3 are transient sources that were discovered in 1977 with SAS-3 (Clark et al. 1978) which were both initially very bright at a level of $``$1.0 $`\times `$ 1038 ergs s-1 (2 $``$ 11 keV) for assumed distances of 71 kpc and remained in outburst for approximately one month (Clark, Li, & van Paradijs 1979). Both sources were also detected with HEAO-1 A-2 experiment (Marshall et al. 1979). However, they have since exhibited little activity. SMC X-2 was not detected in an Einstein IPC survey (Seward & Mitchell 1980) and was seen in only one of two observations made with the ROSAT PSPC detector (Kahabka & Pietsch 1996, Haberl & Sasaki 2000) at a luminosity of 2.7 $`\times `$ 1037 ergs s-1 (0.15 $``$ 2.4 keV). SMC X-3 was also not detected by Seward & Mitchell (1980) but was seen in one ROSAT HRI observation (Haberl & Sasaki 2000). An optical counterpart for SMC X-2 was proposed by Crampton, Hutchings, & Cowley (1978) which was resolved into a pair of stars, “A” and “B”, with the fainter object (“B”), a main-sequence Be star, identified as the likely counterpart by Murdin, Morton, & Thomas (1979). Infra-red observations of SMC X-2 during a quiescent period are reported by Coe et al. (1997). In none of the previous X-ray observations of SMC X-2 or SMC X-3 were any pulsations detected although they are expected to be present as most Be star X-ray sources are X-ray pulsars. Another bright X-ray pulsar candidate within the SMC is the little studied source H0107-750 which also has a Be star as its proposed optical counterpart (Whitlock & Lochner 1994 and references therein).

Although no other SMC X-ray pulsars apart from SMC X-1 were known for many years there has recently been an explosion in the discovery of SMC X-ray pulsars. Hughes (1994) discovered a 2.7s pulsar, and three pulsars were found to be active in a single RXTE observation of the vicinity of SMC X-3 (Corbet et al. 1998). Further RXTE, ASCA and SAX observations have now lead to the identification of over 20 X-ray pulsars in the SMC (e.g. Haberl & Sasaki 2000, Yokogawa et al. 2000, and references therein) and their variability and several optical identifications suggests that at least the vast majority are transient Be star systems. The number of known Be star X-ray binaries in the SMC is thus now starting to rival the number known in the Galaxy (e.g. Bildsten et al. 1997).

In this paper we report the detection of activity from the vicinity of SMC X-2 in early 2000 with the Rossi X-ray Timing Explorer All-Sky Monitor (ASM). Subsequent observations with the RXTE Proportional Counter Array (PCA) then revealed the presence of 2.37s pulsations. Contemporaneous optical photometry of the optical counterpart are compared to the initial discovery photometry. Note that in this paper we use an assumed distance of 65 kpc to the SMC for consistency with several recent publications (e.g. Kahabka & Pietsch 1996, Haberl & Sasaki 2000) which is somewhat smaller than the 71 kpc adopted by Clark et al. (1978).

## 2 Observations and Analysis

### 2.1 RXTE All-Sky Monitor

The All-Sky Monitor (Levine et al. 1996) on board RXTE consists of three Scanning Shadow Cameras, sensitive to X-rays in an energy band of approximately 2-12 keV, which perform sets of 90 second pointed observations (“dwells”) so as to cover $``$80% of the sky every $``$90 minutes. The analysis presented here makes use of daily averaged light curves constructed from the flux measured in individual dwells. The Crab produces approximately 75 counts s-1 in the ASM over the entire energy range and ASM observations of blank field regions away from the Galactic center suggest that background subtraction is accurate to a systematic level about 0.1 counts s-1 (Remillard & Levine 1997). The RXTE All Sky Monitor provides light curves of several hundred sources including SMC X-2.

An examination of the light curve obtained from the ASM on SMC X-2, shown in Figure 1, suggested that the source may have increased in brightness in approximately January 2000. The source was still very faint and the possible increase in brightness was only visible in a very heavily smoothed and rebinned version of the ASM light curve. Because of this possible increase in flux we requested target of opportunity observations of the vicinity of SMC X-2 with the RXTE PCA.

### 2.2 Proportional Counter Array

The Proportional Counter Array (PCA; Jahoda et al. 1996) consists of five Proportional Counter Units (PCUs) sensitive to photons with energies between 2 - 60 keV with a total effective area of $``$6500 cm2. For the entire PCA across the complete energy band the Crab produces a count rate of about 13,000 counts s-1. The PCA spectral resolution at 6 keV is approximately 18% and the collimator gives a field of view of 1 full width half maximum. Due to instrumental problems not all PCUs are always operated and the observations reported here use varying numbers of PCUs.

Before the first set of scheduled observations of SMC X-2 was performed we noted the serendipitous detection of a source consistent with the location of SMC X-2 during slews over the SMC on 2000 April 9. In order to localize the position of the X-ray emission from this region more precisely, and so determine whether it was indeed consistent with the location of SMC X-2, we performed additional short slews in right ascension and declination over the position of SMC X-2 on 2000 April 12 together with short “stares”. Additional pointed observations were made between 2000 April 22 21:04 to April 23 02:43 and 2000 May 2 (06:09 to 06:48), May 3 (06:05 to 06:47 and 09:34 to 09:58) and May 5 (16:42 to 17:16). A log of the observations is given in Table 1. Data reduction used standard HEAsoft software, and background subtraction was performed using the model appropriate for faint sources.

### 2.3 Optical

The optical counterpart of SMC X-2 was observed with the South African Astronomical Observatory 1.0m telescope on 2000 April 14, two days after RXTE observation #2 was performed. Data were collected using the SITe4 CCD giving a field of approximately 5 arcminutes and a detector scale of 0.3 arcseconds per pixel. Observations were made through standard Johnson BV and Cousins RC filters. The data were reduced using IRAF and Starlink software, and the instrumental magnitudes were corrected to the standard system using E region standards. The resulting BVRC magnitudes are shown in Table 2, and all have an associated 1 $`\sigma `$ error of $`\pm `$0.02.

The results may be compared with the photometry of Murdin et al. (1979) which are also given in Table 2. These authors do not quote any errors on their results, but given that they only present their values to one decimal place they could well be $`\pm `$0.1. Hence the values for Star A (not the proposed candidate) are essentially in agreement between the two epochs.

## 3 Results

### 3.1 ASM

The ASM light curve (Fig. 1) suggests that there was an outburst from SMC X-2 that occurred between roughly 2000 January to April. However, the faintness of the source makes it difficult to determine precise start and stop times for this outburst. In addition, there were no SMC X-2 observations between approximately March 6 to April 1 (MJD 51609 to 51635). To obtain a crude characterization of the outburst in a quantitative way we fitted a model consisting of a constant flux plus a symmetric triangular profile outburst to the ASM light curve. This fit yields an outburst starting on MJD 51555.6 $`\pm `$ 12.2 (2000 January 13) and reaching a peak 47 $`\pm `$ 6 days later (February 28). However, since the gap in the ASM observations started at close to the fitted peak time, the time of real peak flux may actually have been later than this. The fitted maximum count rate is 1.0 $`\pm `$ 0.25 counts s-1 which, using the spectral parameters determined with the PCA (Section 3.2) corresponds to absorbed and unabsorbed 2 $``$ 10 keV luminosities of $``$1.5 and $``$1.7 $`\times `$ 1038 ergs s-1 respectively at 65 kpc. The peak flux during this outburst was thus at least comparable to, and may even have been higher than, the $``$0.8 $`\times `$ 1038 (d/65 kpc)2 ergs s-1 reported by Clark et al. (1978). Although the end of the outburst based on the fit of this simple model corresponds to roughly MJD 51649 (2000 April 15) the PCA observations (see next section) show activity after this time.

In an attempt to determine the orbital period of SMC X-2 we searched for periodic modulation in the ASM light curve by calculating power spectra but find no modulation at an amplitude $`>`$ 0.1 ASM counts s-1 ($``$1 mCrab) for periods longer than 0.5 days.

### 3.2 PCA

From the slews made over the location of SMC X-2 in observations #1 and #2 we obtain a source location of R.A. = 00h54m53s, decl. = $``$7338′ with a 90% confidence error circle of 3′ radius. This is consistent with the location of SMC X-2 at R.A. = 00h54m33s.4, decl. = $``$7341′ 04″(J2000; Murdin et al. 1979).

During the pointed observations we searched for pulsations by calculating power spectra from the PCA light curve and in observations #2 and #3 pulsations were very strongly detected at periods of approximately 2.37 seconds. In order to refine our estimates of the period we performed a pulse time arrival analysis. For each observation a template was constructed by folding the light curve on the period found from the power spectrum. The light curves were then divided into sections each of approximately equal duration with breaks at the ends of “good time” intervals and these sections of the light curve were also folded on this period. The phase shifts between the folded profiles and the template were then calculated by cross-correlation. A linear fit was made to the phase shifts found and this was used to re-determine the period and the associated one $`\sigma `$ errors. In this way pulse periods of 2.371532 $`\pm `$ 0.000002 s and 2.371861 $`\pm `$ 0.000003 s were measured from observations #2 and #3 respectively. The PCA light curve from observation #3 folded on the pulse period in four energy bands is plotted in Figure 2.

For the two pointed observations when the source was detected (#2 and #3) we fit the background subtracted spectrum between 3 to 15 keV with an absorbed power-law model. The results of these fits are given in Table 3. The unabsorbed 2 – 10 keV luminosities implied by these fits are $``$3.3 $`\times `$ 1037 (d/65 kpc)2 ergs s-1 and $``$2.7 $`\times `$ 1037 (d/65 kpc)2 ergs s-1 respectively.

In the last two observations made on May 2-3 (#4) and May 5 (#5) no pulsations and no significant flux were detected. Given the large field of view of the PCA and the possibility of contaminating sources the most stringent upper limit on the source flux probably comes from the lack of detection of any pulsations rather than the nominally background subtracted spectrum. Based on the lack of pulsations, and scaling from the two detections in observations #2 and #3, we place upper limits on the source flux of $``$0.38 $`\times `$ 10-11 ergs cm-2 s-1 for observation #4 and $``$0.6 $`\times `$ 10-11 ergs cm-2 s-1 for observation #5. These flux limits corresponds to luminosities of $`<`$ $``$1.8 $`\times `$ 1036 (d/65 kpc)2 ergs s-1 and $`<`$ $``$2.9 $`\times `$ 1036 (d/65 kpc)2 ergs s-1 respectively.

### 3.3 Optical

From Table 2 it can be seen that the proposed optical counterpart of SMC X-2 (Star B) was approximately 0.5 magnitudes fainter in April 2000 than in 1978 and somewhat redder, which is likely indicative of less circumstellar disk emission. However, when Murdin et al. performed their optical observations on 1978 June 15 the source had already faded below the detectable X-ray level: Clark et al (1979) had observed SMC X-2 again on 1977 December 7-15 - just 2 months after their original detection and failed to detect it with a 3$`\sigma `$ upper limit of 1.0 $`\times `$ 1037 ergs s-1. In contrast we measured a flux of $``$3.3 $`\times `$ 1037 ergs s-1 just two days before our optical observations. The X-ray and optical brightness of SMC X-2 may thus not be correlated in a simple fashion. For neither the 1978 or the 2000 outburst do we have optical photometry obtained at the peak of the outbursts and it is possible that the optical counterpart may have been brighter during the peak.

## 4 Discussion

Given the large number of X-ray pulsars in the SMC can we be certain that the new pulsar that we have detected is indeed SMC X-2? The pulse period is rather short and we can make use of the result of Stella, White & Rosner (1986) who noted an anti-correlation between pulse period and maximum X-ray luminosity. SMC X-2 at its peak observed with SAS-3 and the RXTE ASM was rather bright at $`>`$ $``$1038 ergs s-1 which, together with the pulse period that we measure, would be consistent with the relationship found by Stella et al.

The vicinity of SMC X-2 may also be less densely populated with active pulsars than the region around SMC X-3. After emission had been detected from the vicinity of SMC X-3 with RXTE subsequent analysis revealed the presence of three pulsars (Corbet et al. 1998), while ASCA observations showed that none of these were in fact SMC X-3 itself. However, in the case of SMC X-2, we only detect one pulsar and, in addition, its location is consistent at the $``$3 arc minute level with SMC X-2. The brightening observed in the ASM light curve for SMC X-2 is also consistent with the source being not more than a few arcminutes from SMC X-2. Further evidence still for the identification of the 2.37s pulsar with SMC X-2 comes from the detection of this source by Torii et al. (2000) with the Gas Imaging Spectrometer experiments on board ASCA. The GIS is an imaging instrument which gives source locations with 90% errors of 1 to 2 arc minutes radius (Tashiro et al. 1995). In the ASCA GIS survey of the SMC (Yokogawa et al. 2000) 12 X-ray pulsars and 8 candidate X-ray pulsars were identified. We thus obtain an approximate estimate of the chance of finding any X-ray pulsar or candidate within 2 arc minutes of SMC X-2 in a snap-shot survey of a few percent. Another independent estimate of the probability that the source that we have observed is SMC X-2 comes from the serendipitous RXTE slew detection (observation #1). The slews across the SMC that occur during routine RXTE observations provide a monitor for bright new sources within this entire region. There are currently only five sources known in the SMC with luminosities $``$1038 ergs s-1: SMC X-1, SMC X-2, SMC X-3, H0107-750, and the 31 second pulsar XTE J0111.2-7317 (Chakrabarty et al. 1998). The probability of finding a new source with a comparable brightness at a position consistent with one of these known bright sources within the approximately 2 $`\times `$ 2 extent of the SMC is less than 1%.

Together, the comparable luminosity observed with RXTE and SAS-3 and the positional coincidence at the few arc minute level, strongly suggest that the source we have observed is indeed SMC X-2. The evidence would be strengthened further still if the position of the 2.37 s pulsar could be localized more precisely or evidence found for pulsations at this period in the archival SAS-3 or HEAO-1 observations.

The outburst we have observed with RXTE could be due to either changes in the circumstellar envelope of the Be star, associated with orbital variability, or both of these effects. We note that no other outburst of a similar size is apparent in the ASM light curve which would argue for a short lived change in the circumstellar envelope. If SMC X-2 has properties consistent with the general correlation between orbital and pulse periods exhibited by most Be/neutron star binaries (Corbet 1986) then an orbital period of order 20 days to within a factor $``$2 might be expected. However, no evidence for any periodic variability is found in the ASM light curve.

During the course of our observations the pulse period of SMC X-2 increased and was longer still at 2.37230 $`\pm `$ 0.00004 s when ASCA observed on April 25-26 (Torii et al. 2000). This change in pulse period may be due to orbital motion and would correspond to a velocity change of $``$100 km s-1. While there are currently too few pulse period measurements to place any strong constraints on the orbital parameters, if SMC X-2 returns again to an active state it is possible that additional observations may yield the orbital period, either from monitoring the pulse period or from periodic outbursts that the ASM could observe.

The pulse period that we find for SMC X-2 is among the shortest currently known for high mass X-ray binaries (e.g. Bildsten et al. 1997) especially for those systems which contain Be star primaries. As noted above, this short pulse period is consistent with the high source luminosity. SMC X-3, which also exhibited a high maximum luminosity, may also be expected to have a relatively short pulse period. Previously Clark et al. (1978, 1979) argued that the high luminosities of SMC X-1, X-2 and X-3 implied that the luminosity distribution of the SMC X-ray sources is shifted toward higher luminosities compared to Galactic systems. However, now that the transient X-ray pulsar population of the SMC is known to be so large this shift to higher luminosities is called into question (see e.g. Haberl & Sasaki 2000) as the earlier observations were only detecting the high-luminosity tip of the large SMC X-ray pulsar population. Extensive studies of these objects down to low flux levels are required to fully establish their overall properties as a class.

This paper made use of quick look data provided by the RXTE ASM team at MIT and GSFC. We thank our colleagues in the RXTE team at GSFC and MIT for useful discussions on several aspects of RXTE data analysis.

Figure Captions


format: text

text:
# Assisted stellar suicide: the wind-driven evolution of the recurrent nova T Pyxidis

## 1 Introduction

The recurrent nova T Pyx is the most prolific of all recurrent novae, with a mean outburst recurrence time of only 22 years. But even away from outburst, T Pyx’s observational characteristics are remarkable. Most importantly, its quiescent bolometric luminosity is $`L_{\mathrm{bol}}\stackrel{>}{}10^{36}`$ ergs s-1 (Patterson et al. 1998). This is much higher than expected for a short-period cataclysmic variable (CV) like T Pyx, whose orbital period of $`P_{\mathrm{orb}}=1.8`$ hr places it well below the period gap.

The implied disagreement with theoretical expectations is severe. For reasonable estimates of T Pyx’s system parameters, the standard model of CV evolution predicts mass transfer to be driven by gravitational radiation at a rate of $`\dot{M}_{\mathrm{gr}}4\times 10^{11}\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$. The corresponding accretion luminosity is only $`L9\times 10^{32}`$ erg s-1. Thus the gap between T Pyx’s observed luminosity and the theoretically predicted value spans at least three orders of magnitude.

In order to account for T Pyx’s extreme luminosity without having to postulate a similarly extreme mass-transfer rate for the system, it has been suggested that, even in quiescence, the system may be powered predominantly by nuclear burning on the surface of the white dwarf (WD), rather than by accretion (e.g. Webbink et al. 1987; Patterson et al. 1998). However, even this can only reduce the luminosity gap by about an order of magnitude (c.f. Section 2.2). Thus regardless of the actual mode of radiation energy release, T Pyx’s extreme luminosity implies an accretion rate $`\dot{M}_{\mathrm{acc}}\stackrel{>}{}(440)\times 10^9\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$. This lower limit is in line with estimates based on theoretical nova models: in order to reproduce T Pyx’s short outburst recurrence time, an accretion rate $`\dot{M}_{\mathrm{acc}}\stackrel{>}{}10^8\mathrm{M}_{\mathrm{}}\mathrm{yr}^1`$ appears to be required (Livio & Truran 1992; Prialnik & Kovetz 1995; Contini & Prialnik 1997).

Clearly, T Pyx poses a fundamental challenge to our understanding of CV evolution. Put concisely, the question is this: why is T Pyx’s accretion rate so abnormally high relative to “ordinary” CVs? The answer we propose in this Letter is that T Pyx has left the standard CV evolutionary track completely and is instead currently evolving as a wind-driven supersoft X-ray source.

## 2 T Pyx as a Wind-driven Supersoft X-ray Source

### 2.1 Motivation

T Pyx’s observational characteristics in quiescence are clearly extreme among CVs. However, as recently pointed out by Patterson et al. (1998) they are actually a good match to the properties of the supersoft x-ray sources (SSSs). Despite this empirical connection, the standard model for binary SSSs – thermal-timescale mass transfer from a companion initially more massive than the WD – cannot work for T Pyx: the shortest orbital period found in the detailed study of this phase by Deutschmann (1998) is about 6 hr, much longer than T Pyx’s 1.8 hr. In fact, King et al. (2000) find that the two shortest-period systems among known binary SSSs – RX J0537.7$``$7034 ($`P_{\mathrm{orb}}3.5`$ hr) and SMC 13 ($`P_{\mathrm{orb}}4.1`$ hr) – are already difficult to account for within the thermal-timescale mass-transfer framework, even allowing for non-conservative evolution.

We therefore propose that T Pyx is instead a member of an (even more) exotic class of SSSs – the wind-driven supersoft sources. The existence of this class has recently been suggested by van Teeseling & King (1998, hereafter VK98; see also King & van Teeseling 1998, hereafter KV98). These authors show that if the secondary star in a close binary system like T Pyx is strongly irradiated by soft x-rays, a powerful wind may be driven from its surface. Under certain conditions, the mass and angular momentum loss in this wind can dominate the binary evolution, which in turn can drive mass transfer at a rate comparable to the wind mass-loss rate. If the wind-driven mass-transfer rate is high enough to power the irradiating source, a self-sustaining, stable, wind-driven state is created. All that is required in this scenario is an event to trigger the wind-driven evolution.

### 2.2 Theoretical Background

We take as our starting point Equation 7 in VK98, which gives the expected mass-loss rate from a low-mass, near-main-sequence secondary star irradiated by intense soft x-ray emission from the vicinity of the primary as

$$\dot{M}_{w2}3\times 10^7\varphi \frac{r_2}{a_{11}}(m_2L_{37}\eta _s)^{1/2}\mathrm{M}_{\mathrm{}}\mathrm{yr}^1.$$
(1)
In this equation, $`m_2=M_2/M_{\mathrm{}}`$ and $`r_2=R_2/R_{\mathrm{}}`$ are the mass and radius of the secondary, $`a_{11}=a/(10^{11}\mathrm{cm})`$ is the binary separation and $`L_{37}=L/(10^{37}\mathrm{erg}\mathrm{s}^1)`$ is the luminosity of the irradiating source. The factor $`\eta _s`$ measures the efficiency of the irradiating spectrum in producing the wind, and $`\varphi `$ measures the area of the mass-losing regions on the secondary’s surface, relative to $`\pi R_2^2`$.

The irradiating luminosity is most conveniently parameterized in terms of the luminosity expected for steady shell burning, $`L_{\mathrm{sb}}`$, which is

$$L_{\mathrm{sb}}=2.9\times 10^{37}\left(\frac{\dot{M}_{\mathrm{acc}}}{10^7\mathrm{M}_{\mathrm{}}\mathrm{yr}^1}\right)\mathrm{erg}\mathrm{s}^1$$
(2)
(Iben 1982). Defining $`L_{\mathrm{sb},37}=L_{\mathrm{sb}}/(10^{37}\mathrm{erg}\mathrm{s}^1)`$, we therefore write the irradiating luminosity as

$$L_{37}=\eta _aL_{\mathrm{sb},37}=2.9\eta _a\left(\frac{\dot{M}_{\mathrm{acc}}}{10^7\mathrm{M}_{\mathrm{}}\mathrm{yr}^1}\right)$$
(3)
and use the numerical factor $`\eta _a1`$ as a measure of the efficiency with which the accreted material is actually converted into luminosity. Thus $`\eta _a=1`$ for steady shell burning, whereas $`\eta _a=L_{\mathrm{acc}}/L_{\mathrm{sb}}\stackrel{<}{}0.3`$ for accretion onto a massive WD in the absence of any nuclear processing ($`L_{\mathrm{acc}}=GM_1\dot{M}_{\mathrm{acc}}/R_1`$). Combining Equations 1 and 3 yields

$$\begin{array}{c}\dot{M}_{\mathrm{acc}}=g\dot{M}_{w2}\hfill \\ 1.2\times 10^6g^2\varphi ^2\eta _a\eta _s\left(\frac{q^{5/2}}{1+q}\right)^{2/3}m_1\mathrm{M}_{\mathrm{}}\mathrm{yr}^1.\hfill \end{array}$$
(4)
Here, $`g=\dot{M}_{\mathrm{acc}}/\dot{M}_{w2}`$ is the dimensionless accretion rate (measured in units of the wind mass-loss rate), $`m_1=M_1/M_{\mathrm{}}`$, and we have used Paczyński’s (1971) approximation for $`R_{L2}/a`$ to recast $`r_2/a_{11}`$ in terms of $`q`$. Here, $`R_{L2}`$ denotes the volume-averaged Roche-lobe radius of the secondary. If a stationary, stable, wind-driven state exists, the accretion rate in it must be given by Equation 4.

In order to determine $`g`$ as a function of the system parameters, one needs to calculate the mass-transfer rate due to Roche-lobe overflow (RLOF) in a semi-detached close binary system whose evolution is driven by a stellar wind from the secondary. This calculation is described in detail by VK98 and KV98. Briefly, such a system will quickly settle in a stationary state with $`\dot{R_{L2}}/R_{L2}=\dot{R_2}/R_2`$, provided such a state exists and is stable. The stationarity condition is sufficient to find $`g`$ as

$$g=\frac{(6\beta _2+2q)(5+3\zeta )(1+q)}{(1+q)(5+3\zeta 6q)},$$
(5)
where $`\zeta `$ is the effective mass-radius index of the secondary (describing its reaction upon mass loss), and $`\beta _2`$ is the specific angular momentum of the escaping wind material, measured relative to the specific orbital angular momentum of the secondary.

Equation 5 gives the dimensionless RLOF mass-transfer rate in the stationary wind-driven state. It is identical to Equation 30 of VK98, except that we have retained the explicit dependence on $`\beta _2`$ (VK98 set $`\beta _2=1`$ throughout). As shown by KV98, the stationary solution is stable for system parameters appropriate to T Pyx and $`\dot{M}_{w2}L^{1/2}\dot{M}_{\mathrm{acc}}^{1/2}`$. In deriving Equation 5, it has been assumed that all of the material in the stellar wind from the secondary escapes. Also, the effects of mass loss from the primary have been ignored. In reality, T Pyx may undergo significant mass ejection during its nova eruptions, and the long term average $`\dot{M}_{w1}`$ may be comparable to $`\dot{M}_{w2}`$. However, we have verified that even strong, episodic mass loss from the primary has virtually no impact on our results for T Pyx in quiescence, unless the specific angular momentum carried away by the nova ejecta is extremely high (much higher than that of the primary). In that case, the angular momentum loss associated with nova eruptions would further accelerate (and could possibly even dominate) the binary evolution. We will return to this possibility in Section 3.

We finally consider the orbital period derivative of such a wind-driven system. If this is entirely due to stationary mass loss and/or transfer, we can combine Kepler’s 3rd law with Paczyński’s (1971) approximation for $`R_{L2}`$ and differentiate logarithmically to obtain

$$\frac{\dot{P}_{\mathrm{orb}}}{P_{\mathrm{orb}}}=\frac{3\zeta 1}{2}\frac{\dot{M}_2}{M_2}=\frac{(3\zeta 1)(1+g)}{2}\frac{\dot{M}_{w2}}{M_2}$$
(6)
(c.f. Equation 34 in VK98). However, the period derivative measured in real systems may differ from this, since the timescale associated with such measurements is much shorter than the timescale on which the RLOF rate can adjust itself ($`\tau _{\mathrm{RLOF}}\frac{H}{R_2}\frac{M_2}{\dot{M}_{w2}}10^4`$ yr, where $`H/R_210^4`$ is the scale height in the atmosphere of the secondary near the inner Lagrangian point). Thus observationally determined period changes could, for example, be due to fluctuations of $`\dot{M}_{w2}`$ or $`\beta _2`$ on timescales shorter than $`\tau _{\mathrm{RLOF}}`$.

### 2.3 Application to T Pyx

We are now ready to apply the wind-driven evolution scenario to T Pyx. For definiteness in our numerical estimates, we will adopt main-sequence-based values for the mass and radius of T Pyx’s secondary: $`M_2=0.12M_{\mathrm{}}`$, $`R_2=0.17R_{\mathrm{}}`$. These follow from Kepler’s law, the orbital period, Patterson’s (1998) power law approximation to the the M-dwarf mass-radius relationship of Clemens et al. (1998) and Paczyński’s (1971) approximation for $`R_{L2}R_2`$ of a Roche-lobe filling secondary star. In addition, we will use $`M_1=1.2M_{\mathrm{}}`$ as an estimate of the white dwarf mass (Contini & Prialnik 1997). We therefore take the mass ratio in T Pyx to be $`q=M_2/M_1=0.1`$. We note from the outset that even this rather low value for $`q`$ could be an overestimate, since T Pyx’s wind-driven evolution may already have reduced $`M_2`$ below its main-sequence value. Our assumption that the secondary is still close to the main sequence amounts to saying that wind-driving has only just begun.

We begin by noting that, for a low-mass secondary undergoing adiabatic mass loss, we may take $`\zeta 1/3`$ (VK98). Next, we consider two extreme estimates for the angular momentum loss parameter $`\beta _2`$. To obtain a lower limit, we note that the stellar wind material will carry away at least the specific angular momentum of the secondary, in which case $`\beta _2=1`$ (this is the case considered by VK98). On the other hand, the stellar wind material may extract angular momentum from the binary system by frictional processes. An upper limit to the amount of specific angular momentum that is likely to be extracted this way is given by the specific angular momentum of particles escaping through the outer Lagrangian points, which is $`j_21.65a^2\mathrm{\Omega }`$ (Sawada et al. 1984; Flannery & Ulrich 1977; Nariai 1975). Here $`\mathrm{\Omega }`$ is the angular velocity of the binary system, and, in our notation, the corresponding angular momentum loss parameter is $`\beta _2=2`$. These estimates, together with $`q=0.1`$, yield $`g0.5`$ ($`\beta _2=1`$) and $`g2`$ ($`\beta _2=2`$). Substituting these values back into Equation 4, along with $`m_1=1.2`$, we obtain

$$\dot{M}_{\mathrm{acc}}(\mathrm{T}\mathrm{Pyx})\{\begin{array}{cc}7\times 10^9\varphi ^2\eta _s\eta _a\mathrm{M}_{\mathrm{}}\mathrm{yr}^1\hfill & \beta _2=1\hfill \\ & \\ 1\times 10^7\varphi ^2\eta _s\eta _a\mathrm{M}_{\mathrm{}}\mathrm{yr}^1\hfill & \beta _2=2\hfill \end{array}$$
(7)
for T Pyx if it has settled into a stationary, wind-driven state. The corresponding luminosity is given by Equation 3 as

$$L(\mathrm{T}\mathrm{Pyx})\{\begin{array}{cc}2\times 10^{36}\varphi ^2\eta _s\eta _a^2\mathrm{erg}\mathrm{s}^1\hfill & \beta _2=1\hfill \\ & \\ 3\times 10^{37}\varphi ^2\eta _s\eta _a^2\mathrm{erg}\mathrm{s}^1\hfill & \beta _2=2.\hfill \end{array}$$
(8)
Thus wind-driving can indeed account for T Pyx’s extreme accretion rate and luminosity, provided that the various efficiency factors in Equations 7 and 8 are not too far from unity.

VK98 have argued that $`\eta _s1`$ in SSSs, since soft x-rays are absorbed well above the photosphere and should therefore be quite efficient at driving the wind. The factor $`\varphi `$ is just the area of the mass-losing regions on the secondary divided by $`\pi R_2^2`$, so that mass loss from the entire front hemisphere would correspond to $`\varphi =2`$. We may therefore also expect $`\varphi 1`$, even if the secondary is partially shielded by an optically thick accretion disk. However, the most interesting parameter in this context is $`\eta _a`$. As noted in Section 2.2, energy release by accretion yields $`\eta _a\stackrel{<}{}0.3`$. The upper limit in this inequality corresponds to a Carbon-core WD of maximum mass ($`1.4`$ M) and minimum radius ($`0.002R_{}`$; Hamada-Salpeter 1961). For a 1.2 $`M_{}`$ WD on the same mass-radius relation, we have $`\eta 0.1`$. Thus gravitational energy release alone can only meet the system’s luminosity requirements if the WD is even more massive than we have assumed and $`\beta _2>1`$.

The alternative is that nuclear processing continues in T Pyx even in quiescence. This would imply that nuclear burning in T Pyx occurs both quasi-steadily (in quiescence) and explosively (during outbursts), with the former taking place at a rate slightly below the accretion rate. From an empirical point of view, this does not seem unreasonable: as noted by Patterson et al. (1998), some SSSs in the LMC and M31 are recurrent (see Kahabka 1995), and the galactic novae GQ Mus and V1974 remained luminous soft x-ray sources for several years after their nova eruptions (Shanley et al. 1995; Krautter et al. 1996). Symbiotic stars (SySs) provide another interesting point of comparison. Several classical SySs, such as Z And, exhibit erupting behaviour even though though their quiescent luminosities ($`L10^3L_{}`$) suggest nuclear processing as the dominant power source (Mürset, Nussbaumer, Schmid & Vogel 1991). In addition, the luminosities of symbiotic novae (a distinct class from the classical SySs) remain high for decades after their outbursts (Mürset & Nussbaumer 1994).

From a theoretical point of view, the situation is somewhat more difficult. The problem is that the rate at which steady nuclear processing can proceed is not, in principle, a tunable parameter. Based on 1-dimensional models, explosive and steady processing are generally expected to be mutually exclusive regimes, whose dividing line is a function of accretion rate and white dwarf mass (Iben 1982). Thus steady burning is not expected to occur for accretion rates less than $`\dot{M}_{\mathrm{crit}}1.32\times 10^7M_{\mathrm{WD}}^{3.57}`$. The accretion rates predicted by our wind-driven model are below this line (albeit by only a factor of 2.5 for $`\beta =2`$ and $`M_{\mathrm{WD}}=1.2M_{}`$).On the other hand, nuclear processing on real WDs is unlikely to be well-described by spherically symmetric models. For example, thermonuclear runaways (TNRs) triggered by localized temperature or pressure fluctuations may not always spread and evolve into global TNRs (Shara 1982; Shankar & Arnett 1994). Quasi-steady quiescent burning in T Pyx could therefore conceivably be the collective outcome of many successive localized TNRs. The individual mini-eruptions would not necessarily be obvious observationally, provided they are sufficiently small and frequent. These localised TNRs could process non-degenerate material more slowly than it accretes, leading to the build-up of a more and more massive non-degenerate envelope and, eventually, to a global TNR.

The preceding is a highly speculative scenario, and we do not mean to endorse it too strongly. We have outlined it mainly to provide a specific illustration of the general idea that (quasi-)steady and explosive nuclear processing might take place in the same object. This general idea is not new. It was first suggested by Webbink et al. (1986) and again by Patterson et al. (1998), both times without any specific model in mind. The main achievement of the wind-driving mechanism is that, given a high radiative efficiency, T Pyx’s high mass accretion rate and luminosity can be accounted for self-consistently. We do, of course, acknowledge that the requirement $`\eta _a1`$ implies a fair amount of fine-tuning, regardless of whether one invokes an extremely massive WD or quiescent nuclear burning. But some theoretical fine-tuning seems reasonable for a system like T Pyx, whose short orbital period, high luminosity and ability to produce nova eruptions are a unique combination among CVs and SSSs. This is not to say that the wind-driven evolutionary channel must be narrow: most other wind-driven systems may be characterized by somewhat higher accretion rates than T Pyx and may therefore be steady SSSs. Observationally, such steady, wind-driven SSSs may nevertheless be rare, since their evolutionary timescales would be even shorter than T Pyx’s.

We finally turn to the orbital period derivative. On substituting the values for $`g`$ and $`\dot{M}_{w2}=g\dot{M}_{\mathrm{acc}}`$ derived above into Equation 8 and inverting, we find that the expected timescale for period increase due to stationary wind-driven mass loss and mass transfer in T Pyx is

$$\frac{P_{\mathrm{orb}}}{\dot{P}_{\mathrm{orb}}}\left(\mathrm{T}\mathrm{Pyx}\right)\{\begin{array}{cc}5\times 10^6\mathrm{yr}\hfill & \beta _2=1\hfill \\ & \\ 7\times 10^5\mathrm{yr}\hfill & \beta _2=2.\hfill \end{array}$$
(9)
Patterson et al. (1998) derived $`P_{\mathrm{orb}}/\dot{P}_{\mathrm{orb}}3\times 10^5`$ yr from the periodic dip in T Pyx’s optical light curve, but more recent timings suggest a slower rate of change. In any case, as noted in the discussion following Equation 6, the observed $`\dot{P}_{\mathrm{orb}}`$ may not be a valid estimate of the stationary (long-term average) value.

## 3 Discussion and Conclusions

We have shown that the abnormally high luminosity of the short-period recurrent nova T Pyx can be explained if the system belongs to the class of wind-driven SSSs. T Pyx’s accretion rate is then much higher than that of an ordinary CV, because the system’s evolution is dominated by a radiation-induced wind from the secondary star. The wind-driven evolution is self-sustaining, because the wind-induced accretion rate is sufficient to power the luminosity that excites the wind.

Our model predicts that T Pyx is an intrinsically bright soft x-ray source. T Pyx was observed serendipitously as part of the ROSAT all-sky survey, but not detected (Verbunt et al. 1997). This is not yet in conflict with a SSS model for T Pyx, because of the limited sensitivity of these observations and the relatively high column density towards the system ($`N_\mathrm{H}2.1\times 10^{21}`$ cm-2; Dickey & Lockman 1990111Data available online at http://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl.). However, a sensitive, pointed observation with Chandra or XMM could yield a detection and provide a direct estimate of the WD luminosity and temperature. For example, assuming a distance of 3 kpc (e.g. Patterson et al. 1998), the flux produced by a blackbody with $`L_{\mathrm{WD}}=10^{37}`$ ergs/s and $`T_{\mathrm{WD}}=2.4\times 10^5`$ K is below the ROSAT detection limit, but would provide around 1000 counts in a 30 ksec observation with XMM.

If T Pyx is, in fact, a wind-driven SSS, its current evolutionary timescale is $`\tau _{\mathrm{evol}}=|M_2/\dot{M}_2|=\dot{P}_{\mathrm{orb}}/P_{\mathrm{orb}}10^6`$ yrs. Clearly, this is an extremely short-lived evolutionary state. So how will T Pyx’s evolution actually end? We can think of at least three possible outcomes. First, the mass ratio of the system may be driven so low that irradiation and wind-driving cease to be effective. Second, the wind-driven evolution may proceed until the secondary completely evaporates, leaving T Pyx as an isolated WD. Third, if the total mass of the system is somewhat higher than we have assumed, the WD may reach the Chandrasekhar limit and explode in a Type Ia supernova.

The likelihood of the last outcome depends on whether the amount of material that is ejected in each of T Pyx’s recurrent nova outbursts is significantly less than that which is accreted during each quiescent interval. Theoretically, this is plausible (e.g. Prialnik & Kovetz 1995; Livio & Truran 1992). The question might also be addressed observationally: the total mass accreted during quiescence is $`\mathrm{\Delta }M_{\mathrm{acc}}\dot{M}_{\mathrm{acc}}\tau _{\mathrm{rec}}120\times 10^7\mathrm{M}_{\mathrm{}}`$. The mass ejected during an outburst can be inferred by comparing the pre- and post-outburst orbital periods, although this is sensitive to the details of the ejection process (Livio 1991; Livio et al. 1991). Such a comparison could also distinguish between wind-driven and nova-driven evolution scenarios. As noted in Section 2.2, if the specific angular momentum of the mass ejected during T Pyx’s nova eruptions is very high ($`j_{\mathrm{nova}}\stackrel{>}{}a^2\mathrm{\Omega }>>j_1`$), the angular momentum loss associated with nova eruptions could, in principle, itself drive T Pyx’s high mass transfer rate. In the context of this model, the difference between pre- and post-outburst orbital periods must be relatively large. By contrast, the wind-driving model makes no such prediction.

Given the strong possibility that the system will destroy itself on a short timescale, T Pyx represents an evolutionary channel by which short-period CVs may be removed from the general CV population. The existence of such a channel is desirable from a theoretical point of view. As most recently pointed out by Patterson (1998), the observed CV population shows a significant dearth of short-period systems relative to the predictions of standard evolutionary scenarios. Whether T Pyx’s ”assisted stellar suicide” represents a sufficiently wide channel in this context is an important question: do most CVs eventually enter a wind-driven phase or is T Pyx really a unique system? The answer clearly depends on how wind-driving is actually triggered, an issue we have not addressed in the present work. One possibility is that triggering occurs during a period of residual nuclear burning in the aftermath of a nova eruption, though this would have to be limited to systems satisfying other constraints (e.g. low $`q`$ and/or high $`M_{\mathrm{WD}}`$).

###### Acknowledgements.

We are grateful to Jim Truran for a clear, concise and cogent referee’s report. Support for CK was provided by NASA through Hubble Fellowship grant HF-01109 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555. ARK gratefully acknowledges a PPARC Senior Fellowship.


text:
# 1 Graphs for K⁺→π⁺e⁺e⁻ through a virtual photon. (a) and (b) are long-distance graphs, (c) is a short-distance graph and (d) is a pion loop term. In each graph, the blob denotes the weak (strangeness-changing) vertex and the wavy line is an off-shell photon.

Studying the scale and $`q^2`$ dependence of $`K^+\pi ^+e^+e^{}`$ decay

H. Burkhardt***Also at Shell Centre, School of Education, University, Nottingham NG8 1BB, England and J. LoweAlso at Physics Department, University, Birmingham B15 2TT, England

Physics Department, University of New Mexico, Albuquerque, NM 87131, USA

G. Eilam and M.D ScadronPermanent address: Physics Department, University of Arizona, Tucson, AZ 85721, USA

Physics Department, Technion, 32000 Haifa, Israel

## Abstract

We extract the $`K^+\pi ^+e^+e^{}`$ amplitude scale at $`q^2=0`$ from the recent Brookhaven E865 high-statistics data. We find that the $`q^2=0`$ scale is fitted in excellent agreement with the theoretical long-distance amplitude. Lastly, we find that the observed $`q^2`$ shape is explained by the combined effect of the pion and kaon form-factor vector-meson-dominance $`\rho `$, $`\omega `$ and $`\varphi `$ poles, and a charged pion loop coupled to a virtual photon$`e^+e^{}`$ transition.

PACS numbers:13.20.Eb and 13.25.Es

I. INTRODUCTION

The decay $`K^+\pi ^+e^+e^{}`$ has been studied theoretically for many years. There is general agreement that the process is dominated by the long-distance “bremsstrahlung” graphs, in which a virtual photon is radiated by the charged kaon or the pion, away from the strangeness-changing weak vertex. However, qualitative agreement with the available experimental evidence on the branching ratio and the form factor (i.e. the $`q^2`$ dependence) has not previously been achieved. Here, the virtual photon mass, $`q^2`$, is given by $`(p_++p_{})^2`$, where $`p_\pm `$ are the four-momenta of the $`e^\pm `$. Also, most calculations have several unknown parameters, making firm predictions difficult. Experimentally, the situation has been dramatically improved by the recent availability, from Brookhaven experiment E865 , of accurate data on both the branching ratio and the $`q^2`$ dependence. In this paper, we fit these data to about 5% in a calculation based on a specific model with no free parameters (except for relative signs), which fits other related data.

II. AMPLITUDE AT $`q^2=0`$

The E865 experiment  at Brookhaven had measured a new value of the branching ratio for $`K^+\pi ^+e^+e^{}`$. The value depends on the form of the extrapolation to $`q^2=0`$. We choose to extrapolate with a form factor quadratic in $`q^2`$, which gives $`BR(K^+\pi ^+e^+e^{})=(2.99\pm 0.06)\times 10^7`$ (see Table 1, column 2 of Ref. ). From this, one extracts the invariant amplitude defined by

$`_{K\pi ee}=A(q^2)(p_K+p_\pi )^\mu \overline{u}_e\gamma _\mu v_{\overline{e}}`$ (1)
giving

$`|A^{\mathrm{exp}}(0)|=(4.00\pm 0.18)\times 10^9\mathrm{GeV}^2,`$ (2)
found from the 3-body phase-space integral for the rate

$`\mathrm{\Gamma }(K^+\pi ^+e^+e^{})=0.2020{\displaystyle \frac{m_K^5A(0)^2}{3(4\pi )^3}}=(1.59\pm 0.03)\times 10^{23}\mathrm{GeV}.`$ (3)
As shown in , a quadratic fit provides a better agreement than a linear one for the same data. We will return to the $`q^2`$ dependence later.

The long-distance (LD) chiral low-energy model, used to calculate the bremsstrahlung graphs of Fig. 1(a) and (b), predicts

$`A_{\mathrm{LD}}=e^2\left|{\displaystyle \frac{\pi ^+H_WK^+}{m_{K^+}^2m_{\pi ^+}^2}}\right|\left|{\displaystyle \frac{F_{\pi ^+}(q^2)F_{K^+}(q^2)}{q^2}}\right|,`$ (4)
where $`F_{\pi ^+}(q^2)`$ and $`F_{K^+}(q^2)`$ are the pion and kaon electromagnetic form factors. At $`q^2=0`$, Eq. (4) becomes

$`A_{\mathrm{LD}}(0)=e^2\left|{\displaystyle \frac{\pi ^+H_WK^+}{m_{K^+}^2m_{\pi ^+}^2}}\right|\left|{\displaystyle \frac{dF_{\pi ^+}}{dq^2}}{\displaystyle \frac{dF_{K^+}}{dq^2}}\right|_{q^2=0}.`$ (5)
We now present the evaluation of the matrix element in Eq. (5), describing how we obtain $`3.9\times 10^9\mathrm{GeV}^2`$ for $`A_{\mathrm{LD}}(0)`$, reasonably close to the experimental result in Eq. (2). The matrix element $`\pi ^+H_WK^+`$ is well established; in Ref.  its value was deduced theoretically and was confirmed by comparison with values from ten measured kaon decays, $`K_S2\pi ^0`$, $`K3\pi `$, $`K_{L,S}2\gamma `$ and $`K_L\pi ^02\gamma `$, all of which are consistent with $`\pi ^+H_WK^+`$ = $`\pi ^0H_WK_L3.5\times 10^8\mathrm{GeV}^2`$. Specifically the $`K_S2\pi ^0`$ rate $`\mathrm{\Gamma }`$ gives

$$2\pi ^0H_WK_S=m_K\sqrt{16\pi \mathrm{\Gamma }/p}=(37.1\pm 0.2)\times 10^8\mathrm{GeV},$$
(6)
where $`p`$ is the three-momentum of a $`\pi ^0`$ in the $`K_S`$ rest frame. Using partially conserved axial currents (PCAC), this predicts for the pion decay constant $`f_\pi 93`$ MeV the LD scale

$$\pi ^+H_WK^+f_\pi 2\pi ^0H_WK_S3.5\times 10^8\mathrm{GeV}^2.$$
(7)
For $`F_{\pi ^+}(q^2)`$, there are many experimental measurements  for both positive and negative $`q^2`$, including our region $`0<q^2<0.125\mathrm{GeV}^2`$. The experimental data close to the region of relevance here are shown in Fig. 2. As expected from Vector Meson Dominance (VMD), these data are well described by a $`\rho `$ pole; the curve in Fig. 2 is the $`\rho `$-pole, i.e. VMD prediction,

$$F_{\pi ^+}(q^2)=\left(1\frac{q^2}{m_\rho ^2}\right)^1\frac{dF_{\pi ^+}}{dq^2}|_{q^2=0}=\frac{1}{m_{\rho ^0}^2}=1.69\mathrm{GeV}^2.$$
(8)
The pion charge radius $`r_\pi \sqrt{<r_\pi ^2>}`$, where $`r^2=6dF(q^2)/dq^2|_{q^2=0}`$, which in VMD equals $`0.628`$ fm, is in agreement with the experimental value  $`r_\pi =(0.63\pm 0.01)`$ fm. Since both $`F_\pi (q^2)`$ and $`r_\pi `$ agree well with the data, the $`\rho `$-pole expression can either be regarded as an input to our calculation or as an interpolating function between the experimental points. Note that up to a small G–parity violation, $`\omega `$ and $`\varphi `$ do not contribute to $`F_{\pi ^+}`$.

For the kaon case, no such comparable data exist, so we need a model for $`F_{K^+}(q^2)`$. In view of the success of the VMD picture for $`F_{\pi ^+}(q^2)`$, we use the same model for $`F_{K^+}(q^2)`$, resulting in the $`\rho `$, $`\omega `$ and $`\varphi `$ pole structure

$`F_{K^+}(q^2)`$ $`=`$ $`{\displaystyle \frac{N}{2}}\left({\displaystyle \frac{g_{\rho ee}}{m_{\rho ^0}^2q^2}}+{\displaystyle \frac{g_{\omega ee}}{m_\omega ^2q^2}}+\sqrt{2}{\displaystyle \frac{g_{\varphi ee}}{m_\varphi ^2q^2}}\right)`$ (9)
$`{\displaystyle \frac{dF_{K^+}}{dq^2}}|_{q^2=0}={\displaystyle \frac{n_\rho }{m_{\rho ^0}^2}}+{\displaystyle \frac{n_\omega }{m_\omega ^2}}+{\displaystyle \frac{n_\varphi }{m_\varphi ^2}}=1.42\mathrm{GeV}^2,`$
where $`g_{\rho ee}=5.03`$, $`g_{\omega ee}=17.06`$, $`g_{\varphi ee}=13.24`$ (derived from the $`e^+e^{}`$ decay widths), and the $`\rho ^0K^+K^{},\omega K^+K^{},\varphi K^+K^{}`$ $`\mathrm{SU}(3)`$ coefficients are $`\frac{1}{2}`$, $`\frac{1}{2}`$, $`\frac{1}{\sqrt{2}}`$, respectively. We deduce the $`\rho KK`$ coupling constant from $`g_{\rho \pi \pi }`$ (obtained using the $`\rho `$ width) via SU(3). We link the normalization $`F_{K^+}(0)=1`$ to $`N=0.037\mathrm{GeV}^2`$ in Eq. (9) leading to $`6n_{\rho ,\omega ,\varphi }=0.94,3.08,1.99`$. We then find that the charged $`K`$ radius is $`r_{K^+}=0.557\mathrm{fm}`$, in very good agreement with the data . The above lead to the prediction $`r_\pi ^2r_K^2=6\left(dF_{\pi ^+}/dq^2dF_{K^+}/dq^2\right)|_{q^2=0}=0.063\mathrm{fm}^2`$, consistent with the experimental value  $`(0.100\pm 0.045)\mathrm{fm}^2`$. Substituting into Eq. (5) the above values for the derivatives of $`F(q^2)`$ at $`q^2=0`$ (see Eqs. 8 and 9), together with Eq. (7) predict

$$\left|A_{\mathrm{LD}}^{\mathrm{VMD}}(0)\right|=3.9\times 10^9\mathrm{GeV}^2.$$
(10)
Folding in the $`\rho `$ width into the form factors and their derivatives has a negligible effect for the range of $`q^2`$ relevant to us here, i.e. $`q^2`$ between $`0`$ and $`0.125\mathrm{GeV}^2`$.

Our assumptions for $`F_{\pi ^+}(q^2)`$ and $`F_{K^+}(q^2)`$ are consistent with a specific model, which has the VMD and $`r_{\pi ,K}`$ structure built in, namely the quark-level linear $`\sigma `$ model (L$`\sigma `$M) as discussed in Refs. . The pion electromagnetic form factor, obtained from a $`u,d`$ quark triangle graph for a $`\pi ^+`$ probed by a photon (with pseudoscalar $`\pi qq`$ coupling $`g\gamma _5`$), predicts

$$F_{\pi ^+}(q^2)=i4N_cg^2_0^1𝑑x\frac{d^4p}{(2\pi )^4}\frac{1}{[p^2m_q^2+x(1x)q^2]^2},$$
(11)
nonperturbatively normalized to $`F_{\pi ^+}(0)=1`$ via the quark-loop version of the pion decay constant $`f_\pi `$, combined with the quark-level Goldberger-Treiman relation $`f_\pi g=m_q`$. To extract the $`q^2`$ dependence from Eq. (11) while integrating out the quark momentum $`p`$, one can differentiate Eq. (11), using $`r_\pi ^2(q^2)=6dF_{\pi ^+}(q^2)/dq^2`$ for $`g^2=(2\pi )^2/N_c`$ and $`r_{\pi ^+}=1/m_q0.63`$ fm :

$$r_{\pi ^+}^2(q^2)=6_0^1𝑑xx(1x)[m_q^2x(1x)q^2]^1.$$
(12)
Performing the above integration analytically and then making a Taylor series expansion keeping only the leading terms with $`yq^2/m_q^2`$ , one finds

$$m_q^2r_{\pi ^+}^2(q^2)=1+\frac{y}{5}+\frac{y^2}{23\frac{1}{3}}+\frac{y^3}{105}+....$$
(13)
Since $`m_qm_N/3`$, so $`m_\rho ^2/m_q^26`$, one can approximately express Eq. (13) in VMD language for the form factor itself in the low $`q^2`$ region as in Eq. (8). Lastly, the coupling $`g_{\rho \pi \pi }=6.04`$ is taken from the measured $`\rho `$ width, with the L$`\sigma `$M correction to VMD predicting  $`g_{\rho \pi \pi }/g_{\rho ee}=6/5`$, very close to the data ratio $`6.04/5.03`$. Then from $`\mathrm{SU}(3)`$ one gets $`g_{\rho K^+K^{}}=3.02`$. A test of this approach is to extract $`g_{\varphi K^+K^{}}=4.58`$ from the measured partial width of $`\varphi K^+K^{}`$ , from which we find, using $`\mathrm{SU}(3)`$ that $`g_{\rho K^+K^{}}=3.24`$. This value agrees to better than $`10\%`$ with its value $`3.02`$ above.

In the work of Ref. , the short-distance (SD) term was estimated to cause about a 20% reduction in the LD term. But more recent work by Dib, Dunietz and Gilman , based on a much heavier top-quark mass and also including QCD corrections, suggests that the SD term is much smaller. We therefore neglect the SD term of fig. 1(c) in this work.

Our predicted scale, then, is the LD amplitude of Eq. (10), which agrees with the data, Eq. (2), to about 5%.

III. FORM FACTOR

To extract the form factor for $`K^+\pi ^+e^+e^{}`$ from Eq. (4), we need the $`q^2`$ dependence of the electromagnetic form factors, $`F_{\pi ^+}(q^2)`$ and $`F_{K^+}(q^2)`$. In applying Eq. (4), we take $`F_{\pi ^+}(q^2)`$ from the data in Fig. 2, parameterised using the $`\rho `$ pole (Eq. (8)) and for $`F_{K^+}(q^2)`$, we use the equivalent expression, Eq. (9). We emphasise that we do not rely heavily on VMD for $`F_{\pi ^+}(q^2)`$; we use essentially the empirical values. The main reason to show the VMD fit to $`F_{\pi ^+}(q^2)`$ is to establish that this theory fits the data well, so it is reasonable to use it to obtain $`F_{K^+}(q^2)`$, to apply in Eq. (4).

With $`F_{\pi ^+}(q^2)`$, $`F_{K^+}(q^2)`$ and $`\pi ^+H_WK^+`$ taken as above, the LD part, Eq. (4), is completely determined. The resulting form factor, $`A_{\mathrm{LD}}`$, rises with $`q^2`$, in agreement with the data, but gives only about $`30\%`$ of the observed rise in the form factor, $`F^2`$, between $`q^2=0.03`$ and 0.12 $`\mathrm{GeV}^2`$. This is shown by the dashed line in Fig. 3.

Another possible contribution has been discussed by Ecker et al.  and by D’Ambrosio et al. , namely the charged pion loop term, Fig. 1(d). This term is derived in Refs.  and  using dimensional regularisation. We use here the expression for $`W^{\pi \pi }`$ from Ref.  without any additional polynomial. Also, we do not include terms from Refs.  and  other than the loop term; these contributions are calculated explicitly in our amplitude $`A_{\mathrm{LD}}`$. Evaluation of the pion loop term requires a knowledge of the $`K^+\pi ^+\pi ^+\pi ^{}`$ amplitude. In Refs.  and  this is taken from experiment. Ref.  shows that a current algebra-PCAC approach for $`K^+\pi ^+\pi ^+\pi ^{}`$ agrees with $`K_{3\pi }`$ data within 5%, so the $`q^2`$-dependent part of the loop in Refs.  and  is, in fact, compatible with the methods used here. The relative sign of the pion loop term and the dominant part of the amplitude, which may be considered as a parameter, was already established by the experiment .

Adding the amplitude, $`A_{\pi \mathrm{loop}}`$, from the pion loop term to $`A_{\mathrm{LD}}`$ gives the form factor shown as the solid line in Fig. 3. This agrees quite well with the data. The pion loop term makes a negligible contribution to the amplitude at $`q^2=0`$, so it does not disturb the agreement with experiment of the scale, $`A(0)`$, discussed in Sec. II and resulting in Eq. (10).

IV. CONCLUSIONS

In summary, we have shown that the recent Brookhaven data for $`K^+\pi ^+e^+e^{}`$ are in qualitative agreement, both in the $`q^2=0`$ scale and the form factor, with a calculation in which two processes are included. The experimental amplitude, $`|A^{\mathrm{exp}}(0)|=(4.00\pm 0.18)\times 10^9\mathrm{GeV}^2`$ (Eq. (2)), is about $`5\%`$ around the simple VMD bremsstrahlung prediction of $`\left|A_{\mathrm{LD}}(0)\right|=3.9\times 10^9\mathrm{GeV}^2`$ in Eq. (10), having dropped the SD term according to Ref. .

This agreement is as good as can be expected from such a simple model. In particular folding in a $`K^{}`$ pole through $`K^+\pi ^+K^0\pi ^+e^+e^{}`$ and a small, about $`10\%`$ SD contribution to the amplitude , may modify our prediction slightly. Note that as more terms are added, one increases the danger of double counting. Each of the graphs discussed here has been calculated before, though all have never been taken together in one calculation. This level of agreement is reasonable in view of the simplicity and combination of the models used.

There are no free parameters in our treatment, except for the relative sign between the pion loop and the rest of the amplitude and a relative sign between the couplings in Eq. (9); all quantities are taken from experiment, or from models that are known to fit other low-energy kaon or pion processes well. We do not attempt to estimate theoretical errors on our calculation, since these are to some extent dependent on one’s point of view. For example, the pion electromagnetic form factor, shown in Fig. 2, could equally well be regarded as taken from experiment, with its associated errors, or from the VMD prediction which is essentially free from error within the model. Similarly, the value of $`\pi ^+H_WK^+`$ is derived from theoretical arguments, but is confirmed by the fit to no less than ten other kaon decays . Either of these sources could be regarded as the origin of our value for $`\pi ^+H_WK^+`$ and hence as the source of error in it.

Concerning the recent paper  on $`K^+\pi ^+ee`$, we agree with its recommended form factor structure (as in our Eq. (4); compare with its Eqs. (2), (7) and (8)). However we disagree about the scale of Eqs. (7) and (8) in  which seems more akin to the small $`\mathrm{\Delta }I=3/2`$ $`K\pi \pi `$ tree level transition rather than to the much larger $`K^0\pi \pi `$ $`\mathrm{\Delta }I=1/2`$ LD transition, as in our Eqs. (6) and (7).

As for Chiral Perturbation Theory (ChPT) applied to $`K^+\pi ^+ee`$ , the required $`𝒪(p^6)`$ expression contains about 100 parameters. Moreover, the underlying $`\mathrm{\Delta }I=1/2`$ LD scale as in our Eq. (7) disappears for $`N_C\mathrm{}`$ .

ACKNOWLEDGEMENTS

We are grateful to T. Goldman, G.J. Stephenson Jr., P. Singer and F. Krauss for useful discussions. We acknowledge support from the US DOE. GE would like to thank the US-Israeli Binational Science Foundation and the Fund for the Promotion of Research at the Technion for partial support. MDS appreciates partial support from the Technion.


text:
# Simulated ecology-driven sympatric speciation

## 1 Introduction

The theory of evolution studies processes that are extremely complex and whose characteristic times usually extend over many human life spans. Because of this inherent complexity, these studies do not lend themselves readily to the establishment of discrete, falsifiable hypotheses. For this reason, observations and experimental results emanating from evolutionary studies often elicit broadly differing interpretations . A large number of mathematical models have been applied to a variety of evolutionary problems, as well as to the related fields of population dynamics and evolutionary ecology, composing a unified approach under the name of evolutionary population ecology . The success of many of these models, together with the above-mentioned difficulties in the interpretation of observational data, has triggered the development of an integrated computational effort for the study of these problems. Computer simulations of the behavior of virtual populations, immersed in well-defined habitats and evolving under a representation of the dynamics of evolution, have generated a wealth of reliable data with which to test our models and the assumptions they rely on.

Our approach in this paper is to use a computational representation to study the genetical patterns generated by the order-disorder conflict between selective pressure and mutation accumulation in the presence of an environment that favors particular phenotype configurations. Our model is a variation on the popular Penna model , a bit-string representation of genetic dynamics for age-structured populations that has been successfully used to address a variety of problems in the field. It provides a simple metaphor for evolutionary dynamics in terms of the mutation accumulation theory. In essence, each individual is represented by a double strand of inherited genetic material, subject to the addition of harmful mutations at birth. Its life span is determined by the position of the active deleterious mutations in the genome, while non-mutated genes have no effect on the individual adaptation to its environment.

To this underlying age-structured genome we add a non-structured portion, which can be seen as representing genes that do not have an age-dependent action. The configuration of this new portion of the genome is translated into a phenotype which determines, through its match with an environment-dependent ideal, the probability of survival when the number of harmful mutations is not critical. Models with similar characteristics have been studied recently, but under a different perspective . In our case, the phenotype is translated into a fitness function merely by counting the number $`n`$ of positions in the bit-string that are set to a given value, say $`1`$, and the individual fitness is a function of this number $`f(n)`$. This strategy simulates a single phenotypical trait which is encoded in a number of genes. We then study the distribution of phenotypes of the equilibrium population, identifying its main patterns.

Beyond the theoretical interest of unveiling the features generated by evolutionary dynamics on the phenotype of the population, our model also allows us to address a much more challenging problem, namely that of sympatric speciation. Speciation involves the division of a species on an adaptive peak, so that each part moves onto a new adaptive peak without either one going against the upward force of natural selection. This process is readily envisioned if a species becomes subdivided by a physical barrier, whereby each part experiences different mutations, population fluctuations and selective forces, in what is called the allopatric model of speciation. In contrast, conceiving the division of a single population and radiation onto separate peaks without geographical isolation, in what is called sympatric speciation, is intuitively more difficult . Through which mechanism can a single population of interbreeding organisms be converted into two reproductively isolated segments in the absence of spatial barriers or hindrances to gene exchange?

As an answer to this question, a general mechanism for sympatric speciation has been modeled with a simple one-locus two-allele dynamics that shows how a stable polymorphism can exist in a heterogeneous environment with two niches, even when the adults form a single randomly-mating population . Increasingly realistic models that explicitly incorporate density-dependent competition for resources among phenotypes have since appeared in the literature . The present paper is another link in this chain. We show that a stable polymorphism is generated when the environment switches from a one to a two-niche configuration in a model with age-structured genomes and competition among phenotypes, as a necessary result of the dynamics of evolution.

The establishment of a stable polymorphism is the first step towards sympatric speciation. The second is the subsequent evolution of reproductive isolation between the populations in the two niches. The key question here is whether a single gene difference could produce selective coefficients large enough to drive reproductive isolation. It is generically thought that this isolation requires a genetic association between traits conferring adaptation to a niche and traits causing mating preferences . Alternatively, one can consider assortative mating of phenotypes; since a gene for assortative mating can increase to fixation sympatrically, speciation could be the outcome of sexual selection alone. This is what our model shows: once we introduce a gene which describes the mating preferences of an individual, it is lead by the dynamics to a fixation pattern in which alike phenotypes mate exclusively.

As a metaphor to the results obtained in our simulations, we compare them qualitatively with the observed dynamics of the Galapagos finches. Growing evidence for the onset of stable polymorphism in nature comes from years of observations of seasonal morphological variations in the population of ground finches in the Galapagos Islands. The assumptions that selection, mediated through rainfall and its effects on food supply, can have a dramatic impact on finch phenotypes, and that much of the morphological variation in these birds is genetically inherited are reasonably well established since the field work of P.R. Grant and his collaborators .

## 2 The Penna model with phenotype selection

In this section we will briefly describe the main characteristics of the standard age-structured Penna model, together with the features we added to it to represent phenotypic selection. For a recent review on the applications of the model we direct the reader to Ref.  and references therein.

For sexual populations, the genome of each organism is represented by two computer words. In each word, a bit set to one at a locus corresponds to a mutated deleterious allele; a “perfect” strand would be composed solely of zeros. The effect of this mutation may be felt by the individual at all ages equal to or above the numerical order of that locus in the word. As an example, a bit set to $`1`$ at the second position of one of the bit-strings means that a harmful effect may become present in the life history of the organism to which it corresponds after it has lived for two time periods. The diploid character of the genome provides for a coverage against the effectiveness of these harmful mutations. Thus, a mutation in a position of one of the strands is felt as harmful only if this is a homozygote locus or if the mutated allele is dominant. In a homozygote locus, a harmful mutation must be present in both strings at the same position to be effective. The concept of dominance, on the other hand, relates to loci in the genome in which a harmful mutation in just one strand is enough to make it affect the organism’s life. The life span of an individual is controlled by the number of effective mutations active at any instant in time. This number must be smaller than a specified threshold to keep the individual alive; it dies as soon as this limit is reached.

An individual may also die because of intra-specific competition for the finite resources of the environment or because of the action of predators. In the standard model, these constraints are taken care of by the so-called Verhulst factor. This is a logistic-type term that introduces a mean-field random death probability, independent of the quality of the genome, and in the simulations its main purpose is to limit the size of the population. Its usual expression is $`V(P)=P/C`$, where $`P`$ is the total population at some time step and $`C`$ is a parameter of the simulation, traditionally called the carrying capacity of the environment. This name, inherited from simple logistic models , is somewhat misleading in the context of more sophisticated models; in those simple models the equilibrium population ends up to be $`C`$, thus justifying its name. This is not the case for models that include a more detailed microscopic description of the evolutionary interactions, where the equilibrium population is a complex function of many different parameters.

Since really random deaths in nature can hardly play any significant role in population dynamics, this concept has already been criticized in the literature, and a couple of alternate implementations analyzed . Rather than random, one should expect that deaths caused either by predation or because of intra-specific competition should also have selective value. The probability of their occurrence should be dependent on the fitness of the individual to the environment, and this fitness a function of the match between the individual’s phenotypic expression of genetically acquired traits and a phenotype ideally adapted to its habitat. The problems that such a microscopic description of selective dynamics entails are extremely complex and are far from being solved. The present paper reports on an attempt to address this issue in the context of an age-structured model. In addition to the usual double strand of genes with age-dependent expression, we introduce an extra pair for each individual in the population. In this computer representation, each genome may be considered as composed of $`64`$ loci - in a $`32`$-bit word machine - of which half have age-dependent expression. The dynamics of reproduction and mutations, to be described in what follows, is the same for both the age-structured and the new strings. Meiosis, phenotypic expression controlled by homozygose or dominance, and the introduction of random mutations at birth affect equally all loci. The new portion of the phenotype, however, represents some individual trait, such as the size of the beak in a ground finch, for example, that may have selective value. We chose a multi loci representation of a single trait for simplicity, and claim that there is no loss of generality in this choice for our present purposes. The selective value of this phenotype is expressed by a fitness function $`f(n)`$, chosen to depend on a single variable $`n`$: the number of $`1`$’s expressed in the non-structured portion of the phenotype, either because of homozygose ($`11`$ alleles) or the presence of $`10`$ ($`01`$) configuration in a locus where the $`1`$ allele is dominant. This fitness function encapsulates the selective value of a particular phenotype and restricts the multi-dimensional space of the interaction between the individual and its environment to one single dimension. In this discretized mathematical model, the probability of death by intra-specific competition at each time step is now given by $`V(P)/f(n)`$. The interaction with other individuals is still mediated by the mean-field $`V(P)`$, but its effect is no longer uniform. It depends on the particular genetically acquired configuration of each individual and, although stochastic, escapes from the biologically unmotivated randomness mentioned above.

Reproduction is modeled by the introduction of new genomes in the population. Each female becomes reproductive after having reached a minimum age, after which she generates a fixed number of offspring at the completion of each period of life. The meiotic cycle is represented by the generation of a single-stranded cell out of the diploid genome, now composed of two functionally different portions. To do so, each string of the parent genome is cut at a randomly selected position, the same for both strings and independently for each functional part, and the left part of one is combined with the right part of the other, thus generating two new combinations of the original genes - see Figure 1. The selection of one of these completes the formation of the haploid gamete coming from the mother. For mating, a male is randomly selected in the population and undergoes the same meiotic cycle, generating a second haploid gamete out of his genome. The two gametes, one from each parent, are now combined to form the genome of the offspring. Each of its strands was formed out of a different set of genes. The next stage of the reproduction process is the introduction of $`M_a`$ independent mutations in the newly generated genetic strands of the age-structured part of the genome. In this kind of model it is normal to consider only the possibility of harmful mutations, because of their overwhelming majority in nature, and this is what is done for the age-structured portion. For the non-structured portion, mutations are added independently of their harmful or beneficial effect, with some probability $`M_p`$ per string. The gender of the newborn is then randomly selected, with equal probability for each sex. Figure 1 is a schematic representation of the entire reproductive process.

The passage of time is represented by the reading of a new locus in the age-structured portion of the genome of each individual in the population, and the increase of its age by one. After having accounted for the selection pressure of a limiting number of effective harmful mutations and the stochastic action of the fitness function, females that have reached the minimum age for reproduction generate a fixed number of offspring. The simulation runs for a pre-specified number of time steps, at the end of which averages are taken over the population. Typically, measures are taken for the age structure of the population, such as the number of individuals and probability of survival and death by genetic causes for each age group, as well as for the genetic composition distribution. We are particularly interested in the equilibrium distribution of the phenotypes in the population, also averaged at the end of each run.

## 3 Simulation results

We will relate our results qualitatively to the observations cited in the introduction made on the evolutionary studies of the Galapagos finches. Thus, we will have in mind a picture of our population as composed of ground finches, feeding primarily on seeds. The phenotype trait represented by our model describes the size and strength of their beaks, and the fitness measures their adequacy to the kinds and sizes of seeds available in the ecology.

Let us start with results concerning the equilibrium configuration of the single trait that characterizes the fitness of an individual’s phenotype. For this study, the fitness function is

$$f(n)=1\frac{32n}{A};$$
(1)
$`n`$ is the number of bits set to $`1`$ in an individual’s phenotype, and $`A`$ is a boost factor that can be used to control the intensity of the selective pressure; this intensity is a decreasing function of $`A`$. This fitness function expresses an environment where the ideal phenotype would be composed entirely of bits $`1`$. For instance, it can represent an environment in which the availability of edible seeds is a decreasing function of their size, while the number $`n`$ is a measure of the size of the beak of an individual. In this ecology, a well-fit individual would have a large beak, allowing him to take advantage of the distribution of seeds size. The simulation starts with an ill-fit population, in which all the individuals have the smallest beak possible; in our representation, their phenotypes are homozygotes $`00`$ in all loci. Alternatively, this can be expressed as a distribution of beak sizes delta-like peaked at $`n=0`$. Figure 2 shows the equilibrium distribution reached after $`\mathrm{50\; 000}`$ Monte Carlo steps, for a weak selective pressure, corresponding to a broad distribution of seed sizes and, in our representation, to a large value of the boost $`A=128`$. We plot the distribution of beak sizes for various values for the number of loci $`d`$ in which the well-fit ($`1`$) allele is dominant in the phenotype. In all cases, the distribution shows a bell-shaped Gaussian-like aspect in which the peak position and amplitude increase with $`d`$, while the width shows a decreasing pattern. This result is easy to understand: because of the availability of large seeds, represented here by a fitness function that increases with $`n`$, the population is driven to a pattern of larger beaks. The combined effects of mutation pressure and heterozygose avoids its collapse into a delta-like distribution from which it evolved; the width of the distribution increases with the mutation rate. In a strongly selective ecology, represented by a smaller value of $`A`$, the pattern is essentially the same. The peaks are shifted to larger values of $`n`$ and the widths decrease.

The shift in the distributions caused by stronger selection is shown in Figure 3. We compare results obtained for a fixed value of $`d=0`$, where the allele $`1`$ is recessive in all loci. When there is no activation of phenotype selection, which can be translated into a uniform distribution of seed sizes, the broad pattern of beak sizes is a result of random genetic drift. If $`d=16`$, the equilibrium distribution peaks at $`n=16`$. When phenotype selection is activated, the equilibrium distribution tries to follow the availability of seed sizes and its peak moves to larger values of $`n`$, while its width shrinks. It is worthwhile to remark that even in a very strong selective ecology ($`A=33`$), the distribution still retains some broadness, representing a diversity of morphologies also observed in real populations.

We examine next the pattern of homo/heterozygose generated in the phenotype of the population. Figure 4 was obtained for a medium range selective pressure $`A=48`$. The circles indicate what is the dominant allele in each locus. One can immediately identify that the selection pressure pushes the loci in which the allele $`0`$ is dominant to homozygose $`11`$, while on sites where $`1`$ is dominant hetero and homozygose are in fact indistinguishable.

The establishment of a stable polymorphism is the theme of Figure 5. In the simulation, the initial population was again a small sized beak one, now immersed in an ecology with a broad distribution of edible seed sizes available, peaking at middle-sized seeds. The fitness function that represents this situation is

$`f(n)=1{\displaystyle \frac{16n}{A}},n<16`$ (2)
$`=1{\displaystyle \frac{n16}{A}},n16`$ (3)
The population evolves for $`\mathrm{20\; 000}`$ Monte Carlo steps, when a snapshot of the phenotype distribution is taken. In agreement with the result shown in Figure 2, this distribution is again bell-shaped, with its peak now located at $`n=16`$, corresponding to middle-sized beaks. Because mutations can both increase or decrease the beak size, as opposed to the case where the peak of the fitness function was at the largest possible extreme, and because the number of loci where each allele is dominant is the same, corresponding to a value $`d=16`$ for the dominance parameter, there is no inherent bias to the equilibrium distribution; its peak can sit at the same position as the one for the fitness function, as Figure 5 shows. After $`\mathrm{20\; 000}`$ time steps, there is a sudden change in the pattern of seed availability, perhaps because of a variation in the rainfall regime caused by some global climate oscillation, whose effect is to decrease the availability of middle-sized seeds. The fitness function that expresses this new pattern is, for instance,

$`f^{}(n)=1{\displaystyle \frac{A16+n}{A}},n<16`$ (4)
$`=1{\displaystyle \frac{An+16}{A}},n16`$ (5)
The evolution of the population resumes with this new fitness function, and with intra-specific competition within each niche. This is represented in the model by letting the mean-field $`V(P)`$ take into account only the individuals that feed on a particular kind of seed. Thus, finches with small (large) beaks, represented by having $`n<(>)16`$, and feeding on small (large) seeds, will compete only against those feeding on the same niche. To be specific, the simulation now computes at each time step $`P_m=\mathrm{\Sigma }_{population}\theta (16n)`$ and $`P_M=\mathrm{\Sigma }_{population}\theta (n16)`$, where $`\theta (k)`$ is the step function with $`\theta (0)=1/2`$. To determine the probability of death by intra-specific competition for an individual, its beak size, which is proportional to the number $`n`$ of $`1`$-bits set in its phenotype, is used to decide in which niche it is competing. If $`n<(>)16`$, its death probability will be $`V(P_m)/f^{}(n)`$ ($`V(P_M)/f^{}(n)`$). These elements alone force evolution to give rise to a polymorphism, shown in Figure 5 as a resulting $`2`$-peaked equilibrium distribution of beak sizes. This polymorphism is reversible: if, in a subsequent time step, the pattern of availability of edible seeds reverts to its original configuration, so does also the distribution of beak sizes. This is in complete agreement with the field observations of stable polymorphism in the Galapagos ground finches .

If such a stable polymorphism is not challenged by new climate variations for some time, it can induce speciation and generate two reproductively isolated segments, even without geographical barriers. To simulate this effect, we introduced in our model a gene that determines mating selection pattern. The initial population has a panmictic behavior, expressing no preferences related to the phenotype of the partner. When a female is ready to mate, she chooses a partner according to the expression of this gene. A randomly selected male in the population, to be accepted as a partner, has to either feed on the same niche, in which case the mating selection gene becomes irrelevant, or, if he feeds on a different niche, both parents have to be non-selective in their mating preferences for reproduction to occur. The offspring inherit the mating preferences of either the mother or the father, randomly selected at birth, and this gene can also suffer a mutation in either direction with some probability. The resulting distribution of this gene among the population is shown in table I. After the double peaked polymorphism is established, and if the distribution of edible seeds does not change for some time, reproductive isolation sets in.

## 4 Conclusions

The addition of a multi-locus non-structured part to the age-structured genome of the Penna model, subjected to the same representation of evolution dynamics, allows a substantial broadening of its scope. The role of genes with age-independent effects can now be studied through the patterns that they generate in evolutionary populations. A particularly interesting side benefit is to get rid of the unmotivated random effect of the Verhulst dagger and substitute it by inter-specific competition with outcome based on the genetic load of the individuals. With the help of this broadened model we studied the distribution of a single phenotypical trait among the population when it evolves on a single-peaked linear fitness landscape. The simulations exhibited its Gaussian-like character with mean and width that depend in a clear way on the various parameters of the modified model. The distribution of homo and heterozygote loci of the non-structured genome over the population also shows patterns that are easy to explain and to relate to the parameters. The present study opens the path for a variety of new applications of age-structured microscopic models as it establishes their basic properties.

As a hint of the scope of these applications, we undertook the simulation of the effects of a varying pattern of food supply on the phenotype distribution . Our motivation were the reports from field observations of the variability of the beak morphology of the ground finch in the Galapagos Islands as a response to changes in seed availability due to seasonal oscillations in rainfall regime. The results of the simulations show a striking parallelism with those observations, in that a polymorphism is generated as soon as the fitness landscape switches to a double-peaked pattern. Since modern theories on speciation seem to indicate that polymorphism is a necessary first step, we proceeded by introducing a gene that controls mating preferences. Mutation and evolution under a continued double-peaked fitness landscape were enough to generate reproductive isolation and speciation as a result. Speciation without geographic isolation has proven to be a difficult concept in evolution theory so far, but our simulations support a scenario in which its onset arises as a necessary evolutionary consequence of a continued double-peaked fitness landscape.

## Acknowledgments

We thank D. Stauffer for a critical reading of the manuscript. Research by J.S.S.M. was supported as a Visiting Fellow by CIRES, University of Colorado at Boulder. SMO thanks the Brazilian agencies CNPq, CAPES and FAPERJ for partial financial support.


text:
# Asymptotically free theories based on discrete subgroups Work supported in part by Schweizerischer Nationalfonds Talk presented by F. Niedermayer

## 1 Introduction

Here we summarize some results of a recent paper  where discretized version of 2d O(3) sigma model and of 4d SU(N) (N=2,3) gauge models were studied near their critical points. Here we discuss only the 2d sigma model. For further results and a more extensive list of references we also refer to that paper.

We consider models with discretized spin variables taking values from the vertices of regular polyhedra. The symmetry group is accordingly reduced from O(3) to the corresponding discrete subgroup. At strong coupling (i.e. large temperature), the fluctuations are large, hence the effects of restricting the spins to a discrete set are expected to be small. However, at sufficiently small coupling the discrete system has to be frozen. Decreasing $`g`$ from the strong coupling regime one expects a second order phase transition at some value $`g_c`$. We investigate the quantum field theory obtained in the limit $`gg_c`$ and suggest that for sufficiently large subgroups (icosahedron and dodecahedron) the resulting theory is the same as the one given by the original O(3) model in the $`g0`$ limit.

The latter being asymptotically free (AF) this sounds surprising. Our conjecture that the discrete model is AF should be understood in the sense that the physical running coupling goes to zero as the corresponding momentum scale goes to infinity. Here we study only the question whether the discrete model is equivalent to the original O(3) model (which we can establish only numerically, within our numerical precision, of course). The title refers to the standard wisdom that the O(3) model is AF – a belief debated for long time by Patrascioiu and Seiler (for references see ). Here we do not study this question.

The discrete subgroups were introduced  in the SU(2) Yang-Mills theory as an approximation thought to be valid until the Wilson loops measured at the same bare coupling start to disagree. In contrast to this we compare the theories with discrete and continuous groups at the same (large) correlation length. That a discrete symmetry of the action can be enhanced to the continuous group at the criticality is known for the case of the 2d XY-model .

The investigations reported in ref.  and here, were inspired by the works of Patrascioiu and Seiler on the dodecahedron model  where they observed that it behaves as the O(3) model. We investigated the conjecture that the discrete and continuous models are in the same universality class to a high precision, down to O(0.1%).

We consider here discretized spin models when the spin vectors point to the vertices of a regular polyhedron embedded into the sphere $`S_2`$. If not stated otherwise, we consider the standard nearest neighbour (nn) action. The spin models with small subgroups of O(3) (tetrahedron, octahedron and cube, with 4,6 and 8 directions, respectively) are not equivalent to the original O(3) model. The tetrahedron is equivalent to a $`q=4`$ Potts model, while the cube to 3 independent Ising models. Besides the standard O(3) action we also compared some of the results to those for the O(3) fixed point action .

We have chosen physical quantities which can be measured with high precision: the finite size scaling function (FSSF) , the renormalized zero momentum 4-point function $`g_R`$ and the same coupling $`g_R(z)`$ in a finite physical volume. It is the latter quantity which can be measured with the highest precision.

## 2 The finite-size scaling function

We consider an $`L\times L`$ periodic box and measure the second moment correlation length $`\xi (L)`$. At the same inverse coupling $`\beta `$ we also measure the quantity $`\xi (2L)`$ on a $`2L\times 2L`$ lattice, and study the ratio $`\xi (2L)/\xi (L)`$ as a function of $`\xi (L)/L`$. The techique of finite size scaling (FSS) was used very effectively by Lüscher, Weisz and Wolff to obtain the running coupling in the O(3) model. Instead of the (exponential) correlation length defined in a strip $`L\times \mathrm{}`$ used in  we use the second moment correlation length defined in a square box by

$$\xi (L)=\frac{1}{2\mathrm{sin}(\pi /L)}\sqrt{\frac{G_2(0)}{G_2(k_0)}1},$$
(1)
where $`G_2(k)`$ is the 2-point spin correlation function in Fourier space and $`k_0=(2\pi /L,0)`$. The finite-size scaling function for this quantity has been measured by Caracciolo et al. in the same model for a large set of $`\xi (L)/L`$ values. We have chosen to measure around the value $`\xi (L)/L0.4`$.

Fig. 1 shows the FSSF for O(3) and different subgroups. Points with symbols of the same shape belong to the same subgroup at different values of $`L`$, with the largest symbol corresponding to the largest $`L`$. While the tetrahedron, octahedron and cube deviate strongly from the O(3) curve and move away with increasing $`L`$, the icosahedron and dodecahedron results are close to O(3), and at largest $`L`$ lie on the curve within the small statistical errors.

## 3 The renormalized zero momentum 4-point coupling $`g_\mathrm{R}`$

Define the quantity $`g_\mathrm{R}(z)`$ as

$$g_\mathrm{R}(z)=\left(\frac{L}{\xi (L)}\right)^2\left(1+\frac{2}{N}\frac{(𝐌^2)^2}{𝐌^2^2}\right),$$
(2)
where $`\xi (L)`$ is the second moment correlation length, eq. (1), $`z=L/\xi (L)`$, $`N=3`$ for O(3) and $`𝐌`$ is the magnetization,

$$M^a=\underset{x}{}S^a(x).$$
(3)
The coupling $`g_\mathrm{R}=g_\mathrm{R}(\mathrm{})`$ is defined as the infinite volume limit $`z\mathrm{}`$, where first the continuum limit $`\xi \mathrm{}`$ is taken for fixed $`z`$. For the O(3) $`\sigma `$-model $`g_\mathrm{R}`$ has been calculated using the form-factor bootstrap method: $`g_\mathrm{R}=6.770(17)`$ in agreement with the MC result $`g_\mathrm{R}=6.77(2)`$ .

We have measured $`g_\mathrm{R}(z)`$ for the icosahedron and dodecahedron models for $`z6`$ and extrapolated to infinite volume using the finite size formula of ref. . Figure 2 gives the deviation from the O(3) result as the function of $`\xi `$. For comparison, the coupling from the cubic group is $`g_\mathrm{R}^{\mathrm{cubic}}=1/3g_\mathrm{R}^{\mathrm{Ising}}`$ with $`g_\mathrm{R}^{\mathrm{Ising}}=14.6975(1)`$ from ref., i.e. for this case the deviation from O(3) is 1.87(2).

As Figures 1 and 2 show, the two largest non-Abelian subgroups behave qualitatively differently from the small subgroups. The results suggest that the icosahedron and dodecahedron models are in the same universality class as the O(3) non-linear $`\sigma `$-model.

## 4 The renormalized coupling in a finite volume

The quantity $`g_\mathrm{R}(z)`$ can be measured with much better accuracy for smaller z values as for large ones. In addition, one avoids the uncertainty by taking the infinite volume limit $`z\mathrm{}`$. We have chosen to measure in the vicinity of the arbitrary value $`z_0=2.32`$.

Note that both $`g_\mathrm{R}(z)`$ and $`z`$ are measured in the same run. Since they turn out to be strongly correlated the statistical fluctuations can be strongly reduced by considering the quantity $`g_\mathrm{R}(z)c(zz_0)`$ with an appropriate value of $`c`$. Interpolating in $`zz_0`$ gives $`g_\mathrm{R}(z_0)`$. By measuring this quantity we could compare results for different models to a high precision.

Table 1 summarizes the results for the three actions considered.

Figure 3 shows $`g_\mathrm{R}(z_0)`$ as the function of $`a/L`$ for the three actions considered.

The striking feature here is that even the standard O(3) action shows a linear, O($`a/L`$), cut-off effect, while the perturbative approach predicts an O($`(a/L)^2`$) behaviour (up to $`\mathrm{log}(L/a)`$ factors).

Below we present for the O(3) nn action fits with different functional forms, together with the corresponding value of $`\chi ^2/\mathrm{dof}`$:

$`3.1201(2)1.224(9){\displaystyle \frac{a}{L}},`$
$`3.1205(4)1.28(5){\displaystyle \frac{a}{L}}+1.1(1.1)\left({\displaystyle \frac{a}{L}}\right)^2,`$
$`3.1207(6)1.1(1){\displaystyle \frac{a}{L}}0.04(4){\displaystyle \frac{a}{L}}\mathrm{log}{\displaystyle \frac{L}{a}},`$
$`3.1168(3)+103(5)({\displaystyle \frac{a}{L}})^240.4(1.4)\left({\displaystyle \frac{a}{L}}\right)^2\mathrm{log}{\displaystyle \frac{L}{a}},`$
The corresponding values of $`\chi ^2/\mathrm{dof}`$ are 2.1, 2.2, 2.5 and 4.1, respectively. (The fits for the FP action and dodecahedron are given in .)

One concludes that an $`a/L`$ term is needed to describe the cut-off effects. (Note also the large coefficients in the $`a^2`$ fit!)

## 5 On the cut-off effects

For bosonic models in any order of perturbation theory the cut-off effects are given by O($`a^2`$) terms, up to logarithmic corrections . Although not proven beyond perturbation theory, this form of cut-off effects is generally assumed when one extrapolates to the continuum limit. Therefore an O($`a`$) is rather surprising. Note that already in ref.  such behaviour for $`g_\mathrm{R}`$ at $`z5`$ was preferred. However, due to larger errors this form was not compelling.

Because the Ansatz for the cut-off corrections can significantly modify the predictions in the continuum limit, and because a non-standard behaviour of the lattice artifacts can have other theoretical implications, it is important to investigate the related questions further. In particular, it will be interesting to study other non-perturbative quantities as the LWW coupling  in the O(3) sigma model. Obviously, the same questions in d=4 gauge theories are even more relevant.