Issue 
A&A
Volume 633, January 2020



Article Number  A16  
Number of page(s)  24  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201936584  
Published online  23 December 2019 
A 3D shortcharacteristics method for continuum and line scattering problems in the winds of hot stars
^{1}
LMU München, Universitätssternwarte, Scheinerstr. 1, 81679 München, Germany
email: levin@usm.unimuenchen.de, levin@usm.lmu.de
^{2}
Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
Received:
27
August
2019
Accepted:
8
October
2019
Context. Knowledge about hot, massive stars is usually inferred from quantitative spectroscopy. To analyse nonspherical phenomena, the existing 1D codes must be extended to higher dimensions, and corresponding tools need to be developed.
Aims. We present a 3D radiative transfer code that is capable of calculating continuum and line scattering problems in the winds of hot stars. By considering spherically symmetric test models, we discuss potential error sources, and indicate advantages and disadvantages by comparing different solution methods. Further, we analyse the ultraviolet (UV) resonance line formation in the winds of rapidly rotating O stars.
Methods. We consider both a (simplified) continuum model including scattering and thermal sources, and a UV resonance line transition approximated by a twolevelatom. We applied the shortcharacteristics (SC) method, using linear or monotonic Bézier interpolations, for which monotonicity is of prime importance, to solve the equation of radiative transfer on a nonuniform Cartesian grid. To calculate scattering dominated problems, our solution method is supplemented by an accelerated Λiteration scheme using newly developed nonlocal operators.
Results. For the spherical test models, the mean relative error of the source function is on the 5 − 20% level, depending on the applied interpolation technique and the complexity of the considered model. All calculated line profiles are in excellent agreement with corresponding 1D solutions. Close to the stellar surface, the SC methods generally perform better than a 3D finitevolumemethod; however, they display specific problems in searchlightbeam tests at larger distances from the star. The predicted line profiles from fast rotating stars show a distinct behaviour as a function of rotational speed and inclination. This behaviour is tightly coupled to the wind structure and the description of gravity darkening and stellar surface distortion.
Conclusions. Our SC methods are ready to be used for quantitative analyses of UV resonance line profiles. When calculating optically thick continua, both SC methods give reliable results, in contrast to the alternative finitevolume method.
Key words: radiative transfer / methods: numerical / stars: massive / stars: rotation / stars: winds / outflows
© 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” stellarmass 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 GravitationalWave Observatory (aLIGO), Abbott et al. 2016), the formation of heavy stellarmass 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 singlestar evolution, however, requires comparatively moderate mass loss rates (e.g. in low metalicity environments, or massloss 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; WMBASIC: 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. udDoula & Owocki 2002; udDoula 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 DavidUraz et al. (2019), for instance. Additionally, (nonradial) 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 linedriving, 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), multiD 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 socalled formal solution), and iterating the updated sources and sinks until convergence (Cannon 1973).
In order to obtain the formal solution in multiD, the following two major methods exist, each having specific advantages and disadvantages. Firstly, within the finitevolume method (FVM, see Adam 1990 for the first application in the context of radiative transfer problems), the calculation volume is discretized into finite subvolumes. 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 “longcharacteristics” method (LC, Jones & Skumanich 1973; Jones 1973), and the “shortcharacteristics” 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 finitevolume method from Paper I, when appropriate. As a first application to nonspherical 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 timeindependent equation of radiative transfer (e.g. Mihalas 1978), formulated in the observer’s frame:
I_{ν} describes the specific intensity for a given frequency ν and direction n, and χ_{ν} and S_{ν} are the opacities (in units of cm^{−1}) and source functions for continuum and line processes. For our following tests, we either consider a pure continuum case with source function
and thermalization parameter ϵ_{C}, or a line (within an optically thin continuum) approximated by a twolevelatom (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:
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:
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, . V is the local velocity vector in units of , and
is the ratio between the local thermal velocity (accounting for microturbulent velocities) and the fiducial thermal velocity. This description of the profile function allows for a depthindependent frequency grid (see also Paper I).
Finally, we express the continuum and line opacities in terms of the Thomsonopacity, χ_{Th} = n_{e}σ_{e}, and depthindependent input parameters, k_{C}, k_{L}:
where the linestrength k_{L} is related to the frequency integrated opacity (e.g. Paper I, Eqs. (10) and (12), and with = (cm s)^{−1}). 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 timeconsuming 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 longcharacteristics 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 nonuniform, 3dimensional 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 nonmonotonic velocity fields, we use the observer’s frame formulation.
SC methods have been successfully implemented already for 3D nonLTE (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
where dτ := χds is the opticaldepth 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 and , while local quantities are denoted either as or simply q_{ijk}. For a given ray segment, we then obtain:
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 yzplanes, 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 coordinateaxes and are defined in Sect. 3.3. 

Open with DEXTER 
with upwind opticaldepth 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 opticaldepth increments, Δτ_{u} and Δτ_{d} (see below), are required. To calculate the source contribution, the source function is commonly approximated by first or secondorder 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 nonlinear, because its elements now explicitly depend on the stratification of source functions (via corresponding interpolation/integration coefficients). Within any Λiteration scheme, this nonlinearity can lead to oscillations. In extreme cases, “flipflop 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^{1}, 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 datapoints 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^{0}, t^{1}, t^{2} terms:
with downwind opticaldepth 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 , , and 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 flipflop 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 alwaysconvergent 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:
with
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:
where Δs_{u}, Δs_{d} describe the path lengths of the upwind and downwind intervals, respectively, and , 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 socalled 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 . 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 secondorder (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 sspace (Eqs. (B.7)–(B.9)):
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. 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 intergrid points, such that the (local) upwind, current, and downwind quantities are now described by the indices [ℓ − 1, ℓ, ℓ + 1], we obtain:
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:
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)}, , 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
where θ is the colatitude (measured from the Cartesian zaxis), and ϕ is the azimuth (measured from the xaxis), the distances from a considered grid point to the neighbouring xy, xz, and yzplanes are calculated from trigonometry and yield:
with α, β, γ set to ±1 for directionvector 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 , 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 (index2) to (index) to determine upwind quantities, while downwind quantities are calculated from (index1) 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 xyplane at the upwind side. The 2D Bézier interpolation for the upwind point then consists of three 1D Bézier interpolations along the xaxis 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 yaxis at the upwind xcoordinates. With the 2D Bézier interpolation given by Eq. (C.1), we find for each required quantity q_{u, d}:
where the coefficients w^{(ijk)} and 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 (nonzero) 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 nonlinear Λ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 xzplane). Firstly, the considered ray originates from the stellar surface (direction n_{1} in Fig. 2). In this case, we use a corehalo approximation and set , with 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 socalled “ghost points”. Secondly, the considered ray originates from a plane spanned by grid points that are partially located inside the star (direction n_{2} in Fig. 2). Then, the interpolation is performed as in Sect. 3.3, using again representative estimates at the corepoints.
Fig. 2. Boundary conditions for rays propagating in the xzplane at ylevel (j) with three different directions n_{1}, n_{2}, n_{3}, and upwind points u_{1}, u_{2}, u_{3}. For point u_{1}, the intensity is set to 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 u_{2} 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 n_{3}, the unknown intensity inside the core, , is directed inwards. Such situations occur only for ray directions (nearly) parallel to the spatial grid, and thus are relatively seldom. 

Open with DEXTER 
Inside the core, we define and set and all velocity components to zero, where is the inward directed intensity which needs to be specified only in rare situations (Fig. 2, direction n_{3}). 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^{2}. 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:
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 stepfunction 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 GaussLegendre integration performs (slightly) better. However, the corresponding directions are always clustered in certain regions since the nodes of the GaussLegendre quadrature are fixed. Additionally, the GaussLegendre 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 GaussLegendre integration should not be used (see Table 1).
Mean relative error (defined in Sect. 4.2) of the mean intensity for a zeroopacity model, obtained using the Lebedev, GaussLegendre, and trapezoidal integration method, the latter with nodes from Lobel & Blomme (2008).
We have tested a multitude of other quadrature schemes, including trapezoidal and (pseudo)Gaussian rules on triangles, and the socalled 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 noncore 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 Lebedevintegration 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 finitevolume method.
To obtain the scattering integral, we apply the trapezoidal rule for the frequency integration. The scattering integral then reads:
with and 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 nonlinear 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:
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:
with the nth unit vector e_{n}, and matrix indices m, n related to the 3D indices (i, j, k) by
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 matrixelement describes the effect of a nonvanishing source function at grid point n onto grid point m. We emphasize that Eq. (29) holds only for precalculated 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 operatorsplitting methods (Cannon 1973). Within the ALI, the Λoperator is written as
where the first term is an appropriately chosen ALO acting on the new source function, , and the second term acts on the previous one, . 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:
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 multilevel atoms, we emphasize that Λ and Λ^{(A)} may change within the ALIcycle 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:
Equation (33) is solved to obtain a new source function (see below). Since, however, has been optimized only to ensure monotonicity in a specific step k − 1 (based on source function ), the iteration scheme can oscillate due to oscillations in and . Even worse, the new source function might become negative. To overcome these problems, nonlinear 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 ) 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, approaches . 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 matrixelements to be calculated, and the resulting convergence speed. While the inversion of a diagonal ALO reduces to a simple division, a multiband 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 JacobiIteration for sparse systems (see Paper I). For 3D calculations, a multiband 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 26neighbour 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 multilevel atoms (for which the occupation numbers and thus opacities might change during the iteration)^{3}.
In this paper, we analyse the convergence speed of the ALI for a diagonal ALO given by Eq. (D.14), a “directneighbour” (DN)ALO given by Eqs. (D.5), (D.11), (D.13)–(D.15), (D.17), (D.23), and a “nearestneighbour” (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 DNALO, since the diagonal and directneighbour elements depend on several other neighbours through the inclusion of downwind interpolations. Since, however, the downwindintegration 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 thirdorder 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 NNALO, in parallel to the formal solution requires 20%, 30%, and 40% of the total computation time.
Finally, we have implemented the Ngextrapolation (Ng 1974; Olson et al. 1986) to improve the convergence speed further. As in Paper I, the Ngacceleration 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 nonuniform gridspacing 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 (nearestneighbour) Λ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 h and 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, s, and 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^{4}. 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. Searchlightbeam test
A first test of our 3D SC methods is the searchlightbeam 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 , 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 xzplane. 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. Searchlightbeam test for direction n = (1, 0, 1) and a typical grid with N_{x} = N_{y} = N_{z} = 133 grid points. Upper panel: contour plot of the specific intensity as calculated with the SC method using Bézier interpolations in the xzplane (cf. Paper I, Fig. 3, for the finitevolume 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_{*}. 

Open with DEXTER 
Along the beam centre, the SC methods perfectly reproduce the exact solution, whereas the FVM solution decreases significantly due to the finite gridcell 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 to zero for rays intersecting the stellar surface or not, due to the corehalo approximation. As a consequence, almost all interpolations (and interpolation schemes) overestimate the specific intensity^{5}. For optically thick models (where 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 gridcell 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 upwindinterpolation 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. 

Open with DEXTER 
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:
For stellar and wind parameters, R_{*} = 19 R_{⊙}, Ṁ = 5 × 10^{−6} M_{⊙} yr^{−1}, β = 1, v_{min} = 10 km s^{−1}, v_{∞} = 2000 km s^{−1}, the density stratification and the velocity field are completely determined. For the considered scattering problems (ϵ_{C} = ϵ_{L} = 10^{−6}), 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 linestrengths k_{L} = [1, 10^{3}, 10^{5}]. To minimize the computation time, we use a microturbulent velocity v_{micro} = 100 km s^{−1} throughout this paper. Such a large velocity dispersion mimicks the effects of multiply nonmonotonic velocity fields resulting from the linedriven 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, directneighbour, nearestneighbourALO with the Ngextrapolation switched on or off) have been applied. We display the continuum and line calculations for k_{C} = [10, 100] and k_{L} = [10, 10^{5}], 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^{−3} between subsequent iteration steps, emphasizing that a truly converged solution is only found when the curve describing subsequent relative corrections is sufficiently steep^{6}. 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 k_{C} = 10 and k_{L} = 10, the lower row has been calculated using k_{C} = 100 and k_{L} = 10^{5}. Different acceleration techniques have been applied, where the NNALO is implemented only for the SC method. 

Open with DEXTER 
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 multiband ALOs coupled to a 1DSC and 3DLC 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^{5}, 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 to and when using the diagonal, DN, and NNALO within the SClin and SCbez method, respectively. The Ngextrapolation 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 DNALO yield a solution within and iteration steps. Again, the NNALO with the Ngextrapolation scheme switched on performs best, and reduces the number of iteration steps until convergence to ≲20.
In total, we conclude that a NNALO together with the Ngextrapolation 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 testcases, 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 NNALO together with the Ngextrapolation 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
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 k_{C} = [1, 10, 100] from left to right. Bottom panel: line source function with ϵ_{L} = 10^{−6}, and k_{L} = [10^{0}, 10^{3}, 10^{5}] from left to right. 

Open with DEXTER 
Mean and maximum relative errors of the FVM and SClin, SCbez methods, when applied to spherically symmetric test models.
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 zaxis 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 zeroopacity 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^{3}) and strong (k_{L} = 10^{5}) line, obtained from the converged source functions from above.
Fig. 7. Emergent flux profiles of an intermediate (k_{L} = 10^{3}, top panel) and strong (k_{L} = 10^{5}, 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 zeroopacity model. 

Open with DEXTER 
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 opticaldepth 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 shortcharacteristics solution scheme is required to enable a quantitative interpretation of ultraviolet (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 nonspherical 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 nonradial lineforces into their 2D radiationhydrodynamic 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 socalled κeffect). This effect becomes particularly important when the local effective temperature drops below the bistability jump temperature^{7}. 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 socalled ωmodel by Espinosa Lara & Rieutord 2011, which basically results in a slower decrease of effective temperature with colatitude than obtained from the von Zeipel theorem), these authors predict either a “singlewind regime” (with enhanced polar mass loss) or a “twowind regime” (with enhanced mass loss at latitudes where the effective temperature drops below the bistability 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 VH1 code^{8} developed by J. M. Blondin and coworkers. 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 nonradial 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^{−1}), Fig. 8 shows corresponding density contours in the xzplane. The zaxis is aligned with the rotation axis. To explicitly show the prolate wind structure, we have scaled the density by the density resulting from the nonrotating (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 nonrotating 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 xzplane (z being the rotation axis) for a prototypical rotating O star with L_{*} = 10^{6} L_{⊙}, M_{*} = 52.5 M_{⊙}, R_{p} = 19 R_{⊙}, and v_{rot} = 432 km s^{−1} (corresponding to Ω = 0.9). The density has been scaled by values from the nonrotating 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. 

Open with DEXTER 
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 Ω. 

Open with DEXTER 
Specific parameters used and obtained for the rotating wind models.
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):
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:
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)):
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 Thomsonaccelerations. Additionally, we note that our stellar models are well below the Eddington limit (Γ = 0.5). Thus, the critical angular velocity is simply given by . Instead of using Eq. (35) in our final implementation, we approximated the stellar surface by a spheroid with semimajor axes a = b = R_{eq} and semiminor 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 , with the effective temperature as a function of colatitude. For a given luminosity of the star, L_{*}, we obtain (see also Petrenz & Puls 1996):
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^{−1}, and calculated the frequency integrated opacity from Eq. (8) for an intermediate and a strong line with linestrength parameter k_{L} = 10^{3} and k_{L} = 10^{5}. To obtain the source function, we applied the 3D SClin method and set ϵ_{L} = 10^{−6}. 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 xaxes have been normalized to an (arbitrarily chosen) terminal velocity v_{∞} = 3000 km s^{−1}. 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 k_{L} = 10^{3} and k_{L} = 10^{5}, respectively. The inclination angle has been set to sin(i) = [0, 0.707, 1] from left to right. 

Open with DEXTER 
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 nonrotating model. 

Open with DEXTER 

(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 equatoron 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^{−1}). 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 (k_{L} = 10^{3}) and strong (k_{L} = 10^{5}) line, respectively. 

Open with DEXTER 
6. Summary and conclusions
In this study, we have presented a 3D shortcharacteristics method tailored for the solution of continuum and linescattering 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 nonuniform 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 resonanceline transition (approximated by a twolevelatom) 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 multilevel 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 nonlocal approximate Λoperators (ALOs), together with applying the Ngextrapolation method for subsequent iterations. With increasing complexity of the ALO (i.e. from a purely diagonal ALO to a nearestneighbour ALO including also the 26 neighbouring terms), the rate of convergence is improved. When applying the NNALO, 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 finitevolume 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 linestrength 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 longcharacteristics 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 turnaround 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 nonspherical 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 VH1 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 (nonrotating) spherically symmetric wind (obtained using the same stellar parameters), rotating stars should either show relatively low terminal velocities (for equatoron observers) or deeper absorption troughs (for poleon 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 nonrelativistic velocity fields), such as magnetic winds or colliding winds in close binaries, can be investigated.
In this paper, different interpolation schemes are tested by considering simplified (though physically relevant) continuum and line scattering problems. We emphasize that our code will be further developed to enable the solution of more complex, multilevel problems in 3D. For such problems, highly accurate interpolation schemes are required to describe the variation of the mean intensities along a ray.
For a given grid point, the number of neighbouring grid points that can be used for extrapolation is not a priori clear, and depends on the shape of the stellar surface, and the considered direction of the ray. Indeed, there are 64 special cases that would have to be implemented explicitly. This is computationally not feasible.
The 1D solution for the continuum transport is found from the Rybickialgorithm (combined with the solution of the moment equations using variable Eddington factors, see, e.g. Mihalas 1978). To calculate the line, a comovingframe raybyray solution scheme in pzgeometry is applied, ensuring convergence with a diagonal ALO.
In contrast, the number of photons entering and leaving a given grid cell is (nearly) conserved within the FVM by definition. This statement, however, is not completely true for the FVM as formulated by Adam (1990), Lobel & Blomme (2008), and Hennicker et al. (2018), since all these authors apply an (averaged) upwind approximation.
For linearly convergent iteration schemes, the steepness of the convergence curve is described by the relative corrections in subsequent iterations steps, Δ_{k}/Δ_{k − 1} = const. = :q. To obtain a solution within a reasonable amount of computation time, we may demand that q ≲ 0.8, corresponding to a reduction of relative errors by a factor of 10^{−3} every 30th iteration step.
The jump temperature is theoretically motivated by a stronger radiative linedriving due to lower ionization stages of iron for T_{eff} ≲ T_{jump} ≈ 25 kK (Vink et al. 1999). More recently, Petrov et al. (2016) predicted a somewhat lower jump temperature, T_{jump} ≈ 20 kK.
S_{p} ≠ 0 only when considering the grid point p ↔ (ijk). Then, , and only when the upwind point is located on the stellar surface (Eq. (D.14)).
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/91. NDK acknowledges funding from the KU Leuven C1 grant MAESTRO C16/17/007. JOS acknowledges support from FWO Odysseus grant D4545GH9218. Finally, we thank J. M. Blondin and coworkers for providing the VH1 code.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Adam, J. 1990, A&A, 240, 541 [NASA ADS] [Google Scholar]
 Ahrens, C., & Beylkin, G. 2009, Proc. R. Soc. A: Math. Phys. Eng. Sci., 465 [NASA ADS] [CrossRef] [Google Scholar]
 Auer, L. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 3 [NASA ADS] [Google Scholar]
 Beentjes, C. H. L. 2015, Quadrature on a Spherical Surface (Oxford: University of Oxford) [Google Scholar]
 Bjorkman, J. E., & Cassinelli, J. P. 1993, ApJ, 409, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Busche, J. R., & Hillier, D. J. 2005, AJ, 129, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Cannon, C. J. 1973, ApJ, 185, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Cherepashchuk, A. M. 1976, Sov. Astron. Lett., 2, 138 [NASA ADS] [Google Scholar]
 Collins, II., G. W. 1963, ApJ, 138, 1134 [NASA ADS] [CrossRef] [Google Scholar]
 Cranmer, S. R., & Owocki, S. P. 1995, ApJ, 440, 308 [NASA ADS] [CrossRef] [Google Scholar]
 DavidUraz, A., Erba, C., Petit, V., et al. 2019, MNRAS, 483, 2814 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dufton, P. L., Dunstall, P. R., Evans, C. J., et al. 2011, ApJ, 743, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
 Espinosa Lara, F., & Rieutord, M. 2011, A&A, 533, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, G., & Espinosa Lara, F. 2019, A&A, 625, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Georgiev, L. N., Hillier, D. J., & Zsargó, J. 2006, A&A, 458, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Koesterke, L., & Hamann, W.R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H. 1992, J. Quant. Spec. Radiat. Transf., 47, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hauschildt, P. H., Störzer, H., & Baron, E. 1994, J. Quant. Spec. Radiat. Transf., 51, 875 [NASA ADS] [CrossRef] [Google Scholar]
 Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288 [NASA ADS] [CrossRef] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Holzreuter, R., & Solanki, S. K. 2012, A&A, 547, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton: Princeton University Press) [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jones, H. P. 1973, ApJ, 185, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, H. P., & Skumanich, A. 1973, ApJ, 185, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Keszthelyi, Z., Wade, G. A., & Petit, V. 2017, in The Lives and DeathThroes of Massive Stars, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, IAU Symp., 329, 250 [NASA ADS] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P. B., & Olson, G. L. 1988, J. Quant. Spec. Radiat. Transf., 39, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., CerrutiSola, M., & Perinotto, M. 1987, ApJ, 314, 726 [NASA ADS] [CrossRef] [Google Scholar]
 Leenaarts, J., & Carlsson, M. 2009, in The Second Hinode Science Meeting: Beyond DiscoveryToward Understanding, eds. B. Lites, M. Cheung, T. Magara, J. Mariska, & K. Reeves, ASP Conf. Ser., 415, 87 [NASA ADS] [Google Scholar]
 Lobel, A., & Blomme, R. 2008, ApJ, 678, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 1999, A&A, 347, 185 [NASA ADS] [Google Scholar]
 Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
 Marcolino, W. L. F., Bouret, J.C., Sundqvist, J. O., et al. 2013, MNRAS, 431, 2253 [NASA ADS] [CrossRef] [Google Scholar]
 Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009, AJ, 137, 3358 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W. H. Freeman and Co.) [Google Scholar]
 Ng, K.C. 1974, J. Chem. Phys., 61, 2680 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spec. Radiat. Transf., 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1996, ApJ, 472, L115 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., Gayley, K. G., & Cranmer, S. R. 1998, in Properties of Hot Luminous Stars, ed. I. Howarth, ASP Conf. Ser., 131, 237 [NASA ADS] [Google Scholar]
 Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petit, V., Owocki, S. P., Wade, G. A., et al. 2013, MNRAS, 429, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, V., Keszthelyi, Z., MacInnis, R., et al. 2017, MNRAS, 466, 1052 [NASA ADS] [CrossRef] [Google Scholar]
 Petrenz, P., & Puls, J. 1996, A&A, 312, 195 [NASA ADS] [Google Scholar]
 Petrenz, P., & Puls, J. 2000, A&A, 358, 956 [NASA ADS] [Google Scholar]
 Petrov, B., Vink, J. S., & Gräfener, G. 2016, MNRAS, 458, 1999 [NASA ADS] [CrossRef] [Google Scholar]
 Prilutskii, O. F., & Usov, V. V. 1976, Sov. Astron., 20, 2 [NASA ADS] [Google Scholar]
 Puls, J. 1991, A&A, 248, 581 [NASA ADS] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Rivero González, J. G., Puls, J., Najarro, F., & Brott, I. 2012, A&A, 537, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rybicki, G. B., & Hummer, D. G. 1978, ApJ, 219, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schwarz, H. R. 1997, Numerische Mathematik (Stuttgart: B. G. Teubner) [CrossRef] [Google Scholar]
 Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., udDoula, A., Owocki, S. P., et al. 2012, MNRAS, 423, L21 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., & Owocki, S. P. 2002, ApJ, 576, 413 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., Owocki, S. P., & Townsend, R. H. D. 2008, MNRAS, 385, 97 [NASA ADS] [CrossRef] [Google Scholar]
 van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066 [NASA ADS] [CrossRef] [Google Scholar]
 Vanbeveren, D. 1991, A&A, 252, 159 [NASA ADS] [Google Scholar]
 Vath, H. M. 1994, A&A, 284, 319 [NASA ADS] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181 [NASA ADS] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Walborn, N. R., Sana, H., Taylor, W. D., SimónDíaz, S., & Evans, C. J. 2012, in Proceedings of a Scientific Meeting in Honor of Anthony F. J. Moffat, eds. L. Drissen, C. Robert, N. StLouis, & A. F. J. Moffat, ASP Conf. Ser., 465, 490 [NASA ADS] [Google Scholar]
 Zsargó, J., Hillier, D. J., & Georgiev, L. N. 2006, A&A, 447, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
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 inputgrid coordinates in an optimum way. When the input grid uses spherical coordinates (r, Θ, Φ), we define a joint probability distribution
where f(r) and g(Θ) are the probability density functions derived from the distribution of the input coordinates, and is the Jacobian determinant. Since we consider only axisymmetric atmospheres in this paper, we use the xcoordinate distribution also for the ycoordinates. 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 xcoordinates 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 rcoordinates, and a uniformly distributed polar angle. The corresponding distributions of xcoordinates are calculated within our grid construction procedure (see text). Large values of the probability density functions correspond to a high resolution of x and rcoordinates. 

Open with DEXTER 
Since the final number of core and noncore 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 noncore points, and because the total number of used points is memorylimited, 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 renormalization 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 inputgrid 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_{0}, b_{1}, b_{2} (see Fig. B.1), and have the following useful properties:
Fig. B.1. Bézier curves (solid lines) for three given points b_{0}, b_{1}, b_{2}. The blue, red, and green lines represent the resulting curves for different control points b_{1}. The straight connections of the control points with the data points are indicated by the dotted lines. 

Open with DEXTER 

(i)
The boundary points, b_{0} and b_{2} are reproduced exactly by the Bézier curve.

(ii)
The straight connections (b_{1}−b_{0}) and (b_{2}−b_{1}) define the tangent lines of the curve at b_{0} and b_{2}, respectively.

(iii)
Any point on the Bézier curve is located in the convex hull of b_{0}, b_{1}, b_{2}.
In a 2D plane described by coordinates x and y, the quadratic Bézier curve is parameterized as:
with t ∈ [0, 1], and b_{0} = (x_{0}, y_{0}), b_{1} = (x_{1}, y_{1}), b_{2} = (x_{2}, y_{2}). With Eq. (B.1), the properties (i)–(iii) can be exploited to construct a monotonic interpolation scheme by identifying b_{0},b_{2} with two given data points (x_{0}, f_{0}), (x_{2}, f_{2}), and defining b_{1} as a free (and tunable) parameter. Thus, b_{1} 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ézierinterpolation 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):
with (x_{c}, f_{c}) the control point. The abscissa of the control point, x_{c}, can be chosen arbitrarily (at least in principal). To obtain a secondorder interpolation scheme, however, x_{c} needs to be located at the centre of the datapoint’s abscissae^{9}, and is therefore set to x_{c} = (x_{i − 1} + x_{i})/2. Then, the quadratic Bézier interpolation scheme is given by:
where t has been determined from the definition of x_{c} and Eq. (B.2). Since the straight connection of the control point (x_{c}, f_{c}) with the data point (x_{i}, f_{i}) defines the tangent line of the Bézier curve at this data point, f_{c} is calculated as
with Δx_{i} = x_{i} − x_{i − 1}. The unknown derivative at x_{i} needs to be approximated. Using also the information from the next data point, (x_{i + 1}, f_{i + 1}), and assigning a weight ω to the forward and backward derivatives (obtained from finite differences), we find
with Δx_{i + 1} = x_{i + 1} − x_{i}. With a proper choice of ω, we can adjust the Bézier curve to our needs by shifting the control point up or down. For instance, setting ω = Δx_{i + 1}/(Δx_{i} + Δx_{i + 1}) results in the unique parabola connecting the three given data points, while ω = 1 would yield the linear interpolation. To avoid overshoots and negative function values, we demand that the Bézier curve shall be monotonic in the interval [x_{i − 1}, x_{i}]. Noting that monotonicity is obtained when the control point is located in the interval f_{c} ∈ [f_{i − 1}, f_{i}], corresponding ωvalues should lie in between the following limits:
where the superscript [i − 1, i] denotes that ω corresponds to the interpolation scheme in the left interval, [x_{i − 1}, x_{i}]. In the final implementation, we avoid the division by zero if f_{i} = f_{i − 1} or f_{i} = f_{i + 1}, of course. Our standard interpolation is then performed as follows. At first, we calculate ω such that we obtain the unique parabola connecting all three data points. Secondly, if ω lies outside the allowed limits from Eq. (B.5) and (B.6), we adjust ω to yield monotonic interpolations. In Fig. B.2, we display the monotonic Bézier curve resulting from a ωparameter calculated by means of Eq. (B.5), together with linear and quadratic interpolations (the latter connecting the three data points). Since monotonicity is always obtained for ω ∈ [ω_{i − 1}, ω_{i}], we can define even stricter limits in order to avoid oscillations during the iteration scheme, by setting ω = 1 to obtain purely linear interpolations, for instance (see Sect. 3.6).
Fig. B.2. Different interpolation techniques for a set of three data points at xcoordinates indicated by the dotted vertical lines. The solid and dashed lines correspond to the interpolation in the different intervals [x_{i − 1}, x_{i}] and [x_{i}, x_{i + 1}], respectively. Linear interpolations, quadratic interpolations (connecting all three data points), and a monotonic Bézier curve (with ω calculated from Eq. (B.5) in the interval [x_{i − 1}, x_{i}]) are indicated in red, blue, and green. Since the quadratic interpolation is already monotonic in the interval [x_{i}, x_{i + 1}], the monotonic Bézier curve coincides with the dashed, blue line. Control points are indicated with coloured asterisks. 

Open with DEXTER 
To calculate the elements of the (approximate) Λmatrix, the interpolation coefficients are required. Combining Eqs. (B.3) and (B.4) then gives:
with
B.2. Interval [x_{i}, x_{i + 1}]
The interpolation formula for the right interval [x_{i}, x_{i + 1}] uses the same data points as above. Since the value of the control point needs to be calculated at a different xcoordinate, x_{c} = (x_{i + 1} + x_{i})/2, we cannot simply substitute indices. Using
we obtain for this interval:
with
and
The corresponding Bézier curves for different ωparameters (ω = 1 for linear and ω = Δx_{i}/(Δx_{i} + Δx_{i + 1}) for continuous quadratic interpolations) are also shown in Fig. B.2. We note that the Bézier interpolation gives a continuous function in the complete interval [x_{i − 1}, x_{i + 1}] only for those ωvalues that define the parabola connecting all three data points.
Appendix C: 2D Bézier interpolation
To interpolate upwind and downwind quantities, a 2D interpolation scheme is required. Figure C.1 displays the geometry for a 2D rectangular area, with grid points indicated by the black dots. With this setup, we perform a 2D Bézier interpolation by applying three 1D Bézier interpolations along the xaxis on each ylevel at (j − 1), (j), (j + 1), followed by another 1D Bézier interpolation along y at the desired xcoordinate. Within the cyan shaded interval, we obtain with the 1D Bézier interpolation given by Eqs. (B.13)–(B.16):
Fig. C.1. 2D interpolation for upwind or downwind quantities required in the cyan shaded area. The 2D Bézier interpolation consists of three 1D interpolations to obtain the values at the desired xcoordinate (indicated by red dots), followed by a 1D interpolation along the ycoordinate using the obtained values at the red dots. 

Open with DEXTER 
where the subscripts of the interpolation coefficients indicate the coordinate used for each 1D interpolation. We note that all upwind and downwind interpolations are performed in the upper right interval of a given surface, in order to obtain a simple representation of the Λmatrix elements.
Appendix D: ALO coefficients
In this section, we derive the Λmatrix coefficients used to construct the approximate Λoperator. We note that the obtained matrix elements can also be used for any other (2nd or lower order) interpolation scheme using the same geometry, with different interpolation coefficients though.
For a source function set to unity at grid point (ijk) and zero everywhere else, we consider all 27 points ranging from (i − α, j − β, k − γ) to (i + α, j + β, k + γ). The corresponding matrix coefficients are derived from Eq. (29), using the discretized equation of radiative transfer, Eqs. (12)/(20), with upwind and downwind interpolations defined by Eqs. (22) and (23). We further consider only the Λ_{Ων}operator, since the integration over frequency and/or solid angle is straightforward. Each Λmatrix element then corresponds to the intensity (resulting from S_{ijk} = 1) at a considered grid point p (not necessarily identical to (ijk)), and consists of an emission term (defined by the interpolated source functions and opticaldepth steps at the corresponding upwind, current, and downwind points), and the irradiation from the upwind point (defined by the upwind intensity and upwind opticaldepth step). The upper and lower panel of Fig. D.1 show an example in 2D considering the points (i − α, j − β) and (i, j − β) for a source function S_{ij} = 1.
Fig. D.1. 2D example for calculating the Λ_{Ω, ν}matrix elements at a grid point (i − α, j − β) (upper panel) and (i, j − β) (lower panel). The matrix elements correspond to the intensity at the considered grid points calculated for a source function S_{ij} = 1 and zero everywhere else. For such a configuration, the downwind source function is interpolated from grid points indicated with the green dots, while the upwind source function and upwind intensity are obtained from the red dots (for simplicity we here assume linear interpolations for determining upwind and downwind quantities). We emphasize that the upwind intensity vanishes only when considering the grid point (i − α, j − β). 

Open with DEXTER 
In the following, we sketch the derivation of the (3D) matrix element for the first neighbour, and only present the solution for the remaining ones. To save space, we skip the indices Ω, ν. The m, nth Λelement is written as , with matrix indices n, m calculated from Eq. (30). While n corresponds to the 3D indices of the local grid point (S_{ijk} = S_{n} = 1), m represents the neighbouring point, (i − α, j − β, k − γ). Applying Eq. (29) to the specific intensity at point m, we obtain:
with boundary contribution Φ_{B}, nth unit vector e_{n}, and , , the Kroneckerδ for all possible , , and coordinates, respectively. S_{u} and S_{d} are the upwind and downwind source functions corresponding to a considered short characteristic at grid point p ↔ (i − α, j − β, k − γ), S_{p} is the source function at the grid point^{10}, I_{u} is the upwind intensity, and a, b, c, d are the integration coefficients for this particular short characteristic. All upwind and downwind quantities are to be interpolated from neighbouring grid points. We use the notation w, , , to identify different interpolation coefficients corresponding to the upwind source function, upwind intensity, and downwind source function, respectively. Using Eqs. (22) and (23) to interpolate upwind and downwind quantities, we find:
with the upwind intensity interpolated from the same points as the upwind source function, and a compact notation for the interpolation coefficients (with skipped superscripts). Since only S_{ijk} = 1 (and zero everywhere else), and because the upwind intensity vanishes (for this particular grid point, see Fig. D.1 for an example in 2D), we finally obtain:
The matrix element for a point (i − α, j − β, k − γ) with a nonvanishing source function at point (ijk) is thus solely given by the integration coefficient c_{i − α, j − β, k − γ} from the discretized equation of radiative transfer multiplied with the interpolation coefficient for the downwind source function of point I (corresponding to grid point (ijk), see Fig. 1). The other neighbours are obtained analogously, without vanishing incident intensities, however. Accounting also for the interpolation of upwind source functions and intensities when necessary, we find:
Since the integration and interpolation coefficients need to be calculated only once at each considered grid point (here denoted by (u, v, w), to avoid confusion), we obtain the Λmatrix coefficients by substituting indices. For Eq. (D.1), we find:
and proceed analogously for all other elements in Eqs. (D.2)–(D.27). Thus, the ALO can be calculated in parallel to the formal solution scheme.
All Tables
Mean relative error (defined in Sect. 4.2) of the mean intensity for a zeroopacity model, obtained using the Lebedev, GaussLegendre, and trapezoidal integration method, the latter with nodes from Lobel & Blomme (2008).
Mean and maximum relative errors of the FVM and SClin, SCbez methods, when applied to spherically symmetric test models.
All Figures
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 yzplanes, 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 coordinateaxes and are defined in Sect. 3.3. 

Open with DEXTER  
In the text 
Fig. 2. Boundary conditions for rays propagating in the xzplane at ylevel (j) with three different directions n_{1}, n_{2}, n_{3}, and upwind points u_{1}, u_{2}, u_{3}. For point u_{1}, the intensity is set to 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 u_{2} 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 n_{3}, the unknown intensity inside the core, , is directed inwards. Such situations occur only for ray directions (nearly) parallel to the spatial grid, and thus are relatively seldom. 

Open with DEXTER  
In the text 
Fig. 3. Searchlightbeam test for direction n = (1, 0, 1) and a typical grid with N_{x} = N_{y} = N_{z} = 133 grid points. Upper panel: contour plot of the specific intensity as calculated with the SC method using Bézier interpolations in the xzplane (cf. Paper I, Fig. 3, for the finitevolume 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_{*}. 

Open with DEXTER  
In the text 
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. 

Open with DEXTER  
In the text 
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 k_{C} = 10 and k_{L} = 10, the lower row has been calculated using k_{C} = 100 and k_{L} = 10^{5}. Different acceleration techniques have been applied, where the NNALO is implemented only for the SC method. 

Open with DEXTER  
In the text 
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 k_{C} = [1, 10, 100] from left to right. Bottom panel: line source function with ϵ_{L} = 10^{−6}, and k_{L} = [10^{0}, 10^{3}, 10^{5}] from left to right. 

Open with DEXTER  
In the text 
Fig. 7. Emergent flux profiles of an intermediate (k_{L} = 10^{3}, top panel) and strong (k_{L} = 10^{5}, 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 zeroopacity model. 

Open with DEXTER  
In the text 
Fig. 8. Contours of the density in the xzplane (z being the rotation axis) for a prototypical rotating O star with L_{*} = 10^{6} L_{⊙}, M_{*} = 52.5 M_{⊙}, R_{p} = 19 R_{⊙}, and v_{rot} = 432 km s^{−1} (corresponding to Ω = 0.9). The density has been scaled by values from the nonrotating 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. 

Open with DEXTER  
In the text 
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 Ω. 

Open with DEXTER  
In the text 
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 k_{L} = 10^{3} and k_{L} = 10^{5}, respectively. The inclination angle has been set to sin(i) = [0, 0.707, 1] from left to right. 

Open with DEXTER  
In the text 
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 nonrotating model. 

Open with DEXTER  
In the text 
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 (k_{L} = 10^{3}) and strong (k_{L} = 10^{5}) line, respectively. 

Open with DEXTER  
In the text 
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 rcoordinates, and a uniformly distributed polar angle. The corresponding distributions of xcoordinates are calculated within our grid construction procedure (see text). Large values of the probability density functions correspond to a high resolution of x and rcoordinates. 

Open with DEXTER  
In the text 
Fig. B.1. Bézier curves (solid lines) for three given points b_{0}, b_{1}, b_{2}. The blue, red, and green lines represent the resulting curves for different control points b_{1}. The straight connections of the control points with the data points are indicated by the dotted lines. 

Open with DEXTER  
In the text 
Fig. B.2. Different interpolation techniques for a set of three data points at xcoordinates indicated by the dotted vertical lines. The solid and dashed lines correspond to the interpolation in the different intervals [x_{i − 1}, x_{i}] and [x_{i}, x_{i + 1}], respectively. Linear interpolations, quadratic interpolations (connecting all three data points), and a monotonic Bézier curve (with ω calculated from Eq. (B.5) in the interval [x_{i − 1}, x_{i}]) are indicated in red, blue, and green. Since the quadratic interpolation is already monotonic in the interval [x_{i}, x_{i + 1}], the monotonic Bézier curve coincides with the dashed, blue line. Control points are indicated with coloured asterisks. 

Open with DEXTER  
In the text 
Fig. C.1. 2D interpolation for upwind or downwind quantities required in the cyan shaded area. The 2D Bézier interpolation consists of three 1D interpolations to obtain the values at the desired xcoordinate (indicated by red dots), followed by a 1D interpolation along the ycoordinate using the obtained values at the red dots. 

Open with DEXTER  
In the text 
Fig. D.1. 2D example for calculating the Λ_{Ω, ν}matrix elements at a grid point (i − α, j − β) (upper panel) and (i, j − β) (lower panel). The matrix elements correspond to the intensity at the considered grid points calculated for a source function S_{ij} = 1 and zero everywhere else. For such a configuration, the downwind source function is interpolated from grid points indicated with the green dots, while the upwind source function and upwind intensity are obtained from the red dots (for simplicity we here assume linear interpolations for determining upwind and downwind quantities). We emphasize that the upwind intensity vanishes only when considering the grid point (i − α, j − β). 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.