format: text

text:
# The DIRTY Model. I. Monte Carlo Radiative Transfer Through Dust

## 1 Introduction

Dust radiative transfer models are crucial to interpreting many astronomical observations ranging from reflection nebulae with a single star to galaxies comprised of stars, gas, and dust arranged in complex geometries. Almost all observations are affected by the presence of dust intrinsic to the system under study or present in the foreground. For dust present in the foreground as is the case for observations of single stars, the effects of dust can be accurately taken into account using a screen geometry. In this geometry, the effects are well characterized by an extinction curve which is determined by the amount and intrinsic properties (e.g., sizes and materials) of the dust (Clayton et al., 2000). However, when the dust is near or mixed with the system being studied, its effects become dependent not only on the intrinsic properties of the dust grains, but on the geometry of the stars and dust in the system as well (Witt, Thronson, & Capuano, 1992; Witt & Gordon, 1996, 2000). Since the effects of mixing stars and dust are non-trivial, a radiative transfer model is required. Using such a model produces a deeper understanding of the physical properties of the system being studied. This applies to a wide variety of astrophysical systems, from reflection nebulae to galaxies. There are a variety of numerical methods other than Monte Carlo techniques which can be used to solve the radiative transfer through dust. One example is the doubling method (Zubko & Laor, 2000). While the other methods can achieve faster computations for specific cases, Monte Carlo techniques are the most flexible and allow the study of a wider range of astrophysical systems with arbitrary geometries.

There exist many Monte Carlo dust radiative transfer models which have been constructed to study a variety of problems. In most models, the computational needs of the Monte Carlo technique have been reduced by exploiting symmetries and/or making approximations. This necessarily restricts the variety of systems to which a particular model can be successfully applied. Reflection nebulae and circumstellar shells have been successfully approximated as spherically and azimuthally symmetric systems (Witt, 1977; Witt et al., 1982; Voshchinnikov & Karjukin, 1994; Calzetti et al., 1995; Witt & Gordon, 1996). Models for bipolar nebulae, associated with young stellar objects and active galactic nuclei, usually exploit the azimuthal symmetry of such systems (Yusef-Zadeh, Morris, & White, 1984; Whitney & Hartmann, 1992; Fischer, Henning, & Yorke, 1994; Wood et al., 1998; Young, 2000). The general behavior of dust in galaxies has been studied using the concept of spherical galactic environments (Witt, Thronson, & Capuano, 1992; Witt & Gordon, 2000). Spiral galaxies have been studied with models using the approximation of disk galaxies (Bianchi, Ferrara, & Giovanardi, 1996; Wood & Jones, 1997; Kuchinski et al., 1998).

The importance of the multi-phase nature of dust distributions has motivated the construction of models which can handle arbitrary distributions of dust. Models which can handle such dust distributions include those of Wolf, Fischer, & Pfau (1998), Városi & Dwek (1999), Wood & Reynolds (1999), Bianchi et al. (2000), and our model. More recently, the importance of the re-emission of the energy absorbed by dust has been acknowledged and increasingly added to some radiative transfer models. The dust re-emission can be broken into three components; thermal equilibrium emission (large particles), non-equilibrium thermal emission (small particles), and the aromatic features (PAH-like particle emission between 3.3 and $``$20 $`\mathrm{\mu m}`$). The inclusion of the dust re-emission to radiative transfer models has been done with a range of approximations. The easiest approximation is to include dust in thermal equilibrium (Városi & Dwek, 1999; Wolf, Henning, & Stecklum, 1999; Bianchi, Davies, & Alton, 2000). Including the non-equilibrium emission and aromatic feature emission is much harder and, as a result, requires a significant amount of CPU time. Due to the complexity of the non-equilibrium and aromatic feature emission, some models use a simplified (Silva et al., 1998) and/or empirical approach to including them. The model described in this paper and a companion paper (Misselt et al., 2000a) can handle arbitrary dust distributions and fully includes all three dust emission components self-consistently.

The DIRTY model (DustI Radiative Transfer, Yeah!) grew out of spherically symmetric models which were used to study reflection nebulae (Witt et al., 1982; Gordon et al., 1994; Calzetti et al., 1995). Using Monte Carlo techniques, these models computed the radiative transfer through spherical shells of dust illuminated by point sources. The model used by Gordon et al. (1994) allowed for multiple illuminating stars arbitrarily distributed throughout a spherical nebula and the images of such a system could be constructed for arbitrary lines-of-sight. As such, the Gordon et al. (1994) model is the direct ancestor of the DIRTY model. It is important to note that the algorithm used by Gordon et al. (1994), while giving the correct total scattered light value, gave too flat a distribution of scattered light across the system. The correct algorithm for efficient computation of the appearance of a system for a particular line-of-sight is described by Yusef-Zadeh, Morris, & White (1984).

The motivation for the DIRTY model grew out of the conference “Dust Survival in Interstellar/Intergalactic Media” held at the Space Telescope Science Institute in 1994. One of the major concerns at this conference was that all Monte Carlo radiative transfer models at that time assumed a smooth distribution of dust, yet the interstellar medium is known to be clumpy over a large range of size scales (Colomb, Pöppel, & Heiles, 1980; Scalo, 1990; Rosen & Bregman, 1995). We constructed the DIRTY model to answer this concern. The main change was to move from spherical shells to rectangular cells as the basic unit of the dust density distribution. Using the simple system of a spherical nebula with a central illuminating star, Witt & Gordon (1996) explored the differences between the radiative transfer in smooth and clumpy 3-dimensional dust distributions. The implications of clumpy dust for galaxies were discussed by Witt & Gordon (2000) using the concept of spherical galactic environments. The DIRTY model has been applied to starburst galaxies (Gordon et al., 1997, 2000), spiral galaxies (Kuchinski et al., 1998), the UW Cen reflection nebula (Clayton et al., 1999), and the nucleus of M33 (Gordon et al., 1999).

The DIRTY model was constructed to compute the radiative transfer of photons from arbitrary distributions of photon emitters through arbitrary distributions of dust. This model self-consistently models the scattering and absorption of photons through dust and the re-emission of photons from dust. The algorithm for the radiative transfer, including polarization, is based on the work of Witt (1977), Yusef-Zadeh, Morris, & White (1984), and Code & Whitney (1995). The details of this algorithm are the subject of this paper. The re-emission of energy absorbed by the dust is the subject of a companion paper (Misselt et al., 2000a). The radiative transfer of the re-emitted photons is handled by the same algorithm as the stellar photons. This requires an iterative procedure to account for dust absorption and scattering of re-emitted photons and their subsequent re-emission. The details of the algorithms for the dust emission and iterative procedure are given by Misselt et al. (2000a).

## 2 Monte Carlo Radiative Transfer

The goal of the DIRTY model is to compute the radiative transfer of photons through arbitrary distributions of dust. The lack of symmetries motivated us to use Monte Carlo techniques which are based on the work of Witt (1977), Yusef-Zadeh, Morris, & White (1984), and Code & Whitney (1995). The forward scattering nature of dust grains (Gordon et al., 1994, 1997) and the observed multi-phase characteristics of the distribution of dust strongly argue for radiative transfer models based on Monte Carlo techniques. Such techniques rely on the probabilistic understanding of the interaction of a photon with a dust grain and sufficient computing power to evaluate this interaction for many paths through the distribution of the dust. This approach allows us to efficiently compute what a system of stars and dust will look like for any line-of-sight. In addition to the algorithm described in §2.1-2.3, the model requires the physical description of the dust grain properties, the dust distribution, the distribution of photon emitters, and a pseudo-random number generator.

The physical description of the dust grain properties is given by the wavelength dependent behavior of the dust optical depth, albedo, and scattering phase function. The optical depth determines the probability of a photon interaction with a dust grain, the albedo determines the probability of the dust grain scattering (or absorbing) the photon, and the scattering phase function gives the probability of the photon being scattered at a specific angle as well as the change in the polarization state of the photon as a result of the scattering. We have taken the physical description of the dust properties from work by Clayton et al. (2000) and Kim, Martin, & Hendry (1994). Clayton et al. (2000) used MEM techniques to model the size distribution of dust grains for dust extinction curves in the Milky Way, Large Magellanic Cloud, and Small Magellanic Cloud. This work provides albedos and scattering phase functions appropriate for the range of known interstellar dust extinction curves.

The DIRTY model allows for arbitrary distributions of dust limited only by the amount of memory available to store the distribution. This is accomplished by representing the dust distribution with a 3-dimensional grid. Each grid cell represents a region of uniform dust density. Thus, any dust geometry can be represented with the smallest scale of inhomogeneity being the size of a grid cell. The physical dimensions of each axis of the grid are specified by separate 1-dimensional arrays. This allows the physical dimensions of one axis to vary making it possible to efficiently model both cubical geometries as well as thin pancake-like geometries. In addition, this geometry allows for the representation of multi-phase or clumpy dust distributions. In the simplest form, this is a two-phase medium with high-density clumps and a low density inter-clump region. This type of dust distribution can be described using three parameters: the fractional filling factor of the high density dust ($`ff`$), the density ratio between the low and high density dust ($`k_2/k_1`$), and the size of a grid cell in comparison to the system size ($`1/N`$, where $`N`$ is the number of grid cells on an axis). The differences in the radiative transfer between homogeneous and two-phase dust distributions is discussed at length by Witt & Gordon (1996).

The tracking of the photons through the 3-dimensional grid is simple and quite efficient. A photon’s location and direction of travel is uniquely defined by its x, y, and z positions and the $`u`$, $`v`$, and $`w`$ direction cosines. When a photon enters a grid cell, the side through which it exits can be quickly computed in the following manner. Using the x, y, and z dimensions of the cell and the direction cosines of the photon, the distance the photon would have to travel in each of the x, y, and z directions to exit the cell through a plane parallel to the wall of the cell in those directions can be computed. The photon exits through the side which results in the shortest distance traveled. The photon’s position can be stepped by this distance and ancillary information updated (e.g., optical depth traveled).

Just as our model permits the use of arbitrary dust distributions, it also allows photons to be emitted from arbitrary source distributions. These can be a single star, a distribution of stars, or any definable surface (e.g., an accretion disk). The sources can emit isotropically or only in specific directions. We use the “ran2” pseudo-random number generator which is described by Press, Teukolsky, Vetterling, & Flannery (1992) and has a period greater than $`2\times 10^{18}`$.

### 2.1 Radiative Transfer Algorithm - A Photon’s Life

The transfer of photons through dust using Monte Carlo techniques is historically based on following a single photon through a distribution of dust. In this paradigm, the chance of the photon being absorbed or scattered is determined using a random number generator with the results dictating the life of that photon. Therefore, a particular photon might exit the system without interacting with a dust grain, be absorbed on its first scattering (or subsequent scatterings), or scatter one or more times before exiting the system. This is not a particularly efficient method of computing the direct, scattered, or absorbed light from a system as each photon can contribute to only one of these quantities. In addition, each photon exits the system in a particular direction which is not necessarily the direction from which one wants to observe the system. The efficiency of the Monte Carlo method can be increased greatly by assigning the photon a weight which allows every photon to contribute to the absorbed energy, the direct light, and the scattered light for a particular line-of-sight (Witt, 1977; Yusef-Zadeh, Morris, & White, 1984).

The algorithm we use for the radiative transfer in the DIRTY model is based on the work of Witt (1977), Yusef-Zadeh, Morris, & White (1984), and Code & Whitney (1995). In the following description of our algorithm, we will use the concept of the life of a photon. To track each photon’s state, we tabulate its position, direction of travel (via direction cosines $`u`$, $`v`$, and $`w`$), weight, and polarization (via Stokes I, Q, U, and V values).

A photon is born according to an input source distribution and its initial direction is usually chosen from an isotropic distribution (e.g., eq. 5-7 of Witt (1977)). The assumption of isotropic emission can be modified according to the particular problem being studied. For example, a mask with different sized holes at specific locations was applied in the DIRTY model to the central star of the nebula surrounding the R Corona Borealis star UW Cen (Clayton et al., 1999). In this case, only the photons emitted in the direction of the holes in the mask were allowed to move to the next step. The initial weight of the photon is usually

$$W_0^\alpha =L^\alpha ,$$
(1)
where $`L^\alpha `$ is the luminosity associated with the $`\alpha `$th photon and $`\alpha `$ runs from 1 to $`N`$. Usually $`L^\alpha =L/N`$ where $`L`$ is the luminosity of the source distribution at the modeled wavelength and $`N`$ is the number of photons in the model run. The value of $`L^\alpha `$ can be varied depending on the particular problem being studied to optimize the model running time (Yusef-Zadeh, Morris, & White, 1984). The photon is initially assumed to be unpolarized (i.e. $`S_0^\alpha =(I_0^\alpha ,Q_0^\alpha ,U_0^\alpha ,V_0^\alpha )=(1,0,0,0)`$).

The flux corresponding to the part of the photon which escapes from the system in the direction of the observer is

$$F_0^\alpha =W_0^\alpha e^{\tau (obs)_0^\alpha }\frac{1}{4\pi d^2},$$
(2)
where $`\tau (obs)_0^\alpha `$ is the optical depth along the path from the birth position of the photon to the surface of the system in the direction of the observer and $`d`$ is the distance to the system being modeled. In general, $`\tau (obs)_0^\alpha `$ is different for each photon emitted unless the system being modeled has a single central point source embedded in a spherical, homogeneous dust distribution.

The first scattering of the photon is forced to insure that every photon contributes to the scattered light (Cashwell & Everett, 1959). The optical depth to the first scattering site is

$$\tau _1^\alpha =\mathrm{ln}\left[1\xi \left(1e^{\tau _s^\alpha }\right)\right]$$
(3)
where $`\xi `$ is a random number between 0 and 1 and $`\tau _s^\alpha `$ is the optical depth to the surface of the nebula in the direction the $`\alpha `$th photon is traveling. The weight of the photon after the first scattering is then

$$W_1^\alpha =W_0^\alpha a\left(1e^{\tau _s^\alpha }\right)$$
(4)
where $`a`$ is the dust albedo. The fraction of the photon which is absorbed at the scattering site is

$$A_1^\alpha =W_0^\alpha (1a)\left(1e^{\tau _s^\alpha }\right).$$
(5)
The flux corresponding to the fraction of the photon which is scattered at the first scattering toward the observer is

$$F_1^\alpha =W_1^\alpha e^{\tau (obs)_1^\alpha }\mathrm{\Phi }(\theta (obs)_1^\alpha )\frac{1}{4\pi d^2},$$
(6)
where $`\tau (obs)_1^\alpha `$ is the optical depth from the first scattering site along the direction towards the observer, $`\mathrm{\Phi }(\theta )`$ is the scattering phase function, and $`\theta (obs)_1^\alpha `$ is the angle between the direction the photon was traveling before the scattering and the direction towards the observer. This angle is easily calculated for the $`n`$th scattering using

$$\mathrm{cos}(\theta (obs)_n^\alpha )=(u_{n1}^\alpha u_{obs}+v_{n1}^\alpha v_{obs}+w_{n1}^\alpha w_{obs})$$
(7)
where $`u_{n1}`$, $`v_{n1}`$, and $`w_{n1}`$ are the direction cosines of the photon before it scatters for the $`n`$th time and $`u_{obs}`$, $`v_{obs}`$, $`w_{obs}`$, are the direction cosines from the $`n`$th scattering site towards the observer. This equation is equivalent to eq. 17 of Yusef-Zadeh, Morris, & White (1984).

The polarization state of the flux corresponding to the fraction of the photon which reaches the observer from the $`n`$th scattering is calculated using

$$S_n^\alpha =L(\pi i_2^\alpha )R(\theta (obs)_n^\alpha )L(i_1^\alpha )S_{n1}^\alpha $$
(8)
where $`S_n^\alpha `$ is the new Stokes vector, $`L(i_1^\alpha )`$ rotates $`S_{n1}^\alpha `$ from the reference frame to the scattering plane, $`R(\theta )`$ is the scattering matrix, and $`L(\pi i_2^\alpha )`$ rotates the new Stokes vector from the scattering plane to the reference frame. The reference frame for the Stokes vector is set to be the z-axis. The rotation matrix is

$$L(\mathrm{\Psi })=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}(2\mathrm{\Psi })& \mathrm{sin}(2\mathrm{\Psi })& 0\\ 0& \mathrm{sin}(2\mathrm{\Psi })& \mathrm{cos}(2\mathrm{\Psi })& 0\\ 0& 0& 0& 1\end{array}\right]$$
(9)
where $`\mathrm{\Psi }`$ is the rotation angle. The last two equations are eqs. 2-3 of Code & Whitney (1995).

The scattering matrix, $`R(\theta )`$, is the $`4\times 4`$ Mueller matrix and is also known as the phase matrix (Bohren & Huffman, 1983; Mishchenko, Hovenier, & Travis, 2000). For this work, we have used calculations of $`R(\theta )`$ by Clayton et al. (2000) for size distributions of spherical dust grains. As a result, there are only 8 nonzero elements of $`R(\theta )`$, with $`R_{11}=R_{22}`$, $`R_{12}=R_{21}`$, $`R_{33}=R_{44}`$, and $`R_{34}=R_{43}`$. The angles $`i_1^\alpha `$ and $`i_2^\alpha `$ are calculated from the knowledge of the direction the photon is traveling and $`\theta (obs)_n^\alpha `$ using spherical geometry. A useful figure for visualizing the calculation of $`i_1^\alpha `$ and $`i_2^\alpha `$ is Fig. 1 of Code & Whitney (1995). The Stokes polarization vector is renormalized after each scattering so that $`I_n^\alpha =1`$. As a result, the flux contained in the Q polarization component is just $`Q_n^\alpha F_n^\alpha `$. Similar equations give the U and V polarization component fluxes.

The direction into which the photon is scattered at the $`n`$th scattering is described by the spherical angles $`\theta _n^\alpha `$ and $`\varphi _n^\alpha `$. The $`\theta _n^\alpha `$ angle is determined from the the dust scattering phase function. We use either a Henyey-Greenstein phase function or the full scattering matrices. The Henyey-Greenstein phase function is

$$\mathrm{\Phi }(\mathrm{cos}(\theta ),g)=\frac{1g^2}{4\pi \left[1+g^22g\mathrm{cos}(\theta )\right]^{3/2}}$$
(10)
which is eq. 3 of Witt (1977) and where $`g=\mathrm{cos}\theta `$ is the scattering phase function asymmetry. For the Henyey-Greenstein phase function, the cosine of the scattering angle is calculated using

$$\mathrm{cos}\theta _n^\alpha =\frac{1}{2g}\left[(1+g^2)\left(\frac{1g^2}{1g+2g\xi }\right)^2\right]$$
(11)
which is eq. 19 of Witt (1977). For the full scattering matrices, the scattering angle is determined from numerically inverting the function

$$\mathrm{\Phi }(\theta _n^\alpha )=I_{n1}^\alpha R_{11}(\theta _n^\alpha )+Q_{n1}^\alpha R_{22}(\theta _n^\alpha ),$$
(12)
where $`R_{11}(\theta _n^\alpha )`$ and $`R_{22}(\theta _n^\alpha )`$ are elements of the scattering matrix and $`I_{n1}^\alpha `$ and $`Q_{n1}^\alpha `$ are elements of the photon’s Stokes vector. Unlike Code & Whitney (1995) who approximate $`\mathrm{\Phi }(\theta )`$ as $`R_{11}(\theta )`$, we use the polarization dependent scattering phase function. The $`\varphi _n^\alpha `$ angle of scattering is determined using

$$\varphi _n^\alpha =\pi \left(2\xi 1\right).$$
(13)
Given the direction into which the photon is scattered ($`\theta _n^\alpha `$ and $`\varphi _n^\alpha `$), the new direction cosines can be computed from the old direction cosines (e.g., eq. 22 of Witt (1977)). The polarization state of the photon after the $`n`$th scattering is computed using eq. 8 with $`\theta _n^\alpha `$ replacing $`\theta (obs)_n^\alpha `$.

The subsequent scatterings are not forced. By forcing the first scattering we insure that every photon contributes to the scattered flux. By not forcing the subsequent scatterings we insure that the calculation is only done for the scatterings which contribute significantly to the scattered flux. The optical depth to the $`n`$th scattering, $`\tau _n^\alpha `$, is determined using

$$\tau _n^\alpha =\mathrm{ln}\xi $$
(14)
which holds for $`n2`$. If $`\tau _n^\alpha `$ is greater than the optical depth to the surface in the direction the photon is traveling, the photon escapes the system. If not, the photon is scattered in a direction determined using eqs. 11-13. The weight of the photon after the $`n`$th scattering is

$$W_n^\alpha =aW_{n1}^\alpha $$
(15)
which holds for $`n2`$. The fraction of the photon which is absorbed at the scattering site is

$$A_n^\alpha =(1a)W_{n1}^\alpha $$
(16)
which holds for $`n2`$. The fraction of the photon which is scattered towards the observer is

$$F_n^\alpha =W_n^\alpha e^{\tau (obs)_n^\alpha }\mathrm{\Phi }(\theta (obs)_n^\alpha )\frac{1}{4\pi d^2}$$
(17)
which holds for $`n1`$. The total number of scatterings a photon undergoes is limited to be less than $`n_{max}`$. For low to moderate optical depths, we set to $`n_{max}=50`$ to avoid spending significant computation time on photons whose weight is very small as $`W_n^\alpha a^n`$. It is important to remember that there are times (e.g., high $`\tau `$ values) when the value of $`n_{max}`$ should be set to a higher value. For example, if it is important to know the energy absorbed in very high optical depth parts of a dust distribution the value of $`n_{max}`$ should be increased.

### 2.2 Integrated Quantities

While the previous section has described in detail the life of a single photon, the results we are interested in are integrated quantities describing the system being modeled. These integrated quantities are images of the system (e.g., direct and scattered light images) and the 3-dimensional matrix of the energy which was absorbed by the dust.

Images of the system are constructed by projecting the birth position (for the direct light image) or the $`n`$th scattering position (for the scattered light image) of the photon onto the sky. In the DIRTY model, we embed the dust density grid in a sphere making the sky a plane which is tangent to the sphere in the direction of the observer. This plane is divided into rectangular pixels and the images are built up as each photon is run through the model. Special care must be applied when constructing the images which give the polarization state of the scattered flux, specifically the Q and U images. For the radiative transfer of the photon, the polarization state has been referenced to the z-axis of the system. As part of the construction of these images, the polarization vector of the photon must be rotated so that it is referenced to the y-axis of the image plane.

The direct light surface brightness image is constructed using

$$I_D(i,j)=\frac{1}{\mathrm{\Omega }_{pixel}}\underset{\alpha }{}F_0^\alpha $$
(18)
where $`F_0^\alpha `$ is the direct flux for the $`\alpha `$th photon (see eq. 2), $`\mathrm{\Omega }_{pixel}`$ is the surface area of pixel $`(i,j)`$ in steradians, and the sum over $`\alpha `$ is done only for photons whose birth position fall within pixel $`(i,j)`$ when projected onto the sky.

The scattered light surface brightness image is constructed using

$$I_S(i,j)=\frac{1}{\mathrm{\Omega }_{pixel}}\underset{\alpha }{}\underset{n}{}F_n^\alpha $$
(19)
where $`F_n^\alpha `$ is the scattered flux for the $`n`$th scattering of the $`\alpha `$th photon (see eqs. 6 and 17), and the sums over $`\alpha `$ and $`n`$ are done only for photons and scatterings of that photon whose scattering sites fall within pixel $`(i,j)`$ when projected onto the sky. The images which describe the polarization state of the scattered light image are constructed by first rotating the $`S_n^\alpha `$ Stokes vector from the coordinated system referenced to the z-axis to a coordinate system referenced to the y-axis of the image plane. This is done using

$$S_{}^{}{}_{n}{}^{\alpha }=L(\beta ^\alpha )S_n^\alpha $$
(20)
where $`\beta ^\alpha `$ is the angle which rotates the Stokes vector between the two coordinate systems. The Q image, and similarly the U and V images, are constructed using

$$Q(i,j)=\frac{1}{\mathrm{\Omega }_{pixel}}\underset{\alpha }{}\underset{n}{}Q_{}^{}{}_{n}{}^{\alpha }F_n^\alpha $$
(21)
where the sums over $`\alpha `$ and $`n`$ are done over the same limits as eq. 19. Of course, the sum of the direct and scattered light images gives the image one would observe at a telescope.

The 3-dimensional absorbed energy matrix is calculated using

$$A(i,j,k)=\underset{\alpha }{}\underset{n}{}A_n^\alpha $$
(22)
where the sum over $`\alpha `$ and $`n`$ are only done over scatterings which happen in cell $`(i,j,k)`$ and $`A_n^\alpha `$ is the energy which is absorbed the $`n`$th scattering site of the $`\alpha `$th photon (see eqs. 5 and 16). The 3D absorbed energy matrix is what is needed to compute the dust emission spectrum (Misselt et al., 2000a).

### 2.3 Uncertainties

Due to our use of photon weights, the uncertainties in output quantities (i.e., scattered intensity, polarization, etc.) can not be computed directly from the square root of the number of photons as is usually done when non-weighted Monte Carlo techniques are used. The ability to calculate uncertainties is crucial when using Monte Carlo techniques as the accuracy of model results are dependent on the number of photons run. If such uncertainties can be calculated during the model run, they can be used to dynamically set the number of photons needed in a model run to achieve a user input accuracy. We have adopted a simple technique for computing the upper limits on the uncertainties in the output quantities. If the output quantity is $`X`$, then

$$X=\underset{\alpha =1}{\overset{N}{}}\underset{n}{}x_n^\alpha =M\overline{x}$$
(23)
where $`x_n^\alpha `$ is the contribution of the $`n`$th scattering of the $`\alpha `$th photon to $`X`$, $`M`$ is the total number of number photons or scatterings in the model run, and $`\overline{x}`$ is the average contribution each photon or scattering makes to $`X`$. In the case of the direct flux, the sum over $`n`$ is dropped and $`M=N`$. The uncertainty in $`X`$ is then

$$\sigma _X=X\frac{\sigma _x}{\overline{x}}$$
(24)
where $`\sigma _x`$ is the standard deviation of $`\overline{x}`$. The value of $`\sigma _x`$ is calculated using

$`\sigma _x^2`$ $`=`$ $`{\displaystyle \frac{1}{M(M1)}}{\displaystyle \underset{\alpha =1}{\overset{N}{}}}{\displaystyle \underset{n}{}}(x_n^\alpha \overline{x})^2`$ (25)
$`=`$ $`{\displaystyle \frac{1}{M1}}\left(\overline{x^2}\overline{x}^2\right).`$
Again, in the case of the direct flux, the sum over $`n`$ is dropped and $`M=N`$.

The uncertainty $`\sigma _X`$ is only an upper limit on the true uncertainty since this quantity measures both the Monte Carlo noise associated with running a finite number of photons and the intrinsic variation in the quantity. For example, the total scattered flux from a nebula can be computed by directly summing every photon’s scattered weight. The uncertainty, $`\sigma _X`$, in this quantity has a contribution from the intrinsic variation of the scattered flux across the nebula as well as the Monte Carlo noise. An upper limit which would be closer to the true uncertainty can be computed by computing the point-by-point uncertainties in an image of the nebula. The point-by-point uncertainties can be calculated using

$`\sigma _X(i,j)^2`$ $`=`$ $`{\displaystyle \frac{1}{M(i,j)(M(i,j)1)}}{\displaystyle \underset{\alpha }{}}{\displaystyle \underset{n}{}}(x_n^\alpha \overline{X(i,j)})^2`$ (26)
$`=`$ $`{\displaystyle \frac{1}{M(i,j)1}}\left(\overline{X(i,j)^2}\overline{X(i,j)}^2\right).`$
where $`M(i,j)`$ is the number of scatterings contributing to output quantity in pixel $`(i,j)`$, the sum over $`\alpha `$ and $`n`$ is only done for those photons which contribute to the output quantity in pixel $`(i,j)`$, and $`X(i,j)`$ is calculated using eq. 18, 19, or 21. In the case of the direct surface brightness image, the sum over $`n`$ is dropped and $`M(i,j)`$ is the number of photons contributing to pixel $`(i,j)`$.

The uncertainty in the quantity $`X`$ can then be calculated using

$$\sigma _X^2=\underset{(i,j)}{}\sigma _X(i,j)^2.$$
(27)
Calculating the uncertainty with eq. 27 will result in a lower value of the uncertainty because the intrinsic variation of the quantity $`X(i,j)`$ across the nebula will not be included in the uncertainty calculation. This calculation of the uncertainty is still only an upper limit as the uncertainty at a specific point in a nebula will still include a contribution from photons having intrinsically different weights due to having emerged from different depths in the nebula.

We have tested this method of calculating uncertainties by running the same model multiple times, but with a different random number seed each time the model was run. This allowed us to compute the true uncertainty in an output quantity directly from the variation of the quantity between model runs. By comparing the true value of the uncertainty with the values computed using eqs. 25 and 27, we were able to evaluate the accuracy of our simple method of estimating uncertainties in output quantities. For this testing, we chose to use a sphere with a homogeneous dust distribution, a $`\tau _V=1`$, V band Milky Way dust grain properties, and a central illuminating star (see Fig. 2a). Other optical depths give similar results. We ran the model 100 times and varied the total number of photons between $`10^2`$ and $`10^5`$. Figure 1 displays the uncertainties in the scattered flux and Q component of the polarized flux as a function of the number of photons run. The general trend is for the uncertainty calculated using eq. 25 or eq. 27 to underestimate the actual uncertainty by smaller amounts as the number of photons run increases. This is due to small number statistics, especially when the uncertainty was calculated using eq. 27. For a large number of photons (e.g., $`10^5`$), the uncertainty calculated using eq. 25 or 27 is a very good estimate of the actual uncertainty. The reason is that the majority of the scattered flux and Q component of the polarized flux comes from the central region of the nebula (see Fig. 2a) and the intrinsic variation of the scattered flux in the central region is small. This implies that the uncertainty calculated from eq. 25 or 27 is dominated by Monte Carlo noise for this model. There are model systems where this will not be the case and, as a result, care must be taken calculating the uncertainty using the method outlined above.

### 2.4 Comparison with Other Models

We tested the results of the DIRTY model against those produced by Monte Carlo radiative transfer models which do not weight photons and models which use the Witt (1977) photon weighting. These models include ones which we have coded as well ones others have coded (J. Bjorkman 1999, private communication; K. Wood 1999, private communication). For computational reasons, these models are usually restricted to spherically symmetric systems with smoothly varying radial dust distributions. Our two main test cases were for $`\tau =1`$ and $`\tau =10`$. We adopted an albedo of 0.6 and a scattering phase function asymmetry of 0.6. In all cases, the Yusef-Zadeh, Morris, & White (1984) photon weighting method produced statistically similar results to the other two weighting methods. In addition, we computed the wavelength dependence of the polarization for active galactic nucleus models similar to those used by Manzini & di Serego Alighieri (1996) and found qualitative agreement with their results. Quantitative agreement is more difficult to test as we used a different dust grain model than Manzini & di Serego Alighieri (1996).

Figure 2 illustrates the images produced by the DIRTY model. Figure 2a shows how a spherical nebula with a central illuminating star would look in the V band assuming Milky Way type dust with a homogeneous distribution and a radial $`\tau _V=1`$. Figure 2b shows how a biconical nebula or active galactic nucleus inclined by an angle of 30$`\mathrm{°}`$ would look in the V band assuming Milky Way type dust with a homogeneous distribution and a $`\tau _V=1`$. Figure 3 displays the spectral energy distribution before and after the inclusion of dust in a simple starburst system. This is a good illustration of how the dust redistributes energy from the ultraviolet to the infrared. In addition, the three components (thermal equilibrium, thermal non-equilibrium, and aromatic feature emissions) of the dust emission spectrum are shown. Starburst systems are investigated in detail in Misselt et al. (2000a).

## 3 Summary and Discussion

We have presented the DIRTY radiative transfer model in this paper and a companion paper (Misselt et al., 2000a). This model uses Monte Carlo techniques to compute the radiative transfer of photons emitted from arbitrary stellar distributions through arbitrary dust distributions. The weighting scheme used in the DIRTY model results in the efficient computation of the appearance of such systems from any viewing direction. The dust re-emission includes equilibrium thermal emission, non-equilibrium thermal emission, and the aromatic feature emission. The dust emission is computed self-consistently with the dust absorption and scattering using an iterative technique. The details of the dust emission are described by Misselt et al. (2000a).

The DIRTY model can be used to study a wide range of astrophysical objects due to its flexibility. The general behavior of multi-phase dust has been studied using the DIRTY model in reflection nebula environments (Witt & Gordon, 1996) and galactic environments (Witt & Gordon, 2000). The very non-symmetric reflection nebula surrounding the R CrB star UW Cen was successfully modeled using the DIRTY model (Clayton et al., 1999). The optical depth of spiral galaxies was investigated using an early version of the DIRTY model (Kuchinski et al., 1998). The finding that the lack of a 2175 Å depression in the spectra of UV-selected starburst galaxies was not due to radiative transfer effects but due to the dust grain properties was accomplished with the DIRTY model (Gordon et al., 1997). The ultraviolet through near-IR spectral energy distribution of the M33 nucleus was modeled using the DIRTY model and discovered to be a post-starburst stellar population surrounded by a shell of Milky Way-like dust (Gordon et al., 1999). The general behavior of the infrared emission of starburst galaxies has been studied by Misselt et al. (2000a). The full ultraviolet through far-infrared spectral energy distributions of a handful of starburst galaxies was studied by Misselt, Gordon, & Clayton (2000b). The range of objects which the DIRTY model has already been applied to is a good illustration of the usefulness of such a radiative transfer model.

In addition to continuing to use the DIRTY model to explore the nature of dust and its affects on a variety of astrophysical systems, we plan to continue to improve the model itself. One such improvement would be to modify the current way the dust distribution is stored to allow for grid points to be subdivided independently allowing for a larger range of size scale to be efficiently modeled (Wolf, Henning, & Stecklum, 1999). For example, this would allow for modeling the physically small region near the accretion disk in a QSO while simultaneously modeling the large scattering regions in the jets of the QSO. Another improvement would be to include electron scattering which would improve the application of the DIRTY model to the study of QSOs. Finally, the code could be parallelized to allow larger dust distribution grids to be studied in shorter amounts of time.

As there are always more interesting astrophysical systems to be studied than time in the day, anyone interested in using the DIRTY model should contact the authors for possible collaboration.

We thank Jon Bjorkman and Kenny Wood for discussions on Monte Carlo radiative transfer, especially on the subject of the proper weighting of photons for particular lines-of-sight and how to compute the uncertainties in output quantities. We also thank the referee, Bruce Draine, for his comments which substantially improved the clarity of this paper. Support for this work was provided by NASA through LTSA Grants NAGW-3168 and NAG5-7933 and the ATP Grant NAG5-9203. KAM gratefully acknowledges financial support from the Louisiana Space Consortium through NASA grant NGT5-40035.


text:
# The Density Probability Distribution Function in Turbulent, Isothermal, Magnetized Flows in a slab geometry

## 1 Introduction

A fundamental feature of compressible turbulence is the formation of density fluctuations, a property of central interest in astrophysics, as density excesses (“clouds”) in the turbulent interstellar medium (ISM) exhibit numerous statistical properties over a range of sizes whose origin is not yet well understood. Although a full understanding of star formation requires knowledge of the full (multiple-point) statistics in order to determine mean densities as a function of region size, a first step towards this goal is the description and physical understanding of the one-point statistics, or probability density function (“PDF”) of the mass density field.

Previous studies have shown that the density PDF in turbulent compressible flows has a lognormal form in the isothermal case, but exhibits a power-law tail at densities larger (smaller) than the mean for flows with polytropic exponents smaller (larger) than unity. It should be noted that some of those works referred to purely hydrodynamic flows, while the rest referred to magneto-hydrodynamic (MHD) flows.

In particular, ref. \[\] presented a heuristic model for the development of the lognormal and power-law PDFs, based essentially on the behavior of the speed of sound with density in those types of flows. In this sense, the density PDF is a diagnostic for the dependence of the pressure with density. However, such model applied only to non-magnetic flows, while the real ISM is most likely magnetized to a significant extent. In this paper we present preliminary results on the density PDF of turbulent magnetized flows in “1+2/3” dimensions, as a first attempt to characterize the behavior of magnetic pressure with density in the fully turbulent regime. Previous works have focused on the pressure produced by weakly nonlinear Alfvén waves, but here we consider fully turbulent regimes with arbitrarily large magnetic fluctuation amplitudes.

## 2 The Non-Magnetic Case

The form of the density PDF in non-magnetic flows and its relation to the effective equation of state of the system has been understood in terms of a heuristic model by Passot and Vázquez-Semadeni, herafter PVS98 (see also ref. \[\]). As in ref. \[\], this model idealizes the generation of turbulent density fluctuations as a “multi-step” process, in which the local density at any given point in the flow is the result of a series of jumps due to the continuous passage of shock waves. The generation of densities is an iterative multiplicative process, in which every new density is obtained through a jump from the previous one, giving for the final density $`\rho _\mathrm{f}`$ in terms of the intermediate steps $`\rho _i`$: $`\rho _\mathrm{f}=\rho _0_i\left(\rho _{i+1}/\rho _i\right)`$. In the isothermal case, the density jump is given by $`M^2`$, where $`M`$ is the Mach number of the shock and, because the speed of sound is constant, a given velocity jump is always characterized by the same Mach number, regardless of the local density, so that individual jumps can be regarded as independent, but extracted from the same distribution. In terms of the variable $`s=\mathrm{ln}\rho `$, the iterative process is additive so that the Central Limit Theorem can be applied in the limit of a large number of jumps, leading to a normal distribution for $`s`$ and thus a lognormal distribution for $`\rho `$, as observed in the numerical experiments mentioned above.

Concerning the variance of the distribution, PVS98 have suggested (and confirmed numerically), through an analysis of the shock and expansion waves in the system, that for a large range of Mach numbers the typical size of the logarithmic jump is expected to be $`\sigma _sM_{\mathrm{rms}}`$, where $`M_{\mathrm{rms}}`$ is the rms Mach number. The mean $`s_0`$ of the distribution can be directly evaluated from the mass conservation condition $`\rho =_{\mathrm{}}^+\mathrm{}e^sP(s)𝑑s=1`$, where $`P(s)`$ is the PDF of $`s`$, yielding $`s_0=\sigma _s^2/2`$. The isothermal model PDF for $`s`$ thus reads

$$P(s)ds=\frac{1}{\sqrt{2\pi \sigma _s^2}}\mathrm{exp}\left[\frac{(ss_o)^2}{2\sigma _s^2}\right]ds,$$
(1)
with $`\sigma _s^2=\beta M_{\mathrm{rms}}^2`$ and $`\beta `$ a proportionality constant of order unity.

In the general polytropic case where the pressure $`P`$ behaves as $`P\rho ^\gamma `$, with $`\gamma `$ the polytropic exponent, the local Mach number $`M(s)`$ at a density $`\rho =e^s`$ is related to the one at the mean density ($`M`$) by $`M(s)=M\mathrm{exp}^{(1\gamma )s/2}`$, suggesting an ansatz where the PDF keeps the same dependence on $`M`$, provided the above replacement is made. After relocating the term in $`s_0`$ from inside the exponential function to the normalization constant, the model PDF for the polytropic case reads

$$P(s;\gamma )ds=C(\gamma )\mathrm{exp}\left[\frac{s^2e^{(\gamma 1)s}}{2M^2}\alpha (\gamma )s\right]ds.$$
(2)
This equation shows that when $`(\gamma 1)s<0`$, the PDF asymptotically approaches a power law, while in the opposite case it decays faster than a lognormal. Thus, for $`0<\gamma <1`$, the PDF approaches a power law at large densities ($`s>0`$), and at low densities ($`s<0`$) for $`\gamma >1`$.

## 3 The Magnetic Case

In what follows we restrict ourselves to the isothermal case ($`\gamma =1`$) and to a propagation along a uniform ambient magnetic field $`B_\mathrm{o}`$ and concentrate on the deviations from the corresponding lognormal PDF induced by the magnetic field.

### 3.1 Numerical Results

The numerical simulations solve the MHD equations using a pseudo-spectral method in a slab geometry (“1+2/3D”) (variability is considered only with respect to the spatial variable $`x`$ for the three components of the velocity and magnetic fields). A resolution of 2048 grid points is used allowing to handle large enough Mach and Reynolds numbers with only regular second-order viscosity. A random acceleration is applied on the $`y`$\- and $`z`$-components of the velocity field, i.e. transverally to $`B_\mathrm{o}`$, thus inducing Alfvén waves into the flow. Since these waves have finite amplitudes, they in turn induce fast and slow magnetosonic waves. The simulations are evolved over long times (several tens to a few hundred crossing times at the rms flow velocity) in order to obtain meaningful statistics for the density PDFs.

The simulations are essentially characterized by the sonic and Alfvénic Mach numbers, respectively defined as $`M_\mathrm{s}u/c_\mathrm{s}`$ and $`M_\mathrm{a}u/v_\mathrm{A}`$, where $`u`$ is the rms velocity, $`c_\mathrm{s}`$ the sound speed and $`v_\mathrm{A}B_\mathrm{o}^2/\rho `$ the Alfvén speed. In order to investigate the effect of the magnetic field exclusively, we consider two runs (denoted I and II) with approximatively the same value of $`M_\mathrm{s}`$ (4.06 and 3.74 respectively) but very different values of $`M_\mathrm{a}`$ (0.36 and 1.66 resp.). These quantities are integrated over the duration of the run, excluding the first few time units in order to avoid including the uniform-density initial conditions. The rms relative field fluctuations $`\delta B/B`$ are for runs I and II respectively 0.75 and 6.65.

In fig. 2 we show the PDFs for the two magnetic runs, together with lognormal fits (dashed lines). It is clearly seen that, contrary to what one might expect, the high-$`M_\mathrm{a}`$ run, which should be closest to the pure-hydro case, exhibits a strong deviation from lognormality, while the low-$`M_\mathrm{a}`$ PDF is quite close to a lognormal. This implies that the case with the weakest mean field is the one with the strongest effect of magnetic pressure on the density PDF.

### 3.2 Interpretation in Terms of Simple MHD waves

A preliminary interpretation of the behavior reported in the previous section can be given in terms of the so-called “simple MHD waves” (see, e.g., ref. \[\]). These are the finite-amplitude equivalent of linear MHD waves, and have the same three well-known modes: Alvén, slow and fast. A simple wave refers to a solution that depends only on a single variable, any combinations of the spatial ($`x`$) and temporal ($`t`$) independent variables. However, as mentioned above, simple waves have finite amplitudes, and in general they develop shocks in a finite time.

A number of properties of simple waves has been reported by Mann (see also Jeffrey and Taniuti). Among them, most relevant for the present discussion are the propagation velocities of the three modes, and the relations between the density and magnetic field fluctuations. The Alfvén mode with speed $`V`$ such that $`V^2=v_{\mathrm{A}}^{}{}_{x}{}^{2}`$ is not associated with density fluctuations. The fast and slow speed are given by

$$V_\pm ^2=\frac{v_{\mathrm{A}}^{}{}_{}{}^{2}+c_{\mathrm{s}}^{}{}_{}{}^{2}}{2}\left\{1\pm \left[1\frac{4c_{\mathrm{s}}^{}{}_{}{}^{2}v_{\mathrm{A}}^{}{}_{x}{}^{2}}{(v_{\mathrm{A}}^{}{}_{}{}^{2}+c_{\mathrm{s}}^{}{}_{}{}^{2})^2}\right]^{1/2}\right\}$$
(3)
and the relation between the magnitude of the magnetic field $`b`$ and the density is given by

$$\frac{db}{d\rho }=\frac{V^2c_{\mathrm{s}}^{}{}_{}{}^{2}}{b},$$
(4)
where $`v_{\mathrm{A}}^{}{}_{x}{}^{2}=b_x^2/\rho `$, $`v_\mathrm{A}^2=b^2/\rho `$, and $`c_\mathrm{s}`$ is the sound speed. An important point to note is that eqs. (3) and (4) imply a positive correlation between $`b`$ and $`\rho `$ for the fast mode ($`V_+`$) and an anticorrelation for the slow mode ($`V_{}`$),.

Inserting eq. (3) in eq. (4), one can find the density dependence of the magnetic pressure. Excluding the case where the $`\beta `$ of the plasma is of order unity with at the same time small to moderate field distorsions, this dependence simplifies and reads (general numerical solutions have been presented by Mann):

$`{\displaystyle \frac{b^2}{b_x^2}}`$ $``$ $`a_1a_2{\displaystyle \frac{c_\mathrm{s}^2\rho }{b_x^2}}(\mathrm{slow}\mathrm{mode})`$ (5)
$`b^2`$ $``$ $`2a_3\rho ^2(\mathrm{fast}\mathrm{mode}),`$ (6)
where $`b^2=b_x^2+b_{}^2`$, with $`b_{}^2=b_y^2+b_z^2`$ being the magnitude of the magnetic field fluctuation. Note that in our 1+2/3 geometry, $`b_x(=B_\mathrm{o})`$ is constant. The quantities $`a_i`$ denote integration constants.

From these relations, it can be seen that the magnetic pressure, $`b^2`$, scales roughly linearly with the density (albeit inversely) in the case of the slow mode, but quadratically in the case of the fast mode. Note also that, in the case of small field fluctuations (most often the case if $`\beta c_\mathrm{s}^2/(b_x^2/\rho )1`$) the coefficient of $`\rho `$ in eq. (5) is small, implying that, for the slow mode, large density fluctuations can occur even for small variations of $`b`$. Thus, the slow mode is expected to dominate density fluctuation production in the case of small field fluctuations. Instead, at large enough field fluctuations, the quadratic dependence of $`b^2`$ on $`\rho `$ of the fast mode eventually overwhelms the linear dependence of the slow mode (and also of the thermal pressure), so the fast mode is expected to dominate the density fluctuation production.

The dependence of the magnetic pressure on the density seems to describe the observations of sec. 3.1 adequately because, according to the results of the non-magnetic case, a pressure that depends linearly on the density produces a lognormal density PDF, while one with a higher-than-linear dependence produces a near power law at low densities.

## 4 Conclusions

### 4.1 Summary

In this paper we have reviewed our previous results on the development of the density PDF as a consequence of the effective (polytropic) equation of state, and applied them to an understanding of the PDF in the magnetic case.

In the context of the model presented in ref. \[\], isothermal non-magnetic flows have lognormal mass density PDFs while in the non-isothermal cases, a power-law tail develops at high Mach number, at either $`\rho >\rho `$ or $`\rho <\rho `$ depending on whether $`\gamma <1`$ or $`\gamma >1`$, respectively.

In the magnetic case, we considered only isothermal cases and a propagation parallel to the ambient magnetic field, but reported that a deviation from the lognormal PDF occurs when the magnetic fluctuations are large, a case expected when the mean field is small. We interpreted this effect in terms of so-called “simple”, finite-amplitude nonlinear MHD waves. Indeed, for the slow mode of these waves, the magnetic pressure $`b^2`$ behaves linearly with the density $`\rho `$, while for the fast mode, the magnetic pressure behaves quadratically with $`\rho `$. Together with the fact that for the slow wave the density is weighted by a small factor, this suggests that the production of density fluctuations is dominated by the slow waves at small field deviations (i.e., $`B_\mathrm{o}`$ large) giving a lognormal PDF again, while for weak mean fields (large field deviations) the quadratic density dependence of $`b^2`$ produces an excess at small densities in the PDF, corresponding to the $`\gamma >1`$ case of the polytopic description. A more detailed study, including in particular the case of non-parallel propagation is in progress.

### 4.2 Astrophysical Implications

The results of this paper have a number of implications in the astrophysical context. First, we have seen that the main parameter determining the effect of the magnetic field appears to be the magnetic fluctuation amplitude, $`\delta B/B`$, rather than $`\beta `$ alone, which is the parameter most frequently used to characterize MHD flows. In other words, it is important to have a knowledge of the relative importance of the turbulent fluctuations (as given by the Alfvénic Mach number) in addition to simply the ratio of thermal to magnetic pressures (as given by $`\beta `$).

Second, as far as the density PDF is concerned, the effect of the magnetic field is most notorious when the field fluctuations are large (generally, when the mean field is weak). This suggests that the limit of vanishing field does not approach the non-magnetic case, and in this sense the latter is singular. This may imply that it is inadequate to model the ISM as a non-magnetic flow even if the degree of magnetization is very low, as has started to be recently contended at the level of molecular clouds (see, e.g., ref. \[\]).

## Acknowledgements

Part of this work was completed while the authors enjoyed the hospitality of the Institute for Theoretical Physics at the University of California at Santa Barbara, as participants of the Astrophysical Turbulence Program. This work has received partial funding from grants Conacyt (México) 27752-E to E.V.-S, the CNRS National Program “PCMI” to T.P. and NSF PHY94-07194.


text:
# Applying Machine Translation to Two-Stage Cross-Language Information Retrieval

## 1 Introduction

The number of machine readable texts accessible via CD-ROMs and the World Wide Web has been rapidly growing. However, since the content of each text is usually provided in a limited number of languages, the notion of information retrieval (IR) has been expanded so that users can retrieve textual information (i.e., documents) across languages. One application, commonly termed “cross-language information retrieval (CLIR)”, is the retrieval task where the user presents queries in one language to retrieve documents in another language. Thus, as can be predicted, CLIR needs to standardize queries and documents into a common representation, so that monolingual IR techniques can be applied. From this point of view, existing CLIR can be classified into three approaches.

The first approach translates queries into the document language , while the second approach translates documents into the query language . The third approach projects both queries and documents into a language-independent representation by way of thesaurus classes  and latent semantic indexing .

Although extensive comparative experiments among different approaches in a rigorous manner are difficult and expensive, a few cases can be found in past CLIR literature.

Oard  compared the query and document translation methods. For the purpose of English-German CLIR experiments, he used the 21 English queries and SDA/NZZ German collection consisting of 251,840 newswire articles, contained in the TREC-6 CLIR collection. Then, he showed that the MT-based query translation with the Logos system was more effective than various types of dictionary-based query translation methods, and that the MT-based document translation method further outperformed the MT-based query translation method. Those findings were salient especially when the length of queries was large.

McCarley  conducted English/French bidirectional CLIR experiments, where the 141,656 AP English documents and 212,918 SDA French documents in the TREC-6 and TREC-7 collections were used, and applied a statistical MT method to both query and document translation methods. He showed that the relative superiority between query and document translation methods varied depending on the source and target language pair. To put it more precisely, in his case, the quality of French-English translation was better than that of English-French translation, for both query and document translations.

In addition, he showed that a hybrid method, where the relevance degree of each document (i.e., the “score”) is the mean of those obtained with query and document translation methods, outperformed methods based on either query or document translation, irrespective of the source and target language pair. Possible rationales include that since machine translation is not an invertible operation, query and document translations mutually enhance the possibility that query terms correspond to appropriate translations in documents.

To sum up, the MT-based document translation approach is potentially effective in terms of retrieval accuracy. Besides this, since retrieved documents are mostly in a user’s non-native language, the document translation approach is significantly effective for browsing and interactive retrieval.

However, a major drawback of this approach is that the full translation on large-scale collections is prohibitive in terms of computational cost. In fact, Oard , for example, spent approximately ten machine-months in translating the SDA/NZZ collection. This problem is especially crucial in the case where the number of user languages is large, and documents are frequently updated as in the Web. Although a fast MT method  was proposed, this method is currently limited to MT within European languages, which are relatively similar to one another.

In view of the above discussions, we propose a method to minimize the computational cost required for the MT-based document translation, which is fundamentally twofold. First, we translate the query into the document language, and retrieve a fixed number of top-ranked documents (one thousand, for example). Second, we machine translate those documents into the query language, and then re-rank those documents based on the score, combining those individually obtained with query and document translation methods. Consequently, it is expected that the retrieval accuracy is improved with a minimal MT cost.

From a different perspective, our method can be classified as a two-stage retrieval principle. However, in the monolingual two-stage IR, the second stage usually involves re-calculation of term weights and local feedback so as to increase the number of relevant documents in the final result , and that in the case of existing two-stage CLIR, multiple stages are used to improve the quality of query translation .

Section 2 describes our two-stage CLIR system, where we elaborate mainly on the MT-based re-ranking method. Section 3 then evaluates the performance of our system, using the NACSIS test collection , which consists of 39 Japanese queries and approximately 330,000 technical abstracts in English and Japanese.

## 2 System Description

### 2.1 Overview

Figure 1 depicts the overall design of our Japanese/English bidirectional CLIR system, in which we combined query and document translation modules with a monolingual retrieval system. In this section, we explain the retrieval process based on this figure.

First, given a query in the source language (S), a query translation is performed to output a translation in the target language (T). In this phase, we use two alternative methods. The first method is the use of an MT system, for which we use the Transer Japanese/English MT system.111Developed by NOVA, Inc. This MT system uses a general bilingual dictionary consisting of 230,000 entries, and 19 optional technical dictionaries, among which a computer terminology dictionary consisting of 100,000 entries is combined with our system.

However, since in most cases, queries consist of a small number of keywords and phrases, word/phrased-based translation methods are expected to be comparable with MT systems, in terms of query translation. Thus, for the second method, we use the Japanese/English phrase-based translation method proposed by Fujii and Ishikawa , which uses general/technical dictionaries to derive possible word/phrase translations, and resolves translation ambiguity based on statistical information obtained from the target document collection. In addition, for words unlisted in dictionaries, transliteration is performed to identify phonetic equivalents in the target language.

Second, the monolingual retrieval system searches a collection for documents relevant to the translated query, and sorts them according to the degree of relevance (i.e., the score), in descending order. For English documents, we use the SMART system , where the augmented TF$``$IDF term weighting method (“atc”) is used for both queries and documents, and the score is computed based on the similarity between the query and each document in a term vector space. For Japanese documents, we implemented a retrieval system based on the vector space model.

Consequently, only the top $`N`$ documents are selected as an intermediate retrieval result, where $`N`$ is a parametric constant.

Third, the top $`N`$ documents are translated into the source language. Note that unlike the query translation phase, we use solely the Transer MT system, because translations are aimed primarily at human users, and thus the phrase-based translation method potentially degrades readability of retrieval results.

Finally, the $`N`$ documents translated are re-ranked according to the new score. To accomplish this task, we compute the similarity score between the source query (submitted by the user) and each translated document in the term vector space, as performed in the first retrieval stage. We then compute the new score by averaging those obtained independently with English and Japanese monolingual similarity computations. We will elaborate on this process in Section 2.2.

Note that by decreasing the value of $`N`$, we can decrease the computational cost required for machine translation. However, this also decreases the number of relevant documents contained in the top $`N`$ set, and potentially dilutes the effectiveness of the re-ranking. For example, in an extreme case where the top $`N`$ set contains no relevant document, the re-ranking procedure does not change the retrieval accuracy.

The re-ranking procedure is similar to McCarley’s hybrid method , in the sense that his method also combines scores obtained with query and document translations. However, unlike McCarley’s method, which needs to translate the entire document collection prior to the retrieval, in our method the overhead for translating documents is minimized and can be distributed to each user. In other words, the second stage can be performed on each client (i.e., users’ computers or Web browsers). In fact, there are a number of commercial Web browsers combined with MT systems, and thus it is feasible to additionally introduce the re-ranking function to those browsers. Besides this, we can easily replace the MT system with a newer version or those for other language pairs.

### 2.2 MT-based Re-ranking Method

First, given the top $`N`$ documents retrieved and translated into the source language, we first compute the similarity score between each document and the source query provided by the user. Following the vector space model, both queries and documents are represented by a vector consisting of statistical factors associated with indexed terms (i.e., term weights).

In conventional retrieval systems, documents are indexed to produce an inverted file, prior to the retrieval, so that documents containing query terms can efficiently be retrieved even from a large-scale collection. However, in the case of our re-ranking process, since (a) the number of target documents is limited, and (b) real-time indexing degrades the time efficiency, we prefer to use a simple pattern matching method, instead of the inverted file.

For term weighting, we tentatively use a variation of TF$``$IDF , as shown in Equation (1).

$$\begin{array}{ccc}TF\hfill & =\hfill & 1+\mathrm{log}(f_{t,d})\hfill \\ \multicolumn{3}{c}{}\\ IDF\hfill & =\hfill & \mathrm{log}\frac{N}{n_t}\hfill \end{array}$$
(1)
Here, $`f_{t,d}`$ denotes the frequency that term $`t`$ appears in document $`d`$. Note that unlike the common IDF formula, $`N`$ denotes the number of documents retrieved in the first stage (see Section 2.1), and $`n_t`$ denotes the number of documents containing term $`t`$, out of $`N`$ documents.

One may argue that since in our case where the number of target documents is considerably smaller than that of the entire collection, a different term weighting method is needed. For example, the IDF formula proposed for large-scale document collections may be less effective for a limited number of documents. However, a preliminary experiment showed that the use of IDF marginally improved the performance obtained without IDF. On the other hand, since the preliminary experiment showed that the use of document length considerably degraded the performance, we compute the similarity between the query and each document, as the inner product (instead of the cosine of the angle) between their associated vectors.

Thereafter, for each document, we combine two similarity scores obtained in English-English and Japanese-Japanese retrieval processes. We shall call them $`ESIM`$ and $`JSIM`$, respectively. Since those two similarity scores have different ranges, we use a geometric mean, instead of an arithmetic mean, as shown in Equation (2).

$`SIM`$ $`=`$ $`ESIM^\alpha JSIM^\beta `$ (2)
Here, $`SIM`$ is the final similarity score with which we re-rank the top $`N`$ documents, and $`\alpha `$ and $`\beta `$ are parametric constants used to control the degree to which $`ESIM`$ and $`JSIM`$ affect the computation of $`SIM`$. However, in the case where either $`ESIM`$ or $`JSIM`$ is zero, the value of $`SIM`$ always becomes zero, disregarding the value of the other similarity score. To avoid this problem, in such a case we arbitrarily assign the value 0.0001 to either $`ESIM`$ or $`JSIM`$ that takes zero.

Possible factors to set values of $`\alpha `$ and $`\beta `$ include the quality of Japanese-English and English-Japanese translations. In the case where the quality of one of the translations is considerably lower, $`\alpha `$ and $`\beta `$ must be properly set so as to decrease the effect of the similarity score through the lower quality translation. Generally speaking, the quality of English-Japanese translation is higher than that of Japanese-English translation, because morphological and syntactic analyses for Japanese are usually more crucial than those for English. However, we empirically set $`\alpha =\beta =1`$, that is, we consider $`ESIM`$ and $`JSIM`$ equally in the re-ranking process.

## 3 Experimentation

### 3.1 Methodology

We investigated the performance of several versions of our system in terms of Japanese-English CLIR, where each system outputs the top 1,000 documents, and the TREC evaluation software was used to calculate non-interpolated average precision values.

For the purpose of our experiments, we used the official version of the NACSIS test collection . This collection consists of 39 Japanese queries and approximately 330,000 documents (in either a combination of English and Japanese or either of the languages individually), collected from technical papers published by 65 Japanese associations for various fields.

Each document consists of the document ID, title, name(s) of author(s), name/date of conference, hosting organization, abstract and keywords, from which titles, abstracts and keywords were indexed by the SMART system. We used as target documents 187,081 entries that are in both English and Japanese.

Each query consists of the query ID, title of the topic, description, narrative and list of synonyms, from which we used only the description. Figure 2 shows example descriptions (translated into English by one of the authors).

The NACSIS collection was produced for a TREC-type (CL)IR workshop held by NACSIS (National Center for Science Information Systems, Japan) in 1999.222See http://www.rd.nacsis.ac.jp/~ntcadm/workshop/work-en.html for details of the NACSIS workshop. In this workshop, each participant was allowed to submit more than one retrieval result using different methods. However, at least one result had to be gained with only the description field in queries. According to experimental results reported in the proceedings of the workshop , in the case where only the description field was used, average precision values ranged from 0.021 to 0.182.

Relevance assessment was performed based on the pooling method . To put it more precisely, candidates for relevant documents were first pooled by multiple retrieval systems (primarily systems that participated in the NACSIS workshop). Thereafter, for each candidate document, human expert(s) assigned one of three ranks of relevance, that is, “relevant”, “partially relevant” and “irrelevant”. The average number of candidate documents pooled for each query is 2,509, among which the number of relevant and partially relevant documents are approximately 21 and 6, respectively. In our experiments, we did not regard “partially relevant” documents as relevant ones, because interpretation of “partially relevant” is not fully clear to the authors. Note that since the NACSIS collection does not contain English queries, we cannot estimate a baseline for Japanese-English CLIR performance using English-English IR.

In the following two sections, we will show experimental results in terms of the first and second stages (i.e., query translation methods and the MT-based re-ranking method), respectively.

### 3.2 Evaluation of Query Translation Methods

The primal objective in this section is to compare the effectiveness of the phrase-based translation method proposed by Fujii and Ishikawa  and one based on the Transer MT system, in terms of Japanese-English query translation. While the former method is aimed solely at words and phrases, the MT system can also be used for full sentences. In addition, since both methods are, to some extent, complementary to each other, we theoretically gain a query expansion effect, combining query terms translated by individual methods. In view of those above factors, we compared the following query translation methods:

* the use of the Transer MT system for full sentences contained in the description field (“MTS”),
* the use of the Transer MT system for content words and phrases extracted from the description field, for which the ChaSen morphological analyzer  was used (“MTP”),
* the phrase-based translation method applied to the same words and phrases as used for the MTP method (“PBT”),
* the use of query terms obtained with both MTP and PBT, where terms outputed by both methods are considered to appear twice in the query (“MPBT”).

Table 1 shows the non-interpolated average precision values, averaged over the 39 queries, for different query translation methods listed above. The second column denotes the average number of query terms provided with each translation method, some of which were potentially discarded as stopwords by the SMART system. The third column denotes average precision values for different query translation methods. We will explain the fourth and fifth columns in Section 3.3.

Looking at this table, one can see that while two MT-based methods, that is, MTS and MTP, were quite comparable in performance, and that PBT outperformed both of them. In the case of PBT, the transliteration successfully identified English equivalents for katakana words unlisted in the word dictionary, such as “coraboreishon (collaboration)” and “mobairu (mobile)”, which the MT-based methods failed to translate. Another reason was due to the difference in dictionaries used. Generally speaking, PBT tended to output technical words more than the MT-based methods. For example, for Japanese phrases “fukusuu-deeta” and “sekitsui-doubutsu”, PBT outputed “multiple data” and “craniate”, while MTS/MTP outputed “more than one data” and “vertebrate”, respectively. Note that this effect was evident partially because the NACSIS collection consists of technical documents. In addition, MPBT further improved the performance of PBT. Although the difference between PBT and MPBT was marginal, it is worth utilizing both the MT-based and phrase-based methods, if available, for query translation.

To validate those above results in a thorough manner, we used the non-parametric Wilcoxon matched-pairs signed-test for statistical testing (at the 5% level), which investigates whether the difference in average precision is meaningful or simply due to chance . We found that differences in average precision values for pairs “MTP versus MTS”, “MPBT versus MTS”, and “MPBT versus MTP” were significant, although for other pairs, we could not obtain sufficient evidence to conclude a statistical significance. To sum up, we concluded that in query translation, a combination of MT-based and phrase-based translation methods was more effective than a method relying solely on the MT system.

### 3.3 Evaluation of the MT-based Re-ranking Method

First, we consider Table 1 again, where the fourth column “MT” denotes the average precision values for each query translation method, combined with the MT-based re-ranking method. Throughout our experimentation in this paper, the best average precision value by an automatic method was 0.2013 (i.e., one obtained by PBT combined with the MT-based re-ranking method), which is relatively high, when compared with average precision values reported in the NACSIS workshop (ranging from 0.021 to 0.182).

For each query translation method, the improvement in average precision from one without the re-ranking, which is generally noticeable, is indicated in parentheses. In fact, we used the Wilcoxon test again, as conducted in Section 3.2, and confirmed that every improvement was statistically significant. To sum up, the MT-based re-ranking method we proposed was generally effective, irrespective of the query translation method combined, in terms of CLIR performance.

Second, we conducted an error analysis for queries for which the re-ranking method degraded the average precision, and found that roughly two thirds of errors were due to ambiguity in the document translation. For example, the English word “library” was often incorrectly translated into “raiburari (library as a software)”, whereas the original query was intended to “toshokan (library as an institution)”.

Third, to estimate the upper bound of the re-ranking method, as denoted in the fifth column “HT”, we used as human translations Japanese documents comparable to English ones in the NACSIS collection. By comparing the results of “MT” and “HT”, one can see that MT systems with a higher quality, if available, are expected to further improve our CLIR system. In fact, when we manually corrected inappropriate translations in translated documents, such as “library (raiburari/toshokan)” above, the average precision of “MT” became almost equivalent to that of “HT”.

Noted that when combined with the re-ranking method, differences among query translation methods in average precision were relatively overshadowed. In the case of “MT”, the Wilcoxon test showed that differences in only pairs “MPBT versus MTS” and “MPBT versus MTP” were significant, while in the case of “HT”, none of the differences were identified as significant.

Fourth, we investigated how the number of documents retrieved in the first stage (i.e., the value of $`N`$ in Section 2) affected the performance of the re-ranking method. As discussed in Section 2.1, in real world usage, one has to consider the trade-off between the retrieval accuracy (i.e., average precision in our case) and overhead required for the document translation.

Table 2 shows the results, where average precision values in the column “1,000” correspond to those in Table 1. By comparing average precision values for each of four query translation methods (i.e., MTS, MTP, PBT and MPBT) and those suffixed with “+MT” and “+HT” in Table 2, one can see that the re-ranking methods were effective, irrespective of the number of documents retrieved. In other words, it is expected that we can minimize the overhead in translating documents, without decreasing the retrieval accuracy.

Table 3 shows CPU time (sec.) required for the document translation and re-ranking procedures, averaged over four different query translation methods. In the case of $`N=1,000`$, the total CPU time was approximately three minutes, which is perhaps not tolerable for a real-time usage. However, for small values of $`N`$ (e.g., 50 and 100), the CPU time was more acceptable and practical, maintaining the improvement of retrieval accuracy.

## 4 Conclusion

Reflecting the rapid growth in utilization of machine readable texts, cross-language information retrieval (CLIR) has variously been explored in order to facilitate retrieving information across languages.

In brief, existing CLIR systems are classified into three approaches: (a) translating queries into the document language, (b) translating documents into the query language, and (c) representing both queries and documents in a language-independent space. Among these approaches, the second approach, based on machine translation, is effective in terms of retrieval accuracy and user interaction. However, the computational cost in translating large-scale document collections is prohibitive.

To resolve this problem, we proposed a two-stage CLIR method, in which we first used a query translation method to retrieve a fixed number of documents, and then applied machine translation only to those documents, instead of the entire collection, to improve the document ranking.

Through Japanese-English CLIR experiments using the NACSIS collection, we showed that our two-stage method significantly improved average precision values obtained solely with query translation methods. We also showed that our method performed reasonably, even in the case where the number of retrieved documents was relatively small.

## Acknowledgments

The authors would like to thank NOVA, Inc. for their support with the Transer MT system, and Noriko Kando (National Institute of Informatics, Japan) for her support with the NACSIS collection.


text:
# The excitonic collapse of higher Landau level fractional quantum Hall effect

## Abstract

The scarcity of the fractional quantum Hall effect in higher Landau levels is a most intriguing fact when contrasted with its great abundance in the lowest Landau level. This paper shows that a suppression of the hard core repulsion in going from the lowest Landau level to higher Landau levels leads to a collapse of the energy of the neutral excitation, destabilizing all fractional states in the third and higher Landau levels, and almost all in the second Landau level. The remaining fractions are in agreement with those observed experimentally.

Electrons in the lowest Landau level exhibit the spectacular phenomenon of the fractional quantum Hall effect (FQHE). Including the several new fractions observed in the ultra-low temperature measurements in very high quality samples, there now exists evidence for more than 50 fractions in the lowest Landau level. The essential phenomenon is remarkably insensitive to the detailed form of the repulsive interaction, as manifested by the fact that it is quite robust to perturbations arising from Landau level mixing, finite thickness, or the nature of the transverse confinement.

The system of electrons restricted to a higher Landau level (LL) deviates from that in the lowest LL only through the short-distance matrix elements for the Coulomb interaction. Consequently, one would expect the FQHE to be not too sensitive to the Landau level index either. From this point of view, it is astounding that the fractional quantum Hall effect (FQHE) is so rare in higher Landau levels. For example, whereas 10 members of each of the sequences $`\nu =n/(2n+1)`$ and $`\nu =n/(2n1)`$ have been observed in the lowest LL, only the first one has been observed in the second LL, and none at all in the third and higher LLs. In fact, the only fractions outside the lowest LL for which clear experimental evidence exists, in the form of reasonably well quantized plateaus, are $`\nu =1/3`$ and $`\nu =1/2`$ in the second LL. (That the latter has no analog in the lowest LL further underscores the striking difference between the lowest and the higher LL physics.) Hartree-Fock variational (Koulakov, Fogler, and Shklovskii), exact diagonalization (Rezayi, Haldane, and Yang ), and experimental (Lilly et al., Du et al. ) studies make a compelling case that a bubble crystal or an anisotropic stripe phase is favored over the FQHE in higher LLs.

The goal of this work is to start from the FQHE end of the problem and ask by what mechanism is the FQHE destroyed in higher Landau levels. Of course, there can be a mundane origin for the disappearance of the FQHE, e.g., sufficiently strong disorder, but of interest to us here is the possible intrinsic instability of FQHE in higher Landau levels. To this end, we will start by assuming an incompressible FQHE state and then investigate its stability to quantum fluctuations. A necessary condition for FQHE is that (in the absence of disorder) the energy gap to all excitations remain positive definite. A vanishing of the gap signals an instability of the assumed ground state. It will be shown that the suppression of the short-distance Coulomb matrix elements in higher Landau levels leads to a collapse of the energy of the neutral exciton, destabilizing most of the incompressible FQHE states in higher LLs.

The investigations will be carried out within the framework of the composite fermion (CF) theory of the FQHE, which has proven successful in capturing subtle instabilities of the FQHE. At small filling factors ($`\nu 1/9`$), the neutral composite-fermion exciton becomes gapless. In another example, a recent theoretical study  has revealed that the Fermi sea of composite fermions is unstable to Cooper pairing at $`\nu =5/2`$ but not at $`\nu =1/2`$, giving insight into the experimental observation of the FQHE at $`\nu =5/2`$ but not at $`\nu =1/2`$.

A composite fermion, 2pCF, is the bound state of an electron and $`2p`$ quantum mechanical vortices. Of relevance to experiment are three flavors of composite fermions carrying two, four, and six vortices, namely 2CFs, 4CFs, and 6CFs. According to the composite fermion theory, the interacting electrons at the LL filling factor $`\nu _0\frac{n}{2pn+1}`$ transform into weakly interacting composite fermions with an effective filling $`\nu _0^{}=n`$. The ground state here corresponds to $`n`$ filled CF-LLs and the neutral excitation to a particle-hole pair of composite fermions, called the CF exciton. The explicit, parameter-free, lowest-LL form for the microscopic wave functions for the fully polarized CF ground state and the CF exciton can be found in the literature and will not be repeated here. In the lowest Landau level, these wave functions have been found to be accurate in tests against exact diagonalization results available for small systems. They are not as good in higher Landau levels for 2CFs, but should still be a reasonable starting point provided the ground state is an incompressible state. Interestingly, the wave functions are quantitatively accurate for 4CFs in the second LL.

We will study filling factors $`\nu =2s+\nu _0`$. Here, the quantity $`s=0,1,\mathrm{}`$ denotes the Landau level index, with the lowest $`s`$ LLs having both spin states filled, thus contributing $`2s`$ to the filling factor. As a simplification we will take these $`2s`$ LLs to be completely inert (neglecting by the lower Landau levels) and work only with the electrons in the topmost partially filled LL, which will be taken to be fully spin polarized; this is a valid approximation in the limit of high magnetic fields when LL mixing is negligible. Note that only the electrons in the topmost partially filled LL capture vortices to turn into composite fermions; the electrons in the lower, filled LLs remain unchanged.

The wave functions $`\mathrm{\Psi }`$ are most easily constructed within the lowest LL. To calculate energies in higher LLs one may promote them to higher LLs by application of the Landau level raising operator. However, this procedure is technically rather cumbersome. We instead proceed by working with an effective interaction in the lowest LL which mimics the Coulomb interaction in a higher LL. The interparticle interaction in any given Landau level is completely specified by its Haldane pseudopotentials, $`V_m^{(s)}`$, defined via $`V(r_ir_j)=_mV_m^{(s)}P_m^{(s)ij},`$ where $`V_m^{(s)}`$ is the interaction energy of two particles in the $`s`$th LL in the relative angular momentum $`m`$ state, and $`P_m^{(s)ij}`$ is the corresponding projection operator. Following Park et al., we map the problem of a given interaction in the $`s=1`$ Landau level into that of an effective interaction $`V^{}(r)`$ in the lowest ($`s=0`$) Landau level, chosen so that the two have the same pseudopotentials, i.e., $`V_m^{(1)}=V_m^{(0)}`$. We implement this strategy in an approximate scheme by taking the following convenient form for the effective potential:

$$V^{}(r)=\frac{e^2}{ϵ}\left(\frac{1}{r}+a_1^{}e^{\alpha _1^{}r^2}+a_2^{}r^2e^{\alpha _2^{}r^2}\right).$$

The parameters $`a_1^{}`$, $`a_2^{}`$, $`\alpha _1^{}`$, and $`\alpha _2^{}`$ are fixed by requiring that the first three to four odd pseudopotentials of the effective potential in the lowest LL match exactly the corresponding pseudopotentials of the Coulomb potential in the second Landau level. Only the odd pseudopotentials are relevant due to the antisymmetric form of the spatial part of the wave function. The remaining higher order pseudopotentials are asymptotically correct because $`V^{}(r)\frac{e^2}{ϵr}`$ for large $`r`$. The $`s=2`$ LL can be similarly treated, and the corresponding effective potential will be denoted by $`V^{\prime \prime }(r)`$. The pseudopotentials for the Coulomb interaction in the $`s=1`$ and $`s=2`$ LLs and for the corresponding effective potentials $`V^{}(r)`$ and $`V^{\prime \prime }(r)`$ in the $`s=0`$ LL are shown in Fig. (1), demonstrating that the lowest LL problem with the effective interaction is an excellent approximation to the higher LL problem with Coulomb interaction. With this interaction, we then evaluate the energy of the CF exciton using the quantum Monte Carlo method developed earlier for the lowest LL wave functions. Because the pseudopotential are matched for the planar geometry, our results in the spherical geometry are meaningful only for sufficiently large systems; unphysical results may be obtained for small systems due to finite size effects. We will work below with systems containing as many as $`N=66`$ particles, using an efficient updating method discussed earlier.

We first consider 2CFs, corresponding to FQHE at $`\nu =2s+\frac{n}{2n+1}`$. Fig. (2) shows the energy of the CF exciton with $`n=1`$, 2, and 3 in the three lowest LLs ($`s=0`$, 1, and 2), for 66 particles. The energies are quoted in units of $`e^2/ϵł_0`$, where $`l_0=\sqrt{\mathrm{}c/eB}`$ is the magnetic length at $`\nu `$.

All FQHE is unstable in $`s=2`$ for 2CFs. We expect that this would remain the case in still higher LLs. The absence of FQHE in $`s2`$ is consistent with earlier theoretical studies, which have made a convincing case for either a bubble crystal or a stripe phase in third or higher ($`s2`$) LLs. At half filling ($`\nu _0=1/2`$), the wave vector of the stripe phase was estimated to be $`ql_02.4/\sqrt{2s+1}`$, which is also approximately equal to the reciprocal lattice vector associated with the bubble crystal away from half filling. Even though the instability occurs for wave vectors below $`ql_01`$, it is unfortunately not possible to determine from the dispersions shown in Fig. (2) a single wave vector for the instability, which precludes us from ascertaining from our method the reciprocal lattice vector of the true charge density wave ground state. Note that a comparison between the energies of the 1/3 FQHE state and the Wigner crystal state in the third LL does not indicate a lack of FQHE here, which is understandable in view of the fact that the instability is into a more complicated bubble crystal.

Surprisingly, we find the FQHE to be unstable also in the second LL. The only exception is $`\nu _0=1/3`$, which corresponds to $`\nu =7/3`$ in experiments (and, of course, other states related to it by symmetry). We have confirmed that $`\nu =7/3`$ remains stable in the thermodynamic limit by extrapolating the $`kl_00`$ energy to the $`N^10`$ limit. Given that, by its very design, our approach is expected to overestimate the strength of a FQHE state, the results provide strong evidence against FQHE at 12/5 for the pure Coulomb interaction. There is often a minimum observed in $`\rho _{xx}`$ in the vicinity of $`\nu =12/5`$, which may suggest an incipient FQHE state here (although no plateau has been observed yet). Can a FQHE state here be stabilized by LL mixing or finite thickness effects? LL mixing is a weak effect at large magnetic fields, and is expected to have the opposite effect, because it effectively screens the short range part of the interaction, thereby further weakening a given FQHE state. (It is possible, however, that at relatively low magnetic fields, when the screening is strong, the modification of the interaction due to screening by the lower Landau level might alter this conclusion. We have not investigated this question.) We have considered the effect of the transverse width of the electron wave function by employing the Fang-Howard variational form for the effective interaction, appropriate for the triangular well geometry. The instability is weakened, but not eliminated for typical paramters.

The instability in the second Landau level (Fig. 2) appears to occur at a small wave vector, which suggests a uniform compressible ground state (although a charge density wave state with a large lattice spacing can obviously not be ruled out on account of the finite size of our study). The nature of the compressible state in the second LL is not fully understood at present.

What about the other flavors of composite fermions? There are theoretical indications  that these are more stable in higher LLs than the 2CFs. In order to explore this issue further, we have computed the dispersion of the CF exciton for the 4CFs and 6CFs at $`\nu =2s+1/5,2s+2/9,2s+3/13`$ and $`2s+1/7`$ for $`s=0,`$ 1, and 2, shown in Fig. (3). Indeed, these states are more stable in higher Landau levels, with the 4CF states having the largest roton gap (in units of $`e^2/ϵl_0`$) in the second LL and the 6CF states in the third. Unlike for 2CFs, FQHE for 4CFs and 6CFs survives in the second and third LLs. As mentioned earlier, the trial wave function $`\mathrm{\Psi }`$ is quite close to the actual 1/5 ground state in the second LL, in fact more accurate than in the lowest LL, implying that the dispersion shown in Fig. (3) is quantitatively reliable. The 2/7 state is analogous to 1/5: whereas 1/5 is obtained from 1/3 by attachment of two additional vortices, 2/7 is similarly obtained from 2/3. There exists preliminary experimental evidence for both 1/5 and 2/7 in the second LL. There is an indication for 1/7 state in the lowest LL  but none yet in higher LLs. The observation of the 4CF and 6CF states in higher LLs is complicated by their rather small energy gaps (because, for a given density, $`e^2/ϵl_0`$ is much smaller here than at the corresponding lowest LL fraction), as well as their close proximity to strong integral quantum Hall plateaus.

Why is the FQHE for 2CFs unstable in higher Landau levels but relatively stable to changes in effective interaction within the lowest LL? The reason is that in order for the FQHE to occur, the interaction must not only be repulsive but it must have a sufficiently strong hard-core repulsion at short distances. It is useful to characterize the hard-core nature of the repulsion for a given interaction by the ratio $`V_1^{(s)}/V_3^{(s)}`$, given in Fig. (1). It appears that for $`V_1/V_31.3`$, most FQHE is unstable, with the exception of the 1/3 state which is only marginally stable. In contrast, when finite thickness is taken into account in the lowest LL, all pseudopotentials are suppressed more or less uniformly, and the ratio remains above 1.4-1.5 for typical experimental parameters. For 4CFs the pseudopotential $`V_1`$ is ineffective (provided it is not too small to cause an instability) because the wave functions for 4CFs are approximately given by the 2CF wave functions multiplied by a power of Jastrow factor that eliminates the unit relative angular momentum. Similarly, $`V_1`$ and $`V_3`$ are not relevant for 6CFs. Therefore, the appropriate ratios are $`V_3/V_5`$ and $`V_5/V_7`$ for 4CFs and 6CFs, respectively. These, however, increase in going from the lowest to the higher LLs, as shown in Fig. (1), thus explaining why 4CFs and 6CFs have a qualitatively different LL-index dependence as compared to the 2CFs. There is a close correspondence between the appropriate ratio and the stability, as a comparison of Figs. (1) and (3) shows, with 4CFs (6CFs) being strongest in the $`s=1`$ ($`s=2`$) LL.

In conclusion, our study provides an insight into the paucity of FQHE in higher Landau levels, as well as the qualitative difference between the stabilities of composite fermions of various flavors in different Landau levels. Most FQHE states of 2CFs in higher Landau levels are intrinsically unstable due to a collapse of the neutral excitation, which in turn is triggered by a softening of the short range part of the interaction.

This work was supported in part by the National Science Foundation under grant no. DMR-9986806. We are grateful to the Numerically Intensive Computing Group led by V. Agarwala, J. Holmes, and J. Nucciarone, at the Penn State University CAC for assistance and computing time with the LION-X cluster. We thank J.P. Eisenstein, V.J. Goldman, E. Rezayi, and H.L. Stormer for helpful conversations.


text:
# 1 The disputed diagram.

Comment on “Complete One-Loop Analysis of the Nucleon’s Spin Polarisabilities”

In a recent letter,  Gellas et al. published results for the spin polarisabilities of the nucleon calculated to next-to-leading order (NLO) in heavy baryon chiral perturbation theory (HB$`\chi `$PT). Their results differ from those in Refs. . Although this is due to different definitions and not to discrepancies in computing Feynman diagrams, it has caused some confusion in the community. In this comment we show that it is the definition of Ref. , and not that of Gellas et al., which should be compared with the values extracted from dispersion relations (DR’s) .

The disagreement stems from the standard procedure of excluding the nucleon pole graphs from the definition of nucleon polarisabilities. In particular, in fixed-$`t`$ dispersion relations each of the six independent Compton amplitudes is split into two parts, a pole piece and a contribution from an integral over the imaginary part of the amplitude above the $`\pi N`$ cut; the latter is obtained from photoproduction data. The polarisabilities $`\gamma _i`$ are then defined as the lowest terms in the energy expansion of the second (non-pole) piece of the real part of the amplitude. The pole occurs when the intermediate nucleon is on shell. The full pole contribution can be calculated using Dirac spinors with the on-shell $`\gamma N`$ coupling, and depends only on the nucleon charge, mass and magnetic moment. There are no “off-shell” or “sideways” form-factors involved, and so it is gauge-invariant, Lorentz-covariant, and satisfies the DR’s. The low energy theorems (LETs) for Compton scattering are all pole terms.

The difference between the results of Refs.  is in the treatment of the contribution from Fig. 1. If this is expanded in powers of the photon energy, $`\omega `$, the leading term is needed to renormalise the bare value of the anomalous magnetic moment in the leading order (LO) pole terms, as might be expected. (Less obviously, many other graphs also have leading terms which contribute to this correction.) The disagreement between Gellas et al. and us is in the treatment of the $`\omega ^3`$ terms from Fig. 1. We include them in the polarisabilities, they do not.

It is clear, however, that as the pole terms in the DR’s depend only on the on-shell couplings, any part of the contribution from Fig. 1 (or other loop diagrams) which is not purely a renormalisation of the mass or magnetic moment is not part of the pole term. Thus it should be included in the definition of the polarisabilities if comparison is to be made with the values extracted from DR’s. Indeed by Cauchy’s theorem the polarisability contributions from the amplitude of Fig. 1 are related to an integral over the its imaginary part, which arises from intermediate on-shell $`\pi N`$ states—just as in the similar graph where the second photon couples to the nucleon inside rather than outside the pion loop.

For forward scattering, the polarisability $`\gamma _0`$, can be directly related to the spin-dependent total inelastic cross section via a GDH-style DR. Clearly all diagrams contribute to this. However as the pole contribution beyond the LET vanishes for forward (and backward) angles, $`\gamma _0`$ will be the same whether defined this way or via fixed-$`t`$ DR’s.

We stress that this is not a problem of relativistic versus heavy-baryon formulations. If however, like Gellas et al., one wishes to propose a new definition of polarisabilities excluding one-particle reducible graphs—polarisabilities which could not then, as we have argued, be directly compared with values extracted from DR’s—then care must be taken that the definition is Lorentz and gauge invariant. In HB$`\chi `$PT the procedure of Gellas et al. violates Lorentz invariance to the order in the expansion in powers of $`1/m_N`$ to which they are working . Its gauge invariance has also not been demonstrated.

It may be worth stressing here that while HB$`\chi `$PT will not give the pole at any finite order, it can reproduce the expansion of the amplitude in powers of $`1/m_N`$. The true pole does not occur at physically realisable values of the energy and momentum transfer, nor is it so close as to make the expansion badly behaved. Thus HB$`\chi `$PT can calculate the physical amplitude in the vicinity of zero energy without any problems of principle. (Issues of chiral convergence or the effects of the delta are irrelevant to this discussion, though they are important in practice.)

It is true that there are significant discrepancies between the polarizabilities from our NLO HB$`\chi `$PT calculations and the DR analyses . Since the NLO corrections are of the same magnitude as the LO values however, there is no reason to expect that the series has converged.

Michael C. Birse, Xiangdong Ji and Judith A. McGovern

> Department of Physics and Astronomy, University of Manchester, M13 9PL, U.K.
> Department of Physics, University of Maryland, College Park, MD20742

PACS numbers: 12.39Fe 13.60Fz 11.30Rd