© The Authors 2023

1 Introduction

Radiative transfer calculations are essential in stellar spectroscopy. They are the basic building block of spectral synthesis from model atmospheres, which are used to infer the temperature, chemical composition, magnetic fields, and many other quantities. In optically thick stellar atmospheres, detailed calculations outside local thermodynamical equilibrium (LTE), or non-LTE (NLTE), are typically solved iteratively (e.g. Cannon 1973; Rybicki & Hummer 1991). To ensure convergence of these schemes, we need to accurately evaluate the radiation field at all grid points (known as a formal solution). Radiation is sensitive to changes in optical depth. Numerically, the formal solution is most stable when there are no large jumps in optical depth between the grid points. However, this is often not the case in model atmospheres, whose grids are optimised to solve the hydrostatic or magnetohydrodynamic equations, and not the transport of radiation.

To increase the accuracy in radiative transfer in stellar atmospheres, different strategies are used. The brute-force approach would be to interpolate the model to a much finer grid, but this raises many problems. First, it can dramatically increase the computational effort because we need to integrate over potentially many more points. Second, iterative NLTE schemes need to propagate a change in the solution, which is linearly proportional to the number of depth points (see e.g. Olson et al. 1986; Štěpán et al. 2022). Therefore, it is highly desirable to solve the formal solution without substantially increasing the number of grid points. One approach is to use higher-order interpolation (e.g. Auer 1976; Trujillo Bueno 2003; Holzreuter & Solanki 2012; de la Cruz Rodríguez & Piskunov 2013; Janett et al. 2019), although these approaches still struggle in coarse grids or discontinuous media (Steiner et al. 2016). Another approach, used in some 1D codes (e.g. in the current version of the MULTI code; Carlsson 1986) is to interpolate from the original grid to a grid that resolves the gradients in optical depth and other associated quantities better. Yet another approach is to use multi-grid iterations, where solving the radiation on a coarser grid is used to remove the high-frequency components of the error of fine grid solutions (e.g. Fabiani Bendicho et al. 1997; Štěpán & Trujillo Bueno 2013; Bjørgen & Leenaarts 2017).

With time-dependent 3D model atmospheres being widely used, it becomes necessary to carry out radiative transfer calculations in 3D. Compared to the classical 1D model atmospheres, these calculations are computationally much more demanding and require efficient algorithms. Assuming the same number of grid points N in the three spatial dimensions, together with a number of angles N_Ω and number of wavelengths N_λ, Štěpán et al. (2022) estimated that the solution time of a 3D NLTE problem is about O(N⁴N_ΩN_λ). This is a great challenge for current problems in solar physics, in particular, for spectral lines that form in the chromosphere. Here, NLTE calculations are essential, and so is spatial resolution, given that the chromo-spheric plasma is organised into much finer scales than in the photosphere. This leads to exploding computational costs, needing many iterations to obtain the converged radiation field, and with computational costs often nearing millions of CPU hours for a single simulation snapshot (Sukhorukov & Leenaarts 2017; Bjørgen & Leenaarts 2017).

The aim of this work is to explore a different approach to solving 3D NLTE problems in stellar atmospheres: irregular grids. In 1D, using a grid that is better suited to the radiative transfer problem can improve the convergence of NLTE problems (see e.g. Pereira & Uitenbroek 2015). However, this approach cannot be directly applied to the 3D case in regular grids because we would need an optimised grid for every ray. By relaxing the requirement of a regular grid and placing grid points freely in 3D space, we can optimise the grid to the radiative transfer problem.

Irregular grids are not new in astrophysical 3D radiative transfer problems, in particular, for cosmological simulations and clumpy media (e.g. Ritzerveld & Icke 2006; Paardekooper et al. 2010). The SIMPLEX algorithm of Ritzerveld & Icke (2006) was designed to solve continuum radiative transfer in irregular grids created from Voronoi tessellation (Dirichlet 1850; Voronoi 1908). Brinch & Hogerheijde (2010) used a similar approach in the LIME code, but extended it to include line opacity and to solve NLTE problems. In these algorithms, rays travel only between Delaunay lines, which are the lines that connect the cell centres of Voronoi cells. Camps et al. (2013) developed another approach in which rays are allowed to travel on a straight path through a 3D Voronoi grid. All these previous studies find that using irregular grids can give accurate results (or results that are comparable to existing methods) using fewer grid points. For their problem, Camps et al. (2013) find that 3D radiative transfer through a Voronoi grid is about three times slower than on a regular grid with the same number of cells, but that fewer cells are required to reach an accurate result. This outweighs the slower calculation times. Bruls et al. (1999) developed methods for solving radiative transfer on unstructured grids and applied them to estimate the radiative heating in stellar atmospheres, but their approach is limited to LTE and 2D geometry. To the best of our knowledge, no study has employed irregular grids in 3D NLTE problems in optically thick stellar atmospheres. The motivation for our work is therefore to investigate whether methods like this could be of value in stellar atmospheres, in particular, to speed up the convergence of 3D NLTE schemes.

2 Methods

2.1 NLTE radiative transfer

We aim to solve the time-independent radiative transfer equation, which for a given ray can be written as $$\frac{\text{d}{I}_{\lambda }}{\text{d}{\tau }_{\lambda }}=S_{\lambda }-I_{\lambda },$$(1)

where S_λ is the source function, and τ_λ is the optical depth. The formal solution for the intensity I_λ can be written as $$I_{\lambda }\left({\tau }_{\lambda }\right)=I_{\lambda }\left(0\right){\text{e}}^{-{\tau }_{\lambda }}+{\displaystyle {\int }_0^{{\tau }_{\lambda }}{S}_{\lambda }\left({\tau }_{\lambda }^{\prime }\right){\text{e}}^{-\left({\tau }_{\lambda }-{\tau }_{\lambda }^{\prime }\right)}\text{d}{\tau }_{\lambda }^{\prime }}.$$(2)

In an NLTE problem, S_λ will depend on I_λ via the angle-averaged intensity J_λ. Using the two-level atom simplification, this can be expressed as $$S_{\lambda }\left(J_{\lambda },T\right)=\left(1-{\epsilon }_{\lambda }\right)J_{\lambda }+{\epsilon }_{\lambda }{B}_{\lambda }\left(T\right),$$(3)

where ε_λ is the photon destruction probability, and B_λ(T) is the Planck function. Numerically, the Λ iteration procedure involves iteratively solving Eq. (3) until S_λ is consistent with J_λ. For a given estimate of S_λ, we need to solve Eq. (2) for multiple directions to obtain J_λ, which can be expressed by the Λ operator, $$J_{\lambda }=\text{Λ}\left[S_{\lambda }\right].$$(4)

This basic scheme is not very efficient and struggles in cases with strong scattering (ε_λ < 10⁻³). Accelerated schemes for Λ iteration (ALI) have been developed that are much more efficient (e.g. Cannon 1973; Rybicki & Hummer 1991). However, given the exploratory nature of our analysis, we adopted the simple Λ iteration because it is faster to implement, and we restricted our problem to cases in which ε_λ > 10⁻³ by artificially modifying the model atom.

Efficient solutions of the 3D NLTE problem can involve methods to speed up the formal solution and/or reduce the number of iterations needed to reach a converged solution. In this work, we address the latter.

2.2 Model atom and extinction sources

To test radiative transfer on irregular grids, we made use of a simple two-level plus continuum model atom to keep the computational complexity down. We built our simplified model atom based on the hydrogen atom, but with some important differences. First, we used the first two levels of hydrogen, so there is only one spectral line (Lyman-α-like), and two bound-free transitions. Second, we artificially increased the collisional excitation and ionisation rates, originally derived from Johnson (1972), by a factor of 10⁹, to decrease the amount of scattering so that we were able to reduce the number of Λ iterations and use the simpler Λ-iteration scheme. In our model, the lowest value of the destruction probability was ε_λ ≈ 0.09.

We adopted the level energies and statistical weights shown in Table 1, along with an ƒ-value of 0.41641 and a spontaneous deexcitation rate A_ul of 4.7 × 10⁸ s⁻¹, both taken from the NIST database (Kramida et al. 2021). The line and continuum extinction were computed using the Transparency.jl package (Pereira 2021). We calculated the intensity for a total of 90 wavelength points, 50 for the bound-bound transition, and 20 for each of the bound-free transitions. The line profile was computed assuming complete redistribution (CRD) and using standard recipes for van der Waals broadening (Unsöld 1955), linear Stark broadening (Sutton 1978), and quadratic Stark broadening (Gray 2005). For the continuum extinction, we included the following processes: H⁻ free-free based on Stilley & Callaway (1970), H⁻ bound-free following Geltman (1962), H free-free as per Mihalas (1978), ${\text{H}}_2^{+}$ free-free and bound-free according to Bates (1952), Rayleigh scattering, and Thomson scattering.

Fig. 1 Surface temperature and vertical magnetic field from the Bifrost simulation used. The horizontal box size is 6 × 6 Mm2.

2.3 Model atmosphere

As our 3D model atmosphere, we used a single snapshot of a simulation ran with the Bifrost code (Gudiksen et al. 2011). Bifrost solves the magnetohydrodynamic equations on a Cartesian grid, and the simulation employed here included thermal conduction, radiation using a multi-group opacity scheme (Nordlund 1982; Skartlien 2000) including scattering (Hayek et al. 2010), and a recipe for NLTE radiative losses in the higher layers (Carlsson & Leenaarts 2012). This simulation did not include the effects of hydrogen non-equilibrium ionisation.

The simulation had 256 × 256 × 430 grid points, or about 28 million cells. These translate into a physical size of 6 × 6 Mm² in the horizontal direction and 8.7 Mm in the vertical direction. This means that the horizontal spatial pixel size is 23.4 × 23.4 km². The vertical grid spacing varied with height, and its average was 20.3 km.

The magnetic field configuration of the simulation was chosen to mimic a quiet coronal hole, with a mean unsigned magnetic field at the surface of 0.5 mT and a mostly unipolar configuration. We display the surface temperature and vertical magnetic field in Fig. 1.

2.4 Constructing the irregular grid

Our main goal was to investigate whether irregular grids can be used to speed up detailed 3D NLTE radiative transfer calculations. Therefore, a critical step in the process was to construct an irregular grid that is better suited for radiation. There are different ways to build irregular grids, but we employed one of the most general cases: a 3D Voronoi diagram (Voronoi 1908).

A Voronoi diagram is calculated from a set of control points called sites. A Voronoi diagram is the set of all cells associated with each site, where a cell is defined as the region that is closer to its control point than to any other control point. The Voronoi cells give us two useful properties for radiative transfer: cell neighbours, and Delaunay lines. Neighbours are sites whose Voronoi cells have a common border, and Delaunay lines are lines connecting neighbouring sites (Delaunay 1934). Our ray tracing algorithm traces rays along Delaunay lines.

The first step in building a Voronoi diagram from a 3D model atmosphere is to decide where to place the sites, or cell centres, and how many sites to use. As we discuss below, this will have a critical effect on the resulting radiation and speed of convergence. We wish to place a higher density of points in regions that are more important for radiative transfer and therefore follow the distribution of a given quantity. For a given number of sites, we built the Voronoi diagram using the Monte Carlo technique of rejection sampling. In brief, this involves associating a probability distribution with a reference distribution of the quantity we wish to sample, determining how likely it is to find a control point in a given area. For a given target distribution of sites ƒ(r), the process is as follows:

1. We randomly generated a proposed position r_prop and obtain the proposal density p_prop = ƒ(r_prop).

2. A random rejection density p_rej was drawn from a uniform distribution.

3. The proposal density was compared with the rejection density. If p_prop ≥ p_rej, r_prop was accepted as a control point (site). Otherwise, r_prop was rejected.

4. We repeated the steps above until the target number of sites was reached.

In our work, the target function ƒ is a quantity that aims to optimise the computational grid. It should have a high value in areas that are important for formation, and lower values in less important areas. We experimented with different quantities from which to sample the location of the sites, but we focus here only on the three most promising ones: the continuum extinction ${\alpha }_{\lambda }^c$ at a given wavelength (typically the line core), the log number density of protons under LTE log($N_{\text{H\hspace{0.17em}II}}^{\text{LTE}}$) and the quantity log(N_H)⁻²T^−2/5, a combination of the total hydrogen number density N_H and temperature.

When the sites of the irregular grid were chosen, we interpolated the atmospheric quantities from the regular grid of the simulation into the new sites using trilinear interpolation. The quantities we interpolated were the temperature, the total hydrogen number density, the electron density, and the velocity vector.

To construct a Voronoi diagram of the Bifrost simulation, we used the open-source library voro++ (Rycroft 2009). This provided the location of the sites and a list of neighbours for all cells, from which we calculated the Delaunay lines that we used for the ray tracing. In Fig. 2, we show an example Voronoi grid constructed from our Bifrost simulation, using log(N_H) to determine the site density. For clarity, this example contains far fewer cells than the original simulation. Still, the site density clearly drops in the higher layers that map the solar corona because its density is much lower.

Fig. 2 Illustration of a Voronoi grid sampled from the Bifrost simulation. To clarify some features of the grid, the number of points is significantly lower than the cells in the original grid. This example grid was constructed to follow the log hydrogen density, log(NH). The white spheres represent the sites, and red lines indicate the edges between cell borders.

2.5 Rays on irregular grids

The standard recipes for long characteristics (Jones & Skumanich 1973; Jones 1973) or short characteristics (Kunasz & Auer 1988) methods are not directly applicable on irregular grids. To calculate radiation through an irregular grid, we need a method that does not require even partitioning between grid points. Different approaches can be used. Camps et al. (2013) developed a method for tracing rays through straight paths in irregular grids. This has the advantage of allowing us to trace straight rays through the whole domain for any direction, but requires additional computations to find the walls of each cell and entry and exit points of the ray. Ritzerveld & Icke (2006) used a different approach, in which rays were only allowed to travel along Delaunay lines. This approach is also followed in other codes such as LIME (Brinch & Hogerheijde 2010) and SIMPLEX2 (Paardekooper et al. 2010). This approach is simpler and requires fewer interpolations, but a downside is that rays can only travel along the predetermined directions that connect cell centres. However, as Paardekooper et al. (2010) demonstrated in their Fig. 5, the direction of each segment can be chosen, such as the overall direction (after the ray has travelled through several cells), which remains very close to a target direction. In this work, we followed the approach in which the rays travel only through Delaunay lines because it is computationally simpler. Article number, page 4 of 10

Our approach to ray tracing was inspired by the short-characteristics (SC) method of Kunasz & Auer (1988). We wished to calculate the intensity at every site, both for upward-travelling rays starting at the bottom boundary and for downward-travelling rays starting at the top boundary. While LIME and SIMPLEX2 use a Monte Carlo approach to sample different directions, we instead adopted a fixed-angle quadrature, which we used to integrate the intensity in different directions and obtain J_λ. Štěpán et al. (2020) and Jaume Bestard et al. (2021) have developed optimised angle quadratures for 3D radiative transfer in stellar atmospheres, and here we adopted the quadrature optimised for unpolarised radiative transfer with n_rays = 12 and L = 7 from Jaume Bestard et al. (2021).

For each ray in the quadrature, we computed the intensity by tracing rays along the Delaunay lines that lie closest to the characteristic angle. A simple sketch of the ray tracing method is shown in Fig. 3.

In some parts of our atmosphere, the irregular grid can have sparse sampling, with large cells. In these cases, limiting the transport to a single Delaunay line can lead to errors in the formal solution. To better sample the grid, we implemented a method for approximating the intensity by using not only a single Delaunay line, but the two Delaunay lines that have the smallest angle from the characteristic ray, $$I_i=w_1{I}_{i,1}+w_2{I}_{i,2},$$(5)

where we designed the weights w_i to favour Delaunay lines that have a smaller angle to the target quadrature angle, and normalised them, $$w_i=\frac{{\upsilon }_i^r}{{\displaystyle {\sum }_{j=k_1,k_2}{\upsilon }_j^r}},$$(6)

where v_i is the dot product between the direction of the ray and the direction of the Delaunay line, and the exponent r is treated as a free parameter. After some experimentation, we find that r = 7 gives reasonable results, but overall, the method is not very sensitive to the particular choice of r.

The transport of radiation can be designed in different ways (e.g. combining the intensity from additional neighbours or adopting different weights), but we find that our algorithm works well to correct the discrepancy between a quadrature angle and the angles of the Delaunay lines. In an irregular grid with sites drawn from a 3D Poisson distribution, a site has 15.54 neighbours on average (van de Weygaert 1994). With our 3D atmosphere with ≈300000 sites drawn according to the hydrogen density, the average number of neighbours was 15.44. Consequently, every site should have some Delaunay lines whose direction differs from a quadrature ray by only a small amount.

To compute I_λ along each Delaunay line, we solved the integral in the formal solution (2) by using linear interpolation of S_λ. This accuracy can be improved by using higher-order schemes (e.g. de la Cruz Rodriguez & Piskunov 2013; Janett et al. 2019). However, given the exploratory nature of this work and several of the approximations already used, the additional accuracy of higher-order methods appeared to be not worth the additional computational complexity.

Our algorithm for computing I_λ along different rays was used to obtain J_λ, which we then used in the Λ-iteration to obtain the source function and the level populations. After the iteration converged, the end product should be the emergent I_λ at a given direction. To obtain I_λ, we could have used the same ray tracing algorithm with the converged level populations. However, this would create undesirable effects in this type of simulation. A side-effect of placing the number of sites in the locations that are most important for the transport of (optically thick) radiation is that the upper layers of the simulation box, covering the corona, have far fewer cells. This can be seen, to some extent, in Fig. 2. This means that in the irregular grid, I_λ at the top of the simulation box, comes from a small number of sites. This results in poor spatial resolution. Therefore, instead of using the irregular grid to compute the final I_λ, we instead interpolated the converged level populations to the original grid of the simulation and then calculated a formal solution using long characteristics. The choice of interpolation mapping the irregular grid back to the original mesh is important for the intensity calculations. We used inverse distance interpolation with the two nearest neighbours because it reduces noise from the random sampling of the grid and preserves sharp features in intensity. The interpolation allowed us to obtain I_λ at the original grid spacing of the simulation. For the tests performed here, we compared only vertical rays (µ = 1), but any viewing angle can be used with this method.

Fig. 3 Two-dimensional representation of ray tracing in an irregular grid. We wish to calculate the intensity at the red point, and the thick black arrow is the direction of the ray (characteristic). In our algorithm, the incoming intensity is a combination of I0,1 and I0,2, the intensity coming from the two sites with the smallest angles between the Delaunay lines and the direction of the ray, θ1 and θ2.

Fig. 4 Emerging intensities from searchlight beam tests. Left: regular grid with short characteristics and linear interpolation. Right: irregular grid. The box we used had a size of 1 × 1 × 1 m3 and 5l3 grid points. The beam was circular, with a radius of 0.2 m and inclination of θ = 20° (µ ≈ 0.94).

2.6 Implementation

Our implementation of the methods described above is open source and publicly available (Udnæs 2023)¹ Our code is written in the programming language Julia (Bezanson et al. 2017). Its relatively simple architecture makes use of Julia’s thread parallelism, which can run across different cores in a shared memory architecture. The code is not domain decomposed and does not make use of more advanced parallel architectures.

For a fair comparison with existing radiative transfer methods, our code also included an implementation of the short-characteristics scheme of Kunasz & Auer (1988) that we ran on a regular grid. This allowed us to run exactly the same functions to calculate extinction, collisional and radiative rates, Λ-iteration, statistical equilibrium, and the final formal solution. The only difference is the method for computing J_λ, which can be selected as regular grid or irregular grid. Hereafter, when we refer to calculations on the “regular grid”, we mean calculations performed using short characteristics in a regular grid, using linear interpolation to obtain the simulation quantities at the necessary locations.

3 Results

3.1 Searchlight beam test

Our first test of our ray tracing in irregular grids was a searchlight beam test. In this test, we created a box in which the extinction α_λ and source function S_λ were zero everywhere and propagated a radiation beam from the bottom to the top of the box at a given angle. We chose a box size of 1 × 1 × 1 m³ and a total of 5l³ sites. From the centre of the bottom boundary, we propagated a circular beam with I_beam = 1 kWm⁻² nm⁻¹ sr⁻¹, and a radius of 0.2 m.

For this test only, the Voronoi cells were randomly sampled, and not sampled from a simulation quantity. If we had sampled the sites from a quantity such as mass density, the density of sites at the top of the box (where we assess the results) would be too low, and it would be a poor test of the diffusion introduced by the algorithm. Therefore, the 5l³ sites were sampled from a uniform random distribution to ensure a similar density throughout the box.

In Fig. 4, we show the results of a beam with an inclination of 20° from the vertical (µ ≈ 0.94). We compare them with the results from the regular grid. Ray tracing on the irregular grid compares very favourably with results from the regular grid. The amount of diffusion introduced by the irregular grid solver is comparable with that of short characteristics with linear interpolation. For this setup, the peak intensity of the beam is reduced to about 70% of the original value in both cases. For shallower rays, the diffusion in regular grids increases (de Vicente et al. 2021), and again we find comparable results for the irregular grid.

3.2 Radiative transfer on irregular grids in LTE

The next step was a calculation of intensity in LTE. Strictly speaking, this was no test of the irregular grid solver because no radiation was calculated with the irregular solver. Instead, it was a test of the procedure to construct an irregular grid and to interpolate quantities from the regular to irregular grid and vice versa.

A NLTE problem revolves around iterating the source function until the radiation field is consistent with the atomic level populations. The populations from a converged Λ-iteration are the end product that is then used to carry out spectral synthesis. In this LTE test, our procedure was to start by computing the LTE populations for our hydrogen-like atom in LTE in the original grid of the Bifrost simulation. Next, we built an irregular grid and interpolated the populations into to the irregular grid. As no iteration took place, the populations were then interpolated back into the original regular grid. Finally, we computed the emerging continuum intensity at 500 nm for μ = 1 from the regular grid using the interpolated populations. The source function was computed in LTE using only quantities in the original grid, but the extinction and optical depth make use of the interpolated populations because we used them as hydrogen populations, and therefore, they were key in setting the H⁻ populations, whose bound-free and free-free contributions dominate the continuum extinction.

For this test, we built the irregular grids using two methods: by sampling α₅₀₀, the continuum extinction at 500 nm, and by sampling the vertical temperature gradient. We built grids with a different number of sites, from 5 × 10⁵ to 3 × 10⁶ sites. For comparison, the original grid has 2.8 × 10⁷ grid points.

We show the results of the continuum intensity in LTE in Fig. 5. For this case, it is obvious that sampling the irregular grid from α₅₀₀ produces much better results than sampling from the vertical gradient of temperature because the noise from the interpolation is much less noticeable. With the grid sampled by α₅₀₀, the granulation pattern, overall contrast, and absolute intensities become close to the result of the regular grid already with only 5 × 10⁵ sites (almost two orders of magnitude fewer grid points than the original grid). With 3 × 10⁶ sites, the results from the irregular grid are virtually indistinguishable from the regular grid. This is encouraging as it shows that the irregular grid can sample the continuum-forming region very well, allowing calculations with nearly ten times fewer grid points.

3.3 Radiative transieron irregular grids in NLTE

Our final and most extensive step was to do a full NLTE calculation using our simplified model atom with one spectral line and two bound-free edges. With the collisional rates artificially scaled up, we reduced the amount of scattering and ensured that the problem was solvable by using the simple Λ-iteration. Since our code for these experiments has only simple parallelism, running the full 3D NLTE problem is still computationally expensive.

To avoid lengthy running times, we did not run the experiment at the native resolution of the Bifrost simulation (or equivalent number of sites in the irregular case) because this greatly increased not only the time per iteration, but also the total number of iterations necessary to achieve convergence. Instead, we ran the regular grid at one-half resolution, which we defined as taking every second point from the original grid in all three directions. This means that the half-resolution grid has eight times fewer grid points, or about 3.5 million points. For the corresponding irregular grid test, we used three million sites. In addition, we also ran tests using one-third and one-fourth resolution, and the same number of sites in the irregular cases, but their main purpose was to time the differences.

We experimented with different irregular grids that sampled the distribution of three quantities: log₁₀(α_cont), the logarithm of the continuum extinction at the line core; log($N_{\text{H\hspace{0.17em}II}}^{\text{LTE}}$), the logarithm of the proton density under LTE; and log(N_H)⁻²T^−2/5, where N_H is the total hydrogen density, and T is the temperature. The rationale for these three choices was as follows. The continuum extinction at the line core is derived mostly from quantities computed in LTE such as the H⁻ populations and is proportional to the mass density. Following the results from the LTE radiative transfer tests, this results in a grid that is optimised for the line continuum wavelengths (but not necessarily for the line profile). The proton density under LTE follows the mass density, but also the electron density and temperature via the Saha-Boltzmann distribution, and it is more sensitive to changes in the upper layers than α_cont· Finally, log(N_H)⁻²T^−2/5 was an ad hoc formulation found by a little experimentation, with the goal of placing a denser grid in the regions most critical to the line formation. Ideally, an irregular grid would be built by sampling the quantities that are most relevant for radiation at all the wavelengths considered, for example the line extinction or level populations. However, this approach is moot because the values of these quantities are not known before the NLTE calculations, so that we are restricted to the quantities that are available at the start of the calculations (or that are derived from these quantities).

In Fig. 6, we compare the average density of sites or grid points per height as a function of height in the simulation. We also plot histograms of the height of the optical depth unity for the line core and line wing, which is a proxy for the formation height at these wavelengths. The regular grid has a variable height resolution and therefore has a higher resolution in the photosphere and chromosphere than in the corona. The irregular grid sampling log₁₀(α_cont) has a higher density of sites in the photosphere, while log($N_{\text{H\hspace{0.17em}II}}^{\text{LTE}}$) follows a similar curve, but has a lower resolution in the deeper layers and a better resolution from about 2 Mm to the top of the box. The grid following log(N_H)⁻²T^−2/5 has the lowest resolution of all in the photosphere, but a larger number of sites between 1 and 3 Mm, covering the region of formation of our modelled spectral line.

We find that the irregular grid sampled by log(N_H)⁻T^−2/5 gives the most accurate results when compared with the regular grid. In Fig. 7, we show emergent intensity images at μ = 1 for both the regular and irregular grid at three wavelengths: the line wing, the line core, and the 500 nm continuum, and in Fig. 8 we show the μ = 1 spectral profiles. For the sake of brevity, we omitted the results from irregular grids sampled by log₁₀(α_cont) and log($N_{\text{H\hspace{0.17em}II}}^{\text{LTE}}$).

Figure 7 shows that the results are very similar overall in shape and intensity level, but that the irregular grid shows some spurious results of points with very high intensities in the line wing and core (black dots). These points are more visible in the line core images, but are still a small subset of the total number of points. They are also visible in Fig. 8 as the very broad and strong profiles. Nevertheless, most profiles agree very well between the irregular and regular grid calculations. The occurrence of the spurious strong profiles seems to be related to how well the irregular grid covers the line-forming region. They are much more numerous when we use fewer sites (e.g. 10⁵ instead of 3 × 10⁶), and they are also more numerous when we use irregular grids sampled from log₁₀(α_cont) or log($N_{\text{H\hspace{0.17em}II}}^{\text{LTE}}$) (not shown).

The 500 nm continuum images of Fig. 7 also show that the irregular grid produces a noisy map. This is because the irregular grid sampling log(N_H)⁻²T^−2/5 has far fewer sites in the photosphere, and the noise recedes with the grids built by sampling log₁₀(α_cont) or log($N_{\text{H\hspace{0.17em}II}}^{\text{LTE}}$).

In addition to the accuracy of the results, an important aspect of this analysis is to determine whether irregular grids can save computer time when compared to calculations in regular grids. For the computational costs of running a single iteration, we list in Table 2 the time it took to run a single iteration on the regular and irregular grids for different numbers of grid points. The regular grid experiments were performed at one-half, one-third, and one-quarter resolution, and for each case the corresponding irregular grids had exactly the same number of grid points (sites). We performed the experiment using five cores of an AMD EPYC 7763 2.45 GHz CPU and averaged the times over three iterations. There was no difference in time per iteration from the different irregular grids.

We plot in Fig. 9 the convergence rate versus the number of iterations for the three irregular grids and the regular grid to show how many iterations were needed to achieve convergence. These 3D NLTE problems are usually assumed to converge at a level around 10⁻³. The grid sampling log(N_H)⁻²T^−2/5 converges in about the same number of iterations as the regular grid, while the other two converge with about one-third fewer iterations.

Fig. 5 Disk-centre continuum intensity at 500 nm from the Bifrost simulation. Each panel was computed using a different grid. The bottom right panel shows the regular grid, and the other three panels show irregular grids with a different number of sites and constructed sampling distributions of dT/dz and α500, the continuum extinction at 500 nm. All calculations assumed LTE.

Fig. 6 Distributions of the vertical density of the grid points for different grids. The regular grid (black) was taken at half resolution, meaning every second grid point. The different irregular grids (blue, cyan, and yellow) were sampled by the quantity indicated in the legend and had three million sites. The shaded histograms (red) show the τ = 1 heights at the line core and the wing at −27 km s−1 from the core.

4 Discussion

Our main aim with this work was to determine whether it is feasible to use irregular grids to speed up 3D NLTE radiative transfer, just like optimised height grids can be used to speed up ID NLTE calculations. We employed one of the most general irregular grids, based on a 3D Voronoi diagram. In this first exploration, we employed a simplified approach with an algorithm that is not massively parallel and runs only on shared-memory systems. To further reduce the complexity of the code, we also used the simpler Λ-iteration scheme, which breaks down in the strong scattering regimes. To reduce computation times, we used a simplified hydrogen-like model with two bound states and one continuum and assumed complete redistribution. To make the problem tractable given these simplifications, we artificially increased the collisional rates.

We focus on the spectral profiles of the line from our model atom, a modified Lyman-α line. Because we reduced the amount of scattering by increasing the collisional rates, the source function is closer to LTE. Compared to Lyman-α, our line is formed in hotter regions and lacks the central reversal. This explains why in both absolute intensity units and in brightness temperature it is consistently brighter than Lyman-α from observations or simulations (e.g. Schmit et al. 2017).

Given the assumptions above, our experiments still give us valuable insight. Our searchlight beam test showed that our ray tracing algorithm in irregular grids produces good results, showing a diffusion that is comparable to short characteristics with linear interpolation. The test of synthesising continuum intensities in LTE probed the effects of interpolating to a Voronoi grid and back to the original grid, and it gave two key results. First, it demonstrated that it is crucial to build the irregular grid on a quantity relevant for the transport of radiation, with α₅₀₀ being the most relevant quantity here. Second, it showed that a properly optimised irregular grid can facilitate the calculation of continuum radiation with about an order of magnitude fewer grid points than the original grid. However, we recall that continuum radiation comes from a relatively shallow region, so that similar gains can probably be achieved with regular grids by limiting the z range. A much more stringent test is the NLTE formation of our spectral line, where our results are more mixed.

When we solved the 3D NLTE problem for our model atom, we had to compute radiation for 90 wavelengths, covering both the spectral line and bound-free edges, and solve the statistical equilibrium equations. This test is more stringent on the irregular grid because the range of formation heights is much broader, so that an ideal irregular grid should provide good resolution for wavelengths from the line core to the continuum, which span a formation region of more than 4 Mm. Finding a relevant quantity from which to sample the irregular grid in this case is a much more challenging task. Especially because many of the relevant quantities (e.g. line extinction) are not known a priori, before we solved the NLTE problem. The design of the irregular grid has dramatic effects on the resulting populations and line profiles in this case as well.

Irregular grids that sample the distribution of continuum extinction or LTE proton density lead to unsatisfactory results: the emergent line profiles show many instances with excessive intensities that show up noticeably in intensity maps. This is perhaps to be expected from Fig. 6: the grid resolution in these cases is much lower than the regular grid in the regions in which the line is formed. It appears that the spurious line profiles arise from large jumps in the populations when interpolating from coarser grids. To address this issue, we built an irregular grid sampling the seemingly arbitrary quantity log(N_H)⁻²T^−2/5, which led to a much higher density of sites in the line-forming region. With this irregular grid, we still obtained a few spurious profiles (as seen in Figs. 7–8), but they are far less prevalent, and the profiles and intensity level agree well with the results from the regular grid overall. A side-effect from this particular grid is that the photospheric radiation is adversely affected by noise from the interpolation because this grid has a lower resolution in the photosphere. Despite these issues, our results show that it is possible to use irregular grids for 3D NLTE radiative transfer. The next important question is therefore whether this leads to faster run times than regular grids.

We find that one iteration on a irregular grid is about four times slower than one iteration in the regular grid. This is comparable with the results of Camps et al. (2013), who reported that radiative transfer in Voronoi grids was about three times slower than on octree grids (and radiation transport across octree grids is already more complex than our regular grid). However, we find that our most reliable irregular grid reached convergence in about the same number of iterations as the regular grid (≈90 iterations to achieve a 10⁻³ level). This means that in this case, irregular grids are much slower to solve the 3D NLTE problem. Some irregular grids can converge faster, in about two-thirds of the regular grid iterations, but the lower number of iterations is not enough to balance the price penalty of iterations that are four times slower. Tellingly, the irregular grids that converge faster have a lower density of grid points in the line-forming region, which leads us to conclude that they converge faster precisely because there are fewer grid points through which the changes to the source function need to be propagated. These faster-converging grids also have the more serious problem of leading to additional spurious line profiles.

Fig. 7 Disk-centre intensity for three different wavelengths calculated with the half-resolution regular grid (top panels) and with an irregular grid with three million sites sampled from log(NH)−2T−2/5 (bottom panels). The blue wing position is at −27kms−1 from the line centre.

Fig. 8 Individual line profiles for the disk-centre intensity calculated from the half-resolution regular grid and an irregular grid with three million sites sampled from log(NH)−2T−2/5. Every black line corresponds to the brightness temperature from a column in the atmosphere. The spatial average over all columns is represented by the red line. The dashed blue lines indicate the wavelength positions of the blue wing at 121.558 nm (or −27 kms−1) and line centre λ0 at 121.568 nm.

5 Conclusions

We conclude that irregular grids alone do not provide the optimal solution for efficient 3D NLTE calculations. With our setup, they are about four times slower and lead to less accurate results. The design of an irregular grid that lends itself to radiative transfer across a large height range is not trivial and affects both the reliability of the results and the speed of convergence. It may be possible to build a better-performing irregular grid sampling other quantities, but this work goes beyond our exploratory analysis. NLTE radiative transfer in irregular grids is still affected by the curse of dimensionality, just like regular grids, where going to 3D and increasing the spatial resolution leads to slower convergence.

Irregular grids might still be attractive through their potential to achieve results similar to a regular grid with fewer grid points. Based on our results for the LTE continuum intensity, we speculate that that the NLTE results from our irregular grid with 3 × 10⁶ sites, when interpolated back to the original simulation grid, could give comparable results to a calculation on the full resolution regular grid with ≈3 × 10⁷ grid points. We did not test this due to the long running times of the full resolution in a regular grid. These regular grid full-resolution runs should be about eight times slower than our half-resolution runs, and in these cases, it is possible that the irregular grid running times would outweigh the four times slower iterations.

Fig. 9 Λ-iteration convergence for irregular grids with three million sites constructed from the labelled quantities, and the half-resolution regular grid of the Bifrost atmospheric model.

Acknowledgements

We would like to thank Mats Carlsson for helpful suggestions, fruitful discussions, and for providing the Bifrost atmosphere used in this paper. We would also like to thank Jaime de la Cruz Rodriguez for helpful discussions and suggestions for improvements. This work has been supported by the Research Council of Norway through its Centers of Excellence scheme, project number 262622. We kindly acknowledge the computational resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway

References

All Tables

All Figures

	Fig. 1 Surface temperature and vertical magnetic field from the Bifrost simulation used. The horizontal box size is 6 × 6 Mm2.
In the text

	Fig. 2 Illustration of a Voronoi grid sampled from the Bifrost simulation. To clarify some features of the grid, the number of points is significantly lower than the cells in the original grid. This example grid was constructed to follow the log hydrogen density, log(NH). The white spheres represent the sites, and red lines indicate the edges between cell borders.
In the text

Fig. 3 Two-dimensional representation of ray tracing in an irregular grid. We wish to calculate the intensity at the red point, and the thick black arrow is the direction of the ray (characteristic). In our algorithm, the incoming intensity is a combination of I0,1 and I0,2, the intensity coming from the two sites with the smallest angles between the Delaunay lines and the direction of the ray, θ1 and θ2.

In the text

	Fig. 4 Emerging intensities from searchlight beam tests. Left: regular grid with short characteristics and linear interpolation. Right: irregular grid. The box we used had a size of 1 × 1 × 1 m3 and 5l3 grid points. The beam was circular, with a radius of 0.2 m and inclination of θ = 20° (µ ≈ 0.94).
In the text

	Fig. 5 Disk-centre continuum intensity at 500 nm from the Bifrost simulation. Each panel was computed using a different grid. The bottom right panel shows the regular grid, and the other three panels show irregular grids with a different number of sites and constructed sampling distributions of dT/dz and α500, the continuum extinction at 500 nm. All calculations assumed LTE.
In the text

Fig. 6 Distributions of the vertical density of the grid points for different grids. The regular grid (black) was taken at half resolution, meaning every second grid point. The different irregular grids (blue, cyan, and yellow) were sampled by the quantity indicated in the legend and had three million sites. The shaded histograms (red) show the τ = 1 heights at the line core and the wing at −27 km s−1 from the core.

In the text

	Fig. 7 Disk-centre intensity for three different wavelengths calculated with the half-resolution regular grid (top panels) and with an irregular grid with three million sites sampled from log(NH)−2T−2/5 (bottom panels). The blue wing position is at −27kms−1 from the line centre.
In the text

Fig. 8 Individual line profiles for the disk-centre intensity calculated from the half-resolution regular grid and an irregular grid with three million sites sampled from log(NH)−2T−2/5. Every black line corresponds to the brightness temperature from a column in the atmosphere. The spatial average over all columns is represented by the red line. The dashed blue lines indicate the wavelength positions of the blue wing at 121.558 nm (or −27 kms−1) and line centre λ0 at 121.568 nm.

In the text

	Fig. 9 Λ-iteration convergence for irregular grids with three million sites constructed from the labelled quantities, and the half-resolution regular grid of the Bifrost atmospheric model.
In the text