© ESO 2019

1. Introduction

Understanding of hot, massive stars is a basic prerequisite for interpreting fundamental properties of our Universe. Already during their lifetimes, such stars influence the evolution of galaxies via feedback of ionizing radiation and radiatively driven winds, and they enrich the interstellar medium (ISM) with metals. But also, their deaths are of large impact. Depending on initial mass and mass loss, OB stars either explode as supernovae or end up as “heavy” stellar-mass black holes (Heger et al. 2003). Obviously, supernova explosions enrich the ISM with metals even further. Additionally, the associated shock fronts possibly trigger star formation, resulting in a new generation of (massive) stars.

With the advent of gravitational wave (GW) observations (e.g. the black hole merger GW150914 observed at the advanced Laser Interferometric Gravitational-Wave Observatory (aLIGO), Abbott et al. 2016), the formation of heavy stellar-mass black holes becomes of key interest. Since massive stars are frequently found to be members of multiple star systems (see, e.g. Mason et al. 2009; Sana et al. 2013), they might explain the occurrence of GW events in the correct mass range. At least the formation of heavy black holes from single-star evolution, however, requires comparatively moderate mass loss rates (e.g. in low metalicity environments, or mass-loss quenching by magnetic fields, Petit et al. 2017; Keszthelyi et al. 2017).

Current knowledge of OB stars is inferred from quantitative spectroscopy, that is, by comparing observed spectra with synthetic ones, the latter being obtained from numerically modelling their stellar atmospheres (photosphere + wind). State of the art atmospheric modelling is performed by assuming spherical symmetry (e.g. CMFGEN: Hillier & Miller 1998; PHOENIX: Hauschildt 1992; POWR: Gräfener et al. 2002; WM-BASIC: Pauldrach et al. 2001; FASTWIND: Puls et al. 2005; Rivero González et al. 2012).

Several effects may alter the geometry of hot star atmospheres, affecting both the photospheric and wind lines. For instance, binary interaction (see, e.g. Vanbeveren 1991; de Mink et al. 2013 for a discussion about evolutionary aspects, and Prilutskii & Usov 1976; Cherepashchuk 1976; Stevens et al. 1992 for the effects of colliding winds) influences the line formation in such objects. Another example is the phenomenon of magnetic winds (e.g. ud-Doula & Owocki 2002; ud-Doula et al. 2008; Petit et al. 2013), for which the resonance line formation has been extensively discussed in Marcolino et al. (2013), Hennicker et al. (2018), and David-Uraz et al. (2019), for instance. Additionally, (non-radial) pulsations may impact the geometry of hot star atmospheres.

In this paper, we focus on (rapidly) rotating stars. In addition to polar velocity components that affect the radiative line-driving, such stars and their winds are influenced by centrifugal forces and gravity darkening. For a detailed discussion, we refer to Sect. 5.

To analyse these effects, and to distinguish between different theories (resulting in, e.g. prolate vs. oblate wind structures), multi-D radiative transfer codes that include a treatment of (arbitrary) velocity fields are required. In hot star winds, a key challenge is the implementation of scattering processes. Such problems are most easily handled by an accelerated Λ-iteration scheme, that is, by calculating the radiation field for known sources and sink terms (the so-called formal solution), and iterating the updated sources and sinks until convergence (Cannon 1973).

In order to obtain the formal solution in multi-D, the following two major methods exist, each having specific advantages and disadvantages. Firstly, within the finite-volume method (FVM, see Adam 1990 for the first application in the context of radiative transfer problems), the calculation volume is discretized into finite sub-volumes. All required physical values at the cell centres are then considered as suitable averages within each cell. This way, a relatively simple solution scheme can be derived, which, however, is only of low order. Thus, the FVM breaks down for large optical depths (Hennicker et al. 2018, hereafter Paper I). Nevertheless, it is still useful for qualitative interpretations (see, e.g. Lobel & Blomme 2008 for an application to corotating interaction regions in the winds of rotating stars). Secondly, within the characteristics methods, the formal solution is obtained by exact integration of the equation of radiative transfer along a ray. One distinguishes between the “long-characteristics” method (LC, Jones & Skumanich 1973; Jones 1973), and the “short-characteristics” method (SC, Kunasz & Auer 1988). The LC method follows a particular ray from the boundary to the considered grid point, whereas the SC method considers rays only from cell to cell. All required quantities at each upwind point are obtained from interpolation. While the LC method is computationally more expensive than the SC method, the latter introduces numerical errors from the interpolation scheme. To date, few 3D codes that apply one or the other technique already exist, such as PHOENIX/3D (Hauschildt & Baron 2006 and other publications in this series) and IRIS (Ibgui et al. 2013), which, however, are either proprietary, or do not include a suitable Λ-iteration scheme. For additional information on other codes, we refer to Sect. 1 of Paper I.

In this paper, we present a newly developed 3D SC code, which is capable of treating arbitrary velocity fields, and includes an accelerated Λ-iteration scheme to handle continuum and line scattering problems. In Sect. 2, we discuss the basic assumptions of the code. The numerical methods are described in Sect. 3. In Sect. 4, we perform a detailed error analysis by comparing our 3D results for spherically symmetric atmospheres with corresponding 1D solutions. Additionally, we present comparisons with the 3D finite-volume method from Paper I, when appropriate. As a first application to non-spherical winds, we calculated UV resonance line profiles for different models of rotating stars in Sect. 5, and discuss the implications. Our conclusions are summarized in Sect. 6.

2. Basic assumptions

In this paper, we use the same basic assumptions as in Paper I, that is we aim to solve the time-independent equation of radiative transfer (e.g. Mihalas 1978), formulated in the observer’s frame:

$$\begin{array}{c}\hfill \mathit{n}\mathbf{\nabla }{I}_{\nu }(\mathit{r},\mathit{n})={\chi }_{\nu }(\mathit{r},\mathit{n})(S_{\nu }(\mathit{r},\mathit{n})-I_{\nu }(\mathit{r},\mathit{n})).\end{array}$$(1)

I_ν describes the specific intensity for a given frequency ν and direction n, and χ_ν and S_ν are the opacities (in units of cm⁻¹) and source functions for continuum and line processes. For our following tests, we either consider a pure continuum case with source function

$$\begin{array}{c}\hfill S_{\mathrm{C}}=(1-{ϵ}_{\mathrm{C}})J_{\nu }+{ϵ}_{\mathrm{C}}{B}_{\nu },\end{array}$$(2)

and thermalization parameter ϵ_C, or a line (within an optically thin continuum) approximated by a two-level-atom (TLA). A generalization to a multitude of lines, coupled via corresponding rate equations, is one of our next objectives. Within the TLA approach, the line source function reads:

$$\begin{array}{cc}\hfill S_{\mathrm{L}}& =(1-{ϵ}_{\mathrm{L}})\overline{J}+{ϵ}_{\mathrm{L}}B\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill {ϵ}_{\mathrm{L}}& ={\displaystyle \frac{{ϵ}^{\prime }}{1+{ϵ}^{\prime }}},{ϵ}^{\prime }={\displaystyle \frac{C_{\mathrm{ul}}}{A_{\mathrm{ul}}}}[1-exp(-\frac{h\nu }{k_{\mathrm{B}}T})]·\hfill \end{array}$$(4)

$\overline{J}$ is the “scattering integral” (see Eq. (25)), and C_ul and A_ul describe the collisional rate and Einstein coefficient for spontaneous emission, respectively. In this paper, ϵ_L is considered as input parameter. The profile function Φ_x is approximated by a Doppler profile:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_x={\displaystyle \frac{1}{\sqrt{\pi }\delta }}exp[-({\displaystyle \frac{x_{\mathrm{cmf}}}{\delta }}{)}^2]={\displaystyle \frac{1}{\sqrt{\pi }\delta }}exp[-{\left({\displaystyle \frac{x_{\mathrm{obs}}-\mathit{n}·\mathit{V}}{\delta }}\right)}^2],\end{array}$$(5)

where x_cmf and x_obs denote the frequency shift from line centre in the comoving and observer’s frame, in units of a fiducial Doppler width, $\mathrm{\Delta }{\nu }_{\text{D}}^{\ast }={\nu }_{\text{0}}{v}_{\text{th}}^{\ast }/c$. V is the local velocity vector in units of $v_{\text{th}}^{\ast }$, and

$$\begin{array}{c}\hfill \delta ={\displaystyle \frac{1}{v_{\mathrm{th}}^{\ast }}}\sqrt{{\displaystyle \frac{2k_{\mathrm{B}}T}{m_{\mathrm{A}}}}+v_{\mathrm{micro}}^2}\end{array}$$(6)

is the ratio between the local thermal velocity (accounting for micro-turbulent velocities) and the fiducial thermal velocity. This description of the profile function allows for a depth-independent frequency grid (see also Paper I).

Finally, we express the continuum and line opacities in terms of the Thomson-opacity, χ_Th = n_eσ_e, and depth-independent input parameters, k_C, k_L:

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{C}}& =k_{\mathrm{C}}·{\chi }_{\mathrm{Th}}\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill {\chi }_{\mathrm{L}}& =k_{\mathrm{L}}·{\chi }_{\mathrm{Th}}·{\mathrm{\Phi }}_x=:{\overline{\chi }}_{\mathrm{L}}{\mathrm{\Phi }}_x,\hfill \end{array}$$(8)

where the line-strength k_L is related to the frequency integrated opacity ${\overline{\chi }}_0=\mathrm{\Delta }{\nu }_{\text{D}}^{*}{\overline{\chi }}_{\mathrm{L}}=k_{\text{L}}\mathrm{\Delta }{\nu }_{\text{D}}^{*}{\chi }_{\text{Th}}$ (e.g. Paper I, Eqs. (10) and (12), and with $[{\overline{\chi }}_0]$ = (cm s)⁻¹). In the following, we present (efficient) numerical tools for solving the coupled Eqs. (1)–(3). Since a direct solution is computationally prohibitive for 3D atmospheres due to limited memory capacity, we apply an accelerated Λ-iteration (ALI) scheme that overcomes the convergence problems of the classical Λ-iteration for optically thick, scattering dominated atmospheres by using an appropriate approximate Λ-operator (ALO).

3. Numerical methods

The most time-consuming part of the ALI scheme consists of calculating the formal solution. In Paper I, we have shown that a formal solution obtained using the FVM suffers from various numerical inaccuracies related to numerical diffusion and to the order of accuracy, the latter influencing the solution particularly in the optically thick regime. To avoid these errors, we implement an integral method along short characteristics. When compared with a long-characteristics solution scheme, the computation time becomes reduced by roughly a factor of N/2, with N the number of spatial grid points (per dimension), since the LC method integrates the equation of radiative transfer along the complete path from the boundary to a considered grid point. Thus, LC methods become feasible only on massively parallelized architectures.

We follow the same approach as in Paper I, and solve the equation of radiative transfer on a non-uniform, 3-dimensional Cartesian grid. In contrast to curvilinear coordinate systems, the direction vector n then becomes constant with respect to the spatial grid, thus avoiding the otherwise required angular interpolation of upwind intensities and a complicated bookkeeping of intensities (and corresponding integration weights) for different directions. Furthermore, the sweep through the spatial domain is considerably simplified, and can be performed grid point by grid point along the x, y, and z coordinates, since the intensities (required for the upwind interpolation) are always known on the previous grid layer. To enable a straightforward implementation of non-monotonic velocity fields, we use the observer’s frame formulation.

SC methods have been successfully implemented already for 3D non-LTE (NLTE) radiative transfer problems in cool stars (e.g. Vath 1994; Leenaarts & Carlsson 2009; Hayek et al. 2010; Holzreuter & Solanki 2012). These codes, however, are mostly designed for planar geometries, and only account for subsonic and slightly supersonic velocity fields. For scattering problems including highly supersonic velocity fields, there exist, to our knowledge, only the 2D codes by Dullemond & Turolla (2000; planar/spherical), van Noort et al. (2002; planar/spherical/cylindrical), Georgiev et al. (2006) and Zsargó et al. (2006; spherical). The only 3D SC code including arbitrary velocity fields, IRIS (Ibgui et al. 2013), has also been formulated for planar geometries, and lacks the implementation of a Λ-iteration scheme thus far. As has been shown in all these studies, the final performance of the SC method crucially depends on the choice of the applied interpolation schemes.

3.1. The discretized radiative transfer equation along a ray

The equation of radiative transfer along a given direction can be written as

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}\tau }}=S-I,\end{array}$$(9)

where dτ := χds is the optical-depth increment along a ray segment ds. Here and in the following, we suppress the notation for the frequency dependence, and explicitly distinguish between continuum and line only when appropriate. Equation (9) is integrated along a ray propagating through a current grid point p with Cartesian grid indices (ijk) and corresponding upwind point u^(ijk). The geometry for a 3D Cartesian grid is shown in the upper panel of Fig. 1. In the following, upwind and downwind quantities corresponding to a considered grid point (ijk) are indicated by $q_{\text{u}}^{(\mathit{ijk})}$ and $q_{\text{d}}^{(\mathit{ijk})}$, while local quantities are denoted either as $q_{\text{p}}^{(\mathit{ijk})}$ or simply q_ijk. For a given ray segment, we then obtain:

$$\begin{array}{cc}\hfill I_{\mathit{ijk}}& =I_{\mathrm{u}}^{(\mathit{ijk})}{\mathrm{e}}^{-({\tau }_{\mathrm{p}}-{\tau }_{\mathrm{u}})}+{\displaystyle {\int }_{{\tau }_{\mathrm{u}}}^{{\tau }_{\mathrm{p}}}{\mathrm{e}}^{-({\tau }_{\mathrm{p}}-\tau )}S(\tau )\mathrm{d}\tau }\hfill \\ \hfill & =I_{\mathrm{u}}^{(\mathit{ijk})}{\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}+{\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}{\displaystyle {\int }_0^{\mathrm{\Delta }{\tau }_{\mathrm{u}}}{\mathrm{e}}^{t}S(t+{\tau }_{\mathrm{u}})\mathrm{d}t,}\hfill \end{array}$$(10)

Fig. 1. Geometry of the SC method for a particular ray with direction n propagating from the upwind point u(ijk) to a considered grid point p(ijk). The downwind point d(ijk) is required to set the slope of a Bézier curve representing the opacities and source functions along the ray. The middle and lower panel display all possible downwind and upwind intersection surfaces for a short characteristic at a grid point p(ijk). For rays intersecting the xy-, xz-, and yz-planes, the 2D Bézier interpolation is obtained from given quantities at grid points located in the cyan, red, and magenta shaded surfaces, respectively. The coordinate system is indicated at the upper left, where α, β, γ determine the direction of the coordinate-axes and are defined in Sect. 3.3.

with upwind optical-depth increment Δτ_u := τ_p − τ_u ≥ 0, and t := τ − τ_u. For the SC solution scheme, the location of the reference point, τ = 0, plays no role, since only the optical-depth increments, Δτ_u and Δτ_d (see below), are required. To calculate the source contribution, the source function is commonly approximated by first- or second-order polynomials (Kunasz & Auer 1988; van Noort et al. 2002), Bézier curves (Hayek et al. 2010; Holzreuter & Solanki 2012; Auer 2003) or cubic Hermite splines (Ibgui et al. 2013). While the 2nd- (and higher) order methods reproduce the diffusion limit in optically thick media, they suffer from overshoots and need to be monotonized with some effort to ensure that any interpolated quantity remains positive between two given grid points. Monotonicity is usually obtained by manipulating the interpolation scheme whenever overshoots occur. Thus, the actual interpolation crucially depends on the specific stratification of the considered quantity (e.g. the source function). The Λ-operator then becomes non-linear, because its elements now explicitly depend on the stratification of source functions (via corresponding interpolation/integration coefficients). Within any Λ-iteration scheme, this non-linearity can lead to oscillations. In extreme cases, “flip-flop situations” (Holzreuter & Solanki 2012, their Appendix A) may occur, which do not converge at all.

For the source contribution, we implement both a linear approximation as the fastest and most stable method (monotonicity is always provided), and a quadratic Bézier approximation (see Appendix B) for higher accuracy¹, which allows us to preserve monotonicity in a rather simple way. The Bézier interpolation is constructed from two given data points and one control point, the latter setting the slope of the interpolating curve. The control point is located at the centre of the data-points abscissae, with the ordinate calculated by accounting for the information of a third data point to yield the parabola intersecting all three data points. Whenever overshoots occur, the value of the control point will be manipulated to ensure monotonicity (see Fig. B.2). The corresponding formulation is given in Appendix B, Eqs. (B.7)–(B.10). Applying these equations to describe the behaviour of the mean intensities along the optical path, and identifying the indices (i − 1), (i), (i + 1) with the upwind, current, and downwind points, we find, after reordering for the t⁰, t¹, t² terms:

$$\begin{array}{cc}\hfill S(t+{\tau }_{\mathrm{u}})=& S_{\mathrm{u}}^{(\mathit{ijk})}+[{\displaystyle \frac{(\omega -2)}{\mathrm{\Delta }{\tau }_{\mathrm{u}}}}{S}_{\mathrm{u}}^{(\mathit{ijk})}\hfill \\ \hfill & +{\displaystyle \frac{(1-\omega )\mathrm{\Delta }{\tau }_{\mathrm{u}}+(2-\omega )\mathrm{\Delta }{\tau }_{\mathrm{d}}}{\mathrm{\Delta }{\tau }_{\mathrm{u}}\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{S}_{\mathrm{p}}^{(\mathit{ijk})}+{\displaystyle \frac{(\omega -1)}{\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{S}_{\mathrm{d}}^{(\mathit{ijk})}]·t\hfill \\ \hfill & +[{\displaystyle \frac{(1-\omega )}{\mathrm{\Delta }{\tau }_{\mathrm{u}}^2}}{S}_{\mathrm{u}}^{(\mathit{ijk})}+{\displaystyle \frac{(\omega -1)(\mathrm{\Delta }{\tau }_{\mathrm{u}}+\mathrm{\Delta }{\tau }_{\mathrm{d}})}{\mathrm{\Delta }{\tau }_{\mathrm{u}}^2\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{S}_{\mathrm{p}}^{(\mathit{ijk})}\hfill \\ \hfill & +{\displaystyle \frac{(1-\omega )}{\mathrm{\Delta }{\tau }_{\mathrm{u}}\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{S}_{\mathrm{d}}^{(\mathit{ijk})}]·t^2,\hfill \end{array}$$(11)

with downwind optical-depth increment, Δτ_d = τ_d − τ_p ≥ 0. The parameter ω defines the ordinate of the control point (Eq. (B.4)). Within the Bézier interpolation, we emphasize that ω may explicitly depend on $S_{\text{u}}^{(\mathit{ijk})}$, $S_{\text{p}}^{(\mathit{ijk})}$, and $S_{\text{d}}^{(\mathit{ijk})}$ to ensure monotonicity, and not solely on the grid spacing. A major advantage of this parameterization is that we can globally define a minimum allowed ω that can be adapted during the iteration process. The flip-flop situations discussed above can then be avoided by gradually increasing ω_min towards unity (ω ≡ 1 corresponds to a linear interpolation), that is, by suppressing the curvature of the Bézier interpolation. This way, we can construct an always-convergent iteration scheme, though with the drawback of using less accurate interpolations.

Integrating Eq. (10) together with a source function described by Eq. (11), we obtain the discretized equation of radiative transfer:

$$\begin{array}{c}\hfill I_{\mathit{ijk}}=a_{\mathit{ijk}}{S}_{\mathrm{u}}^{(\mathit{ijk})}+b_{\mathit{ijk}}{S}_{\mathrm{p}}^{(\mathit{ijk})}+c_{\mathit{ijk}}{S}_{\mathrm{d}}^{(\mathit{ijk})}+d_{\mathit{ijk}}{I}_{\mathrm{u}}^{(\mathit{ijk})},\end{array}$$(12)

with

$$\begin{array}{cc}& a_{\mathit{ijk}}:=e_0+{\displaystyle \frac{\omega -2}{\mathrm{\Delta }{\tau }_{\mathrm{u}}}}{e}_1+{\displaystyle \frac{1-\omega }{\mathrm{\Delta }{\tau }_{\mathrm{u}}^2}}{e}_2\hfill \\ \hfill & b_{\mathit{ijk}}:={\displaystyle \frac{(1-\omega )\mathrm{\Delta }{\tau }_{\mathrm{u}}+(2-\omega )\mathrm{\Delta }{\tau }_{\mathrm{d}}}{\mathrm{\Delta }{\tau }_{\mathrm{u}}\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{e}_1+{\displaystyle \frac{(\omega -1)(\mathrm{\Delta }{\tau }_{\mathrm{u}}+\mathrm{\Delta }{\tau }_{\mathrm{d}})}{\mathrm{\Delta }{\tau }_{\mathrm{u}}^2\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{e}_2\hfill \\ \hfill & c_{\mathit{ijk}}:={\displaystyle \frac{\omega -1}{\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{e}_1+{\displaystyle \frac{1-\omega }{\mathrm{\Delta }{\tau }_{\mathrm{u}}\mathrm{\Delta }{\tau }_{\mathrm{d}}}}{e}_2\hfill \\ \hfill & d_{\mathit{ijk}}:={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}\hfill \\ \hfill & e_0:={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}{\displaystyle {\int }_0^{\mathrm{\Delta }{\tau }_{\mathrm{u}}}{\mathrm{e}}^t\mathrm{d}t=1-{\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}}\hfill \\ \hfill & e_1:={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}{\displaystyle {\int }_0^{\mathrm{\Delta }{\tau }_{\mathrm{u}}}t{\mathrm{e}}^t\mathrm{d}t=\mathrm{\Delta }{\tau }_{\mathrm{u}}-e_0}\hfill \\ \hfill & e_2:={\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\mathrm{u}}}{\displaystyle {\int }_0^{\mathrm{\Delta }{\tau }_{\mathrm{u}}}{t}^2{\mathrm{e}}^t\mathrm{d}t=\mathrm{\Delta }{\tau }_{\mathrm{u}}^2-2e_1.}\hfill \end{array}$$

The calculation of the upwind and downwind Δτ-steps proceeds similarly, where now the opacity is integrated using the Bézier interpolation. Using Eqs. (B.3), (B.4) for the upwind interval, and Eqs. (B.11), (B.12) for the downwind interval, one easily obtains:

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\tau }_{\mathrm{u}}& ={\displaystyle {\int }_{\mathrm{u}}^{\mathrm{p}}\chi (s)\mathrm{d}s={\displaystyle \frac{\mathrm{\Delta }{s}_{\mathrm{u}}}{3}}({\chi }_{\mathrm{u}}+{\chi }_{\mathrm{c}}^{[\mathrm{u},\mathrm{p}]}+{\chi }_{\mathrm{p}})}\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill \mathrm{\Delta }{\tau }_{\mathrm{d}}& ={\displaystyle {\int }_{\mathrm{p}}^{\mathrm{d}}\chi (s)\mathrm{d}s={\displaystyle \frac{\mathrm{\Delta }{s}_{\mathrm{d}}}{3}}({\chi }_{\mathrm{p}}+{\chi }_{\mathrm{c}}^{[\mathrm{p},\mathrm{d}]}+{\chi }_{\mathrm{d}}),}\hfill \end{array}$$(14)

where Δs_u, Δs_d describe the path lengths of the upwind and downwind intervals, respectively, and ${\chi }_{\text{c}}^{[\text{u},\text{p}]}$, ${\chi }_{\text{c}}^{[\text{p},\text{d}]}$ refer to the opacity at the control points in each interval.

3.2. Grid refinement

Since the opacity of a line transition depends on the velocity field through the Doppler effect, regions of significant opacity may become spatially confined in a highly supersonic wind with strong acceleration. Thus, a grid refinement along the short characteristic might be required to correctly account for all so-called resonance zones. Because the profile function is approximated by a Doppler profile and rapidly vanishes for |x_cmf/δ|≳3, a resonance zone is here defined by a region where x_cmf/δ ∈ [ − 3, 3].

A numerically sufficient condition to resolve all such resonance zones along a given ray is to demand that |Δx_cmf|/δ = |ΔV_proj|/δ ≲ 1/3 if a resonance zone lies in between the points [u^(ijk), p], where ΔV_proj is the projected velocity step along the ray in units of $v_{\text{th}}^{\ast }$. Assuming a linear dependence of the projected velocities on the ray coordinate s, this condition directly translates to an equidistant refined spatial grid along the ray. For short ray segments (as is mostly the case within our calculations), neglecting the second-order (curvature) terms of the projected velocity influences the solution only weakly. The line source function on the refined grid is obtained by Bézier interpolation in s-space (Eqs. (B.7)–(B.9)):

$$\begin{array}{c}\hfill S_{\mathrm{L}}(s_{\ell })=S_{\ell }={\stackrel{\sim }{a}}_{\ell }{S}_{\mathrm{u}}^{(\mathit{ijk})}+{\stackrel{\sim }{b}}_{\ell }{S}_{\mathrm{p}}^{(\mathit{ijk})}+{\stackrel{\sim }{c}}_{\ell }{S}_{\mathrm{d}}^{(\mathit{ijk})},\end{array}$$(15)

where the index ℓ refers to the points on the refined grid, and u^(ijk), p^(ijk), d^(ijk) describe the original geometry of the short characteristic. ${\overline{\chi }}_{\mathrm{L}}$ is obtained analogously, and the required Δτ_ℓ steps are calculated with the trapezoidal rule, for simplicity. Contrasted to the Sobolev method (which also assumes a linear velocity law along the ray segment, e.g. Rybicki & Hummer 1978), our grid refinement procedure explicitly accounts for variations of the opacity and the source function.

Using Eq. (12) for the inter-grid points, such that the (local) upwind, current, and downwind quantities are now described by the indices [ℓ − 1, ℓ, ℓ + 1], we obtain:

$$\begin{array}{cc}\hfill I_{\ell }& =I_{\ell -1}{\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\ell }}+(a_{\ell }{\stackrel{\sim }{a}}_{\ell -1}+b_{\ell }{\stackrel{\sim }{a}}_{\ell }+c_{\ell }{\stackrel{\sim }{a}}_{\ell +1})S_{\mathrm{u}}^{(\mathit{ijk})}\hfill \end{array}$$(16)

$$\begin{array}{cc}& +(a_{\ell }{\stackrel{\sim }{b}}_{\ell -1}+b_{\ell }{\stackrel{\sim }{b}}_{\ell }+c_{\ell }{\stackrel{\sim }{b}}_{\ell +1})S_{\mathrm{p}}^{(\mathit{ijk})}\hfill \end{array}$$(17)

$$\begin{array}{cc}& +(a_{\ell }{\stackrel{\sim }{c}}_{\ell -1}+b_{\ell }{\stackrel{\sim }{c}}_{\ell }+c_{\ell }{\stackrel{\sim }{c}}_{\ell +1})S_{\mathrm{d}}^{(\mathit{ijk})}\hfill \end{array}$$(18)

$$\begin{array}{cc}& :=I_{\ell -1}{\mathrm{e}}^{-\mathrm{\Delta }{\tau }_{\ell }}+{\stackrel{\sim }{\alpha }}_{\ell }{S}_{\mathrm{u}}^{(\mathit{ijk})}+{\stackrel{\sim }{\beta }}_{\ell }{S}_{\mathrm{p}}^{(\mathit{ijk})}+{\stackrel{\sim }{\gamma }}_{\ell }{S}_{\mathrm{d}}^{(\mathit{ijk})}.\hfill \end{array}$$(19)

For a number of N_ref refinement points (including the upwind and current point) within the interval [u^(ijk), p^(ijk)], the intensity at point p^(ijk) is finally given by:

$$\begin{array}{cc}\hfill I_{\mathit{ijk}}=& I_{\mathrm{u}}^{(\mathit{ijk})}{\mathrm{e}}^{-{\sum }_{m=2}^{N_{\mathrm{ref}}}\mathrm{\Delta }{\tau }_m}+S_{\mathrm{u}}^{(\mathit{ijk})}{\displaystyle \sum _{m=2}^{N_{\mathrm{ref}}}{\stackrel{\sim }{\alpha }}_m{\mathrm{e}}^{-{\sum }_{n=m+1}^{N_{\mathrm{ref}}}\mathrm{\Delta }{\tau }_n}}\hfill \\ \hfill & +S_{\mathrm{p}}^{(\mathit{ijk})}{\displaystyle \sum _{m=2}^{N_{\mathrm{ref}}}{\stackrel{\sim }{\beta }}_m{\mathrm{e}}^{-{\sum }_{n=m+1}^{N_{\mathrm{ref}}}\mathrm{\Delta }{\tau }_n}+S_{\mathrm{d}}^{(\mathit{ijk})}\sum _{m=2}^{N_{\mathrm{ref}}}{\stackrel{\sim }{\gamma }}_m{\mathrm{e}}^{-{\sum }_{n=m+1}^{N_{\mathrm{ref}}}\mathrm{\Delta }{\tau }_n},}\hfill \end{array}$$(20)

where the upwind and current points always correspond to the indices m = 1 and m = N_ref, respectively, and the sum over m is performed over N_ref − 1 intervals. The discretized radiative transfer equation for the refined grid obviously has the same form as for the standard short characteristic (Eq. (12)), with different coefficients though.

3.3. Upwind and downwind interpolations

To solve the discretized equation of radiative transfer, the opacities χ_C(u, d), ${\overline{\chi }}_{\mathrm{L}(\text{u},\text{d})}$, source functions S_C(u, d), S_L(u, d), and velocity vectors V_(u, d), are required at the upwind and downwind points, together with the incident intensity, I_u. We emphasize that the subscript C describes continuum quantities, and should not be confused with the subscript c denoting the control points of the interpolation scheme.

All required quantities are obtained from a 2D Bézier interpolation (see Appendix C) on the surfaces that intersect with a given ray. The intersection surfaces depend on the considered direction and the size of the upwind and downwind grid cells. For a given direction

$$\begin{array}{c}\hfill \mathit{n}=\left(\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right)=\left(\begin{array}{c}sin\theta cos\varphi \\ sin\theta sin\varphi \\ cos\theta \end{array}\right),\end{array}$$(21)

where θ is the co-latitude (measured from the Cartesian z-axis), and ϕ is the azimuth (measured from the x-axis), the distances from a considered grid point to the neighbouring xy-, xz-, and yz-planes are calculated from trigonometry and yield:

$$\begin{array}{cc}& \mathrm{\Delta }{s}_{\mathit{xy}}^{(\mathrm{u})}={\displaystyle \frac{z_k-z_{k-\gamma }}{n_z}}\mathrm{\Delta }{s}_{\mathit{xy}}^{(\mathrm{d})}={\displaystyle \frac{z_{k+\gamma }-z_k}{n_z}}\hfill \\ \hfill & \mathrm{\Delta }{s}_{\mathit{xz}}^{(\mathrm{u})}={\displaystyle \frac{y_j-y_{j-\beta }}{n_y}}\mathrm{\Delta }{s}_{\mathit{xz}}^{(\mathrm{d})}={\displaystyle \frac{y_{j+\beta }-y_j}{n_y}}\hfill \\ \hfill & \mathrm{\Delta }{s}_{\mathit{yz}}^{(\mathrm{u})}={\displaystyle \frac{x_i-x_{i-\alpha }}{n_x}}\mathrm{\Delta }{s}_{\mathit{yz}}^{(\mathrm{d})}={\displaystyle \frac{x_{i+\alpha }-x_i}{n_x}},\hfill \end{array}$$

with α, β, γ set to ±1 for direction-vector components n_x, n_y, n_z ≷ 0, respectively. The intersection surface on the upwind and downwind side are then found at the minimum of $\mathrm{\Delta }{s}_{\mathit{xy}}^{(\text{u},\text{d})},\mathrm{\Delta }{s}_{\mathit{xz}}^{(\text{u},\text{d})},\mathrm{\Delta }{s}_{\mathit{yz}}^{(\text{u},\text{d})}$, and the corresponding coordinates are easily calculated.

For each surface, the interpolation requires nine points within the corresponding plane (see Fig. 1 and Eq. (C.1)). In each considered plane, we generally use grid points running from (index-2) to (index) to determine upwind quantities, while downwind quantities are calculated from (index-1) to (index+1). Such a formulation greatly simplifies the calculation of the Λ-matrix elements (see Appendix D). In Fig. 1, we show an example for a ray intersecting the xy-plane at the upwind side. The 2D Bézier interpolation for the upwind point then consists of three 1D Bézier interpolations along the x-axis using the points (J_u, K_u, L_u), (D_u, E_u, F_u), (M_u, N_u, O_u), followed by another 1D Bézier interpolation along the y-axis at the upwind x-coordinates. With the 2D Bézier interpolation given by Eq. (C.1), we find for each required quantity q_u, d:

$$\begin{array}{cc}\hfill q_{\mathrm{u}}^{(\mathit{ijk})}=& w_{\mathrm{A}}^{(\mathit{ijk})}{q}_{i-2\alpha ,j-\beta ,k-2\gamma }+w_{\mathrm{B}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-\beta ,k-2\gamma }+w_{\mathrm{C}}^{(\mathit{ijk})}{q}_{i,j-\beta ,k-2\gamma }\hfill \\ \hfill & +w_{\mathrm{D}}^{(\mathit{ijk})}{q}_{i-2\alpha ,j-\beta ,k-\gamma }+w_{\mathrm{E}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-\beta ,k-\gamma }+w_{\mathrm{F}}^{(\mathit{ijk})}{q}_{i,j-\beta ,k-\gamma }\hfill \\ \hfill & +w_{\mathrm{G}}^{(\mathit{ijk})}{q}_{i-2\alpha ,j-\beta ,k}+w_{\mathrm{H}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-\beta ,k}+w_{\mathrm{I}}^{(\mathit{ijk})}{q}_{i,j-\beta ,k}\hfill \\ \hfill & +w_{\mathrm{J}}^{(\mathit{ijk})}{q}_{i-2\alpha ,j-2\beta ,k-\gamma }+w_{\mathrm{K}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-2\beta ,k-\gamma }+w_{\mathrm{L}}^{(\mathit{ijk})}{q}_{i,j-2\beta ,k-\gamma }\hfill \\ \hfill & +w_{\mathrm{M}}^{(\mathit{ijk})}{q}_{i-2\alpha ,j,k-\gamma }+w_{\mathrm{N}}^{(\mathit{ijk})}{q}_{i-\alpha ,j,k-\gamma }+w_{\mathrm{O}}^{(\mathit{ijk})}{q}_{i,j,k-\gamma }\hfill \\ \hfill & +w_{\mathrm{P}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-2\beta ,k-2\gamma }+w_{\mathrm{Q}}^{(\mathit{ijk})}{q}_{i-\alpha ,j,k-2\gamma }\hfill \\ \hfill & +w_{\mathrm{R}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-2\beta ,k}+w_{\mathrm{S}}^{(\mathit{ijk})}{q}_{i-\alpha ,j,k}+w_{\mathit{ijk}}{q}_{\mathit{ijk}}\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill q_{\mathrm{d}}^{(\mathit{ijk})}=& {\stackrel{\sim }{w}}_{\mathrm{A}}^{(\mathit{ijk})}{q}_{i-\alpha ,j+\beta ,k-\gamma }+{\stackrel{\sim }{w}}_B^{(\mathit{ijk})}{q}_{i,j+\beta ,k-\gamma }+{\stackrel{\sim }{w}}_{\mathrm{C}}^{(\mathit{ijk})}{q}_{i+\alpha ,j+\beta ,k-\gamma }\hfill \\ \hfill & +{\stackrel{\sim }{w}}_{\mathrm{D}}^{(\mathit{ijk})}{q}_{i-\alpha ,j+\beta ,k}+{\stackrel{\sim }{w}}_{\mathrm{E}}^{(\mathit{ijk})}{q}_{i,j+\beta ,k}+{\stackrel{\sim }{w}}_{\mathrm{F}}^{(\mathit{ijk})}{q}_{i+\alpha ,j+\beta ,k}\hfill \\ \hfill & +{\stackrel{\sim }{w}}_{\mathrm{G}}^{(\mathit{ijk})}{q}_{i-\alpha ,j+\beta ,k+\gamma }+{\stackrel{\sim }{w}}_{\mathrm{H}}^{(\mathit{ijk})}{q}_{i,j+\beta ,k+\gamma }+{\stackrel{\sim }{w}}_{\mathrm{I}}^{(\mathit{ijk})}{q}_{i+\alpha ,j+\beta ,k+\gamma }\hfill \\ \hfill & +{\stackrel{\sim }{w}}_{\mathrm{J}}^{(\mathit{ijk})}{q}_{i-\alpha ,j-\beta ,k+\gamma }+{\stackrel{\sim }{w}}_{\mathrm{K}}^{(\mathit{ijk})}{q}_{i,j-\beta ,k+\gamma }+{\stackrel{\sim }{w}}_{\mathrm{L}}^{(\mathit{ijk})}{q}_{i+\alpha ,j-\beta ,k+\gamma }\hfill \\ \hfill & +{\stackrel{\sim }{w}}_{\mathrm{M}}^{(\mathit{ijk})}{q}_{i-\alpha ,j,k+\gamma }+{\stackrel{\sim }{w}}_{\mathrm{N}}^{(\mathit{ijk})}{q}_{i,j,k+\gamma }+{\stackrel{\sim }{w}}_{\mathrm{O}}^{(\mathit{ijk})}{q}_{i+\alpha ,j,k+\gamma }\hfill \\ \hfill & +{\stackrel{\sim }{w}}_{\mathrm{P}}^{(\mathit{ijk})}{q}_{i+\alpha ,j-\beta ,k-\gamma }+{\stackrel{\sim }{w}}_{\mathrm{Q}}^{(\mathit{ijk})}{q}_{i+\alpha ,j,k-\gamma }\hfill \\ \hfill & +{\stackrel{\sim }{w}}_{\mathrm{R}}^{(\mathit{ijk})}{q}_{i+\alpha ,j-\beta ,k}+{\stackrel{\sim }{w}}_{\mathrm{S}}^{(\mathit{ijk})}{q}_{i+\alpha ,j,k},\hfill \end{array}$$(23)

where the coefficients w^(ijk) and ${\stackrel{\sim }{w}}^{(\mathit{ijk})}$ refer to the upwind and downwind interpolations corresponding to a considered point (ijk). Depending on the intersection surface, ten out of these 19 coefficients are set to zero. For the upwind interpolation, we have already included the local coefficient (ijk), which is only required when boundary conditions need to be specified (Sect. 3.4). We note that all (non-zero) interpolation coefficients may depend on the specific values of a considered quantity at the given grid points, via the interpolation parameter ω to ensure monotonicity. As in Sect. 3.1, also these monotonicity constraints result in non-linear Λ-operators.

3.4. Boundary conditions

Since the inner boundary is usually not aligned with the 3D Cartesian grid (e.g. a spherical star at the origin), the upwind (and downwind) interpolations need to be adapted near the stellar surface. For the upwind point, the following two situations may occur (see Fig. 2 for an example in the xz-plane). Firstly, the considered ray originates from the stellar surface (direction n₁ in Fig. 2). In this case, we use a core-halo approximation and set $I_{\text{u}}=I_{\text{c}}^{+}=B_{\nu }(T_{\text{rad}})$, with $I_{\text{c}}^{+}$ the emergent intensity from the core, and T_rad the radiation temperature. Unless explicitly noted, we assume T_rad = T_eff throughout this work. All other quantities are obtained from trilinear interpolation using the points (E_u, F_u, H_u, I_u, N_u, O_u, S_u, p^(ijk)) in Fig. 1, where representative estimates need to be defined at the core points (those points that are located inside the star). In hydrodynamic simulations, the analogue of these points are so-called “ghost points”. Secondly, the considered ray originates from a plane spanned by grid points that are partially located inside the star (direction n₂ in Fig. 2). Then, the interpolation is performed as in Sect. 3.3, using again representative estimates at the core-points.

Fig. 2. Boundary conditions for rays propagating in the xz-plane at y-level (j) with three different directions n1, n2, n3, and upwind points u1, u2, u3. For point u1, the intensity is set to $I_{\text{u}}=I_{\text{c}}^{+}$ and all remaining quantities are obtained by bilinear interpolation from points (i − 1, j, k − 1), (i, j, k − 1), (i − 1, j, k), and (ijk). The required quantities at point u2 are found from Bézier interpolation using the values at (i − 2, j, k − 1), (i − 1, j, k − 1), (i, j, k − 1). The (unknown) quantities inside the core are indicated by red dots, and need to be reasonably approximated (see text). For direction n3, the unknown intensity inside the core, $I_{\text{c}}^{-}$, is directed inwards. Such situations occur only for ray directions (nearly) parallel to the spatial grid, and thus are relatively seldom.

Inside the core, we define $I_{\text{c}}^{+}=S_{\text{L}}=S_{\text{C}}=B_{\nu }(T_{\text{rad}})$ and set $I_{\text{c}}^{-}$ and all velocity components to zero, where $I_{\text{c}}^{-}$ is the inward directed intensity which needs to be specified only in rare situations (Fig. 2, direction n₃). The opacities inside the star are found by extrapolation from the known values outside the star. While this procedure certainly introduces errors (e.g. by over- and underestimating the upwind source function in optically thin and thick cases, respectively), it is still favourable to extrapolating all values directly onto the stellar surface, mainly due to performance reasons². In addition to the error introduced by the predefined values inside the core, the calculation of Δs_r is a certain issue, where Δs_r is the distance of the current grid point to the stellar surface. Since the radiative transfer near the stellar surface is (in most cases) very sensitive to the path length of a considered ray, Δs_r needs to be known exactly. Depending on the shape of the surface, Δs_r can be calculated analytically, or needs to be determined numerically. A numerical solution, however, might be time consuming and should be avoided when possible. Downwind quantities are always calculated from Eq. (23), using the estimates at the core points as defined above when necessary.

3.5. Angular and frequency integration

To obtain the mean intensity at each grid point in the atmosphere, we solve the discretized equation of radiative transfer for many directions and numerically integrate via:

$$\begin{array}{c}\hfill J_{\mathit{ijk}}={\displaystyle \frac{1}{4\pi }}{\displaystyle \int I_{\mathit{ijk}}\mathrm{d}\mathrm{\Omega }=\sum _l{w}_l{I}_{\mathit{ijk}}({\mathrm{\Omega }}_l),}\end{array}$$(24)

where w_l is the integration weight corresponding to a considered direction Ω_l = (θ_l, ϕ_l). The angular integration is particular challenging for optically thin atmospheres, since in such situations each (spatial) grid point is illuminated by the stellar surface, and the distribution of intensities I_ijk(θ, ϕ) becomes a 2D step-function in the θ − ϕ-plane (if no upwind interpolation errors were present). Depending on the considered position, the shape of I_ijk(θ, ϕ) greatly varies. Thus, elaborate integration methods are required to resolve the star and its edges at any point of the atmosphere.

Lobel & Blomme (2008) use the trapezoidal rule with a decreasing number of polar grid points at higher latitudes to reasonably distribute the direction vectors on the unit sphere. For the 3D FVM, we have shown in Paper I that a Gauss-Legendre integration performs (slightly) better. However, the corresponding directions are always clustered in certain regions since the nodes of the Gauss-Legendre quadrature are fixed. Additionally, the Gauss-Legendre integration should only be applied when the distribution of intensities can be described by high order polynomials, that is, when I_ijk(θ, ϕ) is smoothed out (e.g. by numerical diffusion). When numerical diffusion errors are suppressed (e.g. by using elaborate upwind interpolation schemes), the Gauss-Legendre integration should not be used (see Table 1).

We have tested a multitude of other quadrature schemes, including trapezoidal and (pseudo)-Gaussian rules on triangles, and the so-called Lebedev quadrature (see, e.g. Ahrens & Beylkin 2009; Beentjes 2015, and references therein). The Lebedev quadrature is optimized to exactly integrate the spherical harmonics up to a certain degree, with a (nearly) optimum distribution of direction vectors on the unit sphere. In Table 1, we summarize the errors for an optically thin atmosphere using different integration methods. The incident intensities have been obtained exactly (i.e. by setting I_ijk = I_c and I_ijk = 0 for core and non-core rays, respectively), or from the 3D SC method using linear interpolations. Considering the SC solution scheme, the solution has only been slightly improved (if at all) when doubling the angular grid resolution from N_Ω ≈ 1000 to N_Ω ≈ 2000, for all applied integration methods. We note that the mean relative error does not converge to zero due to the upwind interpolation scheme (see Sect. 4.1). For the exact solution of the optically thin radiative transfer, the Lebedev-integration method performs best, and is therefore used within all our calculations. When calculating line transitions, the location of resonance zones depends on the considered direction. Thus, we generally use N_Ω = 2030 direction vectors to ensure that no resonance zone has been overlooked. The corresponding angular resolution is typical when calculating 3D radiative transfer problems in extended stellar atmospheres. For instance, Lobel & Blomme (2008) used N_Ω = 6400 within their 3D finite-volume method.

To obtain the scattering integral, we apply the trapezoidal rule for the frequency integration. The scattering integral then reads:

$$\begin{array}{c}\hfill {\overline{J}}_{\mathit{ijk}}={\displaystyle \frac{1}{4\pi }}{\displaystyle \int \mathrm{d}\mathrm{\Omega }{\int }_{x_{\mathrm{obs}}^{(\min )}}^{x_{\mathrm{obs}}^{(\max )}}\mathrm{d}x{I}_{\mathit{ijk}}{\mathrm{\Phi }}_x^{(\mathit{ijk})}=\sum _l{w}_l\sum _x{w}_x{I}_{\mathit{ijk}}{{\mathrm{\Phi }}_x}^{(\mathit{ijk})},}\end{array}$$(25)

with $x_{\mathrm{obs}}^{(\min )}$ and $x_{\mathrm{obs}}^{(\max )}$ the required frequency shift in the observer’s frame obtained from the maximum absolute velocity occurring in the atmosphere (see also Paper I), and w_x the corresponding frequency integration weight. To resolve the profile function at each point in the atmosphere, we demand that |Δx_obs|/δ ≲ 1/3. Since the profile function depends on the ratio of fiducial to actual thermal width, the fiducial velocity should be set to the minimum thermal velocity present in the atmosphere.

3.6. Λ-iteration

In Sect. 3.1, we have already noted that the Λ-operator becomes non-linear due to monotonicity constraints (implemented by the interpolation parameter ω). In this section, we present a suitable workaround, beginning with a recapitulation of some fundamental ideas.

3.6.1. Λ-matrix elements

With the discretized equation of radiative transfer, Eqs. (12)/(20), and the upwind and downwind quantities obtained from Eqs. (22) and (23), the intensity at each spatial, angular, and frequency grid point can be calculated for a given source function. We use the standard Λ-formalism to write the formal solution of the intensity, mean intensity, and scattering integral as:

$$\begin{array}{cc}\hfill I& ={\mathrm{\Lambda }}_{\mathrm{\Omega },\nu }[S_{\mathrm{C},\mathrm{L}}]\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill J& ={\mathrm{\Lambda }}_{\nu }[S_{\mathrm{C}}]\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill \overline{J}& =\mathrm{\Lambda }[S_{\mathrm{L}}],\hfill \end{array}$$(28)

with subscripts Ω and ν defining the dependence of the Λ-operator on direction and frequency, respectively. In the following, we focus on the line transport. The continuum can be derived analogously. When all interpolation parameters ω have been determined (for a given stratification of source functions and intensities), the Λ-operator is an affine operator described by the Λ-matrix and a constant displacement vector Φ_B representing the propagation of boundary conditions (for a detailed discussion, see Paper I, Puls 1991, and references therein). The Λ-matrix elements can then be obtained by:

$$\begin{array}{c}\hfill {\mathrm{\Lambda }}_{m,n}={\overline{J}}_m({\mathit{S}}_{\mathrm{L}}={\mathit{e}}_n,{\mathbf{\Phi }}_{\mathrm{B}}=0),\end{array}$$(29)

with the nth unit vector e_n, and matrix indices m, n related to the 3D indices (i, j, k) by

$$\begin{array}{c}\hfill m=i+N_x(j-1)+N_x{N}_y(k-1),\end{array}$$(30)

where N_x and N_y denote the number of spatial grid points of the x and y coordinate, respectively. Equation (30) simply transforms a data cube to a 1D array. The m, nth matrix-element describes the effect of a non-vanishing source function at grid point n onto grid point m. We emphasize that Eq. (29) holds only for pre-calculated interpolation parameters ω, obtained from an already known stratification of source functions.

3.6.2. Accelerated Λ-iteration

The classical Λ-iteration scheme is defined by calculating a formal solution for a given source function using Eq. (28), followed by the calculation of a new source function by means of Eq. (3). For optically thick, scattering dominated atmospheres, however, this iteration scheme suffers from severe convergence problems (see Fig. 5 for the convergence behaviour of spherically symmetric test models). To overcome these problems, we apply an accelerated Λ-iteration scheme based on operator-splitting methods (Cannon 1973). Within the ALI, the Λ-operator is written as

$$\begin{array}{c}\hfill \mathrm{\Lambda }={\mathrm{\Lambda }}^{(\mathrm{A})}+(\mathrm{\Lambda }-{\mathrm{\Lambda }}^{(\mathrm{A})}),\end{array}$$(31)

where the first term is an appropriately chosen ALO acting on the new source function, $S_{\mathrm{L}}^{(k)}$, and the second term acts on the previous one, $S_{\mathrm{L}}^{(k-1)}$. For the converged solution, this scheme becomes an exact relation. Using also, and in analogy to the exact Λ-operator, an affine representation for the approximate one, Λ^(A)[S] = Λ^* ⋅ S + Φ_B (cf. above), and evaluating Λ^(A) at the previous iteration step, k − 1, we obtain:

$$\begin{array}{cc}\hfill {\mathit{S}}_{\mathrm{L}}^{(k)}& =\mathit{\zeta }·{\overline{\mathit{J}}}^{(k)}+\mathbf{\Psi }\hfill \\ \hfill & \approx \mathit{\zeta }·{\mathrm{\Lambda }}_{k-1}^{(\mathrm{A})}[{\mathit{S}}_{\mathrm{L}}^{(k)}]+\mathit{\zeta }·({\mathrm{\Lambda }}_{k-1}-{\mathrm{\Lambda }}_{k-1}^{(\mathrm{A})})[{\mathit{S}}_{\mathrm{L}}^{(k-1)}]+\mathbf{\Psi }\hfill \\ \hfill & =\mathit{\zeta }·({\mathbf{\Lambda }}_{k-1}^{\ast }{\mathit{S}}_{\mathrm{L}}^{(k)}+{\mathbf{\Phi }}_{\mathrm{B}}^{(k-1)}+{\overline{\mathit{J}}}^{(k-1)}-{\mathbf{\Lambda }}_{k-1}^{\ast }{\mathit{S}}_{\mathrm{L}}^{(k-1)}-{\mathbf{\Phi }}_{\mathrm{B}}^{(k-1)})+\mathbf{\Psi }.\hfill \\ \hfill \end{array}$$(32)

Here, the iteration indices k − 1 and k are indicated as sub- or superscripts, ζ := 1 − ϵ_L is a diagonal matrix, and Ψ := ϵ_L ⋅ B_ν(T) is the thermal contribution vector. From Eq. (32), it is obvious that the Φ_B terms cancel. For multi-level atoms, we emphasize that Λ and Λ^(A) may change within the ALI-cycle due to the variation of opacities (induced by the subsequently updated occupation numbers). Furthermore, both operators might also change even for the simplified TLA approach considered in this paper, since the corresponding matrix elements depend on the source functions via the interpolation parameters ω_k − 1, ω_k. Rearranging terms, we find:

$$\begin{array}{c}\hfill (\mathbf{1}-\mathit{\zeta }·{\mathbf{\Lambda }}_{k-1}^{\ast }){\mathit{S}}_{\mathrm{L}}^{(k)}=\mathit{\zeta }·({\overline{\mathit{J}}}^{(k-1)}-{\mathbf{\Lambda }}_{k-1}^{\ast }·{\mathit{S}}_{\mathrm{L}}^{(k-1)})+\mathbf{\Psi }.\end{array}$$(33)

Equation (33) is solved to obtain a new source function ${\mathit{S}}_{\mathrm{L}}^{(k)}$ (see below). Since, however, ${\mathbf{\Lambda }}_{k-1}^{*}$ has been optimized only to ensure monotonicity in a specific step k − 1 (based on source function ${\mathit{S}}_{\mathrm{L}}^{(k-1)}$), the iteration scheme can oscillate due to oscillations in ${\mathbf{\Lambda }}_k^{*}$ and ${\mathbf{\Lambda }}_{k-1}^{*}$. Even worse, the new source function might become negative. To overcome these problems, non-linear situations need to be avoided (by providing almost constant Λ^*-matrices over subsequent iteration steps). The following approach has proved to lead to a stable and convergent scheme: We apply purely linear interpolations (ω_k − 1 = ω_k = 1 and thus ${\mathbf{\Lambda }}_k^{*}={\mathbf{\Lambda }}_{k-1}^{*}$) in the first four iteration steps to obtain an already smooth stratification of source functions. Additionally, we globally define a minimum allowed interpolation parameter and demand that ω >  ω_min. Then, ω becomes constant (namely ω = ω_min) in (most) critical situations, and again, ${\mathbf{\Lambda }}_k^{*}$ approaches ${\mathbf{\Lambda }}_{k-1}^{*}$. Whenever negative source functions or oscillations occur within the iteration scheme, ω_min is gradually increased to one. With this approach, we obtain an always convergent iteration scheme, with a formal solution obtained by using linear interpolations only in most challenging cases.

3.6.3. Constructing the ALO

The rate of convergence achieved by the accelerated Λ-iteration scheme increases with the number of Λ-matrix elements included in the ALO. To minimize the computation time of the complete procedure, the choice of the ALO is always a compromise between the number of matrix-elements to be calculated, and the resulting convergence speed. While the inversion of a diagonal ALO reduces to a simple division, a multi-band ALO needs to be inverted with some more effort. When taking only the nearest neighbours into account, however, the ALO becomes a sparse matrix, and can be efficiently inverted by applying the Jacobi-Iteration for sparse systems (see Paper I). For 3D calculations, a multi-band ALO is required to obtain a rapidly converging iteration scheme (see Hauschildt & Baron 2006 and Paper I for solutions obtained with the LC method and the FVM, respectively), and thus implemented within our 3D SC framework. To calculate the corresponding Λ-matrix elements (including all upwind and downwind interpolations), we extend the procedure developed by Olson & Kunasz (1987) and Kunasz & Olson (1988). A detailed derivation is given in Appendix D. Equations (D.1)–(D.27) correspond to the exact Λ-matrix elements for a local point and its 26 neighbours, and thus should give an excellent rate of convergence when included in the ALO (see, e.g. the 26-neighbour ALO of PHOENIX/3D, Hauschildt & Baron 2006). Furthermore, all elements can be calculated in parallel to the formal solution. This property becomes important when the ALO varies during the iteration scheme, that is, when applying monotonic Bézier interpolations (as discussed above), or when accounting for multi-level atoms (for which the occupation numbers and thus opacities might change during the iteration)³.

In this paper, we analyse the convergence speed of the ALI for a diagonal ALO given by Eq. (D.14), a “direct-neighbour” (DN)-ALO given by Eqs. (D.5), (D.11), (D.13)–(D.15), (D.17), (D.23), and a “nearest-neighbour” (NN)-ALO obtained from all Eqs. (D.1)–(D.27). We note that only a moderate improvement of the computation time can be expected when using the diagonal or DN-ALO, since the diagonal and direct-neighbour elements depend on several other neighbours through the inclusion of downwind interpolations. Since, however, the downwind-integration weight is generally negative, neglecting these terms will overestimate the considered matrix elements, possibly resulting in a divergent iteration scheme. On the other hand, when using purely linear interpolations (for the source contribution and upwind interpolations), the calculation of the ALO is greatly simplified since all coefficients c_ijk and w_A, w_B, w_C, w_D, w_G, w_J, w_K, w_L, w_M, w_P, w_Q, w_R vanish. For third-order upwind/downwind interpolations as used in IRIS (Ibgui et al. 2013), the calculation of the ALO coefficients becomes computationally prohibitive at some point. Considering both interpolation techniques used in this paper, the calculation of the diagonal, DN-, and NN-ALO, in parallel to the formal solution requires 20%, 30%, and 40% of the total computation time.

Finally, we have implemented the Ng-extrapolation (Ng 1974; Olson et al. 1986) to improve the convergence speed further. As in Paper I, the Ng-acceleration is applied in every fifth iteration step in order to use independent extrapolations.

3.7. Parallelization and timing

To minimize the computation time of our 3D code, we have implemented an elaborate grid construction procedure using a non-uniform grid-spacing that still enables a reasonably high spatial resolution (see Appendix A). Furthermore, when calculating line transitions, we have parallelized the code using OPENMP as in Paper I. The parallelization is implemented over the frequency grid. We note that OPENMP creates a local copy of the 3D arrays representing the intensity and the (nearest-neighbour) Λ-matrix. With 27 Λ-matrix elements (per spatial grid point) included for the ALO calculations, the (spatial) resolution becomes therefore memory limited. A typical resolution of N_x  =  N_y  =  N_z  =  93, however, is still feasible, and gives reasonable results. For the models calculated in Sect. 4.2, typical computation times are $t_{\mathrm{SC}}^{(\mathrm{linear})}\approx 2$ h and $t_{\mathrm{SC}}^{(\mathrm{B}\acute{\mathrm{e}}\mathrm{zier})}\approx 6$ h per iteration when applying the 3D SC methods on an INTEL XEON X5650 (2.67 GHz) machine with 16 CPUs. As a reference, the FVM from Paper I “only” required roughly 50 min. A more meaningful comparison, however, is the computation time per iteration, per CPU, and per angular and frequency grid point using the same spatial grids for all methods. For an equidistant grid with N_x  =  N_y  =  N_z  =  71 grid points, we find computation times of t_FVM ≈ 0.037 s, $t_{\mathrm{SC}}^{(\mathrm{linear})}\approx 0.138$ s, and $t_{\mathrm{SC}}^{(\mathrm{B}\acute{\mathrm{e}}\mathrm{zier})}\approx 0.448$ s. Thus, the computation times of the 3D SC methods using linear/Bézier interpolations are increased by a factor of roughly four/twelve, when compared to the 3D FVM. These differences originate from the computationally more challenging upwind/downwind interpolations, the integration of the discretized equation of radiative transfer on (possibly) refined grids along a given ray, and from the calculation of an ALO including 26 neighbouring elements (instead of the 6 direct neighbours as used for the 3D FVM).

4. Spherically symmetric models

With the numerical tools developed in Sect. 3, we are able to tackle 3D continuum and line scattering problems for arbitrary velocity fields. In the following, we discuss the performance of the code when applied to spherically symmetric test models. We compare the solutions obtained from the 3D SC method using linear and Bézier interpolations (hereafter denoted by SClin and SCbez, respectively), with those obtained from the 3D FVM and from accurate 1D solvers⁴. Although we are using the same models as in Paper I, the FVM solutions might differ slightly from those presented in this paper, since we are using a different grid, optimized for the 3D SC method.

4.1. Searchlight-beam test

A first test of our 3D SC methods is the searchlight-beam test (e.g. Kunasz & Auer 1988). Within this test, we set the opacity to zero and consider the illumination of the atmosphere by a central star for a single direction. For consistency (cf. Paper I), the direction vector corresponds to θ = 45°, ϕ = 0°. Since the discretized equation of radiative transfer, Eqs. (12)/(20), reduces to $I_{\mathit{ijk}}=I_{\mathit{ijk}}^{(u)}$, this test extracts the effects of the applied interpolation schemes for the upwind intensity. The upper panel of Fig. 3 shows the propagation of the specific intensity scaled by I_c in the xz-plane. Due to the upwind interpolation, the beam emerging from the stellar core becomes widened. These interpolation errors are connected with numerical diffusion, and could only be avoided by applying a LC method. To obtain a quantitative measure of this effect, the lower and middle panel of Fig. 3 display the specific intensity along the given direction at the centre of the beam, and the specific intensity through a circular area perpendicular to the ray direction as a function of impact parameter p. The corresponding exact solutions are given by a constant and rectangular function, respectively.

Fig. 3. Searchlight-beam test for direction n = (1, 0, 1) and a typical grid with Nx = Ny = Nz = 133 grid points. Upper panel: contour plot of the specific intensity as calculated with the SC method using Bézier interpolations in the xz-plane (cf. Paper I, Fig. 3, for the finite-volume method). Middle panel: specific intensity through the perpendicular area indicated by the straight line in the upper panel. The blue, red, and green profiles correspond to the FVM, SClin, and SCbez methods, respectively. The dashed line indicates the theoretical profile. Bottom panel: as middle panel, but along the centre of the searchlight beam. We note that the SC methods reproduce the exact solution at the centre of the beam, whereas the FVM solution decreases significantly for r ≳ 2.5 R*.

Along the beam centre, the SC methods perfectly reproduce the exact solution, whereas the FVM solution decreases significantly due to the finite grid-cell size (see Paper I). Considering the intensity through the perpendicular area, both SC methods perform better than the FVM, with slight advantages of the SCbez method when compared with the SClin method. Within the 3D SC methods, however, energy conservation is violated for our zero opacity models, because the (nominal) specific intensity jumps from $I_{\mathrm{c}}^{+}$ to zero for rays intersecting the stellar surface or not, due to the core-halo approximation. As a consequence, almost all interpolations (and interpolation schemes) overestimate the specific intensity⁵. For optically thick models (where $I_{\mathrm{c}}^{+}$ at the core plays a negligible or minor role), this effect should decrease though. The associated error can be quantified by calculating the corresponding flux, that is, by integrating the specific intensity for a given direction over a corresponding perpendicular area (defined as a circle with virtually infinite radius). We emphasize that the flux as defined here constitutes the most demanding test case, and should not be confused with the flux density (i.e. the first moment of the specific intensity). Figure 4 shows the resulting fluxes (normalized by the nominal value) for searchlight beams with different directions defined by θ = 45° and ϕ ∈ [0° ,90° ]. For different directions ϕ, the searchlight beams propagate through different domains of the spatial grid with accordingly different grid-cell sizes. Due to the distinct behaviour of numerical diffusion errors within these domains, the total flux varies as a function of ϕ. Overall, the total fluxes for the SClin and SCbez methods are larger than theoretically constrained, whereas the FVM gives (despite a small error) reasonable results. This effect is largest in regions far from the star and for diagonal directions. Thus, particularly in these regions, also the mean intensities (for optically thin atmospheres) are expected to be overestimated. The same problem arises when calculating line transitions, since photons may freely propagate over large distances before a resonance region is hit. Numerical diffusion errors can only be avoided by increasing the grid resolution, or using higher order upwind-interpolation methods.

Fig. 4. Photon flux as a function of direction angle ϕ (with fixed θ = 45°) through corresponding perpendicular areas, and with the opacity set to zero. The central distance of all areas to the stellar surface has been set to 2 R*. The same colour coding as in Fig. 3 has been used.

4.2. Spherically symmetric stellar winds

In this section, we test the performance of the 3D SC method when applied to spherically symmetric, stationary atmospheres. For consistency, the same test models as in Paper I have been calculated, with the wind described by a β-velocity law and the continuity equation:

$$\begin{array}{cc}& v(r)=v_{\infty }(1-b{\displaystyle \frac{R_{\ast }}{r}}{)}^{\beta }\hfill \\ \hfill & b=1-({\displaystyle \frac{v_{\min }}{v_{\infty }}}{)}^{1/\beta }\hfill \\ \hfill & \rho (r)={\displaystyle \frac{\dot{M}}{4\pi r^{2}v(r)}}·\hfill \end{array}$$

For stellar and wind parameters, R_* = 19 R_⊙, Ṁ = 5 × 10⁻⁶ M_⊙ yr⁻¹, β = 1, v_min = 10 km s⁻¹, v_∞ = 2000 km s⁻¹, the density stratification and the velocity field are completely determined. For the considered scattering problems (ϵ_C = ϵ_L = 10⁻⁶), effects of the temperature stratification are negligible. The continuum and (frequency integrated) line opacities have been calculated from Eqs. (7) and (8), with the electron density derived for a completely ionized H/He plasma with helium abundance N_He/N_H = 0.1. We have calculated three different continuum models by scaling the opacity with k_C = [1, 10, 100], respectively. These models correspond to an optically thin, marginally optically thick, and optically thick atmosphere, with radial optical depths τ_r = [0.17, 1.7, 17]. The line transport has been calculated for a weak, intermediate, and strong line, with line-strengths k_L = [1, 10³, 10⁵]. To minimize the computation time, we use a microturbulent velocity v_micro = 100 km s⁻¹ throughout this paper. Such a large velocity dispersion mimicks the effects of multiply non-monotonic velocity fields resulting from the line-driven instability (Paper I and references therein). The atomic mass has been set to m_A = 12 m_p. Finally, the emergent intensity from the stellar core is calculated from T_rad = T_eff = 40 kK.

4.2.1. Convergence behaviour

Figure 5 shows the maximum relative corrections of the mean intensity (left panel) and scattering integral (right panel) after each iteration step. Different methods (SClin, SCbez, and FVM), and different acceleration techniques (classical Λ-iteration, and diagonal-, direct-neighbour-, nearest-neighbour-ALO with the Ng-extrapolation switched on or off) have been applied. We display the continuum and line calculations for k_C = [10, 100] and k_L = [10, 10⁵], respectively, using N_x = N_y = N_z = 93 spatial and N_Ω = 974 angular grid points (to save computation time when calculating the slowly converging classical Λ-iteration). In the following, we only discuss the convergence behaviour of the SClin and SCbez methods, as the FVM has already been analysed in Paper I. We usually stop the iteration scheme when the maximum relative corrections become less than 10⁻³ between subsequent iteration steps, emphasizing that a truly converged solution is only found when the curve describing subsequent relative corrections is sufficiently steep⁶. For instance, the classical Λ-iteration (with Λ^* = 0) fails to converge for strong scattering lines (Fig. 5, lower right panel), since the relative corrections become almost constant in each iteration step (“false convergence”, cf. Hubeny & Mihalas 2014).

Fig. 5. Convergence behaviour of the standard spherically symmetric model calculated with the 3D FVM (blue) and the 3D SC method using linear (red) or Bézier (green) interpolations. Left and right panels: convergence behaviour of the continuum and line transfer for ϵC = ϵL = 10−6, respectively. While the upper row displays the optically thin models with kC = 10 and kL = 10, the lower row has been calculated using kC = 100 and kL = 105. Different acceleration techniques have been applied, where the NN-ALO is implemented only for the SC method.

Overall, and as expected, the number of iterations needed to obtain the converged solution is decreasing with increasing number of matrix elements used to define the ALO (see also Hauschildt et al. 1994; Hauschildt & Baron 2006 for multi-band ALOs coupled to a 1D-SC and 3D-LC formal solution scheme, respectively). In most cases, the convergence of the SClin is faster than that of the SCbez method, because the interpolation scheme is intrinsically more localized (with stronger weights assigned to local Λ-matrix elements). The FVM always performs best, since only the direct neighbours directly influence the formal solution within this method.

For parameters k_C, k_L ≤ 10 (Fig. 5, upper panels), all ALOs yield a converged solution within N_iter ≈ 10 iteration steps. When applying the SCbez method, the first peak results from switching the linear interpolations to Bézier interpolations at the fifth iteration step. This peak is less pronounced for parameters k_C, k_L = 100, 10⁵, since the maximum relative corrections are still relatively large in the first few iteration steps.

For the optically thick continuum model (Fig. 5, lower left panel), the number of iterations until convergence is reduced from $N_{\mathrm{iter}}^{\mathrm{diag}}\approx 75$ to $N_{\mathrm{iter}}^{\mathrm{DN}}\approx 65$ and $N_{\mathrm{iter}}^{\mathrm{NN}}\approx 45$ when using the diagonal, DN-, and NN-ALO within the SClin and SCbez method, respectively. The Ng-extrapolation significantly reduces N_iter further, and is required to obtain the converged solution in ≲20 iteration steps.

For the strong line, the convergence behaviour becomes slightly improved, because the radiative transfer is localized to (narrow) resonance regions. Thus, already the diagonal and DN-ALO yield a solution within $N_{\mathrm{iter}}^{\mathrm{diag}}\approx 50$ and $N_{\mathrm{iter}}^{\mathrm{DN}}\approx 35$ iteration steps. Again, the NN-ALO with the Ng-extrapolation scheme switched on performs best, and reduces the number of iteration steps until convergence to ≲20.

In total, we conclude that a NN-ALO together with the Ng-extrapolation is required for the SClin and SCbez methods, in order to obtain a fast convergence behaviour of the iteration scheme. This ALO also performs excellently for extreme test-cases, that is, for optically thick continua and strong lines in scattering dominated atmospheres.

4.2.2. Continuum and line solution

In the following, we discuss the errors resulting from the upwind and downwind interpolations, and from the integration of the discretized radiative transfer equation. We apply the NN-ALO together with the Ng-extrapolation to ensure convergence. When calculating the line, we used N_x = N_y = N_z = 93 spatial grid points. Since the continuum transport has been calculated at only one frequency point, we applied a higher grid resolution with N_x = N_y = N_z = 133 for such problems. Figure 6 shows the continuum and line solutions together with corresponding relative errors obtained for the spherically symmetric model when calculated with the FVM, SClin, and SCbez methods, and compared to the “exact” 1D solution. The mean and maximum relative errors are shown for different regions in Table 2, where the mean and maximum relative errors of any quantity are defined throughout this paper by

$$\begin{array}{cc}& \overline{\mathrm{\Delta }q}:={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{|q_i-q_i^{(\mathrm{exact})}|}{q_i^{(\mathrm{exact})}}}\hfill \\ \hfill & \mathrm{\Delta }{q}_{\max }:=\underset{\forall i\in [1,N]}{max}{\displaystyle \frac{|q_i-q_i^{(\mathrm{exact})}|}{q_i^{(\mathrm{exact})}}},\hfill \end{array}$$

Fig. 6. Solutions for the standard spherically symmetric models as calculated with the 3D FVM (blue) and 3D SC methods using linear (red) or Bézier (green) interpolations, compared to an accurate 1D solution (black solid line). The dots represent the solutions at specific grid points (with different latitudes and longitudes), where only a subset of all grid points is displayed to compress the image. Corresponding errors are indicated at the bottom of each chart. Top panel: mean intensity for the continuum transfer as a function of radius, with ϵC = 10−6, and kC = [1, 10, 100] from left to right. Bottom panel: line source function with ϵL = 10−6, and kL = [100, 103, 105] from left to right.

with N the number of grid points within the considered region.

For the continuum models, the solutions obtained from the 3D SC methods are superior to the solution obtained from the FVM in most cases. Particularly near the stellar surface (at r ≲ 3 R_*), both SC methods are in good agreement with the 1D solution (see Fig. 6, top panel, and bottom of each chart for the radial dependence of the relative errors). When considering the most challenging problem of optically thick, scattering dominated atmospheres, the mean relative errors of the SClin and SCbez method for the complete calculation region are on the 20- and 10%-level, respectively. For such models, the FVM breaks down due to the order of accuracy (see Paper I), and a (high order) SC method is indeed required to solve the radiative transfer with reasonable accuracy. For marginally optically thick continua, the mean relative errors of the SClin and SCbez methods are on the 5%-level and below. While the FVM allows for a qualitative interpretation of the radiation field for such models, the SC methods should be used for quantitative discussions. The optically thin model calculations give mean relative errors of the order of 5% for all methods, with the maximum relative error being lowest for the FVM. Since, additionally, the computation times of the SClin and SCbez methods are typically highest (see Sect. 3.7), the FVM is to be preferred when calculating optically thin continua. We note that all errors originate from the interplay between upwind and downwind interpolations of opacities, source functions, and intensities, and from the integration of the discretized radiative transfer equation. Numerical diffusion and the associated violation of energy conservation influences the converged solution particularly in the optically thin regime.

The mean relative errors for the line transition are of the order of 5 − 10%, with increasing accuracy from strong to weak lines, and slight advantages of the SCbez method when compared to the SClin method. The radial stratification of relative errors for each considered line is shown in the bottom panel of Fig. 6, bottom of each chart. While the FVM gives largest errors near the stellar surface (at r ≲ 3 R_*), both SC techniques are in excellent agreement with the 1D solution in such regions. At larger radii, however, the SC solutions are generally overestimated when compared to the 1D solution due to numerical diffusion errors. The distinct behaviour of the applied solution schemes in different atmospheric regions finally determines the quality of emergent flux profiles.

4.2.3. Emergent flux profile

The converged source functions are used to calculate the emergent flux profile using the same postprocessing LC solver as in Paper I. Based on Lamers et al. (1987), Busche & Hillier (2005), and Sundqvist et al. (2012), this method solves the equation of radiative transfer in cylindrical coordinates with the z-axis being aligned with the line of sight of the observer’s direction under consideration. All quantities required on the rays are found by trilinear interpolation from the 3D grid, and the equation of radiative transfer is integrated using linear interpolations. To extract the error resulting from the FVM and SC methods alone, we interpolated the “exact” 1D solution onto our 3D grid and calculated the reference profile using the same technique. The continuum has been calculated from a zero-opacity model given by the unattenuated illumination from the projected stellar disc. Then, the differences of line profiles are exclusively related to the differences of line source functions. Figure 7 shows the line profiles with corresponding absolute and relative errors for the intermediate (k_L = 10³) and strong (k_L = 10⁵) line, obtained from the converged source functions from above.

Fig. 7. Emergent flux profiles of an intermediate (kL = 103, top panel) and strong (kL = 105, bottom panel) line. The blue, red, and green curves correspond to the solution of the FVM, SClin, and SCbez methods, respectively. The reference profile (black solid line) has been derived from the “exact” 1D source function interpolated onto the 3D Cartesian grid. Corresponding relative and absolute errors are shown at the bottom of each chart. For all profiles, the continuum level has been determined from a zero-opacity model.

The line profiles are in good agreement with the 1D solution for both applied SC methods, with slight advantages of the SCbez method when compared to the SClin. Major (relative) deviations arise particularly at large frequency shifts on the blue side due to the enlarged source functions in corresponding resonance regions (i.e. at large radii in front of the star). At such frequency shifts, however, the line profile is mainly controlled by absorption, and the absolute error remains small. At low frequency shifts, the emission peak becomes slightly overestimated, particularly when considering the strong line. The corresponding resonance regions are mainly located near to the star (at low absolute velocities) and in the whole plane perpendicular to the line of sight (with low projected velocities). For the intermediate line, the emission from this plane at large radii only plays a minor role due to relatively small optical-depth increments along the line of sight. Thus, both SC methods are in excellent agreement with the 1D reference profile. With increasing line strength, however, the emission from the outer wind region contributes significantly to the line formation, and the discrepancies between the 1D and the SClin/SCbez methods become more pronounced. For all test calculations, the Bézier method performs best, closely followed by the SClin method, and (far behind) the FVM.

With Fig. 7 and the argumentation from above, we conclude that (at least) a short-characteristics solution scheme is required to enable a quantitative interpretation of ultra-violet (UV) resonance line profiles, where both the linear and Bézier interpolation techniques perform similarly well. The less accurate (however computationally cheaper) FVM can still be applied for qualitative discussions.

5. Rotating winds

As a first application of the 3D SC method to non-spherical atmospheres, we consider the UV resonance line formation in the winds of (fast) rotating O stars. Fast rotation has two immediate consequences on the stellar geometry and wind structure. Firstly, the surface of any rotating star becomes distorted, with R_eq/R_pole approaching 3/2 for rotational speeds near the critical velocity (Collins 1963 assuming a Roche model, and the critical velocity defined by Eq. (36) for Ω = 1). Secondly, the emergent flux depends on the (local) effective gravity (corrected for the centrifugal acceleration), and thus, decreases towards the equator (“gravity darkening”, see von Zeipel 1924; Maeder 1999; Maeder & Meynet 2000 for uniform and shellular rotation, respectively). The first attempt to model the winds of fast rotating OB stars was made by Bjorkman & Cassinelli (1993). These authors considered a purely radial line force, and neglected gravity darkening and the surface distortion. Within these approximations, a “wind compressed disc” is formed in the equatorial plane. Cranmer & Owocki (1995) and Owocki et al. (1996) included the effects of non-radial line-forces into their 2D radiation-hydrodynamic simulations, and showed that the formation of the disc becomes suppressed due to a small, but significant polewards acceleration, giving rise to an associated polar velocity component that prevents the formation of a disc. When also accounting for gravity darkening (i.e. a decreased radial acceleration in equatorial regions), Owocki et al. (1996) further showed that a prolate wind structure develops, with decreased equatorial mass loss and velocity (see also the review by Owocki et al. 1998). Maeder (1999) proposed that an oblate wind structure might still be possible, when accounting for a polar variation of the ionization equilibrium induced by gravity darkening (the so-called κ-effect). This effect becomes particularly important when the local effective temperature drops below the bi-stability jump temperature⁷. Petrenz & Puls (2000) extended the hydrodynamic calculations from above by allowing for spatially varying line force multipliers, and showed that no major differences from the prolate wind structure arise, at least for OB stars above T_eff ≳ 20 kK with an optically thin Lyman continuum. Recently, Gagnier et al. (2019) reinvestigated the effects of rotation using 2D ESTER models (see Rieutord et al. 2016 for a description of this code). Using a different implementation of gravity darkening (consistent with the so-called ω-model by Espinosa Lara & Rieutord 2011, which basically results in a slower decrease of effective temperature with co-latitude than obtained from the von Zeipel theorem), these authors predict either a “single-wind regime” (with enhanced polar mass loss) or a “two-wind regime” (with enhanced mass loss at latitudes where the effective temperature drops below the bi-stability jump temperature). To understand which of the different models represents reality best (in different temperature regimes), one needs to compare synthetic profiles with observations. In this respect, investigating the effects of prolate and oblate wind structures is particularly important to distinguish between different theories.

As a consequence of the distinct wind structure resulting from a particular model, the wind lines of rotating stars are expected to differ as a function of rotational speed and inclination. To predict UV resonance line profiles of fast rotating stars, we calculated the source function of a prototypical resonance line transition including the effects of gravity darkening and surface distortion for models with different rotational velocities. As a first step, we used a wind description that is consistent with the prolate wind model. For all calculations, we applied the SClin method.

5.1. Wind model

To obtain a model for the structure of rotating winds, we applied a 2D version of the VH-1 code⁸ developed by J. M. Blondin and co-workers. Our model includes the effects of gravity darkening and surface distortion (see below). Using a 1D input model derived from radiation driven wind theory including finite cone angle corrections (Castor et al. 1975; Pauldrach et al. 1986) for the first time step, the radiation hydrodynamic equations (accounting for non-radial line forces) are solved until a (quasi) stationary solution is obtained (see Cranmer & Owocki 1995; Owocki et al. 1996 for the description of the line force). Assuming azimuthal symmetry, the resulting 2D density and velocity structure is then used as input for our 3D SC code. Table 3 summarizes specific parameters used and obtained for our model calculations. While the surface integrated mass flux, Ṁ, becomes only slightly increased with increasing rotational speed, the polar (equatorial) terminal velocities are significantly enhanced (reduced). For the fastest rotating model (v_rot = 432 km s⁻¹), Fig. 8 shows corresponding density contours in the xz-plane. The z-axis is aligned with the rotation axis. To explicitly show the prolate wind structure, we have scaled the density by the density resulting from the non-rotating (spherically symmetric) model, as a function of distance from the stellar surface. For different rotational speeds, Fig. 9 displays the density and radial velocity along the polar axis and along an (arbitrarily defined) axis in the equatorial plane. When compared with the spherically symmetric wind, the densities of the rotating models are enhanced in polar regions, and become reduced along the equator. Further, the radial velocity along the polar axis remains nearly the same, except in regions far from the star, where the terminal velocity of all rotating models becomes enhanced with increasing v_rot. We note that one would expect clearer differences of the (radial) velocity fields for different rotation rates, due to different accelerations induced, particularly, by the different radiative fluxes resulting from gravity darkening, and, though to a lesser extent, by the specific density structure and rotational velocities. Such differences can be barely observed within our simulations, most presumably because the wind structure in polar regions has not completely settled to a stationary state at the last time steps. In contrast, the radial velocity in equatorial regions is significantly reduced at all distances, when compared to the non-rotating wind, and the deviations from spherical symmetry become more pronounced with increasing rotational velocity. Although we have averaged the hydrodynamic simulations over the last 20 time steps, the atmospheric structure still suffers from small numerical artefacts.

Fig. 8. Contours of the density in the xz-plane (z being the rotation axis) for a prototypical rotating O star with L* = 106 L⊙, M* = 52.5 M⊙, Rp = 19 R⊙, and vrot = 432 km s−1 (corresponding to Ω = 0.9). The density has been scaled by values from the non-rotating model with Ω = 0. While the thick solid line corresponds to the surface of the (distorted) star, the dashed line would correspond to a spherical surface. Additionally, the velocity vectors are displayed.

Fig. 9. Density and radial velocity as a function of distance from the stellar surface in polar (first and third panel) and equatorial (second and fourth panel) regions, for different rotation parameters Ω.

To calculate the stellar surface distortion, we consider the gravitational potential of the star accounting for the effects of centrifugal forces. Under reasonable assumptions, we can approximate this potential by a Roche model (e.g. Collins 1963, see also Cranmer & Owocki 1995):

$$\begin{array}{c}\hfill \mathrm{\Phi }(r,\mathrm{\Theta })=-{\displaystyle \frac{G{M}_{\ast }}{r}}-{\displaystyle \frac{{\omega }^2{r}^2{sin}^2(\mathrm{\Theta })}{2}},\end{array}$$(34)

with angular velocity ω. The surface of the star is defined on equipotential lines and can be calculated by setting Φ(R_p, Θ = 0) = Φ(R_*(Θ),Θ), with R_p the polar radius. Solving the resulting cubic equation, one finds:

$$\begin{array}{c}\hfill R_{\ast }(\mathrm{\Omega },\mathrm{\Theta })={\displaystyle \frac{3R_{\mathrm{p}}}{\mathrm{\Omega }sin(\mathrm{\Theta })}}cos\left({\displaystyle \frac{\pi +{cos}^{-1}(\mathrm{\Omega }sin(\mathrm{\Theta }))}{3}}\right),\end{array}$$(35)

with Ω = ω/ω_crit the ratio of the actual to the critical (“breakup”) angular velocity. Defining the rotational speed at the equator v_rot as input parameter, one easily obtains (cf. Cranmer & Owocki 1995, their Eq. (27)):

$$\begin{array}{cc}& \mathrm{\Omega }={\displaystyle \frac{v_{\mathrm{rot}}}{R_{\mathrm{eq}}}}{\displaystyle \frac{1}{{\omega }_{\mathrm{crit}}}}\hfill \end{array}$$(36)

$$\begin{array}{cc}& R_{\mathrm{eq}}={\displaystyle \frac{R_{\mathrm{p}}}{1-v_{\mathrm{rot}}^2{R}_{\mathrm{p}}/2G{M}_{\ast }}},\hfill \end{array}$$(37)

with equatorial radius R_eq. Following Maeder & Meynet (2000), we use the actual stellar mass to calculate the equatorial radius and critical velocity without correcting for Thomson-accelerations. Additionally, we note that our stellar models are well below the Eddington limit (Γ = 0.5). Thus, the critical angular velocity is simply given by ${\omega }_{\mathrm{crit}}={(8G{M}_{*}/27R_{\mathrm{p}}^3)}^{1/2}$. Instead of using Eq. (35) in our final implementation, we approximated the stellar surface by a spheroid with semi-major axes a = b = R_eq and semi-minor axis c = R_p. Such a formulation greatly simplifies the calculation of the intersection of a given ray with the stellar surface (required for the boundary conditions, see Sect. 3.4). For the most extreme test case considered here (Ω = 0.9), the maximum error on R_*(Θ) due to this approximation is well below the 2%-level, and rapidly decreases for lower rotational velocities.

To calculate the intensity emerging from the stellar core, we set $I_{\mathrm{c}}^{+}(\mathrm{\Theta })=B_{\nu }(T_{\text{eff}}(\mathrm{\Theta }))$, with the effective temperature as a function of co-latitude. For a given luminosity of the star, L_*, we obtain (see also Petrenz & Puls 1996):

$$\begin{array}{cc}& T_{\mathrm{eff}}(\mathrm{\Theta })=[{\displaystyle \frac{L_{\ast }}{2\pi {\sigma }_{\mathrm{B}}\mathrm{\Sigma }}}{|\mathit{g}|}^{4{\beta }_Z}{]}^{1/4}\hfill \\ \hfill & \mathrm{\Sigma }={\displaystyle {\int }_0^{\pi }{|\mathit{g}|}^{4{\beta }_Z}{\displaystyle \frac{R_{\ast }^2(\mathrm{\Theta })sin(\mathrm{\Theta })}{-g_r/|\mathit{g}|}}\mathrm{d}\mathrm{\Theta },}\hfill \end{array}$$(38)

with σ_B the Stefan Boltzmann constant, and the surface integrated effective gravity Σ derived from g(Θ) =  − ∇Φ(R_*(Θ),Θ). The parameter β_Z describes the gravity darkening law in terms of T_eff(Θ)∝|g(Θ)|^β_Z. As originally formulated by von Zeipel (1924), β_Z = 1/4. Though β_Z might be significantly lower (e.g. Domiciano de Souza et al. 2014; Gagnier et al. 2019), for simplicity we nevertheless used β_Z = 1/4. As long as we assume constant ionization fractions, the effect of β_Z on the line profiles will be minor anyway.

5.2. Line formation

For our test models, we used v_micro = 100 km s⁻¹, and calculated the frequency integrated opacity from Eq. (8) for an intermediate and a strong line with line-strength parameter k_L = 10³ and k_L = 10⁵. To obtain the source function, we applied the 3D SClin method and set ϵ_L = 10⁻⁶. The resulting (normalized) line profiles are shown in Fig. 10 for different rotational velocities and inclination angles. Additionally, we display the continuum flux used for normalization in Fig. 11. Due to gravity darkening and the surface distortion, the continuum depends on the rotation rate and inclination, with largest fluxes found for high rotation rates and low inclinations (resulting from the higher temperatures in polar regions and a larger projected stellar disc). In Figs. 10 and 11, the x-axes have been normalized to an (arbitrarily chosen) terminal velocity v_∞ = 3000 km s⁻¹. The behaviour of the line profiles can be qualitatively explained with the hydrodynamic structure:

Fig. 10. Predicted emergent flux profiles for the rotating star models with Ω = [0, 0.5, 0.7, 0.9] (see Table 3). Upper and lower panels: intermediate and strong line with kL = 103 and kL = 105, respectively. The inclination angle has been set to sin(i) = [0, 0.707, 1] from left to right.

Fig. 11. Continuum fluxes for different inclinations sin(i) = [0, 0.707, 1] from left to right, and different rotation parameters (using the same colour coding as in Fig. 10). The continuum fluxes have been scaled by the corresponding flux obtained from the non-rotating model.

• (i)

  For pole on observers (sin(i) = 0, left panel of Fig. 10), the absorption column in front of the star is enhanced with increasing rotational velocity due to the larger densities (and opacities) in polar regions. Thus, the absorption trough (of unsaturated lines) becomes more pronounced. The absorption edge of the intermediate lines is found at slightly lower velocities than expected from the hydrodynamic simulations, because the optical depths of the corresponding resonance regions are too low to efficiently contribute to the absorption. When considering the strong lines, the optical depth is larger, and the absorption edge is more consistent with the actual terminal velocity. For both applied line strength parameters, the emission peak becomes weaker with increasing rotation rate, particularly at low projected velocities on the red side of the line centre (for negative x_obs). This effect can be partly explained by the reduced emission from the equatorial plane, due to the lower densities in these regions. More importantly, however, is the increased continuum flux that mainly determines the (normalized) height of the emission peak.

• (ii)

  When increasing the inclination towards equator-on observers (sin(i) = 1, right panel of Fig. 10), the behaviour is reversed. For such directions, the continuum plays an only minor role, since the lower (equatorial) effective temperatures of the rotating models are nearly compensated by the enlarged projected stellar disc. With increasing rotation parameter, the absorption trough of the intermediate line becomes reduced and shifted towards lower terminal velocities, consistent with the hydrodynamical model. When considering the strong line, the absorption becomes saturated, and only the impact of the different terminal velocities can be observed. Additionally, and for both line strengths, the absorption slightly extends towards the red side, because of (negative) projected line of sight velocities near the stellar surface induced by rotation. For the fastest rotating model with Ω = 0.9, this effect becomes suppressed due to an increased emission from the (dense) prolate wind structure. This latter effect is only moderate for lower rotational speeds. Based on the current hydrodynamic wind structure, we would therefore expect to observe either rather low terminal velocities or relatively deep absorption troughs for fast rotating stars, and we are able, at least in principle, to check the theory by comparing our synthetic spectra with (past or future) UV observations. This point, however, is beyond the scope of this paper.

Finally, if the projected rotational velocity is known (e.g. from photospheric lines), one might even estimate the actual rotational velocity from UV resonance lines. This latter point becomes clear from Fig. 12, where we predict the line profiles of models with different rotational speed for a given v sin (i) (set here to 200 km s⁻¹). Since, at least for the intermediate line, the profile shapes differ, sin (i) could be derived if v sin (i) was known. Of course, such constraints will become feasible only if the underlying models correctly describe the wind structure (including possibly varying ionization stages) of rotating stars.

Fig. 12. Predicted emergent flux profiles for the rotating star models with Ω = [0.5, 0.6, 0.7, 0.8, 0.9], and different inclination angles such that v sin (i) = 200 km s−1 for all models. The upper and lower panels display the intermediate (kL = 103) and strong (kL = 105) line, respectively.

6. Summary and conclusions

In this study, we have presented a 3D short-characteristics method tailored for the solution of continuum- and line-scattering problems in the winds of hot stars. To obtain the formal solution, we have implemented a purely linear interpolation scheme (for calculating upwind quantities and for the solution of the radiative transfer equation along a ray), as well as a second order, monotonic, Bézier technique. We use Cartesian coordinates with a non-uniform grid spacing to ensure a reasonable spatial resolution in important regions (i.e. where velocity and/or density gradients are large). As a first step towards full NLTE radiative transfer models, we consider a single resonance-line transition (approximated by a two-level-atom) assuming an optically thin background continuum, whereas for pure continuum problems we use the thermalization parameter, ϵ_C, and split the source function into a scattering and a thermal part. A generalization (including multi-level atoms) is planned for future applications.

To calculate strong scattering lines and optically thick, scattering dominated continua, we have implemented an accelerated Λ-iteration scheme using different non-local approximate Λ-operators (ALOs), together with applying the Ng-extrapolation method for subsequent iterations. With increasing complexity of the ALO (i.e. from a purely diagonal ALO to a nearest-neighbour ALO including also the 26 neighbouring terms), the rate of convergence is improved. When applying the NN-ALO, the converged solution is generally found within ≲20 iteration steps even for the most challenging test cases.

We have estimated the error of the applied methods in different regimes by calculating spherically symmetric test models within our 3D SC framework, and with a 3D finite-volume method. To our knowledge, this is the first study, where different 3D solution schemes for spherical problems have been compared, and their precision explored. When rated against the solution obtained from (accurate) 1D solvers, we found a mean relative error for the converged continuum source function of roughly 5 − 10% and 5 − 20% when using Bézier and linear interpolations, respectively. Particularly for optically thick continua, the (first order) FVM method breaks down, and a (high order) SC or LC method is required to accurately solve the radiative transfer. When considering the solution of the line source function for different line-strength parameters, the mean relative errors of both SC methods are on the 10%-level and below, with slight advantages of the Bézier technique compared to purely linear interpolations. The resulting synthetic line profiles are calculated with a long-characteristics postprocessing routine using the previously calculated converged source functions. The SC method using Bézier interpolations almost perfectly matches the 1D reference profiles for all our models (i.e. for weak and strong lines). When linear interpolations are used, we obtain slight deviations originating mainly in the outer wind regions. In contrast, the 3D FVM always overestimates the emission. Nevertheless, all methods have their own advantages and disadvantages, particularly when also accounting for the computation time (with fastest turn-around times for the FVM method). Thus, the 3D FVM method should be used for qualitative interpretations and for finding (rough) estimates of the parameters of interest, while the SC methods are to be preferred when aiming at a quantitative analysis of line profiles, and for optically thick continua.

As a first application of the 3D SC code to non-spherical problems, we considered the resonance line formation in the winds of (fast) rotating O stars. Assuming azimuthal symmetry, the hydrodynamic structure for a prototypical O star with different rotation rates has been calculated by means of the 2D VH-1 code. We have included the effects of surface distortion and gravity darkening into our 3D radiative transfer framework. Given the hydrodynamic models, we are able to predict the shape of line profiles for different rotational speeds and inclination angles. When compared with a (non-rotating) spherically symmetric wind (obtained using the same stellar parameters), rotating stars should either show relatively low terminal velocities (for equator-on observers) or deeper absorption troughs (for pole-on observers). The latter effect, however, would only be observable when considering intermediate (i.e. unsaturated) lines. Additionally, we showed that the line profile shapes vary as a function of rotational speed at a given v sin (i). Thus, assuming that v sin (i) was known (e.g. from photospheric line modelling), one could estimate the rotational speed, though with a rather large uncertainty, since the differences of the line profiles are only moderate. We emphasize that other effects (such as clumping or a flatter gravity darkening law) may additionally shape the line profiles. When analysing UV observations of fast rotating stars, the 3D SC code developed in this work certainly will help to understand the manifestations of various (aforementioned) effects, and to distinguish between different theoretical predictions (e.g. prolate vs. oblate wind structures). Ideal testbeds for future investigations of fast rotating winds are the stars VFTS102 (O9 Vnnne, Dufton et al. 2011) and VFTS285 (O7.5 Vnnn, Walborn et al. 2012), both rotating at nearly their critical velocity.

Finally, we note that our tools are, of course, not limited to rotating stars. Indeed, almost any kind of stellar wind that deviates from spherical symmetry (with non-relativistic velocity fields), such as magnetic winds or colliding winds in close binaries, can be investigated.

Acknowledgments

We thank our referee, Dr. Achim Feldmeier, for many helpful comments and suggestions, that improved the first version of this manuscript considerably. LH gratefully acknowledges support from the German Research Foundation, DFG, under grant PU 117/9-1. NDK acknowledges funding from the KU Leuven C1 grant MAESTRO C16/17/007. JOS acknowledges support from FWO Odysseus grant D4545-GH9218. Finally, we thank J. M. Blondin and co-workers for providing the VH-1 code.

References

Appendix A: Grid construction

In this section, we describe the grid construction procedure used within our 3D code. Generally, we assume the wind structure (i.e. density, velocity field, and temperature) to be given by an input model obtained from hydrodynamic simulations or external (semi)-analytic calculations. Since the input grid is not necessarily compatible with our 3D SC solver, and to minimize interpolation errors when calculating upwind and downwind quantities, we construct an own grid that uses the distribution of the input-grid coordinates in an optimum way. When the input grid uses spherical coordinates (r, Θ, Φ), we define a joint probability distribution

$$\begin{array}{c}\hfill h(x,z)=f(r)g(\mathrm{\Theta })|J|,\end{array}$$(A.1)

where f(r) and g(Θ) are the probability density functions derived from the distribution of the input coordinates, and $|J|=\sqrt{x^2+z^2}$ is the Jacobian determinant. Since we consider only axisymmetric atmospheres in this paper, we use the x-coordinate distribution also for the y-coordinates. To calculate the probability density functions for x and z, we simply marginalize h(x, z) over z and x, respectively. The discretized coordinates are finally determined by demanding that the probabilities of selecting a (continuous) coordinate in each (discrete) interval shall be the same. Figure A.1 shows the probability density functions of the x-coordinates for two different input distributions of the radial grid. Here, the polar angle Θ has been assumed to be uniformly distributed for both examples.

Fig. A.1. Probability density functions of radial (dashed) and x (solid) coordinates for different spherical and Cartesian grids. In this example, two spherical grids are given in 2D as input to our 3D code, with uniformly (black) or logarithmically (red) distributed r-coordinates, and a uniformly distributed polar angle. The corresponding distributions of x-coordinates are calculated within our grid construction procedure (see text). Large values of the probability density functions correspond to a high resolution of x and r-coordinates.

Since the final number of core and non-core points depends on the slope of the probability density of the radial grid, yielding in worst cases a much larger number of core points than non-core points, and because the total number of used points is memory-limited, we define two input parameters N_core and N_{non − core} to keep control on the final grid. For all test calculations (including those presented in Sect. 4.2), the best solution has always been found for a number of N_core/N_{non − core} ∈ [0.25, 0.5]. An explicit choice of N_core and N_{non − core} corresponds to a re-normalization of the probability density function in the regimes x, z ∈ [0, R_*] and x, z ∈ [R_*, R_max], where R_max defines the border of the computational domain. We note that the same procedure can be used for an input grid given in Cartesian coordinates, with the probability density function of the input-grid coordinates derived directly from the corresponding (discrete) distribution.

Appendix B: 1D Bézier interpolation

In this section, we discuss an interpolation technique using quadratic Bézier curves (e.g. Auer 2003; Schwarz 1997). Such curves are generally constructed from three given points b₀, b₁, b₂ (see Fig. B.1), and have the following useful properties:

Fig. B.1. Bézier curves (solid lines) for three given points b0, b1, b2. The blue, red, and green lines represent the resulting curves for different control points b1. The straight connections of the control points with the data points are indicated by the dotted lines.

• (i)

  The boundary points, b₀ and b₂ are reproduced exactly by the Bézier curve.

• (ii)

  The straight connections (b₁−b₀) and (b₂−b₁) define the tangent lines of the curve at b₀ and b₂, respectively.

• (iii)

  Any point on the Bézier curve is located in the convex hull of b₀, b₁, b₂.

In a 2D plane described by coordinates x and y, the quadratic Bézier curve is parameterized as:

$$\begin{array}{c}\hfill \mathit{b}(t)=\left(\begin{array}{c}x(t)\\ y(t)\end{array}\right)={(1-t)}^2{\mathit{b}}_0+2t(1-t){\mathit{b}}_1+t^2{\mathit{b}}_2,\end{array}$$(B.1)

with t ∈ [0, 1], and b₀ = (x₀, y₀), b₁ = (x₁, y₁), b₂ = (x₂, y₂). With Eq. (B.1), the properties (i)–(iii) can be exploited to construct a monotonic interpolation scheme by identifying b₀,b₂ with two given data points (x₀, f₀), (x₂, f₂), and defining b₁ as a free (and tunable) parameter. Thus, b₁ is commonly named control point, and is “only” required to set the slope of the Bézier curve. To reproduce the underlying function best, and to preserve monotonicity of the resulting curve, the control point should be chosen with care.

In the following, we present a Bézier-interpolation technique for an interval x ∈ [x_i − 1, x_i], given three data points, (x_i − 1, f_i − 1), (x_i, f_i), (x_i + 1, f_i + 1). The interpolation formulas corresponding to the interval x ∈ [x_i, x_i + 1] are given in Appendix B.2.

B.1. Interval [x_i − 1, x_i]

A quadratic Bézier curve in the interval [x_i − 1, x_i] is given from Eq. (B.1):

$$\begin{array}{c}\hfill \left(\begin{array}{c}x(t)\\ f\left(x\left(t\right)\right)\end{array}\right)={(1-t)}^2\left(\begin{array}{c}{x}_{i-1}\\ f_{i-1}\end{array}\right)+2t(1-t)\left(\begin{array}{c}{x}_{\mathrm{c}}\\ f_{\mathrm{c}}\end{array}\right)+t^2\left(\begin{array}{c}{x}_i\\ f_i\end{array}\right),\end{array}$$(B.2)