text:
# Search for Rare and Forbidden Charm Meson Decays 𝐷⁰→𝑉ℓ⁺ℓ⁻ and ℎℎℓℓ

## Abstract

We report results of a search for flavor-changing neutral current (FCNC), lepton-flavor, and lepton-number violating decays of the $`D^0`$ (and its antiparticle) into 3 and 4-bodies. Using data from Fermilab charm hadroproduction experiment E791, we examine modes with two leptons (muons or electrons) and either a $`\rho ^0`$, $`\overline{K}^0`$, or $`\varphi `$ vector meson or a non-resonant $`\pi \pi `$, $`K\pi `$, or $`KK`$ pair of pseudoscalar mesons. No evidence for any of these decays is found. Therefore, we present branching-fraction upper limits at 90$`\%`$ confidence level for the 27 decay modes examined (18 new).

The E791 Collaboration has previously reported limits on rare and forbidden dilepton decays of the charged charm $`D`$ meson . Such measurements probe the SU(2)$`\times `$U(1) Standard Model of electroweak interactions in search of new mediators and couplings . We extend the methodology to 27 dilepton decay modes of the neutral $`D`$ meson. The modes are resonant $`D^0V\mathrm{}^+\mathrm{}^{}`$ decays, where $`V`$ is a $`\rho ^0`$, $`\overline{K}^0`$, or $`\varphi `$, and non-resonant $`D^0hh\mathrm{}\mathrm{}`$ decays, where $`h`$ is a $`\pi `$ or $`K`$. The leptons are either muons or electrons. Charge-conjugate modes are implied. The modes are lepton flavor-violating (e.g., $`D^0\rho ^0\mu ^+e^{}`$), or lepton number-violating (e.g., $`D^0\pi ^{}\pi ^{}\mu ^+\mu ^+`$), or flavor-changing neutral current decays (e.g., $`D^0\overline{K}^0e^+e^{}`$). Box diagrams can mimic FCNC decays, but only at the $`10^{10}`$ to $`10^9`$ level . Long range effects (e.g., $`D^0\overline{K}^0\rho ^0,\rho ^0e^+e^{}`$) can occur at the $`10^6`$ level . Numerous experiments have studied rare decays of charge -1/3 strange quarks. Charge 2/3 charm quarks are interesting because they might couple differently .

The data come from measurements with the Fermilab E791 spectrometer. A total of $`2\times 10^{10}`$ events were taken with a loose transverse energy requirement. These events were produced by a 500 GeV/$`c`$ $`\pi ^{}`$ beam interacting in a fixed target consisting of five thin, well-separated foils. Track and vertex information came from “hits” in 23 silicon microstrip planes and 45 wire chamber planes. This information and the bending provided by two dipole magnets were used for momentum analysis of charged particles. Kaon identification was carried out by two multi-cell Čerenkov counters that provided $`\pi /K`$ separation in the momentum range $`660`$ GeV/$`c`$ . We required that the momentum-dependent light yield in the Čerenkov counters be consistent for kaon-candidate tracks, except for those in decays with $`\varphi K^+K^{}`$, where the narrow mass window for the $`\varphi `$ decay provided sufficient kaon identification (ID).

Electron ID was based on transverse shower shape plus matching wire chamber tracks to shower positions and energies in an electromagnetic calorimeter . The electron ID efficiency varied from 62$`\%`$ below 9 GeV/$`c`$ to 45$`\%`$ above 20 GeV/$`c`$. The probability to misidentify a pion as an electron was $`0.8\%`$, independent of pion momentum.

Muon ID was obtained from two planes of scintillation counters. The first plane (5.5 m $`\times `$ 3.0 m) of 15 counters measured the horizontal position while the second plane (3.0 m $`\times `$ 2.2 m) of 16 counters measured the vertical position. There were about 15 interaction lengths of shielding upstream of the counters to filter out hadrons. Data from $`D^+\overline{K}^0\mu ^+\nu __\mu `$ decays  were used to choose selection criteria for muon candidates. Timing information from the smaller set of muon scintillation counters was used to improve the horizontal position resolution. Counter efficiencies, measured using muons originating from the primary target, were found to be $`(99\pm 1)\%`$ for the smaller counters and $`(69\pm 3)\%`$ for the larger counters. The probability for misidentifying a pion as a muon decreased with momentum, from about 6$`\%`$ at 8 GeV/$`c`$ to $`(1.3\pm 0.1)\%`$ above 20 GeV/$`c`$.

Events with evidence of well-separated production (primary) and decay (secondary) vertices were selected to separate charm candidates from background. Secondary vertices were required to be separated from the primary vertex by greater than $`12\sigma __L`$, where $`\sigma __L`$ is the calculated resolution of the measured longitudinal separation. Also, the secondary vertex had to be separated from the closest material in the target foils by greater than $`5\sigma __L^{}`$, where $`\sigma __L^{}`$ is the uncertainty in this separation. The vector sum of the momenta from secondary vertex tracks was required to pass within $`40\mu `$m of the primary vertex in the plane perpendicular to the beam. Finally, the net momentum of the charm candidate transverse to the line connecting the production and decay vertices had to be less than 300 MeV/$`c`$. Decay track candidates were required to pass approximately 10 times closer to the secondary vertex than to the primary vertex. These selection criteria and kaon identification requirements were the same for both the search mode and for its normalization signal (discussed below). The mass ranges used for the resonant masses were: $`\left|m_{\pi ^+\pi ^{}}m_{\rho ^0}\right|<150`$ MeV/$`c^2`$, $`\left|m_{K^{}\pi ^+}m_{\overline{K}^0}\right|<55`$ MeV/$`c^2`$, and $`\left|m_{K^+K^{}}m_\varphi \right|<10`$ MeV/$`c^2`$.

To determine our selection cuts we used a “blind” analysis technique. Before the selection criteria were finalized, all events having masses within a window $`\mathrm{\Delta }M_S`$ around the mass of the $`D^0`$ were “masked” so that the presence or absence of any potential signal candidates would not bias our choice of selection criteria. All criteria were then chosen by studying events generated by a Monte Carlo (MC) simulation program and background events, outside the signal windows, from real data. The criteria were chosen to maximize the ratio $`N_{MC}/\sqrt{N_B}`$, where $`N_{MC}`$ and $`N_B`$ are the numbers of MC and background events, respectively, after all selection criteria were applied. The data within the signal windows were unmasked only after this optimization. We used asymmetric windows for the decay modes containing electrons to allow for the bremsstrahlung low-energy tail. The signal windows were: $`1.83<M(D^0)<1.90`$ GeV/$`c^2`$ for $`\mu \mu `$ and $`1.76<M(D^0)<1.90`$ GeV/$`c^2`$ for $`ee`$ and $`\mu e`$ modes.

We normalize the sensitivity of our search to topologically similar hadronic 3-body (resonant) or 4-body (non-resonant) decays. One exception to this is the case of $`D^0\rho ^0\mathrm{}^\pm \mathrm{}^{}`$ where we normalize to nonresonant $`D^0\pi ^+\pi ^{}\pi ^+\pi ^{}`$ because there is no published branching fraction for $`D^0\rho ^0\pi ^+\pi ^{}`$. Table I lists the normalization mode used for each signal mode and the fitted number of data events ($`N_{\mathrm{Norm}}`$).

The upper limit for each branching fraction $`B_X`$ is calculated using the following formulae:

$$B_X=\frac{N_X}{N_{\mathrm{Norm}}}\frac{\epsilon _{\mathrm{Norm}}}{\epsilon _X}\times B_{\mathrm{Norm}};\frac{\epsilon _{\mathrm{Norm}}}{\epsilon _X}=\frac{f_{\mathrm{Norm}}^{\mathrm{MC}}}{f_X^{\mathrm{MC}}}.$$
(1)
$`N_X`$ is the 90$`\%`$ confidence level (CL) upper limit on the number of decays for the rare or forbidden decay mode $`X`$ and $`B_{\mathrm{Norm}}`$ is the normalization mode branching fraction obtained from the Particle Data Group . $`\epsilon _{\mathrm{Norm}}`$ and $`\epsilon _X`$ are the detection efficiencies while $`f_{\mathrm{Norm}}^{\mathrm{MC}}`$ and $`f_X^{\mathrm{MC}}`$ are the fractions of Monte Carlo events that are reconstructed and pass the final selection criteria, for the normalization and decay modes respectively.   TABLE I.:  Normalization modes used.    Decay Mode Normalization Mode $`N_{\mathrm{Norm}}`$    $`D^0\rho ^0\mathrm{}^\pm \mathrm{}^{}`$ $`D^0\pi ^+\pi ^{}\pi ^+\pi ^{}`$ 2049$`\pm `$53  $`D^0\overline{K}^0\mathrm{}^\pm \mathrm{}^{}`$ $`D^0\overline{K}^0\pi ^+\pi ^{}`$ 5451$`\pm `$72  $`D^0\varphi \mathrm{}^\pm \mathrm{}^{}`$ $`D^0\varphi \pi ^+\pi ^{}`$ 113$`\pm `$19  $`D^0\pi \pi \mathrm{}\mathrm{}`$ $`D^0\pi ^+\pi ^{}\pi ^+\pi ^{}`$ 2049$`\pm `$53  $`D^0K\pi \mathrm{}\mathrm{}`$ $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$ 11550$`\pm `$113  $`D^0KK\mathrm{}\mathrm{}`$ $`D^0K^+K^{}\pi ^+\pi ^{}`$ 406$`\pm `$41

The MC simulations use Pythia/Jetset  as the physics generator and model the effects of resolution, detector geometry, magnetic fields, multiple scattering, interactions in the detector material, detector efficiencies, and the analysis selection criteria. The efficiencies for the normalization modes varied from approximately $`0.2\%`$ to $`1\%`$ depending on the mode, and the efficiencies for the search modes varied from approximately $`0.05\%`$ to $`0.34\%`$. We take muon and electron ID efficiencies from data.

Monte Carlo studies show that the experiment’s acceptances are nearly uniform across the Dalitz plots, except that the dilepton identification efficiencies typically drop to near zero at the dilepton mass threshold. While the loss in efficiency varies channel by channel, the efficiency typically reaches its full value at masses only a few hundred MeV/$`c^{\mathrm{\hspace{0.17em}2}}`$ above the dilepton mass threshold. We use a constant weak-decay matrix element when calculating the overall detection efficiencies.

The 90$`\%`$ CL upper limits $`N_X`$ are calculated using the method of Feldman and Cousins  to account for background, and then corrected for systematic errors by the method of Cousins and Highland . In these methods, the numbers of signal events are determined by simple counting, not by a fit. All results are shown in Fig. 1 and listed in Table II. Upper limits are determined using the number of candidate events observed and expected number of background events within the signal region.

Background sources that are not removed by the selection criteria discussed earlier include decays in which hadrons (from real, fully-hadronic decay vertices) are misidentified as leptons. These misidentified leptons can come from hadronic showers reaching muon counters, decays-in-flight, and random overlaps of tracks from otherwise separate decays (“accidental” sources). In the case where kaons are misidentified as pions or leptons, candidate masses shift below signal windows. However, we remove these events to prevent them from influencing our background estimate, which is partially obtained from the mass sidebands (see discussion of $`N_{\mathrm{Cmb}}`$ below). To remove these events prior to the selection-criteria optimization, we reconstruct all candidates as each of the non-resonant normalization modes and test whether the masses are consistent with $`m_{D^0}`$. If so, we remove the events, but only if the number of kaons in the final state differs from that of the search mode. We do not remove events having the same number of kaons, as the loss in acceptance for true signal events would be excessive.

There remain two sources of background: hadronic decays where pions are misidentified as leptons ($`N_{\mathrm{MisID}}`$) and “combinatoric” background ($`N_{\mathrm{Cmb}}`$) arising primarily from false vertices and partially reconstructed charm decays. The background $`N_{\mathrm{MisID}}`$ arises from the normalization modes. To estimate the rate for misidentifying $`\pi \pi `$ as $`\mathrm{}\mathrm{}`$, for all but the $`D^0K^{}\pi ^+\mathrm{}^+\mathrm{}^{}`$ modes, we assume all $`D^0K^{}\pi ^+\mathrm{}^+\mathrm{}^{}`$ candidates observed (after subtracting combinatoric background estimated from mass sidebands) result from misidentification of $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$ decays and count the number of $`D^0K^{}\pi ^+\mathrm{}^+\mathrm{}^{}`$ decays passing the final selection criteria. We then divide by twice the number of $`D^0K^{}\pi ^+\pi ^{}\pi ^+`$ normalization events with $`K^{}\pi ^+\mathrm{}^+\mathrm{}^{}`$ mass within $`\mathrm{\Delta }M_S`$ boundaries (twice because there are two possible $`\pi ^+`$ misidentifications).

From this procedure, the following misidentification rates were obtained: $`r_{\mu \mu }=(3.4\pm 2.4)\times 10^4`$, $`r_{\mu e}=(4.2\pm 1.4)\times 10^4`$, and $`r_{ee}=(9.0\pm 6.2)\times 10^5`$. For modes in which two possible pion combinations can contribute, e.g., $`D^0K^{}\pi ^+\mu ^\pm \mu ^{}`$, we use twice the above rate; and for $`D^0\pi ^+\pi ^{}\pi ^+\pi ^{}`$, where there are 4 possible combinations, we use 4 times this rate in calculating $`D^0\pi ^+\pi ^{}\mathrm{}^+\mathrm{}^{}`$. Using these rates, we estimate the numbers of misidentified candidates, $`N_{\mathrm{MisID}}^V\mathrm{}\mathrm{}`$ and $`N_{\mathrm{MisID}}^{hh\mathrm{}\mathrm{}}`$, in the signal windows as follows:

$$N_{\mathrm{MisID}}^{hh\mathrm{}\mathrm{}}=r_{_{\mathrm{}\mathrm{}}}\times N_{\mathrm{Norm}}^{hh\pi \pi }\mathrm{and}N_{\mathrm{MisID}}^V\mathrm{}\mathrm{}=r_{_{\mathrm{}\mathrm{}}}\times N_{\mathrm{Norm}}^{V\pi \pi },$$
(2)
where $`N_{\mathrm{Norm}}^{hh\pi \pi }`$ and $`N_{\mathrm{Norm}}^{V\pi \pi }`$ are the numbers of normalization hadronic decay candidates in the signal windows.

To calculate the upper limits for the $`D^0K^{}\pi ^+\mathrm{}^+\mathrm{}^{}`$ modes, we set $`N_{\mathrm{MisID}}`$ to zero as we do not have an independent estimate of the misidentification rates. This results in conservative upper limits. If we had used the misidentification rates from our previous, 3-body decay study , then our limits for the three $`D^0K^{}\pi ^+\mathrm{}^+\mathrm{}^{}`$ modes would be lower by about a factor of two.

To estimate the combinatoric background $`N_{\mathrm{Cmb}}`$ within a signal window $`\mathrm{\Delta }M_S`$, we count events having masses within an adjacent background mass window $`\mathrm{\Delta }M_B`$, and scale this number ($`N_{\mathrm{\Delta }M_B}`$) by the relative sizes of these windows: $`N_{\mathrm{Cmb}}=(\mathrm{\Delta }M_S/\mathrm{\Delta }M_B)\times N_{\mathrm{\Delta }M_B}`$. To be conservative in calculating our 90$`\%`$ confidence level upper limits, we take combinatoric backgrounds to be zero when no events are located above the mass windows. Table II shows the numbers of combinatoric background, misidentification background, and observed events for all 27 modes.   TABLE II.:  E791 90$`\%`$ confidence level (CL) upper limits on the number of events and branching fraction limits ($`\times 10^5`$). Previously published limits  for the nine $`D^0V\mathrm{}^+\mathrm{}^{}`$ modes are 23, 10, 4.9, 118, 14, 10, 41, 5.2, and 3.4 $`\times \mathrm{\hspace{0.17em}10}^5`$.     Mode (Est. BG)  Sys.  E791  $`D^0`$ $`N_{\mathrm{Cmb}}`$ $`N_{\mathrm{MisID}}`$ $`N_{\mathrm{Obs}}`$ Err. $`N_X`$ Limit    $`\pi ^+\pi ^{}\mu ^+\mu ^{}`$ 0.00 3.16 2 11$`\%`$ 2.96 $`3.0`$  $`\pi ^+\pi ^{}e^+e^{}`$ 0.00 0.73 9 12$`\%`$ 15.2 $`37.3`$  $`\pi ^+\pi ^{}\mu ^\pm e^{}`$ 5.25 3.46 1 15$`\%`$ 1.06 $`1.5`$  $`K^{}\pi ^+\mu ^+\mu ^{}`$ 3.65 0.00 12 11$`\%`$ 15.4 $`35.9`$  $`K^{}\pi ^+e^+e^{}`$ 3.50 0.00 6 15$`\%`$ 7.53 $`38.5`$  $`K^{}\pi ^+\mu ^\pm e^{}`$ 5.25 0.00 15 12$`\%`$ 17.3 $`55.3`$  $`K^+K^{}\mu ^+\mu ^{}`$ 2.13 0.17 0 17$`\%`$ 1.22 $`3.3`$  $`K^+K^{}e^+e^{}`$ 6.13 0.04 9 18$`\%`$ 9.61 $`31.5`$  $`K^+K^{}\mu ^\pm e^{}`$ 3.50 0.17 5 17$`\%`$ 6.61 $`18`$  $`\rho ^0\mu ^+\mu ^{}`$ 0.00 0.75 0 10$`\%`$ 1.80 $`2.2`$  $`\rho ^0e^+e^{}`$ 0.00 0.18 1 12$`\%`$ 4.28 $`12.4`$  $`\rho ^0\mu ^\pm e^{}`$ 0.00 0.82 1 11$`\%`$ 3.60 $`6.6`$  $`\overline{K}^0\mu ^+\mu ^{}`$ 0.30 1.87 3 24$`\%`$ 5.40 $`2.4`$  $`\overline{K}^0e^+e^{}`$ 0.88 0.49 2 25$`\%`$ 4.68 $`4.7`$  $`\overline{K}^0\mu ^\pm e^{}`$ 1.75 2.30 9 24$`\%`$ 12.8 $`8.3`$  $`\varphi \mu ^+\mu ^{}`$ 0.30 0.04 0 33$`\%`$ 2.33 $`3.1`$  $`\varphi e^+e^{}`$ 0.00 0.01 0 33$`\%`$ 2.75 $`5.9`$  $`\varphi \mu ^\pm e^{}`$ 0.00 0.05 0 33$`\%`$ 2.71 $`4.7`$  $`\pi ^{}\pi ^{}\mu ^+\mu ^+`$ 0.91 0.79 1 9$`\%`$ 2.78 $`2.9`$  $`\pi ^{}\pi ^{}e^+e^+`$ 0.00 0.18 1 11$`\%`$ 4.26 $`11.2`$  $`\pi ^{}\pi ^{}\mu ^+e^+`$ 2.63 0.86 4 10$`\%`$ 5.18 $`7.9`$  $`K^{}\pi ^{}\mu ^+\mu ^+`$ 2.74 3.96 14 9$`\%`$ 15.7 $`39.0`$  $`K^{}\pi ^{}e^+e^+`$ 0.88 1.04 2 16$`\%`$ 4.14 $`20.6`$  $`K^{}\pi ^{}\mu ^+e^+`$ 0.00 4.88 7 11$`\%`$ 7.81 $`21.8`$  $`K^{}K^{}\mu ^+\mu ^+`$ 1.22 0.00 1 17$`\%`$ 3.27 $`9.4`$  $`K^{}K^{}e^+e^+`$ 0.88 0.00 2 17$`\%`$ 5.28 $`15.2`$  $`K^{}K^{}\mu ^+e^+`$ 0.00 0.00 0 17$`\%`$ 2.52 $`5.7`$

The sources of systematic errors in this analysis include: errors from the fit to the normalization sample $`N_{\mathrm{Norm}}`$; statistical uncertainty on the selection efficiencies, calculated for Monte Carlo simulated events, for both $`f_{\mathrm{Norm}}^{\mathrm{MC}}`$ and $`f_X^{\mathrm{MC}}`$; uncertainties in the calculation of misidentification background; and uncertainties in the relative efficiency for each mode, including lepton tagging efficiencies. These tagging efficiency uncertainties include: 1) muon counter efficiencies from hardware performance; and 2) the fraction of signal events (based on simulations) that would remain outside the signal window due to bremsstrahlung tails. Also, for the $`D^0\rho ^0\mathrm{}^+\mathrm{}^+`$ modes, an additional systematic error is included because we are using $`D^0\pi ^+\pi ^{}\pi ^+\pi ^{}`$ as the normalization mode since there is no published branching fraction for $`D^0\rho ^0\pi ^+\pi ^{}`$. The sums, taken in quadrature, of these systematic errors are listed in Table II.

In summary, we use a “blind” analysis of data from Fermilab experiment E791 to obtain upper limits on the dilepton branching fractions for 27 flavor-changing neutral current, lepton-number violating, and lepton-family violating decays of $`D^0`$ mesons. No evidence for any of these 3 and 4-body decays is found. Therefore, we present upper limits on the branching fractions at the 90$`\%`$ confidence level. Four limits represent significant improvements over previously published results. Eighteen of these modes have no previously reported limits.

We thank the staffs of Fermilab and participating institutions. This research was supported by the Brazilian Conselho Nacional de Desenvolvimento Científico e Tecnológico, CONACyT (Mexico), the Israeli Academy of Sciences and Humanities, the U.S.-Israel Binational Science Foundation, and the U.S. National Science Foundation and Dept. of Energy. The Universities Research Assn. operates Fermilab for the U.S. Dept. of Energy.


format: text

text:
# What is the Highest Plausible Redshift of Luminous Quasars?

## 1. Introduction

Quasars have long been believed to be powered by the accretion of gas onto supermassive black holes (Salpeter 1964; Zel’dovich 1964; Lynden-Bell 1969). Dynamical studies indicate that remnant black holes (BHs) indeed reside in the quiescent nuclei of most nearby galaxies (e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000, and references therein), and imply that BH formation is a generic consequence of galaxy formation. Although the activity of bright quasars peaks at $`z2.5`$, quasars are known to exist at higher redshifts. Approximately 200 quasars have been found so far at $`z>4`$ (the majority of which are in the Palomar Sky Survey111See http://astro.caltech.edu/~george/z4.qsos for a comprehensive list of known $`z>4`$ quasars.), including the record–holding bright quasar SDSS 1044-0125 at $`z=5.80`$, recently discovered by the Sloan Digital Sky Survey (SDSS; Fan et al. 2000).

Supermassive BHs can grow out of low-mass “seed” BHs through accretion or mergers (Rees 1984; Barkana & Loeb 2000a). In popular hierarchical models for structure formation, the first collapsed gaseous objects have low masses. It is natural to identify these galaxies with baryonic masses just above the cosmological Jeans mass, $`10^4\mathrm{M}_{}`$, as the first sites where seed BHs may form (Larson 2000). Here we postulate that supermassive BHs are the merger products of individual BHs that grow out of these seeds through gas accretion. The natural $`e`$–folding timescale for the growth of a single seed can be written as

$$t_{\mathrm{acc}}\frac{M_{\mathrm{bh}}}{\dot{M}_{\mathrm{bh}}}=4\times 10^7\left(\frac{ϵ}{0.1}\right)\eta ^1\mathrm{yr},$$
(1)
where $`ϵL_{\mathrm{bol}}/\dot{M}_{\mathrm{bh}}c^2`$ is the radiative efficiency for a mass accretion rate $`\dot{M}_{\mathrm{bh}}`$, and $`\eta L_{\mathrm{bol}}/L_{\mathrm{Edd}}`$ is the bolometric accretion luminosity in Eddington units, $`L_E=4\pi GM_{\mathrm{bh}}c\mu _em_p/\sigma _T`$. Here $`\sigma _T`$ is the Thomson cross-section and $`\mu _e=1.15`$ is the mean atomic weight per electron. The growth of a $`10^9M_{}`$ BH out of a stellar mass seed requires about $`\mathrm{ln}(10^9M_{}/10M_{})=18.4`$ $`e`$–folding times or $`7\times 10^8(ϵ/0.1)\eta ^1`$ yrs. This time is shorter than the age of the universe at $`z=5.8`$ if $`(ϵ/0.1)\eta ^11`$. However, the time available for the growth of galaxies is shorter than the entire age of the universe at the corresponding redshift.

Turner (1991) recognized that exceptionally bright, high–redshift quasars may yield interesting constraints on cosmic structure formation. Motivated by the recent discovery of the SDSS 1044-0125 quasar at $`z=5.80`$ and the emergence of a concordance cosmological model for structure formation222Throughout this Letter we adopt a $`\mathrm{\Lambda }`$CDM cosmology with $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_\mathrm{b}=0.045`$, $`h=0.7`$, $`\sigma _{8h^1}=0.9`$, and $`n=1`$., we revisit this problem. In this Letter, we derive new constraints from the observed properties of known high-redshift quasars, and assess the potential for tightening these constraints with future quasar surveys.

## 2. Cosmological Growth of Quasar Black Holes

Our constraints derive from the existence of a BH of a particular mass $`M_{\mathrm{bh}}`$ at redshift $`z`$, within a galaxy halo of total mass $`M_{\mathrm{halo}}`$. In the following sub-sections we first express these masses in terms of observed quantities, and then describe the corresponding theoretical modeling of the BH growth.

### 2.1. Black Hole Mass

For the high-redshift quasars of interest, we can only infer the BH mass indirectly from the observed quasar luminosity (although reverberation mapping provides constraints in other cases; see Kaspi et al. 2000). Given the luminosity in the observed band, we first apply a bolometric correction to obtain the total luminosity $`L_{\mathrm{bol}}`$ (Elvis et al. 1994). For example, the approximate spectral slope of $`F_\nu \nu ^1`$ implies that only $`1\%`$ of the bolometric luminosity of SDSS 1044-0125 is in the SDSS $`z^{}`$ band. The BH mass then follows from the relation, $`L_{\mathrm{bol}}=\eta L_{\mathrm{Edd}}=1.4\times 10^{38}\eta (M_{\mathrm{bh}}/\mathrm{M}_{})\mathrm{erg}\mathrm{s}^1`$. In the case of SDSS 1044-0125, we infer $`M_{\mathrm{bh}}=3.4\times 10^9\eta ^1\mathrm{M}_{}`$. There are no obvious signs of beaming or lensing in this source. However, if present, both of these effects would reduce our inferred BH mass or increase our inferred constraint on $`\eta `$ by correcting the relation between the apparent and the intrinsic luminosity of the quasar. In §3 we will explore the sensitivity of our results to such a correction.

### 2.2. Black Hole Growth

In order to model the growth of BHs in our adopted $`\mathrm{\Lambda }`$CDM cosmology, we rely on the merger history of dark matter halos in the extended Press–Schechter (EPS) formalism (Press & Schechter 1974; Bond et al. 1991; Lacey & Cole 1993). We first need to determine the mass of the halo in which an observed quasar resides. We estimate $`M_{\mathrm{halo}}`$ based on the halo abundance within the EPS formalism. The number of halos with mass $`>M_{\mathrm{halo}}`$, in a survey of solid angle $`\mathrm{\Delta }\mathrm{\Omega }`$, probing a redshift range $`\mathrm{\Delta }z`$, is given by $`N(>M_{\mathrm{halo}})=\mathrm{\Delta }\mathrm{\Omega }\mathrm{\Delta }z\times (dV/d\mathrm{\Omega }dz)\times _{M_{\mathrm{halo}}}^{\mathrm{}}𝑑M[dn/dM(z,M)]`$, where $`dV/d\mathrm{\Omega }dz`$ is the cosmological volume element, and $`dn/dM`$ is the number of halos per halo mass per unit volume. The halo mass for a quasar of luminosity $`L`$ is found by equating the observed number of quasars brighter than $`L`$ to $`f(z)\times N`$, where $`f(z)`$ is the fraction of halos hosting active quasars. Because of the exponential shape of the EPS mass function for rare halos, our inferred halo mass depends only logarithmically on $`f(z)`$, and we conservatively adopt $`f(z)1`$. In the case of SDSS 1044-0125, one bright quasar was found within a $`600`$ deg2 survey area at redshift $`z=5.8\pm 0.5`$; from this we infer $`M_{\mathrm{halo}}1.1\times 10^{13}\mathrm{M}_{}`$, with a corresponding velocity dispersion of $`\sigma 470\mathrm{km}\mathrm{s}^1`$. It is interesting that the values of $`M_{\mathrm{bh}}`$ and $`\sigma `$ we obtain for SDSS 1044-0125 are in good agreement with their relation measured in local galaxies (Gebhardt et al. 2000; Ferrarese & Merritt 2000).

We next compute the growth of the central BH mass based on the assembly of its host halo. Given a halo of total mass $`M_{\mathrm{halo}}`$ at a redshift $`z`$, the EPS formalism specifies its average merger history at higher redshifts. Every branch of this merger tree represents a progenitor of the parent halo, whose mass is continuously growing through accretion or mergers with other halos. To keep our model simple, we assume that every building block of the original halo develops a seed BH of mass $`M_{\mathrm{seed}}`$ as soon as it acquires a minimum mass, $`M_{\mathrm{min}}`$, corresponding to a velocity dispersion $`\sigma _{\mathrm{min}}`$. The physical motivation for this choice is that radiative cooling and feedback processes which likely determine whether or not a massive BH forms in a collapsed halo, depend directly on $`\sigma `$ (e.g., Haehnelt et al. 1998). The mass of each seed BH is assumed to grow exponentially by accretion, $`M_{\mathrm{bh}}(t)=\mathrm{exp}[\mathrm{\Delta }t(z)/t_{\mathrm{acc}}]M_{\mathrm{seed}}`$, where $`\mathrm{\Delta }t`$ is the time elapsed between the formation time of the seed BH and the redshift $`z`$. We assume that eventually all massive BHs merge together to form a single supermassive BH at the center of the parent halo. As long as the BH mergers are completed prior to redshift $`z`$, there is no need to specify when these mergers took place. The mass of the resulting BH in the parent halo is the sum of the individual BHs, each of which has grown by a different amount,

$$M_{\mathrm{bh}}(z,M_{\mathrm{halo}})=M_{\mathrm{seed}}_{\mathrm{}}^z𝑑z^{}\frac{dN_{\mathrm{prog}}}{dz^{}}\mathrm{exp}\left[\frac{\mathrm{\Delta }t(z,z^{})}{t_{\mathrm{acc}}}\right],$$
(2)
where $`N_{\mathrm{prog}}(z^{})`$ is the total number of seeded progenitors at redshift $`z^{}>z`$,

$$N_{\mathrm{prog}}(z^{})=_{M_{\mathrm{min}}}^{M_{\mathrm{halo}}}𝑑M\frac{dP(z,z^{},M_{\mathrm{halo}},M)}{dM}.$$
(3)
Here $`dP(z,z^{},M_{\mathrm{halo}},M)`$ is the number of progenitors in the mass range between $`M`$ and $`M+dM`$ at a redshift $`z^{}`$ for a halo whose mass at redshift $`z`$ is $`M_{\mathrm{halo}}`$ \[see Lacey & Cole 1993, Eq. (2.15)\]. Note that we adopt the most optimistic assumptions regarding the BH growth; in reality the growth of each seed BH may be limited by its fuel reservoir. Also, the central BHs may not get incorporated into a single supermassive BH during galaxy mergers (Hut & Rees 1992); however, for exponential growth, the final BH mass tends to be dominated anyway by the first seed to have formed.

In summary, our model for the assembly of BHs has four free parameters: $`\sigma _{\mathrm{min}}`$, $`M_{\mathrm{seed}}`$, $`ϵ`$, and $`\eta `$. Although the values of these parameters are a’priori uncertain, their “natural” choices are as follows: (i) if the seed BHs are the remnants of massive stars, then their characteristic mass is $`M_{\mathrm{seed}}10\mathrm{M}_{}`$; (ii) a necessary requirement for star–formation is efficient cooling; in a metal–poor primordial gas with no $`\mathrm{H}_2`$ molecules, this implies a minimum virial temperature of $`T>\mathrm{\hspace{0.33em}10}^4`$K for progenitor halos with seed BHs, corresponding to a minimum velocity dispersion $`\sigma _{\mathrm{min}}>\mathrm{\hspace{0.33em}10}\mathrm{km}\mathrm{s}^1`$; (iii) if the gas accretes to a non–rotating BH through a steady thin disk, then the radiative efficiency, $`ϵ0.06`$; and (iv) for high fueling rates bright quasars would naturally shine close to their limiting luminosity, implying $`\eta 1`$. Variations on these parameter values will be discussed in § 3 and § 4.

## 3. Results

Given the parameters $`\sigma _{\mathrm{min}}`$, $`M_{\mathrm{seed}}`$, $`ϵ`$, and $`\eta `$, our model yields the mass $`M_{\mathrm{bh}}`$ of the supermassive BH at the center of a halo of mass $`M_{\mathrm{halo}}`$ at redshift $`z`$. In the following, we derive constraints on these four parameters, by requiring that $`M_{\mathrm{bh}}`$ equals the value estimated for observed high–redshift quasars.

### 3.1. Illustrative Example

In order to illustrate the BH growth process in our model, we show in Figure 3 the evolution of various quantities for the SDSS 1044-0125 quasar. In this example, we have assumed $`M_{\mathrm{seed}}=10\mathrm{M}_{},ϵ=0.1,\eta =1`$, and $`\sigma _{\mathrm{min}}=50\mathrm{km}\mathrm{s}^1`$. As we find by numerical integration of equations (2) and (3), this combination yields the required BH mass of $`M_{\mathrm{bh}}3.4\times 10^9\mathrm{M}_{}`$ at $`z=5.8`$. The top left panel in Figure 3 shows the number of progenitors of the parent halo with a velocity dispersion $`>\sigma _{\mathrm{min}}`$. For reference, the bottom left panel shows the corresponding minimum progenitor mass. As the redshift increases, the number of progenitors increases, peaking at $`z9`$, and then decreasing again as the typical progenitors are broken up into halos with $`\sigma <50\mathrm{km}\mathrm{s}^1`$. The top right panel shows the contribution of progenitors from each redshift to the final BH mass, and demonstrates that the bulk of the BH mass is contributed by seeds from $`z15`$. This redshift is considerably higher than the peak at which most progenitors form. The increased time available between higher redshifts and $`z=5.8`$ (shown explicitly in the bottom right panel) makes the contribution from the first few progenitors dominant.

### 3.2. Constraints from SDSS 1044-0125

Next, we search for general combinations of $`M_{\mathrm{seed}}`$, $`ϵ`$, $`\eta `$, and $`\sigma _{\mathrm{min}}`$ that yield the required BH mass for the SDSS quasar. We solve equation (2) with its left hand side set to $`3.4\times 10^9\eta ^1\mathrm{M}_{}`$, and the halo mass set to $`1.1\times 10^{13}\mathrm{M}_{}`$ at $`z=5.8`$, as discussed above. For a given pair of $`(M_{\mathrm{seed}},\sigma _{\mathrm{min}})`$, we then find $`\eta `$ as a function of $`ϵ`$ by a Newton–Rhapson method.

The results of this procedure are displayed in Figure 3.2 in a fiducial model with $`(M_{\mathrm{seed}},\sigma _{\mathrm{min}})=(10\mathrm{M}_{},10\mathrm{km}\mathrm{s}^1)`$ (solid line) and its variants. The dotted and short–dashed curves show results when $`\sigma _{\mathrm{min}}`$ is either increased to $`50\mathrm{km}\mathrm{s}^1`$ or decreased to $`1\mathrm{km}\mathrm{s}^1`$, respectively. The high value applies if the limiting mass for a progenitor halo is determined by feedback from the UV background following the reionization epoch (Thoul & Weinberg 1996; Navarro & Steinmetz 1997; Haiman, Madau & Loeb 1999); and the low value applies before reionization if sufficient $`\mathrm{H}_2`$ exists to allow cooling in small halos (Haiman, Abel & Rees 2000). The long–dashed curve demonstrates the effect of increasing the typical seed mass to $`10^3\mathrm{M}_{}`$, corresponding to remnants of very massive metal-free Population III stars (VMOs; Bond, Arnett & Carr 1984; see also Schneider et al. 2000). The dot–short–dashed curve describes the case where the observed emission from SDSS 1044-0125 is beamed into a tenth of its sky, so that its true BH mass is reduced by a factor of 10 and the halo abundance is increased by the same factor. Finally, the dot–long–dashed curve shows results for the fiducial model in the case of the exceptionally bright $`z=4.6`$ quasar PSS1347+4956. We have searched through the list of known $`z>4`$ quasars and found this object (BH mass of $`5\times 10^{10}\mathrm{M}_{}`$, and a host halo mass of $`2.7\times 10^{13}\mathrm{M}_{}`$ ($`\sigma =570\mathrm{km}\mathrm{s}^1`$, assuming no beaming or lensing) to provide the second strongest constraint after SDSS 1044-0125. A handful of other bright $`z>4`$ quasars follow closely behind.

Figure 3.2 leads to two interesting conclusions. First, the set of rather standard values, $`M_{\mathrm{seed}}=10\mathrm{M}_{}`$, $`\sigma _{\mathrm{min}}=50\mathrm{km}\mathrm{s}^1`$, $`ϵ=0.1`$, and $`\eta =1`$, yields the required BH mass of SDSS 1044-1215. The corresponding progenitor halos have cooling times much shorter than the dynamical time, and could withstand photo-ionization heating by the intergalactic UV background (e.g. Thoul & Weinberg 1996; Navarro & Steinmetz 1997). In such halos there is no obvious obstacle to the formation of massive stars, which could leave behind the required $`10\mathrm{M}_{}`$ seed BHs. Second, the inferred values of $`\eta `$ are relatively high, suggesting that BHs must radiate close to their Eddington limit. Lower values of $`\eta `$ are allowed if $`ϵ<\mathrm{\hspace{0.33em}0.05}`$; however, such low values of $`ϵ`$ would not account for the total energy output of quasars given the observed mass density of their remnants (e.g. Ho & Kormendy 2000; Salucci et al. 1999; Fabian 2000). The radiative efficiency may in fact approach a value as high as $`ϵ=0.42`$ for a maximally rotating Kerr BH.

### 3.3. Constraints from Future Quasar Surveys

In anticipation of future results from the SDSS survey as well as from deeper surveys of high–$`z`$ quasars, we ask: what is the highest plausible value of $`M_{\mathrm{bh}}`$ that should be anticipated for a quasar around a redshift $`z`$ in a given area on the sky? To answer this question, we compute the mass of the central BH in the most massive dark halo that could be found in a representative survey area of 10 deg2 and redshift bin $`\mathrm{\Delta }z=1`$ around different redshifts. For our set of fiducial parameter values we take $`M_{\mathrm{seed}}=10\mathrm{M}_{}`$, $`\sigma _{\mathrm{min}}=10\mathrm{km}\mathrm{s}^1`$, $`ϵ=0.06`$, and $`\eta =1`$. In Figure 3.2, the solid curve in the lower panel shows the maximum allowed quasar BH mass in this model as a function of redshift. In the upper panel, we show the corresponding apparent $`K(AB)`$ magnitude with the corresponding K–correction (Elvis et al. 1994). In our fiducial model, equation (2) leads to an unphysically large BH mass, in excess of the observed halo mass. In Figure 3.2 we do not allow the BH mass to exceed $`(\mathrm{\Omega }_b/\mathrm{\Omega }_0)M_{\mathrm{halo}}`$. This additional constraint appears as the break occurring at $`z9`$; BHs below this redshift are assumed conservatively to have consumed all the available gas in their host halo. Figure 3.2 also shows the maximum allowed BH mass and the corresponding minimum apparent magnitude, in variants of our fiducial model: examples in which $`\sigma _{\mathrm{min}}`$ is either increased to $`50\mathrm{km}\mathrm{s}^1`$ or decreased to $`1\mathrm{km}\mathrm{s}^1`$ (dotted and short–dashed curves), and a model in which the typical seed mass is increased to $`10^3\mathrm{M}_{}`$ (long–dashed curves). The dot–long–dashed curves show the results when we allow BHs to consume only 1% of the available gas in its host halo. Finally, the dot–short–dashed curves that bracket our fiducial model assume survey areas of 0.1 deg2 and 1000 deg2.

We find that the survey area enters only logarithmically into our constraints. However, the allowed BH mass and the resulting apparent magnitude, are strong functions of redshift. At its planned sensitivity, the Next Generation Space Telescope (NGST) will have a detection threshold of $``$32 mag in the $`15\mu `$m range (with $`<\mathrm{\hspace{0.33em}3}`$ hours of integration and S/N=10). The assumed survey size of 0.1 deg2 can be covered by 23 images taken by NGST in a total observation time of 2.3 days. The proposed PRIME survey333A Small Explorer Mission selected by NASA for study. See a description at http://spacescience.nasa.gov/codesr/smex. plans to map an area of $`10`$ deg2 to $`K=27`$. At a limiting magnitude of $`K20(30)`$, future surveys can yield strong constraints on the growth of supermassive black holes, as they can map out the “brightest quasar” envelope shown in Figure 3.2 out to $`z10(15)`$.

## 4. Conclusions

Existing data on high-redshift quasars implies that BHs as massive as $`3\times 10^9\mathrm{M}_{}`$ were assembled when the universe was only a tenth of its present age. Figure 2 shows that the massive BH inferred for SDSS 1044-0125 at $`z=5.8`$ can grow in hierarchical galaxy formation models with plausible parameter values for the initial seed mass ($`10\mathrm{M}_{}`$), the minimum velocity dispersion of collapsed objects that harbor such a seed ($`50\mathrm{km}\mathrm{s}^1`$), the radiative efficiency ($`6\%`$), and the luminosity in Eddington units ($`1`$). Figure 3 illustrates the upper envelope of plausible luminosity values for high-redshift quasars that might be found in future surveys.

There are several caveats to the constraints we derived. First, beaming or lensing may affect the apparent magnitude of some of the brightest and highest redshift quasars (Barkana & Loeb 2000b). Second, we have assumed that $`\eta `$ and $`ϵ`$ maintain the same values during the luminous quasar phase and the main growth phase of the BH mass. It is possible, however, that $`ϵ`$ is negligibly small during the early growth phase of the BH, and high during the luminous quasar phase (e.g. Haehnelt & Rees 1993). This would be equivalent to a corresponding increase in the seed BH mass in our model. Multiple seeds per halo which eventually coalesce are also equivalent to a single massive seed with the sum of their masses. Finally, although a recent cosmological simulation finds good agreement with the Press–Schechter mass function for the lowest halo masses of interest here at $`z=10`$ (Jang-Condell & Hernquist 2000); the accuracy of this ansatz still remains to be tested over a wider range of redshifts and halo masses, where there might be systematic deviations (Jenkins et al. 2000; Sheth & Tormen 2000).

Recently, various modifications to the standard $`\mathrm{\Lambda }`$CDM model have been proposed because of a potential conflict between theory and observations on small spatial scales (e.g., Kamionkowski & Liddle 2000; Bode, Ostriker & Turok 2000; Barkana, Haiman & Ostriker 2000, and references therein). Our constraints are expected to tighten significantly in such models, since they generically reduce small scale power and thus eliminate potential hosts of seed BHs at high redshifts. Future surveys of high-redshift quasars will test these models.

We thank Xiaohui Fan, Michael Strauss, and Pat Hall for useful discussions. ZH was supported by NASA through the Hubble Fellowship grant HF-01119.01-99A, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555. This work was supported in part by NASA grants NAG 5-7039, 5-7768, and NSF grants AST-9900877, AST-0071019 for AL.


text:
# Two-loop Perturbative Quark Mass Renormalization from Large 𝛽 Monte Carlo

## 1 INTRODUCTION

As actions for lattice calculations become more complicated, short distance calculations for normalizing the operators in the action become more and more difficult. It is useful to have as many tools as possible to attack the challenge. Several types of automated calculation of Feynman diagrams and nonperturbative methods have been developed. Another alternative is the use of direct calculation of short distance quark–gluon Green functions by Monte Carlo methods. Last year this method was applied to a determination of the additive mass renormalization for static quarks to third order. In this paper, we present a preliminary estimate of the mass renormalization for heavy Wilson quarks.

## 2 SIMULATION

There are several potential problems in naively performing simulations at large $`\beta `$ to extract perturbative information. One is the problem associated with infra-red divergences coming from the zero modes. The second is the phenomenon of tunneling between the various vacua of $`SU(3)`$. The third is finite volume effects, since the physical sizes of space-time used in the simulations are quite small. All of these effects may be removed or reduced by using twisted boundary conditions for the fields. These special boundary conditions remove the zero modes completely and suppresses tunneling for sufficiently large $`\beta `$. The remaining finite size effects can be removed by extrapolating to the infinite volume limit.

### 2.1 Twisted Boundary Conditions

We use a 3-way periodic twist (spatial directions only) for the gauge fields defined by,

$`U_\mu (x+L_j\widehat{e}_j)`$ $`=`$ $`\mathrm{\Omega }_jU_\mu (x)\mathrm{\Omega }_j^1`$ (1)
$`\mathrm{\Omega }_x\mathrm{\Omega }_y`$ $`=`$ $`\eta \mathrm{\Omega }_y\mathrm{\Omega }_x`$ (2)
where $`\eta =e^{2\pi i/3}`$ and $`\mathrm{\Omega }_z=\mathrm{\Omega }_y\mathrm{\Omega }_x^{}`$. The $`\mathrm{\Omega }`$’s are constant $`SU(3)`$ matrices.

### 2.2 Gauge Fixing

We must fix the gauge in order to extract energies from quark propagators. Our gauge fixing procedure is as follows:
(1) fix to Coulomb gauge, then
(2) fix the temporal links by maximizing,

$$tr\left[Z(t)\underset{\stackrel{}{x}}{}U_t(\stackrel{}{x},t)Z(t+1)\right]$$

where $`Z(t)`$ are elements of the center of $`SU(3)`$.

### 2.3 Parameters

The following values of the coupling were chosen: $`\beta =9,10,12,15,20,40`$ and $`60`$. Lattice sizes of $`L_s=\{4,6,8,10,12,14,16\}`$ and $`L_t=16`$ were used at each coupling to allow for an extrapolation to the infinite volume limit. In this preliminary study, we have accumulated $`5002000`$ configurations separated by 100 sweeps. The errors after binning did not show significant residual correlations.

## 3 STATIC QUARK SELF-ENERGY

We test our gauge-fixed configurations by calculating the first two coefficients of the static quark self-energy and compare to analytic and previous high $`\beta `$ Monte-Carlo calculations where no gauge fixing was performed. We also compare results from the Polyakov loop and the static quark propagator.

### 3.1 Polyakov Loop

The tadpole-improved static quark self-energy (in a finite volume) was calculated from the measured value of the Polyakov loop according to,

$$E_0(L)=\frac{1}{L}\mathrm{l}n\left[\frac{1}{3}\mathrm{𝑅𝑒𝑇𝑟}\left\{\underset{t=1}{\overset{L}{}}\frac{U_t(t)}{u_0}\right\}\right]$$
(3)
where $`u_0`$ is the mean link from the average plaquette. Our fitting procedure to extract the renormalization constants is as follows:

* Fit the measured values at fixed $`L`$ to a polynomial in $`\alpha _V(q^{})`$ to obtain $`c_1(L)`$ and $`c_2(L)`$, i.e. $`E_0(L)/\alpha =c_1(L)+\alpha c_2(L)+\alpha ^2c_3(L)`$.
* Extrapolate the coefficient $`c_1(L)`$ using $`c_1(L)=c_1+X_1/L`$ to get the one-loop coefficient, $`c_1`$.
* Extrapolate $`c_2(L)`$ with $`c_2(L)=c_2+X_2/L+\beta _0X_1ln(L)/L`$ where $`\beta _0=11/2\pi `$ to get the two loop coefficient, $`c_2`$.

The functional form and the coefficient of the logarithm comes from the anticipation of the (running) Coulomb interaction with the image charges with some effective distance $`L/X_1`$.

We also perform the fit in the reverse order, i.e. extrapolate to infinite volume at each $`\beta `$ and then fit to a cubic polynomial in $`\alpha _V`$. In the infinite volume extrapolation, we use $`E_0(L)=E_0+Y_1/L+\beta _0X_1ln(L)/L`$ where $`X_1`$ is the coefficient from the previous fit. The results were in perfect agreement. The renormalized coupling in the V-scheme, $`\alpha _V`$, was evaluated at $`q^{}=0.84a^1`$.

### 3.2 Static Quark Propagator

On the same set of lattices, we have measured the static quark propagator averaged over the spatial sites. These were then fit to single exponentials starting from timeslice 2 for which we get a good fit. The perturbative coefficients were extracted in the exact same manner as in the previous section.

### 3.3 Comparison

We compare our results with the analytical calculations of Heller and Karsch and with previous measurements of the first two coefficients in Table 1. Our measurements on gauge fixed configurations with both gauge invariant and covariant operators are in good agreement with those in the literature. These results give us confidence in our gauge fixing procedure.

## 4 WILSON QUARK REST MASS

We now apply similar methods to the Wilson quark mass renormalization. In order to accommodate twisted boundary conditions consistently for propagating quarks, we enlarge the quark internal degrees of freedom to $`N_{spin}\times N_{color}\times N_{smell}`$ where the smell group is $`SU(N_{smell}=3)`$. The boundary condition for the quark field reads,

$$\mathrm{\Psi }(x+L_j\widehat{e}_j)=\mathrm{\Omega }_j\mathrm{\Psi }(x)\mathrm{\Omega }_j^1e^{i\pi /N_{color}}$$
(4)
where the phase factor ensures anti-periodicity. The momenta of the quarks are then given by,

$$p_i=\frac{\pi }{N_{color}L_i}(2n_i+1)$$
(5)
$$n_xn_y+n_z=mN_{color}$$
(6)
where $`m`$ and $`n_i`$ are integers. The form of the constraint on $`n_i`$ is due to our particular choice of the twist matrix in the z-direction.

We have measured the ground state only, $`(n_x=0,n_y=0,n_z=0)`$, for which the color-smell wavefunction is proportional to the identity matrix. The extrapolation to the infinite volume limit will then give the rest mass, $`E(\stackrel{}{p}=\stackrel{}{0})`$,

$$M_1=M_1^{[0]}+\alpha M_1^{[1]}+\alpha ^2M_1^{[2]}+\alpha ^3M_1^{[3]}+\mathrm{}$$
(7)
where $`M_1^{[0]}=ln(1+am_0)`$, $`am_0=1/(2\kappa )4`$ and $`M_1^{[1]}`$ is calculated in Ref. . As explained in Ref. , the parameters in the Wilson and clover actions should be taken as mass-dependent when the masses are large. Here we calculate the mass renormalization for $`\kappa =0.1`$.

The ground state momenta was projected out both at the source and the sink. Only the spin components 0 and 2 (set equal to each other) were used in the simulation.

### 4.1 Fitting Procedure

We have fitted the quark propagator to a $`cosh`$ starting at time-slice 3 for which we get a good $`\chi ^2`$. The tree-level mass, $`M_1^{[0]}(L)`$, was then subtracted from the raw data so as to cancel some of the finite volume dependence of the rest mass. The extraction of the one-loop coefficient, $`M_1^{[1]}`$, and the two-loop coefficient, $`M_1^{[2]}`$, was done using the same procedure as the static quarks. However, we have found that a $`1/L^2`$ and $`1/L^3`$ term was required to get a good fit to the one-loop coefficient (see Fig. 2). The second order coefficient was extracted with a $`1/L`$ and $`ln(L)/L`$ term.

As in the case of the static energy, we fit the masses in the reverse order to check for consistency. The infinite volume extrapolation is done with $`1/L,1/L^2,1/L^3,ln(L)/L`$ with the coefficient of $`ln(L)/L`$ fixed as before. We show one example in Fig. 1. The fit to a third order polynomial in $`\alpha _V`$ is shown in Fig. 4 where the value of $`q^{}`$ was obtained in Ref. . The results are tabulated in Table 2.

## 5 SUMMARY

It was shown that the one-loop Wilson quark mass renormalization constant is reproducable with a combined extrapolation and statistical error of $`0.5\%`$ from large $`\beta `$ Monte-Carlo simulations of the quark propagator. Twisted boundary conditions ensured a perturbative expansion around the trivial vacuum and the absence of zero modes. The gauge fixing procedure was checked to reproduce the static quark self-energy using both the Polyakov loop and the static quark propagator.

Our preliminary result for the two-loop Wilson quark mass renormalization (for $`\kappa =0.1`$) is $`7(1)(0.5)`$. We have a crude result for the three-loop term (see Fig. 4), but it is poorly determined and sensitive to the details of the fit procedure. For now, we simply take it as an estimate of the systematic error in the determination of the two-loop term. A more detailed study of the systematic errors is now in progress including the infinite volume extrapolation, dependence on $`q^{}`$ and truncation of the series. The heavy quark mass renormalization for improved actions, such as clover and NRQCD, are also currently under investigation.

## ACKNOWLEDGEMENTS

I would like to thank A. Kronfeld, P. Lepage, P. Mackenzie, N. Shakespeare, J. Simone and H. Trottier for many useful discussions. Fermilab is operated by University Research Association, Inc. for the U. S. Department of Energy.


text:
# Small strange stars and marginally stable orbit in Newtonian theory

## Abstract

It is shown that for very rapidly rotating low mass strange stars the marginally stable orbit is located above the stellar surface. This effect is explained by the very important role of the oblateness of the rotating strange star. The comparison with some “academic” examples is presented. This feature is purely Newtonian in its nature and has nothing to do with relativistic marginally stable orbit. The effect is very large and cannot be treated in a perturbative way. It seems that strange stars as a very dense self-bound objects are the only possibility in Nature to represent these toy models.

Recently marginally stable orbits have been extensively studied as a crucial input of the most popular model explaining kilohertz quasi-periodic oscillations (QPOs) in X-ray binaries (see  for a review). Besides neutron stars as a possible source of these QPOs hypothetical strange stars have been considered. Strange stars, objects built of the mixture of u, d, s quarks, in contrast to neutron stars have no minimum mass; their size extends continuously from that of large nuclei (strangelets) up to the maximum mass configurations of about 2 solar masses (for a recent review, see ). It has been showed  that the marginally stable orbit is located above the surface of the strange star for relatively high mass even at very high rotation rates, which is the main difference with respect to neutron stars . It was suggested  that, for configurations close to the Keplerian one, defined as a mass-shedding limit at the equator, the very large deformation of the strange star leads to the significant increase of the radius of the marginally stable orbit, resulting in a gap between this orbit and the stellar surface. Of course in these considerations, performed for relatively massive ($`M>\mathrm{M}_{}`$) stars, the GR effects seem to play a dominant role. General relativity predicts that the radius of the marginally stable orbit for nonrotating star of mass $`M`$ is equal to $`8.8M/\mathrm{M}_{}`$ km, which is comparable to the radius of a strange star. It is well known that quadrupole moment of the star could also destabilize the orbital motion , but it seemed that this effect is not very large in real physical situations, at least as far as neutron stars are concerned. We show that for objects such as self-bound strange stars the deformations due to the rotation can be so large that the change in Newtonian gravitational potential could result in the existence of a marginally stable orbit, even when GR corrections are completely negligible.

The existence of a marginally stable orbit corresponds to the conditions $`\mathrm{d}V/\mathrm{d}r=0`$ and $`\mathrm{d}^2V/\mathrm{d}r^2=0`$ where $`V(r)`$ is the effective potential (for the radial part of motion) for a particle orbiting the central body at a fixed angular momentum per unit mass $`l`$. The resulting equations are discussed in standard textbooks in the case of the Schwarzschild geometry (see e.g. ). In first approximation, the source of the marginally stable orbit is the term in $`V`$ proportional to $`1/r^3`$ with the same sign as the gravitational potential (i.e. opposite sign to centrifugal force). In general relativity it is a natural result of constants of geodesic motion on circular orbits.

In Newtonian theory this term is generated by the quadrupole moment of the star $`Q`$; in the equatorial plane the gravitational force per unit mass of the orbiting particle is

$$f_g=\frac{GM}{r^2}\frac{3Q}{2r^4}.$$
(1)
With the above convention $`Q>0`$ for an oblate star and $`Q<0`$ for a prolate one. The relativistic definition of $`Q`$ can be found in Ref. .

One can take into account the effect of the oblateness and of the angular momentum of the star using standard methods, i.e. via expansions of the metric functions in powers of $`1/r^2`$. In the framework of GR this leads to the following formula for marginally stable orbit:

$$\begin{array}{ccc}\frac{R_{\mathrm{ms}}}{3R_\mathrm{g}}\hfill & =& \frac{1}{2}(1+\sqrt{1+q_{\mathrm{GR}}/6})\hfill \\ & & \left(\frac{2}{3}\right)^{3/2}j\sqrt{1+\frac{9}{16}(\sqrt{1+q_{\mathrm{GR}}/6}1)},\hfill \end{array}$$
(2)
where $`q_{\mathrm{GR}}=Qc^4/G^2M^3`$, $`j=Jc/GM^2`$, $`R_\mathrm{g}=2GM/c^2`$, $`J`$ being the total angular momentum and we linearized the dependence in $`j`$. In the case of $`Q=0`$ Eq. (2) reduces to the formula derived by Kluzniak and Wagoner .

In nonrelativistic limit we have

$$R_{\mathrm{ms}}\frac{3}{2}\left(R_\mathrm{g}+\sqrt{R_\mathrm{g}^2+\frac{2Q}{3M}}\right)\sqrt{\frac{3Q}{2M}}.$$
(3)
In Newtonian theory the marginally stable orbit can be determined exactly for some toy models. We will present here results for one- and two-dimensional cases.

In one dimension we consider a bar with given constant linear density $`\rho _l`$, i.e. mass of the bar is $`M=2\rho _lR`$, where $`R`$ is the “radius” of the bar equal to the half of its length. The gravitational attraction on the line which coincides with the bar at the distance $`r`$ from the bar center is equal to $`f_g=GM/(r^2R^2)`$. Of course we have to assume that this field does not change during orbital motion of the particle, i.e. the bar rotates with the same angular velocity as the orbiting particle. The radius of the marginally stable orbit is then

$$R_{\mathrm{ms}}^{\mathrm{bar}}=\sqrt{3}R1.73R.$$
(4)
In two dimensions the corresponding body is an infinitesimally thin disk of constant surface density $`\rho _s`$ with mass $`M=\pi \rho _sR^2`$. The gravitational force in the plane of the disk at the distance $`r`$ is then given by the hypergeometric function $`{}_{}{}^{2}F_{1}^{}`$

$$f_g=\frac{GM}{r^2}\frac{1}{(1+R/r)^{3/2}}{}_{}{}^{2}F_{1}^{}(\frac{3}{4},\frac{5}{4},2;q^2),$$
(5)
where $`q=2rR/(r^2+R^2)`$. The form of the effective potential \[the integral of Eq. (5) plus a centrifugal term\] results in the marginally stable orbit located at the radius

$$R_{\mathrm{ms}}^{\mathrm{disk}}1.2885R$$
(6)
The angular momentum of the particle at this orbit is then equal to $`l\sqrt{1.784GMR}`$ (per particle mass).

It should be mentioned that in these examples expansions of the gravitational potential in powers of $`(R/r)^2`$ would be completely unjustified. For example in the case of a disk the expansion of Eq. (5) leads to the formula:

$$f_g=\frac{GM}{r^2}(1+\frac{3}{8}(R/r)^2+\frac{15}{64}(R/r)^4+\mathrm{}).$$
(7)
Using only the second term (quadrupole order) we obtain the radius of marginally stable orbit

$$R_{\mathrm{ms}}^{\mathrm{disk}}[Q]=\sqrt{0.375}R0.61R$$
(8)
i.e more than twice smaller than the exact value.

In previous papers  the marginally stable orbit for rotating strange stars have been determined in the framework of General Relativity and the authors concluded that the main difference with respect to neutron stars is the gap between this orbit and the stellar surface for high rotation rates.

Rotating strange stars are very oblate for large rotation rates and the Keplerian configurations for a sequence of stars with given baryon number is achieved not due to the increase of the rotation frequency $`\nu _{\mathrm{rot}}`$, but due to the deformation of the star . It was suggested that this large oblateness is the true reason which cause the radius of the marginally stable orbit to be above the stellar surface.

In this paper we try to prove this statement by considering strange stars for which relativistic effects are very small. Thus we will focus on the results for small strange stars where the oblateness is crucial for the properties of these objects and in particular, is the source of the marginally stable orbit. The effects of General Relativity, which are very important in massive dense stars, are here negligible. However presented results has been obtained in the framework of GR using the numerical code developed by Gourgoulhon et al. .

We calculate rotating bare strange stars (without crust) using the strange matter model considered as a basic one in previous papers (and called SQM1 in ).

In Fig. 2 we present the equatorial radius of the star $`R_{\mathrm{eq}}`$ with fixed baryon mass $`M=0.005\mathrm{M}_{}`$ and the radius of marginally stable orbit. We choose such a small object in order to extract Newtonian effects - which for massive stars are contaminated (and for slowly rotating stars dominated) by relativistic ones. We see that for high rotation rates there exist a gap between stellar surface and the marginally stable orbit. The curve $`R_{\mathrm{ms}}^{\mathrm{disk}}`$ is defined by Eq. (6) with $`R=R_{\mathrm{eq}}`$. The curve $`R_{\mathrm{ms}}^{\mathrm{quad}}`$ corresponds to the quadrupole approximation of $`R_{\mathrm{ms}}`$ given by Eq. (2) with $`Q`$ calculated in the framework of GR \[Eq. (7) of \]. We see that this quadrupole approximation is wrong and as in the case of disk the resulting radius is about twice smaller than the true value. This conclusion confirms the results of Shibata and Sasaki , who found necessary to include higher multipole moments in the marginally stable orbit determination in the case of polytropic EOS. They concluded that the quadrupole moment is as important as the spin. In our model the effect of the rotational deformation of the star is larger because the matter is self-bound at very high density and the stars are more oblate. Due to these two reasons for fast rotating strange stars we are “closer” to the disk approximation than to the expansion up to the quadrupole term (see Fig. 2). The value of the quadrupole moment for small strange stars close to the Keplerian configuration is $`Q_{\mathrm{max}}=0.17MR^2`$ which is quite close to the “disk” value – $`0.25MR^2`$.

The behavior of $`R_{\mathrm{eq}}`$ and $`R_{\mathrm{ms}}`$ for a star rotating faster and faster, presented in Fig. 2, is typical for small mass stars self-bound at high density. For higher masses (larger than $`1.12\mathrm{M}_{}`$ in this model) the radius of marginally stable orbit is larger than $`R`$ at $`\nu _{\mathrm{rot}}=0`$ and for low rotation rates decreases due to the relativistic effects and then increases due to the oblateness (see  for figures and discussion). Of course the point at which $`R_{\mathrm{ms}}(\nu _{\mathrm{rot}})`$ gets larger than the equatorial radius $`R_{\mathrm{eq}}`$ depend on the mass of the star under consideration.

In the next figures we present the regions of stellar parameters and rotation rates for which marginally stable orbit is located above the stellar surface.

In Fig. 3 we show the boundaries resulting from the Keplerian frequency and from the condition $`R_{\mathrm{eq}}=R_{\mathrm{ms}}`$. We see that in $`MR`$ plane the region with $`R_{\mathrm{ms}}<R_{\mathrm{eq}}`$ looks rather small and even for low-mass stars there exist relatively large area in which $`R_{\mathrm{eq}}<R_{\mathrm{ms}}`$. It is also clearly visible that the stars with masses between $`1.12`$ and $`1.28\mathrm{M}_{}`$ enter the region with $`R_{\mathrm{ms}}<R_{\mathrm{eq}}`$ during slowing down or acceleration although for high and low rotation rates $`R_{\mathrm{ms}}>R_{\mathrm{eq}}`$. These stars for fixed baryon mass are represented by nearly horizontal (vertical) lines in top (bottom) panel of Fig. 3 (“nearly” because gravitational mass increases a little due to the rotation).

In Fig 4 we present limits on the frequency of marginally stable orbit taking into account the conditions $`R_{\mathrm{eq}}<R_{\mathrm{ms}}`$ which in principle can be supported by the observations of QPOs .

There exists a quite large region in the $`\nu _{\mathrm{ms}}`$-$`M`$ plane, where this condition is fulfilled. It should be mentioned that because $`\nu _{\mathrm{ms}}`$ in non-monotonic function of $`J`$ for a star with given baryon mass, one region in $`\nu _{\mathrm{ms}}`$-$`M`$ plane is not uniquely defined. Namely in double hatched area there exist stellar configurations with $`R_{\mathrm{eq}}<R_{\mathrm{ms}}`$ (slowly rotating) as well as stars for which $`R_{\mathrm{eq}}>R_{\mathrm{ms}}`$. For nonrotating stars we have here $`R_{\mathrm{eq}}<R_{\mathrm{ms}}`$ but if such a star rotates faster and faster it moves upwards, turns back at the solid line and moves downwards having $`R_{\mathrm{eq}}>R_{\mathrm{ms}}`$ as long as it do not leave the region between solid lines. One has similar situation (non-uniqueness) in the case of high mass stars ($`M>1.28\mathrm{M}_{}`$), which turn back at long-dashed line but the condition $`R_{\mathrm{eq}}<R_{\mathrm{ms}}`$ is fulfilled all the time. In non hatched region all stars have $`R_{\mathrm{eq}}<R_{\mathrm{ms}}`$.

From Fig. 5 we see that the oblateness of the star, the quadrupole moment $`Q`$ and the ratio of the kinetic to the potential energy $`T/W`$ for the rotating stellar configurations defined by the equality $`R_{\mathrm{eq}}=R_{\mathrm{ms}}`$ are highly correlated and almost linearly dependent. Although the frequency at which this condition is fulfilled is approximately constant for $`M<1\mathrm{M}_{}`$ (0.9-1 kHz, Fig 3) the stars with lower mass are much more oblate at this rotation rate, having proportionally larger $`T/W`$ and $`Q`$. For these configurations large value of $`R_{\mathrm{ms}}`$ results entirely from oblateness of the star. Albeit the role of the oblateness for the location of the marginally stable orbit is most pronouncing in the case of small strange stars close to the Keplerian velocity, also for massive stars rotating with relatively high frequency the deformation of the star is crucial for the existence of the marginally stable orbit.

In this paper we showed that the gap between innermost stable orbit and the surface of the star can be produced by the oblateness of the star in the framework of Newtonian theory. Due to the compactness and very high surface density strange stars are natural candidates for this effect to be a real possibility in Nature. This is due to two reasons. First: the stars are very oblate, and the gravitational force is far from the attraction by a point mass. Second: the stars are self-bound at very high density and in practice an oblate small strange star is quite close to a constant density disk. In contrast to bare strange stars, for oblate neutron stars (and also for strange stars with crust) the low-density region is relatively large making the radius of the star larger than $`R_{\mathrm{ms}}`$ (or in other words the radius of corresponding “disk” is much smaller than the stellar radius).

###### Acknowledgements.

We are very grateful to P. Haensel for careful reading of manuscript and helpful comments and suggestions. This research was partially supported by the KBN grants No. 2P03D.014.13. The numerical calculations have been performed on computers purchased thanks to a special grant from the SPM and SDU departments of CNRS. Main results reported in the present note were obtained during visit of one of the authors (J.L.Z) at DARC, Observatoire de Paris- Meudon, in July 2000, within the framework of the CNRS/PAN program Jumelage Astrophysique.


text:
# The Extended X-ray Halo of the Crab-like SNR G21.5-0.9

## 1 Introduction

Of the $`225`$ known supernova remnants (SNRs) in our galaxy no more than $`5\%`$ are classified as Crab-like remnants (Green green00 (2000)). Prominent members of this class of SNR, also known as plerions, include the Crab Nebula, CTB 87, and 3C58. The characteristic centre-filled radio and X-ray morphologies of the Crab-like SNRs is thought to be due to the presence of an active pulsar which powers a bright synchrotron nebula. In the case of the Crab nebula, 33 ms pulsations testify to the presence of the central spinning neutron star, although in many Crab-like systems there is no direct evidence for pulsed emission. Typically a Crab-like SNR exhibits a flat power-law spectrum in the radio regime which eventually steepens at shorter wavelengths to join smoothly to a similarly featureless power-law X-ray continuum. This non-thermal form for the radio to X-ray spectrum is a consequence of the continuous injection by the pulsar of high energy electrons, which suffer radiation and adiabatic losses as they diffuse through the nebula; however the details of the process are undoubtedly complex (e.g. Reynolds & Chanan reynolds84 (1984)). A feature of the Crab-like remnants, which distinguishes them from both shell-like and composite-type SNRs (the latter typically exhibiting a filled-centre X-ray morphology within a radio shell), is the lack of any evidence for an outer shell structure marking the progress of the blastwave from the original supernova explosion.

G21.5-0.9 is a SNR with many of the characteristics of a Crab-like remnant (e.g. Becker & Szymkowiak becker81 (1981); Asaoka & Koyama asaoka90 (1990)). To date there has been no detection of pulsed emission in either the radio or X-ray bands (e.g. Biggs & Lyne biggs96 (1996); Slane et al. slane00 (2000)). However, recent X-ray observations made by Chandra have pinpointed the probable location of the pulsar, namely a very compact central core on a scale of $`2^{\prime \prime }`$ embedded within a more extended ($`30^{\prime \prime }`$ radius) synchrotron nebula (Slane et al. slane00 (2000)). A faint extended halo was also detected which might correspond to the outer “shell” formed from the expanding ejecta and the passage of the supernova-driven blastwave (Slane et al. slane00 (2000)). In this interpretation G21.5-0.9 would be best described as a composite remnant. Here we report further recent X-ray observations of G21.5-0.9 made by XMM-Newton which help to clarify the nature of the extended X-ray halo in this object.

## 2 XMM-Newton observations

G21.5-0.9 was observed as a calibration target during orbits 60-62 and 64-65 of the XMM-Newton mission. Here we focus on measurements from orbit 60 (date: 2000/04/07; start time: 12.35 UT; end time 22.57 UT) in which the SNR was observed on-axis by the EPIC MOS and PN cameras (Turner et al. turner01 (2001); Strüder et al. struder01 (2001)) giving a total accumulated exposure time of $`30`$ ks. For these observations all three EPIC cameras were operated in the standard full-frame mode with the medium filter selected. Attitude information files are not yet available but all the indications are that a stable pointing was maintained for the duration of the observation111Our nominal position for the centre of G21.5-0.9 is (RA, Dec) = 18 hr 33 min 33.8 s, $`10^{}34^{}6^{\prime \prime }`$ (J2000); (l, b) = $`21.50^{}`$, $`0.88^{}`$..

The recorded events were screened with the XMM Science Analysis Software (SAS) to remove known hot pixels and other bad data and pre-processed using the latest CCD gain values. X-ray events corresponding to patterns 0–12 for the two MOS cameras (similar to grades 0-4 in ASCA) were used, whereas for the PN only pattern 0 events (single pixel events) were accepted. Investigation of the full-field count-rate revealed a number of flaring events in the non-cosmic component of the background during the observation. As a final step in the data screening we identified a period when the background was essentially “quiescent”, with full-field background rates in the MOS and PN cameras of $`2.5`$ and $`6`$ count/s. Selecting events only in this low-background period lead to effective exposure times of 13800 s and 11800 s for the MOS and PN cameras respectively.

## 3 Spatial Analysis

Given the observed soft X-ray cut-off in the spectrum of G21.5-0.9 (see §4), we utilize the 1–10 keV energy band in the spatial analysis. The raw image of G21.5-0.9 from the PN dataset binned into $`4^{\prime \prime }`$ pixels is shown in Fig. 1 (left panel). The bright synchrotron nebula at the centre of this SNR is the most prominent feature of the image, although the central core is clearly surrounded by a more extended low surface-brightness emission region. In addition a point source is evident $`2^{}`$ to the south-west of the centre of the remnant.

To avoid the complications due to the chip boundaries in the PN data we have based our spatial analysis on the data from the MOS cameras (since the $`10.9^{}\times 10.9^{}`$ field of view of the central chip in the MOS cameras comfortably encompasses the full extent of the SNR). Broad band images were constructed for both MOS 1 and 2 cameras utilizing $`1^{\prime \prime }`$ pixels. The MOS 1 and 2 images were then aligned (by eye) using the nebula core and the south-western point source as fiducial marks and co-added. Fig. 1 (right panel) shows the resultant image after applying a Gaussian smoothing filter with a width $`\sigma `$ = 2 pixels (FWHM of $`4.7^{\prime \prime }`$). At a qualitative level there is clearly excellent agreement between the images from the PN and MOS cameras. The bright centre of the synchrotron nebula (i.e. the region within $`r30^{\prime \prime }`$ of the centre) shows deviations from circular symmetry, in particular there is an indentation in the north-west quadrant. At a slightly lower level of surface brightness a curious “spur” is evident which appears initially to track northwards away from the core region but then at $`r100^{\prime \prime }`$ sweeps round in an arc into the north-eastern quadrant of the remnant. These spatial features match rather well to similar structure seen in the higher spatial resolution images from Chandra (Slane et al. slane00 (2000)). One of the strengths of the EPIC cameras on XMM-Newton is, of course, the sensitivity afforded to extended, relatively low-surface X-ray emission. This is evident in the images of Fig. 1 in which a low surface brightness “halo” is well delineated out to $`r150^{\prime \prime }`$. This halo exhibits a near perfect circular symmetry apart from the hint of a “linear termination” along its northern rim (note that out-of-time events, namely events recorded during the period when charge is transferred to the readout node, partly mask this feature in the PN image).

In order to investigate more quantitatively the spatial extent of the low surface brightness halo, we have calculated the azimuthally-averaged radial profile of the X-ray emission (Fig. 2). The surface brightness away from the centre of the synchrotron nebula drops rapidly out to $`r50^{\prime \prime }`$. A plateau region is then reached corresponding to the low surface brightness halo. A sharp further decline in the surface brightness then sets in at $`r130^{\prime \prime }`$, with the signal eventually falling below the detection threshold at $`r180^{\prime \prime }`$. In this paper we define the remnant’s outer radius to be $`150^{\prime \prime }`$ at which point the surface brightness is roughly one-fifth of the plateau value.

## 4 Spectral Analysis

We first examined the spectral properties of G21.5-0.9 by plotting the softness ratio image shown in Fig. 3, which shows that the core of the nebula is significantly harder than the low-surface brightness halo.

We have investigated the spectral characteristics of the core of G21.5-0.9 using data from both the PN and MOS cameras. The X-ray spectrum of the core of the nebula was obtained using a circular extraction region of radius $`48^{\prime \prime }`$ centred on the position of peak surface brightness. Background spectra were extracted from representative regions greater than $`200^{\prime \prime }`$ away from the SNR. The source spectra were binned to a minimum of 20 counts per spectral channel, in order to apply $`\chi ^2`$ minimisation techniques. The source and corresponding background spectra were then analysed using the xspec v11.01 and the most recent response matrices available from the EPIC team. The spectral modelling assumed a power-law continuum (photon spectral index, $`\mathrm{\Gamma }`$, and normalisation, $`A`$) subject to soft X-ray absorption in cool gas along the line of sight (with an equivalent hydrogen column density, $`N_\mathrm{H}`$).

Fig. 4 shows the observed count rate spectra for the PN, MOS 1 and MOS 2 detectors together with the corresponding best-fitting model spectra. The spectral response of the two MOS cameras is so similar that the MOS 1 and MOS 2 count rate spectra almost totally overlap. The figure illustrates that the observed count rate in the PN camera is roughly a factor of two higher than in a single MOS camera at 2 keV, rising to a factor of four by 8 keV. The best-fitting spectral parameters obtained for separate fits to the PN and MOS spectra are summarised in Table 1, where the errors are quoted at a 90% confidence level. The simple absorbed power-law model provides a very good fit to all three data sets. In addition the agreement between the spectral parameters derived independently from the PN, MOS 1 and MOS 2 datasets is excellent.

The spectral variations revealed in Fig. 3 were investigated by extracting spectra from the PN dataset from a set of annular regions centred on the point of peak surface brightness222In this analysis we excluded a $`16^{\prime \prime }`$ radius region centred on the south-western point source.. The resulting spectra were then fitted as before with a simple absorbed power-law model. Initially we allowed $`N_\mathrm{H}`$ to vary as a free parameter between the different spectral datasets but, in the event, excellent fits were obtained with $`N_\mathrm{H}`$ fixed at the value obtained for the core region (Table 1). The derived photon spectral index shows a steady increase with radius from a value of $`1.63\pm 0.04`$ at the centre to $`2.45\pm 0.06`$ at the outer edge of the low surface brightness halo (Fig. 5). We note that the rate of spectral softening with radius will be slightly underestimated due to spreading caused by the instrument PSF (see Fig.2); however the impact on the derived spectral indices is likely to be small ($`\mathrm{\Delta }\mathrm{\Gamma }<0.1`$). Fig. 6 further illustrates the variation of the X-ray spectrum from the core through to the outer halo region.

We next investigated whether there is any evidence for thermal emission from the halo component, focussing on the PN spectrum extracted from an outer ($`r=96^{\prime \prime }144^{\prime \prime }`$) annular region. As a first step we fitted a two component model, comprising a power-law and a (solar abundance, MEKAL) thermal component. The result was a modest improvement in the $`\chi ^2`$ with respect to the power-law model (see Table 2) which the F-test shows to be marginally significant for two additional parameters. However, this putative ($`1`$ keV) thermal component accounts for only $`1\%`$ of the observed flux. The lack of prominent line features in the X-ray spectrum of the outer halo is evident in Fig. 6 which shows the fitting residuals to a simple power-law spectral model. One possibility is that the line emission is suppressed due to the fact that shock-heated X-ray emitting gas in the halo of G21.5-0.9 has not had time to reach a state of ionization equilibrium. When we fit a non-equilibrium ionization model (the NEI model in xspec) we obtain the spectral parameters listed in Table 2. The NEI model provides only a rather poor fit to the data when the column density is constrained to the value measured for the SNR core. A rather better result is obtained when $`N_\mathrm{H}`$ is allowed to vary freely (the last entry in Table 2) but nevertheless the fit is not as good as that obtained for the absorbed power-law model.

We also examined the spectrum of the brightest part of the spur feature (located directly to the north of the core) but unfortunately the signal-to-noise ratio was too low to distinguish between thermal and power-law spectral forms.

The total integrated X-ray flux of G21.5-0.9 within a radius of $`144^{\prime \prime }`$, corrected for the line-of-sight absorption, is $`6.1\times 10^{11}`$ erg cm-2 s-1 in the 2–10 keV band. Assuming a distance of 5 kpc based on HI absorption measurements (Davelaar, Smith & Becker davelaar86 (1986)), this implies an X-ray luminosity of $`1.8\times 10^{35}\mathrm{erg}\mathrm{s}^1`$.

Finally we have briefly considered the spectral properties of the south-west point source. We obtained an on-source spectrum from the PN dataset using an extraction cell of $`16^{\prime \prime }`$ radius. A complementary background spectrum was taken (on the same CCD) using a cell of the same dimension sited equidistant from the centre of G21.5-0.9. Spectral fitting indicated a relatively hard spectrum ($`\mathrm{\Gamma }=1.61\pm 0.34`$ or $`\mathrm{kT}7`$ keV) and a column density of $`1.03\pm 0.22\times 10^{22}\mathrm{cm}^2`$. As the latter is about a factor two lower than the derived $`N_\mathrm{H}`$ of G21.5-0.9, we can conclude that this point source is a foreground object with respect to the SNR. If we estimated its distance as $`3`$ kpc (i.e. roughly half the distance to the SNR in line with the ratio of $`N_\mathrm{H}`$ for the two sources) then the inferred 2–10 keV luminosity is $`4\times 10^{32}\mathrm{erg}\mathrm{s}^1`$. This point source is positionally coincident with the emission-line star SS397.

## 5 X-ray Timing Analysis

We have carried out a preliminary search for pulsed X-ray emission using EPIC data from the full observation (i.e. including periods of relatively high background rate). The photon events in the 1–10 keV band from a region of $`8^{\prime \prime }`$ radius centred on the core of G21.5-0.9 were extracted for all three cameras. Next the event times were barycentric corrected using the standard SAS task “Barycen” and the latest available “Reconstructed Orbit” file. Power spectra were then calculated for the combined data from all three cameras and the PN dataset alone. No evidence was found for a significant periodic signal. In terms of the amplitude of an underlying sinusoidal signal, the detection threshold was $`5\times 10^2`$ ct/s in the frequency range $`5\times 10^4`$ – 0.17 Hz for the PN + MOS combination (the upper frequency limit being set by the frame time of the MOS CCDs). Corresponding values for the PN-only data were $`3\times 10^2`$ ct/s in the frequency range $`5\times 10^4`$ – 2 Hz. In terms of the pulsed fraction (of the core emission) the limits are roughly $`3.5\%`$ and $`5.5\%`$ respectively.

## 6 Discussion

In a recent paper reporting Chandra observations, Slane et al. (slane00 (2000)) consider the possibility that the low-surface brightness halo surrounding the central synchrotron nebula in G21.5-0.9 might be the shell formed from the ejecta and blastwave of the original supernova explosion. The present XMM-Newton observations help delineate the spatial extent and morphology of this component. Crucially the new spectral information demonstrates that the extended halo has a spectral form devoid of any significant line features. If the halo emission is thermal then the lack of line emission implies that the plasma is far from ionization equilibrium. Our NEI modelling confirms a low ionization state with $`n_\mathrm{e}t3\times 10^8\mathrm{cm}^3\mathrm{s}`$. The temperature derived for the continuum (bremsstrahlung) emission is $`45`$ keV which is rather hot, even for a very young SNR. Assuming the distance to G21.5-0.9 is 5 kpc (Davelaar, Smith & Becker davelaar86 (1986)) and taking a maximum shell expansion velocity of $`10,000\mathrm{km}\mathrm{s}^1`$, we can estimate the (minimum) elapse time since the shock heating of the bulk of the X-ray emitting gas in the outer halo to be $`100`$ yrs. Applying the above constraint on $`n_\mathrm{e}t`$ then gives $`n_\mathrm{e}<0.1\mathrm{cm}^3`$. However, the observed X-ray luminosity of the outer halo, if interpreted as bremsstrallung radiation, sets a requirement for a significantly higher electron density, $`n_\mathrm{e}0.5\mathrm{cm}^3`$. There would appear therefore to be considerable difficulties in interpreting the extended X-ray halo in G21.5-0.9 as a thermal shell.

There are known examples of SNR shells which exhibit non-thermal X-ray spectra due to local particle acceleration (e.g. SN 1006, Koyama et al. koyama95 (1995); G347.3-0.5, Slane et al. slane99 (1999)). However, in the case of G21.5-0.9 the lack of any limb brightening, the smooth transition in spectral index throughout the remnant and the remarkable circular symmetry, all suggest an interpretation of the outer halo as an extension of the central synchrotron nebula. We conclude therefore, on the basis of the present observations, that G21.5-0.9 is a true Crab-like system rather than a composite (centre-filled plus thermal shell) object.

The observed spectral softening with radius in G21.5-0.9 most likely reflects the impact of synchrotron radiation losses on very high energy electrons as they flow from the region of the termination shock around the pulsar to the edges of the nebula. For example, if the nebula magnetic field is $`0.4`$ mG (Slane et al. slane00 (2000) ), then the synchrotron lifetime of the $`20`$ TeV electrons producing 5 keV X-rays is only $`2`$ yrs. Similar spectral trends due to “synchrotron burn-off” are seen in other Crab-like SNR, such as 3C58 (Torii et al. torii00 (2000)). In fact the variation in photon spectral index from a value of 1.63 at the centre of the G21.5-0.9 to 2.45 at the edge of its halo is amazingly similar to the index range seen in the Crab nebula (Willingale et al. willingale01 (2001)). In the case of the Crab nebula the spatial and spectral distribution of the optical to hard X-ray continuum has been explained in terms a magnetohydrodynamic flow of plasma and magnetic flux in a particle dominated wind (e.g. Kennel & Coroniti kennel84 (1984)). For G21.5-0.9, the origin of halo which represents a “plateau” in the surface brightness extending from $`r50^{\prime \prime }130^{\prime \prime }`$ remains unclear, but could possibly be the result of an unusual magnetic field geometry (see below).

One interesting contrast is that the radio extent of the Crab nebula is roughly a factor four greater than that observed in X-rays, whereas for G21.5-0.9 the dimensions of the X-ray halo exceeds that of the radio synchrotron nebula (which has a major axis of $`90^{\prime \prime }`$; Fürst et al. furst88 (1988)) by a similar factor. Slane et al. (slane00 (2000)) note that no radio “shell” is detectable to a (1$`\sigma `$) threshold of $`4\times 10^{21}\mathrm{W}\mathrm{m}^2\mathrm{Hz}^2\mathrm{sr}^1`$ at 1 GHz. If we make the assumption that the ratio of the core to halo flux is the same in the X-ray and radio regimes, then we predict the average surface brightness of an underlying radio halo (between $`r=60^{\prime \prime }130^{\prime \prime }`$) to be roughly 5 times this threshold. The inference is that very deep radio observations may well reveal the presence of an extended radio halo in G21.5-0.9.

The differences in the spatial morphology of the Crab and G21.5-0.9 could be due to the viewing orientation. Specifically the near circular symmetry of G21.5-0.9 might be explained if our line of sight is reasonably aligned with the spin axis of the resident pulsar, so as to give a near face-on view of a putative inner torus. This hypothesis might also explain the non-detection of pulsed X-ray emission from the core of G21.5-0.9, but is hard to reconcile with the axisymmetric structure seen in high frequency radio maps and the near radial distribution of the magnetic field lines (Fürst et al. furst88 (1988)). In this setting the spur-like feature which emanates northwards from the core could a “plerionic” wisp associated with a complex inner magnetic field and/or jet structure.

###### Acknowledgements.

The results presented in this paper are based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). We thank the whole XMM-Newton team for the hard work and dedication which underlies the success of the mission.