Method

In the following discussion we use notation of Papers I and II. The basic framework and the methods used for the formal solution and the solution of the scattering problem via operator splitting are discussed in detail in Papers I and II and will thus not be repeated here.

In the following we assume (without restriction) that we have periodic boundary conditions in the $x$ and $y$ coordinates, and for the $z$ coordinate that the `bottom' (large optical depth) is at $z=z_{\rm min}$ and the `top' (interface to empty space) is at $z=z_{\rm max}$. The implementation of the periodic boundary conditions within our framework is simple: We use a ``full characteristics'' approach that completely tracks a set of characteristics of the radiative transfer equation from the outer boundary through the computational domain to their exit voxel and takes care that each voxel is hit by at least one characteristic per solid angle. One characteristic is started on each boundary voxel (in this case these are the planes $z=z_{\rm min}$ and $z=z_{\rm max}$) and then tracked until it leaves the other boundary. The direction of a bundle of full characteristics is determined by a set of solid angles $(\theta,\varphi)$ which correspond to a normalized momentum space vector $(p_x,p_y,p_z)$. The periodic boundary conditions are simply implemented as a wrap-around (e.g., passing $x_{\rm max}$ for $p_x>0$ wraps around to $x_{\min}$) and continuing of the characteristic until it leaves at the $z$ boundary. Characteristics with very small $\vert p_z\vert$ would require a large number of wrap-arounds (and eventually would lead to infinitely long characteristics), therefore, we limit the number of wrap-arounds per voxel to a prescribed value, typically around 16 (tests have shown that larger values do not affect the results, values as small as 4 are usable in plane-parallel tests). The code is parallelized as described in Paper II.