Computation of $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$

As demonstrated by Olson et al. (1987) and Olson & Kunasz (1987), the coefficients $\alpha$, $\beta$, and $\gamma$ can be used to construct diagonal and tri-diagonal $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operators for 1D radiation transport problems. In fact, up to the full $\Lambda$ matrix can be constructed by a straightforward extension of the idea (Hauschildt & Baron, 2004; Hauschildt et al., 1994). These non-local $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operators not only lead to excellent convergence rates but they also avoid the problem of false convergence that is inherent in the $\Lambda$ iteration method and can also be an issue for diagonal (purely local) $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operators. Therefore, it is highly desirable to implement a non-local $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ also for the 3D case. The tri-diagonal operator in the 1D case is simply a nearest neighbor $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ that considers the interaction of a point with its two direct neighbors. In the 3D case, the nearest neighbor $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ considers the interaction of a voxel with the (up to) $3^3-1=26$ surrounding voxels (this definition considers a somewhat larger range of voxels that a strictly face-centered view of just 6 nearest neighbors). This means that the non-local $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ requires the storage of 27 (26 surrounding voxels plus local, i.e., diagonal effects) times the total number of voxels $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ elements. This can be reduced, for example if one considers only the voxels that share a face to a total of 7 elements for each voxel. However, this would ignore the effect of characteristics that pass 'diagonally' through a voxel and would therefore lead to a slower convergence rate.

The construction of the $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operator proceeds in the same way as discussed in Hauschildt (1992). In the 3D case, the 'previous' and 'next' voxels along each characteristic must be known so that the contributions can be attributed to the correct voxel. Therefore, we use a data structure that attaches to each voxel its effects on its neighbors. The scheme can be extended trivially to include longer range interactions for better convergence rates (in particular on larger voxel grids), however the memory requirements to simply store $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ ultimately scales like $n^3$ where $n$ is the total number of voxels. The storage requirements can be reduced by, e.g., using $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$'s of different widths for different voxels. Storage requirements are not so much a problem if a domain decomposition parallelization method is used and enough processors are available. Below we will show some results for test cases with larger operators.

We describe here the general procedure of calculating the $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ with arbitrary bandwidth, up to the full $\Lambda$-operator, for the method in spherical symmetry (Hauschildt et al., 1994). The construction of the $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ is described in Olson & Kunasz (1987), so that we here summarize the relevant formulae. In the method of Olson & Kunasz (1987), the elements of the row of $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ are computed by setting the incident intensities (boundary conditions) to zero and setting $S(i_x,i_y,i_z)=1$ for one voxel $(i_x,i_y,i_z)$ and performing a formal solution analytically.

We describe the construction of $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ using the example of a single characteristic. The contributions to the $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ at a voxel $j$ are given by

$$\Lambda_{i,j} = 0 \quad {\rm for\ } i<j-1$$ (22)

$$\Lambda_{j-1,j} = \gamma_{j-1} \quad {\rm for\ } i=j-1$$ (23)

$$\Lambda_{j,j} = \Lambda_{j-1,j}\exp(-\Delta\tau_{j-1}) + \beta^k_{j} \quad{\rm for\ } i=j$$ (24)

$$\Lambda_{j+1,j} = \Lambda_{j,j}\exp(-\Delta\tau_{j}) + \alpha_{j+1} \quad{\rm for\ } i=j+1$$ (25)

$$\Lambda_{i,j} = \Lambda_{i-1,j}\exp(-\Delta\tau_{i-1}) \quad{\rm for\ } j+1 < i$$ (26)

These contributions are computed along a characteristic, here $i$ labels the voxels along the characteristic under consideration. These contributions are integrated over solid angle with the same method (either deterministic or through the Monte-Carlo integration) that is used for the computation of the $J$. For a nearest neighbor $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$, the process of Eq. 26 is stopped with $i=j+1$, otherwise it is continued until the required bandwidth has been reached (or the characteristic has reached an outermost voxel and terminates).