format: text

text:
# Diffuse X-ray emission from the planetary nebula NGC 7009

## 1 Introduction

It has been commonly accepted that the shell structure of a planetary nebula (PN) is formed by the current fast stellar wind sweeping up the previous asymptotic giant branch (AGB) wind (e.g. Kwok, 1983). In this interacting-stellar-winds model, a PN will comprise a central cavity filled with shocked fast wind at temperatures of 10$`{}_{}{}^{6}10^8`$ K, a dense shell of swept-up AGB wind at 104 K, and an outer envelope of unperturbed expanding AGB wind. The shocked fast wind has a high temperature but does not emit appreciably in X-rays because its density is low. At the interface between the shocked fast wind and the dense nebular shell, heat conduction lowers the temperature and mass evaporation raises the density of the shocked fast wind, producing optimal conditions for X-ray emission (Weaver et al., 1977). Thus, soft X-ray emission from the shocked fast wind in a PN interior is expected to peak near the inner wall of the nebular shell.

The $`ROSAT`$ X-ray satellite had been used to search for diffuse X-ray emission from PNe with limited success (Guerrero et al., 2000). It is only recently that the Chandra X-ray Observatory has unambiguously detected diffuse X-ray emission in the central cavities of three young PNe: BD +303639 (Kastner et al., 2000), NGC 6543 (Chu et al., 2001), and NGC 7027 (Kastner et al., 2001). In the case of NGC 6543 where the X-ray emission is well-resolved, a clear limb-brightening is observed. The plasma temperatures implied by their X-ray spectra are a few $`10^6`$ K. Chandra observations of the evolved PN NGC 7293 (the Helix Nebula), on the other hand, does not show any diffuse X-ray emission. These results suggest that the hot gas content of a PN diminishes as the PN ages. To confirm this hypothesis, PNe at intermediate evolutionary stages need to be observed.

We have examined all pointed and serendipitous observations of PNe in the $`ROSAT`$ archive and found that NGC 7009, the Saturn Nebula, hosts an X-ray source with spectral characteristics of plasma emission. NGC 7009 is more evolved than the three PNe with diffuse X-ray emission because of its larger nebular size and dynamical age. Consequently, NGC 7009 is an excellent candidate not only for diffuse X-ray emission, but also as a PN at intermediate evolutionary stages. We have obtained X-ray observations of NGC 7009 with the XMM-Newton Observatory, and found that the X-ray emission is indeed extended. This paper reports our analysis of the diffuse X-ray emission from NGC 7009.

## 2 Observations

NGC 7009 was observed with the XMM-Newton Observatory in Revolution 255 on 2001 April 30 using the EPIC/MOS1, EPIC/MOS2, and EPIC/pn CCD cameras. The two EPIC/MOS cameras were operated in the Full-Frame Mode for a total exposure time of 39.4 ks, while the EPIC/pn camera was operated in the Extended Full-Frame Mode for a total exposure time of 31.6 ks. The EPIC/pn observations started about 4 ks later than the EPIC/MOS observations because of a larger overhead required by the EPIC/pn, but both observations concluded at the same time. For all observations, the Medium filter was used. At energies below 1 keV, the point spread function (PSF) of the EPIC instruments has a small core with FWHM $`6^{\prime \prime }`$ but extended wings with a half energy width of $`15^{\prime \prime }`$. For the EPIC/MOS and EPIC/pn CCDs, the pixel sizes are 1$`\stackrel{}{.}`$1 and 4$`\stackrel{}{.}`$1, and the energy resolutions are $`60`$ eV and $`100`$ eV at 1 keV, respectively (Dahlem 1999; Strüder et al. 2001).

We received the XMM pipeline products which included the raw and processed data of all instruments. These data were further processed using the XMM-Newton Science Analysis Software (SAS version 5.2.0) and the calibration files from the Calibration Access Layer available on 2002 January 24.

The event files were further screened to eliminate bad events, such as those due to charged particles, and periods of high background. For the EPIC/MOS observations, only events with CCD patterns 0–12 (similar to $`ASCA`$ grades 0–4) were selected; for the EPIC/pn observation, only events with CCD pattern 0 (single pixel events) were selected. To assess the background rate, we binned the counts over 25 s time intervals for each instrument in the 10–12 keV energy range, where the counts were dominated by the background. The time intervals with high background, i.e. count rates $`0.4`$ cnts s-1 for the EPIC/MOS or $`1.5`$ cnts s-1 for the EPIC/pn, were discarded. The resulting exposure times are 35.0 ks, 35.4 ks, and 26.2 ks for the EPIC/MOS1, EPIC/MOS2, and EPIC/pn observations, respectively. A total of $`480\pm 25`$ cnts, $`470\pm 25`$ cnts, and $`1,610\pm 45`$ cnts from NGC 7009 were detected in the EPIC/MOS1, EPIC/MOS2, and EPIC/pn observations, respectively.

## 3 Spatial Analysis

We have analyzed the spatial distribution of the X-ray emission from NGC 7009 and compared it to the optical nebula. As the data analysis is fairly complex, we present images at various stages of the processing in Figure 1. First, the event files from the three EPIC instruments are merged together to construct a raw EPIC image with a pixel size of 1$`\stackrel{}{.}`$5 in the 0.15–1.25 keV energy band (Fig. 1a). As shown in §4, this energy band includes all photons detected from NGC 7009. The raw EPIC image is then adaptively smoothed using Gaussian profiles with FWHM ranging from 5′′ to 17′′. A grey-scale presentation of this adaptively smoothed EPIC image is shown in Fig. 1b. The contours from this smoothed EPIC image are overplotted on an H$`\alpha `$ image of NGC 7009 taken with the Hubble Space Telescope (HST) Wide Field Planetary Camera 2 (WFPC2) and presented in Fig. 1c. Finally, to improve the spatial resolution, this smoothed EPIC image is deconvolved using the Lucy-Richardson algorithm and a theoretical EPIC PSF provided in the XMM-Newton Calibration Access Layer. The contours extracted from the PSF-deconvolved EPIC image are overplotted on the HST WFPC2 H$`\alpha `$ image in Fig. 1d.

The X-ray emission from NGC 7009 is clearly extended. Its spatial distribution in the adaptively smoothed EPIC image (Fig. 1b and 1c) can be described by an elliptical Gaussian with FWHM $`23^{\prime \prime }\times 14^{\prime \prime }`$ and a major axis along the position angle $`70^{}`$. There is an excellent correspondence between the X-ray emission and the nebular shell of NGC 7009. The HST WFPC2 H$`\alpha `$ image of NGC 7009 shows a bright $`25^{\prime \prime }\times 10^{\prime \prime }`$ inner shell surrounded by a $`25^{\prime \prime }\times 20^{\prime \prime }`$ envelope, with two ansae extending from the tips of the major axes of the inner elliptical shell out to $`25^{\prime \prime }`$ from the central star. The position and orientation of the 50% contour of the smoothed EPIC image matches almost perfectly the inner shell of the nebula (Fig. 1c), indicating that the X-ray emission peaks within the central cavity. The PSF-deconvolved EPIC image further shows that the majority of the X-ray emission originates from the central cavity. As the EPIC PSF has extended wings and our PSF-deconvolution is limited by the number of counts and the validity of the theoretical PSF adopted, it is most likely that all diffuse X-ray emission from NGC 7009 is confined within the inner shell, similar to what is observed in BD +303639 and NGC 6543. Furthermore, NGC 7009 has a sharp-edged inner shell similar to that of NGC 6543, which implies an expansion driven by the pressure of the hot gas in the central cavity (Chu et al., 2001). Therefore, we suggest that the hot, X-ray-emitting gas in NGC 7009 is also confined within its central cavity, but the detailed distribution of the hot gas cannot be determined from the currently available X-ray observations. In particular, it is not possible to address whether X-ray emission from its central star is present or not.

## 4 Spectral Analysis

Spectra of NGC 7009 have been extracted from the event files of the three XMM-Newton EPIC cameras separately because they have different instrumental response functions. A $`50\mathrm{}\times 40\mathrm{}`$ elliptical source aperture is used to include all X-ray emission detected from NGC 7009. An elliptical annular region exterior to and $`10`$ times as large as the source aperture is used to assess the background level. The background-subtracted spectra of NGC 7009 are shown in Fig. 2-left.

The EPIC/pn spectrum peaks at $`0.55`$ keV and drops abruptly to a faint plateau between 0.7 and 1.0 keV, but stays essentially flat from 0.5 keV down to 0.3 keV. The EPIC/MOS spectra also peak at $`0.55`$ keV, but have smaller numbers of counts because of lower sensitivities. The detailed differences between the EPIC/pn and EPIC/MOS spectra reflect the different instrumental responses among these cameras and the low S/N ratios in the latter spectra. Neither the EPIC/MOS spectra nor the EPIC/pn spectrum shows any detectable emission at energies greater than 1.25 keV.

The EPIC spectra of NGC 7009 are clearly inconsistent with a power-law or blackbody spectrum. Thus, the X-ray emission from NGC 7009 must be thermal plasma emission and the spectral peak at $``$ 0.55 keV corresponds to the He-like triplet of O vii at $`0.57`$ keV. We can model the observed spectra in order to determine the physical conditions of the X-ray-emitting plasma and the amount of foreground absorption. We have adopted the thin plasma emission model of Raymond & Smith (1977) and the nebular abundances determined from optical and UV spectra of NGC 7009. The abundances of He, C, N, O, Ne, and S relative to hydrogen by number have been measured to be 0.13, 2.9 10-4, 2.2 10-4, 4.3 10-4, 7.2 10-5, and 1.9 10-5, respectively (Kwitter & Henry 1998). Accordingly, we have adopted nebular abundances, relative to the solar values (Anders & Grevesse, 1989), of 1.33, 0.8, 1.9, 0.7, 0.6, and 1.2 for He, C, N, O, Ne, and S, respectively, and 1.0 for the other elements. We have assumed solar abundances for the foreground interstellar absorption, and adopted absorption cross-sections from Balucinska-Church & McCammon (1992).

The spectral fits are carried out by folding the Raymond-Smith model spectrum through the corresponding EPIC/MOS1, EPIC/MOS2, and EPIC/pn response matrices, and fitting the modeled spectra to the observed EPIC spectra in the 0.3–1.5 keV energy range using the $`\chi ^2`$ statistics. The best-fit model, overplotted on the EPIC/pn and EPIC/MOS spectra in Fig. 2-left, has a plasma temperature of $`T`$ = 1.8 106 K (or $`kT`$ = 0.152 keV), an absorption column density of $`N_\mathrm{H}`$ = 8.1 1020 cm-2, and a volume emission measure of 1.9 10$`{}_{}{}^{54}d_{}^{2}`$ cm-3 where $`d`$ is the distance to NGC 7009 in kpc111 The distance estimates to NGC 7009 range from 0.5 to 2.1 kpc. In the following we will adopt a distance of 1.2 kpc..

The goodness of the spectral fits is illustrated by the plot of reduced $`\chi ^2`$ of the fits as a function of $`kT`$ and $`N_\mathrm{H}`$ shown in Fig. 2-right. The 99% confidence contour spans $`T`$ = 1.6–2.0 106 K (or $`kT`$ = 0.14–0.17 keV), $`N_\mathrm{H}`$ = 2–17 1020 cm-2, and volume emission measure 1.5–2.4 10$`{}_{}{}^{54}d_{}^{2}`$ cm-3. Assuming that the H I column density is similar to the absorption column density and adopting a typical gas-to-dust ratio (Bohlin, Savage, & Drake 1978), the corresponding visual extinction is $`A_V`$ = $`0.44_{0.3}^{+0.5}`$ mag, and the logarithmic extinction at the H$`\beta `$ line is $`c_{\mathrm{H}\beta }`$ = $`0.15_{0.15}^{+0.22}`$. This extinction is consistent with other estimates of extinction ($`c_{\mathrm{H}\beta }=0.24\pm 0.06`$, Lame & Pogge 1996; $`c_{\mathrm{H}\beta }=0.10_{0.08}^{+0.06}`$, Rubin et al. 2001). The observed (absorbed) X-ray flux of NGC 7009 in the 0.3–1.5 keV energy range is 7.2 10-14 ergs cm-2 s-1, and the unabsorbed (intrinsic) X-ray flux is 1.7 10-13 ergs cm-2 s-1. The total X-ray luminosity in the 0.3–1.5 keV energy range is $`L_\mathrm{X}=`$ 2.0 10$`{}_{}{}^{31}d_{}^{2}`$ ergs s-1 = 3 1031 ergs s-1 for a distance of 1.2 kpc.

## 5 Discussion

XMM-Newton EPIC observations of NGC 7009 have shown that the diffuse X-ray emission is confined within the central cavity. The temperature of the X-ray-emitting gas, 1.8 106 K, is below that expected for the shocked fast wind from its central star which has a terminal velocity $`2,800`$ km s-1 (Cerruti-Sola & Perinotto, 1989; Iping & Sonneborn, 2001). Similar situation is observed in the other PNe with diffuse X-ray emission: the plasma temperatures implied by the X-ray spectra of BD +303639 and NGC 6543 are several times lower than the expected post-shock temperature for their fast stellar wind (Kastner et al., 2000; Chu et al., 2001). This temperature difference is expected if heat conduction is prevalent at the interface between the shocked fast wind and the dense nebular shell (Soker, 1994; Zhekov & Perinotto, 1996).

In order to determine the role played by the hot gas in the nebular dynamics of NGC 7009, it is possible to derive the rms density and pressure of the X-ray-emitting gas and compare them to those of the surrounding cool nebular shell. Assuming a prolate ellipsoidal central cavity, the volume occupied by the X-ray-emitting gas is $``$6.7 10$`{}_{}{}^{51}(d/1.2)_{}^{3}ϵ`$ cm3, where $`ϵ`$ is the filling factor. From this volume and the aforementioned volume emission measure we derive an rms density of 20 $`(d/1.2)^{1/2}ϵ^{1/2}`$ cm-3, and a total mass for the X-ray-emitting gas of 8 10$`{}_{}{}^{5}(d/1.2)_{}^{5/2}ϵ^{1/2}M_{}`$. This rms density and the plasma temperature of 1.8 106 K imply that the thermal pressure of the X-ray-emitting gas is $``$4 10$`{}_{}{}^{9}(d/1.2)_{}^{1/2}ϵ^{1/2}`$ dynes cm-2. The nebular shell surrounding the hot gas has a density of 5,000 cm-3, a temperature of 10,500 K (Hyung & Aller 1995; Rubin et al. 2001), and thus a thermal pressure of $``$10-8 dynes cm-2. The thermal pressure of the nebular shell is comparable or higher than that of the X-ray-emitting gas unless the filling factor $`ϵ`$ is $``$0.2.

The X-ray luminosity of NGC 7009 and the density and pressure of the hot gas in its interior are the lowest among the four PNe for wich diffuse X-ray emission has been detected thus far. These differences may be associated with the nebular evolution. NGC 7009 has a large nebular shell, 0.145$`\times `$0.06 pc in size, and dynamical age, $`1,700`$ yrs for an expansion velocity of 17 km s-1 (Bohigas, López, & Aguilar, 1994). Therefore, it appears to be in the most advanced evolutionary stage among these four PNe. At a more advanced evolutionary stage, NGC 7293 does not have any detectable diffuse X-ray emission (Guerrero et al. 2001). The Chandra and XMM-Newton observations of these five PNe suggest that PNe are brightest in X-rays when they are young and NGC 7009 may be at an evolutionary stage when the diffuse X-ray emission starts to decline. This evolutionary effect is expected as the stellar wind diminishes and the central cavity grows. More observations of PNe at different evolutionary stages are needed to evaluate the duration of the pressure-driven phase in PNe and its effects on the nebular dynamics. As the evolved PNe are expected to have low X-ray surface brightness, the high sensitivity of XMM-Newton is essential in detecting such PNe.

###### Acknowledgements.

We thank Steve Snowden for advice on the XMM-Newton data reduction, Guillermo García-Segura for enlightening discussion on the X-ray emission from PNe, and Rosina Iping, George Sonneborn, and Robert H. Rubin for useful discussion on the physical properties of NGC 7009. We also thank the referee, Joel Kastner, for his constructive comments. This research was supported by the NASA grant NAG 5-10042.


text:
# On the Theory of Magnetotransport in a Periodically Modulated Two-Dimensional Electron Gas

## Abstract

A semiclassical theory based on the Boltzmann transport equation for a two-dimensional electron gas modulated along one direction with weak electrostatic or magnetic modulations is proposed. It is shown that oscillations of the magnetoresistivity $`\rho _{||}`$ corresponding to the current driven along the modulation lines observed at moderately low magnetic fields, can be explained as classical geometric resonances reflecting the commensurability of the period of spatial modulations and the cyclotron radius of electrons.

Theoretical studies of magnetotransport properties of a two-dimensional electron gas (2DEG) spatially modulated along a certain direction due to periodic electrostatic or magnetostatic fields are lasting more than a decade. At present the theory of magnetotransport in modulated 2DEG is well developed and the most of effects observed in such systems at low magnetic fields have been explained both by quantum mechanical (in a semiclassical limit) one ; two ; three ; four ; five ; six ; seven  and classical eight ; nine ; ten ; eleven ; twelve  transport calculations, giving consistent results. One of a few exceptions is an effect of oscillations of the resistivity component $`\rho _{||}`$ which corresponds to the current driven in parallel with the lines of modulation. These oscillations were observed along with the commensurability oscillations of another resistivity component $`\rho _{}`$ corresponding to the current driven across the modulation lines, and have the same period and a phase opposite to those of $`\rho _{}`$.

The oscillations of $`\rho _{||}`$ had been immediately explained as an effect of a pure quantum nature which originates from the quantum oscillations of the electron density of states (DOS) in the applied magnetic field $`𝐁`$ one ; three . However, actual results obtained in these works one ; three  are not consistent with these conclusions. This can be easily demonstrated, basing on the results of Ref.three  which are reproduced here (Fig.1). It is clearly seen in Fig.1 that both $`\rho _{}`$ and $`\rho _{||}`$ exhibit oscillations of the same period at low (unquantizing) magnetic fields where quantum oscillations of DOS are not resolved, and only effects controlled with classical mechanisms could survive.

It is very well known that periods of quantum oscillations and semiclassical commensurability oscillations (geometric resonance) are different. Their ratio equals $`2\pi /\lambda k_F`$ where $`\lambda `$ is the space period of modulations, and $`k_F`$ is the Fermi wave vector of electrons thirteen . Therefore, the observed coincidence on periods of the low-field Weiss oscillations of the resistivity component $`\rho _{}`$ and the weaker anti-phase oscillations of $`\rho _{||}`$ gives an evidence that for both resistivity components these oscillations are geometric oscillations. Such oscillations arise in modulated 2DEG due to the periodic reproduction of commensurability between the diameters of electron cyclotron orbits and the period of applied modulations $`\lambda .`$ Consequently, we can treat the low-field oscillations of $`\rho _{||}`$ as an effect basically following from classical mechanisms, and develop a theory of this effect based on Boltzmann transport equation. The results for the magnetic field dependence of $`\rho _{||}`$ given by the classical transport calculations have to be in agreement with the results obtained within the quantum mechanical approach in the low-field limit, as well as in the case of Weiss oscillations of the magnetoresistivity component $`\rho _{}.`$

The existing classical theory of magnetotransport in modulated 2DEG fails to provide a description of the low-field oscillations of $`\rho _{||},`$ and concludes that this magnetoresistivity component is not affected by the modulations fourteen . This gives a motivation for the present work. Its purpose is to demonstrate that the most important characteristic features of the low-field oscillations of the resistivity $`\rho _{||}`$ can be qualitatively reproduced within a semiclassical transport theory. We also briefly analyze where the earlier semiclassical theory might have erred. To simplify the following calculations the effects of anisotropy in electron scattering are neglected, and the relaxation time approximation is used. It is also assumed that the external magnetic field is moderately weak so that the electrons cyclotron radius $`R`$ is considerably smaller than their mean free path $`l`$ but larger that the period of modulations $`\lambda ,`$ and $`R>>\sqrt{l\lambda }`$ which provides preferred conditions for observation of commensurability oscillations of transport coefficients of the 2DEG. To facilitate a comparison of the present results with the previous work we intentionally use a theoretical formalism as close as possible to that adopted in the current semiclassical theory eight ; nine ; ten ; eleven .

We consider first electrostatic modulation with a single harmonic of period $`\lambda =2\pi /g`$ along the $`y`$ direction given by $`\mathrm{\Delta }E(y)=dV(y)/dy.`$ The screened modulation potential $`V(y)`$ is parametrized as $`eV(y)=\epsilon E_F\mathrm{sin}(gy)`$ where $`E_F`$ is the Fermi energy of the 2DEG. We examine weak modulations, so that $`|ϵgl|<<1.`$ Electron transport coefficients for this case were first calculated by Beenakker eight . He started from a linearized Boltzmann transport equation for the electron distribution function $`\mathrm{\Phi }(y,\psi )`$ of the form

$$D[\mathrm{\Phi }]+C[\mathrm{\Phi }]=𝐄𝐯$$
(1)
where $`\psi `$ is the angular coordinate of the electron cyclotron orbit; $`𝐄`$ is a homogeneous electric field applied to the 2DEG alongside with the modulating field $`\mathrm{\Delta }E(y);`$ the collision term $`C[\mathrm{\Phi }]`$ is written in the relaxation time approximation with the relaxation towards the local equilibrium distribution, namely:

$$C[\mathrm{\Phi }]=\frac{1}{\tau }\left(\mathrm{\Phi }(y,\psi )\frac{1}{2\pi }_0^{2\pi }\mathrm{\Phi }(y,\psi )𝑑\psi \right),$$
(2)
and the drift term in the case of electrostatic modulations applied along the $`y`$ direction equals

$$D[\mathrm{\Phi }]=v(y)\mathrm{sin}\psi \frac{\mathrm{\Phi }}{y}+(v^{}(y)\mathrm{cos}\psi +\mathrm{\Omega })\frac{\mathrm{\Phi }}{\psi }.$$
(3)
Here, $`\mathrm{\Omega }`$ is the electron cyclotron frequency and the electron velocity vector $`𝐯(y)`$ has a direction $`𝐮(\psi )=(\mathrm{cos}\psi ,\mathrm{sin}\psi ),`$ and the magnitude $`v(y)=v_F\sqrt{1+ϵ\mathrm{sin}(gy)}`$ where $`v_F`$ is the Fermi velocity in the unmodulated 2DEG. The Eq.(1) with the collision and drift terms of the form (2) and (3) agrees with transport equations of later papers \[9–11\]. The derivation of this equation is presented in detail in earlier works (see e.g. ), and we omit it for brevity.

The electron current density in the 2DEG modulated along the $`y`$ direction also depends on $`y`$ and can be written as

$$𝐣(y)=D_0e^2_0^{2\pi }\frac{d\psi }{2\pi }𝐯(y,\psi )\mathrm{\Phi }(y,\psi ).$$
(4)
Here, $`D_0=m/\pi \mathrm{}^2`$ is the electron DOS on their Fermi surface; $`m,e`$ are the effective mass and the charge of an electron, respectively.

Following the usual way  we write $`\mathrm{\Phi }(y,\psi )`$ as:

$$\mathrm{\Phi }(y,\psi )=\mathrm{\Phi }_0(\psi )\frac{v(y)}{v_F}+\rho _0\tau \chi (y,\psi )$$
(5)
where $`\rho _0`$ is Drude resistivity, and $`\tau `$ is the relaxation time. The homogeneous distribution function $`\mathrm{\Phi }_0(\psi )=\rho _0\tau 𝐯_0𝐣_0`$ describes the linear response of the 2DEG to the field $`𝐄`$ in the absence of modulations, and the function $`\chi (y,\psi )`$ satisfies the transport equation:

$$D[\chi ]+C[\chi ]=v^{}(y)v(y)j_{0y}.$$
(6)
Here, as well as before, $`𝐣_0`$ is the current density for the unmodulated 2DEG.

We have to remark that these Eqs.(1) and (6) include some inaccuracies. Namely, the effect of electrostatic modulations is taken into account simply replacing $`v_F`$ by $`v(y)`$ which is defined above. This is an oversimplification. For systematic consideration we should explicitly introduce into the transport equation an additional term which describes the ’electrochemical’ field arising due to inhomogeneity of electron density caused by modulations. As a result the electric field $`𝐄`$ has to be replaced by $`𝐄(1/e)\mu `$ where $`\mu (y)`$ is the correction to the chemical potential of electrons caused by modulations. The latter provides the local equilibrium of the 2DEG with a stationary electron density $`(n/t=0).`$ The value of this correction $`\mu (y)`$ is determined with the continuity equation which in this case has the form $`𝐣=0.`$

The above inaccuracies do not change the main approximation for the distribution function $`\chi (y,\psi )`$ but they can influence next terms in the expansion of $`\chi (y,\psi )`$ in powers of a small parameter $`(\mathrm{\Omega }\tau )^1`$, as we show below. Therefore, using Eq.(6) in further calculations we have to keep only those terms in the expansion of $`\chi (y,\psi )`$ which give for the current density results consistent with the continuity equation, and all remaining terms in this expansion have to be ignored.

To proceed we expand $`\chi (y,\psi )`$ in a Fourier series in the spatial variable $`y`$ which gives the system of differential equations for the Fourier components.

$`{\displaystyle \frac{d\chi _0}{d\psi }}+{\displaystyle \frac{1}{\mathrm{\Omega }}}C[\chi _0]`$ (7)
$`=`$ $`{\displaystyle \frac{gRϵ}{4}}\mathrm{sin}\psi (\chi _1+\chi _1){\displaystyle \frac{gRϵ}{4}}\mathrm{cos}\psi \left({\displaystyle \frac{d\chi _1}{d\psi }}+{\displaystyle \frac{d\chi _1}{d\psi }}\right);`$
$`\pm `$ $`igR\mathrm{sin}\psi \chi _{\pm 1}+{\displaystyle \frac{d\chi _{\pm 1}}{d\psi }}+{\displaystyle \frac{1}{\mathrm{\Omega }}}C[\chi _{\pm 1}]`$
$`=`$ $`{\displaystyle \frac{v_F}{4}}gRϵj_{0y}+{\displaystyle \frac{gRϵ}{2}}\mathrm{sin}\psi \chi _{\pm 2}`$
$``$ $`{\displaystyle \frac{gRϵ}{4}}\mathrm{cos}\psi \left({\displaystyle \frac{d\chi _0}{d\psi }}+{\displaystyle \frac{d\chi _{\pm 2}}{d\psi }}\right).`$
The equations determining Fourier components $`\chi _{\pm n}(\psi )`$ for the $`n2`$ are given by

$`\pm `$ $`gnR\mathrm{sin}\psi \chi _{\pm n}+{\displaystyle \frac{d\chi _{\pm n}}{d\psi }}+{\displaystyle \frac{1}{\mathrm{\Omega }}}C[\chi _{\pm n}]`$ (8)
$`=`$ $`{\displaystyle \frac{gRϵ}{4}}\mathrm{sin}\psi [n(\chi _{\pm n+1}\chi _{\pm n1})+\chi _{\pm n+1}+\chi _{\pm n1}]`$
$`=`$ $`{\displaystyle \frac{gRϵ}{4}}\mathrm{cos}\psi \left({\displaystyle \frac{d\chi _{\pm n+1}}{d\psi }}+{\displaystyle \frac{d\chi _{\pm n1}}{d\psi }}\right).`$
Solving these equations, and keeping the terms of the order or larger than $`(\epsilon gR)^2`$ we arrive at the following approximation for the distribution function $`\chi (y,\psi ):`$

$`\chi (y,\psi )={\displaystyle \frac{v_F}{2}}ϵglj_0^yQ\left(\mathrm{cos}\psi {\displaystyle \frac{sin\psi }{\mathrm{\Omega }\tau }}\right)`$
$`\times `$ $`\{\mathrm{sin}(gR\mathrm{cos}\psi +gy){\displaystyle \frac{1}{2}}ϵgR`$
$`\times `$ $`[\mathrm{cos}(gR\mathrm{cos}\psi )\mathrm{cos}^2gy{\displaystyle \frac{1}{2}}\mathrm{sin}(gR\mathrm{cos}\psi )\mathrm{sin}2gy]\}`$
where

$$Q=\frac{J_0(gR)}{1J_0^2(gR)},$$
(10)
and $`J_0(gR)`$ is the Bessel function. Using the obtained distribution function we can easily calculate the electron current density given by Eq.(4).

Taking into account only greatest terms in the expansion of $`\chi (y,\psi )`$ in powers of small parameter $`(\mathrm{\Omega }\tau )^1`$ we obtain that only $`j_x,`$ component gets a correction due to the modulations along the $`y`$ direction, whereas the component $`j_y`$ remains equal to $`j_y^0`$ and does not depend on spatial coordinates. This agrees with the continuity equation $`𝐣=0`$ which is necessary to obtain correct results for electron transport coefficients in modulated 2DEG.

To proceed we need to define effective (averaged over the period of modulations) transport coefficients. Here, we define the effective conductivity tensor $`\sigma _{eff}`$ as follows:

$$𝐣<𝐣(y)>\frac{g}{2\pi }_0^{2\pi /g}𝐣(y)𝑑y\sigma _{eff}𝐄.$$
(11)
To justify the adopted definition (11) we point out that the expressions for transport coefficients obtained either with quantum mechanical or classical calculations have to be consistent at low magnetic fields. Quantum mechanical calculations of the magnetoresistivity three ; four ; five ; six ; seven  give the expression which passes over to a classical conductivity tensor averaged over the period of modulations. Therefore, the latter is an accurate semiclassical analog of the conductivity calculated within the proper quantum mechanical approach, and our definition of $`\sigma _{eff}`$ agrees with that of Ref. fifteen . The same definition was used previously by Mirlin and Wolfle eleven .

As a result we obtain:

$$\sigma _{eff}^{xx}=\frac{\sigma _0}{1+(\mathrm{\Omega }\tau )^2}+\frac{1}{4}(\epsilon gR)^2\sigma _0\frac{J_0^2(gR)}{1J_0^2(gR)}$$
(12)
where $`\sigma _0=1/\rho _0`$ is the Drude conductivity of the electron system. Other components of the effective conductivity are not influenced due to modulations within our approximation. The second term in Eq. (12) arises due to modulations. This term exhibits oscillations whose origin is explained below.

It is known that in the presence of both electric and magnetic fields centers of electron cyclotron orbits drift in parallel with the cross product $`𝐄\times 𝐁.`$ The speed of this motion $`v_{dr}`$ equals $`E/B`$ which can be easily shown starting from the equations of motion of an electron subject to these fields. This results in one-dimensional diffusion of electrons along the direction of the drift. The latter brings extra contributions to the current density and changes the conductivity of the system. Here, electrostatic modulation $`\mathrm{\Delta }E(y)`$ is applied along the $`y`$ direction therefore guiding centers are drifting along the $`x`$ axis, and the electron diffusion caused by the drift reveals itself as a corection to $`\sigma _{xx}`$ component of the magnetoconductivity tensor. The second term in Eq.(12) represents this correction.

To estimate $`v_{dr}`$ for the inhomogeneous electric field $`\mathrm{\Delta }E(y)`$ we have to average $`\mathrm{\Delta }E(y)`$ over the cyclotron orbit. The mean value is mostly determined with contributions from the vicinities of two stationary points at the orbit where electron is moving in parallel with the modulation equipotential lines (along the $`x`$ direction), and the field $`\mathrm{\Delta }E(y)`$ is varying slowly eight . Correspondingly, $`v_{dr}^2`$ is enhanced when $`E(y)`$ takes on values of the same sign at both stationary points, and reduced otherwise. So the diffusion correction in the conductivity $`\sigma _{eff}^{xx}`$ which is proportional to $`<v_{dr}^2>`$ ehxibits an oscillatory dependence on the magnetic field, an these oscillations are commensurability oscillations due to geometric resonances in the guiding center drift. This physical mechanism was proposed to explain semiclassical Weiss oscillations in the magnetoresistivity $`\rho _{}`$ eight . The present result twelve  shows that the same mechanism controls semiclassical oscillations of $`\rho _{||}`$ (see Eqs.(13) ans (14) below).

The effective magnetoresistivity tensor is defined here as the inverse of the effective conductivity introduced by Eq.(11), $`\rho _{eff}=\sigma _{eff}^1`$. For the current driven across the modulation lines, the corresponding resistivity is:

$$\rho _{}=\rho _{yy}=\rho _0\left\{1+\frac{1}{4}(ϵgl)^2\frac{J_0^2(gR)}{1J_0^2(gR)}\right\}.$$
(13)
Assuming that the current flows along the modulation lines we get:

$$\rho _{||}=\rho _{xx}=\rho _0\left\{1\frac{1}{4}(ϵgR)^2\frac{J_0^2(gR)}{1J_0^2(gR)}\right\}.$$
(14)
For moderately weak magnetic fields $`(gR>>1)`$ the obtained results (13) and (14) describe oscillations of both magnetoresistivity components periodic in the inversed magnetic field magnitude. The oscillations of $`\rho _{}`$ and $`\rho _{||}`$ have the same period in $`1/B`$ and the opposite phases which corroborates the experimental results of one . The amplitude of the oscillations of $`\rho _{||}`$ is considerably smaller than that of $`\rho _{},`$ and this also agrees with the experiments of Ref. one , and with the results of numerical quantum mechanical calculations of Ref. three  in the limit of a weak magnetic field (see Fig.1 above). The result for the resistivity $`\rho _{}`$ as well agrees with the corresponding results of Refs. eight ; nine ; ten ; eleven  obtained within a classical magnetotransport theory.

However, the expression (14) for $`\rho _{||}`$ differs from the well known result of the current semiclassical theory. To analyze this discrepancy we now calculate the current density taking into account next terms in the expansion of the distribution function (9) in powers of $`(\mathrm{\Omega }\tau )^1.`$ Keeping terms no less than $`(ϵgR/\mathrm{\Omega }\tau )^2`$ we obtain that now all components of the conductivity get grating-induced corrections:

$`\delta \sigma _{xx}`$ $`=`$ $`{\displaystyle \frac{\sigma _0(\mathrm{\Omega }\tau )^2}{1+(\mathrm{\Omega }\tau )^2}}\alpha (y);`$
$`\delta \sigma _{xy}`$ $`=`$ $`\delta \sigma _{yx}={\displaystyle \frac{1}{(\mathrm{\Omega }\tau )}}\delta \sigma _{xx};`$ (15)
$`\delta \sigma _{yy}`$ $`=`$ $`{\displaystyle \frac{1}{(\mathrm{\Omega }\tau )^2}}\delta \sigma _{xx}.`$
Here, $`\alpha (y)`$ is the $`y`$ dependent factor of the order of $`(ϵgR)^2`$. It follows from Eq.(15) that both components of the current density depend on $`y`$ and continuity equation is violated. On these grounds we conclude that only main approximation for $`\chi (y,\psi )`$ can be used in magnetotransport calculations, since all remaining terms in the expansion of $`\chi (y,\psi )`$ include mistake arising due to original inaccuracies described before. This gives an extra justification to our results (13) and (14) for magnetoresistivities.

At the same time the expressions (15) enable us to show the origin of the mistake in the current semiclassical theory of Refs. eight ; nine ; ten ; eleven . With some formal transformations of the transport Eq.(6) we can present $`<\alpha (y)>`$ in the form:

$$<\alpha (y)>=\frac{1}{2\pi }_0^{2\pi }<v(y)\mathrm{sin}\psi G(y,\psi )>𝑑\psi $$
(16)
where $`G(y,\psi )`$ satisfies the equation:

$$D[G]+C[G]=\frac{2v^{}(y)v(y)}{v_F^2}.$$
(17)
This gives for $`\sigma _{eff}`$ and $`\rho _{eff}`$ expressions absolutely consistent with the existing semiclassical theory eight ; nine ; ten ; eleven ; twelve . Hence, these expressions could be obtained starting from the incorrect expression for the distribution function $`\chi (y,\psi )`$ which includes terms influenced due to inaccuracies in the transport Eqs.(1) and (6). This gives grounds to seriously doubt the results of earlier works eight ; nine ; ten ; eleven ; twelve  especially those concerning $`\rho _{||}.`$

So, the present analysis shows that simplified transport Eq.(1) can be successfully used in calculations of the main term in the expansions of transport coefficients in powers of $`(\mathrm{\Omega }\tau )^1.`$ As for next terms in these expansions, we need to modify the transport Eq.(1) (both drift and collision terms) to get them. For that we need to carry out consistent and systematic consideration of effects of internal electrochemical field arising due to grating-induced inhomogeneity of electron density. This is important for redistribution of the electron density at the presence of modulations provides the local equilibrium of the system sixteen .

Considering magnetic modulations we arrive at similar results. Assuming that the magnitude of the magnetic field gets a weak spatially dependent correction $`\mathrm{\Delta }B(y)=\mathrm{\Delta }B\mathrm{sin}(gy)(\mathrm{\Delta }B/B<<1)`$ we can derive the transport equation for the contribution to the distribution function of electrons $`\chi (y,\psi )`$ which arises due to the modulations:

$$D[\chi ]+C[\chi ]=\mathrm{\Delta }\mathrm{\Omega }(y)(v_{0y}j_x^0v_{0x}j_y^0)$$
(18)
where $`\mathrm{\Delta }\mathrm{\Omega }(y)`$ is the correction to the cyclotron frequency due to modulations, and the drift term equals:

$$D[\chi ]=v_{0y}\frac{\chi }{y}+(\mathrm{\Omega }+\mathrm{\Delta }\mathrm{\Omega }(y))\frac{\chi }{\psi }.$$
(19)
Writing out the transport Eq.(18) we have taken into account that the magnetic modulations do not change the magnitude of the electron velocity at the Fermi surface, so that $`v(y)=v_F.`$ Due to the same reason the first term in the expression (5) now does not depend on spatial coordinates and equals to the electron distribution function for the unmodulated 2DEG.

Solving Eq.(18) and keeping terms no less than $`(\mathrm{\Delta }B/B)^2`$ in the solution, we calculate the effective conductivity of the 2DEG at the presence of magnetic modulations following the way used before to analyze the magnetoresistivities of the 2DEG modulated by a weak electrostatic field. As a result we obtain that $`\sigma _{xx}`$ component of the effective conductivity tensor gets a correction due to the modulation:

$$\delta \sigma _{xx}=\sigma _0\left(\frac{\mathrm{\Delta }B}{B}\right)^2\frac{J_1^2(gR)}{1J_1^2(gR)}$$
(20)
where $`J_1(gR)`$ is the Bessel function. As in the case of electric modulations this correction originates from a guiding centers drift. Using Eq.(20) we easily arrive at the expressions for the magnetoresistivities $`\rho _{}`$ and $`\rho _{||}`$ for the current driven across and along the modulation lines:

$`\rho _{}`$ $`=`$ $`\rho _0\left\{1+(\mathrm{\Delta }\mathrm{\Omega }\tau )^2{\displaystyle \frac{J_1^2(gR)}{1J_0^2(gR)}}\right\};`$ (21)
$`\rho _{||}`$ $`=`$ $`\rho _0\left\{1\left(\mathrm{\Delta }\mathrm{\Omega }\tau {\displaystyle \frac{R}{l}}\right)^2{\displaystyle \frac{J_1^2(gR)}{1J_0^2(gR)}}\right\}.`$ (22)
These results describe Weiss commensurability oscillations in the limit of weak magnetic fields $`(gR>>1).`$ As well as for electrostatic modulations, the expression (22) reports weaker (due to the small factor $`(R/l)^2)`$ oscillations of $`\rho _{||}`$ whose phase is opposite with respect to that of commensurability oscillations of $`\rho _{}.`$

Finally, the novel result of the paper is that it gives a qualitative explanation of the low-field oscillations of the magnetoresistivity component $`\rho _{||}`$ in the 2DEG modulated along one direction within a semiclassical approach. It is shown here that these oscillations of $`\rho _{||}`$ at low magnetic fields are commensurability oscillations. At low temperatures when the quantum oscillations of the electron DOS at the Fermi surface are resolved, Shubnikov-de Haas oscillations can be superimposed upon the geometric oscillations of the magnetoresistivity. However, this does not change a classical nature of the effect itself.

Acknowledgments: I thank G.M. Zimbovsky for help with the manuscript.


text:
# Statistical Physics in a Finite Volume with Absolute Conservation Laws

## 1 Introduction

Results from a canonical ensemble that obeys exact conservation of additive charges and isospin can be useful in a variety of applications. For example, the unexpected isospin imbalances of CENTAURO events  could be explained by employing a Gaussian thermal source and restricting the pion emission to isoscalar pairs only. The pion multiplicity distributions become broader as compared to those of independent pion emission, thus resulting in larger fluctuations . This finding could be further substantiated by replacing the Gaussian thermal source by a canonical ensemble and the ad-hoc assumption of isoscalar pairs emission by conservation of total isospin and its projection.

Another potential area of application is the modeling of relativistic heavy ion collisions. Here, hydrodynamical models are used to describe the early stages of the collision, when the density is so high that the particles behave collectively as a fluid. As the collision progresses the density falls, the system breaks up, and hydrodynamics is not applicable anymore. Therefore, the break-up and later stages of the collision have to be described by other means such as, for example, microscopic transport models . When switching to a microscopic transport model from a macroscopic hydrodynamical prescription, particles must be formed from the fluid such that momentum and energy as well as isospin and all other charges are conserved. Currently, most schemes for generating particles are based on the grand canonical ensemble. However, exact conservation within each event and within each computational cell is essential for studies on fluctuations and balance functions .

Partition functions for canonical ensembles conserving charges can be written quite easily, but they are very difficult to compute for large systems because the number of partitions that have to be summed over increases rapidly . This difficulty can be overcome by rewriting the partition function as a recursion relation. These recursion relations were applied to ensembles of classical particles in  and to Fermi systems in .

In this work we add total isospin conservation to the canonical ensemble and derive the appropriate recursion relations. Furthermore, isospin fluctuations in the canonical ensembles, which, for example, are relevant for CENTAURO event studies, are calculated. A Monte Carlo algorithm that employs these recursion relations to generate particles is part of ongoing work and will be presented in the future. As our focus is the development of the computational technique, we show sample calculations for simple systems rather than comparisons to experimental results.

## 2 Recursion Relations with Classical Statistics

### 2.1 Partition Functions

The partition function for a canonical ensemble of $`A`$ particles conserving a charge $`Q`$, where $`Q`$ could be a vector containing several charges, is defined as

$$Z_{A,Q}=\underset{\alpha }{}\alpha |e^{\beta H}|\alpha ,$$
(1)
where $`\alpha `$ labels some quantum state with the appropriate quantum numbers, $`\beta `$ is the inverse temperature, and $`H`$ is the Hamiltonian of the system. Using the single particle partition function

$$\omega _k=\underset{i}{}\mathrm{exp}(ϵ_i^{(k)}/T),$$
(2)
which is simply a sum of Boltzmann weights over all available states $`i`$, the partition function can be written in a more tractable form as a sum over partitions subject to particle number and charge conservation

$$Z_{A,Q}=\underset{{\scriptscriptstyle \nu _ka_k}=A;{\scriptscriptstyle \nu _kq_k}=Q}{}\underset{k=1}{\overset{N}{}}\frac{\omega _k^{\nu _k}}{\nu _k!}.$$
(3)
$`N`$ is the number of particle types, $`\nu _k`$ is the number of particles of type $`k`$, and $`q_k`$ is the charge of one particle. The particle number $`a_k`$ indicates how many times a particle contributes to the main conserved quantity, e.g., in the case of multi-fragmentation, $`a_k`$ would be the mass number of the nuclei.

The computational challenge of Eq. (3) is hidden in the number of partitions to be summed over. Assuming $`a_k=1`$ for all $`k`$, the number of partitions is

$$J(A,N)=\frac{(A+1)!}{(A+2N)!},$$
(4)
which is prohibitive to compute for large $`A`$ or $`N`$. However, summing over all partitions can be bypassed when using recursion relations as derived in

$$Z_{A,Q}=\frac{1}{A}\underset{j}{}a_j\omega _jZ_{Aa_j,Qq_j}.$$
(5)
If intermediate values are stored, Eq. (5) requires a linear runtime in both, $`A`$ and $`N`$.

To conserve isospin as well, a sum over all possible isospin configurations for a given partition, $`\{\nu _k\}`$, and isospin weights $`W(I,M|\{\nu _k\})`$ have to be included in the partition function

$$\mathrm{\Omega }_{A,I,M,Q}=\underset{{\scriptscriptstyle \nu _ka_k}=A;{\scriptscriptstyle \nu _kq_k}=Q}{}\underset{\{\nu _k\}}{}W(I,M|\{\nu _k\})\underset{k=1}{\overset{N}{}}\frac{\omega _k^{\nu _k}}{\nu _k!}.$$
(6)
This partition function can be converted to a recursion relation by inserting $`\frac{1}{A}_{j=1}^Na_j\nu _j=1`$

$$\mathrm{\Omega }_{A,I,M,Q}=\frac{1}{A}\underset{j=1}{\overset{N}{}}a_j\omega _j\underset{{\scriptscriptstyle \nu _ka_k}=A;{\scriptscriptstyle \nu _kq_k}=Q}{}\underset{\{\nu _k\}}{}W(I,M|\{\nu _k\})\frac{\omega _j^{\nu _j1}}{(\nu _j1)!}\underset{kj}{}\frac{\omega _k^{\nu _k}}{\nu _k!}.$$
(7)
The summation indices can be switched

$$\nu _k^{}=\{\begin{array}{cc}\nu _k\hfill & ,kj\hfill \\ \nu _k1\hfill & ,k=j\hfill \end{array}$$
(8)
if the sum over isospin configurations is split:

$$\underset{\{\nu _k\}}{}W(I,M|\{\nu _k\})=\underset{I^{}=|Mm_j|}{\overset{I_j+I}{}}I_jm_j;I^{},Mm_j|IM^2\underset{\{\nu _k^{}\}}{}W(I^{},M|\{\nu _k^{}\}).$$
(9)
Finally, the recursion relation with conserved total isospin turns out to be

$$\mathrm{\Omega }_{A,I,M,Q}=\frac{1}{A}\underset{j}{}a_j\omega _j\underset{I^{}=|Mm_j|}{\overset{I+I_j}{}}I_jm_j;I^{},Mm_j|IM^2\mathrm{\Omega }_{Aa_j,I^{},Mm_j,Qq_j}.$$
(10)
### 2.2 Multiplicity Distributions

We derive a recursion relation for the multiplicity distribution of a set $`B`$ of particle types. The probability of observing $`n`$ particles of a type contained in the set $`B`$ can be expressed in terms of factorial moments

$$P_{B,A,I,M}(n)=\underset{ln}{}F_{B,A,I,M}(l)\frac{(1)^{ln}}{(ln)!n!},$$
(11)
for which recursion relations can be found more easily. Following the procedure in the previous subsection the recursion relation for the factorial moment is found to be

$$F_{B,A,I,M}(n)=\underset{kB}{}\omega _k\underset{I^{}}{}I_km_k;I^{},Mm_k|IM^2\frac{\mathrm{\Omega }_{Aa_k,I^{},Mm_k}}{\mathrm{\Omega }_{A,I,M}}F_{B,Aa_k,I^{},Mm_k}(n1),$$
(12)
which leads to the recursion relation for the multiplicity distribution

$$P_{B,A,I,M}(n)=\frac{1}{n}\underset{kB}{}\omega _k\underset{I^{}}{}I_km_k;I^{},Mm_k|IM^2\frac{\mathrm{\Omega }_{Aa_k,I^{},Mm_k}}{\mathrm{\Omega }_{A,I,M}}P_{B,Aa_k,I^{},Mm_k}(n1).$$
(13)
To test Eq. (13) we implemented a program that included the pseudo-scalar and vector mesons as well as their decays. Figure 1 shows the multiplicity distribution for neutral pions: the system was constrained to be in an isospin singlet at a temperature of $`170`$ MeV and had to produce a total of 20 pions.

For this and all following sample calculations we need to make an assumption about the energy states that are available to the particles and that have to be summed over in the single particle partition function of Eq. (2). We required that the linear momentum $`k`$ in one dimension satisfy the condition

$$kR=\pi \left(n+\frac{1}{2}\right),$$
(14)
where $`n=0,1,2,\mathrm{}`$ and $`R`$ is the size of the system. The half integers on the RHS, instead of the more usual integers, were inspired by the MIT bag model and lead to a reduction of zero-point surface energy effects. This becomes important for systems that are confined to a small volume. By assuming that our model systems are confined to a cube of volume $`V`$, we obtain the energy states

$$ϵ_{n,m,l}=\sqrt{\frac{\pi ^2}{R^2}\left[(n+1/2)^2+(m+1/2)^2+(l+1/2)^2\right]+M^2},$$
(15)
where $`n,m,l=0,1,2,\mathrm{}`$, $`M`$ is the mass of the particle, and $`R=V^{1/3}`$.

The even-odd asymmetry in Fig. 1 for small volumes that vanishes for larger volumes can be understood in the following way. At small volumes and low temperatures the system minimizes the number, not the mass, of particles produced; it prefers to produce heavy particles that decay into multiple pions rather than producing pions or kaons. These heavy particles ($`\omega `$, $`\eta `$, $`\eta ^{}`$) decay into an odd number of pions (and also into an odd number of neutral pions). Because the number of heavy particles has to be even to make 20 pions, as we imposed, the number of neutral pions is strongly favored to be even. At higher volumes or temperatures low mass particles dominate and the even-odd asymmetry vanishes.

The factorial moment from Eq. (12) corresponds to the density, when $`n=1`$. For a system of pions in an isosinglet, particle production is significantly suppressed at low volume, but converges towards the grand-canonical result at large volumes (Fig. 2).

## 3 Recursion Relations with Quantum Statistics

For a canonical ensemble of quantum particles with charge conservation the recursion relation for the partition function was given in :

$$Z_{A,Q}=\frac{1}{A}\underset{k}{}\underset{n=1}{\overset{\mathrm{}}{}}C_n^{(k)}Z_{An,Qq_k},$$
(16)
where $`C_n^{(k)}=(\pm 1)^{n1}_i\mathrm{exp}(nϵ_i^{(k)}/T)`$. To derive a partition function $`\mathrm{\Omega }_{A,I,M,Q}`$ that conserves total isospin as well, we first observe that, if the Hamiltonian of the system is independent of $`M`$, then so is the partition function. Then the partition function conserving total isospin can be related to one that conserves only the projection $`Z_{A,M,Q}=_{IM}\mathrm{\Omega }_{A,I,M,Q}`$. Therefore

$$\mathrm{\Omega }_{A,I,M,Q}=Z_{A,M,Q}Z_{A,M+1,Q}.$$
(17)
To calculate densities the 2-point function is needed

$$a_i^{}a_j=\frac{\delta _{ij}}{Z_A}\underset{n}{}(1)^{n1}\mathrm{exp}(nϵ_i/T)Z_{An},$$
(18)
which, however, can only be used when conserving additive quantities, but not isospin. The 4-point function serves to calculate higher moments

$$a_i^{}a_j^{}a_ka_l=\frac{1}{Z_A}(\delta _{il}\delta _{jk}\pm \delta _{ik}\delta _{jl})\underset{n_i,n_j}{}(1)^{n_i+n_j}\mathrm{exp}(n_iϵ_i/Tn_jϵ_j/T)Z_{An_in_j}.$$
(19)
## 4 Isospin Fluctuations

For classical particles the isospin fluctuation

$$\sigma =\underset{I=0}{}I|(n_++n_{}2n_0)^2|I,$$
(20)
where $`n_+`$, $`n_{}`$, $`n_0`$ are the densities of the respective pions, can be readily computed with the recursion relations derived above. However, for quantum particles we cannot compute the density, let alone the higher moments, when conserving total isospin. Therefore, we need to rewrite Eq. (20) in terms of isospin projection states.

The operator in Eq. (20), itself a product of spherical tensor components, can be written as a decomposition into other spherical tensor components $`A_{JM}`$

$$(n_++n_{}2n_0)^2=6T_{20}T_{20}=6\underset{J,M}{}20;20|JMA_{JM}.$$
(21)
By the Wigner-Eckart theorem only $`A_{00}`$ contributes when contracted between isosinglet states, $`\sigma =6_{I=0}I|A_{00}|I`$, which, after some algebra, becomes

$$\sigma =\frac{6}{5}\left(\underset{M=0}{}\underset{M=1}{}\right)M|\frac{3}{2}n_++\frac{3}{2}n_{}+n_0+\frac{1}{6}\left(n_++n_{}2n_0\right)^2|M.$$
(22)
In Figs. 3 and 4 we compare classical to Bose-Einstein statistics in calculations of isospin fluctuations for a system of pions at $`T=170`$ MeV. The fluctuations are divided by $`2A`$, the expected fluctuation for an independent pion gas, to highlight the difference. When the system is fixed to contain 20 pions (Fig. 3), the isospin fluctuation for classical particles does not depend on the volume, as the relative composition of the classical pion gas does not. Bose-Einstein particles, however, prefer to be all in the same state, when the volume becomes small, thus enhancing the isospin fluctuations. When the number of pions is not fixed, the isospin fluctuations converge towards the value for an independent pion gas as the system increases in size (Fig. 4).

## 5 Conclusions

Recursion relations have been derived for a canonical ensemble in which all additive charges as well as the total isospin are strictly conserved. In the case of classical statistics we derived recursion relations for the partition function and multiplicity distribution. For quantum statistics we obtained a recursion relation for the partition function and density. Work is under way to derive the recursion relation for the multiplicity distribution in the latter case.

These recursion relations can be applied to an arbitrary number of particle types and conserved additive charges. They can be computed correctly within the numerical precision of the computer in a reasonable runtime as our sample calculations show. We can now move forward and apply this technique to physical problems and in Monte Carlo algorithms that generate particles with exact conservation of additive charges and isospin. We also derived a method to quickly calculate isospin fluctuations, which makes use of the recursion relations we obtained.

## Acknowledgements

This work was supported by the National Science Foundation, Grant No. PHY-00-70818.