text:
# A BeppoSAX Observation of KS 1731–260 in its Quiescent State: Constraints on the NS Magnetic Field

## 1 Introduction

The Galactic X–ray source KS 1731–260 was discovered in August 1989 with the imaging spectrometer TTM on the Mir–Kvant observatory (Sunyaev 1989; Sunyaev et al. 1990). During the $`15`$ days long Mir–Kvant observation the source intensity varied from 50 to 100 mCrab in the 2–27 keV band. The presence of three type–I bursts indicated that the compact object is an accreting neutron star (hereafter NS) and that the source distance is about 7 kpc (Muno et al. 2000). The factor of 10 intensity variations displayed by the source suggested that KS 1731–260 is a transient source (e.g. Sunyaev et al. 1990). There were numerous detections of the source at a level of 50–100 mCrab since its discovery in 1989, and the monitoring by the Rossi X-ray Timing Explorer (RXTE) All Sky Monitor (ASM) showed that the source was continuously active until February 2001. Therefore the source appeared to belong to a group of NS soft X-ray transients (hereafter SXTs) that display active states extending for a number of years (as opposed to outbursts decaying on timescales of a few months; see Campana et al. 1998b for a review).

The source spectrum resembles that of other SXTs in outburst; the TTM spectrum taken on August 16–31, 1989 could be fit to a thermal bremsstrahlung with $`kT5.7`$ keV while the spectrum obtained by SIGMA on March 14, 1991 (when the source was somewhat fainter, $`9\times 10^{36}`$ erg/s as opposed to $`1.5\times 10^{37}`$ erg/s) was well modelled by a power law with a photon index of $`2.9`$ extending up to 150 keV at least. The source detection during the ROSAT all-sky survey yielded a fairly accurate measurement of column density, $`N_H=(1.00\pm 0.19)\times 10^{22}`$ cm-2 (Barret, Motch, & Predehl 1998).

Based on an RXTE/PCA observation in July 1996, Smith, Morgan, & Bradt (1997) discovered a nearly coherent signal at $`523.92\pm 0.05`$ Hz (corresponding to a period of 1.91 ms) during a two-second interval close to the beginning of the decay of a type I X-ray burst. This signal, as well as similar signals observed during type I X-ray bursts from about ten low mass X-ray binaries, likely corresponds to the NS rotation frequency (or twice its value; for a review see Strohmayer, Swank, & Zhang 1998). An RXTE observation in August 1, 1996 (Wijnands & van der Klis 1997) revealed two simultaneous quasi periodic oscillations (QPO) in the persistent emission of KS 1731–260 at 898 Hz and 1159 Hz, the frequency separation of which ($`260.3\pm 9.6`$ Hz) is compatible with half the frequency of the nearly-coherent signal observed in the type–I bursts. In the beat frequency model for the kHz QPOs, this corresponds to a NS spinning at a period of $`3.8`$ ms, which might imply that only the first harmonic is detected during type–I bursts (although this case seems to be unlikely, see Muno et al. 2000).

RXTE/ASM data showed that the source entered a quiescent state in February 2001 (upper limit of 2 count/s corresponding to a luminosity of $`<3\times 10^{36}`$ ergs/s). A few months after (on 27 March 2001) KS 1731–260 was observed by Chandra: it was still in quiescence with a 0.5–10 keV bolometric luminosity of $`10^{33}`$ ergs/s (assuming a distance of 7 kpc) and an X-ray spectrum that is well described by a blackbody at a temperature of $`0.3`$ keV (Wijnands et al. 2001a). The same data were also analyzed by Rutledge et al. (2001); they fitted the spectrum using a hydrogen atmosphere model, obtaining an effective temperature of 0.12 keV, an emission area radius of $`10`$ km and a bolometric luminosity of $`2.7\times 10^{33}`$ ergs/s (assuming a distance of 8 kpc). We present here the results of a 200 ks BeppoSAX observation of KS 1731–260 taken in March 2001, just after the beginning of this quiescent state, in which the source was detected at a quiescence luminosity level of $`(7\pm 2)\times 10^{32}`$ ergs/s.

## 2 Observations and Data Analysis

KS 1731–260 was observed with the BeppoSAX Narrow Field Instruments (NFI) from 2001 March 2 23:28 UT to March 5 19:07 UT for 200 ks. The effective exposure time was $`91`$ ks. The NFIs consist of four co-aligned instruments, namely a Low Energy Concentrator Spectrometer, LECS, a thin window position sensitive proportional counter with extended low energy response in the band 0.1–10 keV (Parmar et al. 1997, due to UV contamination problems, the LECS was operated only at satellite night time, resulting in a reduced exposure time with respect to the other instruments); two (originally three) Medium Energy Concentrator Spectrometers, MECS, operating in the 1.3–10 keV band (Boella et al. 1997); a High Pressure Gas Scintillation Proportional Counter, HPGSPC, sensitive in the energy range 7–60 keV (Manzo et al. 1997); a Phoswich Detection System, PDS, sensitive in the energy range 13–200 keV (Frontera et al. 1997), consisting of four independent NaI(Tl)/CsI(Na) phoswich scintillation detectors.

Images from the two MECS were summed and the events with PHA energies between 1.3 and 10 keV were considered. We searched for sources in the MECS field of view (FOV) around the nominal position of KS 1731–260 (Barret et al. 1998). We found a source centered at the position RA: 17h 34m 09.4s, DEC: $`26^{}`$ $`06^{}`$ $`22^{\prime \prime }`$, equinox 2000, and extracted all the photons from a $`1.8^{}`$ radius circle (corresponding to 64% of the point spread function, PSF; Fiore, Guainazzi & Grandi 1999). The size of the extraction region was determined to optimize the signal to noise ratio. The accuracy and stability of the pointing during the BeppoSAX observation of KS 1731–26 has been checked and the satellite attitude has been verified within each orbit by telemetry reconstruction. In particular we have verified that the positions of the stars in the on board startrackers FOV remained constant during all the observation. In principle, the error on the BeppoSAX/MECS position is $`0.55^{}`$ ($`1\sigma `$ confidence level); however, given the possibility of confusion with other unresolved sources (see below), we will consider as a conservative estimation of the error on the source position the radius, $`1.8^{}`$, of the extraction region.

The recently reported source position based on a Chandra observation performed on March 27 2001 and corrected for an offset with respect to the 2MASS (Wijnands et al. 2001b) is RA: 17h 34m 13.47s, DEC: $`26^{}`$ $`05^{}`$ $`18.8^{\prime \prime }`$, equinox 2000. Therefore our derived position differs by $`1.40^{}`$ from the Chandra position, which is still within the $`1.8^{}`$ radius circle. Because of the position of the source close to the Galactic center, the possibility of contamination from diffuse emission and other sources, which might fall within the $`1.8^{}`$ radius circle, must be considered. Indeed in the Chandra observation of KS 1731–260, another X-ray source (probably an X-ray active star) was detected at $`0.5^{}`$ from KS 1731–260 (Wijnands et al. 2001a, 2001b). The MECS FOV, with the source we detected, is shown in Figure 1. The position of the Chandra sources are also indicated in the figure. Therefore the luminosity we will derive for KS 1731–260 should be considered as an upper limit.

With this caveat in mind, we derived the luminosity of the source in the following way. Since standard background subtraction methods based on blank sky observations are not applicable, we extracted background counts from nearby source-free regions. The background was taken from a region within the 8 central arcmins of the MECS FOV. We checked that the number of source photons estimated does not vary, within 20%, using background boxes of different sizes and locations. Of the 280 events detected within the $`1.8^{}`$ radius circle centered on the source position, 210 events ($`75\%`$ of the total) are expected to be background events. This corresponds to a $`3.7\sigma `$ detection. From these data we extracted a light curve; no short term variability is significantly detected during the BeppoSAX observation. The count rate corrected for the PSF is $`(1.20\pm 0.33)\times 10^3`$ counts s-1. Because the statistics in the MECS spectrum is not enough to constrain the spectral model, we considered different spectral models, namely a blackbody (with temperatures ranging between 0.2 and 1 keV) and a power law (with photon indices ranging between 1.5 and 2.9). In all cases the column density $`N_H`$ was fixed at $`1.0\times 10^{22}`$ cm-2, consistent with the value obtained from the ROSAT and Chandra observations (Barret et al. 1998; Wijnands et al. 2001a). Within 10% all the models yield a flux of $`1\times 10^{13}`$ erg cm-2 s-1 in the $`1.310`$ keV band. For a blackbody of temperature 0.3 keV (in agreement with the blackbody fitted to the Chandra data of KS 1731–260, Wijnands et al. 2001a), adopting a distance of 7 kpc (Muno et al. 2000), the corresponding bolometric luminosity is $`(7\pm 2)\times 10^{32}`$ erg s-1. This value is compatible with the $`0.510`$ keV bolometric luminosity of $`10^{33}`$ ergs/s deduced from the Chandra observation (Wijnands et al. 2001a). As a cross check we also extracted the data in the MECS FOV using a circular region of $`4^{}`$ radius (within which $`90\%`$ of the source counts are expected) centered on the nominal position of KS 1731–260. In this case the background was measured from a region of the MECS image $`20^{}`$ away from the FOV center and divided by a factor of 0.85 to take into account the reduction of the effective area with the increase of the distance from the FOV center. In this case the background contributes $`80\%`$ of the total source plus background spectrum. The estimated bolometric (unabsorbed) source luminosity was $`1.4\times 10^{33}`$ ergs/s using a blackbody model. In the following we will adopt as a reasonable estimate for the source luminosity in quiescence $`L_q1\times 10^{33}`$ erg s-1.

At the beginning of August 2001, KS 1731–260 was observed, still in its quiescence state (according to the ASM light curve), at 1.4 GHz with the Parkes radiotelescope to search for radio pulses. Taking into account the possibility, suggested by Muno et al. (2000), that the NS belongs to an extremely narrow binary system, with an orbital period of $`1`$ hr, we subdivided the data into segments of the duration of 35 and 4 min, respectively. The analysis of longer segments allows to reach a better sensitivity if the NS is in a wide orbit. On the other hand, using shorter segments, the time of arrival of the radio pulses is not heavily affected by the strong Doppler effect caused by the orbital motion in a narrow binary system, which could reduce the coherence of the signal. We found that no radio signals with millisecond periods are present in the data above a flux limit of 0.21 mJy for the 35 min integration time, and above 0.60 mJy for the 4 min integration time.

## 3 Discussion

The 90 ks BeppoSAX observation of KS 1731–260 carried out on 2001 March 2–5 led to a detection of the source at a luminosity level of $`10^{33}`$ ergs/s. This is similar to the flux determined in a subsequent observation of the source carried out with Chandra on 2001 March 27 (Wijnands et al. 2001a). These results testify that the quiescent X-ray luminosity and spectrum of KS 1731–260 are close to those determined for other NS SXTs. We therefore assume that the X-ray spectrum of KS 1731–260 in quiescence is similar to that of Aql X-1 and Cen X-4 (see e.g. Campana et al. 1998a), the best studied cases. The spectrum of these two sources could be well fit by a soft thermal component (blackbody temperature of $`0.2`$ keV) plus a power-law component with a photon index $`\mathrm{\Gamma }1.5`$. The blackbody component is usually interpreted as thermal emission from a pure hydrogen NS atmosphere (e.g. Rutledge et al. 1999, 2000), while the power-law component is thought to be due to residual accretion or the interaction of a pulsar wind with matter released by the companion star (see e.g. Campana & Stella 2000 and references therein).

In the following we discuss the properties of KS 1731–260 in relation to the different mechanisms that have been proposed to explain the quiescent X-ray emission of NS SXTs. In summary, there exist three sources of energy which are expected to produce X-ray luminosity in quiescence:

* residual accretion onto the NS surface at very low rates (e.g. Stella et al. 1994);
* rotational energy of the NS converted into radiation through the emission from a rotating magnetic dipole, a fraction of which can be emitted in X-rays (e.g. Possenti et al. 2001; Campana et al. 1998b, and references therein);
* thermal energy, stored into the NS during previous phases of accretion, released during quiescence (e.g. Brown, Bildsten, & Rutledge 1998; Colpi et al. 2001; Rutledge et al. 2001).

We will compare the corresponding expectations with the measured upper limit on the X-ray luminosity. Note that the last process has to be always taken into account, as the source was recently active for more than one year at a luminosity of $`10^{37}`$ erg/s and frequently active, at comparable luminosities, since its discovery in 1989. Therefore we firstly derive some constraints resulting from the processes a) and b), and then we combine these with the constraint derived from process c).

Some constraints for the processes a) and b) can be derived if we assume that the coherent pulsations detected during the type I X-ray burst on 1996 July 14 (Smith et al. 1997) correspond to the NS spin. Indeed the high degree of coherence ($`\nu /\mathrm{\Delta }\nu 900`$) is compatible with broadening induced by the finite duration of the oscillations ($`2`$ s). This strongly supports the interpretation of the frequency of these oscillations ($`\nu =523.92`$ Hz) as the rotational frequency of the NS. Note that the frequency separation between the kHz QPOs in KS 1731–260 is $`260`$ Hz (Wijnands & van der Klis 1997), which might be interpreted as the spin frequency of the NS (Strohmayer et al. 1996; Miller, Lamb, & Psaltis 1998), while the burst frequency would be its first overtone in this case (see, however, Muno et al. 2000). For concreteness in the following we will assume the burst frequency as the NS spin frequency and will discuss differences when needed.

It has been demonstrated that a rotating magnetic dipole in vacuum emits electromagnetic dipole radiation. Moreover a wind of relativistic particles associated with magnetospheric currents along the field lines is expected to arise in a rotating NS (e.g. Goldreich & Julian 1969). Both these processes, powered by the rotational energy of the NS, depend on the angle between the NS magnetic moment and the spin axis and compensate in such a way that the total energy emitted is nearly independent of this angle (Bhattacharya & van den Heuvel 1991). Thus the bolometric luminosity of a rotating NS in vacuum can be calculated according to the Larmor’s formula $`L_{\mathrm{bol}}=(2/3c^3)\mu ^2\omega ^4`$, where $`c`$ is the speed of light, $`\mu =B_\mathrm{s}R_{\mathrm{NS}}^3`$ (where $`B_\mathrm{s}`$ is the surface magnetic field at the magnetic equator and $`R_{\mathrm{NS}}`$ is the NS radius), and $`\omega `$ is the angular frequency of the NS. A small fraction of this luminosity is emitted in the radio band, probably as the results of electron accelerations in ultra-strong electric potentials (gaps). A still open question is the location of the accelerator (gap), in the outer magnetosphere, close to the light cylinder radius ($`r_{\mathrm{LC}}=cP/2\pi `$, i.e. the radius at which an object corotating with the NS attains the speed of light), as proposed for the pulsar emission mechanism by Halpern & Ruderman (1993), or close to the magnetic cap, as originally proposed by Arons (1981) and Daugherty & Harding (1982). The observability radio emission can be strongly affected by the matter surrounding the NS. In particular free-free absorption can strongly reduce the radio flux especially at low radio frequencies (see e.g. Burderi & King 1994, and, more recently, D’Amico et al. 2001). However, independent of the observability of the radio pulsar, this radiative regime certainly occurs when the space surrounding the NS is free of matter up to $`r_{\mathrm{LC}}`$ and the pressure of this radiation can overcome the pressure of the accretion flow, thus preventing further accretion (see Ruderman, Shaham, & Tavani 1989; Illarionov & Sunyaev 1975). Therefore, in the following, we will adopt the hypothesis that once the magnetospheric radius is outside the light cylinder radius (and thus the space surrounding the NS up to the light cylinder is free of matter) the NS emits radiation according to the Larmor’s formula. In this case, irrespective of the amount of radiation emitted in the radio band and of its observability, we will say, for shortness, that the radio pulsar is active.

As mentioned above, there are two possibilities: the first possibility, scenario a), is that the magnetospheric radius is inside the light cylinder radius and therefore the radio pulsar is off. In this case, for a non-zero NS magnetic field, we should have some matter flow toward the NS in order to keep the magnetospheric radius small enough. In this scenario we again have two possibilities: a1) the magnetospheric radius is inside the co-rotation radius (see below for a definition), so accretion onto the NS surface is possible; a2) the magnetospheric radius is outside the co-rotation radius (but still inside the light cylinder radius), so the accretion onto the NS is centrifugally inhibited, but an accretion disk can still be present and emit X-rays. The other possibility, scenario b), is that the magnetospheric radius is outside the light cylinder radius and therefore the radio pulsar is on. In this case X-ray emission can be produced by: b1) reprocessing of part of the bolometric luminosity of a rotating NS (Larmor’s radiation formula) into X-rays in a shock front; b2) the intrinsic X-ray emission of the radio pulsar.

Let us consider scenario a), i.e. a non-zero accretion rate during the quiescent phase, and in particular scenario a1). If the NS has an intrinsic magnetic field, we can estimate an upper limit on the NS magnetic moment. Note that a similar argument was applied to an ASCA observation of the Rapid Burster in quiescence by Asai et al. (1996). Also in that case it was concluded that a highly magnetized and rapidly rotating NS can be excluded. The matter accreting onto the NS forms an accretion disk whose inner radius is truncated at the magnetosphere by the interaction of the accretion flow with the magnetic field of the NS. In this case the magnetospheric radius $`r_\mathrm{m}`$ is a fraction $`\varphi 1`$ (an expression for $`\varphi `$ can be found in Burderi et al. 1998111$`\varphi =0.21\alpha ^{4/15}n_{0.615}^{8/27}m^{142/945}[(L_{37}/ϵ)^{8/7}R_6^{8/7}\mu _{26}^{5/7}]^{4/135}`$, see the text for the definition of the symbols.; for $`L10^{33}`$ ergs/s we get $`\varphi 0.2`$) of the Alfvén radius $`R_\mathrm{A}`$ defined as the radius at which the energy density of the (assumed dipolar) NS magnetic field equals the kinetic energy density of the spherically accreting (free falling) matter:

$$R_\mathrm{A}=2.23\times 10^6R_6^{2/7}m^{1/7}\mu _{26}^{4/7}ϵ^{2/7}L_{37}^{2/7}\mathrm{cm}$$
(1)
(see, e.g., Hayakawa 1985), where $`R_6`$ is the NS radius, $`R_{\mathrm{NS}}`$, in units of $`10^6`$ cm, $`m`$ is the NS mass in solar masses, $`\mu _{26}`$ is the NS magnetic moment in units of $`10^{26}`$ G cm3, $`ϵ`$ is the ratio between the specific luminosity and the specific binding energy ($`L=ϵ\times GM\dot{M}/R_{\mathrm{NS}}`$, $`G`$ is the gravitational constant, $`M`$ is the NS mass, $`\dot{M}`$ is the accretion rate), and $`L_{37}`$ is the accretion luminosity in units of $`10^{37}`$ erg/s, respectively. Actually, accretion onto a spinning, magnetized NS is centrifugally inhibited once the magnetospheric radius is outside the corotation radius, the radius at which the Keplerian frequency of the orbiting matter is equal to the NS spin frequency:

$$r_{\mathrm{co}}=1.5\times 10^6P_3^{2/3}m^{1/3}\mathrm{cm}$$
(2)
where $`P_3`$ is the spin period in ms. The condition to allow accretion then reads $`r_\mathrm{m}/r_{\mathrm{co}}1`$. This requires:

$$\mu _{26}0.5\varphi ^{7/4}R_6^{1/2}ϵ^{1/2}L_{37}^{1/2}m^{1/3}P_3^{7/6}.$$
(3)
Adopting our upper limit, $`L_{37}=10^4`$, for the quiescent luminosity we obtain $`\mu _{26}0.28m^{1/3}`$ for $`P_3=1.91`$ ($`\mu _{26}0.63m^{1/3}`$ for $`P_3=3.82`$), where we have assumed $`ϵ=R_6=1`$ and $`\varphi 0.2`$.

If the magnetospheric radius is outside the corotation radius and still inside the light cylinder radius (scenario a2) the accretion is centrifugally inhibited, but the disk can still emit X-rays. This will increase the upper limit on the magnetic field derived above because in the propeller regime the system will be under-luminous, by a factor $`ϵ_{\mathrm{prop}}=R_{\mathrm{NS}}/(2r_{\mathrm{disk}})`$, for a given accretion rate. In this case, the maximum accretion rate for a given luminosity occurs when the efficiency factor $`ϵ`$ is minimum. This occurs when the inner disk radius is the furthest from the NS surface without allowing the switch-on of the radio pulsar, i.e. for $`r_{\mathrm{disk}}r_{\mathrm{LC}}`$. Assuming $`r_\mathrm{m}=r_{\mathrm{LC}}`$ and taking as upper limit on $`\dot{M}`$ the condition $`L_q=L_{\mathrm{disk}}(r_{\mathrm{LC}})=ϵ_{\mathrm{prop}}GM\dot{M}/R_{\mathrm{NS}}`$, with $`ϵ_{\mathrm{prop}}=R_{\mathrm{NS}}/(2r_{\mathrm{LC}})`$, we can calculate the upper limit on the magnetic moment in this case:

$$\mu _{26}11.7\varphi ^{7/4}L_{37}^{1/2}P_3^{9/4}m^{1/4},$$
(4)
where we have assumed $`R_6=1`$. This gives $`\mu _{26}8.4m^{1/4}`$ for $`P_3=1.91`$ (and $`\mu _{26}39.8m^{1/4}`$ for $`P_3=3.82`$), where we have assumed $`L_{37}=10^4`$ for the quiescent luminosity, and $`\varphi 0.2`$ (in fact, as $`\varphi ϵ_{\mathrm{prop}}^{32/945}`$, $`\varphi =0.180.2`$ is also appropriate in the propeller case).

Let us consider now the scenario b). If there is no accretion in quiescence or the magnetospheric radius falls outside the light cylinder radius, we expect that the radio pulsar switches on. The pulsar radiation pressure may overcome the pressure of the accretion disk, thus determining the destruction of the disk and the ejection of matter from the system (see Burderi et al. 2001). However there is the possibility that a fraction $`\eta `$ of this power can be converted into X-rays in a shock front between the emerging radiation and the circumstellar matter (scenario b1). Typical values for $`\eta `$ are in the range $`0.010.1`$ (e.g. Campana et al. 1998b). In this case $`\eta L_{\mathrm{rad}}L_q`$ or

$$\mu _{26}5.1L_{37}^{1/2}P_3^2\eta ^{1/2}.$$
(5)
In our case, assuming the minimum value for the efficiency in the conversion of the rotational energy into X-rays, $`\eta =0.01`$, which implies the highest value for the magnetic moment, we get $`\mu _{26}1.9`$ for $`P_3=1.91`$ (and $`\mu _{26}7.5`$ for $`P_3=3.82`$). Let us consider the possibility that part of the spin-down energy loss is directly emitted in X-rays (scenario b2). Indeed about 40 out of $`1400`$ known radio pulsars have been detected in the X-ray range so far. In this case a correlation has been observed between the X-ray and spin-down luminosities (e.g. Becker & Trumper 1997; Possenti et al. 2001, and references therein). Considering the empirical relation, derived by Possenti et al. (2001) analyzing a sample of 37 pulsars, between the 2–10 keV luminosity, $`L_{37}`$, and the rate of spin-down energy loss, $`L_{\mathrm{rad}}`$,

$$L_{37}=2.51\times 10^{52}L_{\mathrm{rad}}^{1.31},$$
(6)
we can again calculate an upper limit on the magnetic moment in the case of KS 1731–260, which gives $`\mu _{26}8.7`$ for $`P_3=1.91`$ (and $`\mu _{26}34.6`$ for $`P_3=3.82`$).

In the last two cases, the X-ray emission is expected to be non-thermal (a power-law spectrum could be appropriate). Interestingly Wijnands et al. (2001) found that the Chandra spectra of KS 1731–260 in quiescence could be equally well fitted by a blackbody or a power law plus a blackbody (although the power law is not statistically required). In this second case one might think that the thermal component arises from the NS surface and the non-thermal power law is associated with the shock emission discussed above. Indeed, according to the results of the spectral deconvolution of Wijnands et al. (2001), the fraction of the total luminosity emitted in the power law is $`15\%`$, small as compared to the total luminosity. However, because of the low statistics, the spectral deconvolution is not secure and, in the case of our BeppoSAX observation, no statistically significative spectral analysis is possible. For this reason, in deriving our upper limits on the magnetic moment strength in these scenarios, we have adopted the very conservative assumption that the whole flux detected by BeppoSAX can be ascribed to the non-thermal emission.

An upper limit on the NS magnetic field can also be derived from the upper limit on the radio pulsed emission from KS 1731–260. Radio observations of KS 1731–260, taken at the Parkes radiotelescope while the system was in quiescence, showed that the upper limit on the flux of pulsed radio emission is 0.21 mJy using 35 min integration time and 0.60 mJy using 4 min integration time. Equating the rotational energy lost by a magnetized NS with its magnetic dipole radiation, we obtain:

$$P\dot{P}=9.75\times 10^{24}\mu _{26}^2I_{45}^1$$
(7)
where $`P`$ and $`\dot{P}`$ are respectively the period of the radio pulsar and its derivate and $`I_{45}`$ is the NS momentum of inertia in units of $`10^{45}`$ g cm2. Proszynski & Przyibician (1985), using a sample of 275 pulsars, found an empirical relation between the observed 400 MHz luminosity, $`L_{400}`$, and the pulsar parameters:

$$\mathrm{log}(L_{400})=\frac{1}{3}\mathrm{log}\left(\frac{\dot{P}_{15}}{P^3}\right)+A$$
(8)
where $`L_{400}`$ is expressed in mJy kpc2, $`\dot{P}_{15}`$ is in units of $`10^{15}`$ s s-1 and $`A=1.1`$. Kulkarni, Narayan & Romani (1990) found that the same relation could easily fit a sample of 11 recycled binary pulsars. For KS 1731–260, assuming a distance of 7.0 kpc and a spectral index $`\alpha =1.7`$, typical of millisecond pulsars, we obtain for the 400 MHz luminosity an upper limit of 88.47 mJy kpc2 in the case of 35 min integration, and of 250.25 mJy kpc2 in the case of 4 min integration. Using these values and combining equations 7 and 8, we can again calculate an upper limit to the magnetic field of KS 1731–260. If the NS is in a wide orbit, i.e. if the phase shifts introduced by the doppler effect are negligible over an integration time of 35 min, we obtain an upper limit of $`\mu _{26}0.68`$ for $`P_3=1.91`$ (and $`\mu _{26}2.7`$ for $`P_3=3.82`$), assuming $`I_{45}=R_6=1`$. If, instead, the NS in KS 1731–260 belongs to a narrow binary system and then the 4 min observations are more appropriate, we obtain $`\mu _{26}3.1`$ for $`P_3=1.91`$ (and $`\mu _{26}12.6`$ for $`P_3=3.82`$). Note however that free-free absorption from matter surrounding the system (expected if the mass loss from the companion proceeds at some level during the quiescent phase, see Burderi et al. 2001 for a detailed discussion of this possibility in transient systems) can efficiently hamper the detection of the pulsed radio emission from an active radio pulsar.

We will consider now the process c), i.e. the thermal emission from the NS. Deep crustal heating was first developed by Haensel & Zdunick (1990). It was pointed out by Brown et al. (1998) that this scenario was relevant to accreting NSs in transient systems, and may contribute to the emission at low luminosities. This scenario applies if the system has time to reach a steady state, in which the heat deposited during short and frequent outbursts in the NS is equal to the luminosity radiated in quiescence. Colpi et al. (2001) have explored the thermal evolution of a NS undergoing episodes of accretion, lasting for a time $`t_{\mathrm{out}}`$, separated by periods of quiescence, lasting $`t_{\mathrm{rec}}`$ in the hypothesis that the thermal luminosity in quiescence is dominated by emission from the hot NS core, with which the crust and atmosphere of the NS are in thermal equilibrium. Adopting an outburst luminosity of $`10^{37}`$ ergs/s and a quiescent luminosity of $`10^{33}`$ ergs/s, this model predicts a ratio of the recurrence time $`t_{\mathrm{rec}}`$ to the outburst time $`t_{\mathrm{out}}`$ of $`70`$ for a 1.4 $`M_{}`$ NS and less than 4 for a NS more massive than 1.7 $`M_{}`$. A $`t_{\mathrm{rec}}/t_{\mathrm{out}}70`$ would imply that, on average, the total time the source should spend in outburst in 13 yrs (i.e. since its discovery) should be about 2 months. This is in conflict with time history of this source and the fact that the last outburst episode observed by the ASM lasted much more than one year. As already pointed out by Wijnands et al. (2001a), we are therefore left with two possibilities: i) we are witnessing a peculiar period of hyperactivity of the source, which will be followed by a much longer (up to a thousand years); ii) the NS in this system is quite massive, in agreement with the expectation of the standard conservative recycling scenario (e.g. Bhattacharya & van den Heuvel 1991) that the NS has accreted a significant amount of mass and has been spun-up to a millisecond spin period.

However, Rutledge et al. (2001) found that, since KS 1731-260 has recently experienced an extremely long outburst, the thermal luminosity in quiescence could be dominated by emission from the heated NS crust out of thermal equilibrium with the core. Indeed the prediction of the quiescent luminosity from the cooling crust model agrees (within a factor of a few) with the observed bolometric luminosity in quiescence. In this case the model above cannot be applied to the system, given that the core might be much cooler than indicated by the thermal luminosity in quiescence. This implies that the presence of a massive NS is not required to explain the behavior of this system, and that the estimated long outburst recurrence timescale mentioned above should be considered as a lower limit. Rutledge et al. (2001) also calculated the expected thermal evolution of the crust, which should cool in 1–30 years, and therefore the time evolution of the quiescent luminosity for different crustal conductivities and cooling processes, showing that a monitoring of this source can provide valuable information on the crust and core microphysics.

A further constraint on the magnetic moment can be derived considering that the NS spin period must be longer than the equilibrium spin period, namely the Keplerian period at the magnetospheric radius:

$$P_{\mathrm{eq}}=1.81\times 10^3\varphi ^{3/2}R_6^{3/7}\mu _{26}^{6/7}ϵ^{3/7}L_{37}^{3/7}m^{2/7}\mathrm{s}.$$
(9)
Adopting $`\varphi =0.2`$, and an outburst luminosity of $`L_{37}=1`$ (i.e. the averaged luminosity observed with ASM and previous missions when the source was detected) we obtain $`\mu _{26}16.4m^{1/3}`$ for $`P_3=1.91`$ (and $`\mu _{26}36.8m^{1/3}`$ for $`P_3=3.82`$). As KS 1731–260 is a transient source it is possible that, during quiescence, the NS enters a propeller phase during which an efficient torque could spin down the NS far from the equilibrium period defined above. In this case the magnetic moment is significantly less than the value derived in this paragraph. On the other hand, if the propeller phases are infrequent and/or the spin down torque is unefficient, the NS is expected to be in spin equilibrium and $`\mu _{26}=16.4m^{1/3}`$ ($`\mu _{26}=36.8m^{1/3}`$).

From the discussion above we conclude that, in any case, the NS magnetic field is most probably less than $`10^9`$ Gauss, and less than $`4\times 10^9`$ Gauss in the worst case.

Let us finally discuss the consequences of the spin equilibrium constraint when applied to the proposed scenarios in quiescence. In fact, if we consider the case that the NS in KS 1731–260 is quite massive, the amount of mass accreted, $`0.3`$ $`M_{}`$, is sufficient, in principle, to spin up the NS below 1 ms for soft or moderately stiff equations of state for the NS matter and below 1.5 ms even for the stiffest equations of state (see Burderi et al. 1999). Therefore:

* The NS in quiescence accretes matter onto its surface at very low rates. In this case $`\mu _{26}0.28m^{1/3}`$ ($`\mu _{26}0.63m^{1/3}`$). This, compared with the spin equilibrium constraint, i.e. $`\mu _{26}16.4m^{1/3}`$ ($`\mu _{26}36.8m^{1/3}`$), implies that the NS is far from spinning at equilibrium, i.e. propeller phases must be frequent and the spin down torque very effective, although this quiescent phase does not correspond to a propeller.
* The NS in quiescence is in a propeller phase. In this case $`\mu _{26}8.4m^{1/4}`$ ($`\mu _{26}39.8m^{1/4}`$). This, compared with the spin equilibrium constraint, $`\mu _{26}16.4m^{1/3}`$ ($`\mu _{26}36.8m^{1/3}`$), indicates that the NS is compatible with spinning close to the equilibrium. In this case propeller phases must be rare (although the present quiescent phase does correspond to a propeller) and/or the spin down torque not very effective.
* The NS in quiescence is not accreting matter and the radio pulsar is active; a fraction of the power emitted by the rotating NS is converted into X-rays in a shock front between the emerging radiation and the circumstellar matter. In this case $`\mu _{26}1.9`$ ($`\mu _{26}7.5`$). This, compared with the spin equilibrium constraint, implies that the NS is far from spinning at equilibrium, i.e. again propeller phases must be frequent and the spin down torque very effective, although this quiescent phase does not correspond to a propeller.
* The NS in quiescence is not accreting matter and the radio pulsar is active; in this case the fraction of power emitted by the rotating NS that is converted into X-rays in a shock front between the emerging radiation and the circumstellar matter is negligible, and the X-ray emission is the intrinsic emission from the rotating NS. In this case $`\mu _{26}8.7`$ ($`\mu _{26}34.6`$). This implies that the NS is compatible with spinning close to the equilibrium, i.e. propeller phases must be rare (indeed the present quiescent phase does not correspond to a propeller) and/or the spin down torque not very effective.

In conclusion, scenarios a2) and b2) seems to be the most reasonable, and both indicate that the NS is compatible with spinning close to the equilibrium period.

The authors would like to thank R.E. Rutledge and the anonymous referee for enlightening discussions during the preparation of this work. This work was partially supported by a grant from the Italian Ministry of University and Research (Cofin-99-02-02).


format: text

text:
# At the boundary between biological and cultural evolution: The origin of surname distributions

## 1 Introduction

Biological and cultural features of human populations have been traditionally studied by separate disciplines, but the parallelisms between biological and cultural evolution have been put forward by a number of researchers. Already Darwin (1871) pointed out that “the formation of different languages and of distinct species and the proofs that both have been developed through a gradual process are curiously parallel.”

Cultural traits are transmitted from ancestors to their descendence, in a process analogous to inheritance, and are subject to changes, similar to mutations, by interaction between individuals –such as teaching and imitation. Moreover, they usually fulfill a practical purpose, which amounts to being subject to selection. In fact, they enhance the relationships within human groups, defining social entities comparable to certain biological species and populations. The quantification of cultural traits has been attempted only recently. For example, Cavalli-Sforza et al. (1982) applied concepts from the quantitative theory of biological evolution to construct a theory of cultural evolution. They analyzed forty traits, ranging from political preferences to superstitions. Many of these traits are subject to high mutability, since they are influenced by fashion, individual acquaintances, and personal experience, and one does not expect their quantitative properties to be directly comparable to those of any biological feature. Other traits, on the other hand, are better preserved. Among them we find languages and surnames. Language is essential to integrate the individual to society; surnames are historical –though recent– signs of identity in social groups. Quite early, Galton and Watson (1874) dealt with the problem of the extinction of surnames. This problem is equivalent to that of the extinction of a mutant allele in a population, although this relation was only established half a century later (Fisher, 1922; Haldane, 1927), when the first quantitative approaches to biological evolutionary processes took place.

Comparative methods analogous to biological taxonomy are used to determine the degree of similarity between languages. This method returns the genetic classification of linguistic diversity (see for example Greenberg, 1992; Ruhlen, 1992). Recently, a quantitative study of the taxonomy of languages has been carried out (Zanette, 2001), showing that the distribution of the number of subtaxa within a taxa displays scaling properties, quantitatively similar to those disclosed in biological taxonomy by Burlando (1990, 1993). That is, if $`n`$ is the number of subtaxa belonging to a given taxa –say, the number of languages in the Indo-European family, or the number of species in the genus Canis– the fraction $`p(n)`$ of taxa that have precisely $`n`$ subtaxa scales as $`p(n)n^\beta `$. The exponent is found in the interval $`1\beta 2`$ in both cases. This gives quantitative support to Darwin’s observation on the “equivalence” between the mechanisms behind biological and linguistic evolution. A complementary comparison is that of linguistic abundance, measured as the number of individuals speaking a language, with the number of individuals of a biological species. Again, the frequency as a function of the number of individuals has scaling properties both for languages (Gomes et al., 1999) and for species (Pielou, 1969), see Fig. 1.

Surnames are cultural traits (Cavalli-Sforza & Feldman, 1981) whose transmission bears strong similarity with that of some biological features. They are inherited from one of the parents, usually the father, much in the same way as the Y-chromosome or the mitochondrial DNA. The extinction of a surname and the persistence of a non-recombining neutral allele are equivalent problems. This is not only a mathematical fact, but has also practical implications. Indeed, to assess the multiple or single origin of a surname one can turn to genetic measures, since males sharing the same surname might also share the same haplotype in the nonrecombining segment of the Y chromosome (Sykes & Irven, 2000). In a very large population, the statistical properties of the surname distribution can be strongly correlated with genetic diversity (Barrai et al., 1996), and may even be used to understand the genetic structure of a population (Yasuda et al., 1974). Recent reports on actual populations (Miyazima et al., 2000; Zanette & Manrubia, 2001) show that the distribution of surnames follows the same statistical law observed for languages and biological species. Namely, if $`n`$ is the number of individuals bearing a given surname, the fraction $`p(n)`$ of surnames decreases with $`n`$ as $`p(n)n^\gamma `$ (see Fig. 1). Here, however, the exponent is always $`\gamma 2`$. This paper focuses on a model aimed at predicting this kind of regularities, observed in disparate human populations.

## 2 Theoretical approaches to surname evolution

At the mathematical level, several models have been proposed and analysed in order to identify quantitative properties of surnames evolution –or, equivalently, those of nonrecombining neutral alleles. The main points addressed by these studies are (i) the probability of fixation of a given surname/allele in a closed population (Galton & Watson, 1874; Fisher, 1922; Haldane, 1927; Moran, 1962; Lange, 1981; Rannala 1997; Hull 1998), and (ii) the distribution of the number of individuals bearing the same surname/allele (Kimura & Crow, 1964; Karlin & McGregor, 1967; Fox & Lasker, 1983; Panaretos, 1989; Gale 1990; Consul, 1991; Islam, 1995; Zanette & Manrubia, 2001). Indeed, these two questions cover complementary aspects of the same problem. In (i), one deals with a closed population (no immigrants enter the system), and implicitly assumes that the mutation rate is low enough, such that fixation can indeed occur before a mutant form appears. Suppose that there are $`N`$ individuals in the population. We know from coalescence theory that the time $`n_g`$ (in units of the number of generations) required for a neutral allele to be fixed is of order $`n_gN`$. Now suppose that the mutation rate per generation and per individual is $`r`$. Then $`Mrn_gN`$ is the average number of mutants after $`n_g`$ generations have elapsed. Only if $`M1`$, that is if $`rN^2`$, will the fixation of the allele be possible. As soon as this inequality is violated, a new situation arises, in which both neutral drift and mutation play relevant roles. In this case, a broad distribution of surnames or of genetic diversity are expected. Actually, the value of $`r`$ will be usually fixed by the nature of the problem, while the size $`N`$ of the considered population can increase enough such that $`N^2>1/r`$. Therefore, the statistical behaviour crosses over to the second regime, where the appearance of mutants cannot be discarded. We are then in the assumptions of the models of class (ii).

All the models in class (ii) represent systems where statistical equilibrium is reached. While the neutral drift drives the less frequent surnames out of the system, mutations generate new surnames. For large times, the number of different surnames in the system is much larger than unity. Yasuda and co-workers (Yasuda et al., 1974) used a stochastic model for a population with a fixed number of individuals. At each evolution step, a new individual with the surname of his father is added, but he acquires a new surname with a given probability. In any case, a randomly chosen individual is eliminated in order to keep the population constant. Their analytical results compared successfully to field data despite the restriction of constant population, which forces the elimination of individuals. In this model, the size of the progeny is Poisson distributed.

Other models start with a single individual in the population and new individuals carrying the same surname or a new one are sequentially added (Panaretos, 1989; Consul, 1991; Islam, 1995). These fall in the category of branching processes (Harris, 1963) with an increasing number of individuals in the population. Panaretos (1989) rephrased a model introduced by Simon (1955) in the context of linguistics, as follows. At each step, an individual is chosen to be the father of a newborn and his surname is transmitted with probability $`1\alpha `$. With the complementary probability $`\alpha `$, the newborn is assigned a surname not present in the population. Recently, we have modified Simon’s model by (i) introducing an additional parameter $`\mu `$ which represents the death rate of the individuals in the population, and (ii) allowing for arbitrary initial conditions. Preliminary results have been successfully compared with large data sets (Zanette & Manrubia, 2001). This is the starting point of the present work, where we present new analytical and numerical results for this birth-death model.

We describe the model in the next section, where our analytical results are derived. Questions like the role played by the composition of the founder population and by death events in an exponentially growing population are addressed. We also analyse a system with constant population, which turns out to quantitatively differ from the previous case. In Section 4 we numerically check our analytical results and run computer simulations of the model in situations not amenable to an analytical description. In Section 5 our results are tested against field data. We finish with an overall discussion in the last section.

## 3 Birth-death model for surname evolution

Our model population evolves in discrete steps, each step corresponding to the birth of a new individual. At each time step, moreover, an individual is chosen at random from the whole population and becomes removed with probability $`\mu `$, representing a death event. The total population at step $`s`$, $`P(s)`$, is therefore a stochastic process governed by the evolution equation

$$P(s+1)=P(s)+1w(s),$$
(1)
where $`w(s)`$ is the dichotomic stochastic process

$$w(s)=\{\begin{array}{cc}1\hfill & \text{with probability }\mu ,\hfill \\ 0\hfill & \text{with probability }1\mu .\hfill \end{array}$$
(2)
We show in Appendix A that, averaging over realizations of the stochastic process $`w(s)`$, the average total population $`\overline{P}`$ grows exponentially in time:

$$\overline{P}(t)=P_0\mathrm{exp}[\nu (1\mu )t].$$
(3)
Here, $`\nu `$ is the birth rate per individual and unit time, and the product $`\nu \mu `$ turns out to be the corresponding death rate; $`P_0`$ is the initial population.

We think of the population as divided into groups –the families– within which all individuals bear the same surname. At each birth event, the newborn is assigned a new surname, not previously present in the population, with probability $`\alpha `$. This probability can be seen as a mutation rate for surnames, but could also be interpreted as a measure of migration towards the population of individuals with new surnames. With the complementary probability, $`1\alpha `$, a preexistent individual is chosen at random from the population and becomes the newborn’s father, i.e. his surname is given to the newborn. Surnames are therefore assigned with probabilities proportional to the size of the corresponding families, allowing however for fluctuations inherent to the random distribution of births among families. Consequently, the distribution of surnames is driven by a stochastic multiplicative process (Van Kampen, 1981), modulated in turn by the total population growth. This process is analogous to the mechanism proposed by Simon (1955) to model the frequency distribution of words and city sizes, described by Zipf’s law (Zipf, 1949), among other instances. Our model, in fact, reduces to Simon’s model if mortality is neglected, i.e. for $`\mu =0`$. Note that allowing for death events adds a relevant process in the case of surnames, namely, the possibility that a surname disappears if it is borne by a single individual and the individual dies.

Since neither the probability of becoming father of a newborn nor the death probability depend on the individual’s age, the present model population can be thought of as ageless. As a consequence, if $`\nu `$ and $`\mu `$ are constant, the probability that an individual dies at an age between $`T`$ and $`T+dT`$ is

$$dp(T)=\nu \mu \mathrm{exp}(\nu \mu T)dT,$$
(4)
which implies a life expectancy $`\overline{T}=1/\nu \mu `$. Moreover, the probability that an individual has exactly $`m`$ children with the same family name during its whole life equals

$$p(m)=(1\alpha )\mu [1+(1\alpha )\mu ]^{m1},$$
(5)
giving a fertility $`\overline{m}=1/(1\alpha )\mu `$. The exponential distribution in (5) is in reasonable agreement with actual data collected over relatively short periods, for instance, for the United States (Hull, 1998), but contrasts with the Poissonian distribution found in data integrated over several centuries, as in the case of England (Dewdney, 1986). This discrepancy can presumably be ascribed to long-range variations of real birth and death rates (see Fig. 2).

We are assuming that an individual’s surname is inherited from the father. Consequently, the model –as presented above– describes the evolution of the male population only. The same evolution rules apply however if it is assumed that surnames are transmitted with the same probability by either parent. In this case, there is no sex distinction and the model encompasses the whole population. The real situation is in fact intermediate between these two limiting cases: whereas some societies strongly favor inheritance of the surname either from the father or from the mother, other groups allow for a choice between the two possibilities. In some developed western countries, where the surname is preferentially transmitted by the father, opting for the mother’s surname has been proposed as a method to avoid the persistent extinction of surnames, due to the extremely slow population growth (Legay & Vernay, 1999).

### 3.1 Distribution of families by size. The case of $`\mu =0`$

As stated above, our model reduces to Simon’s model (Simon, 1955) for $`\mu =0`$. In this case, $`w(s)=0`$ for all $`s`$ and the total population evolves deterministically, $`P(s)=P_0+s`$, since exactly one individual is added to the population at each step. The number of families with exactly $`i`$ individuals at step $`s`$, $`n_i(s)`$, grows to $`n_i(s)+1`$ when an individual is added to a family of size $`i1`$. This occurs with probability $`(1\alpha )(i1)/P(s)`$. On the average, thus, the number of families with$`i`$ individuals varies according to

$$\overline{n}_i(s+1)=\overline{n}_i(s)+\frac{1\alpha }{P(s)}\left[(i1)\overline{n}_{i1}(s)i\overline{n}_i(s)\right],$$
(6)
for $`i>1`$. To the families with only one individual, on the other hand, the positive contribution comes from the creation of new surnames with probability $`\alpha `$. Therefore,

$$\overline{n}_1(s+1)=\overline{n}_1(s)+\alpha \frac{1\alpha }{P(s)}\overline{n}_1(s).$$
(7)
Note that, since for $`\mu =0`$ the stochastic process $`w(s)`$ becomes trivial, overlines indicate here average over realizations of the stochastic process by which each newborn’s surname is chosen.

The system of equations (6) and (7) can be completely solved for arbitrary initial conditions. In fact, (7) is an autonomous equation for $`\overline{n}_1(s)`$, whose solution reads

$$\overline{n}_1(s)=\frac{\alpha }{2\alpha }(P_0+s)+\left(n_1(0)\frac{\alpha }{2\alpha }P_0\right)\frac{\mathrm{\Gamma }(P_0)\mathrm{\Gamma }(P_0+s1+\alpha )}{\mathrm{\Gamma }(P_0+s)\mathrm{\Gamma }(P_01+\alpha )},$$
(8)
where $`\mathrm{\Gamma }(z)`$ is the gamma-function (Abramowitz & Stegun, 1970). Then, eqn (6) can be used to recursively obtain $`\overline{n}_i(s)`$ for $`i>1`$. For long times, $`s\mathrm{}`$, $`\overline{n}_1(s)`$ is essentially given by the first term in the right-hand side of eqn (8), which grows linearly with $`s`$, as the total population. The second term, which contains all the information about the initial condition, can be seen to decrease as $`s^{(1\alpha )}`$. In this limit, the recursion from eqn (6) can be immediately solved. Thus, for $`s\mathrm{}`$, we find

$$\overline{n}_i(s)=\frac{\alpha }{2\alpha }\frac{\mathrm{\Gamma }(i)\mathrm{\Gamma }\left(\frac{1}{1\alpha }+2\right)}{\mathrm{\Gamma }\left(\frac{1}{1\alpha }+1+i\right)}(P_0+s).$$
(9)
The asymptotic number of families of a given size $`i`$ turns out to be proportional to the total population. For fixed $`s`$ and sufficiently large values of $`i`$, $`\overline{n}_i(s)`$ decreases as a power law, $`\overline{n}_ii^z`$, with

$$z=1+\frac{1}{1\alpha }.$$
(10)
In the limit of small mutation rate, $`\alpha 0`$, we have $`z2`$. As advanced in the Introduction, this is the value observed in actual surname distributions. This same result was obtained by Simon (1955) for a special initial condition. Elsewhere (Zanette & Manrubia, 2001) we have already discussed the fact that, though the effects of the initial condition in Simon’s model fade out for long times, transients can strongly depend on the initial distribution $`n_i(0)`$ (see also section 4). In our case, this aspect could be relevant in populations where modern surnames have appeared recently, like the Japanese (Miyazima et al., 2000). We shall return to this specific point later.

Equations (6) and (7) imply that, in average, the total number of surnames in the population grows linearly:

$$\overline{N}(s)=\underset{i=1}{\overset{\mathrm{}}{}}\overline{n}_i(s)=N_0+\alpha s,$$
(11)
with $`N_0`$ the initial number of different surnames. As a function of time, the number of surnames increases exponentially:

$$\overline{N}(t)=N_0+\alpha P_0[\mathrm{exp}(\nu t)1].$$
(12)
We have already pointed out that, in contrast, the number of surnames in some real populations at present times is decreasing (Cavalli-Sforza & Feldman, 1981; Legay & Vernay, 1999). Under suitable conditions, adding mortality allows to reproduce this particular behaviour.

### 3.2 Effect of death events

Under the action of mortality, the growth of the total population $`P(s)`$ fluctuates stochastically, according to eqn (1), depending on the occurrence of death events at each evolution step. The evolution of $`\overline{n}_i(s)`$ can be implemented in two substeps, as follows. First, eqns (6) and (7) are used to calculate the intermediate values

$$\overline{n}_i^{}(s)=\overline{n}_i(s)+\frac{1\alpha }{P(s)}\left[(i1)\overline{n}_{i1}(s)i\overline{n}_i(s)\right],$$
(13)
and

$$\overline{n}_1^{}(s)=\overline{n}_1(s)+\alpha \frac{1\alpha }{P(s)}\overline{n}_1(s),$$
(14)
and the population is updated to $`P^{}(s)=P(s)+1`$. Second, the evolution due to mortality is applied. In terms of the random process $`w(s)`$ of eqn (2), we have

$$\overline{n}_i(s+1)=\overline{n}_i^{}(s)+\overline{\left[\frac{w(s)}{P^{}(s)}\right]}\left[(i+1)\overline{n}_{i+1}^{}(s)i\overline{n}_i^{}(s)\right],$$
(15)
and $`P(s+1)=P^{}(s)w(s)`$ \[cf. eqn (1)\].

Naive replacement of the stochastic process $`w(s)`$ by its average value, $`w(s)\mu `$, in order to obtain a deterministic approximation to the problem, would lead to an equation for which positive solutions cannot be insured. In fact, it is possible to give initial conditions for such deterministic equation which lead to negative values of $`\overline{n}_i(s)`$ at sufficiently large $`s`$ and $`i`$. A satisfactory deterministic approximation can however be proposed by assuming that the solution $`\overline{n}_i(s)`$ varies slowly both in $`s`$ and $`i`$. This assumption is generically verified for large populations with smooth initial conditions. In this situation, eqns (13) and (15) admit a continuous approximation in terms of real variables $`z`$ and $`y`$, which replace the integer variables $`s`$ and $`i`$, respectively. Meanwhile, $`\overline{n}_i(s)`$ is replaced by a continuous function $`n(y,z)`$.

The continuous approximation can be analyzed at different truncation orders, as discussed in Appendix B. At the first order, the approximate solution to eqns (13) and (15) reads

$$n(y,z)=\alpha \frac{P_0+(1\mu )z}{1\alpha \mu }y^{1(1\mu )/(1\alpha \mu )}$$
(16)
for $`y<y_T(z)`$, and

$$n(y,z)=y_T^1n(y/y_T(z),0)$$
(17)
for $`y>y_T(z)`$, with

$$y_T(z)=\left(1+\frac{1\mu }{P_0}z\right)^{(1\alpha \mu )/(1\mu )}.$$
(18)
In this first-order continuous approximation, thus, the average number of families $`n(y,z)`$ exhibits two separated regimes. For large values of the size variable, $`y>y_T(z)`$, the distribution is essentially determined by the initial condition. At the first evolution steps, where $`z0`$ and $`y_T1`$, this regime covers practically all the domain of variable $`y`$. As time elapses and $`y_T`$ grows, however, this regime recedes and is replaced for $`y<y_T`$ by a power-law distribution, $`ny^\zeta `$, with

$$\zeta =1+\frac{1\mu }{1\alpha \mu }.$$
(19)
Note that, as the exponent $`z`$ in eqn (10), this new exponent coincides with the observed value, $`\zeta 2`$, in the relevant limit $`\alpha 0`$. This result is independent of the death probability $`\mu `$. For sufficiently long times, the power-law regime will be observed for all the relevant family sizes. All the contribution of the initial condition will become restricted to the domain of the largest families.

The boundary $`y_T(z)`$ between the two regimes, eqn (18), grows as a power of the ratio between the average current population $`P_0+(1\mu )z`$ and the initial population $`P_0`$. For $`\alpha 0`$, the exponent is practically equal to unity. Consequently, for the boundary to reach a given value $`y_0`$, in such a way that all the families with sizes below $`y_0`$ are in the power-law regime, the total population must grow by a factor practically equal to $`y_0`$. Taking into account eqn (3), the position of the boundary as a function of real time reads

$$y_T(t)=\mathrm{exp}[\nu (1\alpha \mu )t].$$
(20)
The transient $`t_0`$ needed for the power-law regime to develop up to a given size $`y_0`$ is therefore logarithmic, $`t_0\mathrm{ln}y_0`$.

The average total number of surnames in the continuous approximation is calculated by analogy with eqn (11):

$$\overline{N}(z)=_1^{\mathrm{}}n(y,z)𝑑y.$$
(21)
Using the above first-order solution for $`n(y,z)`$, we find $`\overline{N}(z)=N_0+\alpha z`$ or, as a function of time,

$$\overline{N}(t)=N_0+\frac{\alpha P_0}{1\mu }\{\mathrm{exp}[\nu (1\mu )t]1\},$$
(22)
as shown in Appendix B \[cf. eqn (12)\]. In consequence, within this approximation, the surname diversity always grows exponentially. We point out, however, that this conclusion is valid for sufficiently smooth distributions and for $`\mu <1\alpha `$, two necessary conditions for the continuous approximation to apply to our problem. It could therefore be argued that having a decreasing number of surnames, as found in some modern developed societies (Cavalli-Sforza & Feldman, 1981; Legay & Vernay, 1999), requires violation of the above conditions. For instance, an initial condition with sharp variations would violate the smoothness condition during the first evolution stages. A death probability $`\mu 1`$ would also threaten the validity of the continuous approximation. This is precisely the case of modern developed societies, where the population growth rate is practically vanishing and, as a matter of fact, reaches occasionally negative values –a situation not described by the present model. In Section 3.3 we treat the special case $`\mu =1`$ and show that, in this limit, the total number of surnames can in fact decrease.

Though the first-order continuous approximation gives a rather rough description of the solution to our model, as a piecewise function with a moving discontinuity, it provides a quite clear qualitative picture of how the solution behaves. The growth of the power-law regime as the initial condition recedes is in good qualitative agreement with the evolution observed in numerical realizations of the model. A better analytical approximation is obtained from the second-order truncation. Second-order derivatives would in fact introduce diffusive-like effects in the variable $`y`$, with the consequent smoothing of the discontinuities of the first-order approximation. The second-order equation, however, cannot be analytically solved for arbitrary initial conditions. Nevertheless, it is possible to give an asymptotic approximation for long times, as

$$n(y,z)=\frac{\alpha P(z)}{1\alpha \mu }\left(2\frac{1\alpha \mu }{1\alpha +\mu }\right)^{\zeta 1}y^1U(\zeta 1,0,2\frac{1\alpha \mu }{1\alpha +\mu }y),$$
(23)
where $`U(a,b,x)`$ is the Kummer’s function (see Appendix B). With respect to the first-order approximation, eqn (16), this solution predicts a lower value of $`n(y,z)`$ for small $`y`$. For larger family sizes, however, it behaves as a power law, $`n(y,z)y^\zeta `$ with exactly the same exponent as in eqn (19). The asymptotic behaviour for large $`y`$ in the power-law regime is therefore not modified. It must be pointed out that, for small values of $`\alpha `$, the profile of $`n(y,z)`$ given by eqn (23) results to be quite independent of the death probability over a considerable range of values of $`\mu `$. Only for $`\mu 1\alpha 1`$, where the exponent $`\zeta `$ strongly depends on $`\mu `$, does $`n(y,z)`$ change sensibly as $`\mu `$ is varied. We take advantage of this feature in Section 5, where the second-order approximation is compared with actual data of surname distributions.

### 3.3 The case of constant population, $`\mu =1`$

As argued above, the limit $`\mu =1`$ is relevant to the discussion of the evolution of surname distributions in some modern developed societies. In fact, this limit corresponds to the case where the population growth rate vanishes, a situation which is closely approached, for instance, in many European countries. In this case, eqn (15) is again deterministic, with $`w(s)=1`$ and $`P^{}(s)=P_0+1`$. The total population at any time is $`P(s)=P_0`$.

The main difference between this case and that of $`\mu <1`$ is that, now, the distribution $`\overline{n}_i(s)`$ becomes independent of time for asymptotically large times. This feature is in agreement with the fact that the asymptotic distribution for $`\mu <1`$ is proportional to the total population $`P(s)`$. The asymptotic surname distribution $`\overline{n}_i^{\mathrm{}}`$ is given, for $`\mu =0`$, by the recurrence equations

$$\overline{n}_{i+1}^{\mathrm{}}=\frac{[(2\alpha )P_02i(1\alpha )]i\overline{n}_i^{\mathrm{}}(1\alpha )(P_0+1i)(i1)\overline{n}_{i1}^{\mathrm{}}}{(i+1)[P_0(1\alpha )(i+1)]}$$
(24)
for $`i>1`$, and

$$\overline{n}_2^{\mathrm{}}=\frac{[(2\alpha )P_02(1\alpha )]\overline{n}_1^{\mathrm{}}\alpha P_0^2}{2[P_02(1\alpha )(i+1)]}.$$
(25)
Note that, since the total population is always $`P_0`$, we have $`\overline{n}_i^{\mathrm{}}=0`$ for $`i>P_0`$, as no family can be larger than $`P_0`$.

The solution reveals two well-defined regimes, depending on how the product $`\alpha P_0`$ compares with unity. For $`\alpha P_0>1`$, the asymptotic distribution behaves as

$$\overline{n}_i^{\mathrm{}}\frac{\alpha P_0}{i}(1\alpha )^{i1}$$
(26)
for a vast range of family sizes. Departures from this behaviour are found very close to $`i=P_0`$ only. For $`i>P_0`$, in fact, the distribution must vanish. We stress the remarkable difference between the exponential stationary distribution (26) and the long-time power-law solution obtained for $`\mu <1`$. In this regime, the stationary total number of surnames is

$$\overline{N}^{\mathrm{}}\frac{\alpha P_0}{1\alpha }|\mathrm{ln}\alpha |.$$
(27)
For $`\alpha P_0<1`$, on the other hand, the distribution behaves as a power law, $`\overline{n}_i^{\mathrm{}}i^1`$, over practically the whole range of family sizes. Note that the exponent of this power law is in agreement with eqn (16) in the limit $`\mu 1`$. For $`iP_0`$, however, the distribution deviates from the power law and exhibits a sharp peak. In the limit $`\alpha 0`$, the distribution becomes an isolated peak at $`i=P_0`$, namely, $`\overline{n}_i^{\mathrm{}}=0`$ for $`i<P_0`$ and $`\overline{n}_i^{\mathrm{}}=1`$ for $`i=P_0`$. Therefore, the total number of surnames is $`\overline{N}^{\mathrm{}}=1`$. This special solution describes the well-known case of a closed population with no surname mutations where, by random drift, a single surname survives for asymptotically large times (Cavalli-Sforza & Feldman, 1981). Notice also that this limit corresponds to the birth-death model introduced by Moran (1962) to study probabilities of fixation of alleles when generations overlap.

We conclude that, for a given population $`P_0`$, the asymptotic number of surnames can be very small if the mutation probability $`\alpha `$ is, in turn, small enough. Consequently, in a steady population with many surnames at the initial stage and with a sufficiently low mutation probability, the number of surnames will decrease towards the stationary value as time elapses. Numerical simulations, discussed in detail in the next section, show that for $`\mu `$ just below unity and small $`\alpha `$, an initial transient where the number of surnames decreases temporarily can be observed. Since $`\mu <1`$, however, the population grows steadily and, as a consequence of mutations, $`\overline{N}`$ will also increase in the long run. The situation in modern developed societies is that the population growth rate has been constantly decreasing, to reach values around zero at present times. Starting from a state with a wealth of surnames –due to the combined effect, a few centuries ago, of a high population growth and the frequent appearance of new surnames– these societies have now reached a regime of almost stationary population where the number of surnames decreases. This situation would be reverted if the growth rate could be maintained above zero during substantial periods. Presently, the only mechanism acting in this direction seems to be immigration, which is in addition an effective source of new surnames.

## 4 Numerical results

In this section we present results of numerical realizations of our model, in order to compare with the analytical approximations presented in Section 3.2, and to illustrate the behaviour within the regimes where such approximations do not hold. Emphasis is put on the role of the initial conditions, the duration of transients, and the evolution of the number of different surnames. Also, we discuss actual surname distributions of modern populations in the light of our results.

### 4.1 Role of initial conditions

First, it is worthwhile to illustrate how the shape of the distribution $`\overline{n}_i(s)`$ depends on the initial condition $`n_i(0)`$ \[see also Zanette & Manrubia (2001)\]. We focus the attention on initial distributions of the form $`n_i(0)=N_0`$ for $`i=i_0`$, and $`n_i(0)=0`$ for $`ii_0`$, in which there are $`N_0`$ families of equal size formed by $`i_0`$ individuals each. Consequently, the initial population is $`P_0=i_0N_0`$. In the following, we denote such an initial condition by the pair $`(i_0,N_0)`$. Figure 3 shows four normalized distributions of family sizes,

$$p_i(s)=\frac{\overline{n}_i(s)}{_i\overline{n}_i(s)},$$
(28)
obtained both from the numerical realization of the stochastic birth-death model and from the iterative solution of the deterministic eqns (6) and (7), for $`\mu =0`$. We find a very good agreement between both methods and, at the same time, clearly realize the relevant role of the initial condition in determining the profile of the distribution for large family sizes. These same features are found for other values of the death probability.

The solution of the first-order approximation to our model, eqns (16) to (18), predicts that in the zone of small family sizes –i.e. for $`y<y_T`$ in the continuous variables– the only dependence on the initial condition appears through the quantity $`P_0`$. This means that the distribution of family sizes is not sensitive to a variation of $`i_0`$ and $`N_0`$ as far as their product is kept constant. Note that the three cases $`(4,5)`$, $`(1,20)`$, and $`(20,1)`$ of Fig. 3 share the same value of $`P_0=20`$. For the parameters of the figure, the crossover between the regions of small and large family sizes, given for the continuous variables by eqn (18), should occur at $`i_T470`$. Though two of the distributions indeed have the same profile –with the expected power-law decay– up to that value, the distribution corresponding to the initial condition $`(1,20)`$ deviates considerably below $`i_T`$. This deviation can be attributed to the fact that the initial condition $`(1,20)`$ corresponds to a quite singular distribution, with a high peak at $`i=1`$. A continuous approximation for such a distribution is arguably expected to give a poor description of the real situation.

### 4.2 Computational measurement of transients

Equation (18) defines, within the first-order approximation, the boundary that separates the asymptotic regime and the zone dominated by the initial condition. Alternatively, for a fixed family size, it can be used to determine the transient before the asymptotic distribution builds up at that size. In order to test this aspect of the first-order approximation, we have devised an independent computational method to evaluate the point at which the cut-off between the two regimes actually takes place. According to eqns (16) and (28), the asymptotic normalized distribution of family sizes is $`p_i^{\mathrm{}}p_i(s\mathrm{})(\zeta 1)i^\zeta `$. In numerical realizations, the distribution is expected to adopt values very close to $`p_i^{\mathrm{}}`$ for large $`i`$ and large enough $`s`$. In the simulations, we fix a certain family size $`i_T`$ and stop the calculation at the step $`s_T`$ when the measured value of $`p_{i_T}`$ satisfies, for the first time, $`|p_{i_T}p_{i_T}^{\mathrm{}}|<\mathrm{\Delta }`$. The results presented in the following correspond to the choice $`i_T=10^2`$ and $`\mathrm{\Delta }=10^1`$. Since, as shown in Section 4.1, the agreement between the numerical simulations of the stochastic birth-death model and the iterative solution of the corresponding average equations is very good, we use this second description to measure the transient $`s_T`$.

We focus the analysis on the dependence with the initial conditions, and keep the values of $`\alpha =10^2`$ and $`\mu =0`$ fixed. As far as $`\alpha 0`$, the results are qualitatively the same for other values of $`\mu `$. We consider initial conditions of the form $`(i_0,N_0)`$, as defined in Section 4.1. Figure 4 shows the transient $`s_T`$ measured from the numerical solution to eqns (6) and (7) according to the criterion introduced above, as a function of $`i_0`$ and $`N_0`$. On the one hand, for small $`N_0`$ and, especially, for small $`i_0`$, $`s_T`$ varies quite irregularly. On the other hand, we see that for larger values of $`i_0`$ and $`N_0`$, say $`i_0>10`$ and $`N_0>5`$, $`s_T`$ exhibits a well-defined linear dependence with both parameters. In this regime, the transient is well approximated by

$$s_T=ai_0N_0.$$
(29)
Linear fitting of $`s_T`$ as a function of $`i_0`$ and $`N_0`$ for the data shown in Fig. 4 yields $`a=237\pm 4`$. Note that eqn (29) implies that, for large $`i_0`$ and $`N_0`$ the transient length is proportional to the initial population $`P_0=i_0N_0`$. This is to be compared with the analogous result within the first-order approximation, eqn (18), which implies

$$s_T=\frac{P_0}{1\mu }\left[i_T^{(1\mu )/(1\alpha \mu )}1\right].$$
(30)
The first-order approximation predicts then a proportionality between $`s_T`$ and $`P_0`$ –independently of the values of $`\alpha `$ and $`\mu `$– in full agreement with the numerical results. We point out, however, that the corresponding proportionality coefficients cannot be directly compared. In fact, the coefficient $`a`$ in eqn (29) is evaluated following the computational definition of the transient, and thus depends on the parameter $`\mathrm{\Delta }`$, which has no correlate in the analytical approach. Note, nevertheless, that for $`\alpha 0`$ and $`\mu =0`$ eqn (30) predicts a coefficient approximately equal to $`i_T=100`$, which is of the same order as the numerical result.

### 4.3 Evolution of surname diversity

As discussed in Section 3, the evolution of the average number $`\overline{N}(s)`$ of different surnames is driven in our model by two competing mechanisms. The diversity increases due to the appearance of new surnames at rate $`\alpha `$, and surnames borne by single individuals disappear when the individual in question dies. Within the first-order continuous approximation to our model, for $`\mu <1`$, the average number of different surnames increases linearly with $`s`$ as $`\overline{N}(s)=N_0+\alpha s`$ (see Section 3.2). In this approximation, the death probability $`\mu `$ plays a role in the variation of diversity as a function of time, eqn (22). On the one hand, increase of the surname diversity in a steadily growing population is generally expected when new surnames are created at a constant rate. On the other, it is possible to conceive special situations where the number of different surnames should temporarily decrease, violating the first-order approximation. Imagine, for instance, that the initial population consists of an ensemble of families with only one individual each. For moderate values of $`\mu `$ and small $`\alpha `$, the initial stage would be characterized by the death of some individuals, with the consequent disappearance of their surnames, and no significative appearance of new surnames. Since the total population grows, however, $`\overline{N}(s)`$ will eventually attain a minimum and, from then on, will increase.

To illustrate this situation, we numerically solve eqns (13) to (15) for the initial condition $`(1,20)`$, with $`20`$ families of one individual each. Note that this is the initial condition for which we detected the largest deviations from the first-order approximation in Section 4.1. Figure 5 shows the evolution of $`\overline{N}(s)`$ for different values of $`\mu `$. As expected, an initial transient where the surname diversity drops is found for $`\mu >0`$. The transient is longer and the minimum in $`N(s)`$ is deeper as $`\mu `$ grows. In all cases, however, the subsequent growth of $`\overline{N}(s)`$ is clearly seen. Note that the slope of this growth depends slightly on $`\mu `$, a feature not predicted by the first-order continuous approximation.

These long transients where (for large $`\mu `$ and suitable initial conditions) the number of different surnames is expected to decrease could be relevant to the description of modern populations with declining diversity (Cavalli-Sforza & Feldman, 1981; Legay & Vernay, 1999). A realistic description of this situation should however take into account that, in recent times, the relative values of birth and death rates in actual populations have considerably changed. The general trend to population growth observed worldwide in the nineteenth century has by now been reversed in many developed areas, such as in Europe, where the total population is practically stationary. In our model, this corresponds to an increment in the value of $`\mu `$ to values close to unity. In particular, the death probability must be allowed to vary with time. Additionally, if the description is expected to encompass the periods where the appearance of new surnames was frequent –specifically, the Middle Ages in the case of Europe– $`\alpha `$ should also change during the evolution.

As a qualitative demonstration of the effect of varying the model parameters with time, we focus the attention on a change of $`\mu `$ from a low value to a high value. For the initial condition $`(1,20)`$ we fix $`\mu =0`$ during the first $`s_0`$ evolution steps. Then, $`\mu `$ is left to increase linearly with $`s`$, such that it reaches unity after $`s_1`$ additional steps. In the numerical calculations shown in Fig. 6 we fix $`s_1=10^3`$ and consider several values of $`s_0`$. For $`s<s_0`$, the number of different surnames increases (cf. Fig. 5). Then, as the death probability grows, $`\overline{N}(s)`$ attains a maximum and begins to decrease. For asymptotically large times, it is expected to approach a stationary value, as predicted for the case $`\mu =1`$ in Section 3.3. Interestingly, $`\overline{N}(s)`$ responds faster to the effects of growing $`\mu `$ for smaller values of $`s_0`$, where the contribution of the initial condition is still important during the variation of the death probability. The asymptotic surname distribution seems to be, in this sense, quite robust to the action of varying $`\mu `$. This feature should be related to the fact that the asymptotic distribution is rather insensitive to the value of $`\mu `$, as discussed at the end of Section 3.2.

## 5 Comparison of the model with field data

Finally, we compare some of the analytical results for our model, derived under several approximations in Section 3.2, with actual data from three modern populations. Specifically, we focus the attention on the second-order approximation, eqn (23), for the asymptotic distribution of surnames in the range of small family sizes, where the effect of the initial condition is negligible. In fact, it is virtually impossible to extract information on the distribution of surnames in historical times from the data presently available. We recall that eqn (23) correctly describes the asymptotic power law $`n(y,z)y^\zeta `$ with $`\zeta 2`$, as observed in real data. The validity of our second-order approximation is therefore to be especially evaluated in the region of very small family sizes, where the distribution differs from the power law.

Our three data sets were obtained from surname counts in telephone books. They correspond to (i) the almost 350,000 different surnames of the whole 1996 Argentine telephone book, (ii) the 6,400 surnames beginning by A in the 1996 Berlin telephone book, and (iii) the surnames in five Japanese cities, with populations ranging between $`2\times 10^3`$ and $`2\times 10^5`$, reported by Miyazima et al. (2000). Let us point out that the three populations involved here have considerably different demographic history. The modern Argentine population has predominantly European ancestors, who immigrated mainly in the period 1880-1915 and just after the World War II. Their surname distribution has therefore to be considered as a combined sample from the countries that contributed the immigrants. From the times of the largest immigration waves to the present, both the birth and death rates have substantially changed. As for Berlin, this city was practically abandoned in the late stages of World War II and subsequently repopulated with a mixture of the ancient inhabitants and newcomers from other cities of Germany. In this case, the surname distribution has consequently resulted from a combination of several German regions. The modern population in Berlin, in addition, presents the particularity of having been artificially separated into two practically immiscible groups during several decades, ending in 1989. Since European surnames appeared mostly during the Middle Ages (Legay & Vernay, 1999), the populations that contributed surnames to both Argentina and Berlin are expected to have developed the asymptotic surname distribution on a substantial range of family sizes. This situation contrasts with the case of modern Japanese surnames, that were originated some $`120`$ years ago (Miyazima et al., 2000). In Japan, moreover, the contribution of immigration in the relevant period should have been considerably less important than in Berlin and Argentina.

Explicit evaluation of eqn (23) requires fixing the values of the mutation rate $`\alpha `$ and the death probability $`\mu `$. We know that the mutation rate, i.e. the probability that a new individual acquires a surname not previously present in the population, is very low for any modern society. Since, as far as $`\alpha 0`$, the limit of eqn (23) for $`\alpha 0`$ is well defined, fixing any sufficiently small value for $`\alpha `$ gives a correct description of $`n(y,z)`$. In our comparison with real data, we have taken $`\alpha =10^3`$. As for the death probability $`\mu `$, unfortunately, we have found it impossible to fix a reliable value. For all real populations, mortality has considerably changed during the periods relevant to the evolution of surname frequencies. Moreover, since $`\mu `$ measures the relative frequency of death and birth events –the latter also including, in real populations, arrival of new individuals by immigration– an evaluation of the death probability should also involve a detailed description of immigration effects. On the other hand, as discussed at the end of Section 3.2, the asymptotic profile of $`n(y,z)`$, given by eqn (23), is practically independent of the death probability as far as $`\mu `$ is not close to unity. We therefore decided to fit the field data with eqn (23) allowing $`\mu `$ to vary in order to get the best approximation in the relevant zone of small family sizes.

Figure 7 shows the three data sets and the corresponding fittings with eqn (23). Fittings have been optimized in the domain of small family sizes, where our analytical approximation is expected to hold. In the case of Argentina, the agreement with real data is excellent up to family sizes above $`i=100`$. Only for $`i>200`$ a systematic deviation is observed, where the analytical result overestimates the frequency. According to our discussion in Section 3.2, this deviation would be the remnant contribution of initial conditions. For the Berlin data, the statistics are poorer. However, it is clear that the analytical approximation fits the data with good precision in the whole range shown here. In the case of the Japanese data the fitting is very good for relatively small family sizes, $`i<10`$, but, on the other hand, noticeable systematic deviations appear for $`i>20`$. This agrees with our expectation on the effect of initial conditions. In the $`120`$ years elapsed from the appearance of modern Japanese family names, the population in the cities from which the present data were obtained has increased by a factor of, at most, order $`10`$. Therefore, according to eqn (18), initial conditions should contribute to the distribution for relatively small family sizes, $`i>10`$, just as observed.

## 6 Discussion

We have analysed a birth-death model with overlapping generations, in order to study several statistical properties of monoparental inheritance in large populations. We have focused the analysis on a cultural trait, namely, the distribution of surnames, taking advantage of the availability of big corpora of real data. However, the model applies generally to biological traits associated with non recombining alleles. Our model has two parameters, namely, the probability that a new individual carries a surname not present in the population, $`\alpha `$, and the average number of death events per birth $`\mu `$. For any value of $`\alpha >0`$ the system attains a broad, stationary distribution of surname diversity. If $`\mu <1`$ the total population grows exponentially in time and so does the total number of different surnames. The marginal case $`\mu =1`$ corresponds to constant total population and, on average, to constant surname diversity for asymptotically large times.

Our analytical results for the stationary distribution of surnames frequency are in good agreement with field data for modern human populations in different countries. Through an analysis of the transient time required for this distribution to reach its asymptotic shape, we have shown that some deviations observed in real data might actually reflect the composition of the founder population. This result has implications in the study of polyphyletism. Indeed, if the same surname can have multiple origins and thus the individuals carrying it are not always phylogenetically related, this will affect the shape of the surname distribution. In particular, it is not difficult to estimate the time when surnames originated in a population (using historical records) or, in the biological counterpart, when a mutant allele first appeared (through molecular clock analysis). Then, the approximate cut-off until which the stationary distribution follows the asymptotic shape can be calculated. If the distribution underestimates the frequency of values larger than the cut-off, the system is mainly polyphyletic. If it overestimates that frequency, then the simultaneous appearance of many individual carrying different surnames took place in the past.

The strong resemblance between the cultural inheritance of the surname and the biological process in which non recombining neutral alleles are passed to offspring has justified to apply results from field data in the former case to the latter (Barrai et al., 1996). In the few cases where data on genetic diversity was available, it was possible to retrieve information on past populations by comparing both sets of data (Sykes & Irven, 2000). A specific example comes from the small island of Tristan da Cunha, where the fact that its $`300`$ inhabitants represent only seven surnames and five mitochondrial lineages reflects without doubt the small size of the founder population (Soodyall et al., 1997).

Our results could be applied to population genetics under some hypotheses. Indeed, we are assuming that the number of different alleles at a given locus is practically infinite (this is analogous to the assumption made by Kimura & Crow (1964), in their model of infinitely many alleles), since the possibility of backward mutations is discarded. Nonetheless, this factor could be accounted for just by lowering the value of $`\alpha `$, because this process implies that the “new” mutant is in fact identical to one of the forms present in the population. Since, as long as $`\alpha `$ is small, our results are not sensibly modified, we could also work with a finite but large number $`R`$ of different surnames (equivalently, large genetic polymorphism), and use the same model as long as the current diversity is lower than $`R`$. We have also assumed that the death rate is constant during the life time of individuals, while it is known that the life expectancy depends not only on the age of individuals, but also varies as a function of time (Vaupel et al., 1998). A more realistic model of human inheritance could be constructed by taking into account the variation of $`\mu `$ along the lifetime of each individual.

The evolution of surname diversity follows truly neutral evolution: family names do not fulfill any practical purpose but identifying the lineage of each individual. They cannot be selected for through natural selection. Hence, as we have shown, their statistical distribution closely follows the predictions of a neutral model of monoparental inheritance. It would be interesting to test our theoretical results against the genetic diversity of a large population sample. Unfortunately, however, data on frequencies of allelic polymorphisms are still scarce to carry out a massive study like the one presented here for actual surname distributions. Moreover, one should make sure that the different haplotypes in the sample do not confer any advantage to the individuals, in which case positive feedbacks and deviations from neutral statistics would be expected. Different tests have been proposed in the literature to detect deviations from neutrality (Nielsen, 2001; Fu & Li, 1993), and a number of neutral haplotypes have been positively identified (Sánchez Mazas et al., 1994; Stenoien, 1999). Hopefully, enough data will be available in a near future to calculate reliable diversity distributions in human populations.

## Acknowledgements

We thank M. Montemurro for kindly suplying the whole surname count from the 1996 Argentine telephone book.

## References

Abramowitz, M. & Stegun, I. A. (1970). Handbook of Mathematical Functions. New York: Dover.

Barrai, I., Scapoli, C., Beretta, M., Nesti, C., Mamolini, E. & Rodriguez - Larralde, A. (1996). Isonymy and the genetic structure of Switzerland. I. The distribution of surnames. Ann. Hum. Biol. 23, 431-455.

Burlando, B. (1990). The fractal dimension of taxonomic systems. J. theor. Biol. 146, 99; (1993). The fractal geometry of evolution. J. theor. Biol. 163, 161.

Cavalli-Sforza, L. L. & Feldman, M. W. (1981). Cultural Transmission and Evolution: A Quantitative Approach. Princeton: Princeton University Press.

Cavalli-Sforza, L. L., Feldman, M. W., Chen, K. H. & Dornbusch, S.M. (1982). Theory and observation in cultural transmission. Science 218, 19-27.

Consul, P. C. (1991). Evolution of surnames. Int. Stat. Rev. 59, 271-278.

Darwin, C. (1871). The Descent of Man and Selection in Relation to Sex. London: Murray.

Dewdney, A. K. (1986). Computer recreations: Branching phylogenies of the Paleozoic and the fortunes of English family names. Sci. Am. 254, 12-16.

Fisher, R.A. (1922). On the dominance ratio. Proc. Roy. Soc. Edin. 42, 321-341.

Fox, W.R. & Lasker, G.W. (1983). The distribution of surname frequencies. Int. Stat. Rev. 51, 81-87.

Fu, Y.X. & Li, W.H.  (1993). Statistical tests of neutrality of mutations. Genetics 133, 693-709.

Gale, J.S. (1990). Theoretical population genetics. Unwin Hyman, London.

Galton, F. & Watson, H.W. (1874). On the probability of the extinction of families. J. Roy. Anthropol. Inst. 4, 138-144.

Gomes, M. A. F., Vasconcelos, G. L., Tsang, I. J. & Tsang, I. R. (1999). Scaling relations for diversity of languages. Physica A 271, 489-495.

Greenberg, J. H. (1992). Preliminaries to a systematic comparison between biological and linguistic evolution. In: The evolution of human languages (Hawkins J. A. and Gell-Mann M., eds) Reading: Addison-Wesley.

Haldane, J. B. S. (1927). A mathematical theory of natural and artificial selection. Part V: Selection and mutation. Proc. Camb. Phil. Soc. 26, 838-844.

Harris, Th. E. (1963). The theory of branching processes. Berlin: Springer.

Hull, D. M. (1998). A reconsideration of Galton’s problem (using a two-sex population). Theor. Pop. Biol. 54, 105-116.

Islam, M. N. (1995). A stochastic model for surname evolution. Biom. J. 37, 119-126.

Karlin, S. & McGregor, J. (1967). The number of mutant forms maintained in a population. Proc. fifth Berkeley Symp. Math. Stat. Prob. 4, 415-438.

Kimura, M. & Crow, J.F. (1964). The number of alleles that can be maintained in a finite population. Genetics 49, 725-738.

Lange, K. (1981). Minimum extinction probability for surnames and favorable mutations. Math. Biosci. 54, 71-78.

Legay, J. M. & Vernay, M. (1999). L’ instructive histoire des noms de familie. Pour la Science n. 255, 58-65.

Lotka, A. J. (1931). Population analysis–the extinction of families–I. Washington Acad. Sci. 21, 377-380; ibid. 453-459.

Miyazima, S., Lee, Y., Nagamine, T. & Miyajima, H. (2000). Power-law distribution of family names in Japanese societies. Physica A 278, 282-288.

Moran, P.A.P. (1962). The statistical processes of evolutionary theory. Clarendon Press, Oxford.

Nielsen, R. (2001). Statistical tests of selective neutrality in the age of genomics. Heredity 86, 641-647.

Panaretos, J. (1989). On the evolution of surnames. Int. Stat. Rev. 57, 161-167.

Pielou, E. C. (1969). An introduction to mathematical ecology. New York: Wiley.

Poore, M. E. D. (1968). Studies in Malaysian rainforest. Ecology 56, 143-196.

Rannala, B. (1997). Gene genealogy in a population of variable size. Heredity 78, 417-423.

Ruhlen, M. (1992). An overview of genetic classification. In: The evolution of human languages (Hawkins, J.A. and Gell-Mann M., eds) Reading: Addison-Wesley.

Sánchez Mazas, A., Butler Brunner, E., Excoffier, L., Ghanem, N., Salem, M.B., Breguet, G., Dard, P., Pellegrini, B., Tikkanen, M.J., Langaney, A., Lefranc, G., & Butler, R. (1994). New data for AG haplotype frequencies in caucasoid populations and selective neutrality of the AG polymorphism. Human Biology 66, 27-48.

Simon, H. A. (1955). On a class of skew distribution functions. Biometrika 42, 425-440.

Soodyall, H., Jenkins, T., Mukherjee, A. du Toit, E., Roberts, D.F. & Stoneking, M. (1997). The founding mitochondrial DNA lineages of Tristan da Cunha Islanders. Am. J. Phys. Anthropol. 104, 157-166.

Solé, R. V. and Alonso, D. (1998). Random walks, fractals, and the origin of rainforest diversity. Adv. Complex Systems 1, 203-220.

Stenoinen, H.K. (1999). Are enzyme loci selectively neutral in haploid populations of nonvascular plants? Evolution 53, 1050-1059.

Sykes, B. & Irven, C. (2000). Surnames and the Y chromosome. Am. J. Hum. Genet. 66, 1417-1419.

van Kampen, N. G. (1981). Stochastic Processes in Physics and Chemistry. Amsterdam: North-Holland.

Vaupel, J. W., Carey, J. R., Christensen, K., Johnson, T. E., Yashin, A. I., Holm, N. V., Iachine, I. A., Kannisto, V., Khazaeli, A. A., Liedo, P., Longo, V. D., Zeng, Y., Manton, K. G. & Curtsinger, J. W. (1998). Biodemographic trajectories of longevity. Science 280, 855-860.

Yasuda, N., Cavalli-Sforza, L. L., Skolnick, M. & Moroni, A. (1974). The evolution of surnames: An analysis of their distribution and extinction. Theor. Pop. Biol. 5, 123-142.

Zanette, D. H. (2001). Statistical regularities in the taxonomic classification of human languages. Adv. Complex Systems 4, 281-286.

Zanette, D. H. & Manrubia, S.C. (2001). Vertical transmission of culture and the distribution of family names. Physica A 295, 1-8.

Zipf, G. K. (1949) Human Behavior and the Principle of Least Effort. Cambridge: Addison-Wesley.

## Appendix A Average evolution of the total population

As explained in the main text, at each step, the total population $`P(s)`$ in our model either increases by $`1`$ with probability $`1\mu `$, due to the occurrence of a birth event and no death event, or remains constant with probability $`\mu `$, due to the occurrence of both a birth and a death event. Equations (1) and (2) quantify the corresponding stochastic process.

Since the total population changes during the evolution, the time interval $`\delta t`$ to be associated with each step –which corresponds to the interval between consecutive births– must also change. In fact, the birth frequency is proportional to the total population, so that $`\delta t`$ is inversely proportional to $`P(t)`$. We write, at step $`s`$,

$$\delta t(s)=\frac{1}{\nu P(s)},$$
(31)
where the frequency $`\nu `$ fixes time units.

In consequence, the real-time average variation in $`P(s)`$ can be obtained from

$$\frac{d\overline{P}}{dt}\overline{\left[\frac{P(s+1)P(s)}{\delta t(s)}\right]}=\nu \overline{[1w(s)]P(s)}=\nu (1\mu )\overline{P},$$
(32)
where overlines indicate average over realizations of the stochastic process $`w(s)`$. To obtain this average evolution equation we have used eqn (1), and have taken into account that $`w(s)`$ and $`P(s)`$ are independent stochastic processes at each step $`s`$, so that $`\overline{w(s)P(s)}=\overline{w}(s)\overline{P}(s)=\mu \overline{P}(s)`$.

If $`\nu `$ is identified with the birth rate per individual and unit time, the product $`\nu \mu `$ is the mortality rate. For constant $`\nu `$ and $`\mu `$, the solution to eqn (32) is

$$\overline{P}(t)=P_0\mathrm{exp}[\nu (1\mu )t],$$
(33)
with $`P_0`$ the initial population.

## Appendix B Continuous approximation to the distribution of family sizes

Under the action of both birth and death events, the evolution of the average number of families of size $`i`$, $`\overline{n}_i`$, is given by eqns (13) and (15). These equations can be approximately solved assuming that the solution varies slowly on $`s`$ and $`i`$, so that these two discrete variables admit replacement by continuous variables $`z`$ and $`y`$, respectively. Accordingly, $`\overline{n}_i(s)`$ is replaced by a continuous function $`n(y,z)`$. The approximation is based on the expansions $`\overline{n}_{i\pm 1}(s)n(y,z)\pm _yn(y,z)+\mathrm{}`$ and $`\overline{n}_i(s+1)n(y,z)+_zn(y,z)+\mathrm{}`$, to be introduced in eqns (13) and (15) at different truncation orders.

### B.1 First-order approximation

To the first order, we obtain the differential equation

$$\frac{n}{z}+\frac{1\alpha \mu }{P_0+(1\mu )z}\frac{}{y}(yn)=0$$
(34)
for $`y>1`$. The contribution for $`y=1`$, given by eqn (14), can be incorporated to (34) either as a boundary condition or as a singular inhomogeneity in the equation itself. The latter yields

$$\frac{n}{z}+\frac{1\alpha \mu }{P_0+(1\mu )z}\frac{}{y}(yn)=\alpha \delta (y1),$$
(35)
where $`\delta (x)`$ stands for the Dirac delta distribution. In the following, we obtain and analyse the solution to this equation.

Introducing the auxiliary variables

$$\xi =\mathrm{ln}\left(1+\frac{1\mu }{P_0}z\right),\eta =\mathrm{ln}y,$$
(36)
eqn (35) becomes

$$\frac{f}{\xi }+\frac{1\alpha \mu }{1\mu }\frac{f}{\eta }=\frac{\alpha }{1\mu }\delta (\eta ).$$
(37)
with $`f(\eta ,\xi )yn(y,z)`$. This is a one-dimensional wave equation, with “spatial” variable $`\eta `$ and “temporal” variable $`\xi `$. It describes shape-preserving advection of the “density” $`f`$ with constant velocity $`v=(1\alpha \mu )/(1\mu )`$, subject to the action of a point source of intensity $`\alpha /(1\mu )`$ at $`\eta =0`$. In our problem, thus, this equation makes sense for $`v>0`$ ($`\mu <1\alpha `$) since $`\eta `$ must be non-negative ($`yi1`$). The general solution for $`\eta 0`$ is given by an arbitrary combination of functions of the form $`f(\eta ,\xi )=f_0(\eta v\xi )`$, where $`f_0`$ is in principle arbitrary. The combination must be chosen in such a way that the boundary and initial conditions are satisfied. In our case, this is achieved with a combination of two such functions. In the original variables, the solution is given piecewise as

$$n(y,z)=\alpha \frac{P_0+(1\mu )z}{1\alpha \mu }y^{1(1\mu )/(1\alpha \mu )}$$
(38)
for $`y<y_T(z)`$, and

$$n(y,z)=y_T^1n(y/y_T(z),0)$$
(39)
for $`y>y_T(z)`$. Here, $`n(y,0)`$ is the initial distribution. The transition point between the two pieces is located at

$$y_T(z)=\left(1+\frac{1\mu }{P_0}z\right)^{(1\alpha \mu )/(1\mu )}.$$
(40)
These expressions give the first-order continuous approximation $`n(y,t)`$ to the distribution of families by size.

The average total number of surnames in the continuous approximation is defined in (21). For the first-order approximation we find

$$\overline{N}(z)=_1^{\mathrm{}}n(y,z)𝑑y=\alpha \frac{P_0+(1\mu )z}{1\alpha \mu }_1^{y_T}y^\zeta 𝑑y+_{y_T}^{\mathrm{}}n(y/y_T,0)\frac{dy}{y_T}.$$
(41)
The last integral is clearly equal to the initial number of surnames, $`N_0`$, as may be immediately realized by the change of variables $`y/y_Ty`$. Explicit calculation of the first integral shows that $`\overline{N}(z)=N_0+\alpha z`$, i.e. exactly the same result as for the case of $`\mu =0`$, eqn (11). As a function of time, we have

$$\overline{N}(t)=N_0+\frac{\alpha P_0}{1\mu }\{\mathrm{exp}[\nu (1\mu )t]1\}.$$
(42)
### B.2 Second-order approximation

Truncation to the second order in the continuous approximation to eqns (13) and (15) yields

$$\frac{n}{z}+\frac{1\alpha \mu }{P_0+(1\mu )z}\frac{}{y}(yn)+\frac{1\alpha +\mu }{2[P_0+(1\mu )z]}\frac{^2}{y^2}(yn)=\alpha \delta (y1).$$
(43)
Unfortunately, this equation cannot be analytically solved for arbitrary initial conditions. However, a particular solution can be found in the form of a separate function, $`n(y,z)=h(y)P(z)=h(y)[P_0+(1\mu )z]`$. Comparison with eqns (9) and (16) suggests that this particular solution will correspond to the long-time asymptotic evolution. It reads

$$n(y,z)=\frac{\alpha P(z)}{1\alpha \mu }\left(2\frac{1\alpha \mu }{1\alpha +\mu }\right)^{\zeta 1}y^1U(\zeta 1,0,2\frac{1\alpha \mu }{1\alpha +\mu }y),$$
(44)
where $`U(a,b,x)`$ is the logarithmic Kummer’s function (Abramowitz & Stegun, 1970).