text:
# Radiative transfer in disc galaxies – III. The observed kinematics of dusty disc galaxies

## 1 Introduction

During the last decade of the past century there has been a vivid discussion about the opacity of spiral galaxies. This discussion was initiated by Disney, Davies & Phillipps (1989) and Valentijn (1990), who countered the conventional view that spiral galaxies are optically thin over their entire optical discs (Holmberg 1958; de Vaucouleurs, de Vaucouleurs & Corwin 1976; Sandage & Tammann 1981). This issue has not yet been settled, although an overall consensus seems to be emerging that spiral galaxies are, in general, at least moderately optically thick within their optical disc, i.e. that they have central face-on optical depth of at least unity. We will not repeat a detailed overview of the different points of view and the large variety of observational tests, and refer to Davies & Burstein (1995), Xilouris et al. (1997), Kuchinski et al. (1998) and Calzetti (2001) for excellent reviews. We just wish to point out an important observational constraint which, in our opinion, is not always given enough attention: the far-infrared (FIR) emission of spiral galaxies. This emission originates from the absorption of starlight by dust grains, which is thermally re-emitted at longer wavelengths. The most recent results show that about 30 per cent of the bolometric luminosity of late-type galaxies is FIR emission by interstellar dust (Popescu & Tuffs 2002). Notice that this figure refers to normal quiescent spiral galaxies; special classes of galaxies such as starburst galaxies emit even higher fractions of their energy at FIR wavelengths. These facts demonstrate quite convincingly that interstellar dust forms an important constituent of spiral galaxies.

The physical processes of absorption and scattering by dust grains should hence be taken into account when interpreting and modelling the observed optical properties of galaxies. We would like to stress that all optical observables will to some extent be affected by dust, including the observed kinematics. Indeed, when interstellar dust grains absorb or scatter photons, the kinematical information contained within these photons, in the form of the Doppler shift, will also be absorbed or scattered. It is important to realize that most of the kinematical data of galaxies are optical measurements. Firstly, stellar kinematics are always measured from optical absorption lines. In principle, the attenuation effects of interstellar dust could be minimized by using the sensitive 12CO absorption feature at 2.29 $`\mu `$m (Gaffney, Lester & Doppmann 1995). Measuring stellar kinematics at near-infrared (NIR) wavelengths is, however, observationally extremely difficult due to the high sky surface brightness. Secondly, the vast majority of spatially resolved rotation curves of spiral galaxies have been measured using optical emission lines (usually the H$`\alpha `$ line). They can also be measured using the Hi or CO lines, where dust extinction is not an issue. However, these methods are observationally much more expensive than the more straightforward H$`\alpha `$ observations.

Currently, the only way to investigate the effects of dust attenuation on the observed kinematics of galaxies is by means of detailed radiative transfer modelling. Many different approaches exist to handle the radiative transfer problem, and various teams have adopted different techniques to study the effect of dust on the photometry and spectral energy distributions (SEDs) of galaxies (e.g. Bruzual, Magris & Calvet 1988; Witt, Thronson & Capuano 1992; Byun, Freeman & Kylafis 1994; Wise & Silva 1996; Bianchi, Ferrara & Giovanardi 1996; Corradi, Beckman & Simonneau 1996; Baes & Dejonghe 2001a). The effect of dust attenuation on the observed kinematics, however, is largely unexplored, and limited to the rotation curves of spiral galaxies. Davies (1990) argued that the effects of dust absorption on the observed optical rotation curve of a disc galaxy can lead to a severe underestimation of the true rotational velocity in the inner regions. This argument was quantified by Bosma et al. (1992), who calculated the effect of dust absorption on the observed rotation curve of edge-on galaxies. The most detailed work on this subject comes from Matthews & Wood (2001), who calculate the effects of dust attenuation on the rotation curve of low surface brightness (LSB) galaxies through Monte Carlo simulations.

We have embarked on a program to systematically investigate the effects of dust on the observed kinematics of galaxies. In previous papers we focused on the effects of dust attenuation on the observed stellar kinematics of elliptical galaxies. Initially, we only considered absorption by dust grains in our models (Baes & Dejonghe 2000; Baes, Dejonghe & De Rijcke 2000). The inclusion of absorption in the calculation of the observed kinematics is a fairly simple procedure: it basically adds a weight to the contribution of each individual star along a line of sight. The radiative transfer problem becomes much more complicated however when scattering is also included. In this case, photons can leave their original path, such that any star can contribute to the observed kinematics in any line of sight. Moreover, not only the stellar velocities but also the individual dust grain velocities should be taken into account in the calculation of the observed kinematics. We argue that the Monte Carlo method is the only radiative transfer modelling method in which kinematical information can be included in an elegant and straightforward way. Using a one-dimensional spherical Monte Carlo code, we found that the effects of dust scattering are fairly important: dust attenuation can serve as an additional or alternative explanation for the stellar kinematical evidence of a dark matter halo around ellipticals (Baes & Dejonghe 2001b, 2002a).

In this paper, we tackle the more complicated problem of investigating the effect of dust attenuation on the observed kinematics of disc galaxies. We have written a new Monte Carlo code (SKIRT, acronym for Stellar Kinematics Including Radiative Transfer) that can handle any geometry of stars and dust. This code is presented in section 2, and tested in section 3. We present our model and the results of our Monte Carlo simulations in section 4. These results are discussed in section 5, and section 6 presents our conclusions.

## 2 Description of the SKIRT code

The basic characteristics of Monte Carlo radiative transfer have been explained at length by various authors (e.g. Cashwell & Everett 1959; Mattila 1970; Witt 1977; Fischer, Henning & Yorke 1994; Bianchi et al. 1996). In essence, a Monte Carlo radiative transfer code follows the life of a very large number of individual photons. A photon is, at each stage in its existence, characterized by various quantities, such as its position $`𝒓`$, propagation direction $`𝒌`$, wavelength $`\lambda `$ and Doppler shift $`u`$ as measured in a frame of rest.111Throughout this paper, we will make no distinction between the measured Doppler shift $`\mathrm{\Delta }\lambda `$ and the corresponding line-of-sight velocity $`u`$, which are related by the expression $`\mathrm{\Delta }\lambda /\lambda _0=u/c`$. After being emitted, photons propagate on straight lines through the interstellar medium until they either interact with a dust grain or leave the galaxy. The various interactions alter the properties of the photon, according to random numbers generated from the appropriate probability functions. When at last the photon escapes from the galaxy, its final properties are recorded. After recording a large numbers of photons in this way, the global observed properties of the system can be calculated.

We will not repeat an in-depth description of the principles and the numerous equations of Monte Carlo radiative transfer here, as they are well described in the articles listed above. Instead, we will focus on a number of aspects which make our code different from the existing ones. These differences are mainly our choice of the dust grid and the possibility to include kinematical information into the radiative transfer calculations.

### 2.1 The emission process

When no kinematical information is included in the radiative transfer calculations, photons are characterized at each moment by a wavelength $`\lambda `$, a position $`𝒓`$ and a direction $`𝒌`$. Initial values for these quantities must be generated randomly from the emitting stellar system, which can be composed of several stellar components.222We describe the code as adopted for radiative transfer calculations in a dusty galaxy. The code can, however, equally well be applied to any other environment. Terms like “stellar component” must therefore be interpreted in a broad way, and can also refer to a single star, emitting gas, an AGN, etc., i.e. any source of photons of high enough energy. The initial emission direction $`𝒌_0`$ can be sampled from the unit sphere, the initial position $`𝒓_0`$ is sampled from the spatial distribution of the stellar component, and the initial wavelength $`\lambda _0`$ is sampled from the SED of the stellar component at the position $`𝒓_0`$. The SKIRT code contains a library with a set of common stellar components (including spherical Jaffe, Hernquist, Plummer and de Vaucouleurs models, and axisymmetric exponential, sech, isothermal disc models), from which it is straightforward to sample a random position (see Appendix A). Another library contains a set of SEDs, including Planck functions and realistic tabulated stellar and galaxy SEDs (Kinney et al. 1996; Pickles 1998). This library also contains so-called monochromatic SEDs, which are described by a Dirac delta function. Sampling a random wavelength from the latter is trivial, because all photons are emitted with the same wavelength.

When the SKIRT code is used to calculate the observed kinematics of a dusty galaxy, photons must also carry with them the kinematic signature of the star that has emitted them, in the form of a Doppler shift. Because this is a one-dimensional quantity, not all kinematical information of the star can be transmitted to the observer, and it is important to consider how the observed Doppler shift relates to the velocity vectors of the star and the scattering dust grains. This problem was discussed in detail in section 2.3 of Baes & Dejonghe (2002a), where it is shown that the general expression for the Doppler shift after $`M`$ scattering events is

$$u_M=𝒗_{}𝒌_0+\underset{i=1}{\overset{M}{}}𝒗_{\text{d}_i}(𝒌_i𝒌_{i1}),$$
(1)
where $`𝒗_{}`$ and $`𝒗_{\text{d}_i}`$ are the velocity vectors of the emitting star and the $`i`$th dust grain respectively. The initial value of the Doppler shift attached to the photon is hence the component of the star in the direction of the original emission, $`u_0=𝒗_{}𝒌_0`$.

To complete the initialization of a photon, we also need to generate a line-of-sight velocity from the appropriate probability distribution. In this case, this function is the spatial line-of-sight velocity distribution (LOSVD) $`\varphi _{}(𝒓_0,𝒌_0,u)`$, which describes the probability for a star at a position $`𝒓_0`$ to have a velocity component $`u`$ in the direction $`𝒌_0`$. In general, the spatial LOSVD is a marginal probability distribution of the stellar distribution function $`F_{}(𝒓,𝒗)`$, which describes the probability density of the stellar component in six-dimensional phase space,

$$\varphi _{}(𝒓,𝒌,u)=\frac{1}{n_{}(𝒓)}F_{}(𝒓,𝒗)\text{d}v__1\text{d}v__2,$$
(2)
where $`n_{}(𝒓)`$ is the stellar number density and the integration covers the entire velocity space perpendicular to the direction $`𝒌`$. Thus, in order to be of use for a kinematical Monte Carlo run, each stellar component must allow the generation of a velocity from its spatial LOSVDs for arbitrary $`𝒓`$ and $`𝒌`$.

The spatial LOSVDs can be calculated analytically for a number of self-consistent spherical models. Apart from isotropic models (e.g. Plummer 1911; Hénon 1959; Jaffe 1983; Hernquist 1990), this is possible for the anisotropic families of Plummer and Hernquist models described by Dejonghe (1987) and Baes & Dejonghe (2002b). For more complicated dynamical models, however, an analytical evaluation of the spatial LOSVDs is not possible. For such models, we approximate the velocity distribution at each position in the galaxy as a local gaussian distribution. This means that at each position $`𝒓`$, the distribution function can be written as

$$\begin{array}{c}F_{}(𝒓,𝒗)=\frac{n_{}(𝒓)}{(2\pi )^{3/2}\sigma _1\sigma _2\sigma _3}\hfill \\ \hfill \times \mathrm{exp}\left[\frac{(v_1\overline{v}_1)^2}{2\sigma _1^2}\frac{(v_2\overline{v}_2)^2}{2\sigma _2^2}\frac{(v_3\overline{v}_3)^2}{2\sigma _3^2}\right].\end{array}$$
(3)
Here $`(\overline{v}_1,\overline{v}_2,\overline{v}_3)`$ and $`(\sigma _1,\sigma _2,\sigma _3)`$ represent the mean velocities and velocity dispersions in the directions $`(𝒆_1,𝒆_2,𝒆_3)`$, the principle axes of the velocity ellipsoid at the position $`𝒓`$. The orientation of the velocity ellipsoid and the values of the mean velocities and velocity dispersions can vary from position to position, such that the distribution function (3) represents a fairly general distribution function. But a distribution function of this form has the advantage that the spatial LOSVDs can be calculated exactly for any direction $`𝒌`$ by applying the formula (2). We find that the spatial LOSVD will also be a gaussian distribution,

$$\varphi _{}(𝒓,𝒌,u)=\frac{1}{\sqrt{2\pi }\sigma _u}\mathrm{exp}\left[\frac{(u\overline{u})^2}{2\sigma _u^2}\right],$$
(4)
with parameters

$$\overline{u}=k_1\overline{v}_1+k_2\overline{v}_2+k_3\overline{v}_3,$$
(5)
$$\sigma _u^2=k_1^2\sigma _1^2+k_2^2\sigma _2^2+k_3^2\sigma _3^2,$$
(6)
where obviously $`(k_1,k_2,k_3)`$ are the components of the vector $`𝒌`$ with respect to the orthonormal reference system $`(𝒆_1,𝒆_2,𝒆_3)`$. If hence, for each stellar component, we specify the orientation and parameters of the velocity ellipsoid at each position $`𝒓`$, we can easily sample line-of-sight velocities from the spatial LOSVDs in an arbitrary direction, and hence initialize the photons’ initial Doppler shifts.

### 2.2 The dust iteration

The life of a single photon can be thought of as a loop, whereby at each iteration (representing a physical process) we must update its position $`𝒓`$, propagation direction $`𝒌`$, wavelength $`\lambda `$ and Doppler shift $`u`$. The initial values are determined randomly at the emission phase, and change at every scattering event. As we only consider coherent scattering, and we are not taking into account the re-emission of photons at longer wavelength333For the present study, we are primarily interested in the attenuation of starlight by interstellar dust at optical wavelengths, where the contribution of re-emitted photons is negligible., the only wavelength change of a photon along its path is due to the varying Doppler shift. These wavelength variations are so tiny that the optical properties of the dust (extinction curve, dust albedo and asymmetry parameter) do not vary. The wavelength of a photon can hence be considered fixed throughout its lifetime. We must therefore only update the three quantities $`𝒓`$, $`𝒌`$ and $`u`$ at each iteration, from their old values ($`𝒓_{i1}\mathrm{}`$) to the new ones ($`𝒓_i\mathrm{}`$). This proceeds in several steps.

The first step in the iteration consists of determining whether an interaction with a dust grain will take place or whether the photon will leave the galaxy. To do this, we sample an optical depth $`\tau _\lambda `$ from an exponential distribution and compare it to $`\tau _{\text{path},\lambda }`$, the total optical depth along the path. When $`\tau _\lambda >\tau _{\text{path},\lambda }`$, the photon will leave the galaxy; in the other case, the photon will interact with a dust grain. This interaction can either be a scattering or an absorption event, which is determined by the scattering albedo. In the former case, the photon will continue its journey through the galaxy, in the latter case, the loop is ended.

The second step in the loop, if the photon is scattered, is the determination of the position of the scattering. Therefore we have to translate the sampled optical depth $`\tau _\lambda `$ to a physical path length $`s`$. With this path length known, the new position is $`𝒓_i=𝒓_{i1}+s𝒌_{i1}`$.

The third step in the iteration is the determination of the new propagation direction $`𝒌_i`$ of the photon. It is found by sampling a direction from the probability density $`p(𝒌_i)=\mathrm{\Phi }_\lambda (𝒌_{i1},𝒌_i)`$, which represents the scattering phase function. Various phase functions are built into the SKIRT code, including the isotropic phase function and the anisotropic Henyey-Greenstein phase function. The sampling of a direction from these phase functions can be performed analytically.

The final step is to update the photon’s Doppler shift. The relative orientation of the propagation directions before and after the scattering event cause a change in Doppler shift \[from equation (1)\]

$$u_iu_{i1}=𝒗_{\text{d}_i}(𝒌_i𝒌_{i1})$$
(7)
To calculate the updated Doppler shift, we need to sample a dust grain velocity from the dust velocity distribution. Therefore, the dust components in the SKIRT library must not only be specified by a spatial distribution, but also by the entire phase space distribution $`F_\text{d}(𝒓,𝒗)`$. Analogous to the stellar distribution function, we assume that the dust velocity field can be described by a trivariate gaussian distribution. Because the dispersions in the interstellar medium are fairly small compared to the rotational velocities however, we restrict ourselves to assuming an isotropic dispersion tensor, such that

$$F_\text{d}(𝒓,𝒗)=\frac{n_\text{d}(𝒓)}{(2\pi )^{3/2}\sigma _\text{d}^3}\mathrm{exp}\left[\frac{(𝒗\overline{𝒗})^2}{2\sigma _\text{d}^2}\right],$$
(8)
with $`\overline{𝒗}`$ and $`\sigma _\text{d}`$ being the mean velocity vector and the velocity dispersion of the dust. Similarly as in the stellar case, we don’t have to sample a full three-dimensional velocity vector for the dust, but only a component in the direction

$$𝒌_i^{}=\frac{𝒌_i𝒌_{i1}}{𝒌_i𝒌_{i1}}.$$
(9)
This means that we just have to sample a velocity from a one-dimensional gaussian distribution with mean velocity $`\overline{u}=\overline{𝒗}𝒌_i^{}`$ and dispersion $`\sigma _u=\sigma _\text{d}`$.

### 2.3 Integration through the dust

The first step in the iteration described in the previous subsection requires the calculation of the total optical depth along a given path. It is found through the integral

$$\tau _{\text{path},\lambda }(𝒓,𝒌)=_0^{\mathrm{}}\kappa _\lambda (𝒓+s𝒌)\text{d}s,$$
(10)
where $`\kappa _\lambda (𝒓)`$ represents the total (absorption plus scattering) opacity of the dust at a position $`𝒓`$. To execute the second step in the iteration, we need to convert an optical depth $`\tau _\lambda `$ into a physical path length $`s(𝒓,𝒌,\tau _\lambda )`$ along a given path, which is done by solving the equation

$$_0^s\kappa _\lambda (𝒓+s^{}𝒌)\text{d}s^{}=\tau _\lambda $$
(11)
for $`s=s(𝒓,𝒌,\tau _\lambda )`$. Only for a few simple geometries, such as constant density or some spherical dust distributions, can these integrations be performed analytically. In general, the calculation of the two quantities $`\tau _{\text{path},\lambda }(𝒓,𝒌)`$ and $`s(𝒓,𝒌,\tau _\lambda )`$ must be performed numerically. These integrations are usually the most time-consuming part of Monte Carlo codes. Therefore, it is important to think about an efficient and flexible way to do this; we have considered two approaches.

#### 2.3.1 The UDD grid

In order to allow a completely arbitrary distribution of dust (including clumpy dust distributions), most modern radiative transfer Monte Carlo codes adopt an approach that can be described as a uniform dust density (UDD) grid. It consists of dividing space into a number of cells, and attaching a uniform dust density to each cell, for example the value at the centre of the cell. This uniform density (and hence opacity) makes it very easy to calculate the optical depth \[equation (10)\] along a given path: simply calculate the distance that the photon runs through each cell it crosses, multiply this distance with the opacity in the cell, and add all the pieces together. The inversion of the equation (11) is not much more difficult. First, look for the cell in which the interaction will take place. Since the opacity increases linearly with path length within a single cell, this inversion of the equation is then straightforward.

In the SKIRT code, a UDD grid is included based on a three-dimensional cartesian grid. The grid cells can be chosen arbitrarily, in order to achieve an efficient grid for a variety of dust distributions. An algorithm is included to construct clumpy dust distributions from an underlying smooth dust distribution, as described in Witt & Gordon (1996) and Bianchi et al. (2000).

#### 2.3.2 The IDD grid

In addition to the UDD grid, we explore a new approach to handle the integration through the dust: an interpolated dust density (IDD) grid. The basis is the same: we divide space into a number of cartesian grid cells, but, instead of attaching a uniform opacity to each cell, we use the correct values of the opacity at the eight border points of the cells, and we apply a simple trilinear interpolation routine to determine the opacity in each of the points within a grid cell. Consider a photon at the position $`𝒓=(x,y,z)`$, a position within the $`(i,j,k)`$’th cell with border points $`x_ixx_{i+1}`$, etc. We define the dimensionless quantities $`p`$, $`q`$ and $`r`$ as

$$p=\frac{xx_i}{x_{i+1}x_i},q=\frac{yy_j}{y_{j+1}y_j},r=\frac{zz_k}{z_{k+1}z_k}.$$
(12)
Since $`p`$ increases linearly from 0 to 1 as $`x`$ moves linearly from $`x_i`$ to $`x_{i+1}`$, and similarly for $`q`$ and $`r`$, the trilinear interpolation rule for the opacity simply reads

$`\kappa _\lambda (𝒓)`$ $`=\left(1p\right)\left(1q\right)\left(1r\right)\kappa _{i,j,k}`$
$`+\left(1p\right)\left(1q\right)r\kappa _{i,j,k+1}`$
$`+\mathrm{}`$
$`+pqr\kappa _{i+1,j+1,k+1},`$ (13)
where $`\kappa _{i,j,k}`$, etc. represent the opacity in the border points. Now consider the fraction of the path that remains in this cell. This path can be parametrized using the expression $`𝒓+s𝒌`$, with the path length $`s`$ as a free parameter. In each of the points along this path (as long as we remain in the same cell), we can calculate the opacity by a similar formula as (13),

$`\kappa _\lambda (𝒓+s𝒌)`$ $`=\left[1p(s)\right]\left[1q(s)\right]\left[1r(s)\right]\kappa _{i,j,k}`$
$`+\left[1p(s)\right]\left[1q(s)\right]r(s)\kappa _{i,j,k+1}`$
$`+\mathrm{}`$
$`+p(s)q(s)r(s)\kappa _{i+1,j+1,k+1},`$ (14)
where the parameters $`p`$, $`q`$ and $`r`$ are now linear functions of the path length $`s`$,

$$p(s)=\frac{x+sk_xx_i}{x_{i+1}x_i},$$
(15)
$$q(s)=\frac{y+sk_yy_j}{y_{j+1}y_j},$$
(16)
$$r(s)=\frac{z+sk_zz_k}{z_{k+1}z_k}.$$
(17)
Along a path within a single cell, the opacity is thus a cubic polynomial in $`s`$, the distance traveled along the path,

$$\kappa _\lambda (𝒓+s𝒌)=\underset{m=0}{\overset{3}{}}a_ms^m,$$
(18)
with constant coefficients $`a_m`$. With this expression for the opacity, the calculation of the optical depth $`\tau _{\text{path},\lambda }(𝒓,𝒌)`$ is very straightforward. Indeed, it is sufficient to calculate, for each cell the path passes through, the two boundary positions. This enables us to calculate the path length $`s`$ within the cell and the coefficients $`a_m`$. Given these coefficients, the portion of the optical depth in that cell is found by substitution into the quartic polynomial

$$_0^s\kappa _\lambda (𝒓+s^{}𝒌)\text{d}s^{}=\underset{m=1}{\overset{4}{}}\frac{a_ms^m}{m}.$$
(19)
Also, the calculation of the path length $`s(𝒓,𝒌,\tau _\lambda )`$ is straightforward. First, we determine the cell in which the optical depths at the border points bracket the value of $`\tau _\lambda `$, in a similar way as described above to calculate the total optical depth along the path. The determination of $`s`$ then comes down to finding the root of a quartic equation. The obvious way of finding this root is by means of a Newton-Raphson iteration, because we can easily compute both the function values and the derivatives of this function. This routine has a very fast quadratic convergence (the number of significant digits approximately doubles with each iteration). Usually, no more than two iterations are necessary.

#### 2.3.3 Comparison of the grids

To compare the accuracy and efficiency of the two approaches, we applied both techniques to a simple spherical dust geometry, characterized by an opacity

$$\kappa _\lambda (𝒓)\left(1+\frac{r^2}{c^2}\right)^{\alpha /2}.$$
(20)
The integrations through the dust, i.e. the calculation of the two functions $`\tau _{\text{path},\lambda }(𝒓,𝒌)`$ and $`s(𝒓,𝒌,\tau _\lambda )`$, can be performed exactly for such a dust component with $`\alpha =2`$, 3 or 5 \[see Baes (2001) for details\]. These exact results were compared with the results from the UDD and IDD grids. In figure 1 we plot the relative error of the calculated values of the opacity, the total optical depth along a path and the path length corresponding to a given optical depth for 1000 randomly generated photons for such a model. For a given number of dust grid cells, the accuracy of the integrations in the IDD grid is much better than those in the UDD grid. However, in terms of computational efficiency, comparing the two grids with the same numbers of cells is not fair, because the integrations in the UDD grid are less complicated and hence faster. We found that the IDD grid integrations are a factor 3 slower in the mean. On the other hand, in order to reach a similar accuracy as the IDD grid, many more cells are necessary in the UDD grid. Apart from slowing down the integrations, this will also increase the memory requirements significantly. After several tests, we found that the IDD grid is computationally preferable above the UDD grid, in particular when the system harbours a large dynamical range of dust densities (such as an exponential profile which is appropriate in disc galaxies). It should be noted however, that the UDD grid still remains a very useful alternative, in particular to include clumpy dust distributions into the SKIRT code.

### 2.4 Optimization of the code

The procedure described at the beginning of the previous section is the most basic version of a Monte Carlo iteration step. When executed as such, however, it would be very inefficient. Throughout the years a number of ’intelligent tricks’ have been presented which increase the computational efficiency considerably. Most of them attach a weight to each photon, which can be altered during its lifetime. We incorporate three such techniques into our SKIRT code.

The first optimization to the basic Monte Carlo scenario is to turn all interactions into scattering events. In reality, of course, the interaction between a photon and a dust grain can either be an absorption or a scattering event; the nature of this event is determined by the albedo $`\omega _\lambda `$ of the dust grains. The problem with absorptions is that the photon is lost, and therefore does not contribute to the observed radiation field anymore. To overcome this problem, we force each interaction to be a scattering event and alter the weight of the photon after each scattering by a factor $`\omega _\lambda `$ in order to compensate for the fraction of absorbed photons.

The second technique is the peeling-off procedure introduced by Yusef-Zadeh et al. (1984). In the basic Monte Carlo iteration, each photon leaves the galaxy at a certain stage, and it can be classified in bins according to its position and direction. This approach has a number of disadvantages. Firstly, if we are interested in only one single observing position, most of the photons are not used. Secondly, if the observer has the misfortunate of being located in a direction in which not many photons leave the galaxy, he might have to wait for a very long time to detect enough photons to obtain good statistics. Thirdly, if one wants to compare the view of a system from two observing positions close to each other, the direction bins must be chosen to be very small, which requires a very large number of photons. The peeling-off technique can be used to tackle these problems. It consists of creating a new photon during each scattering \[resp. emission\] process, which is scattered \[resp. emitted\] exactly in the direction of the observer. This photon will be detected by the observer. Its weight is altered by a factor $`\mathrm{\Phi }(𝒌,𝒌_{\text{obs}})`$ \[resp. 1\] to compensate for the probability of scattering \[resp. emission\] in the direction $`𝒌_{\text{obs}}`$, and a factor $`\mathrm{exp}[\tau _{\text{path},\lambda }(𝒓,𝒌_{\text{obs}})]`$ to compensate for extinction along the line-of-sight.

The third optimization we apply is the principle of forced first scattering (Witt 1977), which is useful in the relatively low optical depth regime. A problem in such systems is that most of the photons leave the galaxy directly without a single interaction with a dust grain, and therefore many photons need to be generated in order to have a good statistics of the scattered radiation. To overcome this problem, we split the emitted photon into two parts. The first part of the photon directly leaves the galaxy and obtains a weight $`\mathrm{exp}(\tau _{\text{path},\lambda })`$ (because we adopt the peeling-off technique, we don’t have to consider this fraction anymore). The other part of the photon with the remaining weight is forced to be scattered at least once, which is achieved by sampling an optical depth $`\tau _\lambda `$ such that $`\tau _\lambda <\tau _{\text{path},\lambda }`$. This technique is called the principle of forced first scattering, but nothing prevents us applying it only to freshly emitted photons: it can equally well be applied to scattered photons (Bianchi et al. 1996). The number of forced scatterings we apply in our calculations is usually 3, but this value can be set arbitrarily according to the studied system.

### 2.5 The detection and data reduction processes

The last phase in the Monte Carlo cycle is the detection phase, whereby the information contained in the photon is transmitted to the observer. Because of the peeling-off procedure, the number and position of the observers can be chosen randomly. At each observing position, the user can choose from several simulated instruments to detect the photons: photometers, low-resolution and high-resolution spectrographs. Each instrument consists of a number of pixels on the plane of the sky and they come in two geometries: rectangular and circular. Each pixel contains a detector, the nature of which depends on the kind of instrument.

The photometers are the most simple instruments, whereby the detector in each pixel is just a counter. Every time a photon hits the pixel and its wavelength is within the detector’s bandwidth, its weight is added to the intensity. At the end of the simulation, the 2D surface brightness image of the system is directly obtained. The photometers are designed such that each photon remembers the number of scatterings it has undergone, and this information is also used by the photometric detectors. This allows us to construct images for each individual scattering (direct emission, photons that have been scattered once, etc.), which can be useful to disentangle the effects of absorption and scattering.

The low-resolution spectrographs work in a very similar way, except that there is not just one detector, but a range of detectors corresponding to different wavelength bins. In this way a data cube is produced (actually a set of data cubes, one for every scattering), very similar to the data cubes obtained by radio or Fabry-Perot observations. The high-resolution spectrographs basically produce similar data cubes, with Doppler shift bins replacing the wavelength bins. Combining the different Doppler shift bins at a given position $`𝒙`$ on the sky, we obtain the LOSVD or line profile $`\varphi _\text{p}(𝒙,u)`$, which describes the entire projected velocity distribution. From the LOSVDs, the 2D mean projected velocity field $`v_\text{p}(𝒙)`$ and projected velocity dispersion field $`\sigma _\text{p}(𝒙)`$ can be calculated.

From a conceptual point of view, high-resolution and low-resolution spectrographs work in a very similar way. Both separate the incoming photons in different bins according to the extra information carried by the photon, which is the wavelength in the former and the Doppler shift in the latter case. From a computational point of view, however, there is an important difference between the two. Indeed, the wavelength of the photon is actively used throughout the Monte Carlo iteration, because the optical properties of the dust grains depend on the photon’s wavelength. On the contrary, a photon’s path does not depend on its Doppler shift. Moreover, the changes in Doppler shift experienced by a photon when it moves through the dust are independent of the initial stellar Doppler shift. Therefore, we can optimize our code in the following way. Instead of initializing the photon’s Doppler shift by generating a $`u_0`$ from the stellar spatial LOSVD, we postpone this contribution of the stellar velocity until later and initialize $`u_0`$ to zero. When the photon leaves the galaxy after $`M`$ scattering events, it will then carry a Doppler shift

$$u_M^\text{d}=\underset{i=1}{\overset{M}{}}𝒗_{\text{d}_i}(𝒌_i𝒌_{i1}),$$
(21)
where we still have to add the stellar Doppler shift in order to obtain the correct value. Instead of generating a single initial stellar Doppler shift $`u_0`$ from the spatial LOSVD $`\varphi _{}(𝒓_0,𝒌_0,u)`$, we split the photon over all velocity bins. Each fraction will then be weighted by a weight factor that represents the probability that the final Doppler shift of the photon will fall within the corresponding bin. In order to fall within a bin with boundaries $`u_k`$ and $`u_{k+1}`$, the photon must have an initial stellar Doppler shift satisfying $`u_ku_0+u_M^\text{d}u_{k+1}`$. The weight factor corresponding to the $`k`$th bin will hence equal

$$w_k=_{u_ku_M^\text{d}}^{u_{k+1}u_M^\text{d}}\varphi _{}(𝒓_0,𝒌_0,u)\text{d}u.$$
(22)
When we substitute the gaussian spatial LOSVD (4) into this expression, we can evaluate this weight factor exactly

$$w_k=\frac{1}{2}\left[\mathrm{erf}\left(\frac{u_{k+1}u_M^\text{d}\overline{u}}{\sqrt{2}\sigma _u}\right)\mathrm{erf}\left(\frac{u_ku_M^\text{d}\overline{u}}{\sqrt{2}\sigma _u}\right)\right].$$
(23)
This method is of the same kind as the optimization techniques described in section 2.4 and also increases the efficiency of the SKIRT code in a significant way.

### 2.6 General architecture and performance

The SKIRT code is implemented in the ANSI standard C++ language and makes optimal use of the object-oriented structure of this language. The output of the code consists of FITS files, which are created by means of the FITSIO and CCFITS packages.

The run time of the SKIRT code with reliable output depends critically on the number of photons, the number of observing directions, the spatial and velocity resolution of the instruments and the number of cells in the dust grid. For the simulations described in section 4, we found that about $`5\times 10^7`$ photons are necessary to obtain accurate results. The code spends 60 per cent of the time moving through the dust grid, while most of the remaining time is spent in the calculation of the spatial LOSVDs. When run on a modern Pentium iv PC, the run time of a single simulation (with 7 observing positions) is of the order of 50 to 100 hours.

## 3 Testing the code

The code was first run without dust included in order to test the projection algorithm and the calculation of the observed kinematics. We used simple non-rotating spherical geometries which allow an analytical expression of the surface brightness and projected velocity dispersion profiles (Dejonghe 1987; Hernquist 1990; Jaffe 1983, Baes & Dejonghe 2002b). Also the surface brightness profiles of face-on and edge-on exponential discs (which can be calculated analytically) are accurately recovered. Moreover, when we attach a simple velocity structure to an exponential disc, with a flat rotation curve and exponentially decreasing velocity dispersions, the observed kinematics can be also calculated analytically in the face-on and edge-on directions. These are recovered very well, proving the correct behaviour of the code.

Next, we compared the results of the SKIRT code with our previous Monte Carlo calculations of the observed kinematics of elliptical galaxies (Baes & Dejonghe 2001b, 2002a). The Monte Carlo code used in these papers also allowed for the calculation of the observed kinematics, but only for spherical models, chosen such that the integrations through the dust could be performed analytically. We found an excellent agreement between these results and those calculated with the new SKIRT code.

Finally, we compared our (photometric) results with the radiative transfer calculations of other authors. In figure 2 we plot the two-dimensional edge-on surface brightness distribution of an Sb galaxy model in the $`B`$ band for various values of the optical depth. The models are the same models as the BT0.3 models of Byun et al. (1994), and they consist of a stellar disc, a stellar bulge and a dust disc. Comparing this figure with their figure 4b, we find an excellent agreement between the two results. Notice that Byun et al. (1994) adopted a completely different approach to solve the radiative transfer problem.

## 4 The observed kinematics of dusty disc galaxies

We apply the SKIRT code to investigate the effect of dust absorption and scattering on the photometry and the stellar kinematics of normal galactic discs. We limit our modelling only to a stellar disc and do not include a bulge, because we are primarily interested in the disc kinematics.

### 4.1 Presentation of the model

Our model consists of a simple axisymmetric double exponential disc galaxy, with stellar emissivity

$$\mathrm{}_\lambda (𝒓)=\frac{L_\lambda }{4\pi h_{}^2z_{}}\mathrm{exp}\left(\frac{R}{h_{}}\right)\mathrm{exp}\left(\frac{|z|}{z_{}}\right).$$
(24)
The luminosity $`L_\lambda `$ is a free scaling parameter in our models, as both the input emissivity and the output light profile scale with the luminosity. For the scale length and scale height we chose the values adopted by Byun et al. (1994) to describe the Milky Way, $`h_{}=4`$ kpc and $`z_{}=350`$ pc. As we focus on the observed kinematics, which are always measured in a very narrow wavelength region, we chose a trivial monochromatic SED centered on the $`V`$ band.

In order to run a kinematical simulation for such a disc galaxy, we also need to provide details on the kinematical structure of the disc. We assume that the stellar velocities at each position in the disc can be represented by a trivariate Gaussian distribution,

$$F_{}(𝒓,𝒗)\mathrm{exp}\left[\frac{v_R^2}{2\sigma _R^2}\frac{(v_\phi \overline{v}_\phi )^2}{2\sigma _\phi ^2}\frac{v_z^2}{2\sigma _z^2}\right].$$
(25)
Such a velocity distribution was first proposed by Schwarzschild (1907) to describe the distribution of stellar velocities in the solar neighbourhood. The mean azimuthal velocity and the three velocity dispersion components can vary from position to position. The entire kinematical structure of the disc will be completely determined if we have expressions for the four functions $`\overline{v}_\phi (𝒓)`$, $`\sigma _R(𝒓)`$, $`\sigma _\phi (𝒓)`$ and $`\sigma _z(𝒓)`$. Since the system is axisymmetric, these quantities will depend only on the cylindrical radius $`R`$ and the height $`z`$. Moreover, because stellar discs are thin and the contribution of photons high above or under the symmetry plane of the galaxy will be quite small due to the vertical exponential decrease, we will assume that these four functions are a function of $`R`$ only.

For the $`\overline{v}_\phi `$ profile, we assume the arctan profile, a simple two-parameter fitting function that fits the observed rotation curves of disc galaxies fairly well (e.g. Courteau 1997),

$$\overline{v}_\phi (𝒓)=\frac{2v_{\text{max}}}{\pi }\mathrm{arctan}\left(\frac{R}{R_0}\right).$$
(26)
We have adopted an amplitude $`v_{\text{max}}=220`$ km/s and a scale radius $`R_0=0.5`$ kpc, which results in a rotation curve that has a steep initial gradient and reaches its flat part fairly quickly. The internal velocity dispersion profiles in disc galaxies are not well constrained. For a thin self-gravitating disc, the Jeans equations yield that the vertical velocity dispersion scales with the square root of the surface density. If we assume a constant mass-to-light ratio, we obtain the following parametrization,

$$\sigma _z(𝒓)=\sigma _{z0}\mathrm{exp}\left(\frac{R}{2h_{}}\right).$$
(27)
As all dispersions in a stellar disc are thought to exist from the gradual stellar heating of an initial cold disc, it makes sense to assume that the velocity dispersions in the other direction are intimately coupled to the vertical dispersion. Therefore, we assume a constant axis ratio of the velocity ellipsoid in our model, such that $`\sigma _R`$ and $`\sigma _\phi `$ have the same spatial variation as $`\sigma _z`$. Notice that such behaviour is found to be in agreement with the sparse existing data on the velocity dispersions in external galaxies (Gerssen et al. 1997, 2000). The only remaining quantities that need to be set are normalizations of the dispersions. We adopted the values $`(\sigma _{R0},\sigma _{\phi 0},\sigma _{z0})=(100,75,50)`$ km/s. At a galactocentric radius of 8 kpc, this gives velocity dispersions of (37,28,18) km/s, which roughly corresponds to the measured velocity dispersions in the solar neighbourhood.

The dust was chosen to be distributed in a similar way as the stellar distribution,

$$\kappa _\lambda (𝒓)=\frac{\tau _\lambda }{2z_\text{d}}\mathrm{exp}\left(\frac{R}{h_\text{d}}\right)\mathrm{exp}\left(\frac{|z|}{z_\text{d}}\right),$$
(28)
whereby $`\tau _\lambda `$ (here) represents the face-on optical depth through the centre of the galaxy,

$$\tau _\lambda =_{\mathrm{}}^{\mathrm{}}\kappa _\lambda (𝒓)\text{d}z.$$
(29)
For the dust scale length and scale height, we again choose the values from Byun et al. (1994), $`h_\text{d}=4`$ kpc and $`z_\text{d}=140`$ pc. Notice that the dust is hence significantly more confined to the plane of the galaxy, a feature which gives rise to the well-known dust lanes in edge-on galaxies. For the velocity field of the dust grains, we assumed an isotropic gaussian distribution,

$$F_\text{d}(𝒓,𝒗)\mathrm{exp}\left[\frac{v_R^2+(v_\phi \overline{v}_\phi )^2+v_z^2}{2\sigma _\text{d}^2}\right].$$
(30)
We assumed the same mean velocity as the stellar distribution, and a velocity dispersion $`\sigma _\text{d}=8`$ km/s, which is a typical value for the dispersion of the cold gas in the interstellar medium (Sofue & Rubin 2001). Finally, for the optical properties of the dust we took the average Galactic values (Gordon, Calzetti & Witt 1997).

These functions and parameters completely determine our model, except for the one free parameter $`\tau _\text{V}`$, the central face-on optical depth in the $`V`$ band (the corresponding optical depth at other wavelengths is determined by the assumed optical properties of the dust). As discussed in the Introduction, the optical depth in spiral galaxies is still a matter of debate. We ran our models for the optical depth values $`\tau _\text{V}=0`$, 0.5, 1, 2, 5 and 10, in order to cover a wide range of possible scenarios.

### 4.2 Modelling results

In this section we present and describe the results of our Monte Carlo simulations. We will not discuss the effects of dust on the photometry, since such results have been discussed at length by other authors, such as Byun et al. (1994), Corradi et al. (1996) and Bianchi et al. (1996). Instead we concentrate on the effects of dust on the observed kinematics, in particular on the mean projected velocity $`v_\text{p}(𝒙)`$ and the projected velocity dispersion $`\sigma _\text{p}(𝒙)`$.

In figures 3 and 4 we compare, at various inclination angles, the 2D mean projected velocity and projected velocity dispersion fields of an optically thin galaxy (left) with the corresponding fields of a dusty galaxy with optical depth $`\tau _\text{V}=1`$ (centre) and $`\tau _\text{V}=5`$ (right) respectively. Notice that the mean projected velocity field is not plotted in the face-on direction, as there is no rotation perpendicular to the galactic plane. To facilitate the interpretation of these 2D kinematical fields, we plotted in figure 5 the major and minor axis kinematical profiles of our models at various inclination angles, explicitly as a function of the optical depth.

#### 4.2.1 Edge-on galaxies

For edge-on galaxies, the effects of dust on both the mean projected velocity and the projected velocity dispersion are very strong, as is obvious from the bottom row panels of figure 3 and 4, and the panels in the right column of figure 5. The projected mean velocity profile along the major axis rises more and more slowly as the optical depth increases, and tends towards a solid body rotation. These trends were expected, and in agreement with the arguments of Davies (1990) and the absorption-only calculations of Bosma et al. (1992). To quantify this change of slope, we have calculated, for each optical depth, the major axis radius at which the mean projected velocity reaches half of its maximum value $`v_\text{p}^{\text{max}}=220\mathrm{sin}i`$ km s-1. These values are tabulated in the second column of table 1, and demonstrate that the effects of dust attenuation are very severe, even for small optical depths. The reason for this shallower slope of the rotation curve is that we mainly see the radiation from the optically thin part of the galaxy, which is more and more restricted to the outer part of the galaxy as the optical depth increases.

In the same way, we can understand the effects of dust attenuation on the projected velocity dispersion profile, where also a prominent dust lane is visible in the contour plot: we mainly see the light (and hence the Doppler shift) of the stars from the outer part of the disc, where the intrinsic velocity dispersion is much smaller.

#### 4.2.2 Galaxies at intermediate inclinations

Bosma et al. (1992) found that the signature of dust absorption on the rotation curve depends strongly on the inclination angle. For inclinations which deviate as few as 5 degrees from purely edge-on, they found that the effects of dust on the rotation curve is already severely reduced, even for optical depths up to $`\tau _\text{V}10`$. This result remains the same if scattering is also taken into account, because the absorption effects of dust strongly dominate in highly inclined galaxies. Table 1 quantitatively demonstrates these results. For the projected velocity dispersion, a similar effect is valid, at least for the major axis dispersion profile. When the entire projected velocity dispersion field (or just the minor axis dispersion profile) is taken into account, the asymmetrical signature of the dust can still be recognized for inclinations down to 80 degrees.

Moving to intermediate inclinations, the effects of dust attenuation on the observed kinematics quickly become very limited. For example, there is hardly any effect noticeable for an inclination of 60 degrees, which can easily be seen by comparing the corresponding panels in figure 3 and 4.

#### 4.2.3 Face-on galaxies

In the light of the previous results, we would expect that the effects of dust attenuation on the observed kinematics of face-on galaxies are completely neglegible. And indeed, on the mean projected velocity there is no effect at all: the projected mean velocity is zero for face-on galaxies, regardless of whether dust is present or not. However, we do see a conspicuous effect of dust attenuation on the projected velocity dispersion: it increases with increasing optical depth. This cannot be an absorption effect, because all stars along a given line of sight in our face-on galaxy discs have the same velocity dispersion $`\sigma _z`$ in the direction of the observer. The effect arises as a result of photons emitted in the central regions in a direction parallel to the plane of the galaxy. When these photons are scattered into a vertical direction, they have a large probability to leave the galaxy. They will then contribute the high-velocity information of the stars in the central regions to the LOSVDs, causing high-velocity wings and hence an increased velocity dispersion.

The effect is similar to what we found for the projected velocity dispersion profile of elliptical galaxies containing a diffuse dust component (Baes & Dejonghe 2002a). Compared to these models, where the signature of the dust is very significant, it is fairly modest for the face-on discs. The reason for this weaker impact is the fact that the velocity dispersion profile in our disc galaxy models have a fairly shallow slope, compared to the elliptical galaxy models of Baes & Dejonghe (2002a). In disc galaxies, the Doppler shifts contributed by the scattered photons will in the mean differ less from the Doppler shifts contributed by the direct photons as in elliptical galaxies.

We could ask ourselves the question whether this conclusion is biased by the chosen star-dust geometry of our disc galaxy models. Indeed, Baes & Dejonghe (2002a) found their results depended strongly on the relative star-dust geometry: the scattering effect was strong only for models in which the dust was more extended than the stars. In our disc galaxy models, the dust distribution has the same scale length as the stars, whereas the scale height is smaller, such that the dust is always sandwiched between the stars. Recently however, evidence has been found for extended distributions of cold dust in spiral galaxies, both from the radiative transfer modelling of the dust lanes in edge-on galaxies (Xilouris et al. 1997, 1998, 1999) and from ISO and SCUBA dust emission (Alton et al. 1999; Davies et al. 1999). In order to investigate whether this would bias our results, we considered a new set of models with an extended dust distribution. They were constructed by increasing the dust scale length or the dust scale height by factors up to three, while keeping the face-on optical depth fixed. Running new face-on simulations for these models, we found that the strength of the scattering effect on the projected velocity dispersion increases only minimally, with deviations of the order of a few per cent. The modest effects found in figure 5 are hence representative for disc galaxies.

## 5 Discussion

### 5.1 Stellar kinematics in disc galaxies

One of the most useful and promising applications of stellar kinematics in disc galaxies is their contribution to mass modelling. In its simplest form, this mass modelling consists of decomposing the observed rotation curve into contributions by the disc and the dark halo (assuming the contributions of a bulge and molecular/atomic gas are negligible). Even in this most simple scenario, there is still a virtually complete degeneracy: the mass-to-light ratio of the disc. Authors have often adopted a solution in which the contribution of the disc is maximized, leading to the so-called maximum disc hypothesis (van Albada & Sancisi 1986). There are however, no convincing arguments to favour a particular value for the stellar mass-to-light ratio, and often scenarios in which the dark matter halo dominates the disc are equally possible (e.g. Lake & Feinswog 1989). Stellar kinematics can, in principle, provide the necessary additional constraint to break this degeneracy. Indeed, for a self-gravitating disc, the vertical velocity dispersion scales linearly with the square root of the surface density of the disc, and can therefore be used to constrain the mass-to-light ratio of the disc.

The most obvious candidates for such an analysis would be face-on galaxies, because their projected velocity dispersion directly reflects the vertical velocity dispersion in the disc. We found that for a dusty disc galaxy, the observed velocity dispersion is affected by scattering. The effects are modest, however, even for models with an extended dust distribution. Alas, face-on galaxies do have strong observational disadvantages. Foremost, their low surface brightness make it very hard to obtain the stellar kinematics out to sufficiently large distances to isolate the disc component. Moreover, the shape and amplitude of the rotation curve cannot be measured directly. The amplitude can be estimated fairly accurate using the Tully-Fisher relation, whereas the shape could be approximated through the (controversial) universal rotation curve formalism (Persic, Salucci & Stel 1996). This is unsatisfactory, however, because there is a large spread in universal rotation curves, and the details of mass modelling usually depend crucially on small scale features of the rotation curve which statistical methods cannot incorporate.

Edge-on galaxies don’t suffer from these observational problems: their surface brightness is much higher, and their apparent rotation curve can be determined directly. However, the interpretation of the kinematical data is non-trivial. Firstly, the vertical velocity dispersion must be linked to the observed dispersion, which is a mixture of contributions from the radial and tangential velocity components. Based on the observed axis ratio of velocity ellipsoid in the solar neighborhood, Bottema (1993) estimated the vertical velocity dispersion in edge-on galaxies by means of the radial velocity dispersion evaluated at one scale length. It has been shown, however, that the velocity ellipsoid axis ratios in disc galaxies can vary substantially (Gerssen et al. 2000). Secondly, our modelling indicates that both the projected mean velocity and the projected velocity dispersion field in edge-on galaxies are strongly affected by interstellar dust.

Galaxies of an intermediate inclination ($`i<80^{}`$), however, are most suitable for such an analysis. Gerssen et al. (1997) have shown that for galaxies with an intermediate inclination, the three components of the velocity ellipsoid can, in principle, be determined by studying the kinematics along major and minor axes. Also, because the rotation curve of inclined galaxies can directly be determined and their surface brightnesses are high enough to allow the measurements of their kinematics in a reasonable time, they form ideal targets for such a study. Moreover, determining the three components of the velocity ellipsoid is very useful to constrain the dynamical history of disc galaxies, because the various proposed mechanisms leave a different imprint on the velocity ellipsoid. Our modelling demonstrates that the observed kinematics of inclined disc galaxies can be directly interpreted and are not biased by dust attenuation. We argue that studying the observed kinematics of a substantial set of optically smooth inclined galaxies could thus greatly contribute to understanding the mass structure and dynamical history of spiral galaxies.

### 5.2 Optical rotation curves in disc galaxies

The results of our modelling are important concerning the derivation of the rotation curve of spiral galaxies. Whereas the outer parts of the rotation curves of spirals were the preferred field of interest for nearly three decades, it is the inner slope of the rotation curve which nowadays stands in the spotlight. To measure the slope in this inner region accurately, a high spatial resolution is required, which makes the usual 21-cm measurements less suitable, except for the most nearby galaxies. A useful alternative are CO rotation curves, which can achieve a high resolution in both the spatial and velocity directions, but these observations are quite costly in observing time. The main alternative remains optical H$`\alpha `$ observations, either long-slit or Fabry-Perot.

Using H$`\alpha `$ rotation curves, several authors (Blais-Ouellette, Amram & Carignan 2001; Borriello & Salucci 2001; de Blok, McGaugh & Rubin 2001; de Blok & Bosma 2002) have recently found shallow slopes for the inner rotation curve of a significant number of (mainly LSB and dwarf) galaxies. Such slopes are not in agreement with the results of high-resolution cosmological N-body simulations, where it is found that dark matter halos have a strong cusp, and hence that rotation curves of dark matter dominated galaxies should be steeply rising (Navarro, Frenk & White 1997; Moore et al. 1997). A number of “solutions” have been proposed to explain this discrepancy, amongst them the suggestion that the internal extinction in the centre of galaxies could be a major source of uncertainty. Our models show that, indeed, an edge-on galaxy with an observed slowly rising rotation curve could in principle be a dusty galaxy with an intrinsically steep rotation curve. We found, however, that this effect cannot be obtained with realistic amounts of dust if the galaxies deviate more than a few degrees from exactly edge-on, confirming the results of Bosma et al. (1992) and Matthews & Wood (2001). In addition, because LSB galaxies are assumed to be dust poor (McGaugh 1994, Tully et al. 1998; Matthews et al. 1999), we have demonstrated that dust effects cannot explain the observed discrepancies.

We do however demonstrate that dust effects are important for galaxies within a few degrees from edge-on, and that dust attenuation must be seriously taken into account in such cases.

### 5.3 The importance of the dust grain velocities

The calculation of the final Doppler shift of a photon in the SKIRT code takes into account the velocity information of the star and the individual velocity vectors of each individual scattering dust grain. This means that for every scattering event, we need to sample a velocity from the dust velocity field, and this is a costly operator. It would be very convenient from a computational point of view, if we could neglect the dust grain velocities. We have done this in our previous Monte Carlo simulations (Baes & Dejonghe 2001b, 2002a), where we studied the observed stellar kinematics of elliptical galaxies. In these models, both stars and dust grains are supported by random motions, and because the dust is colder than the stars, the extra Doppler shifts due to scattering are generally smaller than the stellar Doppler shift. In a spiral galaxy, however, this argument does not hold anymore, because the velocities of both stars and dust grains are dominated by the rotation and are therefore of the same magnitude.

In order to investigate the importance of including the individual dust grain velocities into the code, we ran the same models again but without taking the dust velocities into account. In figure 6 we plot the major and minor axes mean projected velocity and projected velocity dispersion profiles of our models calculated this way. Comparing this figure with the corresponding this figure with figure 5, we find that the effects of dust on the mean projected velocity profiles are significantly overestimated when the dust grain velocities are not taken into account, in particular for inclined galaxies. The differences between the projected velocity dispersion curves in both figures are even worse: when the dust grain velocities are neglected, the apparent effect of dust attenuation is to increase the dispersion at nearly all radii and all inclinations, by factors of up to 200 per cent or more for realistic optical depths.

To understand this behaviour, consider a photon emitted by a star with velocity $`𝒗_{}`$ in the direction $`𝒌_0`$, and assume this photon is scattered by a dust grain with velocity $`𝒗_\text{d}`$ into the direction $`𝒌_1`$, whereafter it manages to escape the galaxy. When we do not take the dust grain velocity into account, the observed Doppler shift of the photon will simply be $`u_0=𝒗_{}𝒌_0`$. On the contrary, the correct observed Doppler shift of the photon reads

$`u_1`$ $`=𝒗_{}𝒌_0+𝒗_\text{d}(𝒌_1𝒌_0)`$ (31)
$`=(𝒗_{}𝒗_\text{d})𝒌_0+𝒗_\text{d}𝒌_1.`$ (32)
In a disc galaxy, the motion of both stars and dust grains is dominated by the same rotation. If the pathlength of the photon between its emission and scattering is short, the velocity vectors of the star and dust grain will therefore be very similar, $`𝒗_{}𝒗_\text{d}`$, such that $`u_1𝒗_{}𝒌_1`$. The new correct Doppler shift is hence a typical line-of-sight velocity of a star in the new propagation direction $`𝒌_1`$ of the photon. Therefore, taking the dust motion into account more or less forces the dust grain to adopt a Doppler shift that is “appropriate” for the new propagation direction. Along a given line-of-sight, the velocity information carried by scattered photons will, in the mean, only modestly deviate from the velocity information carried by photons directly emitted in same direction. If we don’t take the dust grain velocities into account, the photon can contribute, to the observed kinematics in a given direction, a line-of-sight velocity typical for a completely different direction.

In the elliptical galaxy models from Baes & Dejonghe (2002a), this does not make much difference, because the velocities of the stars have similar magnitudes in the different directions. In a rotating disc galaxy however, this makes a huge difference because the stars have hugely different velocities in different directions: the vertical and radial motions are determined by the dispersion and are of the order of a few tens of km/s, whereas the azimuthal motion is dominated by the rotation and is of the order of 200 km/s. If the dust grain velocities are not correctly taken into account, scattered photons can thus cause erroneous extremely-high-velocity wings in the LOSVDs. These give rise to modest but significant errors in the mean projected velocity, and very large errors in the projected velocity dispersion.

As a remark, we want to point out that Matthews & Wood (2001) found that “also including the Doppler shifts arising from the relative bulk motion of the scattering dust particles had a negligible effect on the final rotation curves”. Whereas this statement is indeed true for the mean projected velocity of highly inclined galaxies, we nevertheless judge that this could easily be misinterpreted as a justification to neglect dust velocities altogether. It must be realized that for the calculation of the entire LOSVDs, and for the projected velocity dispersion in particular, the dust velocities do play a crucial role, and must be taken into account.

## 6 Conclusions

In this paper we have presented a novel Monte Carlo code that can take velocity information into account, and can hence be used for kinematical modelling of dusty objects. We applied this code to calculate the observed kinematics of dusty disc galaxies. The main results of this paper can be summarized as follows.

1. A correct inclusion of kinematical information into radiative transfer problems requires the inclusion of the velocities of both the emitting stars and the scattering dust grains. We see no other way of tackling this problem except with Monte Carlo techniques.
2. A new approach is presented to optimize the integration through the dust in a Monte Carlo code, using a trilinear interpolation instead of a constant opacity within each cell. We compared both kinds of grids, and find that, to obtain a similar accuracy, the new approach is more efficient in terms of computation time and memory.
3. The effects of dust attenuation on the kinematics of edge-on disc galaxies are severe. Both the mean projected velocity and the projected velocity dispersion are severely affected, even for modest optical depths. Therefore we strongly advise to always interprete the stellar kinematics and optical rotation curves of edge-on galaxies very cautiously.
4. For galaxies which are more than a few degrees from strictly edge-on, the effects of interstellar dust on the observed kinematics are much weaker. Therefore, we argue that dust attenuation cannot be invoked as a possible mechanism to reconcile the observed slope of LSB galaxies with the predicted CDM cosmological models.
5. Dust attenuation does not affect the kinematics of intermediately inclined galaxies. Such galaxies hence form the ideal targets for stellar kinematical studies, in particular to constrain the mass structure and to study the kinematical history of disc galaxies.
6. The projected velocity dispersion of face-on galaxies increases slightly due to scattering of dust grains into the line of sight. The effects are relatively small, however, even for extended dust distributions.
7. Neglecting the extra Doppler effect caused by the scattering medium results in incorrect projected kinematics, in particular if the velocity components of the stars (and dust grains) in the galaxy differ considerably, such as in a rotating disc. If the dust grain velocities are neglected, the mean projected velocity is significantly underestimated, and the projected velocity dispersion field is completely overestimated.

## Acknowledgements

M. B. gratefully thanks the Fund for Scientific Research for financial support.

## Appendix A Generating positions from an exponential disc galaxy

In this Appendix, we show how the SKIRT code generates a random position from the exponential disc model, i.e. how a random position is drawn from the three-dimensional probability density

$$p(𝒓)\text{d}𝒓\mathrm{exp}\left(\frac{R}{h_{}}\right)\mathrm{exp}\left(\frac{|z|}{z_{}}\right)\text{d}𝒓.$$
(33)
As this probability density is a separable function of $`R`$, $`\phi `$ and $`z`$, we can generate a random position by independently generating each of these three coordinates, whereby the generation of an azimuth is trivial. To generate a random radius and height, we use the transformation technique, which says that a random variable $`x`$ can be drawn from a probability density $`p(x)\text{d}x`$ by generating a uniform deviate $`X`$ and solving the equation

$$X=_{\mathrm{}}^xp(x^{})\text{d}x^{}/_{\mathrm{}}^{\mathrm{}}p(x^{})\text{d}x^{}$$
(34)
for $`x`$. The calculation of a random $`z`$ is easy using this principle and results in

$$z=z_{}\mathrm{sgn}(1+2X)\mathrm{ln}(1|12X|).$$
(35)
This procedure works equally well when we introduce a cut-off, or when a sech or isothermal profile is used instead of the exponential profile to describe the vertical distribution of the stars (e.g. Bianchi et al. 1996).

Applying the same technique to generate a random $`R`$, we can find $`R`$ by solving the equation

$$X=1\left(1+\frac{R}{h_{}}\right)\mathrm{exp}\left(\frac{R}{h_{}}\right).$$
(36)
Bianchi et al. (1996) invert this transcendental equation numerically to find $`R`$. However, this equation can be solved exactly by means of the Lambert function, yielding

$$R=h_{}\left[1W_1\left(\frac{X1}{\text{e}}\right)\right].$$
(37)
The Lambert function, also known as the product log function, is generally defined as the inverse of the function $`wf(w)=w\text{e}^w`$, and $`W_1(z)`$ represents the only real branch of this complex function besides the principle branch. The Lambert function can be computed in a very efficient way by means of Halley iteration; for more information on this function see Corless et al. (1996).


format: text

text:
# Untitled Document

Binding energies of hydrogen-like impurities in a

semiconductor in intense terahertz laser fields

J.L. Nie1,2), W. Xu1),∗, and L.B. Lin2)

1) Department of Theoretical Physics

Research School of Physical Sciences and Engineering

Australian National University, Canberra ACT 0200, Australia

2) Department of Physics, Sichuan University

Chengdu - 610064, The People’s Republic of China

PACS: 71.55.Eq, 78.50.Ge, 78.66.Fd

We present a detailed theoretical study of the influence of linearly polarised intense terahertz (THz) laser radiation on energy states of hydrogen-like impurities in semiconductors. The dependence of the binding energy for ground ($`1s`$) and first excited ($`2s`$) states, $`E_{1s}`$ and $`E_{2s}`$, on intensity and frequency of the THz radiation has been examined for a GaAs-based system. It is found that $`E_{1s}`$, $`E_{2s}`$ and $`E_{2s}E_{1s}`$ decrease with increasing radiation intensity or with decreasing radiation frequency, which implies that an intense THz field can enhance ionisation of dopants in semiconductors. Our analytical and numerical results show that one of the most important results obtained by A.L.A. Fonseca et al \[phys. stat. sol. (b) 186, K57 (1994)\] is incorrect.

1. Introduction

The ionisation of dopants (i.e., donor and acceptor impurities) provides the major source of carriers (i.e., electrons and holes) in semiconductors. The ionised dopants are also the major sources of electron-impurity scattering in semiconductor devices, which determine the transport and optical properties of the device systems at low-temperatures. Hence, the investigation of ionisation of dopants (such as binding energies of donors and acceptors, transition energies and probabilities among different impurity states, etc.) is fundamental in understanding almost all physically measurable properties in semiconductors. In the absence of an intense electro-magnetic (EM) radiation field, the binding energies of dopants and the transition energies among different impurity states in popularly used semiconductor materials are known  and the theoretical approaches to calculate these impurity states are well-documented .

It should be noted that in semiconductor materials such as GaAs, Ge and Si, the binding energies of donor and acceptor impurities are of the order of terahertz ($`10^{12}`$ Hz or THz) photon energies  so that an intense THz radiation can affect strongly the impurity states. With development and application of coherent, high-power, long-wavelength, frequency-tunable and linearly polarised radiation sources such as THz or far-infrared (FIR) free-electron lasers (FELs) , it has now become possible to measure the effect of an intense laser radiation on ionisation and perturbation of dopants (especially shallow impurities) in different semiconductor systems. In recent years, using THz FELs \[such as Free Electron Laser for Infrared eXperiments (FELIX) in The Netherlands and CW FELs at UCSB\] as intense radiation sources, THz-photon-induced impact ionisation in InAs heterostructures , time-resolved shallow donor spectrum in Si-doped GaAs  and Lyman transitions in Be-doped GaAs  have been investigated through, e.g., transport and/or photoconduction measurements. The results obtained experimentally indicate that in the presence of intense laser radiation such as FEL fields, i) impurity states in different semiconductor systems are perturbated by the intensity and frequency of the THz laser fields \[5-7\]; ii) photoconduction experiments are more sensitive than optical measurements for detecting transition energies of impurity states in semiconductors ; and iii) some interesting intense radiation phenomena, such as impact ionisation of dopants  and splitting and broadening of the impurity spectrum , can be observed. In order to understand these fundamentally new experimental findings and to predicate new radiation phenomena, it is essential to know theoretically how an intense laser field affects the binding energies of impurities in semiconductors and is the prime motivation of the present study.

In this paper, we study binding energy of hydrogen-like impurities in bulk semiconductors in the presence of intense THz laser fields. For semiconductor materials such as GaAs, there is an absolute conduction band minimum at $`𝐊=0`$ with $`\mathrm{\Gamma }_6`$ symmetry and the electron-effective-mass and dielectric constant can be considered to be isotropic. Therefore, the impurity states in materials such as GaAs are strictly scaled-down versions of those in the hydrogen atom. The small electron-effective-mass coupled with the large static dielectric constant result in small binding energy for donors in these materials . It has been noticed by atomic physicists in their early work that intense UV radiation can affect strongly the structure of atomic hydrogen . Using theoretical approaches proposed and developed in Ref. , Nunes and co-workers have investigated the influence of intense laser radiation on hydrogen-like impurity states in bulk semiconductors  and in semiconductor quantum well structures . However, in Ref. , only the ground-state energy at high-intensity radiation was presented analytically. We find that it is hard to understand the result given by Eq. (13) in Ref. , because it shows that at high-intensity radiation limit (ie, $`\alpha _0\mathrm{}`$) the ground-state binding energy $`E_0R_y^{}`$ (binding energy in the absence of the radiation) with $`R_y^{}`$ being the effective Rydberg constant. In order to see how impurity states in a semiconductor are perturbated within the frequency and intensity range of the current generation of the THz FELs, we think it is necessary to re-examine the dependence of the impurity binding energy in a semiconductor on frequency and intensity of the THz fields. Moreover, since most of the experimental work in this area has been carried out using photoconduction experiments  which measure the transition energies among different impurity states, it is also necessary to know how an intense THz radiation affects the excited impurity states in a semiconductor. In Section 2, we will first briefly introduce the theoretical approaches used for the calculations then present analytical results of the binding energies for ground and first excited impurity states. The numerical results will be presented and discussed in Section 3 and the conclusions drew from this study will be summarised in Section 4.

2. Analytical results

In this work, we consider the situation where a laser field with a vector potential $`A(t)`$ is applied along the x-direction of a semiconductor and the radiation field is polarised linearly along the z-direction. Under the effective mass approximation, the time-dependent Hamiltonian to describe the electron-impurity system is given as

$$H(t)=\frac{1}{2m^{}}[p_x^2+p_y^2+(p_zeA(t))^2]\frac{e^2}{\kappa R}.$$
$`(1)`$
Here, $`m^{}`$ is the electron-effective-mass, $`p_x=i\mathrm{}/x`$ is the momentum-operator along the x-direction, under the dipole approximation $`A(t)=(F_0/\omega )\mathrm{sin}(\omega t)`$ is the vector potential induced by the radiation field with $`\omega `$ being the radiation frequency and $`F_0`$ the electric field strength of the radiation field, $`𝐑=(𝐫,z)=(x,y,z)`$, and $`e^2/\kappa R`$ is the Couloumb potential induced by electron-impurity interaction with $`\kappa `$ being the dielectric constant. It can be seen that due to time-dependent nature of the radiation field, in principle, we have to solve time-dependent Schrödinger equation: $`[i\mathrm{}/tH(t)]\mathrm{\Psi }(𝐑,t)=0`$ to obtain the electronically bounded impurity states. However, at present, there is no simple analytical solution to this problem . Here we employ a tractable theoretical approach proposed and developed by atomic physicists  to study the binding energy of the impurity states. The key issue of this approach is to find laser-dressed potential energy. It is found that if the coordinate along the z-axis is shifted by $`z^{}=z+z_0(t)`$ with $`z_0(t)=(eF_0/m^{}\omega ^2)\mathrm{cos}(\omega t)`$, the Schrödinger equation becomes

$$i\mathrm{}\frac{\mathrm{\Psi }(𝐑^{},t)}{t}=\left[\frac{p_x^2+p_y^2+p_z^{}^2}{2m^{}}+V(𝐑^{},t)\right]\mathrm{\Psi }(𝐑^{},t),$$
$`(2a)`$
where $`𝐑^{}=(𝐫,z^{})`$ and

$$V(𝐑^{},t)=\frac{e^2}{\kappa \sqrt{r^2+[z^{}z_0(t)]^2}}+\frac{[eA(t)]^2}{2m^{}}.$$
$`(2b)`$
After averaging $`V(𝐑^{},t)`$ over a period of the radiation field, we obtain laser-dressed potential as

$$V(𝐑^{})=\frac{\omega }{2\pi }_0^{2\pi /\omega }𝑑tV(𝐑^{},t)=E_{em}\frac{e^2}{2\kappa }[|𝐑^{}+𝐑_0|^1+|𝐑^{}𝐑_0|^1].$$
$`(3)`$
Here $`E_{em}=(eF_0)^2/4m^{}\omega ^2`$ is the energy shift induced by the radiation field due to dynamical Franz-Keldysh effect  and $`𝐑_0=(0,0,eF_0/m^{}\omega ^2)`$. It should be noted that in atomic physics, $`E_{em}`$ is also called ponderomotive energy . If we replace $`V(𝐑^{},t)`$ in Eq. (2a) by $`V(𝐑^{})`$ given by Eq. (3), the time-dependent problem becomes approximately a time-independent one and the solution is $`\mathrm{\Psi }(𝐑^{},t)=e^{i(E_N+E_{em})t/\mathrm{}}\mathrm{\Psi }_N(𝐑^{})`$, where the wavefunction $`\mathrm{\Psi }_N(𝐑^{})`$ and the energy spectrum $`E_N`$ are determined by $`(H_0E_N)\mathrm{\Psi }_N(𝐑^{})=0`$ with

$$H_0=\frac{p_x^2+p_y^2+p_z^{}^2}{2m^{}}\frac{e^2}{2\kappa }(|𝐑^{}+𝐑_0|^1+|𝐑^{}𝐑_0|^1).$$
$`(4)`$
In contrast to Refs. \[8-10\], in the derivations shown above we have included the terms associated with $`F_0^2`$ which, as can be seen, only contributes to the energy shift induced by the radiation field alone. The main merit of this approach is that now one can convert a time-dependent problem into a time-independent one and, thus, one can use, e.g., variational method to calculate the energy states and the corresponding wavefunctions. It should be noted that the present theoretical approach is a generalization of those proposed by the atomic physicists. As has been pointed out by Ref. , this approach is the lowest order approximation to the atomic problem. Furthermore, the results given by Eqs. (3) and (4) have been documented by Refs.  for the case of the hydrogen-like impurities in a semiconductor.

In the present work, we employ the wavefunctions for a hydrogen atom as the trial functions for ground state (ie, $`1s`$-state) and first excited state (ie, $`2s`$-state) of the donor impurities, which reads

$$\mathrm{\Psi }_{1s}(𝐑^{})=ae^{\beta R^{}/a_0^{}}$$
$`(5a)`$
and

$$\mathrm{\Psi }_{2s}(𝐑^{})=c(1+bR^{}/a_0^{})e^{\gamma R^{}/2a_0^{}},$$
$`(5b)`$
where $`\beta `$ and $`\gamma `$ are variational parameters, $`a_0^{}=\mathrm{}^2\kappa /m^{}e^2`$ is the effective Bohr radius, $`a=(\beta ^3/\pi a_{0}^{}{}_{}{}^{3})^{1/2}`$ is determined by $`<1s|1s>=1`$, and $`b=(\beta +\gamma /2)/3`$ and $`c=\gamma ^{5/2}/[8\pi a_{0}^{}{}_{}{}^{3}(12b^2+6b\gamma +\gamma ^2)]^{1/2}`$ are determined respectively by $`<1s|2s>=0`$ and $`<2s|2s>=1`$. After introducing these trial functions to the variational method, i.e., $`E_N=<N|H_0|N>`$, the binding energies of the $`1s`$-state and $`2s`$-state are obtained as

$$E_{1s}=R_y^{}\left[\frac{2}{\alpha _0^{}}\beta ^2\frac{2}{\alpha _0^{}}(1+\alpha _0^{}\beta )e^{2\alpha _0^{}\beta }\right],$$
$`(6a)`$
and

$$E_{2s}=\frac{R_y^{}}{4}\frac{(4/\alpha _0^{})(24b^2+12b\gamma +2\gamma ^2𝒜e^{\alpha _0^{}\gamma })\gamma ^2(4b^2+2b\gamma +\gamma ^2)}{12b^2+6b\gamma +\gamma ^2},$$
$`(6b)`$
where $`\alpha _0^{}=eF_0/m^{}\omega ^2a_0^{}`$, $`R_y^{}=m^{}e^4/2\mathrm{}^2\kappa ^2`$ is the effective Rydberg constant, and

$$𝒜=24b^2+12b\gamma +2\gamma ^2+\alpha _0^{}\gamma (18b^2+8b\gamma +\gamma ^2)+2\alpha _{0}^{}{}_{}{}^{2}b\gamma ^2(3b+\gamma )+\alpha _{0}^{}{}_{}{}^{3}b^2\gamma ^3.$$

Furthermore, the variational parameters $`\beta `$ and $`\gamma `$ are determined respectively by

$$\beta (1+2\alpha _0^{}\beta )e^{2\alpha _0^{}\beta }=0\mathrm{and}E_{2s}/\gamma =0.$$

The results shown above indicate that the effect of the radiation field on binding energy of the impurities is mainly achieved through a dimensionless factor $`\alpha _0^{}=e^3F_0/\mathrm{}^2\kappa \omega ^2F_0/\omega ^2`$. For high-frequency and/or low-intensity radiation so that $`\alpha _0^{}1`$, we have

$$E_{1s}R_y^{}\left(1\frac{4}{3}\alpha _{0}^{}{}_{}{}^{2}\right)\mathrm{and}E_{2s}\frac{R_y^{}}{4}\left(1\frac{2}{3}\alpha _{0}^{}{}_{}{}^{2}\right),$$
$`(7)`$
which suggests that the binding energy decreases quickly with increasing radiation intensity or with decreasing radiation frequency when $`\alpha _0^{}1`$. When $`\alpha _0^{}=0`$, $`E_{1s}=R_y^{}`$ and $`E_{2s}=R_y^{}/4`$ are well-known results obtained in the absence of the radiation field. For low-frequency and/or high-intensity radiation, entailing $`\alpha _0^{}1`$, we find that there is no simple analytical expression for both $`E_{1s}`$ and $`E_{2s}`$, in contrast to Eq. (13) obtained by Ref. . When $`\alpha _0^{}\mathrm{}`$, we find $`\beta 0`$ and $`\gamma 0`$ so that $`E_{1s}0`$ and $`E_{2s}0`$ (note that a zero binding energy implies that all impurities at this state are ionised). These theoretical results indicate that $`E_{1s}`$ and $`E_{2s}`$ are altered from respectively $`R_y^{}`$ and $`R_y^{}/4`$ at zero field to zero at high-field limit. One of the most important results obtained by Ref.  is that at high-intensity radiation fields, the ground-state binding energy is given by Eq. (13) which shows $`E_{1s}R_y^{}`$ when $`\alpha _0\mathrm{}`$ and $`E_{1s}<R_y^{}`$. Our results here suggest that Eq. (13) obtained by Ref.  is incorrect.

3. Numerical results

Below we present numerical results for semiconductor materials such as GaAs. The material parameters for GaAs taken within the calculations are the effective-electron-mass ratio $`m^{}/m_e=0.0665`$ with $`m_e`$ being the electron rest mass and the static dielectric constant $`\kappa =12.9`$. The dependence of binding energies $`E_{1s}`$ and $`E_{2s}`$ as well as transition energy $`E_{2s}E_{1s}`$ on THz laser radiation fields is shown in Figures 1 - 4. From these results, we see that:

a) in the presence of the radiation fields, $`E_{1s}`$ and $`E_{2s}`$ are altered from respectively $`R_y^{}`$ and $`R_y^{}/4`$ at low-field limit to zero at high-field limit and $`E_{2s}E_{1s}`$ is altered from $`3R_y^{}/4`$ at low-field limit to zero at high-field limit. This confirms that Eq. (13) obtained by Ref.  is incorrect;

b) with increasing radiation intensity and/or decreasing radiation frequency, $`E_{1s}`$, $`E_{2s}`$ and $`E_{2s}E_{1s}`$ decrease;

c) $`E_{2s}`$ depends a bit weakly on the radiation field than $`E_{1s}`$ does, which can be understood by Eq. (7);

d) the strong effect of the radiation field on binding energy and transition energy can be observed at $`\alpha _0^{}1`$; and

e) $`E_{1s}`$, $`E_{2s}`$ and $`E_{2s}E_{1s}`$ depend more strongly on radiation frequency than on radiation intensity because $`\alpha _0^{}F_0/\omega ^2`$.

The results discussed above indicate that in the presence of the intense THz laser fields, the binding energies of the impurity states can be reduced and the impurity spectrum in semiconductors can be shifted significantly by the radiation.

It should be noted that the current generation of the FELs can provide intense THz radiation sources in the frequency and intensity range $`f0.110`$ THz and $`F_00.1100`$ kV/cm  so that the condition $`\alpha _0^{}=e^3F_0/\mathrm{}^2\kappa \omega ^21`$ can be satisfied by most of the popularly used semiconductor materials. We therefore believe that it has now become possible to investigate the effects of the intense laser radiation on ionisation and perturbation of dopants in semiconductor systems by using current generation of THz or FIR FELs.

4. Further remarks

In the present study, our calculations have been performed for shallow-impurities in semiconductor materials such as GaAs, because these impurities are hydrogen-like. It is well known that in semiconductors, deep-impurities may not be hydrogen-like because their states depend strongly on lattice structure of the material systems, such as the symmetry of the crystal potential, electronic band structure, etc. To examine the effect of intense laser radiation on deep-impurity states in a semiconductor, the above mentioned theoretical approach has to be modified.

To our knowledge, except those theoretical results reported by Refs.  and  for semiconductor systems, so far very few theoretical and experimental results have been reported regarding the direct effects of the intense laser radiation on binding energies of electrons in atomic systems (two recent reviews are given in Refs.  and ) or of hydrogen-like impurities in semiconductor systems. The main reason behind this is that when an electronic system (atoms or semiconductors) is subjected to intense laser fields, the investigation of the electronic states (or electronically bounded impurity states) is essentially a time-dependent problem. As a consequence, these states are no longer the eigenstates and transitions among these time-dependent states can occur. Therefore, it is very hard to measure experimentally the binding energy for a certain atomic or impurity state when an intense laser radiation is present. Normally, the results obtained from atomic physics experiments  are contributions from all possible states and transitions among these states.

The analytical and numerical results shown in this paper are what one can get from the simplest theoretical model. As has been pointed out by Ref. , the approach used in the present study is the lowest order approximation to deal with the problem. Under such an approximation, the systems can be described in terms of a set of quasi-eigenstates of atom (or impurity) and laser field and there is no transition among these quasi-eigenstates. Namely, under this approximation, the system is “stable to ionisation”. Thus, the results obtained from this model can be useful in understanding atomic physics phenomena such as free-free transitions and dichotomy of the atomic states . Furthermore, it can be shown  that the conditions under which the approach can be applied are $`\mathrm{}\omega E_{em}`$ and $`R_0^2\omega 1`$. Therefore, the very low-frequency and high-intensity results shown in Figs. 1-4 in this paper may not be the case. To examine the effect of these extreme radiation conditions on states of hydrogen-like impurities in a semiconductor, the contributions from higher order terms (see Eq. (28) in Ref. ) have to be included. The inclusion of these higher order contributions requires considerably further analytical and numerical calculations and, therefore, we do not attempt it in the present study.

5. Conclusions

In this paper, we have re-examined how a linearly polarised intense laser field affects the binding and transition energies of shallow-donor impurities in a bulk semiconductor. The present study has been conducted on the basis of a theoretical approach proposed by Ref.  and documented by Refs. . We have obtained results of binding energies for both ground and first excited impurity states and thus we can look into the effects of the radiation field on transition energy between $`1s`$ and $`2s`$ states, in conjunction with experimental work carried out by photoconductivity measurements . The numerical and analytical results obtained from this work show that Eq. (13) in Ref.  is incorrect.

We have found that the strong influence of the radiation field on binding and transition energies of the shallow impurity states in bulk semiconductors can be observed when the condition $`\alpha _0^{}=e^3F_0/\mathrm{}^2\kappa \omega ^21`$ is satisfied. For most popularly used semiconductor materials such as GaAs, Si and Ge, these radiation conditions have been realised by the current generation of the THz or FIR FELs. We therefore hope that the phenomena studied and predicated in this paper can be verified experimentally.

Acknowledgment

One of us (W.X.) is a Research Fellow of the Australian Research Council (ARC). This work was also supported by the ARC IREX Grant, Visiting Scholar Foundation of Key Laboratory in University of the Ministry of Education, China and National Natural Science Foundation of China. Discussions with P. Fisher and R.A. Lewis (UoW, Australia) are gratefully acknowledged.

References

Corresponding author and electronic mail: wen105@rsphysse.anu.edu.au

 For a review, see, e.g., A.K. Ramdas and S. Rodriguez, Rep. Prog. Phys. 44, 1297 (1981).

 For a early work, see, e.g., J.M. Luttinger and W. Kohn, Phys. Rev. 97, 869 (1955).

 For example, in the absence of a radiation field, the donor binding energy in n-GaAs is 5.4 meV (or 1.3 THz) and for P in Si is 45.6 meV (or 11 THz).

 For recent development of the FELs, see, eg, Free Electron Lasers 1999, edited by J. Feldhaus and H. Weise (Elsevier Science, North Holland, 2000).

 A.G. Markelz, N.G. Asmar, B. Brar and E.G. Gwinn, Appl. Phys. Lett. 69, 3975 (1996).

 P.C.M. Planken, P.C. van Son, J.N. Hovenier, T.O. Klaassen, W. Th. Wenckebach, B.N. Murdin and G.M.H. Knippels, Phys. Rev. B 51, 9643 (1995); Infrared Phys. Technol. 36, 333 (1995).

 R.A. Lewis, L.V. Bradley, and M. Henini, Solid State Commun. 122, 223 (2002).

 M. Gavrila and J.Z. Kaminski, Phys. Rev. Lett. 52, 613 (1984); M. Point, N.R. Walet, M. Gavrila, and J. McCurdy, Phys. Rev. Lett. 61, 939 (1988).

 A.L.A. Fonseca, M.A. Amato, and O.A.C. Nunes, Phys. Stat. Sol. (b) 186, K57 (1994).

 Q. Fanyao, A.L.A. Fonseca, and O.A.C. Nunes, Phys. Stat. Sol. (b) 197, 349 (1996); Superlattices and Microstructures 23, 1005 (1998).

 A.D. Bandrauk and H.Z. Lu, Phys. Rev. A 62, 53406 (2002).

 W. Xu, Phys. Rev. B 57, 15282 (1998).

 See, e.g., K. Burnett, V.C. Reed, and P.L. Knight, J. Phys. B: At. Mol. Opt. Phys. 26, 561 (1993).

 The connection between the electric field strength of a laser field ($`F_0`$) and the laser output power ($`I`$) in vacuum is: $`I=0.5\sqrt{ϵ/\mu }|F_0|^21.32|F_0|^2`$ kW cm-2 where $`F_0`$ is in units of kV/cm.

 For a recent review, see, e.g., M. Protopapas, C.H. Keitel, and P.L. Knight, Rep. Prog. Phys. 60, 389 (1997).

Figure captions

Fig. 1: Binding energies, $`E_{1s}`$ and $`E_{2s}`$, as a function of THz radiation frequency ($`f=\omega /2\pi `$) for different radiation intensities. $`R_y^{}`$ is the effective Rydberg constant and for GaAs $`R_y^{}=5.44`$ meV.

Fig. 2: Transition energy between $`2s`$ and $`1s`$ impurity states, $`E_{2s}E_{1s}`$, as a function of radiation frequency for different radiation intensities.

Fig. 3: $`E_{1s}`$ and $`E_{2s}`$ as a function of radiation intensity for different radiation frequencies.

Fig. 4: $`E_{2s}E_{1s}`$ as a function of radiation intensity for different radiation frequencies.

J.L. Nie et. al. Figure 1.

J.L. Nie et. al. Figure 2.

J.L. Nie et. al. Figure 3.

J.L. Nie et. al. Figure 4.