text:
# Dispersion in the lifetime and accretion rate of T Tauri discs

## 1 Introduction

Observations of young stars in star forming regions such as Taurus Aurigae show that there is a wide dispersion in the lifetimes of circumstellar discs, in the range 1-10 Myr (Strom et al. 1989; Haisch, Lada & Lada 2001). Although disc lifetimes of this order are consistent with the meteoretic evidence from our own Solar System (Russell et al. 1996), the reason for the large dispersion in the disc lifetime remains unknown. The obvious controlling parameter would be the stellar mass, but in fact observations of low mass stars in young clusters show little correlation between mass and disc lifetime (e.g. Rebull et al. 2002).

The simplest disc models also fail to reproduce other properties of circumstellar discs. One problem is to match the rapid dispersal of the disc. An accretion disc evolves on a single, viscous timescale, which in order to reproduce the aforementioned disc lifetimes needs to be of the order of $`10^6\mathrm{yr}`$. Dispersal of the disc on this timescale would lead to the prediction of a substantial class of objects with properties intermediate between accreting Classical T Tauri stars (CTTS) and discless Weak-lined T Tauri stars (WTTS) (Armitage, Clarke & Tout 1999). In fact, few stars are caught in the act of clearing their discs, which constrains the disc dispersal time to be short – of the order of $`10^5\mathrm{yr}`$ (Simon & Prato 1995; Wolk & Walter 1996; Duvert et al. 2000).

A second problem is that disc evolution is largely indifferent to the presence or absence of close binary companions. The viscous timescale increases with the characteristic radius of the disc, so discs hemmed in by close binary companions ought to evolve more rapidly than discs around relatively isolated stars. Although examples are known of binaries that appear to have dissipated their discs at an early epoch (Meyer et al. 1997), surveys suggest that binary companions as close as 10 au generally have little effect on the lifetime or mass accretion rate of the disc surrounding the primary (Simon & Prato 1995; Bouvier, Rigaut & Nadeau 1997; Ghez, White & Simon 1997; White & Ghez 2001).

Recently, Clarke, Gendrin & Sotomayor (2001) have shown that models which include mass loss from the disc lead to a dispersal timescale that is much shorter than the viscous evolution timescale, reproducing the two-timescale behavior required by the observations. In this paper we demonstrate that the same class of models, when combined with a modest range of initial disc parameters, are quantitatively consistent with the wide dispersion in disc lifetimes. In Section 2 we set out the disc model and show that when mass loss is included, the evolution of the accretion rate depends mainly on the initial disc mass and viscous timescale at large radius, and only weakly on the radial dependence of disc parameters such as the surface density and viscosity. In Section 3, we experiment with the range of initial conditions that are required to reproduce the observed age distributions of CTTS and WTTS in Taurus-Aurigae, and show that this range also reproduces the observed scatter in accretion rates in CTTS inferred by Gullbring et. al. (1998). In Section 4, we apply the model to binary star systems, and show that due to the wind mass loss from large radii in the disc, the disc lifetime may increase with decreasing disc binary separation over the binary separation range $`10100`$ au. Section 5 summarizes our conclusions.

## 2 Disc model

### 2.1 Evolution equation

The evolution of the surface density, $`\mathrm{\Sigma }(r,t)`$, of a geometrically thin accretion disc is described by (Lynden-Bell & Pringle 1974; Pringle 1981),

$$\frac{\mathrm{\Sigma }}{t}=\frac{3}{r}\frac{}{r}\left[r^{1/2}\frac{}{r}\left(\nu \mathrm{\Sigma }r^{1/2}\right)\right]+\dot{\mathrm{\Sigma }}_{\mathrm{wind}}\left(t\right),$$
(1)
where the term $`\dot{\mathrm{\Sigma }}_{\mathrm{wind}}`$ represents the rate of mass loss per unit surface area of the disc in the equatorial plane. We assume that any mass lost leaves the disc with Keplerian specific angular momentum, and thus does not drive radial inflow or outflow in the remaining disc gas. For the kinematic viscosity $`\nu `$ we adopt a simple power-law in radius,

$$\nu =\nu _0r^\beta ,$$
(2)
with $`\nu _0`$ and $`\beta `$ constants. Within the context of $`\alpha `$ disc models (Shakura & Sunyaev 1973), the viscosity in a protoplanetary disc is determined by the intrinsic efficiency of angular momentum transport (i.e. the value of $`\alpha `$); by the relative contribution of viscous heating and stellar irradiation; and by the opacity which controls the vertical disc structure. A passive irradiated disc, in which both the central temperature and the surface temperature decline with radius as $`T_cr^{1/2}`$, has a viscosity (assuming a constant $`\alpha `$) which scales as $`\nu r`$. More elaborate calculations of protoplanetary disc structure – though still within the $`\alpha `$ disc formalism – likewise yield $`\beta 1`$ (e.g. Bell et al. 1997; D’Alessio et al. 1998). Only limited observational constraints are available at the (small) radii of interest, though $`\beta =1`$ is consistent with VLA observations of TW Hydrae by Wilner et al. (2000).

In a steady-state, the surface density profile corresponding to a given mass accretion rate $`\dot{M}`$ is,

$$\nu \mathrm{\Sigma }=\frac{\dot{M}}{3\pi }\left(1\sqrt{\frac{r_{\mathrm{in}}}{r}}\right),$$
(3)
where we have assumed standard zero-torque boundary conditions at $`r_{\mathrm{in}}`$, the inner edge of the disc.

Equation (1) is a diffusion equation for the surface density. Provided that the surface density evolution is controlled by angular momentum transport rather than mass loss, the characteristic timescale for evolution is the viscous time,

$$t_\nu =\frac{r^2}{\nu }.$$
(4)
Since the viscous timescale increases with radius, the evolution of a disc described by equation (1) speeds up at late times if the outer edge of the disc is truncated by mass loss.

Standard numerical methods (e.g. Pringle, Verbunt & Wade 1986) are employed to solve equation (1) on a non-uniform radial mesh. We use zero-torque ($`\mathrm{\Sigma }(r_{\mathrm{in}})=0`$) boundary conditions at the inner edge of the disc. The spatial domain is large enough that the outer boundary condition has no influence upon the results. Formally, we set the radial velocity $`v_r=0`$ at this boundary.

### 2.2 Disc mass loss

Observations provide direct evidence for the existence of two types of disc outflow, jets and photoevaporative flows. Of these, photoevaporative flows, which lead to mass loss from the outer regions of the disc, are the most obviously promising for producing two-timescale disc evolution and rapid disc dispersal (Clarke, Gendrin & Sotomayor 2001; Matsuyama, Johnstone & Hartmann 2002). We consider mass loss via photoevaporation exclusively in this paper, while noting that almost identical behaviour is also possible in models of magnetically driven outflows in which mass loss extends beyond the inner few au of the disc (e.g. Ostriker 1997).

Photoevaporative flows are thermal winds that occur when the disc surface is heated by exposure to ultraviolet radiation (Hollenbach et al. 1994; Johnstone, Hollenbach & Bally 1998; Hollenbach, Yorke & Johnstone 2000; Richling & Yorke 2000). The heated gas can escape at radii $`r>r_g`$, where the critical radius $`r_g`$ is given (approximately) by the condition that the sound speed $`c_s`$ of the hot gas exceeds the local escape velocity,

$$r_g\frac{GM}{c_s^2}.$$
(5)
For parameters appropriate to T Tauri stars, $`r_g510\mathrm{au}`$. The flux of photoevaporating radiation can be external, as in Orion where UV radiation from massive stars is responsible (Churchwell et al. 1997; Johnstone, Hollenbach & Bally 1998; Scally & Clarke 2001; Henney et al. 2002). Alternatively, a lower level of mass loss can result from the UV and X-ray flux of the low mass star itself (Shu, Johnstone & Hollenbach 1993).

In Taurus, there are no massive stars which could generate a strong external flux of photoevaporating radiation. We therefore adopt a simple mass loss prescription appropriate for the case where the mass loss is driven by local irradiation (Hollenbach et al. 1994). The mass loss per unit area of the disc, $`\dot{\mathrm{\Sigma }}_{\mathrm{wind}}`$, is a function of radius but not of time111Implicitly, this assumes that the photoevaporating flux is not tied to the accretion rate. Models by Matsuyama, Johnstone & Hartmann (2002) show that photoevaporation driven solely by the star-disc accretion shock fails to disperse the disc quickly enough., and is given by,

$`\dot{\mathrm{\Sigma }}_{\mathrm{wind}}`$ $`=`$ $`0,r<r_g`$
$`\dot{\mathrm{\Sigma }}_{\mathrm{wind}}`$ $``$ $`r^{5/2},r>r_g.`$ (6)
We express the normalization of the mass loss via the parameter $`\dot{M}_{\mathrm{wind}}`$, which is defined as the total mass loss rate if the disc extends to 25 au. The instantaneous mass loss rate will therefore differ from this depending upon the outer disc radius. We take $`r_g=5\mathrm{au}`$.

### 2.3 Model parameters

The model described above is defined by two parameters describing the disc physics ($`\nu _0`$ and $`\beta `$), two describing the wind ($`\dot{M}_{\mathrm{wind}}`$ and $`r_g`$), and two for the initial conditions ($`\dot{M}_{\mathrm{init}}`$ and the initial disc mass)222Equivalently, we could parameterize the initial conditions in terms of the initial mass and angular momentum content of the disc.. Sensible values for most of these parameters can be estimated from existing observations. We adopt $`\dot{M}_{\mathrm{init}}=5\times 10^8M_{}\mathrm{yr}^1`$, since most T Tauri stars have accretion rates of this order or lower (Gullbring et al. 1998). With this initial accretion rate, a disc lifetime of 1-10 Myr mandates an initial disc mass of $`0.1M_{}`$. After 2 Myr of evolution, this initial disc mass has declined (see Table 1) to a few hundredths of a Solar mass, a value which is broadly consistent with estimates of T Tauri disc masses (Beckwith & Sargent 1991). The assumption that mass loss is driven by photoevaporation fixes $`r_g`$ to be 5-10 au, and the requirement that this mass loss drives rapid dispersal of the disc suggests that $`\dot{M}_{\mathrm{wind}}`$ be no more than a couple of orders of magnitude smaller than $`\dot{M}_{\mathrm{init}}`$ (Clarke, Gendrin & Sotomayor 2001).

Determining the parameters for the disc viscosity is more difficult. We have already noted that theoretical models for protoplanetary discs favour values for $`\beta `$ around unity, but no-one would claim that this parameter is as yet accurately determined. We have therefore chosen for our fiducial models combinations of $`\nu _0`$ and $`\beta `$, shown in Table 1 along with the other parameters of the disc model, that yield a disc lifetime of approximately 5 Myr. This choice of disc lifetime fixes the viscous time at a radius of around 10 au, but does not constrain either $`\nu _0`$ or $`\beta `$ individually. Indeed, as shown in Figure 1, disc models that include a disc wind have only a limited period during which the accretion rate drops as a power-law defined by the value of $`\beta `$ in the viscosity law. As a consequence of this, models with different choices of $`\beta `$ can be constructed that yield almost identical evolution of the disc accretion rate, even though they have entirely different surface density profiles. This degeneracy implies that it is not possible to discriminate observationally – even if we had perfect data – between models with different $`\beta `$, using only measurements of the evolution of the accretion rate with time. We cannot, therefore, use the accretion rate evolution measured by Gullbring et al. (1998) to constrain our models in the same way as was done for disc models without mass loss by Hartmann et al. (1998) and Stepinski (1998). In the following discussion, we will leave $`\beta `$ as an unconstrained parameter, and consider the three models outlined in Table 1 as equally valid possible descriptions of the disc.

To illustrate how the results depend upon some of the other model parameters, we show in Fig. 2 how the accretion rate evolves with time for different choices of the initial disc mass, initial accretion rate, and disc mass loss rate. We start with a disc of mass $`M_{\mathrm{disc}}`$, and accretion rate $`\dot{M}_{\mathrm{init}}`$. The disc has an initial surface density profile defined by equation (3) out to the radius where the enclosed mass is $`M_{\mathrm{disc}}`$, and $`\mathrm{\Sigma }=0`$ at greater radii. As is obvious from Fig. 2, the two-timescale behaviour discussed by Clarke, Gendrin & Sotomayor (2001) occurs in a range of disc models which incorporate mass loss. The decline of the accretion rate steepens dramatically in all models at the epoch when mass loss in the disc wind becomes important.

## 3 The disc fraction in Taurus

In a recent study, Haisch, Lada & Lada (2001) determined (using JHKL band photometric surveys)the fraction of stars possessing discs as a function of the previously published mean cluster age. They found an evolution in the disc fraction shown in the lower panel of Figure 4. However, since the spread in stellar ages can be comparable to the mean age for clusters that are only $`10^6`$ yr old, it is not immediately obvious whether the smooth decline in disc fraction with ‘age’ seen in the Haisch, Lada & Lada (2001) sample reflects primarily,

* A true dispersion in the lifetime of circumstellar discs.
* A convolution of a constant disc lifetime with an extended epoch of star formation in each cluster.

To distinguish between these possibilities, it is necessary to date the stars within young clusters on a star by star basis. We present such an analysis below, and show that in Taurus the evidence suggests that for low mass stars there is a significant intrinsic dispersion in disc lifetime.

### 3.1 Data

To study the evolution of the disc fraction with stellar age, we make use of data compiled for a study of star formation in Taurus by Palla & Stahler (2002). The sample of 151 stars includes approximately 10 stars from Wichmann et al. (2000) whose pre-main-sequence status is confirmed by lithium observations, and which fall within the boundaries of the 12CO map that contains the bulk of the Taurus-Auriga population. ROSAT All Sky Survey sources that fall at larger distances have not been included. Of this sample, 67 stars are classified as CTTS, and 57 are WTTS. For 27 stars discovered in more recent surveys we were unable to find published classifications, and these are excluded from further analysis.

For each of the stars in the CTTS and WTTS samples, we have estimated ages using the evolutionary tracks of Palla & Stahler (1999). In general, age estimates for pre-main-sequence stars are subject to considerable uncertainty, due to different assumptions as to the initial conditions and accretion history (Tout, Livio & Bonnell 1999; Hartmann 2001; Baraffe et al. 2002). In the specific case of Taurus, Hartmann (2003) has argued that the ages of the higher mass stars (with masses $`\stackrel{>}{}1M_{}`$) are especially uncertain, due to different treatments of the stellar birthline (e.g. Hartmann, Cassen & Kenyon 1997; Palla & Stahler 1999) and possible foreground contamination. We therefore further restrict the sample under consideration to those stars with effective temperatures $`T_e4350\mathrm{K}`$, leaving a final sample of 57 CTTS and 41 WTTS. With the effective temperature cut, these are all low mass stars.

Figure 3 shows histograms of the age distribution of the low mass CTTS and WTTS in Taurus. The most immediately obvious conclusion to be drawn from the Figure is that there is substantial overlap in the derived ages of the CTTS and the WTTS. However, there is also evidence that the WTTS are significantly older as a class. A KS test reveals that the age distributions of the WTTS and CTTS samples differ at more than the 2$`\sigma `$ level (the probability that the distributions are the same is approximately 0.03).

Although the empirical conclusion that the WTTS are significantly older than the CTTS falls short of statistical proof, the evidence is consistent with our expectations if CTTS turn into WTTS. Indeed, if no such difference were found, one would either have to posit no evolutionary link between the two (thus begging the question of what CTTS should turn into) or else appeal to contrived models in which the lifetime of discs were shorter in the case of stars born more recently. We shall show below that the Taurus data is instead compatible with the simple hypothesis that CTTS turn into WTTS, but that the clock for this process varies from star to star. In Chamaeleon I, by contrast, it has been claimed that the age distribution of CTTS and WTTS is indistinguishable (Lawson, Feigelson & Huenemorder 1996). It is now known, however, that this dataset was missing a number of young CTTS (Persi et al 2000), and a more complete census of this region may not ultimately yield an answer different from the Taurus result reported here.

Figure 4 shows the CTTS fraction in Taurus as a function of stellar age. We have binned the data into 5 bins, evenly spaced in the $`\mathrm{log}`$ of the estimated stellar age. Most of the transition between CTTS and WTTS occurs for ages between 1 Myr and 10 Myr. At the youngest inferred ages, CTTS are dominant, but there is still a significant admixture of WTTS. We note, however, that for such young ages – less than a Myr – theoretical uncertainties in the age are probably comparable to the age itself (Baraffe et al. 2002). Any conclusions we might draw as to the evolution of the CTTS fraction prior to 1 Myr should therefore be treated with caution.

Figure 4 also shows the evolution of the disc fraction as measured by Haisch, Lada & Lada (2001) in a sample of 6 young clusters with estimated ages ranging between 0.3 Myr and 30 Myr. Within the errors, which are substantial (arising both from the limited number of stars, and from the uncertainty in cluster age which is not shown in the Figure), the same trend with age is seen as for the Taurus data.

### 3.2 A model for the dispersion in disc lifetime

A dispersion in the disc lifetime could be caused by several factors, including different disc initial conditions or varying rates of mass loss. For photoevaporation, the mass loss rate is predicted to vary only weakly (as the square root) with the ionizing flux. We have therefore explored whether the data shown in Fig. 4 can instead be explained as arising from different initial conditions for the discs. We focus on the effect of different initial disc masses, as Fig. 2 indicates that the disc lifetime changes markedly as we vary this parameter.

Fig. 5 shows how the disc lifetime varies with the initial disc mass. We define the disc lifetime as the time at which the accretion rate first drops below $`10^{10}M_{}\mathrm{yr}^1`$ (once this epoch is reached, as shown in Figure 1, mass loss in the wind dominates the evolution, and the accretion rate is dropping precipitously). Apart from changing $`M_{\mathrm{disc}}`$, the other parameters are kept fixed at the values shown in Table 1 (i.e. we change the initial conditions for the disc, but keep the parameters describing the disc viscosity and disc wind constant). For all three disc models, the disc lifetime has an approximately power-law dependence upon the initial disc mass, with an index of $`0.5`$.

Figure 6 shows the predicted disc fraction as a function of age, under the assumption that $`\mathrm{log}(M_{\mathrm{disc}})`$ is distributed as a Gaussian around the central values given in Table 1. We find that the general trend in the plot of disc fraction with age can be reproduced provided that there is a 1 sigma dispersion of approximately 0.5 in the $`\mathrm{log}`$ of the initial disc mass. A dispersion in initial disc mass of this magnitude leads to almost all stars losing their discs between 1 Myr and 10 Myr, as required to fit the evolution of the CTTS fraction in Taurus, and the evolution of the disc fraction in the clusters studied by Haisch, Lada & Lada (2001).

Although the models do a reasonable job in reproducing the width of the transition seen in Fig. 6, they fail in that they predict a 100% disc fraction for the youngest stars, those with ages less than a Myr. Indeed, from Fig. 5 we find that a disc lifetime as short as (say) 0.5 Myr requires an initial disc mass that is very small – of the order of few $`\times 10^3M_{}`$. As noted previously, one possibility is simply that the ages of the youngest WTTS have been underestimated. We find, however, that the spatial distribution of these young WTTS, relative to the cluster gas, is more consistent with the distribution of the young CTTS (Palla & Stahler 2002) than of the older stars. Another possibility is that these objects, classified as WTTS on the grounds of their equivalent width in $`H_\alpha `$, have highly variable line emission, and have thus been classed as WTTS even though they are not bona fide discless objects. However, we find that these systems are also lacking significant excesses both in the near and far infrared, thus confirming their discless status. These considerations support the contention that the apparently young, discless stars really are young and discless, and constitute a separate population which have lost their discs at a earlier epoch due to processes that do not affect the majority of stars. An ad hoc model in which 30% of stars lose their discs within a Myr, while the rest evolve in the way which we have modelled, is shown in Fig. 6, and is consistent with the available observations.

If stars form in isolation, then it is hard to understand why a significant fraction of them should lose their discs very early. Rapid early accretion onto the star is expected to cease when the disc first becomes gravitationally stable, at a mass of the order of $`0.1M_{}`$, after which time slower ‘viscous’ evolution continues. Within this framework, initial disc masses that are smaller than average by one to two orders of magnitude are puzzling. In an environment where stars are strongly sub-clustered at birth (Scally & Clarke 2002), however, close encounters between stars can efficiently reduce the disc mass, both directly by unbinding disc material (Clarke & Pringle 1993), and indirectly by reducing the outer disc radius and hastening viscous evolution (Armitage & Clarke 1997). Recent simulations of clustered star formation (Bate, Bonnell & Bromm 2003) show that the required interactions may well be extremely common. In the Bate et al. (2003) calculation, most stars have encounters with another star with a minimum separation distance of a few au or smaller. The young WTTS seen in Taurus may simply be the fraction of those stars which fail to accrete a significant new disc following their last close flyby.

### 3.3 Evolution of the accretion rate of CTTS

With the caveats noted above, the evolution of the disc fraction in Taurus is consistent with a theoretical model in which differences in disc lifetime are ascribed to differences in the initial disc mass. A $`1\sigma `$ dispersion of around 0.5 in the $`\mathrm{log}`$ of the initial disc mass is required to reproduce typical disc lifetimes of 1-10 Myr. Figure 7 shows how this model compares with observations of the accretion rate of CTTS in Taurus. We plot estimates of $`\dot{M}`$ by Gullbring et al. (1998) (as tabulated in Hartmann et al. 1998) against our own revised determinations of the stellar ages. The vertical error bars show Hartmann et al.’s (1998) best estimate of the error in determining the accretion rate, which is 0.6 in $`\mathrm{log}(\dot{M})`$. Against this data, we plot theoretical curves for the $`\beta =3/2`$ disc model. We show the predicted evolution of $`\dot{M}`$ for the central estimate of the initial disc mass of $`0.1M_{}`$, and for the models with initial disc masses differing by $`\pm 1\sigma `$.

From the Figure, it is clear that the $`\pm 1\sigma `$ curves bracket the observed range of CTTS accretion rates at a given age. The dispersion in initial disc mass required to reproduce the observed range in disc lifetimes also produces roughly the correct spread in accretion rates. Since we can reproduce the same result with other choices of $`\beta `$, our interpretation is that the Hartmann et al. (1998) plot of $`\dot{M}`$ against $`t`$, shown in a modified form as Figure 7, tells us more about the variation in the initial conditions for discs than it does about the internal physics governing their evolution.

## 4 The influence of binary companions

Most T Tauri stars are found in binary systems (Leinert et al. 1993; Ghez, Neugebauer & Matthews 1993; Simon et al. 1995), which can affect the evolution of circumstellar discs in several ways. First, gravitational torques from the companion act to truncate the disc (Paczynski 1977), and prevent it from expanding freely. Second, if the binary is sufficiently close, it will limit the reservoir of gas available to be accreted. Both of these effects would be expected to reduce the lifetime of discs in binary systems as compared to single stars. A substantially reduced disc lifetime in binaries is not, however, observed (Simon & Prato 1995; Bouvier, Rigaut & Nadeau 1997; Ghez, White & Simon 1997; White & Ghez 2001).

Fig. 8 shows how binaries alter the predicted disc lifetime in our models. We have assumed that the only effect of the binary on the initial surface density distribution is to truncate the disc at a radius $`r_t=d/3`$, where $`d`$ is the binary separation. The initial accretion rate in the binary models is therefore the same as for the isolated disc models. At $`r_t`$, we model the torques from the companion using a $`v_r=0`$ boundary condition. Parameters for the disc wind and the disc viscosity remain unchanged.

By construction, the disc lifetime for sufficiently wide binaries tends to the observed value of $`6\mathrm{Myr}`$. At smaller separations, all three disc models display similar behaviour. The disc lifetime initially increases, reaching a peak at a separation of 15 au to 100 au, depending upon the model. This increase occurs because the disc, which is prevented from expanding freely, has less area from which mass can be lost in a wind. This compensates for, and initially overwhelms, the tendency for smaller disc to have reduced viscous timescales. As a result, only very close binaries ($`d`$ less than 10 au to 40 au, depending upon the model) are predicted to lead to a reduction in the disc lifetime as compared to the lifetime of an isolated disc.

In addition to asking whether the disc lifetime around the primary is a function of binary separation, we can also ask whether differential evolution of the primary and secondary disc is expected to yield ‘mixed’ binary systems that pair CTTS and WTTS. The naive expectation is that the lower mass component of the binary, whose Roche lobe encloses a smaller disc, should lose its disc and become a WTTS first. Observations, however, show that mixed pairs are rather rare. Using a sample of close binaries (most of which had projected separations between 15 au and 100 au) Hartigan & Kenyon (2002) found that only 25 percent (4 out of 16 systems) of binaries containing at least one CTTS were mixed systems. Combining these results with earlier studies (Brandner & Zinnecker 1997; Prato & Simon 1997; Duchene et al. 1999), which obtained similarly small mixed fractions, there is no evidence for a trend in the fraction of mixed systems with separation. Comparing the observations to disc models that do not include winds, this is only consistent with theoretical expectations if the initial conditions yield more massive discs around the lower mass star (Armitage, Clarke & Tout 1999).

If we assume that the initial disc accretion rate, and the mass loss rate, are the same for both stars in a binary, then the disc models presented in this paper predict a low level of mixed WTTS / CTTS binaries. For a mass ratio $`q=0.5`$, the maximum radii of discs surrounding the primary and secondary differ by around 50 percent (Papaloizou & Pringle 1977). From Figure 8, we then estimate that an initially doubly strong (CTTS + CTTS) binary with separation $`d\stackrel{>}{}10`$ au will appear as a mixed binary for at most 10% (if $`\beta =3/2`$) to 20% (if $`\beta =1`$) of the time. A larger fraction (up to 50%) of mixed binaries is predicted for the $`\beta =1/2`$ disc model. As was the case with the disc lifetime, therefore, only very close binaries or binaries with extreme mass ratios are expected to show significant differential disc evolution.

## 5 Summary

We have shown that protoplanetary disc models which include mass loss from the disc at large radius are generally consistent with observations of accretion in T Tauri stars. In particular, we can simultaneously reproduce both the rapid disc dispersal required by observations (Simon & Prato 1995; Wolk & Walter 1996; Duvert et al. 2000), and the near-independence of the disc lifetime on the binary environment (White & Ghez 2001). The evolution of the disc fraction in young clusters, which we have studied in detail for Taurus, can also be modelled within this framework, provided that we allow for a spread in initial disc masses around a central value of $`0.1M_{}`$. We believe that the required level of variation in the initial conditions of discs is a plausible outcome of the turbulent conditions that give rise to star formation within clusters. Simulations of clustered star formation (e.g. Bate, Bonnell & Bromm 2003) show that interactions between small groups of forming stars and their discs are almost inevitable. In such a dynamic environment, we expect large star to star variations in the initial mass and angular momentum content of circumstellar material. In fact, a significant fraction of the stars formed in the Bate et al. (2003) simulation suffer encounters with other stars that are close enough to truncate circumstellar discs to sizes of au or smaller. Such interactions provide a possible origin for the otherwise puzzling observation of a population of apparently very young Weak lined T Tauri stars seen in the Taurus region.

## Acknowledgements

This work has been supported by the European Community’s Research Training Network under contract HPRN-CT-2000-0155, “Young Stellar Clusters”. We thank the referee, Lee Hartmann, for a prompt and valuable review. CJC gratefully acknowledges support from the Leverhulme Trust in the form of a Philip Leverhulme Prize, and from JILA via their visitors’ programme. PJA thanks the Institute of Astronomy for support from their visitors’ programme.


format: text

text:
# Distribution of Transverse Distances in Directed Animals Sumedhasumedha@theory.tifr.res.in and Deepak Dharddhar@theory.tifr.res.in Department Of Theoretical Physics Tata Institute Of Fundamental Research Homi Bhabha Road, Mumbai 400005 India

## Abstract

We relate $`\varphi (\text{x},s)`$, the average number of sites at a transverse distance x in the directed animals with $`s`$ sites in $`d`$ transverse dimensions, to the two-point correlation function of a lattice gas with nearest neighbor exclusion in $`d`$ dimensions. For large $`s`$, $`\varphi (\text{x},s)`$ has the scaling form $`\frac{s}{R_s^d}f(|\text{x}|/R_s)`$, where $`R_s`$ is the root mean square radius of gyration of animals of $`s`$ sites. We determine the exact scaling function for $`d=1`$ to be $`f(r)=\frac{\sqrt{\pi }}{2\sqrt{3}}\text{erfc}(r/\sqrt{3})`$. We also show that $`\varphi (\text{x}=0,s)`$ can be determined in terms of the animals number generating function of the directed animals.

The directed animals(DA) problem describes the large-scale geometrical structure of subcritical directed percolation clusters, and qualitatively models diverse situations such as trees, river networks and dilute polymers in a flowing solvent . Formally, it corresponds to $`p0`$ limit of the directed percolation problem . The imposed directionality changes the universality class, and the large-scale structure of directed animals is different from isotropic animals. The enumeration of directed site animals in $`d+1`$ dimensions is related to hard-core lattice gas(HCLG) at negative activity with repulsive interactions in $`d`$ dimensions and the Yang-Lee edge problem in $`d`$ dimensions . This equivalence helps to obtain several exact results for the directed animals problem . In particular, the critical exponents $`\theta `$ and the transverse size exponent $`\nu _{}`$ are known exactly in $`d=0,1`$ and $`2`$ . For reviews of available exact results, see . The problem is also related to a recently studied model of quantum gravity , and there is an unexpected relation between the number of distinct eigenvalues for Potts model partition function on strips of width $`w`$ and the number of directed animals with $`w`$ sites in two dimensions .

In this paper, we extend the known relation between the directed site animals enumeration (DSAE) problem to hard-core lattice gases at negative activity to obtain the average number of sites at a given transverse distance x from the origin for $`d+1`$-dimensional directed animals from the density-density correlation function of the lattice gas in $`d`$ dimensions. For $`d=1`$, using the exact generating function for this correlation function, we determine the scaling form for the average number of sites at a given transverse distance in a $`2`$-dimensional directed animal having $`s`$ sites, for large $`s`$. In this case, the average transverse size of the animal scales as $`s^\nu _{}`$, where $`\nu _{}`$ is known to be $`\frac{1}{2}`$ for $`d=1`$. For large $`s`$, the average value of $`q`$-th transverse moment $`<|x|^q>`$ varies as $`R_s^qC_q`$ , where $`R_s`$ is the root mean square radius of gyration of animals of $`s`$ sites. Using the exact scaling function we are able to determine the universal constants $`C_q`$ for all $`q`$.

We consider directed animals on a square lattice. Each site $`(x,t)`$ has two bonds directed outwards towards the sites $`(x+1,t+1)`$ and $`(x1,t+1)`$ (Fig. 1). A directed site animal $`𝒜`$ rooted at the origin is a set of occupied sites including origin, such that for each occupied site $`(x,t)`$ other than the origin, at least one of the two sites $`(x1,t1)`$ and $`(x+1,t1)`$ is also occupied. Fig. 1 shows a directed animal on a square lattice in two dimensions. The number of sites in $`𝒜`$ will be denoted by $`|𝒜|`$ or $`s`$. We define $`n(x|𝒜)`$ as the number of sites of $`𝒜`$ having the transverse coordinate $`x`$. In Fig. 1, $`n(x|𝒜)`$ takes the values $`1,3`$ and $`0`$ for $`x`$ equal to $`4,2,5`$ respectively, with $`|𝒜|=20`$. We shall define the radius of gyration of $`𝒜`$ as $`\sqrt{\frac{1}{s}x^2}`$, where the sum is over all sites of the animal. The animal shown in Fig. 1, has a squared radius of gyration $`21/4`$.

We define $`A(y)`$ as $`y^{|𝒜|}`$ summed over all animals $`𝒜`$ rooted at the origin. If $`A_s`$ is the number of distinct animals having $`s`$ sites, it is easy to see that

$$A(y)=\underset{𝒜}{}y^{|𝒜|}=\underset{s=1}{\overset{\mathrm{}}{}}A_sy^s$$
(1)
For large $`s`$, $`A_s`$ varies as $`\lambda ^ss^\theta `$, where $`\theta `$ is a critical exponent. The radius of gyration $`R_s`$ is expected to vary as $`s^\nu _{}`$, where the exponent $`\nu _{}`$ is related to the animals number exponent $`\theta `$ by the hyperscaling relation $`\theta =d\nu _{}`$ .

For a directed animal of size $`s`$ we define the generating function $`\mathrm{\Psi }(\text{x},y)`$ as

$$\mathrm{\Psi }(\text{x};y)=\underset{𝒜}{}n(\text{x}|𝒜)y^{|𝒜|}$$
(2)
where the summation over $`𝒜`$ is the summation over all animal configurations. This can be written as

$$\mathrm{\Psi }(\text{x};y)=\underset{s}{}\varphi (\text{x},s)A_sy^s$$
(3)
where $`\varphi (\text{x},s)`$ is the value of $`n(\text{x},𝒜)`$ averaged over all animals $`𝒜`$ of size $`s`$.

The DSAE problem in $`d+1`$-dimensions is related to time development of thermal relaxation of a hard core lattice gas (HCLG) with nearest neighbor exclusion on $`d`$ dimensional lattice . On a $`d`$ dimensional body-centered hypercubical lattice (a $`d`$ dimensional body-centered hypercubical lattice is a hypercubical lattice having $`2^{d1}`$ bonds going outward from each site), the dynamics of the lattice gas is defined as follows: The evolution is stochastic discrete -time Markovian. At odd(even) times, all the odd(even) sites are examined in parallel, and if a site has all neighbors empty, its occupation number is set to $`1`$ with probability $`p`$, and to $`0`$ with probability $`(1p)`$. The rates of transitions satisfy detailed balance condition corresponding to the Hamiltonian

$$H=+\mathrm{}\underset{<ij>}{}n_in_j(\text{ln}z)\underset{i}{}n_i$$
(4)
with $`z=p/(1p)`$.

Let $`\rho (p)`$ be the average density of particles in the steady state of this system. In , it was shown that we have

$$A(y)=\rho (p=y)$$
(5)
For the $`d=1`$ nearest-neighbor-exclusion lattice gas, it is straight forward to determine the average density corresponding to chemical activity $`p/(1p)`$. This for square lattice gives

$$A(y)=\frac{1}{2}\left[1\sqrt{\left(\frac{1+y}{13y}\right)}\right]$$
(6)
The animal numbers $`A_s`$, which are the coefficients in the Taylor expansion of $`A(y)`$ for square lattice can be written as

$$A_s=_0^{2\pi }\frac{d\theta }{2\pi }(1+cos\theta )(1+2cos\theta )^{s1}.$$
(7)
For large $`s`$, $`A_s\frac{1}{\sqrt{3\pi }}3^ss^{\frac{1}{2}}`$, which shows that in this case $`\lambda =3`$, and $`\theta =1/2`$.

The derivation of Eq.(5) is easily generalized to the case where the values of $`p`$ at different sites are different. Let the probability that site $`i`$ gets occupied at time $`t+1`$ given that all its neighbors are unoccupied at time $`t`$ be $`p_i`$. Then the rates of this process still satisfy detailed balance condition corresponding to the hamiltonian

$$H=+\mathrm{}\underset{<ij>}{}n_in_j\underset{i}{}(\text{ln}z_i)n_i$$
(8)
where $`z_i=p_i/(1p_i)`$. The probability that site $`i`$ is occupied in the steady state depends on the $`p_j`$’s for all sites $`j`$, and will be denoted by $`\rho _i(\{p_j\})`$. In the corresponding DA problem, we have to define the weight of an animal $`𝒜`$ as product of weights of all occupied sites, the weight corresponding to a site with x-coordinate $`j`$ being $`y_j`$. Then define $`A_i(\{y_j\})`$ be the sum of weights of all animals rooted at $`i`$. Then clearly $`A_i(\{y_j\})`$ is a formal power series in the variables $`\{y_j\}`$. If all $`y_j=y`$, this becomes independent of $`i`$, and reduces to the function $`A(y)`$. For unequal $`y_j`$’s, Eq.(5) becomes

$$A_i(\{y_j\})=\rho _i(\{p_j=y_j\})$$
(9)
Applying the operator $`y_𝐱\frac{}{y_𝐱}`$ on the weight of any particular animal gives us the weight multiplied by the number of occupied sites with transverse coordinate x in the animal. Thus, clearly, we get

$$y_𝐱\frac{}{y_𝐱}A_0(\{y_j\})|_{\{y_j=y\}}=\mathrm{\Psi }(\text{x};y)$$
(10)
Let $`\mathrm{\Omega }(\{z_j\})`$ be the grand partition function for the HCLG given by the Hamiltonian in Eq. (8) with $`z_j=p_j/(1p_j)`$. Then, $`\mathrm{\Omega }(\{z_j\})`$ is a linear function of each of the variables $`z_j`$. Let $`\eta _j`$ be the indicator variable taking value $`1`$ if the site $`j`$ is occupied, and $`0`$ if not. The density-density correlation function of the gas $`G(𝐢,𝐤)`$ is defined as

$$G(𝐢,𝐤;\{z_j\})=<\eta _𝐤\eta _𝐢><\eta _𝐤><\eta _𝐢>=z_𝐤\frac{}{z_𝐤}\rho _𝐢(\{z_j\}).$$
(11)
When $`z_j=z`$ for all $`j`$, the correlation function depends only on $`(𝐢𝐤)`$ and hence can be written as $`G(𝐢𝐤;z)`$.

Now, using Eq.(9), we get

$$\mathrm{\Psi }(\text{x};y)=\frac{1}{1+y}G\left(\text{x};z=\frac{y}{1+y}\right)$$
(12)
This equation is the main result of this paper, and equates the pair correlation function $`G(\text{x};z)`$ of the HCLG with the function $`\mathrm{\Psi }(\text{x};y)`$ which gives the density profile of the DA problem.

For the special case $`d=1`$, it is a simple excercise to calculate $`G(x;z)`$ explicitly, using transfer matrix methods . This gives

$$G(x;z)=\frac{z}{1+4z}\left[\frac{1\sqrt{14z}}{1+\sqrt{1+4z}}\right]^{|x|}$$
(13)
and we get the explicit expression for $`\mathrm{\Psi }(x;y)`$ on square lattice to be

$$\mathrm{\Psi }(x;y)=\frac{y}{(1+y)(13y)}\left[1\sqrt{\frac{13y}{1+y}}\right]^{|x|}\left[1+\sqrt{\frac{13y}{1+y}}\right]^{|x|}$$
(14)
This determines the density profile in the constant fugacity ensemble, where an animal having $`s`$ sites has weight $`y^s`$. However, it is more instructive to look at the profile in the constant-$`s`$ ensemble. This is obtained by looking at the Taylor coefficient of $`y^s`$ in the above equation. This can become rather messy. However, the behavior for large $`s`$ can be determined easily.

In general, for any dimension $`d`$ , for large $`s`$, $`\varphi (\text{x},s)`$ has the scaling form

$$\varphi (\text{x},s)\frac{s}{R_{s}^{}{}_{}{}^{d}}f\left(\frac{|\text{x}|}{R_s}\right)$$
(15)
Usually the argument of the scaling function is defined only upto a multiplicative constant. We have made specific choice for this by using the variable as $`|\text{x}|/R_s`$. The normalization of scaling function $`f(r)`$ is chosen such that is satisfies

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}d^d\text{x}f(|\text{x}|)=1`$ (16)
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}d^d\text{x}|\text{x}|^2f(|\text{x}|)=1`$ (17)
In the DA problem, as the number of animals grow as $`\lambda ^s`$, the series expansions for $`\mathrm{\Psi }(\text{x};y)`$ or $`A(y)`$ in powers of $`y`$ converge for $`y<y_c=1/\lambda `$. For $`y`$ near $`y_c`$, the singular part of the function $`A(y)`$ varies as $`(1y\lambda )^{\theta 1}`$. For the HCLG problem, this corresponds to a singularity in the series for the density $`\rho `$ in powers of activity for the activity $`z_{LY}=(1+\lambda )^1`$ \[Eq. (12)\]. This singularity on the negative real line is the Lee-Yang edge singularity for this problem.

For $`z`$ near the critical value $`z=z_{LY}`$, this correlation function $`G(\text{x};z)`$ is expected to have the scaling form

$$G(\text{x};z=z_{LY}e^ϵ)=cϵ^ag(b|\text{x}|ϵ^\nu )+\mathrm{higher}\mathrm{order}\mathrm{terms}\mathrm{in}ϵ$$
(18)
where $`b`$ and $`c`$ are non-universal, lattice dependent constants. We choose $`b`$ such that $`g(\xi )=exp(\xi )`$ for large $`\xi `$, and $`c`$ is fixed by requiring $`g(0)=1`$.

The scaling function $`g(\xi )`$ tends to a constant limiting value as $`\xi `$ tends to zero, and decreases to zero exponentailly fast as $`\xi `$ tends to infinity. Let the power-series expansion of $`g(\xi )`$ about $`\xi =0`$ be given by

$$g(\xi )=\underset{k=0}{\overset{\mathrm{}}{}}g_k\xi ^k$$
(19)
Substituting this in Eq.(18), we get

$$G(\text{x};z)\underset{k=0}{\overset{\mathrm{}}{}}g_kb^k|\text{x}|^k\left(1\frac{z}{z_{LY}}\right)^{\nu ka}$$
(20)
To get $`\varphi (\text{x},s)`$, we need to determine the coefficient of $`z^s`$ in the above expansion. From the binomial expansion of $`(1\frac{z}{z_{LY}})^{\nu ka}`$ we immediately get

$$G(\text{x};z)\underset{s=0}{\overset{\mathrm{}}{}}\left(\frac{z}{z_{LY}}\right)^s\underset{k=0}{\overset{\mathrm{}}{}}g_kb^k|\text{x}|^k\frac{\mathrm{\Gamma }(a+s\nu k)}{\mathrm{\Gamma }(s+1)\mathrm{\Gamma }(a\nu k)}$$
(21)
For fixed $`k`$, for large $`s`$, we have

$$\frac{\mathrm{\Gamma }(s+a\nu k)}{\mathrm{\Gamma }(s+1)}s^{a1\nu k}$$
(22)
Using this in Eq.(21), we get

$$f(|\text{x}|s^\nu _{})=\underset{k=0}{\overset{\mathrm{}}{}}g_k\left(\frac{|\text{x}|}{s^\nu }\right)^k\frac{1}{\mathrm{\Gamma }(a\nu k)}$$
(23)
where the correlation length exponent $`\nu `$ for the HCLG problem is the same as the transverse size exponent $`\nu _{}`$ for directed animals.

For $`d=1`$, it is easy to determine $`G(x;z)`$ explicitly. In this case, using Eq.(13) we get the scaling form for $`G(x;z)`$ as

$$G(x;ϵ)=\frac{1}{4ϵ}\text{exp}(2x\sqrt{ϵ})$$
(24)
so that $`a=1`$, $`\nu =1/2`$, $`b=2`$ and $`c=1/4`$. The scaling function $`g(\xi )`$ is simply given by

$$g(\xi )=\text{exp}[\xi ]$$
(25)
hence $`g_k=1/\mathrm{\Gamma }[k+1]`$. Using this, and the values of $`a`$ and $`\nu `$ in Eq.(23) we get the leading singular behavior of $`\varphi (x,s)A_s`$ on square lattice to be

$$\varphi (x,s)A_s=\frac{1}{4}y_c^s\underset{k=0}{\overset{\mathrm{}}{}}\frac{(\frac{\sqrt{3}x}{\sqrt{s}})^k}{\mathrm{\Gamma }(k+1)}\frac{1}{\mathrm{\Gamma }(1\frac{k}{2})}$$
(26)
Since $`\mathrm{\Gamma }(1\frac{k}{2})`$ has poles when $`k`$ is an even integer, only the odd terms contribute to the sum. The resulting series is easy to sum explicitly, giving

$$\varphi (x,s)A_s=\frac{1}{4}3^s\text{erfc}\left(\frac{\sqrt{3}x}{2\sqrt{s}}\right)+correctiontoscalingterms$$
(27)
where $`\text{erfc}(x)=\frac{2}{\sqrt{\pi }}_x^{\mathrm{}}e^{x^2}𝑑x`$.

This gives $`\varphi (x,s)`$ for large $`s`$ to be

$$\varphi (x,s)=\frac{\sqrt{3\pi s}}{4}\text{erfc}\left(\frac{\sqrt{3}x}{2\sqrt{s}}\right)+correctiontoscalingterms$$
(28)
From $`\varphi (x,s)`$ we can as well derive the expression for $`q`$th transverse moment of the directed animals. The $`q`$th transverse moment of a cluster of size $`s`$ denoted by $`\mu _{qs}`$ is defined as

$$\mu _{qs}=\underset{i=s}{\overset{s}{}}\varphi (i,s)|i|^q$$
(29)
Using the scaling form of $`\varphi (i,s)`$, for large $`s`$, we see that $`(\mu _{qs}/sR_s^q)`$ is a universal constant. Denoting it by $`C_q`$, we get

$$\frac{\mu _{qs}}{sR_s^q}=C_q=_{\mathrm{}}^{\mathrm{}}d^d\text{x}|\text{x}|^qf(|\text{x}|)$$
(30)
For $`d=1`$, we get

$$C_q=\frac{1}{q+1}(3)^{\frac{q}{2}}\mathrm{\Gamma }\left[1+\frac{q}{2}\right]$$
(31)
and $`R_s=\frac{2}{3}\sqrt{s}`$. Hence we get the scaling function(16) in $`d=1`$ to be

$$f(r)=\frac{\sqrt{\pi }}{2\sqrt{3}}\text{erfc}\left(\frac{r}{\sqrt{3}}\right)$$
(32)
In the entire low-density phase of the HCLG, in any dimension $`d`$, we expect the correlation function $`G(\text{x};z)`$ to have an exponential decay at large $`|\text{x}|`$. But the behavior of the scaling function $`f(r)`$ for large $`r`$ is in general different. Suppose $`\mathrm{ln}f(r)`$ varies as $`r^\alpha `$ for large $`r`$. Putting this behavior in Eq.(3), and using $`A_s\lambda ^s`$ and Eq. (12) , we get

$$G(\text{x};z_{LY}e^ϵ)\underset{s}{}\text{exp}(|\text{x}|^\alpha s^{\nu \alpha }ϵs)$$
(33)
For large $`|\text{x}|`$, the integral can be estimated by steepest descent, and gives $`\mathrm{log}g(\xi )`$ varying as $`\xi ^{\frac{\alpha }{1+\nu \alpha }}`$. Since this should be linear in $`\xi `$, we see that $`\alpha =\frac{1}{1\nu }`$. As a check, we see that in $`d=1`$, $`\nu =1/2`$, and $`f(r)`$ varies as $`\text{exp}(r^2)`$ for large $`r`$.

The case $`x=0`$ is special, in that the density-density correlation function $`G(0,z)`$ is always equal to $`\rho (1\rho )`$ for hard-core lattice gas for any $`d`$ dimensional lattice. Hence, if one knows $`\rho `$ as a function of the activity $`z`$ ( equivalently, in the DA problem, one knows the animal numbers generating function $`A(y)`$), then one can determine $`\mathrm{\Psi }(0;y)`$ in terms of $`A(y)`$ alone. For a $`d+1`$ dimensional DA on a body-centered hypercubical lattice it is given by

$$\mathrm{\Psi }(0;y)=\frac{1}{1+y}A(y)[1+A(y)]$$
(34)
For $`d=1`$, using the $`A(y)`$ from Eq. (6), we get for the square lattice DA problem

$$\mathrm{\Psi }(0;y)=\frac{y}{12y3y^2}$$
(35)
Expanding in powers of $`y`$ we get

$$\varphi (0,s)=\frac{3^s+(1)^{s1}}{4A_s}$$
(36)
where $`A_s`$ is as given in Eq.(7). For large $`s`$, $`\varphi (0,s)`$ varies as $`s^{1/2}`$ as expected.

Similar analysis can be extended to higher dimensions. In general the scaling function $`f(r)`$ tends to a finite value as $`r`$ tends to $`0`$, and hence $`\varphi (0,s)`$ varies as $`s^{1\theta }`$.

For $`d>7`$, mean-field theory becomes asymptotically exact  with $`R_ss^{1/4}`$ and the scaling form of $`f(r)`$ is

$$f(r)=\frac{2}{\mathrm{\Omega }_d}\frac{1}{r^{d2}}\text{exp}(r^2)$$
(37)
where $`\mathrm{\Omega }_d`$ is the surface area of a $`d`$ dimensional unit sphere and is equal to $`\mathrm{\Omega }_d=\frac{2\pi ^{d/2}}{\mathrm{\Gamma }(d/2)}`$.


text:
# Polarized X-Ray Emission from Magnetized Neutron Stars: Signature of Strong-Field Vacuum Polarization

## Abstract

In the atmospheric plasma of a strongly magnetized neutron star, vacuum polarization can induce a MSW-like resonance across which a X-ray photon may (depending on its energy) convert from one mode into the other, with significant changes in opacities and polarizations. We show that this vacuum resonance effect gives rise to an unique energy-dependent polarization signature in the surface emission from neutron stars: for “normal” field strengths ($`10^{12}B7\times 10^{13}`$ G), the plane of linear polarization at photon energy $`E1`$ keV is perpendicular to that at $`E4`$ keV, while for “superstrong” field strengths ($`B7\times 10^{13}`$ G), the polarization planes at different energies coincide. The detection of polarized X-rays from neutron stars can provide a direct probe of strong-field quantum electrodynamics and constrain the neutron star magnetic field and geometry.

Surface emission from isolated neutron stars (radio pulsars and radio quiet neutron stars, including magnetars thompsonduncan95 ) provides an useful probe of the neutron star (NS) interior physics, surface magnetic fields, and composition. The advent of X-ray telescopes in recent years has made the detection and detailed study of NS surface emission a realitypavlovetal02 . On the other hand, X-ray emission from accreting X-ray pulsars has long yielded information on the dynamics of accretion and radiative mechanisms in the strong gravity, strong magnetic field regime meszaros92 .

In the magnetized plasma that characterizes NS atmospheres, X-ray photons propagate in two normal modes: the ordinary mode (O-mode) is mostly polarized parallel to the $`𝐤`$-$`𝐁`$ plane, while the extraordinary mode (X-mode) is mostly polarized perpendicular to the $`𝐤`$-$`𝐁`$ plane, where $`𝐤`$ is the photon wave vector and $`𝐁`$ is the external magnetic field meszaros92 . This description of normal modes applies under typical conditions, when the photon energy $`E`$ is much less than the electron cyclotron energy $`E_{Be}=\mathrm{}eB/(m_ec)=11.6B_{12}`$ keV \[where $`B_{12}=B/(10^{12}\mathrm{G})`$\], $`E`$ is not too close to the ion cyclotron energy $`E_{Bi}=6.3B_{12}(Z/A)`$ eV (where $`Z`$ and $`A`$ are the charge number and mass number of the ion), the plasma density is not too close to the vacuum resonance (see below) and $`\theta _{kB}`$ (the angle between $`𝐤`$ and $`𝐁`$) is not close to zero. Under these conditions, the X-mode opacity (due to scattering and absorption) is greatly suppressed compared to the O-mode opacity, $`\kappa _X(E/E_{Be})^2\kappa _O`$ lodenquaietal74 ; meszaros92 . As a result, the X-mode photons escape from deeper, hotter layers of the NS atmosphere than the O-mode photons, and the emergent radiation is linearly polarized to a high degree (as high as $`100\%`$) gnedinsunyaev74 ; meszarosetal88 ; pavlovzavlin00 . Measurements of X-ray polarization, particularly when phase-resolved and measured in different energy bands, could provide unique constraints on the magnetic field strength and geometry and the compactness of the NS meszarosetal88 ; pavlovzavlin00 ; heyletal03 .

It has long been predicted from quantum electrodynamics (QED) that in a strong magnetic field the vacuum becomes birefringent heisenbergeuler36 ; schwinger51 ; adler71 ; tsaierber75 ; heylhernquist97 . Acting by itself, the birefringence from vacuum polarization is significant (with the index of refraction differing from unity by more than $`10\%`$) only for $`B300B_Q`$, where $`B_Q=m_e^2c^3/(e\mathrm{})=4.414\times 10^{13}`$ G is the critical QED field strength. However, when combined with the birefringence due to the magnetized plasma, vacuum polarization can greatly affect radiative transfer at much smaller field strengths. A “vacuum resonance” arises when the contributions from the plasma and vacuum polarization to the dielectric tensor “compensate” each other gnedinetal78 ; meszarosventura79 ; pavlovgnedin84 ; laiho02 . For a photon of energy $`E`$, the vacuum resonance occurs at the density $`\rho _V9.64\times 10^5Y_e^1B_{12}^2E_1^2f^2\mathrm{g}\mathrm{cm}^3`$, where $`Y_e`$ is the electron fraction, $`E_1=E/(1\mathrm{keV})`$, and $`f=f(B)`$ is a slowly varying function of $`B`$ and is of order unity ($`f=1`$ for $`BB_Q`$ and $`f(B/5B_Q)^{1/2}`$ for $`BB_Q`$; see refs. laiho02 ; holai03 ). For $`\rho >\rho _V`$ (where the plasma effect dominates the dielectric tensor) and $`\rho <\rho _V`$ (where vacuum polarization dominates), the photon modes (for $`EE_{Be}`$, $`EE_{Bi}`$ and $`\theta _{kB}0`$) are almost linearly polarized; near $`\rho =\rho _V`$, however, the normal modes become circularly polarized as a result of the “cancellation” of the plasma and vacuum effects — both effects tend to make the mode linearly polarized, but in mutually orthogonal directions. When a photon propagates in an inhomogeneous medium, its polarization state will evolve adiabatically (i.e. following the $`K_+`$ or $`K_{}`$ curve in Fig. 1) if the density variation is sufficiently gentle. Thus, a X-mode (O-mode) photon will be converted into the O-mode (X-mode) as it traverses the vacuum resonance, with its polarization ellipse rotated by $`90^{}`$ (Fig. 1). This resonant mode conversion is analogous to the Mikheyev-Smirnov-Wolfenstein neutrino oscillation that takes place in the Sunhaxton95 ; bahcall03 . For this conversion to be effective, the adiabatic condition must be satisfied laiho02 ; laiho03  (cf. pavlovgnedin84 )

$$EE_{\mathrm{ad}}=2.52\left(f\mathrm{tan}\theta _{kB}|1u_i|\right)^{2/3}\left(\frac{1\mathrm{cm}}{H_\rho }\right)^{1/3}\mathrm{keV},$$
(1)
where $`u_i=(E_{Bi}/E)^2`$ and $`H_\rho =|dz/d\mathrm{ln}\rho |`$ is the density scale height (evaluated at $`\rho =\rho _V`$) along the ray. For an ionized Hydrogen atmosphere, $`H_\rho 2kT/(m_pg\mathrm{cos}\theta )=1.65T_6/(g_{14}\mathrm{cos}\theta )`$ cm, where $`T=10^6T_6`$ K is the temperature, $`g=10^{14}g_{14}`$ cm s-2 is the gravitational acceleration, and $`\theta `$ is the angle between the ray and the surface normal. In general, the mode conversion probability is given by laiho02 ; laiho03

$$P_{\mathrm{con}}=1\mathrm{exp}\left[(\pi /2)(E/E_{\mathrm{ad}})^3\right].$$
(2)
The probability for a nonadiabatic “jump” is $`(1P_{\mathrm{con}})`$.

Because the two photon modes have vastly different opacities, the vacuum resonance can significantly affect the transfer of photons in NS atmospheres. When the vacuum polarization effect is neglected, the decoupling densities of the O-mode and X-mode photons (i.e., the densities of their respective photospheres, where the optical depth measured from outside is $`2/3`$) are approximately given by (for Hydrogen plasma and $`\theta _{kB}`$ not too close to $`0`$) $`\rho _O0.42T_6^{1/4}E_1^{3/2}G^{1/2}`$ g cm-3 and $`\rho _X486T_6^{1/4}E_1^{1/2}B_{14}G^{1/2}`$ g cm-3, where $`G=1e^{E/kT}`$ laiho02 . There are two different magnetic field regimes: For “normal” magnetic fields,

$$B<B_l6.6\times 10^{13}T_6^{1/8}E_1^{1/4}G^{1/4}\mathrm{G},$$
(3)
the vacuum resonance lies outside both photospheres ($`\rho _V<\rho _O,\rho _X`$); for the magnetar field regime, $`B>B_l`$, the vacuum resonance lies between these two photospheres, i.e., $`\rho _O<\rho _V<\rho _X`$ (the condition $`\rho _V<\rho _X`$ is satisfied for all field strengths and relevant energies and temperatures). These two field regimes yield qualitatively different X-ray polarization signals.

Consider the “normal” field strengths, $`10^{12}\mathrm{G}BB_l`$, which apply to most NSs (see Fig. 2). In this regime, the atmosphere structure and total spectrum can be calculated without including vacuum polarization \[to be more accurate, we require $`(B/B_l)^41`$ for this to be valid\]. For concreteness, we consider emission from a hot spot (magnetic polar cap) on the NS; the magnetic field at the hot spot (with size much smaller than the stellar radius) is perpendicular to the stellar surface. Let the specific intensities of the O-mode and X-mode emerging from their respective photospheres (which lie below the vacuum resonance) be $`I_O^{(0)}`$ and $`I_X^{(0)}`$, which we calculate using our H atmosphere models developed previously holai01 ; holai03 . For a given $`B`$ and $`T_{\mathrm{eff}}`$, both $`I_O^{(0)}`$ and $`I_X^{(0)}`$ depend on $`E`$ and $`\theta _{kB}`$ at emission (the hot spot). As the the radiation crosses the vacuum resonance, the intensities of the O-mode and X-mode become $`I_O=(1P_{\mathrm{con}})I_O^{(0)}+P_{\mathrm{con}}I_X^{(0)}`$ and $`I_X=(1P_{\mathrm{con}})I_X^{(0)}+P_{\mathrm{con}}I_O^{(0)}`$. In calculating $`P_{\mathrm{con}}`$, we use the temperature profile of the atmosphere model to determine the density scale height at the vacuum resonance. We note that in principle, circular polarization can be produced when a photon crosses the vacuum resonance laiho03 , but the net circular polarization is expected to be zero when photons from a finite-sized polar cap are taken into account.

To determine the observed polarization, we must consider propagation of polarized radiation in the NS magnetosphere, whose dielectric property in the X-ray band is dominated by vacuum polarization heyletal03 . As a photon (of a given mode) propagates from the polar cap through the magnetosphere, its polarization state evolves adiabatically following the varying magnetic field it experiences, up to the “polarization-limiting radius” $`r_{\mathrm{pl}}`$, beyond which the polarization state is frozen. We set up a fixed coordinate system $`XYZ`$, where the $`Z`$-axis is along the line-of-sight and the $`X`$-axis lies in the plane spanned by the $`Z`$-axis and $`𝛀`$ (the spin angular velocity vector). The direction of the magnetic field $`𝐁`$ following the photon trajectory in this coordinate system is specified by the polar angle $`\theta _B(s)`$ and the azimuthal angle $`\varphi _B(s)`$ (where $`s`$ is the affine parameter along the ray). The $`X,Y`$-components of the electric field of the photon mode are $`(E_X,E_Y)=(\mathrm{cos}\varphi _B,\mathrm{sin}\varphi _B)`$ for the O-mode and $`(\mathrm{sin}\varphi _B,\mathrm{cos}\varphi _B)`$ for the X-mode. The adiabatic condition requires $`|n_Xn_O|2(\mathrm{}c/E)|d\varphi _B/ds|`$ (where $`n_X,n_O`$ are the indices of refraction of the two modes); this condition breaks down at $`r=r_{\mathrm{pl}}`$, which is much greater than the stellar radius for all parameter regimes of interest in this paper. The observed Stokes parameters (normalized to the total intensity $`I`$) are given by $`Q/I=(2P_{\mathrm{con}}1)p_e\mathrm{cos}2\varphi _B(r_{\mathrm{pl}})`$ and $`U/I=(2P_{\mathrm{con}}1)p_e\mathrm{sin}2\varphi _B(r_{\mathrm{pl}})`$, where $`p_e=(I_X^{(0)}I_O^{(0)})/(I_X^{(0)}+I_O^{(0)})`$ is the “intrinsic” polarization fraction at emission, and $`\varphi _B(r_{\mathrm{pl}})`$ is the azimuthal angle of the magnetic field at $`r=r_{\mathrm{pl}}`$ \[For $`r_{\mathrm{pl}}c/\mathrm{\Omega }`$, or for spin frequency $`70`$ Hz, $`\varphi _B(r_{\mathrm{pl}})\pi +\varphi _B(R)`$\]. We calculate the observed spectral fluxes $`F_I,F_Q,F_U`$ (associated with the intensities $`I,Q,U`$) using the standard procedure, including the effect of general relativity pechenicketal83 .

Figure 3 shows the total flux and polarization “light curves” generated from the hot polar cap of a rotating NS. As the star rotates, the plane of linear polarization rotates. The light curves obviously depend on the geometry (specified by the angles $`\beta `$ and $`\gamma `$). Note that $`F_Q`$ for low-energy ($`E1`$ keV) photons is opposite to that for high-energy photons ($`E4`$ keV), which implies that the planes of linear polarization at low and high energies are orthogonal. This feature also manifests in the phase-averaged linear polarization (see Fig. 4). This is an unique signature of photon mode conversion induced by vacuum polarization.

In the magnetar field regime, $`B>B_l`$, the vacuum resonance lies between the photospheres of the two modes. At low energies (e.g., $`E1`$ keV), no mode conversion occurs at the vacuum resonance, and vacuum polarization makes the X-mode photosphere lie above (i.e. at lower density) its “original” location (i.e. when vacuum polarization is turned off) because the photon opacity exhibits a large spike at the vacuum resonance laiho02 . At high energies (e.g., $`E4`$ keV), the effective X-mode photosphere lies very near the vacuum resonance. For both low and high energies, the emergent radiation is dominated by the X-mode (for $`\theta _{kB}`$ not too close to 0). Therefore the planes of linear polarization at different energies coincide and evolve “in phase” as the star rotates (see Fig. 5). These polarization signals are qualitatively different from the “normal” field regime. Note that in the magnetar field regime, vacuum polarization significantly affects the total spectral flux from the atmosphere laiho02 ; holai03 ; laiho03 : (i) Vacuum polarization makes the effective decoupling density of X-mode photons (which carry the bulk of the thermal energy) smaller, thereby depleting the high-energy tail of the spectrum and making the spectrum closer to black-body; (ii) vacuum polarization suppresses the proton cyclotron line and other spectral line features hoetal03  in the spectra by making the decoupling depths inside and outside the line similar. As we discussed in previous papers holai03 ; hoetal03 , the absence of lines in the observed thermal spectra of several magnetar candidates juettetal02  may be considered as a proof of the vacuum polarization effect at work in these systems. Measurement of X-ray polarization would provide an independent probe of the magnetic fields of these objects.

Finally, we note that although the specific results shown in Figs. 3-5 refer to emission from a hot polar cap on the NS, we expect the vacuum polarization signature (e.g., that the planes of linear polarization at $`E1`$ keV and at $`E4`$ keV are orthogonal for $`B7\times 10^{13}`$ G) to be present in more complicated models (e.g. when several hot spots or the whole stellar surface contribute to the X-ray emission). This is because the polarization-limiting radius (due to vacuum polarization in the magnetosphere) lies far away from the star (see above), where rays originating from different patches of the star experience the same dipole field heyletal03 . Our results therefore demonstrate the unique potential of X-ray polarimetrycostaetal01  in probing the physics under extreme conditions (strong gravity and magnetic fields) and the nature of various forms of NSs.

This work was supported by NASA grant NAG 5-12034 and NSF grant AST 0307252.


text:
# Untitled Document

Unfortunate Terminology

Jerzy Klemens Kowalczyński

Mailing address: Institute of Physics, Polish Academy of Sciences,
Al. Lotników 32/46, 02-668 Warsaw, Poland

E-mail address: jkowal@ifpan.edu.pl

Abstract. The phrase “negative squared rest mass” can sometimes be found in papers on neutrinos and frequently occurs in the tachyonic literature. Consequently, some authors say that “the rest mass of tachyons is imaginary”. Besides, the statement “photons have zero rest mass” is almost common. In terms of relativity, however, the state of rest cannot be reasonably defined for luxons and tachyons, and, therefore, today it does not make sense to speak of rest mass of such objects. It is shown here that the phrases “negative squared mass”, “imaginary mass”, and “photon’s zero mass” result from applying bradyonic dynamical relativistic relations to determine properties of luxons and tachyons; and that this erroneous procedure results from an unfortunate interpretation of kinematical relativistic relations. It is also shown that the use of proper relativistic relations, i.e. luxonic or tachyonic ones, gives a positive quantity having the squared mass dimension. Thus we obtain a nonzero real quantity having the mass dimension (called masslike quantity), which is positive for the photon and may (by intuition – should) be positive for other luxons and for tachyons.

From some papers on measurements of the neutrino mass one can learn that “the squared rest mass of neutrinos is negative”. In the tachyonic literature it is frequently stated that “the squared rest mass of tachyons is negative”, and, consequently, some authors conclude that “the rest mass of tachyons is imaginary”, though in the case of neutrinos I have not met such a heroic author. Besides, the statement “photons have zero rest mass” is almost common. These statements are said to be conclusions from relativity but this is not true.

In relativity the term “rest mass” does not make sense in the case of luxons and tachyons, since the state of rest can be reasonably defined for these objects neither within standard relativity nor in its consistent extensions. This is obvious in the luxonic case since, e.g., the Lorentz transformation is singular for speeds equal to $`c`$. If we were to assume that any tachyon may be at rest, then three independent states of its rest would have to exist since observers who see our spacetime as that having three time dimensions and only one space dimension would then have to be admitted.

As regards the phrases “negative squared mass”, “imaginary mass”, and “photon’s zero mass”, we shall proceed step by step.

Consider the world line $`x^\mu \left(\sigma \right)`$ of a pointlike object. Assume, for simplicity, that the object is free in flat spacetime endowed with the Lorentzian coordinates (i.e. $`x^\mu \left(\sigma \right)`$ is straight), that $`\sigma `$ is the normalized affine parameter of $`x^\mu \left(\sigma \right)`$, and that the signature is, e.g., $`+++`$. Note that in the metric form expressions, $`ds^2=dx_\mu dx^\mu ,ds^2`$ is only a conventional symbol, and therefore it need not be the square of an infinitesimal real quantity. In the case under consideration

$$ds^2=dx^2+dy^2+dz^2c^2dt^2,$$
$`(1)`$
and for $`x^\mu \left(\sigma \right)`$ we have

$$ds^2=k\left(d\sigma \right)^2,$$
$`(2)`$
where $`d\sigma `$ is indeed an infinitesimal real quantity, and where the discrete dimensionless parameter $`k`$ is as follows:

$`k=1`$ in the bradyonic (timelike, subluminal) case,

$`k=0`$ in the luxonic (null, luminal) case, and

$`k=1`$ in the tachyonic (spacelike, superluminal) case.

(If the signature $`+`$ were chosen, then by Eq. (2) we would have $`k=1`$ in the bradyonic case and $`k=1`$ in the tachyonic case.) Dividing Eqs. (1) and (2) by $`\left(d\sigma \right)^2`$ we get

$$k=\left(u^x\right)^2+\left(u^y\right)^2+\left(u^z\right)^2\left(u^t\right)^2,$$
$`(3)`$
where $`u^\mu :=dx^\mu /d\sigma `$ is a four-velocity vector. The kinematical Eq. (3) concerns every type of world lines – timelike, null, and spacelike. The type is determined by $`k`$.

Multiplying Eq. (3) by $`m^2c^2`$, where $`m`$ has the mass dimension (we do not yet determine physical meanings of $`m`$), we get the well-known special relativistic formula for a four-momentum vector $`p^\mu `$:

$$km^2c^2=\left(p^x\right)^2+\left(p^y\right)^2+\left(p^z\right)^2\left(p^t\right)^2𝐩^2c^2E^2,$$
$`(4)`$
where

$$p^\mu :=mcu^\mu ,$$
$`(5)`$
and where $`\left(p^x\right)^2+\left(p^y\right)^2+\left(p^z\right)^2𝐩^2`$ and $`\left(p^t\right)^2c^2E^2`$. If we had $`m=0`$, then by definition (5) we would have no four-momentum, i.e. no object on our world line (not speaking that the multiplication of equations by zero does not make sense). Thus

$$m0.$$
$`(6)`$
If $`m`$ were imaginary, then by definition (5) also the four-momentum components $`p^\mu `$ would be imaginary, which would give us a new physics yet unknown. If we had real $`m<0`$, then by definition (5) we would have opposite senses of the four-vectors $`u^\mu `$ and $`p^\mu `$. Such a situation is yet unknown and today seems strange, though perhaps it will be considered in future. Anyway, we are entitled to put real $`m>0`$ for every type of the objects under consideration (Ockham’s principle!).

The unfortunate phrases have resulted from the fact that some authors have not taken into account the existence of three values of $`k`$ (1, 0, $`1`$) and have applied the bradyonic variants of Eqs. (1)–(4) for luxons and tachyons. The use of proper values of $`k`$ allows to avoid the difficulties.

In the bradyonic case, $`m`$ is the rest mass of our object. In the luxonic case the physical meaning of $`m`$ is not determined in general, though it is so for the photon for which $`m=c^2E=c^2h\nu >0`$. Anyway, the dynamical luxonic relation $`𝐩^2c^2=E^2`$ does not result from the condition $`m=0`$, which is false (inequality (6)), but it does result from the condition $`k=`$ 0, i.e. it is determined at the kinematical level of Eqs. (1)–(3). In the tachyonic case we have yet no operational definition of $`m`$ (for lack of rest), and therefore the term “masslike quantity” has been proposed. (The terms “pseudomass” or “quasimass” are shorter but semantically inferior.)

Additional remarks and amusing details can be found in Section 6 of my prior paper .

 J. K. Kowalczyński, Acta Physica Slovaca 50 (2000) 381, or hep-ph/9911441.


text:
# Spectral representation of the Casimir Force Between a Sphere and a Substrate

## Appendix A Appendix

By using the Green’s function from Eq. (8), one can find the density of states of the system as a function of the spectral variable $`u`$, as $`\rho (u)=\frac{1}{\pi }\mathrm{Im}\left[\mathrm{Tr}𝔾(u)\right].`$ In the particular case, when the sphere is metallic and its dielectric function is described by the Drude model, $`u(\omega )=\omega (\omega +i\tau )/\omega _p^2`$. Then, the density of states as a function of the frequency takes the following explicit form

$$\rho (\omega )=\frac{2w_p}{\pi }\underset{s}{}\sqrt{n_s}\left(\frac{\omega /\tau }{(\omega ^2\omega _p^2n_s)^2+(\omega /\tau )^2}\right),$$
(11)
and the energy can be calculated as $`V=\frac{1}{2}_0^{\mathrm{}}\mathrm{}\omega \rho (\omega )𝑑\omega .`$


text:
# Mathisson’ spinning electron : noncommutative mechanics & exotic Galilean symmetry, 66 years ago

## 1 Introduction

Non-commutative (quantum) mechanics, where the position coordinates satisfy

$$\{X_1,X_2\}=\theta ,$$
(1.1)
has been at the center of recent interest . In the plane and in the non-relativistic context, such theories are closely related to the “exotic” Galilean symmetry associated with the two-fold central extension of the planar Galilei group . A model which provides a physical realization of this symmetry has been presented by Lukierski, Stichel and Zakrzewski  who considered the acceleration-dependent Lagrangian

$$L=\frac{m\dot{\stackrel{}{x}}^2}{2}+\frac{\kappa }{2}\dot{\stackrel{}{x}}\times \ddot{\stackrel{}{x}}.$$
(1.2)
My aim here is to point out that the model of Lukierski et al. can actually be derived from that published by Mathisson in ‘37 , and further discussed by Weyssenhoff and Raabe . Not surprisingly, their theory shows interesting analogies also with Schrödinger’s Zitternde Elektron .

This Note is dedicated to the memory of these outstanding physicists who, with an extreme courage, tried to continue their scientific activity under those terrible years of World War II.

## 2 The Mathisson electron

Two years before the outbreak of World War II, Mathisson  proposed to describe a classical electron with the relativistic equations

$$\begin{array}{cc}m\dot{u}^\alpha +\frac{1}{c^2}S^{\alpha \sigma }\ddot{u}_\sigma =f^\alpha \hfill & \\ \dot{S}^{\alpha \beta }\frac{1}{c^2}S^{\alpha \sigma }\dot{u}_\sigma u^\beta +\frac{1}{c^2}S^{\beta \sigma }\dot{u}_\sigma u^\alpha =0\hfill & \end{array}$$
(2.3)
where $`m`$ is the mass, $`u^\alpha `$ the four-velocity, $`f^\alpha `$ the force; the dot means differentiation w.r.t. proper time. The antisymmetric tensor $`S^{\alpha \beta }`$ represents the spin of the electron and is assumed to satisfy the orthogonality condition

$$S^{\alpha \beta }u_\beta =0.$$
(2.4)
In the rest frame, the spatial components of $`S^{\alpha \beta }`$ form therefore a three-vector $`\stackrel{}{S}`$.

In the non-relativistic limit, $`\stackrel{}{S}`$ becomes a constant of the motion. In the absence of external force, the motion is \[apart of free motion along the direction of $`\stackrel{}{S}`$\], in the plane perpendicular to $`\stackrel{}{S}`$ and satisfies the third-order equation

$$m\ddot{x}_i=\kappa ϵ_{ij}\stackrel{\dot{}\dot{}\dot{}}{x_j},$$
(2.5)
where the new constant $`\kappa `$ has been defined by the Jackiw-Nair Ansatz

$$s=\kappa c^2,$$
(2.6)
$`s=|\stackrel{}{S}|`$ being the length of the spin vector. Eqn. (2.5) is precisely the equation of motion put forward by of Lukierski et al. .

From now on we drop the coordinate parallel to $`\stackrel{}{S}`$ and focus our attention to motion in the plane.

## 3 Conserved Quantities

The equations of motion (2.5) are associated with the Lagrangian (1.2). Then Lukierski et al. derive the conserved quantities associated to the space-time symmetries applying the higher-order version of Noether’s theorem. Let us now reproduce their results using elementary methods.

$``$ An obvious first integral of (2.5) is the momentum,

$$P_i=m\dot{x}_i+\kappa ϵ_{ij}\ddot{x}_j.$$
(3.7)
Eq. (2.5) is in fact $`\dot{P}_i=0`$.

$``$ Multiplying (2.5) by the velocity, $`\dot{\stackrel{}{x}}`$, yields a total time derivative, where we recognize the conserved energy,

$$H=\frac{m\dot{\stackrel{}{x}}^2}{2}+\kappa \dot{\stackrel{}{x}}\times \ddot{\stackrel{}{x}}.$$
(3.8)
$``$ Similarly, taking the vector product of (2.5) with $`\stackrel{}{x}`$ yields the conserved angular momentum,

$$J=m\stackrel{}{x}\times \dot{\stackrel{}{x}}+\frac{\kappa }{2}\dot{\stackrel{}{x}}^2\kappa \stackrel{}{x}\ddot{\stackrel{}{x}}.$$
(3.9)
$``$ A Galilean boost $`\stackrel{}{x}\stackrel{}{x}+\stackrel{}{b}t`$ shifts the momentum as $`\stackrel{}{P}\stackrel{}{P}+m\stackrel{}{b}`$. A rest frame where the momentum vanishes can be found, providing us with the conserved boost vector

$$K_i=mx_it\left(m\dot{x}_i+\kappa ϵ_{ij}\ddot{x}_j\right)+\kappa ϵ_{ij}\dot{x}_j.$$
(3.10)
Somewhat surprisingly, one more conserved quantity can be found.

$``$ the vector product of (2.5) with the acceleration, $`\ddot{\stackrel{}{x}}`$, yields the square of the acceleration,

$$I=\frac{\kappa ^3}{2m^2}\left(\ddot{\stackrel{}{x}}\right)^2,$$
(3.11)
where a constant factor has been included for later convenience.

$``$ curiously, multiplying (2.5) by $`\stackrel{\dot{}\dot{}\dot{}}{\stackrel{}{x}}`$ yields once again the same quantity, namely $`(m/2)\left(\ddot{\stackrel{}{x}}\right)^2=(m/\kappa )^3I`$.

The construction of this new quantity reminds one to that of angular momentum and of energy. Its precise origin will be clarified below.

Let us observe that, owing to the conservation of $`I`$, $`\ddot{\stackrel{}{x}}=0`$ can be consistently required. Then the conserved quantities found above reduce to those of an “elementary exotic particle” studied in .

## 4 Zitterbewegung and center-of-mass decomposition

The equation of motion (2.5) is integrated at once. Putting indeed

$$Q_i=\left(\frac{\kappa }{m}\right)^2ϵ_{ij}\ddot{x}_j,$$
(4.12)
eqn. (2.5) becomes

$$\dot{Q}_i=\frac{m}{\kappa }ϵ_{ij}Q_j,$$
(4.13)
showing that the acceleration rotates uniformly with angular velocity $`m/\kappa `$. Putting $`Q=Q_1+iQ_2`$, $`Q(t)=Q_0e^{i(m/\kappa )t}`$. This is plainly consistent with the conservation of the magnitude of the acceleration, Eq. (3.11). Then

$$X_i=x_i+ϵ_{ij}Q_j$$
(4.14)
moves freely,

$$\ddot{X}_i=0.$$
(4.15)
In conclusion, the motion has been separated into the free motion of the center of mass coordinate $`\stackrel{}{X}`$, combined with the “Zitterbewegung” \[uniform rotation\] of the internal coordinate $`\stackrel{}{Q}`$.

A key feature of Mathisson’s electron is that the internal variable $`\stackrel{}{Q}`$ measures in fact the extent of how much the momentum, $`\stackrel{}{P}`$, differs from \[$`m`$-times\] the velocity, $`\dot{\stackrel{}{x}}`$,

$$\stackrel{}{Q}=\frac{\kappa }{m^2}(m\dot{\stackrel{}{x}}\stackrel{}{P}).$$
(4.16)
Re-writing the conserved quantities in terms of the new coordinates confirms the above interpretation. In fact,

$$\begin{array}{ccc}\stackrel{}{P}& & =m\dot{\stackrel{}{X}}\\ H& =H_{CM}+H_{int}& =\frac{m\dot{\stackrel{}{X}}^2}{2}\frac{m^3}{2\kappa ^2}\stackrel{}{Q}^2\\ J& =J_{CM}+J_{int}& =m\stackrel{}{X}\times \dot{\stackrel{}{X}}+\frac{\kappa }{2}\dot{\stackrel{}{X}}^2+\frac{m^2}{2\kappa }\stackrel{}{Q}^2\\ K_i& & =m(X_i\dot{X}_it)+\kappa ϵ_{ij}X_j\\ I& & =\frac{m^2}{2\kappa }\stackrel{}{Q}^2.\end{array}$$
(4.17)
Mathisson’s electron is hence a composite system. Note that in (4.17) the center of mass behaves precisely as an elementary exotic particle ; the internal coordinate only contributes to the energy and the angular momentum. In fact, $`H_{int}={\displaystyle \frac{m}{\kappa }}I`$ and $`J_{int}=I`$. The new conserved quantity found in (3.11) is hence the internal angular momentum and also the internal energy \[which are linked in a $`2`$-dimensional phase space\].

Let us now observe that the equations of motion (4.13-4.15) are consistent with the Poisson structure associated with the symplectic form

$$\mathrm{\Omega }=\mathrm{\Omega }_{CM}+\mathrm{\Omega }_{int}=dP_idX_i+\frac{\kappa }{2m^2}ϵ_{ij}dP_idP_j+\frac{m^2}{\kappa }ϵ_{ij}dQ_idQ_j.$$
(4.18)
The $`6`$ dimensional phase space is hence the direct sum of the four-dimensional “exotic” phase space of the center of mass with coordinates $`\stackrel{}{X}`$ and $`\stackrel{}{P}`$, with the two-dimensional internal phase space of the $`\stackrel{}{Q}`$, endowed with a canonical symplectic structure.

The Poisson structure can be used to calculate the algebraic structure of the symmetries. Consistently with Lukierski et al. , we find that $`\stackrel{}{P},H,J,\stackrel{}{K}`$, supplemented with the central charges $`m`$ and $`\kappa `$, realize the “exotic” \[two-fold centrally extended\] planar Galilei group. The structure relations of this latter only differ from those of the usual Galilei group in that the Poisson bracket of the boost components yields the “exotic” central charge,

$$\{K_1,K_2\}=\kappa .$$
(4.19)
Similarly, the center-of-mass coordinates have a nonvanishing Poisson bracket,

$$\{X_1,X_2\}=\frac{\kappa }{m^2},\{Q_1,Q_2\}=\frac{\kappa }{m^2}.$$
(4.20)
Both the center-of-mass and the internal coordinates are hence noncommuting, cf. (1.1) with $`\theta =(\kappa /m^2)`$ \[while the original coordinates $`x_i`$ commute\]. This is similar to what happens in the Landau problem where the guiding center coordinates are noncommuting, with $`\theta =1/eB`$.

The additional conserved quantity $`I`$ in (3.11) is actually associated with the internal symmetries of the system. The translations and boosts form indeed an invariant subgroup $`K`$ of the Galilei group. The quotient $`G/K`$, which consists of rotations and time translations, is hence a group that can be made to act separately on the center-of-mass and the internal space. We can, e. g., rotate the internal coordinate $`\stackrel{}{Q}`$ alone and leave the center-of-mass coordinate $`\stackrel{}{X}`$ fixed. This is plainly a symmetry, and the associated conserved quantity is the internal angular momentum $`J_{int}=I`$. (A physical rotation moves both the external and internal coordinates, yielding the total angular momentum in (4.17)). The internal energy arises in a similar way. In conclusion, the non-relativistic limit of the Mathisson electron admits the direct product of the “exotic” Galilei group with the internal rotations and time translations, $`SO(2)\times 𝐑`$, as symmetry. Here the action of the Galilei group is transitive on the submanifolds $`I=const`$ i. e., $`\stackrel{}{Q}^2=const`$.

The same statement is valid for any composite nonrelativistic system, i. e. one upon which the Galilei group acts by symmetries but not transitively .

## 5 Relation to Schrödinger’s Zitternde Elektron

The results of Section 4 remind those Schrödinger derived in his original paper on Zitterbewegung . Schrödinger starts in fact with the Dirac Hamiltonian

$$H=c\stackrel{}{\alpha }\stackrel{}{P}+m^2c^2\beta $$
(5.21)
where $`\stackrel{}{\alpha }`$ and $`\beta `$ denote the usual Dirac matrices. In the Heisenberg picture, the operators satisfy

$$\frac{d\stackrel{}{P}}{dt}=0,\frac{dH}{dt}=0\frac{d\stackrel{}{x}}{dt}=c\stackrel{}{\alpha }.$$
(5.22)
The last equation can be rewritten as $`i\frac{d\stackrel{}{\eta }}{dt}=2H\stackrel{}{\eta },`$ ($`\mathrm{}=1`$), where $`\stackrel{}{\eta }=\stackrel{}{\alpha }cH^1\stackrel{}{P}.`$ This can be integrated as $`\stackrel{}{\eta }(t)=e^{2iHt}\stackrel{}{\eta }_0=\stackrel{}{\eta }_0e^{2iHt}`$, where $`\stackrel{}{\eta }_0`$ is a constant operator. Hence

$`{\displaystyle \frac{d\stackrel{}{x}}{dt}}=c^2H^1\stackrel{}{P}+c\stackrel{}{\eta }_0e^{2iHt},`$
which can again be integrated to yield

$$\stackrel{}{x}(t)=\left\{\stackrel{}{X}_0+c^2H^1\stackrel{}{P}t\right\}+{\scriptscriptstyle \frac{1}{2}}ic\stackrel{}{\eta }_0H^1e^{2iHt}$$
(5.23)
where $`\stackrel{}{X}_0`$ is a constant operator. The structure is clearly the same as in (4.14), with the operator

$$\stackrel{}{X}(t)=\stackrel{}{X}_0+c^2H^1\stackrel{}{P}t$$
(5.24)
representing the freely moving center-of-mass, and the second term describing the internal Zitterbewegung. The precise relation is more subtle, though. Intuitively, dropping the third component and working in the plane, putting $`s=1/2`$ and $`s/c^2\kappa `$ \[which would require the spin to diverge as $`c\mathrm{}`$ rather then remain a constant\], setting $`c\stackrel{}{\alpha }\dot{\stackrel{}{x}}`$ and replacing $`Hmc^2`$, would transform (5.23) formally into (4.14). In fact,

$`\stackrel{}{X}(t)\stackrel{}{X}_0+{\displaystyle \frac{\stackrel{}{P}}{m}}t,mc\stackrel{}{\eta }{\displaystyle \frac{m^2}{\kappa }}\stackrel{}{Q},e^{i2Ht}=e^{i(H/s)t}e^{i(m/\kappa )t}.`$ (5.25)
Note that

$$mc\stackrel{}{\eta }=mc\stackrel{}{\alpha }mc^2H^1\stackrel{}{P}m\dot{\stackrel{}{x}}\stackrel{}{P}.$$
(5.26)
consistently with (4.16).

A distinctive feature of Schrödinger’s Zitternde Elektron is that the center-of-mass coordinates satisfy the nontrivial commutation relation

$$[X_i,X_j]=ic^2E^2ϵ_{ijk}S_k,$$
(5.27)
where $`E=c\sqrt{\stackrel{}{P}^2+m^2c^2}`$ and $`\stackrel{}{S}=(i/4)\stackrel{}{\alpha }\times \stackrel{}{\alpha }`$ is the spin operator. If we assume that the spin is polarized in the third direction, $`S_3=s`$, and we consider the non-relativistic limit $`Emc^2+\stackrel{}{P}^2/2m`$ together with the Ansatz (2.6), we find for the planar components

$$[X_1,X_2]i\frac{s}{c^2m^2}=i\frac{\kappa }{m^2}$$
(5.28)
cf. (1.1) with $`\theta =\kappa /m^2`$. Let us remark that our procedure here is in fact the quantum version of the subtle non-relativistic limit proposed by Jackiw and Nair .

## 6 The relativistic description of Weyssenhoff and Raabe

Mathisson’s classical electron was further elaborated by Weyssenhoff and Raabe in a paper published after the War . They posit the equations

$$\begin{array}{cc}\dot{p}^\alpha =0,p^\alpha =mu^\alpha +\frac{1}{c^2}S^{\alpha \beta }\dot{u}_\beta ,\hfill & \\ \dot{S}^{\alpha \beta }=p^\alpha u^\beta p^\beta u^\alpha ,\hfill & \\ S^{\alpha \sigma }u_\sigma =0,\hfill & \end{array}$$
(6.29)
where $`m=\frac{1}{c^2}u_\beta p^\beta .`$ Eliminating $`p^\alpha `$ yields the relativistic Mathisson equations (2.3) once again. Eqns. (6.29) imply that $`m`$ is constant of the motion, $`\dot{m}=0`$, identified as the rest-mass of the particle. $`S_{\alpha \beta }S^{\alpha \beta }=s^2`$ is also a constant of the motion. They also observe that, owing to $`\dot{p}^\alpha =0`$, the quantity $`M`$ defined by $`p_\alpha p^\alpha =M^2c^2`$ is another constant of the motion. It is worth noting that the position satisfies again a third-order equation analogous to (2.5), namely

$$m\ddot{x}_\alpha =\frac{1}{c^2}S_{\alpha \sigma }\stackrel{\dot{}\dot{}\dot{}}{x}_\sigma .$$
(6.30)
Then Weyssenhoff and Raabe proceed to integrate the free relativistic equations of motion. In a suitable inertial frame (called the proper system) the spatial components, $`P_i`$, of the vector $`p^\alpha `$ can be made to vanish, so that its time component is $`Mc`$. In this frame $`\stackrel{}{S}`$ is constant. The mass is $`m=M/\sqrt{1(\stackrel{}{v}/c)^2}`$ where $`v_i=u_i\sqrt{1(\stackrel{}{v}/c)^2}`$ denotes the three-velocity. Hence the time component of the four-velocity is also constant so that the four-acceleration is proportional to the three-acceleration, $`\stackrel{}{a}=d^2\stackrel{}{x}/dt^2`$. Transforming from proper time to $`t`$, $`\stackrel{}{P}=0`$ reduces finally to

$$M\stackrel{}{v}+\frac{1}{c^2}\stackrel{}{S}\times \stackrel{}{a}=0.$$
(6.31)
The particle moves hence along a circle in the plane perpendicular to $`\stackrel{}{S}`$, with uniform angular velocity

$$\frac{mc^2}{s}\left(1\frac{\stackrel{}{v}^2}{c^2}\right).$$
(6.32)
In a general Lorentz frame, the motion is a superposition of such a motion with a uniform translation.

Our clue is to observe that in the non-relativistic limit these formulæ reduce, with the Jackiw-Nair Ansatz $`s=\kappa c^2`$ cf. (2.6), to those we derived in Section 4.

It is worth mentionning that the equations of Weyssenhoff and Raabe have again and again re-emerged in the course of the years. Consider, for example, (6.29) infive dimensions and for $`s={\scriptscriptstyle \frac{1}{2}}`$. Multiplication of $`p^\beta `$ with $`S_{\alpha \beta }`$ allows us to express the five-vector $`\dot{u}_\alpha `$ as

$$\dot{u}_\alpha =\frac{4}{c^2}S_{\alpha \sigma }p^\sigma $$
(6.33)
which, together with the remaining relations in (6.29) and the constraint $`u_\alpha u^\alpha =1`$, are precisely the equations proposed by Barut and Zanghi  as a “Kaluza-Klein” description of a classical Dirac electron.

## 7 Conclusion

In this Note we have shown that the non-relativistic limit of Mathisson’s classical spinning electron yields the acceleration-dependent model of Lukierski et al. . This latter has non-commuting coordinates and realizes the “exotic” Galilean symmetry.

Our results confirm once again the relation between the relativistic spin and the non-relativistic “exotic” structure, advocated by Jackiw and Nair . Their rule (2.6) is, however, a rather strange one, since it requires the spin to diverge as $`c\mathrm{}`$ so that $`s/c^2`$ remains finite. For this reason, the use of a Dirac equation valid for the fixed velue $`s={\scriptscriptstyle \frac{1}{2}}`$ \[as in Section 5 above\] is clearly illegitimate, and should be replaced by some anyon equation, valid for any real spin $`s`$ .

Another intriguing feature of this procedure is the following. While the relativistic model is associated with an irreducible representation of the Poincaré group, its dequantized & non-relativistic limit, namely the model of Lukierski et al., only carries a reducible representation of the Galilei group : irreducibility is lost in the procedure.

A final remark concerns the spin constraint (2.4) which appears to lie at the very root of the Zitterbewegung. Trading it for

$$S^{\alpha \beta }p_\beta =0$$
(7.34)
would in fact eliminate the Zitterbewegung altogether and lead to models of the type discussed in .

Acknowledgement. I am indebted to Professor J. Lukierski for sending me copies of those old Acta Physica Polonica papers, and also to Professor A. Staruszkiewicz who provided me with some biographical data.