IRIS: a generic threedimensional radiative transfer code
^{1}
LERMA, Observatoire de Paris, CNRS et UPMC,
5 place J. Janssen,
92195
Meudon,
France
email: laurent.ibgui@obspm.fr; chantal.stehle@obspm.fr
^{2}
Steward Observatory, University of Arizona,
933 North Cherry
Avenue, Tucson,
AZ
85721,
USA
email: hubeny@as.arizona.edu
^{3}
Laboratoire J.L. Lagrange, Université de NiceSophia Antipolis,
CNRS, Observatoire de la Côte d’Azur, BP 4229, 06304
Nice,
France
email: thierry.lanz@oca.eu
Received:
30
September
2012
Accepted:
14
November
2012
Context. For most astronomical objects, radiation is the only probe of their physical properties. Therefore, it is important to have the most elaborate theoretical tool to interpret observed spectra or images, thus providing invaluable information to build theoretical models of the physical nature, the structure, and the evolution of the studied objects.
Aims. We present IRIS, a new generic threedimensional (3D) spectral radiative transfer code that generates synthetic spectra, or images. It can be used as a diagnostic tool for comparison with astrophysical observations or laboratory astrophysics experiments.
Methods. We have developed a 3D shortcharacteristic solver that works with a 3D nonuniform Cartesian grid. We have implemented a piecewise cubic, locally monotonic, interpolation technique that dramatically reduces the numerical diffusion effect. The code takes into account the velocity gradient effect resulting in gradual Doppler shifts of photon frequencies and subsequent alterations of spectral line profiles. It can also handle periodic boundary conditions. This first version of the code assumes local thermodynamic equilibrium (LTE) and no scattering. The opacities and source functions are specified by the user. In the near future, the capabilities of IRIS will be extended to allow for nonLTE and scattering modeling.
Results. IRIS has been validated through a number of tests. We provide the results for the most relevant ones, in particular a searchlight beam test, a comparison with a 1D planeparallel model, and a test of the velocity gradient effect.
Conclusions. IRIS is a generic code to address a wide variety of astrophysical issues applied to different objects or structures, such as accretion shocks, jets in young stellar objects, stellar atmospheres, exoplanet atmospheres, accretion disks, rotating stellar winds, cosmological structures. It can also be applied to model laboratory astrophysics experiments, such as radiative shocks produced with high power lasers.
Key words: methods: numerical / radiative transfer
© ESO, 2013
1. Introduction
Essentially all objects in the Universe have a complicated, threedimensional, 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 planeparallel, horizontallyhomogeneous 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 onedimensional problem.
In the last two decades, there was much activity devoted to extend the traditional modeling techniques to treat structures that are truly threedimensional (3D). In fact, detailed 3D hydrodynamic simulations have now become essentially routine. A nonexhaustive 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 longrange 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 timedependent, nonLTE 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 CopenhagenOslo Stagger Codes (Nordlund & Galsgaard 1995; Galsgaard & Nordlund 1996; Hansteen 2004; Hansteen et al. 2007), CO^{5}BOLD (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 frequencyintegrated 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 postprocessing 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 shortcharacteristics 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 threedimensional radiative transfer solver, named IRIS. Analogously to most of the above mentioned techniques, it uses the shortcharacteristics 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 (nonrelativistic) 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 socalled nonLTE (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 shortcharacteristics 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 timeindependent 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 freeflight 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 socalled “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 threedimensional 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 twodimensional 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 shortcharacteristics” method (Olson & Kunasz 1987). All these methods are called the “longcharacteristics” 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 shortcharacteristics 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 nonLTE 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 longcharacteristics scheme is useful for some purposes, such as allowing for an efficient parallel scheme within domain decomposition (Heinemann et al. 2006), we chose the shortcharacteristics scheme. The classical setup of the longcharacteristics 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 N^{3} grid points (and also ≈ N^{3} elementary cells). A longcharacteristics scheme would consider N^{3} 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(N^{4}) operations. In contrast, the shortcharacteristics 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(N^{3}) operations.
To avoid confusion, we mention that when using the longcharacteristics scheme for obtaining an emergent intensity only, it requires also O(N^{3}) operations, because we have N^{2} 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 shortcharacteristics scheme. In some cases, the use of longcharacteristics scheme may even be more computationally advantageous than the shortcharacteristics, and therefore some codes (e.g., ASSeT, SCATE) use the shortcharacteristics scheme when dealing iteratively with scattering, and use the longcharacteristics 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(N^{3}). In astrophysics, the vast majority of approaches is based on the shortcharacteristics 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 longcharacteristics method as well, but in view of our intended future development of the solver to treat scattering and nonLTE situations, the shortcharacteristic scheme is clearly the appropriate choice.
3. 3D shortcharacteristics
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 r_{u} to a current position r_{c}, the integral form of the formal solution of the timeindependent RTE is: (5)where , , n is the unit vector along the straight line (r_{u},r_{c}), τ is the optical depth from r_{u} to an intermediate position r_{u} + sn, with , and τ_{uc} is the optical depth from r_{u} to r_{c}.
We consider a threedimensional Cartesian grid (see Fig. 1) defined in the (O,x,y,z) coordinate system with the unit vectors (e_{x},e_{y},e_{z}). 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: Δx_{cell}(x), Δy_{cell}(y), and Δz_{cell}(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 e_{z}, and by the azimuthal angle ϕ, between e_{x} and the projection of n on the xy 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 I_{c} in a current point M_{c} in direction n. M_{c} is defined by the intersection of the cells marked with indices i, i+1 in xdirection, j, j+1 in ydirection, and k, k+1 in zdirection (see Fig. 1). This is accomplished by using Eq. (5). Following the propagation of the ray that goes through M_{c} in direction n, we define a shortcharacteristic by the line joining the intersection of the ray with the upwind cell, M_{u}, to the current point M_{c} (the upwind cell is defined by the indices (i,j,k) in Fig. 1). I_{c} can be determined if we know the following quantities: the upwind specific intensity I_{u}, 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 shortcharacteristic from M_{u} to M_{c}. 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 shortcharacteristic, typically as loworder polynomials; we choose third degree polynomials. To do so, we need to know (S,χ) in the upwind end point M_{u}, (S_{u},χ_{u}), and in the downwind end point M_{d}, (S_{d},χ_{d}). M_{d} is the intersection of the shortcharacteristic 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 (I_{u},S_{u},χ_{u}) and the downwind quantities (S_{d},χ_{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.
Fig. 1 Shortcharacteristics 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 (e_{x},e_{y},e_{z}). The specific intensity is calculated in the current point M_{c}, for a radiation propagating from the upwind end point M_{u} to the downwind end point M_{d}. The direction n of the ray is defined by the polar angle θ and the azimuthal angle ϕ. The shortcharacteristic is defined by the line joining M_{u} and M_{c}. The transfer equation is solved in its integral form along this shortcharacteristic. The cells around point M_{c} are marked with the following indices: i, i+1 in xdirection, j, j+1 in ydirection, and k, k+1 in zdirection. 

Open with DEXTER 
3.2. Intersections of a shortcharacteristic with the neighboring cell faces
The coordinates of the intersections of a shortcharacteristic with the neighboring cell faces, i.e., the coordinates of the upwind end point M_{u} and the downwind end point M_{d}, are easily determined in the local coordinate system defined by (M_{c},e_{x},e_{y},e_{z}) (see Fig. 1).
Figure 2 displays the upwind cell (top panel) and the downwind cell (bottom panel), as they are defined for the shortcharacteristic 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 f_{x} (respectively f_{y}, f_{z}) as the generic name of a cell face that is perpendicular to xaxis (respectively yaxis, zaxis), and that may be intersected by a shortcharacteristic, i.e., to which M_{u} or M_{d} may belong. In our example, f_{x} in the upwind cell is specified by x = −Δx_{i}, while f_{x} in the downwind cell is specified by x = Δx_{i + 1}, both in the local coordinate system. In the same vein, f_{y} is y = −Δy_{j} in the upwind cell, and y = Δy_{j + 1} in the downwind cell. Finally, f_{z} is z = −Δz_{k} in the upwind cell, and z = Δz_{k + 1} in the downwind cell.
Fig. 2 Upwind cell (top) and downwind cell (bottom), for the shortcharacteristic M_{u}M_{c} 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 (M_{c},e_{x},e_{y},e_{z}). The upwind end point M_{u} belongs to the cell face f_{z}:z = −Δz_{k}. However, depending on (θ,ϕ) in this octant, M_{u} may belong to f_{x}:x = −Δx_{i}, or to f_{y}:y = −Δy_{j}. In the general case, Δx_{i} ≠ Δy_{j}, Δx_{i} ≠ Δz_{k}, Δy_{j} ≠ Δz_{k}, though equalities are possible. The downwind end point M_{d} belongs to f_{z}:z = Δz_{k + 1}. However, depending on (θ,ϕ) in this octant, M_{d} may belong to f_{x}:x = Δx_{i + 1} or to f_{y}:y = Δy_{j + 1}. See Sect. 3.2 for more details. 

Open with DEXTER 
In the example displayed in Figs. 1 and 2, M_{u} belongs to the face f_{z}:z = −Δz_{k}, and M_{d} belongs to the face f_{z}:z = Δz_{k + 1}. However, depending on the values of (θ,ϕ), M_{u} may belong to f_{y}:y = −Δy_{j}, or to f_{x}:x = −Δx_{i}, and M_{d} may belong to f_{y}:y = Δy_{j + 1}, or to f_{x}:x = Δx_{i + 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 Δx_{u} and Δx_{d} as the generic names for the size in xdirection of, respectively, the upwind and the downwind cell, for a given current point M_{c} and a given direction (θ,ϕ). Similarly, let us introduce the analogous notations Δy_{u},Δy_{d} and Δz_{u},Δz_{d} 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 Δx_{u},Δy_{u},Δz_{u}, for any current grid point M_{c}. 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 M_{u(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.
Correspondence between octant, directional angles (θ,ϕ), upwind cell indices, sizes of the upwind cell faces in each direction, and (α_{xu},α_{yu},α_{zu}) factors.
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 shortcharacteristics method requires to define laws of variation of physical quantities within faces of the grid and along the shortcharacteristics. 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 shortcharacteristic’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 (x_{i})_{i = 1,n}, with . In other words, w_{i} = w(x_{i}) for i = 1,2,...,n are given. Our objective is to find an interpolation function that goes through all the data points (x_{i},w_{i}) 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, x_{i} and x_{i + 1} (see Fig. 3). Within the interval [x_{i},x_{i + 1}] , a cubic Hermite polynomial is defined as (see Eq. (2) of Auer 2003): (8)where (9)
Fig. 3 Piecewise interpolation. w(x) is the function to be interpolated between the end points (or nodes) x_{i} and x_{i + 1}. We assume that we know w and its derivatives at the end points: w_{i} = w(x_{i}), w_{i + 1} = w(x_{i + 1}), , . The local interpolation function, H_{i}(x), is a cubic Hermite polynomial. It is a third degree polynomial defined between x_{i} and x_{i + 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 (H_{i}(x))_{i = 1,n − 1} is monotonic in the interpolation interval [x_{i},x_{i + 1}] . More precisely, the sign of H_{i}(x), for x ∈] x_{i},x_{i + 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 [x_{i},x_{i + 1}] is equal the local function H_{i}(x). H(x) is continuous. Moreover, it is easy to demonstrate that its derivative H′(x) is also continuous throughout the whole interval [x_{1},x_{n}] , specifically at each inner node (x_{i})_{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 w_{i}, w_{i + 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 B_{i}(x) the local Bézier polynomial defined in [x_{i},x_{i + 1}] . B_{i}(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 B_{i}(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 B_{i}(x). In this case, the continuity of the derivative B′(x) of the piecewise interpolant B(x), which is defined over [x_{1},x_{n}] 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 (x_{i})_{i = 2,n − 1} from the unique parabola passing through the points (x_{i − 1},w_{i − 1};x_{i},w_{i};x_{i + 1},w_{i + 1}), if this parabola is monotonic in [x_{i − 1},x_{i + 1}] . If not, the node derivative is defined as the common slope shared by two parabolas, which are locally monotonic in [x_{i − 1},x_{i}] and [x_{i},x_{i + 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 shortcharacteristics 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 zaxis, in the sense that we are sweeping the 3D grid gradually, zplane by zplane, in the direction of increasing (or decreasing) z. In an astrophysical context, zaxis 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 N_{x} cells in xdirection, N_{y} cells in ydirection, and N_{z} cells in zdirection, the positions of the grid nodes in each direction may be noted as (x_{0},x_{1},...,x_{Nx}), (y_{0},y_{1},...,y_{Ny}), and (z_{0},z_{1},...,z_{Nz}). The sizes of the cells in each direction are then noted (Δx_{1},...,Δx_{Nx}), (Δy_{1},...,Δy_{Ny}), and (Δz_{1},...,Δz_{Nz}), 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 = x_{0}, y = y_{0}, z = z_{0}.
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 (I_{u},S_{u},χ_{u}) at the intersection of the shortcharacteristic with the upwind face, and (S_{d},χ_{d}) at the intersection of the shortcharacteristic with the downwind face.
Fig. 4 Examples of cubic monotonic interpolations in cell faces intersected by a shortcharacteristic, 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 (S_{u},χ_{u}), (S_{d},χ_{d}). But, for the upwind specific intensity I_{u}, it occurs only when horizontal faces f_{z}:z = ± Δz_{k} or boundary faces are intersected. The second configuration is displayed in the bottom panel. It impacts the determination of I_{u} when vertical nonboundary faces, f_{x}:x = ± Δx_{i} or f_{y}:y = ± Δy_{j}, 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 threedimensional grid. And, by definition of a formal solution, the source function S is known at each grid point. Then, S_{u}, χ_{u}, (respectively S_{d}, χ_{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 f_{z}:z = −Δz_{k} 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. M_{u} is the upwind end point. is the projection of the current point M_{c} on f_{z}, and coincides with a grid point. Quantities defined in the grid points are interpolated in ydirection at the y_{u} position of M_{u}. Then, the interpolated values, in positions A_{i − 2}, A_{i − 1}, A_{i}, A_{i + 1}, are themselves interpolated in xdirection at the x_{u} position of M_{u}, which provides the upwind quantity S_{u}, χ_{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 zplane by zplane, as described above in the first paragraph of this section. Therefore, for each new zplane to be processed, the specific intensity is known at each grid point of the upwind zplane. Accordingly, the procedure described above for S and χ works if M_{u} belongs to a horizontal face. However, the specific intensity is not known at all grid points of a preceding vertical plane (f_{x}:x = ± Δx_{i}, or f_{y}:y = ± Δy_{j}, 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 = z_{k} is the current zplane, z = z_{k − 1} is the upwind one. M_{u} belongs to the vertical face identified by indices j and k, and by f_{x}:x = −Δx_{i}. We still consider a ray propagating in the first octant (increasing x, y, z). The specific intensity is fully known in the preceding z = z_{k − 2} and z = z_{k − 1} planes, where we can perform interpolations in ydirection at the y_{u} position of M_{u}. This provides the values of the specific intensity in points A_{k − 2} and A_{k − 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 = z_{k} plane. Therefore, the specific intensity is known in the three consecutive grid points ((y_{j − 2},z_{k}), (y_{j − 1},z_{k}), (y_{j},z_{k})), but not in (y_{j + 1},z_{k}). It is also not known in the z = z_{k + 1} plane (except in y = y_{j − 2} if this yplane is a boundary). Consequently, we can determine the specific intensity in point A_{k} in z = z_{k} plane, with a cubic Hermite interpolation in ydirection, using the values in ((y_{j − 2},z_{k}), (y_{j − 1},z_{k}), (y_{j},z_{k})). Finally, the interpolated values in A_{k − 2}, A_{k − 1}, A_{k} are themselves interpolated in zdirection at the z_{u} position of M_{u}, which provides the upwind specific intensity I_{u}. 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 ydirection and A_{k} in zdirection, 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 xdirection, for example, and for a propagation in increasing x, ghost nodes are required for determining the upwind quantities at x_{1}, and for determining the downwind quantities at x_{Nx − 1} and x_{Nx}. The ghost nodes (x_{2},x_{1}) and (x_{Nx + 1},x_{Nx + 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 x_{1}, x_{Nx − 1} and x_{Nx}, 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 shortcharacteristic
In the preceding subsection, we have described our procedure to determine the upwind quantities (I_{u},S_{u},χ_{u}), and the downwind quantities (S_{d},χ_{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 shortcharacteristics, of the source function and the opacity. We use a monotonic cubic Hermite polynomial between the upwind endpoint M_{u} and the current grid point M_{c} (see Fig. 1). Since we know the values of S and χ in three points, M_{u}, M_{c}, and the downwind endpoint M_{d}, the derivatives of these quantities at M_{u}, and , are provided by the slopes defined by Eq. (12a), while the derivatives at M_{c}, 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 M_{u} to M_{c}. We obtain: (15)where s_{uc} is the path length from M_{u} to M_{c}.
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, I_{c} → I_{u}, and the optically thick or even opaque path for which τ_{uc} tends to infinity, and, therefore, I_{c} → S_{c}.
5.3. The velocity gradient effect (Doppler shift)
We consider a medium with a nonzero 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 shortcharacteristics method. If the velocity difference between the current point M_{c} and the upwind end point M_{u} (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 M_{u}, 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 M_{c}, 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 M_{u} to M_{c}. In fact, the line goes through its maximum value at this frequency, somewhere between M_{u} and M_{c}, 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 shortcharacteristic 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 M_{u}, the calculation frequency can be located close to the center of one line. And, at M_{c}, 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 M_{u} and M_{c} 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 shortcharacteristic 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 shortcharacteristic. 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 V_{th} 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 angleaveraged 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 (e_{x},e_{y},e_{z}) 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, F_{x}, F_{y}, F_{z}. Since the radiation pressure tensor is symmetric (Mihalas 1978; Mihalas & Mihalas 1984), it is defined by six components, P_{xx}, P_{yy}, P_{zz}, P_{xy}, P_{xz}, P_{yz}, 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 n_{m} = (n_{mx},n_{my},n_{mz}) and the corresponding weights w_{m}, 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 (e_{x},e_{y},e_{z});

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 e_{i} stands for e_{x}, e_{y}, or e_{z}, 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 threedimensional medium with an infinite extension in the horizontal plane (x,y), and a finite extension along the vertical zaxis between its lower boundary z_{0} and its upper boundary z_{Nz}. We assume that this medium has a double periodicity, one in xdirection, and one in ydirection. The boundary conditions are known at z_{0} and z_{Nz}. For example, we may consider a nonirradiated stellar atmosphere with no incoming radiation at the outer surface z_{Nz} and a black body radiation at its inner surface z_{0}. The computational grid ranges from z_{0} to z_{Nz} in vertical direction, from x_{0} to x_{Nx}, and from y_{0} to y_{Ny} in the horizontal plane, so that (x_{Nx} − x_{0}) defines one xperiod and (y_{Ny} − y_{0}) defines one yperiod. Now, the 3D shortcharacteristics 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 shortcharacteristic intersects a lateral boundary, we prolong this characteristic, which becomes a longcharacteristic, 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 zlayer, 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 M_{c} in the plane x = x_{1}, which is adjacent to the boundary plane x = x_{0}. The upwind end point M_{u} of the shortcharacteristic (M_{u},M_{c},M_{d}) 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 shortcharacteristic , the upwind end point belongs to the horizontal boundary plane x = x_{0}, in which the specific intensity is not explicitly defined. The shortcharacteristic is then lengthened, and becomes a longcharacteristic, down to its first intersection, , with a horizontal face. The longcharacteristic 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.
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 xdirection and in ydirection. In xdirection, one period is defined from x_{0} to x_{Nx}. The plane x = x_{0} is a boundary plane. The cell defined by indices (− 1,j,k), drawn with thick lines, is identical to the cell (N_{x},j,k) (not drawn). The upwind specific intensity is known at M_{u} for the shortcharacteristic (M_{u},M_{c},M_{d}), since M_{u} belongs to a horizontal face. On the other hand, the upwind specific intensity is not known at for the shortcharacteristic , 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 zplane 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 = x_{0} and y = y_{0}, therefore the first grid point in which we calculate the specific intensity is (x_{1},y_{1}). Then, the specific intensity is calculated progressively in all the grid points in the zplane along xdirection and ydirection. 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 zplane, provided that the grid points where the specific intensity is calculated span one xperiod and one yperiod. We mark this starting point as (x_{i0},y_{j0}). The indices i_{0} and j_{0} are defined so as to minimize the number of cells that are intersected by the longcharacteristics, thus saving computing time. This way, i_{0} is chosen such that the corresponding upwind cell has the largest size in xdirection, Δx_{i0} = Δx_{max}. Similarly, we choose j_{0} such that Δy_{j0} = Δj_{max}. Then, starting from (x_{i0},y_{j0}), we calculate the specific intensity with the longcharacteristics method at y_{j0}, for the grid points along xdirection that span one xperiod. In the same vein, we use the longcharacteristics method at x_{i0}, for the grid points along ydirection that span one yperiod. Once these two first rows have been treated, we resume the shortcharacteristics method for all the next rows of the current zlayer. Note that here we do not introduce ghost nodes (see Sect. 5.1), since the medium has an infinite extension in xdirection and ydirection.
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 zplane. Both span one xperiod and one yperiod. The original domain is depicted by the rectangle with solid thick lines. The positions of the grid nodes in the horizontal zplane were introduced in Sect. 5.1, (x_{0}, x_{1}, ..., x_{Nx}) and (y_{0}, y_{1}, ..., y_{Ny}). The vertical boundary planes are x = x_{0} and y = y_{0}. The indices i_{0} and j_{0} 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 (x_{i0}, x_{i0 + 1}, ..., x_{i0 + Nx}) and (y_{j0}, y_{j0 + 1}, ..., y_{j0 + 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 longcharacteristics method along xdirection from (x_{i0},y_{j0}) to (x_{i0 + Nx},y_{j0}), and along ydirection from (x_{i0},y_{j0}) to (x_{i0},y_{j0 + Ny}). Then, we resume our usual shortcharacteristics method for the next rows of the current zlayer. 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 longcharacteristics may consist merely of shortcharacteristics. This is the case for a given grid point (x_{i, i ∈ [[i0,Nx]]},y_{j0}), if the corresponding Δx_{i} and Δy_{max} are large enough, so that the related upwind end point intersects the horizontal upwind face. Similarly, this is the case for (x_{i0},y_{j, j ∈ [[j0,Ny]]}), if Δy_{j} and Δx_{max} 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).
Fig. 6 Original domain (solid thick lines) and calculation domain (dotted thick lines) at a given zplane, 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 xdirection and one period in ydirection. The nodes of the original domain are (x_{0},x_{1},...,x_{Nx}) and (y_{0},y_{1},...,y_{Ny}). The nodes of the calculation domain, where the specific intensity is computed, are (x_{i0},x_{i0 + 1},...,x_{i0 + Nx}) and (y_{j0},y_{j0 + 1},...,y_{j0 + Ny}). The indices i_{0} and j_{0} are chosen such that Δx_{i0} = Δx_{max} and Δy_{j0} = Δj_{max}. 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 welltested onedimensional scheme, and the test of the velocity gradient effect applied on one given spectral line.
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 xaxis, 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 shortcharacteristics 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, I_{c} = I_{u}. 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 hardedged beam that propagates throughout a threedimensional 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 100^{3} 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 30^{2} 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 hardedges 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, [(x_{1},y_{1});(x_{2},y_{1});(x_{1},y_{2});(x_{2},y_{2})] , 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 xaxis (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 xysurface integral of the specific intensity at z = 0 equals the the xysurface 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 shortcharacteristics 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 planeparallel models
It is important to check whether IRIS can reproduce the results provided by onedimensional radiative transfer schemes, when the medium is made of homogeneous planeparallel layers, i.e., in the case of a 1D planeparallel 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 zplanes, 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 planeparallel layers. The hydrodynamics structure, provided by Matthias González^{1}, is obtained with the threedimensional 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 threedimensional 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 zplane, all the state parameters (temperature, density, velocity) are independent of x and y. Two postprocessing calculations were performed: one with IRIS applied on the 3D grid, and one with a welltested 1D shortcharacteristics solver applied on a 1D grid that is extracted from the 3D one. Both were done at seventeen frequencies that range from hν = 2 eV (λ = 620 nm) to hν = 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, hν = 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 zplane 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, F_{x}, F_{y}, F_{z}, in the bottom left panel, and the components of the radiation pressure tensor P_{xx}, P_{yy}, P_{zz}, P_{xy}, P_{xz}, P_{yz}, 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 planeparallel structures. The radiative flux is zero in x and y direction, F_{x} = 0, F_{y} = 0, and it depends only on z in zdirection, F_{z}(z,ν). The nondiagonal components of the radiation pressure tensor are all zero: P_{xy} = 0, P_{xz} = 0, P_{yz} = 0. The other three pressure components depend only on z, P_{xx}(z,ν), P_{yy}(z,ν), P_{zz}(z,ν), and we have, for any z and ν, P_{xx}(z,ν) = P_{yy}(z,ν).
These tests demonstrate that IRIS is capable of reproducing 1D planeparallel simulations, and, thereby, of handling periodic boundary conditions.
Fig. 8 Comparison with 1D planeparallel models. The figures shows the monochromatic radiative quantities, at hν = 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 planeparallel 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 (P_{zz}:blue,square;P_{xx} = P_{yy}: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 shortcharacteristic 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 quasianalytically. 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 zdirection. In our example, the material moves toward an observer. The layer is numerically modeled with a single cell between two grid points, z_{1}, the farther position from the observer, and z_{2}, the closer position to the observer. The wind velocities at grid points are V_{1} > 0 and V_{2} > V_{1}. 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 z_{1}. 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 z_{1} to z_{2}: (34)where the absorption coefficient κ(ν,z) is written as (35)where ν_{ctr}(z), defined in Sect. 5.3 by Eq. (16), is the shifted linecenter 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 shortcharacteristic, τ_{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: V_{1} = 2000 km s^{1}, V_{2} = 2250 km s^{1}, and a thermal velocity (constant in our uniform layer) of V_{th} = 25 km s^{1}. With such a choice, any line is Doppler shifted by an amount of 10 times its Doppler width between position z_{1} and position z_{2}, i.e., (ν_{ctr2} − ν_{ctr1})/Δν_{D} = (V_{2} − V_{1})/V_{th} = 10. We also adopt a temperature T = 10^{5} K, and a layer size z_{2} − z_{1} = 10^{9} cm.
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 z_{1} and z_{2}, with wind velocities V_{1} > 0 (the material approaches the observer) and V_{2} > V_{1}. 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 z_{2} toward the observer, as computed by IRIS. This emission spectrum is plotted versus photon energy hν (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 (V_{2} = V_{1}). 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 z_{1}, and at ν_{ctr2} ≈ 80.60 eV at z_{2}. 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 z_{2} toward the observer, as calculated by IRIS. This emission spectrum is displayed versus photon energy hν (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 V_{2} = V_{1}. 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 shortcharacteristic between z_{1} and z_{2} 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 shortcharacteristic. The cyan curve shows the profile obtained with ϵ_{D} = 11, for which we checked that the shortcharacteristic 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 shortcharacteristic 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 threedimensional 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 freeflight 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 postprocessing 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 nonstationary flow. We consider nonrelativistic 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 shortcharacteristics 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 shortcharacteristics 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 shortcharacteristics.
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 (nonrelativistic) nonmonotonic 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 XUV 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. Postprocessing 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.
Acknowledgments
We are grateful to Matthias González (AIM, CEA/DSM/IRFU, CNRS, Université Paris Diderot, 91191 GifsurYvette, 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 planeparallel 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 08BLAN026307. The work of I.H. was supported in part by NSF grant AST0807496. 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
 Auer, L. 1976, J. Quant. Spectroc. Radiat. Transf., 16, 931 [NASA ADS] [CrossRef] [Google Scholar]
 Auer, L. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 3 [Google Scholar]
 Auer, L. H., & Paletou, F. 1994, A&A, 285, 675 [NASA ADS] [Google Scholar]
 Auer, L., Bendicho, P. F., & Trujillo Bueno, J. 1994, A&A, 292, 599 [NASA ADS] [Google Scholar]
 Baron, E., & Hauschildt, P. H. 2007, A&A, 468, 255 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Botnen, A. 1997, Master’s thesis, Inst. Theor. Astrophys. Oslo [Google Scholar]
 Brodlie, K. W. 1980, in Mathematical Methods in Computer Graphics and Design, ed. K. W. Brodlie (London: Academic Press), 1 [Google Scholar]
 Bruls, J. H. M. J., Vollmöller, P., & Schüssler, M. 1999, A&A, 348, 233 [NASA ADS] [Google Scholar]
 Bryan, G. L., & Norman, M. L. 1997 [arXiv:astroph/9710187] [Google Scholar]
 Carlson, B. 1963, in Methods in Computational Physics, eds. B. Alder, S. Fernbach, & M. Rotenberg (Academic Press), 1, 1 [Google Scholar]
 Carlsson, M. 1986, Uppsala Astronomical Observatory Reports, 33 [Google Scholar]
 Carlsson, M. 2008, Phys. Scr. Vol. T, 133, 014012 [NASA ADS] [CrossRef] [Google Scholar]
 Carlsson, M. 2009, Mem. Soc. Astron. Italiana, 80, 606 [NASA ADS] [Google Scholar]
 Castor, J. I., Dykema, P. G., & Klein, R. I. 1992, ApJ, 387, 561 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1950, Radiat. Transf. (London: Oxford University Press) [Google Scholar]
 Cunningham, A. J., Frank, A., Varnière, P., Mitran, S., & Jones, T. W. 2009, ApJS, 182, 519 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, S. W., Stone, J. M., & Jiang, Y.F. 2012, ApJS, 199, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Fabiani Bendicho, P. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 419 [Google Scholar]
 Fabiani Bendicho, P., & Trujillo Bueno, J. 1999, in Solar Polarization, eds. K. N. Nagendra, & J. O. Stenflo, Astrophys. Space Sci. Lib., 243, 219 [Google Scholar]
 Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. 1997, A&A, 324, 161 [Google Scholar]
 Feautrier, P. 1964, C. R. Acad. Sci., 258, 3189 [NASA ADS] [Google Scholar]
 Freytag, B., Steffen, M., & Dorch, B. 2002, Astron. Nachr., 323, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Fritsch, F. N., & Butland, J. 1984, SIAM J. Sci. Stat. Comput., 5, 300 [CrossRef] [MathSciNet] [Google Scholar]
 Galsgaard, K., & Nordlund, Å. 1996, J. Geophys. Res., 101, 13445 [NASA ADS] [CrossRef] [Google Scholar]
 González, M., Audit, E., & Huynh, P. 2007, A&A, 464, 429 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gudiksen, B. V., Carlsson, M., Hansteen, V. H., et al. 2011, A&A, 531, A154 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hansteen, V. H. 2004, in MultiWavelength Investigations of Solar Activity, eds. A. V. Stepanov, E. E. Benevolenskaya, & A. G. Kosovichev, IAU Symp., 223, 385 [Google Scholar]
 Hansteen, V. H., Carlsson, M., & Gudiksen, B. 2007, in The Physics of Chromospheric Plasmas, eds. P. Heinzel, I. Dorotovič, & R. J. Rutten, ASP Conf. Ser., 368, 107 [Google Scholar]
 Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 2008, A&A, 490, 873 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hayek, W., Asplund, M., Collet, R., & Nordlund, Å. 2011, A&A, 529, A158 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heinemann, T., Dobler, W., Nordlund, Å., & Brandenburg, A. 2006, A&A, 448, 731 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hubeny, I. 1988, Comput. Phys. Commun., 52, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Lanz, T. 1995, ApJ, 439, 875 [NASA ADS] [CrossRef] [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., et al. 2012a, in SF2A2012: Proc. Ann. Meet. French Soc. Astron. Astrophys., eds. S. Boissier, P. de Laverny, N. Nardetto, R. Samadi, D. VallsGabaud, & H. Wozniak [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., et al. 2012b, in Astronum 2012, ASP Conf. Ser. [Google Scholar]
 Koesterke, L. 2009, in AIP Conf. Ser. 1171, eds. I. Hubeny, J. M. Stone, K. MacGregor, & K. Werner, 73 [Google Scholar]
 Koesterke, L., Allen de Prieto, C., & Lambert, D. L. 2008, ApJ, 680, 764 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectroc. Radiat. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P. B., & Olson, G. L. 1988, J. Quant. Spectroc. Radiat. Transf., 39, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lathrop, K. D., & Carlson, B. G. 1965, LASL3186 [Google Scholar]
 Leenaarts, J., & Carlsson, M. 2009, in The Second Hinode Science Meeting: Beyond DiscoveryToward Understanding, eds. B. Lites, M. Cheung, T. Magara, J. Mariska, & K. Reeves, ASP Conf. Ser., 415, 87 [Google Scholar]
 Leenaarts, J., Carlsson, M., Hansteen, V., & Rouppe van der Voort, L. 2009, ApJ, 694, L128 [NASA ADS] [CrossRef] [Google Scholar]
 Michaut, C., Stehlé, C., Leygnac, S., Lanz, T., & Boireau, L. 2004, Eur. Phys. J. D, 28, 381 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of radiation hydrodynamics [Google Scholar]
 Mihalas, D., Auer, L. H., & Mihalas, B. R. 1978, ApJ, 220, 1001 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, A. 1982, A&A, 107, 1 [NASA ADS] [Google Scholar]
 Nordlund, Å., & Galsgaard, K. 1995, A 3D MHD Code for Parallel Computers, Tech. Rep., Astronomical Observatory, Copenhagen University [Google Scholar]
 Norman, M. L., & Bryan, G. L. 1999, in Astrophys. Num. Astrophys., eds. S. M. Miyama, K. Tomisaka, & T. Hanawa, Space Sci. Lib., 240, 19 [Google Scholar]
 Norman, M. L., Bryan, G. L., Harkness, R., et al. 2007 [arXiv:0705.1556] [Google Scholar]
 Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spectroc. Radiat. Transfer, 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
 O’Shea, B. W., Bryan, G., Bordner, J., et al. 2004 [arXiv:astroph/0403044] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in Fortran. The art of scientific computing [Google Scholar]
 Steffen, M. 1990, A&A, 239, 443 [NASA ADS] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992a, ApJS, 80, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992b, ApJS, 80, 791 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Mihalas, D., & Norman, M. L. 1992, ApJS, 80, 819 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trujillo Bueno, J., & FabianiBendicho, P. 1995, ApJ, 455, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [Google Scholar]
 Uitenbroek, H. 2006, ApJ, 639, 516 [NASA ADS] [CrossRef] [Google Scholar]
 van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066 [NASA ADS] [CrossRef] [Google Scholar]
 Vögler, A., Shelyag, S., Schüssler, M., et al. 2005, A&A, 429, 335 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
Correspondence between octant, directional angles (θ,ϕ), upwind cell indices, sizes of the upwind cell faces in each direction, and (α_{xu},α_{yu},α_{zu}) factors.
All Figures
Fig. 1 Shortcharacteristics 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 (e_{x},e_{y},e_{z}). The specific intensity is calculated in the current point M_{c}, for a radiation propagating from the upwind end point M_{u} to the downwind end point M_{d}. The direction n of the ray is defined by the polar angle θ and the azimuthal angle ϕ. The shortcharacteristic is defined by the line joining M_{u} and M_{c}. The transfer equation is solved in its integral form along this shortcharacteristic. The cells around point M_{c} are marked with the following indices: i, i+1 in xdirection, j, j+1 in ydirection, and k, k+1 in zdirection. 

Open with DEXTER  
In the text 
Fig. 2 Upwind cell (top) and downwind cell (bottom), for the shortcharacteristic M_{u}M_{c} 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 (M_{c},e_{x},e_{y},e_{z}). The upwind end point M_{u} belongs to the cell face f_{z}:z = −Δz_{k}. However, depending on (θ,ϕ) in this octant, M_{u} may belong to f_{x}:x = −Δx_{i}, or to f_{y}:y = −Δy_{j}. In the general case, Δx_{i} ≠ Δy_{j}, Δx_{i} ≠ Δz_{k}, Δy_{j} ≠ Δz_{k}, though equalities are possible. The downwind end point M_{d} belongs to f_{z}:z = Δz_{k + 1}. However, depending on (θ,ϕ) in this octant, M_{d} may belong to f_{x}:x = Δx_{i + 1} or to f_{y}:y = Δy_{j + 1}. See Sect. 3.2 for more details. 

Open with DEXTER  
In the text 
Fig. 3 Piecewise interpolation. w(x) is the function to be interpolated between the end points (or nodes) x_{i} and x_{i + 1}. We assume that we know w and its derivatives at the end points: w_{i} = w(x_{i}), w_{i + 1} = w(x_{i + 1}), , . The local interpolation function, H_{i}(x), is a cubic Hermite polynomial. It is a third degree polynomial defined between x_{i} and x_{i + 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 
Fig. 4 Examples of cubic monotonic interpolations in cell faces intersected by a shortcharacteristic, 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 (S_{u},χ_{u}), (S_{d},χ_{d}). But, for the upwind specific intensity I_{u}, it occurs only when horizontal faces f_{z}:z = ± Δz_{k} or boundary faces are intersected. The second configuration is displayed in the bottom panel. It impacts the determination of I_{u} when vertical nonboundary faces, f_{x}:x = ± Δx_{i} or f_{y}:y = ± Δy_{j}, 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 
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 xdirection and in ydirection. In xdirection, one period is defined from x_{0} to x_{Nx}. The plane x = x_{0} is a boundary plane. The cell defined by indices (− 1,j,k), drawn with thick lines, is identical to the cell (N_{x},j,k) (not drawn). The upwind specific intensity is known at M_{u} for the shortcharacteristic (M_{u},M_{c},M_{d}), since M_{u} belongs to a horizontal face. On the other hand, the upwind specific intensity is not known at for the shortcharacteristic , 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 
Fig. 6 Original domain (solid thick lines) and calculation domain (dotted thick lines) at a given zplane, 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 xdirection and one period in ydirection. The nodes of the original domain are (x_{0},x_{1},...,x_{Nx}) and (y_{0},y_{1},...,y_{Ny}). The nodes of the calculation domain, where the specific intensity is computed, are (x_{i0},x_{i0 + 1},...,x_{i0 + Nx}) and (y_{j0},y_{j0 + 1},...,y_{j0 + Ny}). The indices i_{0} and j_{0} are chosen such that Δx_{i0} = Δx_{max} and Δy_{j0} = Δj_{max}. See text in Sect. 7 for a discussion. 

Open with DEXTER  
In the text 
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 xaxis, 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 
Fig. 8 Comparison with 1D planeparallel models. The figures shows the monochromatic radiative quantities, at hν = 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 planeparallel 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 (P_{zz}:blue,square;P_{xx} = P_{yy}:red,plus). See Sect. 8.2 for more details. 

Open with DEXTER  
In the text 
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 z_{1} and z_{2}, with wind velocities V_{1} > 0 (the material approaches the observer) and V_{2} > V_{1}. 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 z_{2} toward the observer, as computed by IRIS. This emission spectrum is plotted versus photon energy hν (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 (V_{2} = V_{1}). See Sect. 8.3 for detailed explanations. 

Open with DEXTER  
In the text 