EDP Sciences
Free Access
Issue
A&A
Volume 549, January 2013
Article Number A126
Number of page(s) 16
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201220468
Published online 10 January 2013

© ESO, 2013

1. Introduction

Essentially all objects in the Universe have a complicated, three-dimensional, dynamic structure. Understandably, throughout most of the history of astronomy the observed objects had to be modeled using significant restrictions and approximations. Among them, a set of assumptions about the global geometry of an object always played a pivotal role. For instance, stellar atmospheres were typically treated assuming a plane-parallel, horizontally-homogeneous geometry, which simplified the problem to one spatial dimension. Similarly, stellar winds, novae, planetary nebulae, and other extended sources, were modeled using an assumption of spherical symmetry, which again renders it a one-dimensional problem.

In the last two decades, there was much activity devoted to extend the traditional modeling techniques to treat structures that are truly three-dimensional (3D). In fact, detailed 3D hydrodynamic simulations have now become essentially routine. A non-exhaustive list of 3D astrophysical (radiation)(magneto)hydrodynamics (RHD, MHD or RMHD) codes includes ZEUS (Stone & Norman 1992a,b), ENZO (Bryan & Norman 1997; Norman & Bryan 1999; O’Shea et al. 2004; Norman et al. 2007), RAMSES (Teyssier 2002), HERACLES (González et al. 2007), PLUTO (Mignone et al. 2007), ATHENA (Stone et al. 2008) and its radiation module (Davis et al. 2012), AstroBEAR (Cunningham et al. 2009).

However, since most of the information about an astronomical object is conveyed to a distant observer by its radiation, detailed 3D hydrodynamic models have to be accompanied by adequate 3D radiation transfer solutions that provide a spectroscopic diagnostic information. Compared to 3D hydrodynamic models, the 3D radiation transfer solvers are much more complicated and difficult to treat numerically, because (i) there are many more quantities that describe a radiation field, as compared to the hydrodynamics, due to the directional and spectral dependence of radiation; (ii) a long-range interaction between the radiation field and the plasma arises because of a typically much larger mean free path of a photon compared to a mean free path of massive particles.

Although providing exact time-dependent, non-LTE radiation (magneto)hydrodynamic models of astronomical objects is generally viewed as a mighty goal, such a goal has not been fully achieved yet. However, a large progress was accomplished in recent years, specifically by 3D RMHD codes for stellar atmospheres, such as: the Copenhagen-Oslo Stagger Codes (Nordlund & Galsgaard 1995; Galsgaard & Nordlund 1996; Hansteen 2004; Hansteen et al. 2007), CO5BOLD (Freytag et al. 2002), MURaM (Vögler et al. 2005), Bifrost (Gudiksen et al. 2011). All of these codes solve the 3D radiative transfer by approximating the monochromatic opacities with a small number of mean opacities (the opacity binning scheme of Nordlund 1982). With typically four bins, they can model the frequency-integrated radiative energy losses and gains quite well.

There are also several recently developed computer codes specifically designed to provide 3D radiative transfer solvers, typically for a post-processing of 3D hydrodynamic simulations. The most widely used among them are the following: MUGA (Trujillo Bueno & Fabiani Bendicho 1995; Fabiani Bendicho et al. 1997; Fabiani Bendicho & Trujillo Bueno 1999), MULTI3D (Botnen 1997; Leenaarts & Carlsson 2009; Leenaarts et al. 2009), which is a generalization of the 1D code MULTI (Carlsson 1986), the RH code (Uitenbroek 2001, 2006), a 3D version of the atmospheric code PHOENIX (Hauschildt & Baron 2006, 2008; Baron & Hauschildt 2007), ASSeT (Koesterke et al. 2008; Koesterke 2009), SCATE (Hayek et al. 2011). They differ in their intended range of applications, and in many details of the numerical techniques. They all use some variant of the short-characteristics scheme (see Sect. 3), but differ in details of what is assumed about the behavior of state parameters between the grid points, and how the necessary interpolations are being performed. For details, the reader is referred to the above mentioned papers, and several excellent review papers (Carlsson 2008, 2009).

In this paper, we describe our variant of the three-dimensional radiative transfer solver, named IRIS. Analogously to most of the above mentioned techniques, it uses the short-characteristics scheme. Our solver differs from the previous ones in several respects: it is formulated only in the Cartesian coordinate system, but within this system it is formulated in the most universal way that allows for several variants of the spatial integration and interpolation to be tested and used depending on the actual application; it carefully treats subgriding to allow for line radiation transport in the presence of arbitrary (non-relativistic) velocity fields. Although the current version assumes local thermodynamic equilibrium (LTE) and no scattering, the code is prepared to be relatively straightforwardly extended to treat scattering and departures from LTE, the so-called non-LTE (or NLTE) situations. This will be reported in future papers.

The paper is organized as follows. In Sect. 2 we outline the main physical assumptions of the model solved by IRIS. In Sect. 3 we summarize the short-characteristics method, and its application in a 3D Cartesian grid. In Sect. 4 we describe the piecewise, cubic, locally monotonic, interpolation technique that we use. Then, in Sect. 5 we describe our implementation in IRIS of the techniques presented in Sects. 3 and 4. We also explain the procedure to handle general macroscopic velocity fields. The calculation of the moments of the radiation field and the related angular integration methods are the subjects of Sect. 6. We detail in Sect. 7 the method that we employ to treat media with periodic boundary conditions. Section 8 presents the results of relevant tests. We conclude our paper in Sect. 9 with a summary of the features of IRIS, and with an outline of possible extensions to the code.

2. Radiative transfer: basic definitions and assumptions

The general unpolarized radiative transfer equation (RTE) in the observer’s frame reads (Chandrasekhar 1950; Mihalas 1978; Mihalas & Mihalas 1984) (1)where I(r,n,ν,t) is the specific intensity of radiation at position r, propagating in the direction specified by the unit vector n, at frequency ν, and time t; χ(r,n,ν,t) is the absorption coefficient, η(r,n,ν,t) is the emission coefficient, and c is the speed of light. Although the above quantities depend on time, the current version of IRIS solves the time-independent RTE, for a given hydrodynamics structure at a given instant, i.e., for a given snapshot. In other words, we consider regimes in which the photon free-flight time is small compared to the fluid flow dynamical timescales, so that the radiation field gets fully stabilized before any change occurs in the flow dynamical properties. We assume that the typical flow velocities of the hydrodynamic structures are non relativistic.

As is customary, we introduce the source function, (2)To simplify the notations, we drop the explicit dependence of radiative quantities on the position, direction and time, and denote the dependence on frequency with a subscript ν. The transfer equation, Eq. (1), becomes (3)where s is the path length along the ray in the direction of propagation of the radiation.

The choice of a numerical method depends on the purpose of the simulation. If one is interested only in computing emergent specific intensity from a 3D computational box, and if the opacity and the source function are fully known within the box (the so-called “formal solution”), one can select a number of rays (photon paths) that emerge from the box and cross the whole extent of the box, and then solve the radiative transfer equation along such rays. Obviously, a ray does not generally go through the original grid points. One can choose any discretization of the ray, and then obtain needed values of the opacity and source function at the discretized positions by a three-dimensional interpolation. A more efficient way is to discretize the ray by taking its intersections with the individual planes defined by the grid points; in this case one deals with two-dimensional interpolations.

The problem is thus reduced to a set of 1D problems. Any method capable of solving the transfer equation along a ray can be used here; most popular in the astrophysical applications being the Feautrier method, its 2nd order variant (Feautrier 1964) or the 4th order variant (Auer 1976); the Discontinuous Finite Element method (Castor et al. 1992); or the “1D short-characteristics” method (Olson & Kunasz 1987). All these methods are called the “long-characteristics” scheme in some specific situations. However, the terminology is not used consistently; some studies use this term in a more restricted meaning, namely for using a 1D short-characteristics scheme for solving the transport for a set of rays passing through the whole computational region.

The situation is different when the source function is not known a priori within the computational box. The simplest situation of this sort arises when the opacity is specified in the computational domain and one can even assume LTE, but the continuum scattering is taken into account. Another case is a general non-LTE situation where line scattering also is important. The most complex situation is when dealing with the true radiation hydrodynamics, where one has to pass the values of radiation moments at all grid points to the hydrodynamics part of the code.

Even in the simplest case of LTE with electron (Thomson) scattering, the total source function depends on the mean intensity of radiation, and therefore has to be determined iteratively. The essential point is, however, that even if we are interested only in the emergent radiation, we have to determine the total source function, and therefore the specific intensities, in all grid points to be able to interpolate it into the discretized positions along the rays.

Although the long-characteristics scheme is useful for some purposes, such as allowing for an efficient parallel scheme within domain decomposition (Heinemann et al. 2006), we chose the short-characteristics scheme. The classical setup of the long-characteristics method is to pass a ray through each grid point through the whole computational box, in which case the number of necessary operations is obtained as follows. Let us assume that the computational box is discretized using N points in each direction, so that there are N3 grid points (and also  ≈ N3 elementary cells). A long-characteristics scheme would consider N3 rays. Solving the transfer equation along one ray, for one direction and one frequency, would require O(N) operations, an exact value depending on how exactly are the necessary interpolations performed, assuming that they require O(1) operations. For one direction and one frequency one thus needs O(N4) operations. In contrast, the short-characteristics scheme solves the transfer equation within one cell only, and therefore the total number of operations is proportional to the number of cells, that is O(N3) operations.

To avoid confusion, we mention that when using the long-characteristics scheme for obtaining an emergent intensity only, it requires also O(N3) operations, because we have N2 rays; each going through one grid point on the boundary plane that faces an “observer”, and a transfer solution along each ray requires O(N) operations. In this particular case, the scaling of the computer time is the same as for the short-characteristics scheme. In some cases, the use of long-characteristics scheme may even be more computationally advantageous than the short-characteristics, and therefore some codes (e.g., ASSeT, SCATE) use the short-characteristics scheme when dealing iteratively with scattering, and use the long-characteristics scheme to evaluate the emergent intensity for the final converged model.

We note that there are several possible numerical methods that propagate an information form one edge of the cell to another, whose numbers of operations also scale as the number of cells, O(N3). In astrophysics, the vast majority of approaches is based on the short-characteristics scheme, and this is what we adopt in this paper as well. In this paper we assume LTE and no scattering. In this case, the absorption coefficient χν is a function of the local mass density ρ and the local temperature T, and the source function Sν is equal to the local Planck function Bν(T). With these two assumptions (LTE and no scattering), we could have adopted the long-characteristics method as well, but in view of our intended future development of the solver to treat scattering and non-LTE situations, the short-characteristic scheme is clearly the appropriate choice.

3. 3D short-characteristics

In the context of astrophysical radiative transfer, the method was first used by Mihalas et al. (1978), and later by (Olson & Kunasz 1987; Kunasz & Auer 1988; Kunasz & Olson 1988), and subsequently in many studies – see references below.

3.1. Overview of the method

Let us consider, at time t, the radiation propagating in direction n. The optical depth from position r to position r + sn, where s is the path length between these two positions, is: (4)For a radiation propagating from an upwind position ru to a current position rc, the integral form of the formal solution of the time-independent RTE is: (5)where , , n is the unit vector along the straight line (ru,rc), τ is the optical depth from ru to an intermediate position ru + sn, with , and τuc is the optical depth from ru to rc.

We consider a three-dimensional Cartesian grid (see Fig. 1) defined in the (O,x,y,z) coordinate system with the unit vectors (ex,ey,ez). The grid cells are rectangular boxes with irregular spacing in each direction. The sizes of the cells in each direction depend on the position in this direction: Δxcell(x), Δycell(y), and Δzcell(z). The state parameters, i.e., the mass density, the temperature, and the velocity components in the observer’s frame, are defined at each grid point. They are provided by the (radiation)(magneto)hydrodynamics (RMHD) simulations. A direction of propagation of the radiation, defined by the unit vector n, is specified by the polar angle θ between n and the unit vector ez, and by the azimuthal angle ϕ, between ex and the projection of n on the x-y plane. In order to span the whole 4π sr angular domain, we impose θ ∈  [0]  and ϕ ∈  [0,2π [.

The task is to calculate the specific intensity Ic in a current point Mc in direction n. Mc is defined by the intersection of the cells marked with indices i, i+1 in x-direction, j, j+1 in y-direction, and k, k+1 in z-direction (see Fig. 1). This is accomplished by using Eq. (5). Following the propagation of the ray that goes through Mc in direction n, we define a short-characteristic by the line joining the intersection of the ray with the upwind cell, Mu, to the current point Mc (the upwind cell is defined by the indices (i,j,k) in Fig. 1). Ic can be determined if we know the following quantities: the upwind specific intensity Iu, the source function S(r,n,ν,t) and the absorption coefficient χ(r,n,ν,t) (in order to deduct the optical depths τ and τuc) along the short-characteristic from Mu to Mc. However, S and χ are specified only in the vertices of the cells (grid points). Therefore, we are essentially free to define laws of variation of these quantities along the short-characteristic, typically as low-order polynomials; we choose third degree polynomials. To do so, we need to know (S,χ) in the upwind end point Mu, (Su,χu), and in the downwind end point Md, (Sd,χd). Md is the intersection of the short-characteristic with the downwind cell, which is defined by the indices (i + 1,j + 1,k + 1) in Fig. 1. We also choose third degree polynomials to get the upwind quantities (Iu,Su,χu) and the downwind quantities (Sd,χd), through interpolations from the values in the neighboring grid points. Further, in Sect. 4, we describe and discuss in detail the discretization of Eq. (5), and the mathematical form of the third degree polynomials.

thumbnail Fig. 1

Short-characteristics method illustrated with an example of a 3D irregularly spaced Cartesian grid. We consider all quantities in the observer’s frame. They are defined in the vertices of the cells (grid points). (O,x,y,z) is an example of a coordinate system with the unit vectors (ex,ey,ez). The specific intensity is calculated in the current point Mc, for a radiation propagating from the upwind end point Mu to the downwind end point Md. The direction n of the ray is defined by the polar angle θ and the azimuthal angle ϕ. The short-characteristic is defined by the line joining Mu and Mc. The transfer equation is solved in its integral form along this short-characteristic. The cells around point Mc are marked with the following indices: i, i+1 in x-direction, j, j+1 in y-direction, and k, k+1 in z-direction.

Open with DEXTER

3.2. Intersections of a short-characteristic with the neighboring cell faces

The coordinates of the intersections of a short-characteristic with the neighboring cell faces, i.e., the coordinates of the upwind end point Mu and the downwind end point Md, are easily determined in the local coordinate system defined by (Mc,ex,ey,ez) (see Fig. 1).

Figure 2 displays the upwind cell (top panel) and the downwind cell (bottom panel), as they are defined for the short-characteristic plotted in Fig. 1. In this example, the upwind cell is determined by indices (i,j,k), while the downwind cell is determined by indices (i + 1,j + 1,k + 1). Let us refer to fx (respectively fy, fz) as the generic name of a cell face that is perpendicular to x-axis (respectively y-axis, z-axis), and that may be intersected by a short-characteristic, i.e., to which Mu or Md may belong. In our example, fx in the upwind cell is specified by x = −Δxi, while fx in the downwind cell is specified by x = Δxi + 1, both in the local coordinate system. In the same vein, fy is y = −Δyj in the upwind cell, and y = Δyj + 1 in the downwind cell. Finally, fz is z = −Δzk in the upwind cell, and z = Δzk + 1 in the downwind cell.

thumbnail Fig. 2

Upwind cell (top) and downwind cell (bottom), for the short-characteristic MuMc displayed in Fig. 1. In this example, the direction of propagation of the radiation, n, is in the first octant, which corresponds to (θ,ϕ) ∈  [0/2 [×  [0/2 [. The local coordinate system is defined by (Mc,ex,ey,ez). The upwind end point Mu belongs to the cell face fz:z = −Δzk. However, depending on (θ,ϕ) in this octant, Mu may belong to fx:x = −Δxi, or to fy:y = −Δyj. In the general case, Δxi ≠ Δyj, Δxi ≠ Δzk, Δyj ≠ Δzk, though equalities are possible. The downwind end point Md belongs to fz:z = Δzk + 1. However, depending on (θ,ϕ) in this octant, Md may belong to fx:x = Δxi + 1 or to fy:y = Δyj + 1. See Sect. 3.2 for more details.

Open with DEXTER

In the example displayed in Figs. 1 and 2, Mu belongs to the face fz:z = −Δzk, and Md belongs to the face fz:z = Δzk + 1. However, depending on the values of (θ,ϕ), Mu may belong to fy:y = −Δyj, or to fx:x = −Δxi, and Md may belong to fy:y = Δyj + 1, or to fx:x = Δxi + 1. In addition, the upwind cell and the downwind cell are respectively (i,j,k) and (i + 1,j + 1,k + 1) only for (θ,ϕ) ∈  [0/2 [×  [0/2 [, as exemplified in Figs. 1 and 2. In fact, for each grid point, eight different upwind cells are possible, one different cell for (θ,ϕ) in each of the eight octants in the 4π sr directional domain. Let us refer to Δxu and Δxd as the generic names for the size in x-direction of, respectively, the upwind and the downwind cell, for a given current point Mc and a given direction (θ,ϕ). Similarly, let us introduce the analogous notations Δyuyd and Δzuzd for the corresponding sizes in y and z directions. Table 1 indicates, for each of the eight octants, the corresponding range of (θ,ϕ), the indices of the upwind cell, along with the values of Δxuyuzu, for any current grid point Mc. Table 2 indicates the equivalent quantities for the downwind cell.

The following generic inequalities identify the cell face to which an upwind (respectively downwind) end point belongs: The local coordinates of the upwind (respectively downwind) end point Mu(d) depend on the face it belongs to, as follows: where the factors (αxu(d),αyu(d),αzu(d)), are equal to  ± 1, depending on the octant to which the direction n belongs, as detailed in Tables 1 and 2.

Table 1

Correspondence between octant, directional angles (θ,ϕ), upwind cell indices, sizes of the upwind cell faces in each direction, and (αxu,αyu,αzu) factors.

Table 2

Same as Table 1, but for the downwind cell faces.

4. Piecewise cubic, locally monotonic, interpolation

4.1. The advantages of a piecewise, locally monotonic, interpolant

As mentioned in Sect. 3.1, solving the integral form of the radiative transfer equation with the short-characteristics method requires to define laws of variation of physical quantities within faces of the grid and along the short-characteristics. The linear law is the simplest and fastest method. However, it leads to large numerical diffusion. This means that a sharp beam is significantly dispersed as it propagates throughout a grid, which was discussed in detail by several authors (Kunasz & Auer 1988; Hayek et al. 2010; Davis et al. 2012). Moreover, a second or higher order polynomial is mandatory, in order to recover the diffusion approximation at large optical depth. A parabolic law, however, generates overshoots that may be negative and thus unphysical. This point was discussed in detail by Kunasz & Auer (1988) and Auer & Paletou (1994). An efficient method to overcome this difficulty is to use in each direction (x, y, z and the short-characteristic’s direction) a piecewise continuous, locally monotonic, interpolation function, as proposed by Auer (2003). The monotonicity suppresses any spurious extrema. Hereafter, we describe the method that we have adopted.

4.2. Cubic Hermite polynomial and weighted harmonic mean node derivative

Let us consider a physical quantity that depends on only one variable, w(x), and which is specified only in a set of m discrete positions or nodes (xi)i = 1,n, with . In other words, wi = w(xi) for i = 1,2,...,n are given. Our objective is to find an interpolation function that goes through all the data points (xi,wi) and that preserves the monotonicity of the data. A solution to this problem is provided by a piecewise cubic Hermite polynomial coupled with a weighted harmonic mean defining the derivative at each node, . We detail this solution hereafter.

Among the nodes above, let us choose two consecutive ones, xi and xi + 1 (see Fig. 3). Within the interval  [xi,xi + 1] , a cubic Hermite polynomial is defined as (see Eq. (2) of Auer 2003): (8)where (9)

thumbnail Fig. 3

Piecewise interpolation. w(x) is the function to be interpolated between the end points (or nodes) xi and xi + 1. We assume that we know w and its derivatives at the end points: wi = w(xi), wi + 1 = w(xi + 1), , . The local interpolation function, Hi(x), is a cubic Hermite polynomial. It is a third degree polynomial defined between xi and xi + 1, from the four relations: , , , . The arrows indicate the slopes of H(x) at the end points, or, equivalently, the derivatives , . See Sect. 4.2 for detailed explanations.

Open with DEXTER

By definition, the Hermite polynomial matches the interpolated function w(x) and its derivatives at both ends of the interval: The derivative at an inner node, , is defined by the weighted harmonic mean suggested by Brodlie (1980) and Fritsch & Butland (1984), and brought to the attention of the astrophysics community by Auer (2003): (11a)where The derivatives at the outer nodes, and , are defined by the local slopes: Adopting this definition of the node derivatives, each local cubic Hermite interpolant (Hi(x))i    =    1,n − 1 is monotonic in the interpolation interval  [xi,xi + 1] . More precisely, the sign of Hi(x), for x ∈] xi,xi + 1 [, is equal to the sign of . Therefore, spurious extrema never occur, which guarantees the positivity of the interpolated physical quantities. Let us note H(x) the piecewise interpolation function defined over the whole interval , so that its restriction to each interval  [xi,xi + 1]  is equal the local function Hi(x). H(x) is continuous. Moreover, it is easy to demonstrate that its derivative H′(x) is also continuous throughout the whole interval  [x1,xn] , specifically at each inner node (xi)i    =    2,n − 1.

In addition to the Hermite interpolation, Auer (2003) suggests a second possibility: using a Bézier polynomial as a local interpolant. The latter is quite close to a Hermite polynomial. Both functions match the local end point values wi, wi + 1. However, while the Hermite polynomial matches the end point derivatives , , the Bézier polynomial does not do it necessarily. The reason is as follows. Let us note Bi(x) the local Bézier polynomial defined in  [xi,xi + 1] . Bi(x), and, therefore, its derivative , depend on free parameters, called “control values” (Auer 2003). These values can be adjusted in order to control the variation of Bi(x) in the vicinity of the end points. In particular, they can be restrained within a minimum value and a maximum value (see Eq. (11) of Auer 2003), in order to guarantee the monotonicity of Bi(x). In this case, the continuity of the derivative B′(x) of the piecewise interpolant B(x), which is defined over  [x1,xn]  similarly to H(x), is not necessarily verified at the nodes. In short, whatever the definition of the node derivatives of w(x), , a piecewise Hermite interpolation function, H(x), guarantees that H′(x) matches these node derivatives and that H′(x) is continuous at the nodes. On the other hand, a piecewise Bézier interpolation function B(x), along with an adequate choice of the control values, guarantees to be locally monotonic, but B′(x) does not necessarily match the node derivatives and is not necessarily continuous at the nodes. Now, if we adopt very specific control values (see Eq. (9) of Auer 2003 in the cubic case), we can force the Bezier polynomial to be identical to the Hermite polynomial. And, as we discussed in the paragraph above, with an adequate choice of the definition of the nodes derivatives (Eqs. (11) and (12)), a piecewise Hermite interpolation function can be made locally monotonic.

Finally, let us mention another possible definition of the node derivatives , that ensures the monotonicity of a Hermite polynomial. Steffen (1990) suggests to calculate the slope at inner node (xi)i = 2,n − 1 from the unique parabola passing through the points (xi − 1,wi − 1;xi,wi;xi + 1,wi + 1), if this parabola is monotonic in  [xi − 1,xi + 1] . If not, the node derivative is defined as the common slope shared by two parabolas, which are locally monotonic in  [xi − 1,xi]  and  [xi,xi + 1] . Specific definitions are proposed for the outer nodes derivatives and . Although this possibility suggested by Steffen (1990) satisfies our requirements of monotonicity, we did not implement it in IRIS, because the calculation of such a derivative, for a given node, uses a larger cpu time than the calculation with the weighted harmonic mean formula (Eq. (11a) in this paper may be compared to Eq. (11) of Steffen 1990).

5. Implementation

5.1. Interpolations in cell faces, ghost nodes

For each direction of propagation of the radiation, we apply the 3D short-characteristics method, which consists in solving the integral form of the RTE, Eq. (5), throughout the 3D Cartesian grid points. We repeat this procedure for all the directions specified over a solid angle. The method that we employ for propagating the radiation gives a privileged role to z-axis, in the sense that we are sweeping the 3D grid gradually, z-plane by z-plane, in the direction of increasing (or decreasing) z. In an astrophysical context, z-axis should represent a global direction of the energy transport in an object; that is the direction from the interior to the outer layers of a stellar or a planetary atmosphere, the vertical distance from an accretion disk plane, or the direction of an accretion column in a young stellar object.

Boundary conditions must be defined in appropriate planes, which depend on the direction of propagation of the radiation. Assuming that the grid is defined by Nx cells in x-direction, Ny cells in y-direction, and Nz cells in z-direction, the positions of the grid nodes in each direction may be noted as (x0,x1,...,xNx), (y0,y1,...,yNy), and (z0,z1,...,zNz). The sizes of the cells in each direction are then noted (Δx1,...xNx), (Δy1,...yNy), and (Δz1,...zNz), and they are defined by: (13)In agreement with our notations in Fig. 1, a propagation in the first octant corresponds to directional angles (θ,ϕ) ∈  [0/2 [×  [0/2 [, and, therefore, to a propagation along increasing x, increasing y, and increasing z. Consequently, the boundary specific intensities must be known at the bottom frontiers of the 3D Cartesian grid, defined by the following planes x = x0, y = y0, z = z0.

As explained in Sect. 3.1, calculating the formal solution of the monochromatic RTE in its integral form, at a given grid point and for a given direction of propagation of the radiation, requires the determination of (Iu,Su,χu) at the intersection of the short-characteristic with the upwind face, and (Sd,χd) at the intersection of the short-characteristic with the downwind face.

thumbnail Fig. 4

Examples of cubic monotonic interpolations in cell faces intersected by a short-characteristic, for a ray propagating in the first octant, i.e., along increasing x, y, and z. Two configurations are possible. The first one is displayed in the top panel, where all the quantities are known in the grid points neighboring the cell face. This situation always occurs for the upwind and downwind source functions and opacities (Su,χu), (Sd,χd). But, for the upwind specific intensity Iu, it occurs only when horizontal faces fz:z =  ± Δzk or boundary faces are intersected. The second configuration is displayed in the bottom panel. It impacts the determination of Iu when vertical non-boundary faces, fx:x =  ± Δxi or fy:y =  ± Δyj, are intersected. In this case, the upwind specific intensity is not known in all the grid points neighboring the cell face. See Sect. 5.1 for detailed explanations.

Open with DEXTER

The opacities χ are known at each point of the three-dimensional grid. And, by definition of a formal solution, the source function S is known at each grid point. Then, Su, χu, (respectively Sd, χd) are determined by successive cubic monotonic interpolations in each of the two directions of a given upwind (respectively downwind) horizontal of vertical face. The top panel of Fig. 4 provides an illustration of a sequence of interpolations in a horizontal upwind cell face. This example corresponds to the case shown in the top panel of Fig. 2. The upwind face is identified by indices i and j and belongs to the plane defined by fz:z = −Δzk in the local coordinate system (see Sect. 3.2). The ray propagates in the first octant, i.e., along increasing x, y, and z. The black dots indicate the grid points. Mu is the upwind end point. is the projection of the current point Mc on fz, and coincides with a grid point. Quantities defined in the grid points are interpolated in y-direction at the yu position of Mu. Then, the interpolated values, in positions Ai − 2, Ai − 1, Ai, Ai + 1, are themselves interpolated in x-direction at the xu position of Mu, which provides the upwind quantity Su, χu. The curved arrows indicate the grid points which contribute to the interpolations. Four points are necessary to accomplish each interpolation (see Sect. 4.2 along with Eqs. (8) and (11)).

Unlike opacities and source functions, the upwind specific intensity is not always known in all the grid points neighboring the upwind cell face. The sweeping procedure throughout the 3D Cartesian grid is made z-plane by z-plane, as described above in the first paragraph of this section. Therefore, for each new z-plane to be processed, the specific intensity is known at each grid point of the upwind z-plane. Accordingly, the procedure described above for S and χ works if Mu belongs to a horizontal face. However, the specific intensity is not known at all grid points of a preceding vertical plane (fx:x =  ± Δxi, or fy:y =  ± Δyj, in the local coordinate system), except if it is a boundary plane. The bottom panel of Fig. 4 provides an example of a sequence of interpolations of the specific intensity in a vertical upwind cell face. In this case, z = zk is the current z-plane, z = zk − 1 is the upwind one. Mu belongs to the vertical face identified by indices j and k, and by fx:x = −Δxi. We still consider a ray propagating in the first octant (increasing x, y, z). The specific intensity is fully known in the preceding z = zk − 2 and z = zk − 1 planes, where we can perform interpolations in y-direction at the yu position of Mu. This provides the values of the specific intensity in points Ak − 2 and Ak − 1. The crosses indicate the grid points in which the specific intensity is not known. Since the sweeping follows the direction of propagation of the ray, it is performed along increasing y, in the z = zk plane. Therefore, the specific intensity is known in the three consecutive grid points ((yj − 2,zk), (yj − 1,zk), (yj,zk)), but not in (yj + 1,zk). It is also not known in the z = zk + 1 plane (except in y = yj − 2 if this y-plane is a boundary). Consequently, we can determine the specific intensity in point Ak in z = zk plane, with a cubic Hermite interpolation in y-direction, using the values in ((yj − 2,zk), (yj − 1,zk), (yj,zk)). Finally, the interpolated values in Ak − 2, Ak − 1, Ak are themselves interpolated in z-direction at the zu position of Mu, which provides the upwind specific intensity Iu. For the last two interpolations, since we know the intensities in only three grid points, the derivatives of the cubic Hermite polynomials in the edge points, in y-direction and Ak in z-direction, are provided by the local slopes as defined in Eqs. (12).

We have added at each boundary of each of the three x, y, z dimensions, two layers of ghost nodes that embed the initial 3D Cartesian grid. These ghost nodes are necessary if we want to apply the above interpolation procedure uniformly to all the grid points, which spares us from testing whether the grid points are close to the boundaries or not. In x-direction, for example, and for a propagation in increasing x, ghost nodes are required for determining the upwind quantities at x1, and for determining the downwind quantities at xNx − 1 and xNx. The ghost nodes (x-2,x-1) and (xNx + 1,xNx + 2) are added. Their positions per se do not matter, as we define in these nodes positive linear extrapolations of the state parameters. Then, it is easy to show that we can still use the cubic Hermite polynomial interpolants in x1, xNx − 1 and xNx, and that the derivatives in these nodes, even though calculated with the weighted harmonic mean, are equal to the local slopes defined by Eqs. (12).

5.2. The discretized integrated RTE along a short-characteristic

In the preceding subsection, we have described our procedure to determine the upwind quantities (Iu,Su,χu), and the downwind quantities (Sd,χd). Once we know them, we are in a position to calculate the integral form of the RTE, Eq. (5).

To do so, we define laws of variation, along the short-characteristics, of the source function and the opacity. We use a monotonic cubic Hermite polynomial between the upwind endpoint Mu and the current grid point Mc (see Fig. 1). Since we know the values of S and χ in three points, Mu, Mc, and the downwind endpoint Md, the derivatives of these quantities at Mu, and , are provided by the slopes defined by Eq. (12a), while the derivatives at Mc, and , are calculated with the weighted harmonic mean, Eq. (11). Then, the RTE (Eq. (5)) can be integrated analytically. We obtain the following expression: (14a)where The optical depth τuc is calculated by integrating analytically Eq. (4) from Mu to Mc. We obtain: (15)where suc is the path length from Mu to Mc.

When τuc is small enough (≲5 × 10-2), we use Taylor expansions up to the third order of Eqs. (14). The bottom line is that the RTE can always be calculated, including the two asymptotic limits: the optically thin or even transparent path for which τuc equals zero, and, therefore, Ic → Iu, and the optically thick or even opaque path for which τuc tends to infinity, and, therefore, Ic → Sc.

5.3. The velocity gradient effect (Doppler shift)

We consider a medium with a non-zero macroscopic velocity with respect to the observer’s frame. In this case, the velocity gradient between two positions in the medium causes a Doppler shift of photon frequencies, and therefore a shift of the corresponding spectral lines. This effect must be properly taken into account in an implementation of the short-characteristics method. If the velocity difference between the current point Mc and the upwind end point Mu (see Fig. 1) is large enough, then a simple interpolation scheme may underestimate or overestimate opacities.

To illustrate this point, let us consider a single emission line and a given frequency at which we solve the RTE. The following situation may happen. At Mu, the frequency is located at a wing of the line, quite far from the center of the line, such that the absorption coefficient is very small compared to the line center value. And, at Mc, the line is shifted, such that this same frequency is at the opposite side of the line, in the other wing, far from the center. The law of variation then infers a very weak opacity for this line throughout the path from Mu to Mc. In fact, the line goes through its maximum value at this frequency, somewhere between Mu and Mc, and this is not taken into account. This issue was pointed out by van Noort et al. (2002) who emphasized the fact that the opacity along a short-characteristic may be underestimated by several orders of magnitude. In addition to that example, we provide here a case of overestimation. Let us consider now two lines that are quite distant from each other. We may face the following situation. At Mu, the calculation frequency can be located close to the center of one line. And, at Mc, due to the shift of the lines, this frequency can be close to the maximum of the second line. Consequently, a variation law defined between Mu and Mc overestimates the opacity if the two line centers are distant enough, such that we miss the trough between them.

A common solution of these problems is provided by a subgriding (e.g., van Noort et al. 2002). Let n be the macroscopic velocity in the observer’s frame, and its projection on the direction of photon propagation n. The idea of subgriding is to subdivide a short-characteristic into a set of subintervals, so that the difference in between two points of this subdivision remains small compared to the thermal velocity. In IRIS, we follow this procedure. Specifically, we introduce a dimensionless parameter factor ϵD that provides to the user a control over the subdivision. Let us identify with indices l and l + 1 two points of the subdivision of a given short-characteristic. Let us consider a single line. We note ν0 the transition frequency (emission or absorption of a photon) in the particle’s frame (or comoving frame). ν0 coincides with the center of the line in the fluid rest frame where v = 0, but this frequency is shifted to the value νctr in the observer’s frame where v ≠ 0, as follows, (16)where we consider a Doppler shift of frequencies to first order in , since we assume the velocities to be non relativistic (see Sect. 2). νctr is the center of the line in the observer’s frame. The shift of νctr between position l and position l + 1 is given by (17)The Doppler width associated with this line, in a given position in the medium, is (18)with c being the speed of light, and Vth the thermal velocity. The latter is defined as (19)where k is the Boltzmann constant, m the mass of the particle. ϵD is a control parameter defined such that the Doppler shift of the center of the line between position l and position l + 1 is bounded as follows: (20)The quantity is an average of the local Doppler width between its value in l position and its value in l + 1 position. It is easy to demonstrate that this inequality is independent of the transition frequency in the particle’s frame ν0, and that it is equivalent to the following condition on the gradients of the projected velocities: (21)ϵD is a tunable parameter, whose value depends on the required accuracy of the solution (ideally between 1/3 and 1). The smaller it is, the higher accuracy is achieved; however at the expense of increased computer demands.

6. Radiation moments and angular integration

In addition to computing specific intensity, one can calculate important angle-averaged quantities related to radiation moments, such as the mean intensity J(r,ν,t), the radiation flux vector F(r,ν,t), and the radiation pressure tensor P(r,ν,t). These are defined as where dΩ is the elementary solid angle, and the symbol designates an integration over 4π sr. The flux vector and the radiation pressure tensor can also be written in component form: (25)(26)In Cartesian coordinates, the three directions (1, 2, 3) are the unit vectors (ex,ey,ez) introduced in Sect. 3.1, and the components of the radiation propagation vector n are (27)IRIS computes each of the three components of the radiation flux vector, Fx, Fy, Fz. Since the radiation pressure tensor is symmetric (Mihalas 1978; Mihalas & Mihalas 1984), it is defined by six components, Pxx, Pyy, Pzz, Pxy, Pxz, Pyz, all of them being computed by IRIS.

The angular integrations are numerically performed by using quadratures. A quadrature is a set of M discretized direction vectors nm = (nmx,nmy,nmz) and the corresponding weights wm, which are the numerical representation of, respectively, the direction vector n and the elementary solid angle dΩ. An appropriate quadrature for an integration over the whole 4π sr range should satisfy at least the first two of the three following conditions:

  • 1.

    symmetry: the directions of the quadrature are invariant after rotation of 90° around the three coordinate directions (ex,ey,ez);

  • 2.

    for an isotropic radiation field, the three moments should be numerically exact;

  • 3.

    for an isotropic radiation field, the half first moment should be numerically exact.

The second condition reads as follows:

where δ is the unit tensor. The third condition is optional, and is related to the following issue. One may need to calculate a radiative flux over a surface or at a wall. In this case, the flux is determined by an integration over 2π sr, instead of 4π sr: . The third condition reads then: (29)where ei stands for ex, ey, or ez, and where we consider surfaces or walls perpendicular to one of these three directions.

We have implemented in IRIS the quadratures of type A proposed by Carlson (1963), and whose construction principle is summarized by Bruls et al. (1999). They satisfy conditions 1 and 2, but not condition 3. We also have implemented the quadratures built by Lathrop & Carlson (1965). The latter provide several sets of quadratures that satisfy conditions 1 and 2, but not condition 3, with 24 up to 288 directions. They also provide several sets of another type of quadratures that satisfy the three conditions above, with 24 up to 80 directions.

We also have implemented the Gaussian quadrature (Press et al. 1992). The latter is not appropriate for the calculation of the radiation moments, because it does not satisfy any of the three conditions above and leads to serious non physical asymmetries, unless one chooses a large number of directions (typically, more than 100). However, it is suitable for the integration of the incoming specific intensity over a viewing solid angle, in order to determine the radiative power per unit surface that receives a detector looking at a given object.

7. Periodic boundary conditions

We consider a three-dimensional medium with an infinite extension in the horizontal plane (x,y), and a finite extension along the vertical z-axis between its lower boundary z0 and its upper boundary zNz. We assume that this medium has a double periodicity, one in x-direction, and one in y-direction. The boundary conditions are known at z0 and zNz. For example, we may consider a non-irradiated stellar atmosphere with no incoming radiation at the outer surface zNz and a black body radiation at its inner surface z0. The computational grid ranges from z0 to zNz in vertical direction, from x0 to xNx, and from y0 to yNy in the horizontal plane, so that (xNx − x0) defines one x-period and (yNy − y0) defines one y-period. Now, the 3D short-characteristics method consists in solving the integral form of the RTE by propagating the rays throughout a computational domain, from up to three upwind boundary sides in which the specific intensity is assumed to be known, down up to the three other faces of the domain. For horizontally periodic media, while the vertical boundary conditions are specified explicitly, the lateral boundary conditions are defined implicitly, such that, for any physical quantity f(ν,x,y,z), we have the following relations: Consequently, when the upwind end point of a short-characteristic intersects a lateral boundary, we prolong this characteristic, which becomes a long-characteristic, until it intersects a horizontal face, following the suggestions by Auer et al. (1994) and Fabiani Bendicho (2003). For a given direction of propagation, and for each z-layer, this treatment affects only the rows that are adjacent to the boundary faces.

Figure 5 illustrates our point with an example showing a propagation of the radiation in the first octant, i.e., along increasing x, y, and z. Two characteristics are plotted from the current point Mc in the plane x = x1, which is adjacent to the boundary plane x = x0. The upwind end point Mu of the short-characteristic (Mu,Mc,Md) belongs to a horizontal face. Therefore, the specific intensity can be calculated by interpolations from the values at the neighboring grid points in the upwind horizontal plane. If we consider the short-characteristic , the upwind end point belongs to the horizontal boundary plane x = x0, in which the specific intensity is not explicitly defined. The short-characteristic is then lengthened, and becomes a long-characteristic, down to its first intersection, , with a horizontal face. The long-characteristic can go through several cells (only one cell in the example, for the clarity of the figure) beyond the vertical boundary face, before hitting a horizontal face.

thumbnail Fig. 5

Extract of a 3D Cartesian grid. Unlike the general case presented in Fig. 1, we consider here the particular case of a medium with a double horizontal periodicity, in x-direction and in y-direction. In x-direction, one period is defined from x0 to xNx. The plane x = x0 is a boundary plane. The cell defined by indices (− 1,j,k), drawn with thick lines, is identical to the cell (Nx,j,k) (not drawn). The upwind specific intensity is known at Mu for the short-characteristic (Mu,Mc,Md), since Mu belongs to a horizontal face. On the other hand, the upwind specific intensity is not known at for the short-characteristic , because belongs to a vertical boundary face. Therefore, the characteristic is prolonged up to the first intersection with a horizontal face, . See Sect. 7 for more explanations.

Open with DEXTER

When the medium is not periodic, we start the sweeping of a given z-plane at the grid point that is the closest to the vertical boundary planes. For example, if the ray propagates in the first octant, these planes are x = x0 and y = y0, therefore the first grid point in which we calculate the specific intensity is (x1,y1). Then, the specific intensity is calculated progressively in all the grid points in the z-plane along x-direction and y-direction. Now, as we consider in this section a medium with a double horizontal periodicity, we can start the sweeping at any grid point in the z-plane, provided that the grid points where the specific intensity is calculated span one x-period and one y-period. We mark this starting point as (xi0,yj0). The indices i0 and j0 are defined so as to minimize the number of cells that are intersected by the long-characteristics, thus saving computing time. This way, i0 is chosen such that the corresponding upwind cell has the largest size in x-direction, Δxi0 = Δxmax. Similarly, we choose j0 such that Δyj0 = Δjmax. Then, starting from (xi0,yj0), we calculate the specific intensity with the long-characteristics method at yj0, for the grid points along x-direction that span one x-period. In the same vein, we use the long-characteristics method at xi0, for the grid points along y-direction that span one y-period. Once these two first rows have been treated, we resume the short-characteristics method for all the next rows of the current z-layer. Note that here we do not introduce ghost nodes (see Sect. 5.1), since the medium has an infinite extension in x-direction and y-direction.

The schematic diagram of Fig. 6 clarifies our point, for a radiation propagating in the first octant, therefore along increasing x, y, and z. It shows two computational domains in a given z-plane. Both span one x-period and one y-period. The original domain is depicted by the rectangle with solid thick lines. The positions of the grid nodes in the horizontal z-plane were introduced in Sect. 5.1, (x0, x1, ..., xNx) and (y0, y1, ..., yNy). The vertical boundary planes are x = x0 and y = y0. The indices i0 and j0 are defined above in the preceding paragraph. The RTE is solved in the computational domain depicted by the rectangle with thick dotted lines. The positions of the grid nodes of this domain are (xi0, xi0 + 1, ..., xi0 + Nx) and (yj0, yj0 + 1, ..., yj0 + Ny), with the following correspondence of the grid sizes, derived from the double periodicity of the medium: (31)We also have the following correspondence for the physical quantities f(ν,x,y,z): (32)We employ the long-characteristics method along x-direction from (xi0,yj0) to (xi0 + Nx,yj0), and along y-direction from (xi0,yj0) to (xi0,yj0 + Ny). Then, we resume our usual short-characteristics method for the next rows of the current z-layer. Once calculated, the specific intensity is defined in the original computational domain, with a simple rearrangement of indices as per Eq. (32). Note that some long-characteristics may consist merely of short-characteristics. This is the case for a given grid point (xi,   i ∈ [[i0,Nx]],yj0), if the corresponding Δxi and Δymax are large enough, so that the related upwind end point intersects the horizontal upwind face. Similarly, this is the case for (xi0,yj,   j ∈ [[j0,Ny]]), if Δyj and Δxmax are large enough.

We mention here that it is also possible to handle periodic conditions by iterative methods, as done for instance by van Noort et al. (2002) and Davis et al. (2012). However, such an approach can lead to a slow convergence in the case of optically thin media, as noted by van Noort et al. (2002).

thumbnail Fig. 6

Original domain (solid thick lines) and calculation domain (dotted thick lines) at a given z-plane, for a medium that is periodic in the horizontal (x, y) plane, and for a radiation propagating in the first octant, i.e., along increasing x, y, and z. Both domains span one period in x-direction and one period in y-direction. The nodes of the original domain are (x0,x1,...,xNx) and (y0,y1,...,yNy). The nodes of the calculation domain, where the specific intensity is computed, are (xi0,xi0 + 1,...,xi0 + Nx) and (yj0,yj0 + 1,...,yj0 + Ny). The indices i0 and j0 are chosen such that Δxi0 = Δxmax and Δyj0 = Δjmax. See text in Sect. 7 for a discussion.

Open with DEXTER

8. Tests of IRIS

The code underwent a number of various tests to validate each of its functionalities that we describe in this paper. We focus here on the three most important ones: the searchlight beam test, a comparison with a well-tested one-dimensional scheme, and the test of the velocity gradient effect applied on one given spectral line.

thumbnail Fig. 7

Searchlight beam test. The top left panel of the figure depicts the purpose of this test. A slanted beam, with direction n, enters a zero opacity box at its base, and emerges from it at the top. The objective is to compare the entering beam profile with the emerging one; both are theoretically the same. The bottom panels show the normalized specific intensity of the beam, as a function of positions x and y, at the bottom of the box (left panel), and at the top of the box (right panel). The top right panel displays sectional views of the beam along x-axis, at z = 0 and y = 3 (red curve), at z = 10 and y = 6.8 (green curve). The blue curve is identical to the the green one, but shifted so that its center fits the center of the red curve. See Sect. 8.1 for explanations.

Open with DEXTER

8.1. The searchlight beam test

In the context of astrophysical radiation transport, this test was first proposed by Kunasz & Auer (1988). Several authors then used it to evaluate their multidimensional short-characteristics algorithms (Stone et al. 1992; Auer & Paletou 1994; Fabiani Bendicho 2003; Fabiani Bendicho & Trujillo Bueno 1999; Hayek et al. 2010; Davis et al. 2012). The purpose is to examine how a beam propagates through a transparent medium. Theoretically, the beam crosses the medium without being absorbed or dispersed. The RTE, Eq. (5), is greatly simplified since, in this particular case, Ic = Iu. Note that the numerical counterpart of the RTE, Eq. (14) is consistent with this behavior when the opacity is zero, as pointed out in Sect. 5.2. Therefore, the calculation uses only the upwind interpolation of the specific intensity. This test challenges the capabilities of the piecewise cubic, locally monotonic, interpolation technique (see Sect. 4), as applied to cell face interpolations (see Sect. 5.1).

For comparison, the numerical test was performed with the same parameters as the ones used by Hayek et al. (2010). We consider a hard-edged beam that propagates throughout a three-dimensional zero opacity box along a slanted direction defined by θ = 28.1° and ϕ = 45.0°, as shown in the sketch (top left panel) of Fig. 7. The box is made with 1003 grid points. It extends in each of the three directions (x,y,z) from 0 to 10 in arbitrary units. The beam is made with 302 grid points at the base of the box. It enters the box at the base and emerges from it at the top. The bottom panels of Fig. 7 show the normalized specific intensity of the beam as a function of x and y, at the base of the box (left panel), and at the top of the box (right panel). We verify that the beam profile is very well conserved. There is no absorption: the maximum value of the normalized specific intensity remains equal to one. The dispersion is very limited: the hard-edges are conserved, along with the size of the beam. In addition, the beam exits at the right position of the box. The positions of the four corners of the beam,  [(x1,y1);(x2,y1);(x1,y2);(x2,y2)] , are at the base:  [(1.5,1.5);(4.5,1.5);(1.5,4.5);(4.5,4.5)] , and at the top:  [(5.3,5.3);(8.3,5.3);(5.3,8.3);(8.3,8.3)] , which is consistent with the (θ,ϕ) values specified above and the size of the box along z direction. In order to show in greater detail the dispersion of the beam, we have plotted sectional views of the beam along x-axis (top right panel). The red curve represents the normalized specific intensity as a function of x, at z = 0 (base of the box) and y = 3 (middle of the beam). The green curve represents this quantity at z = 10 (top of the box) and y = 6.8 (middle of the beam). In order to compare the two profiles, we have plotted the blue curve, which is the green one artificially shifted so that its center fits the center of the red curve. Such a superposition shows the symmetry of the beam’s dispersion, which also remains small. We verified that our scheme guarantees the photon conservation: the xy-surface integral of the specific intensity at z = 0 equals the the xy-surface integral at z = 10, to the machine accuracy.

This test represents an excellent validation of the piecewise cubic, locally monotonic, interpolation scheme that we use. Such a technique almost suppresses the numerical diffusion effect of the short-characteristics method, and is more efficient than linear, parabolic (the first one is extremely diffusive, the second one introduces spurious extrema, as shown by Kunasz & Auer 1988), or even monotonic quadratic schemes (Auer & Paletou 1994).

8.2. Comparison with 1D plane-parallel models

It is important to check whether IRIS can reproduce the results provided by one-dimensional radiative transfer schemes, when the medium is made of homogeneous plane-parallel layers, i.e., in the case of a 1D plane-parallel structure. Simulating the radiative transfer for such a medium with IRIS can be achieved by imposing periodic boundary conditions in faces perpendicular to the parallel layers. Assuming that these layers are parallel to the horizontal z-planes, the specific intensity has then the following dependence, I(z,cosθ,ν). Applying the full 3D calculation with IRIS for such a medium perfectly reproduces the radiative results provided by 1D radiative models at any frequency, as we show in the example below.

We calculate the radiative transfer through a radiative shock structure made of homogeneous plane-parallel layers. The hydrodynamics structure, provided by Matthias González1, is obtained with the three-dimensional radiation hydrodynamics code HERACLES (González et al. 2007). The test case is a simulation of an experimental radiative shock generated in a tube full of Xenon assumed to be a perfect gas, with the following upstream conditions: fluid velocity = 60 km   s-1, pressure = 7 bar, temperature = 1 eV. The opacities that we use are from Michaut et al. (2004). Our objective here is focused on the comparison of the 1D and the 3D radiative models; a future paper (Ibgui et al. 2012, in prep.) will present more details of three-dimensional models of radiative shocks.

The 3D computational grid is built with 20 × 20 × 510 points in respectively x, y, and z directions. In each horizontal z-plane, all the state parameters (temperature, density, velocity) are independent of x and y. Two post-processing calculations were performed: one with IRIS applied on the 3D grid, and one with a well-tested 1D short-characteristics solver applied on a 1D grid that is extracted from the 3D one. Both were done at seventeen frequencies that range from  = 2   eV (λ = 620   nm) to  = 494   eV (λ = 2.41   nm). The velocity gradient effect was not taken into account in this test. Figure 8 shows several radiative quantities provided by IRIS and by the 1D solver. By way of example, we select one frequency,  = 296   eV, which corresponds to λ = 4.19   nm.

The top left panel of Fig. 8 displays the specific intensity as a function of z, and for different polar angles θ: 0° (green), 45° (blue), 89° (red), 135° (orange), 180° (cyan). The solid lines are the results calculated with IRIS. The symbols are the results provided by the 1D solver. Note that, for the readability of the figure, we plot the 1D results for only some of the 510 z values. Note also that we plot the 3D results calculated with several azimuthal angles ϕ(°) = (0,45,89,90,135,179,180,225,269,270,315,359). All the curves for these twelve ϕ angles are perfectly superposed, and, therefore, indistinguishable. This demonstrates that the specific intensity is independent of ϕ. In addition, we verified that for any given pair (θ,ϕ), the specific intensity at a given z-plane is the same for all the grid points of this plane, which demonstrates that this quantity is independent of x and y. Last but not least, the 1D symbols perfectly fit with the 3D curves.

The other three panels of Fig. 8 display the moments obtained by angular integration (see Sect. 6): the mean intensity J in the top right panel, the components of the radiation flux vector, Fx, Fy, Fz, in the bottom left panel, and the components of the radiation pressure tensor Pxx, Pyy, Pzz, Pxy, Pxz, Pyz, in the bottom right panel. All the curves show that the 1D and 3D results perfectly agree. We checked this result for all the seventeen tested frequencies. In addition, the 3D results provided by IRIS verify, for any frequency, the following properties that are theoretically valid for 1D plane-parallel structures. The radiative flux is zero in x and y direction, Fx = 0, Fy = 0, and it depends only on z in z-direction, Fz(z,ν). The non-diagonal components of the radiation pressure tensor are all zero: Pxy = 0, Pxz = 0, Pyz = 0. The other three pressure components depend only on z, Pxx(z,ν), Pyy(z,ν), Pzz(z,ν), and we have, for any z and ν, Pxx(z,ν) = Pyy(z,ν).

These tests demonstrate that IRIS is capable of reproducing 1D plane-parallel simulations, and, thereby, of handling periodic boundary conditions.

thumbnail Fig. 8

Comparison with 1D plane-parallel models. The figures shows the monochromatic radiative quantities, at  = 296   eV (λ = 4.19   nm), calculated with IRIS (solid lines) and with a 1D solver (symbols), in the case of a radiative shock structure made of homogeneous plane-parallel layers. The top left panel displays the specific intensity for different polar angles θ: 0° (green, square), 45° (blue, triangle), 89° (red, filled circle), 135° (orange, diamond), 180° (cyan, asterisk), and for a set of twelve azimuthal angles ϕ(°) = (0,45,89,90,135,179,180,225,269,270,315,359). There is no legend for the ϕ values, because all the curves, for a given θ, are perfectly superposed. The top right, bottom left, and bottom right panels show respectively the mean intensity J (black, square), the components of the radiative flux vector F (black, square), and the components of the radiation pressure tensor P (Pzz:blue,square;Pxx = Pyy:red,plus). See Sect. 8.2 for more details.

Open with DEXTER

8.3. Test of the velocity gradient effect

As explained in Sect. 5.3, IRIS can handle the velocity gradient effect on spectra (caused by the Doppler shifts of photon frequencies) by subdividing a short-characteristic in a set of subintervals. This subdivision is controlled by a tunable parameter ϵD, so that the Doppler shift of any spectral line between two subinterval points remains bounded by a fraction ϵD of the local Doppler width of the line, as per Eq. (20).

We illustrate the effect of our treatment in a simple ideal case that can be otherwise calculated quasi-analytically. We consider one plasma layer with a uniform temperature T in a stellar wind region with a velocity gradient, as sketched by Fig. 9, top panel. We focus on one velocity direction, referred to as z-direction. In our example, the material moves toward an observer. The layer is numerically modeled with a single cell between two grid points, z1, the farther position from the observer, and z2, the closer position to the observer. The wind velocities at grid points are V1 > 0 and V2 > V1. We assume that there is only one radiating species. We focus on the radiation of the layer itself, in the direction of the wind toward the observer, and neglect any incoming radiation in z1. Although the situation in a stellar wind is significantly different from LTE, we do assume LTE for the purpose of this test. We also neglect scattering. Therefore, the source function is uniform and equals the Planck function at the layer’s temperature, . We consider one spectral line, for which we assume, for simplicity, a Doppler profile centered at a frequency ν0 in the particle’s frame. With these assumptions, the Doppler width ΔνD is a constant, and the monochromatic specific intensity of the radiation emerging from the layer toward the observer is simply given by: (33)where τ12(ν) is the monochromatic optical depth from z1 to z2: (34)where the absorption coefficient κ(ν,z) is written as (35)where νctr(z), defined in Sect. 5.3 by Eq. (16), is the shifted line-center frequency for a velocity V(z) at position z: νctr(z)    =    ν0   (1 + V(z)/c). The quantity K is a constant that contains all the attributes of the line (population of the lower level, and oscillator strength) and the value of the constant Doppler width.

Now, since we use monotonic laws (cubic Hermite polynomials, see Sect. 5.2) to interpolate physical quantities along a short-characteristic, τ12(ν) may be written as (36)where R is the function relating the shifted transition frequency νctr to the position z within the layer: z = R(νctr), R′ being its derivative, and where we define: (37)In this simplified configuration, the optical depth appears as the convolution of a function A(ν) and a Gaussian function B(ν) (the line profile), both defined as follows: (38a)where For the test, we adopt the following typical values of velocities in winds of early massive stars: V1 = 2000   km   s-1, V2 = 2250   km   s-1, and a thermal velocity (constant in our uniform layer) of Vth = 25   km   s-1. With such a choice, any line is Doppler shifted by an amount of 10 times its Doppler width between position z1 and position z2, i.e., (νctr2 − νctr1)/ΔνD = (V2 − V1)/Vth = 10. We also adopt a temperature T = 105 K, and a layer size z2 − z1 = 109 cm.

thumbnail Fig. 9

Test of the velocity gradient effect. The top panel represents a single plasma layer in a stellar wind, with a velocity gradient. The layer is at uniform temperature T, with one radiating species, assuming LTE and no scattering. The layer is modeled in IRIS with one cell between grid points z1 and z2, with wind velocities V1 > 0 (the material approaches the observer) and V2 > V1. The red thick curve in the bottom left panel displays the monochromatic specific intensity Iν, due to the radiation that emerges from the layer at z2 toward the observer, as computed by IRIS. This emission spectrum is plotted versus photon energy  (eV) or photon wavelength λ   (Å). It is obtained with a parameter ϵD = 0.3. Too large values of ϵD result in incorrect profiles, as shown by the following curves: cyan (ϵD = 11), green (ϵD = 6), orange (ϵD = 5), blue (ϵD = 3), grey (ϵD = 2). The bottom right panel displays the corresponding transmission spectrum. For comparison, the black thick curve portrays the spectral line that would obtain if there was no velocity gradient (V2 = V1). See Sect. 8.3 for detailed explanations.

Open with DEXTER

We consider an artificial line, arbitrarily centered at hν0 = 80.00 eV in the particle’s frame. In the observer’s frame, the line is centered at νctr1 ≈ 80.53 eV at z1, and at νctr2 ≈ 80.60 eV at z2. Note also that we have defined the K value (see Eq. (35)) so that the transmission of the layer at the line center equals 0.5 in the case when there is no velocity gradient. In this case indeed, the monochromatic thermal absorption κ(ν,z) is independent of z, so that K is inferred from e − K(z2 − z1) = 0.5. This way, we avoid the extreme cases of an optically thin or an optically thick layer. The red thick curve in Fig. 9, bottom left panel, displays the monochromatic specific intensity I(ν) that emerges from the layer at position z2 toward the observer, as calculated by IRIS. This emission spectrum is displayed versus photon energy  (eV) or photon wavelength λ   (Å). For comparison, the black thick curve represents the line, centered at 80.53 eV, that would obtain if there was no velocity gradient, i.e., if V2 = V1. The bottom right panel of the figure portrays the transmission, eτ12(ν), of the layer. In the highly simplified configuration of our test, the transmission and the specific intensity are related by the linear Eq. (33). This is why this transmission spectrum looks the same, but inverted, as the emission spectrum. In the presence of a velocity gradient, the absorption of the medium decreases, is shifted (blueward if the material approaches the observer, redward if it recedes from the observer), and spreads over a larger spectral range. This also applies to the emission of the medium.

In this ideal test case, we can explain the shape and the position of the line depicted by the red curve. The optical depth is given by Eq. (36). In addition, the layer is modeled with a unique cell. So, the cubic monotonic Hermite polynomial that defines V(z) along the short-characteristic between z1 and z2 almost coincides with a linear law, since the downwind endpoint is a ghost point in which state parameters are defined by linear extrapolations (see Sect. 5.1). In these specific conditions, z is trivially related to νctr by a linear function R(νctr), so that A(ν) is merely a boxcar function. Then, the optical depth is (39)where νctr1 and νctr2 are defined by Eqs. (37). The calculation of the monochromatic specific intensity with Eq. (33), using the relation (39) above, provides a result that exactly fits the red curve obtained with IRIS. Note that this spectral line is centered at (νctr1 + νctr2)/2 ≈ 80.57 eV.

Last but not least, Fig. 9 (emission spectrum, left, or transmission spectrum, right) shows that the ϵD parameter (cf. Eq. (20)) plays a crucial role in the precision of the calculated line, through the subgriding of the short-characteristic. The cyan curve shows the profile obtained with ϵD = 11, for which we checked that the short-characteristic is not subdivided. This profile is composed of two peaks, which are located in the two extreme positions of the line center, νctr1 and νctr2, with a large trough in the middle. Such a discrepancy between this profile and the correct red one may result in erroneous interpretations of observed spectra when compared with synthetic spectra. As ϵD decreases, the trough is being gradually filled, as shown by the other curves obtained with different ϵD: 6 (green), 5 (orange), 3 (blue), 2 (grey). This test verifies that the line shape starts to stabilize from ϵD = 1. Tiny waves remain in the profile calculated with ϵD = 1. This profile is not shown in the figure because these waves are barely perceptible, so that it almost coincides with the red profile. Our numerical tests indicate that the short-characteristic is subdivided in 14 subintervals for ϵD = 1, whereas it is subdivided in 46 subintervals for ϵD = 0.3. Not surprisingly, increasing the precision involves a larger computational cost. Such a test confirms that the ideal value of ϵD lies somewhere between 1/3 and 1 (see Sect. 5.3).

We stress again that the profile depicted by the red curves in Fig. 9 has been obtained in the ideal case of a single line with a Doppler profile, under the assumption of a unique layer at uniform temperature, in LTE, and with a linear velocity variation within the layer. In the general case, none of these situations holds, which results in a more complex profile. Only a full calculation with our subgriding method can provide the actual profile.

9. Conclusion

In this paper, we have described IRIS, a new generic three-dimensional spectral radiative transfer code, which solves the monochromatic static 3D radiative transfer equation in the observer’s frame, in Cartesian coordinates. We drop the time derivative in the transfer equation because we consider situations in which the dynamical timescales are large compared to the photon free-flight time, such that the studied medium is assumed to be in a frozen state during the time when the radiation field propagates throughout its structure. The code is primarily intended for post-processing any hydrodynamics snapshot, i.e., any structure provided at a given instant by a (radiation)(magneto)hydrodynamics simulation. That way, IRIS can determine, at each instant, the radiation field of a non-stationary flow. We consider non-relativistic flow velocities. Currently, we assume local thermodynamic equilibrium.

IRIS is mainly a diagnostic tool to be used for comparisons between predicted spectra, maps or images one the one hand, and astrophysical observations or laboratory astrophysics measurements on the other hand. It computes the monochromatic specific intensity and optical depth at any position within a hydrodynamics structure, for any direction. Associated with appropriate opacities specified by the user in a dedicated module, it calculates synthetic spectra and maps or images emerging from the studied structure or astrophysical object. It can also calculate, from the specific intensity and through angular integrations, the moments of the radiation field: the mean intensity, the radiation flux vector, and the radiation pressure tensor, at any position within the studied medium.

IRIS works with any 3D Cartesian grid, the latter being nonuniform in each direction. The code uses a short-characteristics solver to determine the formal solution of the radiative transfer equation. We have implemented a very efficient piecewise cubic, locally monotonic, interpolation technique, that considerably reduces the numerical diffusion effects of the short-characteristics method. The latter is used for interpolating any physical quantity in the cell faces of the computational grid, and for defining laws of variation of physical quantities along short-characteristics.

IRIS is able to handle horizontal periodic boundary conditions of a simulation box. This configuration occurs for a medium with an infinite (that is, large compared to the extent of the computational box) extension in its horizontal plane and a double periodicity in this plane, along with a finite extension in its vertical direction, such as a stellar atmosphere or an accretion disk.

Since it is formulated in the observer’s frame, IRIS can deal with any (non-relativistic) non-monotonic macroscopic velocities. The code can be applied to a large number of radiating astrophysical objects or structures, such as accretion shocks or jets in young stellar objects, stellar atmospheres, exoplanet atmospheres, accretion disks, rotating stellar winds, and cosmological structures. IRIS has already been applied to predict X-UV spectra of laboratory generated radiative shocks (see conference proceedings Ibgui et al. 2012a,b; we are also writing a paper on this topic: Ibgui et al. 2012, in prep.).

We envision various extensions in the near future. First, we will implement an iterative method based on the accelerated lambda iteration (ALI) technique to handle scattering and NLTE effects. We will interface IRIS with the NLTE stellar atmosphere code TLUSTY (Hubeny 1988; Hubeny & Lanz 1995), from which we take the treatment of atomic data, and all local physics (opacities, and the scheme for solving the set of statistical equilibrium equations using the preconditioning method to obtain NLTE level populations). Later, extending IRIS to polarized radiative transfer will provide the code with the capability to diagnose magnetic fields.

On the computational side, IRIS is written in Fortran 95. The current version runs on a single processor. We plan to parallelize the code with the message passing interface (MPI) library. Post-processing of (R)(M)HD calculations performed on adaptively refined grids generated with the adaptive mesh refinement (AMR) technique will benefit from 3D nonuniform Cartesian grids in IRIS, though its full implementation will require additional work.


1

AIM, CEA/DSM/IRFU, CNRS, Université Paris Diderot, 91191 Gif-sur-Yvette, France

Acknowledgments

We are grateful to Matthias González (AIM, CEA/DSM/IRFU, CNRS, Université Paris Diderot, 91191 Gif-sur-Yvette, France) for providing us with the 3D hydrodynamics results generated by the HERACLES code. These simulations, which mimic a 1D radiative shock in ideal gas, were used as an input for the radiative transfer code IRIS, in our test to reproduce the radiative field in a 1D plane-parallel medium (see Sect. 8.2). We are also grateful to the referee for a prompt report and insightful suggestions that have improved the paper. The work is supported by French ANR, under grant 08-BLAN-0263-07. The work of I.H. was supported in part by NSF grant AST-0807496. I.H. also acknowledges the travel support from the Observatoire de Paris, and the Visiting Professorship at the Université Pierre et Marie Curie where a part of work was done.

References

All Tables

Table 1

Correspondence between octant, directional angles (θ,ϕ), upwind cell indices, sizes of the upwind cell faces in each direction, and (αxu,αyu,αzu) factors.

Table 2

Same as Table 1, but for the downwind cell faces.

All Figures

thumbnail Fig. 1

Short-characteristics method illustrated with an example of a 3D irregularly spaced Cartesian grid. We consider all quantities in the observer’s frame. They are defined in the vertices of the cells (grid points). (O,x,y,z) is an example of a coordinate system with the unit vectors (ex,ey,ez). The specific intensity is calculated in the current point Mc, for a radiation propagating from the upwind end point Mu to the downwind end point Md. The direction n of the ray is defined by the polar angle θ and the azimuthal angle ϕ. The short-characteristic is defined by the line joining Mu and Mc. The transfer equation is solved in its integral form along this short-characteristic. The cells around point Mc are marked with the following indices: i, i+1 in x-direction, j, j+1 in y-direction, and k, k+1 in z-direction.

Open with DEXTER
In the text
thumbnail Fig. 2

Upwind cell (top) and downwind cell (bottom), for the short-characteristic MuMc displayed in Fig. 1. In this example, the direction of propagation of the radiation, n, is in the first octant, which corresponds to (θ,ϕ) ∈  [0/2 [×  [0/2 [. The local coordinate system is defined by (Mc,ex,ey,ez). The upwind end point Mu belongs to the cell face fz:z = −Δzk. However, depending on (θ,ϕ) in this octant, Mu may belong to fx:x = −Δxi, or to fy:y = −Δyj. In the general case, Δxi ≠ Δyj, Δxi ≠ Δzk, Δyj ≠ Δzk, though equalities are possible. The downwind end point Md belongs to fz:z = Δzk + 1. However, depending on (θ,ϕ) in this octant, Md may belong to fx:x = Δxi + 1 or to fy:y = Δyj + 1. See Sect. 3.2 for more details.

Open with DEXTER
In the text
thumbnail Fig. 3

Piecewise interpolation. w(x) is the function to be interpolated between the end points (or nodes) xi and xi + 1. We assume that we know w and its derivatives at the end points: wi = w(xi), wi + 1 = w(xi + 1), , . The local interpolation function, Hi(x), is a cubic Hermite polynomial. It is a third degree polynomial defined between xi and xi + 1, from the four relations: , , , . The arrows indicate the slopes of H(x) at the end points, or, equivalently, the derivatives , . See Sect. 4.2 for detailed explanations.

Open with DEXTER
In the text
thumbnail Fig. 4

Examples of cubic monotonic interpolations in cell faces intersected by a short-characteristic, for a ray propagating in the first octant, i.e., along increasing x, y, and z. Two configurations are possible. The first one is displayed in the top panel, where all the quantities are known in the grid points neighboring the cell face. This situation always occurs for the upwind and downwind source functions and opacities (Su,χu), (Sd,χd). But, for the upwind specific intensity Iu, it occurs only when horizontal faces fz:z =  ± Δzk or boundary faces are intersected. The second configuration is displayed in the bottom panel. It impacts the determination of Iu when vertical non-boundary faces, fx:x =  ± Δxi or fy:y =  ± Δyj, are intersected. In this case, the upwind specific intensity is not known in all the grid points neighboring the cell face. See Sect. 5.1 for detailed explanations.

Open with DEXTER
In the text
thumbnail Fig. 5

Extract of a 3D Cartesian grid. Unlike the general case presented in Fig. 1, we consider here the particular case of a medium with a double horizontal periodicity, in x-direction and in y-direction. In x-direction, one period is defined from x0 to xNx. The plane x = x0 is a boundary plane. The cell defined by indices (− 1,j,k), drawn with thick lines, is identical to the cell (Nx,j,k) (not drawn). The upwind specific intensity is known at Mu for the short-characteristic (Mu,Mc,Md), since Mu belongs to a horizontal face. On the other hand, the upwind specific intensity is not known at for the short-characteristic , because belongs to a vertical boundary face. Therefore, the characteristic is prolonged up to the first intersection with a horizontal face, . See Sect. 7 for more explanations.

Open with DEXTER
In the text
thumbnail Fig. 6

Original domain (solid thick lines) and calculation domain (dotted thick lines) at a given z-plane, for a medium that is periodic in the horizontal (x, y) plane, and for a radiation propagating in the first octant, i.e., along increasing x, y, and z. Both domains span one period in x-direction and one period in y-direction. The nodes of the original domain are (x0,x1,...,xNx) and (y0,y1,...,yNy). The nodes of the calculation domain, where the specific intensity is computed, are (xi0,xi0 + 1,...,xi0 + Nx) and (yj0,yj0 + 1,...,yj0 + Ny). The indices i0 and j0 are chosen such that Δxi0 = Δxmax and Δyj0 = Δjmax. See text in Sect. 7 for a discussion.

Open with DEXTER
In the text
thumbnail Fig. 7

Searchlight beam test. The top left panel of the figure depicts the purpose of this test. A slanted beam, with direction n, enters a zero opacity box at its base, and emerges from it at the top. The objective is to compare the entering beam profile with the emerging one; both are theoretically the same. The bottom panels show the normalized specific intensity of the beam, as a function of positions x and y, at the bottom of the box (left panel), and at the top of the box (right panel). The top right panel displays sectional views of the beam along x-axis, at z = 0 and y = 3 (red curve), at z = 10 and y = 6.8 (green curve). The blue curve is identical to the the green one, but shifted so that its center fits the center of the red curve. See Sect. 8.1 for explanations.

Open with DEXTER
In the text
thumbnail Fig. 8

Comparison with 1D plane-parallel models. The figures shows the monochromatic radiative quantities, at  = 296   eV (λ = 4.19   nm), calculated with IRIS (solid lines) and with a 1D solver (symbols), in the case of a radiative shock structure made of homogeneous plane-parallel layers. The top left panel displays the specific intensity for different polar angles θ: 0° (green, square), 45° (blue, triangle), 89° (red, filled circle), 135° (orange, diamond), 180° (cyan, asterisk), and for a set of twelve azimuthal angles ϕ(°) = (0,45,89,90,135,179,180,225,269,270,315,359). There is no legend for the ϕ values, because all the curves, for a given θ, are perfectly superposed. The top right, bottom left, and bottom right panels show respectively the mean intensity J (black, square), the components of the radiative flux vector F (black, square), and the components of the radiation pressure tensor P (Pzz:blue,square;Pxx = Pyy:red,plus). See Sect. 8.2 for more details.

Open with DEXTER
In the text
thumbnail Fig. 9

Test of the velocity gradient effect. The top panel represents a single plasma layer in a stellar wind, with a velocity gradient. The layer is at uniform temperature T, with one radiating species, assuming LTE and no scattering. The layer is modeled in IRIS with one cell between grid points z1 and z2, with wind velocities V1 > 0 (the material approaches the observer) and V2 > V1. The red thick curve in the bottom left panel displays the monochromatic specific intensity Iν, due to the radiation that emerges from the layer at z2 toward the observer, as computed by IRIS. This emission spectrum is plotted versus photon energy  (eV) or photon wavelength λ   (Å). It is obtained with a parameter ϵD = 0.3. Too large values of ϵD result in incorrect profiles, as shown by the following curves: cyan (ϵD = 11), green (ϵD = 6), orange (ϵD = 5), blue (ϵD = 3), grey (ϵD = 2). The bottom right panel displays the corresponding transmission spectrum. For comparison, the black thick curve portrays the spectral line that would obtain if there was no velocity gradient (V2 = V1). See Sect. 8.3 for detailed explanations.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.