A 3D short-characteristics method for continuum and line scattering problems in the winds of hot stars
aa r X i v : . [ a s t r o - ph . S R ] O c t Astronomy&Astrophysicsmanuscript no. paper˙v15 c (cid:13)
ESO 2019October 30, 2019
A 3D short-characteristics method for continuum and linescattering problems in the winds of hot stars
L. Hennicker , J. Puls , N. D. Kee , and J. O. Sundqvist LMU M¨unchen, Universit¨atssternwarte, Scheinerstr. 1, 81679 M¨unchen, Germany, e-mail: [email protected] Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, BelgiumReceived 27 August 2019; Accepted 8 October 2019
ABSTRACT
Context.
Knowledge about hot, massive stars is usually inferred from quantitative spectroscopy. To analyse non-spherical 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 hotstars. By considering spherically symmetric test models, we discuss potential error sources, and indicate advantages and disadvantagesby comparing di ff erent solution methods. Further, we analyse the ultra-violet (UV) resonance line formation in the winds of rapidlyrotating O stars. Methods.
We consider both a (simplified) continuum model including scattering and thermal sources, and a UV resonance linetransition approximated by a two-level-atom. We applied the short-characteristics (SC) method, using linear or monotonic B´ezierinterpolations, for which monotonicity is of prime importance, to solve the equation of radiative transfer on a non-uniform Cartesiangrid. To calculate scattering dominated problems, our solution method is supplemented by an accelerated Λ -iteration scheme usingnewly developed non-local operators. Results.
For the spherical test models, the mean relative error of the source function is on the 5 −
20 % percent level, depending onthe applied interpolation technique and the complexity of the considered model. All calculated line profiles are in excellent agreementwith corresponding 1D solutions. Close to the stellar surface, the SC methods generally perform better than a 3D finite-volume-method; however, they display specific problems in searchlight-beam tests at larger distances from the star. The predicted line profilesfrom fast rotating stars show a distinct behaviour as a function of rotational speed and inclination. This behaviour is tightly coupledto 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 opticallythick continua, both SC methods give reliable results, in contrast to the alternative finite-volume method.
Key words.
Radiative transfer – Methods: numerical – Stars: massive – Stars: rotation – Stars: winds, outflows
1. Introduction
Understanding of hot, massive stars is a basic prerequisite for in-terpreting fundamental properties of our Universe. Already dur-ing their lifetimes, such stars influence the evolution of galax-ies via feedback of ionizing radiation and radiatively drivenwinds, and they enrich the interstellar medium (ISM) with met-als. But also, their deaths are of large impact. Depending on ini-tial mass and mass loss, OB stars either explode as supernovaeor end up as ‘heavy’ stellar-mass black holes (Heger et al. 2003).Obviously, supernova explosions enrich the ISM with metalseven further. Additionally, the associated shock fronts possiblytrigger 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 ad-vanced Laser Interferometric Gravitational-Wave Observatory(aLIGO), Abbott et al. 2016), the formation of heavy stellar-mass black holes becomes of key interest. Since massive starsare frequently found to be members of multiple star systems(see, e.g. Mason et al. 2009, Sana et al. 2013), they might ex-plain the occurrence of GW events in the correct mass range. Atleast the formation of heavy black holes from single-star evolu-tion, however, requires comparatively moderate mass loss rates (e.g. in low metalicity environments, or mass-loss quenching bymagnetic fields, Petit et al. 2017, Keszthelyi et al. 2017).Current knowledge of OB stars is inferred from quanti-tative spectroscopy, that is, by comparing observed spectrawith synthetic ones, the latter being obtained from numeri-cally modelling their stellar atmospheres (photosphere + wind).State of the art atmospheric modelling is performed by assum-ing spherical symmetry (e.g. CMFGEN: Hillier & Miller 1998;PHOENIX: Hauschildt 1992; P o WR: Gr¨afener et al. 2002; WM- basic : Pauldrach et al. 2001; FASTWIND: Puls et al. 2005 andRivero Gonz´alez et al. 2012).Several e ff ects may alter the geometry of hot star atmo-spheres, a ff ecting both the photospheric and wind lines. Forinstance, binary interaction (see, e.g. Vanbeveren 1991 andde Mink et al. 2013 for a discussion about evolutionary aspects,and Prilutskii & Usov 1976, Cherepashchuk 1976, Stevens et al.1992 for the e ff ects of colliding winds) influences the line for-mation in such objects. Another example is the phenomenon ofmagnetic winds (e.g. ud-Doula & Owocki 2002, ud-Doula et al.2008, Petit et al. 2013), for which the resonance line forma-tion has been extensively discussed in Marcolino et al. (2013),Hennicker et al. (2018), and David-Uraz et al. (2019), for in-stance. Additionally, (non-radial) pulsations may impact the ge-ometry of hot star atmospheres.
1. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
In this paper, we focus on (rapidly) rotating stars. In addi-tion to polar velocity components that a ff ect the radiative line-driving, such stars and their winds are influenced by centrifugalforces and gravity darkening. For a detailed discussion, we referto Sect. 5.To analyse these e ff ects, and to distinguish between di ff er-ent theories (resulting in, e.g. prolate vs. oblate wind structures),multi-D radiative transfer codes that include a treatment of (ar-bitrary) velocity fields are required. In hot star winds, a keychallenge is the implementation of scattering processes. Suchproblems are most easily handled by an accelerated Λ -iterationscheme, that is, by calculating the radiation field for knownsources and sink terms (the so-called formal solution), and iter-ating the updated sources and sinks until convergence (Cannon1973).In order to obtain the formal solution in multi-D, the fol-lowing two major methods exist, each having specific ad-vantages and disadvantages. Firstly, within the finite-volumemethod (FVM, see Adam 1990 for the first application in thecontext of radiative transfer problems), the calculation vol-ume is discretized into finite sub-volumes. All required phys-ical values at the cell centres are then considered as suitableaverages within each cell. This way, a relatively simple solu-tion scheme can be derived, which, however, is only of loworder. Thus, the FVM breaks down for large optical depths(Hennicker et al. 2018, hereafter paper I). Nevertheless, it is stilluseful for qualitative interpretations (see, e.g. Lobel & Blomme2008 for an application to corotating interaction regions inthe winds of rotating stars). Secondly, within the character-istics methods, the formal solution is obtained by exact in-tegration of the equation of radiative transfer along a ray.One distinguishes between the ‘long-characteristics’ method(LC, Jones & Skumanich 1973, Jones 1973), and the ’short-characteristics’ method (SC, Kunasz & Auer 1988). The LCmethod follows a particular ray from the boundary to the consid-ered grid point, whereas the SC method considers rays only fromcell to cell. All required quantities at each upwind point are ob-tained from interpolation. While the LC method is computation-ally more expensive than the SC method, the latter introducesnumerical errors from the interpolation scheme. To date, few 3Dcodes that apply one or the other technique already exist, such asP hoenix/
3D (Hauschildt & Baron 2006 and other publications inthis series) and IRIS (Ibgui et al. 2013), which, however, are ei-ther proprietary, or do not include a suitable Λ -iteration scheme.For additional information on other codes, we refer to Sect. 1 ofpaper I.In this paper, we present a newly developed 3D SC code,which is capable of treating arbitrary velocity fields, and in-cludes an accelerated Λ -iteration scheme to handle continuumand line scattering problems. In Sect. 2, we discuss the basicassumptions of the code. The numerical methods are describedin Sect. 3. In Sect. 4, we perform a detailed error analysis bycomparing our 3D results for spherically symmetric atmosphereswith corresponding 1D solutions. Additionally, we present com-parisons with the 3D finite-volume method from paper I, whenappropriate. As a first application to non-spherical winds, wecalculated UV resonance line profiles for di ff erent models of ro-tating stars in Sect. 5, and discuss the implications. Our conclu-sions are summarized in Sect. 6.
2. Basic assumptions
In this paper, we use the same basic assumptions as in paper I,that is we aim to solve the time-independent equation of radiative transfer (e.g. Mihalas 1978), formulated in the observer’s frame: n ∇ I ν ( r , n ) = χ ν ( r , n ) (cid:0) S ν ( r , n ) − I ν ( r , n ) (cid:1) . (1) I ν describes the specific intensity for a given frequency ν anddirection n , and χ ν and S ν are the opacities (in units of cm − )and source functions for continuum and line processes. For ourfollowing tests, we either consider a pure continuum case withsource function S C = (1 − ǫ C ) J ν + ǫ C B ν , (2)and thermalization parameter ǫ C , or a line (within an opticallythin continuum) approximated by a two-level-atom (TLA). Ageneralization to a multitude of lines, coupled via correspond-ing rate equations, is one of our next objectives. Within the TLAapproach, the line source function reads: S L = (1 − ǫ L ) ¯ J + ǫ L B (3) ǫ L = ǫ ′ + ǫ ′ , ǫ ′ = C ul A ul " − exp (cid:16) − h ν k B T (cid:17) . (4)¯ J is the ‘scattering integral’ (see Eq. (25)), and C ul and A ul de-scribe the collisional rate and Einstein coe ffi cient for sponta-neous emission, respectively. In this paper, ǫ L is considered asinput parameter. The profile function Φ x is approximated by aDoppler profile: Φ x = √ πδ exp " − (cid:16) x cmf δ (cid:17) = √ πδ exp " − (cid:16) x obs − n · V δ (cid:17) , (5)where x cmf and x obs denote the frequency shift from line centre inthe comoving and observer’s frame, in units of a fiducial Dopplerwidth, ∆ ν ∗ D = ν v ∗ th / c . V is the local velocity vector in units of v ∗ th , and δ = v ∗ th r k B Tm A + v (6)is the ratio between the local thermal velocity (accounting formicro-turbulent velocities) and the fiducial thermal velocity.This description of the profile function allows for a depth-independent frequency grid (see also paper I).Finally, we express the continuum and line opacities in termsof the Thomson-opacity, χ Th = n e σ e , and depth-independent in-put parameters, k C , k L : χ C = k C · χ Th (7) χ L = k L · χ Th · Φ x = : ¯ χ L Φ x , (8)where the line-strength k L is related to the frequency integratedopacity ¯ χ = ∆ ν ∗ D ¯ χ L = k L ∆ ν ∗ D χ Th (e.g. paper I, Eqs. (10) and(12), and with [ ¯ χ ] = (cm s) − ). In the following, we present(e ffi cient) numerical tools for solving the coupled equations (1)and (2) / (3). Since a direct solution is computationally prohibitivefor 3D atmospheres due to limited memory capacity, we applyan accelerated Λ -iteration (ALI) scheme that overcomes the con-vergence problems of the classical Λ -iteration for optically thick,scattering dominated atmospheres by using an appropriate ap-proximate Λ -operator (ALO).
3. Numerical methods
The most time-consuming part of the ALI scheme consists ofcalculating the formal solution. In paper I, we have shown thata formal solution obtained using the FVM su ff ers from various
2. Hennicker et al.: 3D short-characteristics method in the winds from OB stars numerical inaccuracies related to numerical di ff usion and to theorder of accuracy, the latter influencing the solution particularlyin the optically thick regime. To avoid these errors, we imple-ment an integral method along short characteristics. When com-pared with a long-characteristics solution scheme, the compu-tation time becomes reduced by roughly a factor of N /
2, with N the number of spatial grid points (per dimension), since theLC method integrates the equation of radiative transfer alongthe complete path from the boundary to a considered grid point.Thus, LC methods become feasible only on massively paral-lelized architectures.We follow the same approach as in paper I, and solve theequation of radiative transfer on a non-uniform, 3-dimensionalCartesian grid. In contrast to curvilinear coordinate systems, thedirection vector n then becomes constant with respect to the spa-tial grid, thus avoiding the otherwise required angular interpo-lation of upwind intensities and a complicated bookkeeping ofintensities (and corresponding integration weights) for di ff erentdirections. Furthermore, the sweep through the spatial domainis considerably simplified, and can be performed grid point bygrid point along the x , y , and z coordinates, since the intensities(required for the upwind interpolation) are always known on theprevious grid layer. To enable a straightforward implementationof non-monotonic velocity fields, we use the observer’s frameformulation.SC methods have been successfully implemented alreadyfor 3D non-LTE (NLTE) radiative transfer problems in coolstars (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 ac-count for subsonic and slightly supersonic velocity fields.For scattering problems including highly supersonic velocityfields, there exist, to our knowledge, only the 2D codes byDullemond & Turolla (2000) (planar / spherical), van Noort et al.(2002) (planar / spherical / cylindrical), Georgiev et al. (2006) andZsarg´o et al. (2006) (spherical). The only 3D SC code includingarbitrary velocity fields, IRIS (Ibgui et al. 2013), has also beenformulated for planar geometries, and lacks the implementationof a Λ -iteration scheme thus far. As has been shown in all thesestudies, the final performance of the SC method crucially de-pends on the choice of the applied interpolation schemes. The equation of radiative transfer along a given direction can bewritten as d I d τ = S − I , (9)where d τ : = χ d s is the optical-depth increment along a raysegment d s . Here and in the following, we suppress the nota-tion for the frequency dependence, and explicitly distinguish be-tween continuum and line only when appropriate. Eq. (9) is in-tegrated along a ray propagating through a current grid point pwith Cartesian grid indices ( i jk ) and corresponding upwind pointu ( i jk ) . The geometry for a 3D Cartesian grid is shown in the upperpanel of Fig. 1. In the following, upwind and downwind quanti-ties corresponding to a considered grid point ( i jk ) are indicatedby q ( i jk )u and q ( i jk )d , while local quantities are denoted either as q ( i jk )p or simply q i jk . For a given ray segment, we then obtain: I i jk = I ( i jk )u e − ( τ p − τ u ) + Z τ p τ u e − ( τ p − τ ) S ( τ )d τ = I ( i jk )u e − ∆ τ u + e − ∆ τ u Z ∆ τ u e t S ( t + τ u )d t , (10)with upwind optical-depth increment ∆ τ u : = τ p − τ u ≥ t : = τ − τ u . For the SC solution scheme, the locationof the reference point, τ =
0, plays no role, since only theoptical-depth increments, ∆ τ u and ∆ τ d (see below), are re-quired. To calculate the source contribution, the source functionis commonly approximated by first- or second-order polynomi-als (Kunasz & Auer 1988, van Noort et al. 2002), B´ezier curves(Hayek et al. 2010, Holzreuter & Solanki 2012, Auer 2003) orcubic Hermite splines (Ibgui et al. 2013). While the 2nd- (andhigher) order methods reproduce the di ff usion limit in opticallythick media, they su ff er from overshoots and need to be mono-tonized with some e ff ort to ensure that any interpolated quantityremains positive between two given grid points. Monotonicityis usually obtained by manipulating the interpolation schemewhenever overshoots occur. Thus, the actual interpolation cru-cially depends on the specific stratification of the consideredquantity (e.g. the source function). The Λ -operator then be-comes non-linear, because its elements now explicitly depend onthe stratification of source functions (via corresponding interpo-lation / integration coe ffi cients). Within any Λ -iteration scheme,this non-linearity can lead to oscillations. In extreme cases, ‘flip-flop situations’ (Holzreuter & Solanki 2012, their Appendix A)may occur, which do not converge at all.For the source contribution, we implement both a linear ap-proximation as the fastest and most stable method (monotonicityis always provided), and a quadratic B´ezier approximation (seeAppendix B) for higher accuracy , which allows us to preservemonotonicity in a rather simple way. The B´ezier interpolation isconstructed from two given data points and one control point,the latter setting the slope of the interpolating curve. The controlpoint is located at the centre of the data-points abscissae, withthe ordinate calculated by accounting for the information of athird data point to yield the parabola intersecting all three datapoints. Whenever overshoots occur, the value of the control pointwill be manipulated to ensure monotonicity (see Fig. B.2). Thecorresponding formulation is given in Appendix B, Eqs. (B.7)-(B.10). Applying these equations to describe the behaviour ofthe mean intensities along the optical path, and identifying theindices ( i − i ), ( i +
1) with the upwind, current, and down-wind points, we find, after reordering for the t , t , t terms: S ( t + τ u ) = S ( i jk )u + " ( ω − ∆ τ u S ( i jk )u + (1 − ω ) ∆ τ u + (2 − ω ) ∆ τ d ∆ τ u ∆ τ d S ( i jk )p + ( ω − ∆ τ d S ( i jk )d · t + " (1 − ω ) ∆ τ S ( i jk )u + ( ω − ∆ τ u + ∆ τ d ) ∆ τ ∆ τ d S ( i jk )p + (1 − ω ) ∆ τ u ∆ τ d S ( i jk )d · t , (11) In this paper, di ff erent interpolation schemes are tested by consider-ing simplified (though physically relevant) continuum and line scatter-ing problems. We emphasize that our code will be further developed toenable the solution of more complex, multi-level problems in 3D. Forsuch problems, highly accurate interpolation schemes are required todescribe the variation of the mean intensities along a ray. 3. Hennicker et al.: 3D short-characteristics method in the winds from OB stars with downwind optical-depth increment, ∆ τ d = τ d − τ p ≥
0. Theparameter ω defines the ordinate of the control point (Eq. (B.4)).Within the B´ezier interpolation, we emphasize that ω may ex-plicitly depend on S ( i jk )u , S ( i jk )p , and S ( i jk )d to ensure monotonicity,and not solely on the grid spacing. A major advantage of this pa-rameterization is that we can globally define a minimum allowed ω that can be adapted during the iteration process. The flip-flopsituations discussed above can then be avoided by gradually in-creasing ω min towards unity ( ω ≡ I i jk = a i jk S ( i jk )u + b i jk S ( i jk )p + c i jk S ( i jk )d + d i jk I ( i jk )u , (12)with a i jk : = e + ω − ∆ τ u e + − ω ∆ τ e b i jk : = (1 − ω ) ∆ τ u + (2 − ω ) ∆ τ d ∆ τ u ∆ τ d e + ( ω − ∆ τ u + ∆ τ d ) ∆ τ ∆ τ d e c i jk : = ω − ∆ τ d e + − ω ∆ τ u ∆ τ d e d i jk : = e − ∆ τ u e : = e − ∆ τ u Z ∆ τ u e t d t = − e − ∆ τ u e : = e − ∆ τ u Z ∆ τ u t e t d t = ∆ τ u − e e : = e − ∆ τ u Z ∆ τ u t e t d t = ∆ τ − e . The calculation of the upwind and downwind ∆ τ -steps proceedssimilarly, where now the opacity is integrated using the B´ezierinterpolation. Using Eqs. (B.3), (B.4) for the upwind interval,and Eqs. (B.11), (B.12) for the downwind interval, one easilyobtains: ∆ τ u = Z pu χ ( s )d s = ∆ s u χ u + χ [u , p]c + χ p ) (13) ∆ τ d = Z dp χ ( s )d s = ∆ s d χ p + χ [p , d]c + χ d ) , (14)where ∆ s u , ∆ s d describe the path lengths of the upwind anddownwind intervals, respectively, and χ [u , p]c , χ [p , d]c refer to theopacity at the control points in each interval. Since the opacity of a line transition depends on the velocityfield through the Doppler e ff ect, regions of significant opacitymay become spatially confined in a highly supersonic wind withstrong acceleration. Thus, a grid refinement along the short char-acteristic might be required to correctly account for all so-calledresonance zones. Because the profile function is approximatedby a Doppler profile and rapidly vanishes for | x cmf /δ | &
3, a res-onance zone is here defined by a region where x cmf /δ ∈ [ − , ffi cient condition to resolve all such reso-nance zones along a given ray is to demand that | ∆ x cmf | /δ = | ∆ V proj | /δ . / ( i jk ) , p], where ∆ V proj is the projected velocity step along theray in units of v ∗ th . Assuming a linear dependence of the projectedvelocities on the ray coordinate s , this condition directly trans-lates to an equidistant refined spatial grid along the ray. For shortray segments (as is mostly the case within our calculations), ne-glecting the second-order (curvature) terms of the projected ve-locity influences the solution only weakly. The line source func-tion on the refined grid is obtained by B´ezier interpolation in s -space (Eqs. (B.7)-(B.9)): S L ( s ℓ ) = S ℓ = ˜ a ℓ S ( i jk )u + ˜ b ℓ S ( i jk )p + ˜ c ℓ S ( i jk )d , (15)where the index ℓ refers to the points on the refined grid, andu ( i jk ) , p ( i jk ) , d ( i jk ) describe the original geometry of the short char-acteristic. ¯ χ L is obtained analogously, and the required ∆ τ ℓ steps are calculated with the trapezoidal rule, for simplicity.Contrasted to the Sobolev method (which also assumes a lin-ear velocity law along the ray segment, e.g. Rybicki & Hummer1978), our grid refinement procedure explicitly accounts forvariations of the opacity and the source function.Using Eq. (12) for the inter-grid points, such that the (local)upwind, current, and downwind quantities are now described bythe indices [ ℓ − , ℓ, ℓ + I ℓ = I ℓ − e − ∆ τ ℓ + (cid:16) a ℓ ˜ a ℓ − + b ℓ ˜ a ℓ + c ℓ ˜ a ℓ + (cid:17) S ( i jk )u (16) + (cid:16) a ℓ ˜ b ℓ − + b ℓ ˜ b ℓ + c ℓ ˜ b ℓ + (cid:17) S ( i jk )p (17) + (cid:16) a ℓ ˜ c ℓ − + b ℓ ˜ c ℓ + c ℓ ˜ c ℓ + (cid:17) S ( i jk )d (18): = I ℓ − e − ∆ τ ℓ + ˜ α ℓ S ( i jk )u + ˜ β ℓ S ( i jk )p + ˜ γ ℓ S ( i jk )d . (19)For a number of N ref refinement points (including the upwindand current point) within the interval [u ( i jk ) , p ( i jk ) ], the intensityat point p ( i jk ) is finally given by: I i jk = I ( i jk )u e − P N ref m = ∆ τ m + S ( i jk )u N ref X m = ˜ α m e − P N ref n = m + ∆ τ n + S ( i jk )p N ref X m = ˜ β m e − P N ref n = m + ∆ τ n + S ( i jk )d N ref X m = ˜ γ m e − P N ref n = m + ∆ τ n , (20)where the upwind and current points always correspond to theindices m = m = N ref , respectively, and the sum over m isperformed over N ref − ff erent co-e ffi cients though. To solve the discretized equation of radiative transfer, the opac-ities χ C(u , d) , ¯ χ L(u , d) , source functions S C(u , d) , S L(u , d) , and velocityvectors V (u , d) , are required at the upwind and downwind points,together with the incident intensity, I u . We emphasize that thesubscript C describes continuum quantities, and should not beconfused with the subscript c denoting the control points of theinterpolation scheme.All required quantities are obtained from a 2D B´ezier inter-polation (see Appendix C) on the surfaces that intersect with agiven ray. The intersection surfaces depend on the considered di-rection and the size of the upwind and downwind grid cells. Fora given direction n = n x n y n z = sin θ cos φ sin θ sin φ cos θ , (21)
4. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
PSfragreplacements A u B u C u D u E u F u G u H u I u J u K u L u M u N u O u P u Q u R u S u n A d B d C d D d E d F d G d H d I d J d K d L d M d N d O d P d Q d R d S d n p ( ijk ) E u F u N u O u N d O d H d I d α x βyγ z n u ( ijk ) p ( ijk ) d ( ijk ) PSfragreplacements A u B u C u D u E u F u G u H u I u J u K u L u M u N u O u P u Q u R u S u n A d B d C d D d E d F d G d H d I d J d K d L d M d N d O d P d Q d R d S d n p ( ijk ) E u F u N u O u N d O d H d I d α x βyγ z n u ( i jk ) p ( i jk ) d ( i jk ) PSfragreplacements A u B u C u D u E u F u G u H u I u J u K u L u M u N u O u P u Q u R u S u n A d B d C d D d E d F d G d H d I d J d K d L d M d N d O d P d Q d R d S d n p ( ijk ) E u F u N u O u N d O d H d I d α x βyγ z n u ( i jk ) p ( i jk ) d ( i jk ) Fig. 1.
Geometry of the SC method for a particular ray with di-rection n propagating from the upwind point u ( i jk ) to a consid-ered grid point p ( i jk ) . The downwind point d ( i jk ) is required to setthe slope of a B´ezier curve representing the opacities and sourcefunctions along the ray. The middle and lower panel display allpossible downwind and upwind intersection surfaces for a shortcharacteristic at a grid point p ( i jk ) . For rays intersecting the x y -, xz -, and y z -planes, the 2D B´ezier interpolation is obtained fromgiven quantities at grid points located in the cyan, red, and ma-genta shaded surfaces, respectively. The coordinate system is in-dicated at the upper left, where α , β , γ determine the direction ofthe coordinate-axes and are defined in Sect. 3.3.where θ is the co-latitude (measured from the Cartesian z -axis),and φ is the azimuth (measured from the x -axis), the distances from a considered grid point to the neighbouring x y -, xz -, and y z -planes are calculated from trigonometry and yield: ∆ s (u) x y = z k − z k − γ n z ∆ s (d) x y = z k + γ − z k n z ∆ s (u) xz = y j − y j − β n y ∆ s (d) xz = y j + β − y j n y ∆ s (u) y z = x i − x i − α n x ∆ s (d) y z = x i + α − x i n x , with α, β, γ set to ± n x , n y , n z ≷
0, respectively. The intersection surface on the upwindand downwind side are then found at the minimum of ∆ s (u , d) x y , ∆ s (u , d) xz , ∆ s (u , d) y z , and the corresponding coordinates are eas-ily calculated.For each surface, the interpolation requires nine pointswithin the corresponding plane (see Fig. 1 and Eq. (C.1)). Ineach considered plane, we generally use grid points runningfrom (index-2) to (index) to determine upwind quantities, whiledownwind quantities are calculated from (index-1) to (index + Λ -matrix elements (see Appendix D). In Fig. 1, we show an ex-ample for a ray intersecting the x y -plane at the upwind side.The 2D B´ezier interpolation for the upwind point then con-sists of three 1D B´ezier interpolations along the x -axis usingthe points (J u , K u , L u ), (D u , E u , F u ), (M u , N u , O u ), followed byanother 1D B´ezier interpolation along the y -axis at the upwind x -coordinates. With the 2D B´ezier interpolation given by Eq. (C.1),we find for each required quantity q u , d : q ( i jk )u = w ( i jk )A q i − α, j − β, k − γ + w ( i jk )B q i − α, j − β, k − γ + w ( i jk )C q i , j − β, k − γ + w ( i jk )D q i − α, j − β, k − γ + w ( i jk )E q i − α, j − β, k − γ + w ( i jk )F q i , j − β, k − γ + w ( i jk )G q i − α, j − β, k + w ( i jk )H q i − α, j − β, k + w ( i jk )I q i , j − β, k + w ( i jk )J q i − α, j − β, k − γ + w ( i jk )K q i − α, j − β, k − γ + w ( i jk )L q i , j − β, k − γ + w ( i jk )M q i − α, j , k − γ + w ( i jk )N q i − α, j , k − γ + w ( i jk )O q i , j , k − γ + w ( i jk )P q i − α, j − β, k − γ + w ( i jk )Q q i − α, j , k − γ + w ( i jk )R q i − α, j − β, k + w ( i jk )S q i − α, j , k + w i jk q i jk (22) q ( i jk )d = ˜ w ( i jk )A q i − α, j + β, k − γ + ˜ w ( i jk ) B q i , j + β, k − γ + ˜ w ( i jk )C q i + α, j + β, k − γ + ˜ w ( i jk )D q i − α, j + β, k + ˜ w ( i jk )E q i , j + β, k + ˜ w ( i jk )F q i + α, j + β, k + ˜ w ( i jk )G q i − α, j + β, k + γ + ˜ w ( i jk )H q i , j + β, k + γ + ˜ w ( i jk )I q i + α, j + β, k + γ + ˜ w ( i jk )J q i − α, j − β, k + γ + ˜ w ( i jk )K q i , j − β, k + γ + ˜ w ( i jk )L q i + α, j − β, k + γ + ˜ w ( i jk )M q i − α, j , k + γ + ˜ w ( i jk )N q i , j , k + γ + ˜ w ( i jk )O q i + α, j , k + γ + ˜ w ( i jk )P q i + α, j − β, k − γ + ˜ w ( i jk )Q q i + α, j , k − γ + ˜ w ( i jk )R q i + α, j − β, k + ˜ w ( i jk )S q i + α, j , k , (23)where the coe ffi cients w ( i jk ) and ˜ w ( i jk ) refer to the upwind anddownwind interpolations corresponding to a considered point( i jk ). Depending on the intersection surface, ten out of these19 coe ffi cients are set to zero. For the upwind interpolation,we have already included the local coe ffi cient ( i jk ), which isonly required when boundary conditions need to be specified(Sect. 3.4). We note that all (non-zero) interpolation coe ffi cientsmay depend on the specific values of a considered quantity atthe given grid points, via the interpolation parameter ω to en-sure monotonicity. As in Sect. 3.1, also these monotonicity con-straints result in non-linear Λ -operators.
5. Hennicker et al.: 3D short-characteristics method in the winds from OB stars z x
PSfrag replacements n n n p(ijk)u u u i − , j , k − i − , j , k − i , j , k − i − , j , k Fig. 2.
Boundary conditions for rays propagating in the xz -planeat y -level ( j ) with three di ff erent directions n , n , n , and up-wind points u , u , u . For point u , the intensity is set to I u = I + c and all remaining quantities are obtained by bilinear interpo-lation from points ( i − , j , k − i , j , k − i − , j , k ), and( i jk ). The required quantities at point u are found from B´ezierinterpolation using the values at ( i − , j , k − i − , j , k − i , j , k − n , the unknown intensity inside the core, I − c , is directed inwards. Such situations occur only for ray direc-tions (nearly) parallel to the spatial grid, and thus are relativelyseldom. Since the inner boundary is usually not aligned with the 3DCartesian grid (e.g. a spherical star at the origin), the upwind(and downwind) interpolations need to be adapted near the stel-lar surface. For the upwind point, the following two situationsmay occur (see Fig. 2 for an example in the xz -plane). Firstly,the considered ray originates from the stellar surface (direction n in Fig. 2). In this case, we use a core-halo approximation andset I u = I + c = B ν ( T rad ), with I + c the emergent intensity from thecore, and T rad the radiation temperature. Unless explicitly noted,we assume T rad = T e ff throughout this work. All other quan-tities are obtained from trilinear interpolation using the points(E u , F u , H u , I u , N u , O u , S u , p ( i jk ) ) in Fig. 1, where representativeestimates need to be defined at the core points (those points thatare located inside the star). In hydrodynamic simulations, theanalogue of these points are so-called ‘ghost points’. Secondly,the considered ray originates from a plane spanned by grid pointsthat are partially located inside the star (direction n in Fig. 2).Then, the interpolation is performed as in Sect. 3.3, using againrepresentative estimates at the core-points.Inside the core, we define I + c = S L = S C = B ν ( T rad ) and set I − c and all velocity components to zero, where I − c is the inwarddirected intensity which needs to be specified only in rare situa-tions (Fig. 2, direction n ). The opacities inside the star are foundby extrapolation from the known values outside the star. Whilethis procedure certainly introduces errors (e.g. by over- and un-derestimating the upwind source function in optically thin and Table 1.
Mean relative error (defined in Sect. 4.2) of the meanintensity for a zero-opacity model, obtained using the Lebedev,Gauss-Legendre, and trapezoidal integration method, the latterwith nodes from Lobel & Blomme (2008). For ∆ J ex and ∆ J SC ,the incident intensities have been calculated exactly and from theSC method using linear upwind interpolations, respectively. Trapez Legendre Lebedev N Ω ∆ J ex [%] 12.1 8.0 14.3 9.4 11.2 7.7 ∆ J SC [%] 11.3 10.8 10.9 10.9 10.8 10.8 thick cases, respectively), it is still favourable to extrapolatingall values directly onto the stellar surface, mainly due to per-formance reasons . In addition to the error introduced by thepredefined values inside the core, the calculation of ∆ s r is a cer-tain issue, where ∆ s r is the distance of the current grid point tothe stellar surface. Since the radiative transfer near the stellarsurface is (in most cases) very sensitive to the path length of aconsidered ray, ∆ s r needs to be known exactly. Depending on theshape of the surface, ∆ s r can be calculated analytically, or needsto 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), usingthe estimates at the core points as defined above when necessary. To obtain the mean intensity at each grid point in the atmosphere,we solve the discretized equation of radiative transfer for manydirections and numerically integrate via: J i jk = π Z I i jk d Ω = X l w l I i jk ( Ω l ) , (24)where w l is the integration weight corresponding to a considereddirection Ω l = ( θ l , φ l ). The angular integration is particular chal-lenging for optically thin atmospheres, since in such situationseach (spatial) grid point is illuminated by the stellar surface,and the distribution of intensities I i jk ( θ, φ ) becomes a 2D step-function in the θ − φ -plane (if no upwind interpolation errorswere present). Depending on the considered position, the shapeof I i jk ( θ, φ ) greatly varies. Thus, elaborate integration methodsare required to resolve the star and its edges at any point of theatmosphere.Lobel & Blomme (2008) use the trapezoidal rule with a de-creasing number of polar grid points at higher latitudes to rea-sonably distribute the direction vectors on the unit sphere. Forthe 3D FVM, we have shown in paper I that a Gauss-Legendreintegration performs (slightly) better. However, the correspond-ing directions are always clustered in certain regions since thenodes of the Gauss-Legendre quadrature are fixed. Additionally,the Gauss-Legendre integration should only be applied whenthe distribution of intensities can be described by high orderpolynomials, that is, when I i jk ( θ, φ ) is smoothed out (e.g. bynumerical di ff usion). When numerical di ff usion errors are sup- For a given grid point, the number of neighbouring grid points thatcan be used for extrapolation is not a priori clear, and depends on theshape of the stellar surface, and the considered direction of the ray.Indeed, there are 64 special cases that would have to be implementedexplicitly. This is computationally not feasible.6. Hennicker et al.: 3D short-characteristics method in the winds from OB stars pressed (e.g. by using elaborate upwind interpolation schemes),the Gauss-Legendre integration should not be used (see Table 1).We have tested a multitude of other quadrature schemes,including trapezoidal and (pseudo)-Gaussian rules ontriangles, and the so-called Lebedev quadrature (see,e.g. Ahrens & Beylkin 2009, Beentjes 2015, and refer-ences therein). The Lebedev quadrature is optimized to exactlyintegrate the spherical harmonics up to a certain degree, with a(nearly) optimum distribution of direction vectors on the unitsphere. In Table 1, we summarize the errors for an optically thinatmosphere using di ff erent integration methods. The incidentintensities have been obtained exactly (i.e. by setting I i jk = I c and I i jk = N Ω ≈ N Ω ≈ N Ω = N Ω = J i jk = π Z d Ω Z x (max)obs x (min)obs d xI i jk Φ ( i jk ) x = X l w l X x w x I i jk Φ ( i jk ) x , (25)with x (min)obs and x (max)obs the required frequency shift in the ob-server’s frame obtained from the maximum absolute velocityoccurring in the atmosphere (see also paper I), and w x the cor-responding frequency integration weight. To resolve the pro-file function at each point in the atmosphere, we demand that | ∆ x obs | /δ . /
3. Since the profile function depends on the ratioof fiducial to actual thermal width, the fiducial velocity should beset to the minimum thermal velocity present in the atmosphere. Λ -iteration In Sect. 3.1, we have already noted that the Λ -operator becomesnon-linear due to monotonicity constraints (implemented by theinterpolation parameter ω ). In this section, we present a suitableworkaround, beginning with a recapitulation of some fundamen-tal ideas. Λ -matrix elements With the discretized equation of radiative transfer,Eqs. (12) / (20), and the upwind and downwind quantitiesobtained from Eqs. (22) and (23), the intensity at each spatial,angular, and frequency grid point can be calculated for a givensource function. We use the standard Λ -formalism to write theformal solution of the intensity, mean intensity, and scattering integral as: I = Λ Ω ,ν [ S C , L ] (26) J = Λ ν [ S C ] (27)¯ J = Λ [ S L ] , (28)with subscripts Ω and ν defining the dependence of the Λ -operator on direction and frequency, respectively. In the fol-lowing, we focus on the line transport. The continuum can bederived analogously. When all interpolation parameters ω havebeen determined (for a given stratification of source functionsand intensities), the Λ -operator is an a ffi ne operator describedby the Λ -matrix and a constant displacement vector Φ B repre-senting the propagation of boundary conditions (for a detaileddiscussion, see paper I, Puls 1991, and references therein). The Λ -matrix elements can then be obtained by: Λ m , n = ¯ J m ( S L = e n , Φ B = , (29)with the n -th unit vector e n , and matrix indices m , n related to the3D indices ( i , j , k ) by m = i + N x ( j − + N x N y ( k − , (30)where N x and N y denote the number of spatial grid points of the x and y coordinate, respectively. Eq. (30) simply transforms adata cube to a 1D array. The m , n -th matrix-element describesthe e ff ect of a non-vanishing source function at grid point n ontogrid point m . We emphasize that Eq. (29) holds only for pre-calculated interpolation parameters ω , obtained from an alreadyknown stratification of source functions. Λ -iteration The classical Λ -iteration scheme is defined by calculating a for-mal solution for a given source function using Eq. (28), followedby the calculation of a new source function by means of Eq. (3).For optically thick, scattering dominated atmospheres, however,this iteration scheme su ff ers from severe convergence problems(see Fig. 5 for the convergence behaviour of spherically sym-metric test models). To overcome these problems, we apply anaccelerated Λ -iteration scheme based on operator-splitting meth-ods (Cannon 1973). Within the ALI, the Λ -operator is written as Λ = Λ (A) + ( Λ − Λ (A) ) , (31)where the first term is an appropriately chosen ALO acting onthe new source function, S ( k )L , and the second term acts on theprevious one, S ( k − . For the converged solution, this scheme be-comes an exact relation. Using also, and in analogy to the ex-act Λ -operator, an a ffi ne representation for the approximate one, Λ (A) [ S ] = Λ ∗ · S + Φ B (cf. above), and evaluating Λ (A) at theprevious iteration step, k −
1, we obtain: S ( k )L = ζ · ¯ J ( k ) + Ψ ≈ ζ · Λ (A) k − [ S ( k )L ] + ζ · ( Λ k − − Λ (A) k − )[ S ( k − ] + Ψ = ζ · (cid:16) Λ ∗ k − S ( k )L + Φ ( k − + ¯ J ( k − − Λ ∗ k − S ( k − − Φ ( k − (cid:17) + Ψ . (32)Here, the iteration indices k − k are indicated as sub- or su-perscripts, ζ : = − ǫ L is a diagonal matrix, and Ψ : = ǫ L · B ν ( T ) isthe thermal contribution vector. From Eq. (32), it is obvious thatthe Φ B terms cancel. For multi-level atoms, we emphasize that
7. Hennicker et al.: 3D short-characteristics method in the winds from OB stars Λ and Λ (A) may change within the ALI-cycle due to the varia-tion of opacities (induced by the subsequently updated occupa-tion numbers). Furthermore, both operators might also changeeven for the simplified TLA approach considered in this paper,since the corresponding matrix elements depend on the sourcefunctions via the interpolation parameters ω k − , ω k . Rearrangingterms, we find:( − ζ · Λ ∗ k − ) S ( k )L = ζ · ( ¯ J ( k − − Λ ∗ k − · S ( k − ) + Ψ . (33)Eq. (33) is solved to obtain a new source function S ( k )L (see be-low). Since, however, Λ ∗ k − has been optimized only to ensuremonotonicity in a specific step k − S ( k − ), the iteration scheme can oscillate due to oscillations in Λ ∗ k and Λ ∗ k − . Even worse, the new source function might be-come negative. To overcome these problems, non-linear situ-ations need to be avoided (by providing almost constant Λ ∗ -matrices over subsequent iteration steps). The following ap-proach has proved to lead to a stable and convergent scheme:We apply purely linear interpolations ( ω k − = ω k = Λ ∗ k = Λ ∗ k − ) in the first four iteration steps to obtain analready smooth stratification of source functions. Additionally,we globally define a minimum allowed interpolation parameterand demand that ω > ω min . Then, ω becomes constant (namely ω = ω min ) in (most) critical situations, and again, Λ ∗ k approaches Λ ∗ k − . Whenever negative source functions or oscillations occurwithin the iteration scheme, ω min is gradually increased to one.With this approach, we obtain an always convergent iterationscheme, with a formal solution obtained by using linear interpo-lations only in most challenging cases. The rate of convergence achieved by the accelerated Λ -iterationscheme increases with the number of Λ -matrix elements in-cluded in the ALO. To minimize the computation time of thecomplete procedure, the choice of the ALO is always a compro-mise between the number of matrix-elements to be calculated,and the resulting convergence speed. While the inversion of adiagonal ALO reduces to a simple division, a multi-band ALOneeds to be inverted with some more e ff ort. When taking onlythe nearest neighbours into account, however, the ALO becomesa sparse matrix, and can be e ffi ciently inverted by applying theJacobi-Iteration for sparse systems (see paper I). For 3D calcula-tions, a multi-band ALO is required to obtain a rapidly converg-ing iteration scheme (see Hauschildt & Baron 2006 and paper Ifor solutions obtained with the LC method and the FVM, re-spectively), and thus implemented within our 3D SC framework.To calculate the corresponding Λ -matrix elements (includingall upwind and downwind interpolations), we extend the proce-dure developed by Olson & Kunasz (1987) and Kunasz & Olson(1988). A detailed derivation is given in Appendix D. Eqs. (D.1)-(D.27) correspond to the exact Λ -matrix elements for a localpoint and its 26 neighbours, and thus should give an excellentrate of convergence when included in the ALO (see, e.g. the26-neighbour ALO of P hoenix/ . For the simplified continuum and the TLA considered in this paper,the linear interpolation scheme is particularly advantageous in terms of
In this paper, we analyse the convergence speed of the ALIfor a diagonal ALO given by Eq. (D.14), a ‘direct-neighbour’(DN)-ALO given by Eqs. (D.5), (D.11), (D.13), (D.14), (D.15),(D.17), (D.23), and a ‘nearest-neighbour’ (NN)-ALO obtainedfrom all Eqs. (D.1)-(D.27). We note that only a moderate im-provement of the computation time can be expected when us-ing the diagonal or DN-ALO, since the diagonal and direct-neighbour elements depend on several other neighbours throughthe inclusion of downwind interpolations. Since, however, thedownwind-integration weight is generally negative, neglectingthese terms will overestimate the considered matrix elements,possibly resulting in a divergent iteration scheme. On the otherhand, when using purely linear interpolations (for the sourcecontribution and upwind interpolations), the calculation of theALO is greatly simplified since all coe ffi cients c i jk and w A , w B , w C , w D , w G , w J , w K , w L , w M , w P , w Q , w R vanish. For third-orderupwind / downwind interpolations as used in IRIS (Ibgui et al.2013), the calculation of the ALO coe ffi cients becomes com-putationally prohibitive at some point. Considering both inter-polation techniques used in this paper, the calculation of the di-agonal, DN-, and NN-ALO, in parallel to the formal solutionrequires 20%, 30%, and 40% of the total computation time.Finally, we have implemented the Ng-extrapolation (Ng1974, Olson et al. 1986) to improve the convergence speed fur-ther. As in paper I, the Ng-acceleration is applied in every fifthiteration step in order to use independent extrapolations. To minimize the computation time of our 3D code, we have im-plemented an elaborate grid construction procedure using a non-uniform grid-spacing that still enables a reasonably high spatialresolution (see Appendix A). Furthermore, when calculating linetransitions, we have parallelized the code using O pen
MP as inpaper I. The parallelization is implemented over the frequencygrid. We note that O pen
MP creates a local copy of the 3D arraysrepresenting the intensity and the (nearest-neighbour) Λ -matrix.With 27 Λ -matrix elements (per spatial grid point) included forthe ALO calculations, the (spatial) resolution becomes thereforememory limited. A typical resolution of N x = N y = N z = t (linear)SC ≈ t (B´ezier)SC ≈ ntel X eon X5650 (2.67 GH z ) machine with16 CPUs. As a reference, the FVM from paper I ‘only’ requiredroughly 50 minutes. A more meaningful comparison, however, isthe computation time per iteration, per CPU, and per angular andfrequency 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 ≈ . t (linear)SC ≈ . t (B´ezier)SC ≈ . / B´ezier interpolations are increased bya factor of roughly four / twelve, when compared to the 3D FVM.These di ff erences originate from the computationally more chal-lenging upwind / downwind interpolations, the integration of thediscretized equation of radiative transfer on (possibly) refinedgrids along a given ray, and from the calculation of an ALO in-cluding 26 neighbouring elements (instead of the 6 direct neigh-bours as used for the 3D FVM). computation time, since the corresponding ALO remains constant overall iteration steps, and therefore needs to be calculated only once.8. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
4. Spherically symmetric models
With the numerical tools developed in Sect. 3, we are able totackle 3D continuum and line scattering problems for arbitraryvelocity fields. In the following, we discuss the performance ofthe code when applied to spherically symmetric test models. Wecompare the solutions obtained from the 3D SC method usinglinear and B´ezier interpolations (hereafter denoted by SClin andSCbez, respectively), with those obtained from the 3D FVMand from accurate 1D solvers . Although we are using the samemodels as in paper I, the FVM solutions might di ff er slightlyfrom those presented in this paper, since we are using a di ff erentgrid, optimized for the 3D SC method. A first test of our 3D SC methods is the searchlight-beam test(e.g. Kunasz & Auer 1988). Within this test, we set the opacityto zero and consider the illumination of the atmosphere by a cen-tral star for a single direction. For consistency (cf. paper I), thedirection vector corresponds to θ = ◦ , φ = ◦ . Since the dis-cretized equation of radiative transfer, Eq. (12) / (20), reduces to I i jk = I ( u ) i jk , this test extracts the e ff ects of the applied interpola-tion schemes for the upwind intensity. The upper panel of Fig. 3shows the propagation of the specific intensity scaled by I c inthe xz -plane. Due to the upwind interpolation, the beam emerg-ing from the stellar core becomes widened. These interpolationerrors are connected with numerical di ff usion, and could only beavoided by applying a LC method. To obtain a quantitative mea-sure of this e ff ect, the lower and middle panel of Fig. 3 displaythe specific intensity along the given direction at the centre ofthe beam, and the specific intensity through a circular area per-pendicular to the ray direction as a function of impact parameter p . The corresponding exact solutions are given by a constant andrectangular function, respectively.Along the beam centre, the SC methods perfectly reproducethe exact solution, whereas the FVM solution decreases signif-icantly due to the finite grid-cell size (see paper I). Consideringthe intensity through the perpendicular area, both SC methodsperform better than the FVM, with slight advantages of theSCbez method when compared with the SClin method. Withinthe 3D SC methods, however, energy conservation is violatedfor our zero opacity models, because the (nominal) specific in-tensity jumps from I + c to zero for rays intersecting the stellarsurface or not, due to the core-halo approximation. As a con-sequence, almost all interpolations (and interpolation schemes)overestimate the specific intensity . For optically thick models(where I + c at the core plays a negligible or minor role), this ef-fect should decrease though. The associated error can be quanti-fied by calculating the corresponding flux, that is, by integratingthe specific intensity for a given direction over a correspond-ing perpendicular area (defined as a circle with virtually infiniteradius). We emphasize that the flux as defined here constitutesthe most demanding test case, and should not be confused with The 1D solution for the continuum transport is found from theRybicki-algorithm (combined with the solution of the moment equa-tions using variable Eddington factors, see, e.g. Mihalas 1978). To cal-culate the line, a comoving-frame ray-by-ray solution scheme in pz-geometry is applied, ensuring convergence with a diagonal ALO. In contrast, the number of photons entering and leaving a givengrid cell is (nearly) conserved within the FVM by definition. This state-ment, however, is not completely true for the FVM as formulated byAdam (1990), Lobel & Blomme (2008), Hennicker et al. (2018), sinceall these authors apply an (averaged) upwind approximation.
Fig. 3.
Searchlight-beam test for direction n = (1 , ,
1) and atypical grid with N x = N y = N z =
133 grid points. Upper panel:Contour plot of the specific intensity as calculated with the SCmethod using B´ezier interpolations in the xz -plane (cf. paper I,Fig. 3, for the finite-volume method). Middle panel: Specific in-tensity through the perpendicular area indicated by the straightline in the upper panel. The blue, red, and green profiles cor-respond 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 thecentre of the beam, whereas the FVM solution decreases signif-icantly for r & . R ∗ .
9. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. 4.
Photon flux as a function of direction angle φ (with fixed θ = ◦ ) through corresponding perpendicular areas, and withthe opacity set to zero. The central distance of all areas to thestellar surface has been set to 2 R ∗ . The same colour coding as inFig. 3 has been used.the flux density (i.e. the first moment of the specific intensity).Fig. 4 shows the resulting fluxes (normalized by the nominalvalue) for searchlight beams with di ff erent directions defined by θ = ◦ and φ ∈ [0 ◦ , ◦ ]. For di ff erent directions φ , the search-light beams propagate through di ff erent domains of the spatialgrid with accordingly di ff erent grid-cell sizes. Due to the dis-tinct behaviour of numerical di ff usion errors within these do-mains, the total flux varies as a function of φ . Overall, the totalfluxes for the SClin and SCbez methods are larger than theoreti-cally constrained, whereas the FVM gives (despite a small error)reasonable results. This e ff ect is largest in regions far from thestar and for diagonal directions. Thus, particularly in these re-gions, also the mean intensities (for optically thin atmospheres)are expected to be overestimated. The same problem arises whencalculating line transitions, since photons may freely propagateover large distances before a resonance region is hit. Numericaldi ff usion errors can only be avoided by increasing the grid reso-lution, or using higher order upwind-interpolation methods. In this section, we test the performance of the 3D SC methodwhen applied to spherically symmetric, stationary atmospheres.For consistency, the same test models as in paper I have beencalculated, with the wind described by a β -velocity law and thecontinuity equation: v ( r ) = v ∞ (cid:16) − b R ∗ r (cid:17) β b = − (cid:16) v min v ∞ (cid:17) /β ρ ( r ) = ˙ M π r v ( r ) . For stellar and wind parameters, R ∗ = R ⊙ , ˙ M = · − M ⊙ yr − , β = v min =
10 km s − , v ∞ = − , thedensity stratification and the velocity field are completely deter-mined. For the considered scattering problems ( ǫ C = ǫ L = − ),e ff ects of the temperature stratification are negligible. The con-tinuum and (frequency integrated) line opacities have been cal- culated from Eqs. (7) and (8), with the electron density derivedfor a completely ionized H / He plasma with helium abundance N He / N H = .
1. We have calculated three di ff erent continuummodels by scaling the opacity with k C = [1 , , τ r = [0 . , . , k L = [1 , , ]. To minimize the computation time,we use a microturbulent velocity v micro =
100 km s − through-out this paper. Such a large velocity dispersion mimicks the ef-fects of multiply non-monotonic velocity fields resulting fromthe line-driven instability (paper I and references therein). Theatomic mass has been set to m A = m p . Finally, the emergentintensity from the stellar core is calculated from T rad = T e ff =
40 kK.
Fig. 5 shows the maximum relative corrections of the mean in-tensity (left panel) and line source function (right panel) aftereach iteration step. Di ff erent methods (SClin, SCbez, and FVM),and di ff erent acceleration techniques (classical Λ -iteration, anddiagonal-, direct-neighbour-, nearest-neighbour-ALO with theNg-extrapolation switched on or o ff ) have been applied. We dis-play the continuum and line calculations for k C = [10 , k L = [10 , ], respectively, using N x = N y = N z =
93 spa-tial and N Ω =
974 angular grid points (to save computation timewhen calculating the slowly converging classical Λ -iteration). Inthe following, we only discuss the convergence behaviour of theSClin and SCbez methods, as the FVM has already been anal-ysed in paper I. We usually stop the iteration scheme when themaximum relative corrections become less than 10 − betweensubsequent iteration steps, emphasizing that a truly convergedsolution is only found when the curve describing subsequent rel-ative corrections is su ffi ciently steep . For instance, the classical Λ -iteration (with Λ ∗ =
0) fails to converge for strong scatteringlines (Fig. 5, lower right panel), since the relative corrections be-come almost constant in each iteration step (‘false convergence’,cf. Hubeny & Mihalas 2014).Overall, and as expected, the number of iterations neededto obtain the converged solution is decreasing with increasingnumber of matrix elements used to define the ALO (see alsoHauschildt et al. 1994 and Hauschildt & Baron 2006 for multi-band ALOs coupled to a 1D-SC and 3D-LC formal solutionscheme, respectively). In most cases, the convergence of theSClin is faster than that of the SCbez method, because the in-terpolation scheme is intrinsically more localized (with strongerweights assigned to local Λ -matrix elements). The FVM alwaysperforms best, since only the direct neighbours directly influencethe formal solution within this method.For parameters k C , k L ≤
10 (Fig. 5, upper panels), all ALOsyield a converged solution within N iter ≈
10 iteration steps.When applying the SCbez method, the first peak results fromswitching the linear interpolations to B´ezier interpolations at thefifth iteration step. This peak is less pronounced for parameters For linearly convergent iteration schemes, the steepness of the con-vergence curve is described by the relative corrections in subsequentiterations steps, ∆ k / ∆ k − = const . = : q. To obtain a solution within areasonable amount of computation time, we may demand that q . . − every30th iteration step.10. Hennicker et al.: 3D short-characteristics method in the winds from OB stars Fig. 5.
Convergence behaviour of the standard spherically symmetric model calculated with the 3D FVM (blue) and the 3D SCmethod using linear (red) or B´ezier (green) interpolations. The left and right panels show the convergence behaviour of the contin-uum and line transfer for ǫ C = ǫ L = − , 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 = . Di ff erent acceleration techniques have been applied,where the NN-ALO is implemented only for the SC method. k C , k L = , , since the maximum relative corrections arestill relatively large in the first few iteration steps.For the optically thick continuum model (Fig. 5, lower leftpanel), the number of iterations until convergence is reducedfrom N diagiter ≈
75 to N DNiter ≈
65 and N NNiter ≈
45 when usingthe diagonal, DN-, and NN-ALO within the SClin and SCbezmethod, respectively. The Ng-extrapolation significantly reduces N iter further, and is required to obtain the converged solution in .
20 iteration steps.For the strong line, the convergence behaviour becomesslightly improved, because the radiative transfer is localized to(narrow) resonance regions. Thus, already the diagonal and DN-ALO yield a solution within N diagiter ≈
50 and N DNiter ≈
35 iterationsteps. Again, the NN-ALO with the Ng-extrapolation schemeswitched on performs best, and reduces the number of iterationsteps until convergence to . In the following, we discuss the errors resulting from the up-wind and downwind interpolations, and from the integration ofthe discretized radiative transfer equation. We apply the NN-ALO together with the Ng-extrapolation to ensure convergence.When calculating the line, we used N x = N y = N z =
93 spatialgrid points. Since the continuum transport has been calculatedat only one frequency point, we applied a higher grid resolutionwith N x = N y = N z =
133 for such problems. Fig. 6 shows thecontinuum and line solutions together with corresponding rela-tive errors obtained for the spherically symmetric model whencalculated with the FVM, SClin, and SCbez methods, and com-pared to the ‘exact’ 1D solution. The mean and maximum rela-tive errors are shown for di ff erent regions in Table 2, where themean and maximum relative errors of any quantity are definedthroughout this paper by ∆ q : = N N X i = | q i − q (exact) i | q (exact) i ∆ q max : = max ∀ i ∈ [1 , N ] | q i − q (exact) i | q (exact) i ,
11. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. 6.
Solutions for the standard spherically symmetric models as calculated with the 3D FVM (blue) and 3D SC methods usinglinear (red) or B´ezier (green) interpolations, compared to an accurate 1D solution (black solid line). The dots represent the solutionsat specific grid points (with di ff erent latitudes and longitudes), where only a subset of all grid points is displayed to compress theimage. Corresponding errors are indicated at the bottom of each chart. The top panel shows the mean intensity for the continuumtransfer as a function of radius, with ǫ C = − , and k C = [1 , , ǫ L = − , and k L = [10 , , ] from left to right. Table 2.
Mean and maximum relative errors of the FVM andSClin, SCbez methods, when applied to spherically symmetrictest models. The mean relative errors are listed for di ff erent re-gions with r ∈ [ R ∗ , R ∗ ], r ∈ [3 R ∗ , R max ], and r ∈ [ R ∗ , R max ],from top to bottom. ∆ J [%] for r ∈ [ R ∗ , R ∗ ] ∆ S L [%] for r ∈ [ R ∗ , R ∗ ] k C τ r FVM SClin SCbez k L FVM SClin SCbez10 .
17 3.4 6.2 6.6 10 .
17 18 1.2 2.2 10 .
10 3.3 2.2 ∆ J [%] for r ∈ [3 R ∗ , R max ] ∆ S L [%] for r ∈ [3 R ∗ , R max ] k C τ r FVM SClin SCbez k L FVM SClin SCbez10 .
17 5.2 4.4 4.3 10 .
17 22 10 3.7 10 .
12 19 6.2 ∆ J [%] for r ∈ [ R ∗ , R max ] ∆ S L [%] for r ∈ [ R ∗ , R max ] k C τ r FVM SClin SCbez k L FVM SClin SCbez10 .
17 4.4 5.2 5.4 10 .
17 20 6.0 3.0 10 .
11 13 4.7 ∆ J max [%] for r ∈ [ R ∗ , R max ] ∆ S L , max [%] for r ∈ [ R ∗ , R max ] k C τ r FVM SClin SCbez k L FVM SClin SCbez10 .
17 17 51 49 10
27 46 4510 .
17 36 39 14 10
73 50 4510 .
70 65 55 with N the number of grid points within the considered region. For the continuum models, the solutions obtained from the3D SC methods are superior to the solution obtained from theFVM in most cases. Particularly near the stellar surface (at r . R ∗ ), both SC methods are in good agreement with the1D solution (see Fig. 6, top panel, and bottom of each chartfor the radial dependence of the relative errors). When consid-ering the most challenging problem of optically thick, scatteringdominated atmospheres, the mean relative errors of the SClinand SCbez method for the complete calculation region are onthe 20- and 10 %-level, respectively. For such models, the FVMbreaks down due to the order of accuracy (see paper I), and a(high order) SC method is indeed required to solve the radiativetransfer with reasonable accuracy. For marginally optically thickcontinua, the mean relative errors of the SClin and SCbez meth-ods are on the 5 %-level and below. While the FVM allows for aqualitative interpretation of the radiation field for such models,the SC methods should be used for quantitative discussions. Theoptically thin model calculations give mean relative errors of theorder of 5 % for all methods, with the maximum relative errorbeing lowest for the FVM. Since, additionally, the computationtimes of the SClin and SCbez methods are typically highest (seeSect. 3.7), the FVM is to be preferred when calculating opticallythin continua. We note that all errors originate from the inter-play between upwind and downwind interpolations of opacities,source functions, and intensities, and from the integration of thediscretized radiative transfer equation. Numerical di ff usion andthe associated violation of energy conservation influences theconverged solution particularly in the optically thin regime.The mean relative errors for the line transition are of the or-der of 5 −
10 %, with increasing accuracy from strong to weaklines, and slight advantages of the SCbez method when com-pared to the SClin method. The radial stratification of relativeerrors for each considered line is shown in the bottom panel of
12. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. 6, bottom of each chart. While the FVM gives largest er-rors near the stellar surface (at r . R ∗ ), both SC techniquesare in excellent agreement with the 1D solution in such regions.At larger radii, however, the SC solutions are generally over-estimated when compared to the 1D solution due to numericaldi ff usion errors. The distinct behaviour of the applied solutionschemes in di ff erent atmospheric regions finally determines thequality of emergent flux profiles. The converged source functions are used to calculate the emer-gent flux profile using the same postprocessing LC solver as inpaper I. Based on Lamers et al. (1987), Busche & Hillier (2005),Sundqvist et al. (2012), this method solves the equation of ra-diative transfer in cylindrical coordinates with the z -axis beingaligned with the line of sight of the observer’s direction underconsideration. All quantities required on the rays are found bytrilinear interpolation from the 3D grid, and the equation of ra-diative transfer is integrated using linear interpolations. To ex-tract the error resulting from the FVM and SC methods alone,we interpolated the ‘exact’ 1D solution onto our 3D grid andcalculated the reference profile using the same technique. Thecontinuum has been calculated from a zero-opacity model givenby the unattenuated illumination from the projected stellar disc.Then, the di ff erences of line profiles are exclusively related tothe di ff erences of line source functions. Fig. 7 shows the lineprofiles with corresponding absolute and relative errors for theintermediate ( k L = ) and strong ( k L = ) line, obtainedfrom the converged source functions from above.The line profiles are in good agreement with the 1D solu-tion for both applied SC methods, with slight advantages of theSCbez method when compared to the SClin. Major (relative) de-viations arise particularly at large frequency shifts on the blueside due to the enlarged source functions in corresponding res-onance regions (i.e. at large radii in front of the star). At suchfrequency shifts, however, the line profile is mainly controlledby absorption, and the absolute error remains small. At lowfrequency shifts, the emission peak becomes slightly overesti-mated, particularly when considering the strong line. The corre-sponding resonance regions are mainly located near to the star(at low absolute velocities) and in the whole plane perpendicularto the line of sight (with low projected velocities). For the in-termediate line, the emission from this plane at large radii onlyplays a minor role due to relatively small optical-depth incre-ments along the line of sight. Thus, both SC methods are in ex-cellent agreement with the 1D reference profile. With increasingline strength, however, the emission from the outer wind regioncontributes significantly to the line formation, and the discrep-ancies between the 1D and the SClin / SCbez methods becomemore pronounced. For all test calculations, the B´ezier methodperforms best, closely followed by the SClin method, and (farbehind) the FVM.With Fig. 7 and the argumentation from above, we concludethat (at least) a short-characteristics solution scheme is requiredto enable a quantitative interpretation of ultra-violet (UV) reso-nance line profiles, where both the linear and B´ezier interpola-tion techniques perform similarly well. The less accurate (how-ever computationally cheaper) FVM can still be applied for qual-itative discussions.
Fig. 7.
Emergent flux profiles of an intermediate ( k L = , toppanel) and strong ( k L = , bottom panel) line. The blue, red,and green curves correspond to the solution of the FVM, SClin,and SCbez methods, respectively. The reference profile (blacksolid line) has been derived from the ‘exact’ 1D source functioninterpolated onto the 3D Cartesian grid. Corresponding relativeand absolute errors are shown at the bottom of each chart. For allprofiles, the continuum level has been determined from a zero-opacity model.
5. Rotating winds
As a first application of the 3D SC method to non-sphericalatmospheres, we consider the UV resonance line formation inthe winds of (fast) rotating O stars. Fast rotation has two im-mediate consequences on the stellar geometry and wind struc-ture. Firstly, the surface of any rotating star becomes distorted,with R eq / R pole approaching 3 / Ω =
13. Hennicker et al.: 3D short-characteristics method in the winds from OB stars the emergent flux depends on the (local) e ff ective gravity (cor-rected for the centrifugal acceleration), and thus, decreases to-wards the equator (‘gravity darkening’, see von Zeipel 1924, andMaeder 1999, Maeder & Meynet 2000 for uniform and shellu-lar rotation, respectively). The first attempt to model the windsof fast rotating OB stars was made by Bjorkman & Cassinelli(1993). These authors considered a purely radial line force, andneglected gravity darkening and the surface distortion. Withinthese approximations, a ‘wind compressed disc’ is formed in theequatorial plane. Cranmer & Owocki (1995) and Owocki et al.(1996) included the e ff ects of non-radial line-forces into their2D radiation-hydrodynamic simulations, and showed that theformation of the disc becomes suppressed due to a small, butsignificant polewards acceleration, giving rise to an associatedpolar velocity component that prevents the formation of a disc.When also accounting for gravity darkening (i.e. a decreasedradial acceleration in equatorial regions), Owocki et al. (1996)further showed that a prolate wind structure develops, with de-creased equatorial mass loss and velocity (see also the reviewby Owocki et al. 1998). Maeder (1999) proposed that an oblatewind structure might still be possible, when accounting for apolar variation of the ionization equilibrium induced by grav-ity darkening (the so-called κ -e ff ect). This e ff ect becomes par-ticularly important when the local e ff ective temperature dropsbelow the bi-stability jump temperature . Petrenz & Puls (2000)extended the hydrodynamic calculations from above by allowingfor spatially varying line force multipliers, and showed that nomajor di ff erences from the prolate wind structure arise, at leastfor OB stars above T e ff &
20 kK with an optically thin Lymancontinuum. Recently, Gagnier et al. (2019) reinvestigated the ef-fects of rotation using 2D ESTER models (see Rieutord et al.2016 for a description of this code). Using a di ff erent imple-mentation of gravity darkening (consistent with the so-called ω -model by Espinosa Lara & Rieutord 2011, which basicallyresults in a slower decrease of e ff ective temperature with co-latitude than obtained from the von Zeipel theorem), these au-thors predict either a ‘single-wind regime’ (with enhanced po-lar mass loss) or a ‘two-wind regime’ (with enhanced mass lossat latitudes where the e ff ective temperature drops below the bi-stability jump temperature). To understand which of the di ff erentmodels represents reality best (in di ff erent temperature regimes),one needs to compare synthetic profiles with observations. Inthis respect, investigating the e ff ects of prolate and oblate windstructures is particularly important to distinguish between di ff er-ent theories.As a consequence of the distinct wind structure resultingfrom a particular model, the wind lines of rotating stars are ex-pected to di ff er as a function of rotational speed and inclination.To predict UV resonance line profiles of fast rotating stars, wecalculated the source function of a prototypical resonance linetransition including the e ff ects of gravity darkening and surfacedistortion for models with di ff erent rotational velocities. As afirst step, we used a wind description that is consistent with theprolate wind model. For all calculations, we applied the SClinmethod. Fig. 8.
Contours of the density in the xz -plane ( z being the ro-tation axis) for a prototypical rotating O star with L ∗ = L ⊙ , M ∗ = . M ⊙ , R p = R ⊙ , and v rot =
432 km s − (correspond-ing to Ω = . Ω =
0. While the thick solid linecorresponds to the surface of the (distorted) star, the dashed linewould correspond to a spherical surface. Additionally, the veloc-ity vectors are displayed.
Table 3.
Specific parameters used and obtained for the rotatingwind models. For a given stellar luminosity L ∗ = L ⊙ , stellarmass M ∗ = . M ⊙ , and polar radius R p = R ⊙ , rows twoto eight display the rotation parameter Ω , the equatorial radius R eq , the polar and equatorial e ff ective temperature T e ff , p , T e ff , eq ,the total mass loss rate ˙ M , and the polar and equatorial terminalvelocity v ∞ , p , v ∞ , eq , for di ff erent equatorial rotation speeds v rot . v rot [km s − ] 0 210 294 432 Ω R eq [ R p ] 1 1.04 1.09 1.22 T e ff , p [kK] 41.84 42.44 43.07 44.66 T e ff , eq [kK] 41.84 40.61 39.28 35.20˙ M [10 − M ⊙ yr − ] 2.70 2.73 2.79 2.93 v ∞ , p [km s − ] 2781 2989 3255 3159 v ∞ , eq [km s − ] 2781 2651 2556 2273 To obtain a model for the structure of rotating winds, we applieda 2D version of the VH-1 code developed by J. M. Blondin andco-workers. Our model includes the e ff ects of gravity darken-ing and surface distortion (see below). Using a 1D input modelderived from radiation driven wind theory including finite coneangle corrections (Castor et al. 1975 and Pauldrach et al. 1986)for the first time step, the radiation hydrodynamic equations (ac-counting for non-radial line forces) are solved until a (quasi)stationary solution is obtained (see Cranmer & Owocki 1995and Owocki et al. 1996 for the description of the line force).Assuming azimuthal symmetry, the resulting 2D density and ve-locity structure is then used as input for our 3D SC code. Table 3summarizes specific parameters used and obtained for our model The jump temperature is theoretically motivated by a stronger ra-diative line-driving due to lower ionization stages of iron for T e ff . T jump ≈
25 kK (Vink et al. 1999). More recently, Petrov et al. (2016)predicted a somewhat lower jump temperature, T jump ≈
20 kK. http://wonka.physics.ncsu.edu/pub/VH-1/
14. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
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 di ff erent rotation parameters Ω .calculations. While the surface integrated mass flux, ˙ M , becomesonly slightly increased with increasing rotational speed, the po-lar (equatorial) terminal velocities are significantly enhanced(reduced). For the fastest rotating model ( v rot =
432 km s − ),Fig. 8 shows corresponding density contours in the xz -plane.The z -axis is aligned with the rotation axis. To explicitly showthe prolate wind structure, we have scaled the density by thedensity resulting from the non-rotating (spherically symmetric)model, as a function of distance from the stellar surface. Fordi ff erent rotational speeds, Fig. 9 displays the density and ra-dial velocity along the polar axis and along an (arbitrarily de-fined) axis in the equatorial plane. When compared with thespherically symmetric wind, the densities of the rotating mod-els are enhanced in polar regions, and become reduced along theequator. Further, the radial velocity along the polar axis remainsnearly the same, except in regions far from the star, where theterminal velocity of all rotating models becomes enhanced withincreasing v rot . We note that one would expect clearer di ff erencesof the (radial) velocity fields for di ff erent rotation rates, due todi ff erent accelerations induced, particularly, by the di ff erent ra-diative fluxes resulting from gravity darkening, and, though toa lesser extent, by the specific density structure and rotationalvelocities. Such di ff erences can be barely observed within oursimulations, most presumably because the wind structure in po-lar regions has not completely settled to a stationary state at thelast time steps. In contrast, the radial velocity in equatorial re-gions is significantly reduced at all distances, when compared tothe non-rotating wind, and the deviations from spherical symme-try become more pronounced with increasing rotational velocity.Although we have averaged the hydrodynamic simulations overthe last 20 time steps, the atmospheric structure still su ff ers fromsmall numerical artefacts.To calculate the stellar surface distortion, we consider thegravitational potential of the star accounting for the e ff ects ofcentrifugal forces. Under reasonable assumptions, we can ap-proximate this potential by a Roche model (e.g. Collins 1963,see also Cranmer & Owocki 1995): Φ ( r , Θ ) = − GM ∗ r − ω r sin ( Θ )2 , (34)with angular velocity ω . The surface of the star is defined onequipotential lines and can be calculated by setting Φ ( R p , Θ = = Φ ( R ∗ ( Θ ) , Θ ), with R p the polar radius. Solving the resultingcubic equation, one finds: R ∗ ( Ω , Θ ) = R p Ω sin( Θ ) cos (cid:16) π + cos − (cid:0) Ω sin( Θ ) (cid:1) (cid:17) , (35)with Ω = ω/ω crit the ratio of the actual to the critical (‘breakup’)angular velocity. Defining the rotational speed at the equator v rot as input parameter, one easily obtains (cf. Cranmer & Owocki 1995, their Eq. 27): Ω = v rot R eq ω crit (36) R eq = R p − v R p / GM ∗ , (37)with equatorial radius R eq . Following Maeder & Meynet (2000),we use the actual stellar mass to calculate the equatorial ra-dius and critical velocity without correcting for Thomson-accelerations. Additionally, we note that our stellar models arewell below the Eddington limit ( Γ = . ω crit = (8 GM ∗ / R ) / .Instead of using Eq. (35) in our final implementation, we ap-proximated the stellar surface by a spheroid with semi-majoraxes a = b = R eq and semi-minor axis c = R p . Such a for-mulation greatly simplifies the calculation of the intersection ofa given ray with the stellar surface (required for the boundaryconditions, see Sect. 3.4). For the most extreme test case consid-ered here ( Ω = . R ∗ ( Θ ) due to thisapproximation is well below the 2%-level, and rapidly decreasesfor lower rotational velocities.To calculate the intensity emerging from the stellar core, weset I + c ( Θ ) = B ν ( T e ff ( Θ )), with the e ff ective temperature as a func-tion of co-latitude. For a given luminosity of the star, L ∗ , we ob-tain (see also Petrenz & Puls 1996): T e ff ( Θ ) = h L ∗ πσ B Σ | g | β Z i / (38) Σ = Z π | g | β Z R ∗ ( Θ ) sin( Θ ) − g r / | g | d Θ , with σ B the Stefan Boltzmann constant, and the surface inte-grated e ff ective gravity Σ derived from g ( Θ ) = − ∇ Φ ( R ∗ ( Θ ) , Θ ).The parameter β Z describes the gravity darkening law in termsof T e ff ( Θ ) ∝ | g ( Θ ) | β Z . As originally formulated by von Zeipel(1924), β Z = /
4. Though β Z might be significantly lower(e.g. Domiciano de Souza et al. 2014, Gagnier et al. 2019), forsimplicity we nevertheless used β Z = /
4. As long as we as-sume constant ionization fractions, the e ff ect of β Z on the lineprofiles will be minor anyway. For our test models, we used v micro =
100 km s − , and calculatedthe frequency integrated opacity from Eq. 8 for an intermedi-ate and a strong line with line-strength parameter k L = and k L = . To obtain the source function, we applied the 3D SClinmethod and set ǫ L = − . The resulting (normalized) line pro-files are shown in Fig. 10 for di ff erent rotational velocities andinclination angles. Additionally, we display the continuum flux
15. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. 10.
Predicted emergent flux profiles for the rotating star models with
Ω = [0 , . , . , .
9] (see Table 3). The upper and lowerpanels show the intermediate and strong line with k L = and k L = , respectively. The inclination angle has been set tosin( i ) = [0 , . ,
1] from left to right.
Fig. 11.
Continuum fluxes for di ff erent inclinations sin( i ) = [0 , . ,
1] from left to right, and di ff erent rotation parameters(using the same colour coding as in Fig. 10). The continuumfluxes have been scaled by the corresponding flux obtained fromthe non-rotating model.used for normalization in Fig. 11. Due to gravity darkening andthe surface distortion, the continuum depends on the rotation rateand inclination, with largest fluxes found for high rotation ratesand low inclinations (resulting from the higher temperatures inpolar regions and a larger projected stellar disc). In Figs. 10 and11, the x -axes have been normalized to an (arbitrarily chosen)terminal velocity v ∞ = − . The behaviour of the line profiles can be qualitatively explained with the hydrodynamicstructure:(i) For pole on observers (sin( i ) =
0, left panel of Fig. 10), theabsorption column in front of the star is enhanced with increas-ing rotational velocity due to the larger densities (and opacities)in polar regions. Thus, the absorption trough (of unsaturatedlines) becomes more pronounced. The absorption edge of theintermediate lines is found at slightly lower velocities than ex-pected from the hydrodynamic simulations, because the opticaldepths of the corresponding resonance regions are too low to e ffi -ciently contribute to the absorption. When considering the stronglines, the optical depth is larger, and the absorption edge is moreconsistent with the actual terminal velocity. For both applied linestrength parameters, the emission peak becomes weaker with in-creasing rotation rate, particularly at low projected velocities onthe red side of the line centre (for negative x obs ). This e ff ect canbe partly explained by the reduced emission from the equato-rial plane, due to the lower densities in these regions. More im-portantly, however, is the increased continuum flux that mainlydetermines the (normalized) height of the emission peak.(ii) When increasing the inclination towards equator-on ob-servers (sin( i ) =
1, right panel of Fig. 10), the behaviour is re-versed. For such directions, the continuum plays an only minorrole, since the lower (equatorial) e ff ective temperatures of the ro-tating models are nearly compensated by the enlarged projectedstellar disc. With increasing rotation parameter, the absorptiontrough of the intermediate line becomes reduced and shifted to-wards lower terminal velocities, consistent with the hydrody-namical model. When considering the strong line, the absorptionbecomes saturated, and only the impact of the di ff erent termi-
16. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. 12.
Predicted emergent flux profiles for the rotating starmodels with
Ω = [0 . , . , . , . , . ff erent inclina-tion angles such that v sin( i ) =
200 km s − for all models. Theupper and lower panel display the intermediate ( k L = ) andstrong ( k L = ) line, respectively.nal velocities can be observed. Additionally, and for both linestrengths, the absorption slightly extends towards the red side,because of (negative) projected line of sight velocities near thestellar surface induced by rotation. For the fastest rotating modelwith Ω = .
9, this e ff ect becomes suppressed due to an increasedemission from the (dense) prolate wind structure. This latter ef-fect is only moderate for lower rotational speeds. Based on thecurrent hydrodynamic wind structure, we would therefore ex-pect to observe either rather low terminal velocities or relativelydeep absorption troughs for fast rotating stars, and we are able,at least in principle, to check the theory by comparing our syn-thetic 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 ac-tual rotational velocity from UV resonance lines. This latterpoint becomes clear from Fig. 12, where we predict the line pro-files of models with di ff erent rotational speed for a given v sin( i )(set here to 200 km s − ). Since, at least for the intermediate line,the profile shapes di ff er, sin( i ) could be derived if v sin( i ) wasknown. Of course, such constraints will become feasible onlyif the underlying models correctly describe the wind structure(including possibly varying ionization stages) of rotating stars.
6. Summary and conclusions
In this study, we have presented a 3D short-characteristicsmethod tailored for the solution of continuum- and line-scattering problems in the winds of hot stars. To obtain the for-mal solution, we have implemented a purely linear interpolationscheme (for calculating upwind quantities and for the solutionof the radiative transfer equation along a ray), as well as a sec-ond order, monotonic, B´ezier technique. We use Cartesian coor-dinates with a non-uniform grid spacing to ensure a reasonablespatial resolution in important regions (i.e. where velocity and / ordensity gradients are large). As a first step towards full NLTEradiative transfer models, we consider a single resonance-linetransition (approximated by a two-level-atom) assuming an op-tically thin background continuum, whereas for pure continuumproblems we use the thermalization parameter, ǫ C , and split thesource function into a scattering and a thermal part. A general-ization (including multi-level atoms) is planned for future appli-cations.To calculate strong scattering lines and optically thick, scat-tering dominated continua, we have implemented an accelerated Λ -iteration scheme using di ff erent non-local approximate Λ -operators (ALOs), together with applying the Ng-extrapolationmethod for subsequent iterations. With increasing complexityof the ALO (i.e. from a purely diagonal ALO to a nearest-neighbour ALO including also the 26 neighbouring terms), therate of convergence is improved. When applying the NN-ALO,the converged solution is generally found within .
20 iterationsteps even for the most challenging test cases.We have estimated the error of the applied methods in dif-ferent regimes by calculating spherically symmetric test mod-els within our 3D SC framework, and with a 3D finite-volumemethod. To our knowledge, this is the first study, where di ff er-ent 3D solution schemes for spherical problems have been com-pared, and their precision explored. When rated against the so-lution obtained from (accurate) 1D solvers, we found a meanrelative error for the converged continuum source function ofroughly 5 −
10 % and 5 −
20 % when using B´ezier and linearinterpolations, respectively. Particularly for optically thick con-tinua, the (first order) FVM method breaks down, and a (highorder) SC or LC method is required to accurately solve the ra-diative transfer. When considering the solution of the line sourcefunction for di ff erent line-strength parameters, the mean relativeerrors of both SC methods are on the 10 %-level and below, withslight advantages of the B´ezier technique compared to purelylinear interpolations. The resulting synthetic line profiles are cal-culated with a long-characteristics postprocessing routine usingthe previously calculated converged source functions. The SCmethod using B´ezier interpolations almost perfectly matches the1D reference profiles for all our models (i.e. for weak and stronglines). When linear interpolations are used, we obtain slight devi-ations originating mainly in the outer wind regions. In contrast,the 3D FVM always overestimates the emission. Nevertheless,
17. Hennicker et al.: 3D short-characteristics method in the winds from OB stars all methods have their own advantages and disadvantages, par-ticularly when also accounting for the computation time (withfastest turn-around times for the FVM method). Thus, the 3DFVM method should be used for qualitative interpretations andfor finding (rough) estimates of the parameters of interest, whilethe SC methods are to be preferred when aiming at a quantitativeanalysis of line profiles, and for optically thick continua.As a first application of the 3D SC code to non-sphericalproblems, we considered the resonance line formation in thewinds of (fast) rotating O stars. Assuming azimuthal symme-try, the hydrodynamic structure for a prototypical O star withdi ff erent rotation rates has been calculated by means of the 2DVH-1 code. We have included the e ff ects of surface distortionand gravity darkening into our 3D radiative transfer framework.Given the hydrodynamic models, we are able to predict the shapeof line profiles for di ff erent rotational speeds and inclination an-gles. When compared with a (non-rotating) spherically symmet-ric wind (obtained using the same stellar parameters), rotatingstars should either show relatively low terminal velocities (forequator-on observers) or deeper absorption troughs (for pole-on observers). The latter e ff ect, however, would only be ob-servable when considering intermediate (i.e. unsaturated) lines.Additionally, we showed that the line profile shapes vary as afunction of rotational speed at a given v sin( i ). Thus, assumingthat v sin( i ) was known (e.g. from photospheric line modelling),one could estimate the rotational speed, though with a ratherlarge uncertainty, since the di ff erences of the line profiles areonly moderate. We emphasize that other e ff ects (such as clump-ing or a flatter gravity darkening law) may additionally shapethe line profiles. When analysing UV observations of fast ro-tating stars, the 3D SC code developed in this work certainlywill help to understand the manifestations of various (afore-mentioned) e ff ects, and to distinguish between di ff erent theo-retical predictions (e.g. prolate vs. oblate wind structures). Idealtestbeds for future investigations of fast rotating winds are thestars VFTS102 (O9 Vnnne, Dufton et al. 2011) and VFTS285(O7.5 Vnnn, Walborn et al. 2012), both rotating at nearly theircritical velocity.Finally, we note that our tools are, of course, not limited torotating stars. Indeed, almost any kind of stellar wind that de-viates from spherical symmetry (with non-relativistic velocityfields), such as magnetic winds or colliding winds in close bina-ries, can be investigated. Acknowledgements.
We thank our referee, Dr. Achim Feldmeier, for many help-ful comments and suggestions, that improved the first version of this manuscriptconsiderably. LH gratefully acknowledges support from the German ResearchFoundation, DFG, under grant PU 117 / / / References
Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters,116, 061102Adam, J. 1990, A&A, 240, 541Ahrens, C. & Beylkin, G. 2009, in Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, Vol. 465Auer, L. 2003, in Astronomical Society of the Pacific Conference Series, Vol.288, Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner,3Beentjes, C. H. L. 2015, Quadrature on a Spherical SurfaceBjorkman, J. E. & Cassinelli, J. P. 1993, ApJ, 409, 429Busche, J. R. & Hillier, D. J. 2005, AJ, 129, 454Cannon, C. J. 1973, ApJ, 185, 621Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157Cherepashchuk, A. M. 1976, Soviet Astronomy Letters, 2, 138 Collins, II, G. W. 1963, ApJ, 138, 1134Cranmer, S. R. & Owocki, S. P. 1995, ApJ, 440, 308David-Uraz, A., Erba, C., Petit, V., et al. 2019, MNRAS, 483, 2814de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ,764, 166Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569,A10Dufton, P. L., Dunstall, P. R., Evans, C. J., et al. 2011, ApJ, 743, L22Dullemond, C. P. & Turolla, R. 2000, A&A, 360, 1187Espinosa Lara, F. & Rieutord, M. 2011, A&A, 533, A43Gagnier, D., Rieutord, M., Charbonnel, C., Putigny, B., & Espinosa Lara, F.2019, A&A, 625, A88Georgiev, L. N., Hillier, D. J., & Zsarg´o, J. 2006, A&A, 458, 597Gr¨afener, G., Koesterke, L., & Hamann, W.-R. 2002, A&A, 387, 244Hauschildt, P. H. 1992, J. Quant. Spec. Radiat. Transf., 47, 433Hauschildt, P. H. & Baron, E. 2006, A&A, 451, 273Hauschildt, P. H., St¨orzer, H., & Baron, E. 1994, J. Quant. Spec. Radiat. Transf.,51, 875Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003,ApJ, 591, 288Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140Hillier, D. J. & Miller, D. L. 1998, ApJ, 496, 407Holzreuter, R. & Solanki, S. K. 2012, A&A, 547, A46Hubeny, I. & Mihalas, D. 2014, Theory of Stellar AtmospheresIbgui, L., Hubeny, I., Lanz, T., & Stehl´e, C. 2013, A&A, 549, A126Jones, H. P. 1973, ApJ, 185, 183Jones, H. P. & Skumanich, A. 1973, ApJ, 185, 167Keszthelyi, Z., Wade, G. A., & Petit, V. 2017, in IAU Symposium, Vol. 329,The Lives and Death-Throes of Massive Stars, ed. J. J. Eldridge, J. C. Bray,L. A. S. McClelland, & L. Xiao, 250–254Kunasz, P. & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67Kunasz, P. B. & Olson, G. L. 1988, J. Quant. Spec. Radiat. Transf., 39, 1Lamers, H. J. G. L. M., Cerruti-Sola, M., & Perinotto, M. 1987, ApJ, 314, 726Leenaarts, J. & Carlsson, M. 2009, in Astronomical Society of the PacificConference Series, Vol. 415, The Second Hinode Science Meeting: BeyondDiscovery-Toward Understanding, ed. B. Lites, M. Cheung, T. Magara,J. Mariska, & K. Reeves, 87Lobel, A. & Blomme, R. 2008, ApJ, 678, 408Maeder, A. 1999, A&A, 347, 185Maeder, A. & Meynet, G. 2000, A&A, 361, 159Marcolino, W. L. F., Bouret, J.-C., Sundqvist, J. O., et al. 2013, MNRAS, 431,2253Mason, B. D., Hartkopf, W. I., Gies, D. R., Henry, T. J., & Helsel, J. W. 2009,AJ, 137, 3358Mihalas, D. 1978, Stellar atmospheres (2nd edition) (San Francisco:W. H. Freeman and Co., 1978)Ng, K.-C. 1974, J. Chem. Phys., 61, 2680Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf.,35, 431Olson, G. L. & Kunasz, P. B. 1987, J. Quant. Spec. Radiat. Transf., 38, 325Owocki, S. P., Cranmer, S. R., & Gayley, K. G. 1996, ApJ, 472, L115Owocki, S. P., Gayley, K. G., & Cranmer, S. R. 1998, in Astronomical Societyof the Pacific Conference Series, Vol. 131, Properties of Hot Luminous Stars,ed. I. Howarth, 237Pauldrach, A., Puls, J., & Kudritzki, R. P. 1986, A&A, 164, 86Pauldrach, A. W. A., Ho ff mann, T. L., & Lennon, M. 2001, A&A, 375, 161Petit, V., Keszthelyi, Z., MacInnis, R., et al. 2017, MNRAS, 466, 1052Petit, V., Owocki, S. P., Wade, G. A., et al. 2013, MNRAS, 429, 398Petrenz, P. & Puls, J. 1996, A&A, 312, 195Petrenz, P. & Puls, J. 2000, A&A, 358, 956Petrov, B., Vink, J. S., & Gr¨afener, G. 2016, MNRAS, 458, 1999Prilutskii, O. F. & Usov, V. V. 1976, Soviet Ast., 20, 2Puls, J. 1991, A&A, 248, 581Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, Journal of ComputationalPhysics, 318, 277Rivero Gonz´alez, J. G., Puls, J., Najarro, F., & Brott, I. 2012, A&A, 537, A79Rybicki, G. B. & Hummer, D. G. 1978, ApJ, 219, 654Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107Schwarz, H. R. 1997, Numerische Mathematik (B. G. Teubner, Stuttgart 1997)Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265Sundqvist, J. O., ud-Doula, A., Owocki, S. P., et al. 2012, MNRAS, 423, L21ud-Doula, A. & Owocki, S. P. 2002, ApJ, 576, 413ud-Doula, A., Owocki, S. P., & Townsend, R. H. D. 2008, MNRAS, 385, 97van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066Vanbeveren, D. 1991, A&A, 252, 159Vath, H. M. 1994, A&A, 284, 319
18. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. A.1.
Probability density functions of radial (dashed) and x (solid) coordinates for di ff erent spherical and Cartesian grids.In this example, two spherical grids are given in 2D as input toour 3D code, with uniformly (black) or logarithmically (red) dis-tributed r -coordinates, and a uniformly distributed polar angle.The corresponding distributions of x -coordinates are calculatedwithin our grid construction procedure (see text). Large values ofthe probability density functions correspond to a high resolutionof x and r -coordinates. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 1999, A&A, 350, 181von Zeipel, H. 1924, MNRAS, 84, 665Walborn, N. R., Sana, H., Taylor, W. D., Sim´on-D´ıaz, S., & Evans, C. J. 2012, inAstronomical Society of the Pacific Conference Series, Vol. 465, Proceedingsof a Scientific Meeting in Honor of Anthony F. J. Mo ff at, ed. L. Drissen,C. Robert, N. St-Louis, & A. F. J. Mo ff at, 490Zsarg´o, J., Hillier, D. J., & Georgiev, L. N. 2006, A&A, 447, 1093 Appendix A: Grid construction
In this section, we describe the grid construction procedure usedwithin our 3D code. Generally, we assume the wind structure(i.e. density, velocity field, and temperature) to be given by an in-put model obtained from hydrodynamic simulations or external(semi)-analytic calculations. Since the input grid is not necessar-ily compatible with our 3D SC solver, and to minimize interpo-lation errors when calculating upwind and downwind quantities,we construct an own grid that uses the distribution of the input-grid coordinates in an optimum way. When the input grid usesspherical coordinates ( r , Θ , Φ ), we define a joint probability dis-tribution h ( x , z ) = f ( r ) g ( Θ ) | J | , (A.1)where f ( r ) and g ( Θ ) are the probability density functions derivedfrom the distribution of the input coordinates, and | J | = √ x + z is the Jacobian determinant. Since we consider only axisym-metric atmospheres in this paper, we use the x -coordinate dis-tribution also for the y-coordinates. To calculate the probabil-ity density functions for x and z , we simply marginalize h ( x , z )over z and x , respectively. The discretized coordinates are finallydetermined by demanding that the probabilities of selecting a(continuous) coordinate in each (discrete) interval shall be thesame. Fig. A.1 shows the probability density functions of the x -coordinates for two di ff erent input distributions of the radialgrid. Here, the polar angle Θ has been assumed to be uniformlydistributed for both examples. Fig. B.1.
B´ezier curves (solid lines) for three given points b , b , b . The blue, red, and green lines represent the resulting curvesfor di ff erent control points b . The straight connections of thecontrol points with the data points are indicated by the dottedlines.Since the final number of core and non-core points dependson the slope of the probability density of the radial grid, yieldingin worst cases a much larger number of core points than non-corepoints, and because the total number of used points is memory-limited, we define two input parameters N core and N non − core tokeep control on the final grid. For all test calculations (includingthose presented in Sect. 4.2), the best solution has always beenfound for a number of N core / N non − core ∈ [0 . , . N core and N non − core corresponds to a re-normalization ofthe probability density function in the regimes x , z ∈ [0 , R ∗ ] and x , z ∈ [ R ∗ , R max ], where R max defines the border of the computa-tional domain. We note that the same procedure can be used foran input grid given in Cartesian coordinates, with the probabil-ity density function of the input-grid coordinates derived directlyfrom the corresponding (discrete) distribution. Appendix B: 1D B ´ezier interpolation
In this section, we discuss an interpolation technique usingquadratic B´ezier curves (e.g. Auer 2003, Schwarz 1997). Suchcurves are generally constructed from three given points b , b , b (see Fig. B.1), and have the following useful properties:(i) The boundary points, b and b are reproduced exactly bythe B´ezier curve.(ii) The straight connections ( b − b ) and ( b − b ) define thetangent lines of the curve at b and b , respectively.(iii) Any point on the B´ezier curve is located in the convex hullof b , b , b .In a 2D plane described by coordinates x and y , the quadraticB´ezier curve is parameterized as: b ( t ) = x ( t ) y ( t ) ! = (1 − t ) b + t (1 − t ) b + t b , (B.1)with t ∈ [0 , b = ( x , y ), b = ( x , y ), b = ( x , y ).With Eq. (B.1), the properties (i)-(iii) can be exploited to con-struct a monotonic interpolation scheme by identifying b , b with two given data points ( x , f ), ( x , f ), and defining b asa free (and tunable) parameter. Thus, b is commonly named
19. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
Fig. B.2. Di ff erent interpolation techniques for a set of three datapoints at x -coordinates indicated by the dotted vertical lines. Thesolid and dashed lines correspond to the interpolation in the dif-ferent intervals [ x i − , x i ] and [ x i , x i + ], respectively. Linear in-terpolations, quadratic interpolations (connecting all three datapoints), and a monotonic B´ezier curve (with ω calculated fromEq. (B.5) in the interval [ x i − , x i ]) are indicated in red, blue, andgreen. Since the quadratic interpolation is already monotonic inthe interval [ x i , x i + ], the monotonic B´ezier curve coincides withthe dashed, blue line. Control points are indicated with colouredasterisks.control point, and is ‘only’ required to set the slope of the B´eziercurve. To reproduce the underlying function best, and to preservemonotonicity of the resulting curve, the control point should bechosen with care.In the following, we present a B´ezier-interpolation techniquefor an interval x ∈ [ x i − , x i ], given three data points, ( x i − , f i − ),( x i , f i ), ( x i + , f i + ). The interpolation formulas corresponding tothe interval x ∈ [ x i , x i + ] are given in subsection B.2. B.1. Interval [ x i − , x i ]A quadratic B´ezier curve in the interval [ x i − , x i ] is given fromEq. (B.1): x ( t ) f ( x ( t )) ! = (1 − t ) x i − f i − ! + t (1 − t ) x c f c ! + t x i f i ! , (B.2)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 asecond-order interpolation scheme, however, x c needs to be lo-cated at the centre of the data-point’s abscissae , and is thereforeset to x c = ( x i − + x i ) /
2. Then, the quadratic B´ezier interpolationscheme is given by: f ( x ) = (1 − t ) f i − + t (1 − t ) f c + t f i (B.3) t = ( x − x i − ) / ( x i − x i − ) , where t has been determined from the definition of x c andEq. (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 If x c was located at x c = x i − + / x i − x i − ), for instance, one caneasily show that the resulting B´ezier curve never reproduces the unitparabola for any ordinate value f c . B´ezier curve at this data point, f c is calculated as f c = f i − d f d x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x i ∆ x i , with ∆ x i = x i − x i − . The unknown derivative at x i needs tobe approximated. Using also the information from the next datapoint, ( x i + , f i + ), and assigning a weight ω to the forward andbackward derivatives (obtained from finite di ff erences), we find f c = f i − ∆ x i ω f i − f i − ∆ x i + (1 − ω ) f i + − f i ∆ x i + ! , (B.4)with ∆ x i + = x i + − x i . With a proper choice of ω , we can adjustthe B´ezier curve to our needs by shifting the control point upor down. For instance, setting ω = ∆ x i + / ( ∆ x i + ∆ x i + ) resultsin the unique parabola connecting the three given data points,while ω = x i − , x i ]. Noting thatmonotonicity is obtained when the control point is located in theinterval f c ∈ [ f i − , f i ], corresponding ω -values should lie in be-tween the following limits: ω [ i − , i ] i − : = ω ( f [ i − , i ]c = f i − ) = + − f i + − f i f i − f i − ∆ x i ∆ x i + (B.5) ω [ i − , i ] i : = ω ( f [ i − , i ]c = f i ) = − f i − f i − f i + − f i ∆ x i + ∆ x i , (B.6)where the superscript [ i − , i ] denotes that ω corresponds tothe interpolation scheme in the left interval, [ x i − , x i ]. In the fi-nal implementation, we avoid the division by zero if f i = f i − or f i = f i + , of course. Our standard interpolation is then per-formed as follows. At first, we calculate ω such that we obtainthe unique parabola connecting all three data points. Secondly, if ω lies outside the allowed limits from Eq. (B.5) and (B.6), we ad-just ω to yield monotonic interpolations. In Fig. B.2, we displaythe monotonic B´ezier curve resulting from a ω -parameter calcu-lated by means of Eq. (B.5), together with linear and quadraticinterpolations (the latter connecting the three data points). Sincemonotonicity is always obtained for ω ∈ [ ω i − , ω i ], we can de-fine even stricter limits in order to avoid oscillations during theiteration scheme, by setting ω = Λ -matrix, theinterpolation coe ffi cients are required. Combining Eqs. (B.3) and(B.4) then gives: f (cid:0) x ∈ [ x i − , x i ] (cid:1) = ˜ a [ i − , i ] f i − + ˜ b [ i − , i ] f i + ˜ c [ i − , i ] f i + , (B.7)with˜ a [ i − , i ] = + ( ω − x − x i − x i − x i − + (1 − ω ) x − x i − x i − x i − ! (B.8)˜ b [ i − , i ] = (1 − ω ) ∆ x i + (2 − ω ) ∆ x i + ∆ x i + x − x i − x i − x i − + ( ω − ∆ x i + ∆ x i + ∆ x i + x − x i − x i − x i − ! (B.9)˜ c [ i − , i ] = ( ω − ∆ x i ∆ x i + x − x i − x i − x i − − ( ω − ∆ x i ∆ x i + x − x i − x i − x i − ! . (B.10)
20. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
B.2. Interval [ x i , x i + ]The interpolation formula for the right interval [ x i , x i + ] uses thesame data points as above. Since the value of the control pointneeds to be calculated at a di ff erent x-coordinate, x c = ( x i + + x i ) /
2, we cannot simply substitute indices. Using f ( x ) = (1 − t ) f i + t (1 − t ) f c + t f i + (B.11) t : = ( x − x i ) / ( x i + − x i ) f c = f i + ∆ x i + ω f i + − f i ∆ x i + + (1 − ω ) f i − f i − ∆ x i ! , (B.12)we obtain for this interval: f ( x ∈ [ x i , x i + ]) = ˜ a [ i , i + f i − + ˜ b [ i , i + f i + ˜ c [ i , i + f i + , (B.13)with˜ a [ i , i + = ( ω − ∆ x i + ∆ x i x − x i x i + − x i − ( ω − ∆ x i + ∆ x i x − x i x i + − x i ! (B.14)˜ b [ i , i + = − ω ∆ x i + ( ω − ∆ x i + ∆ x i x − x i x i + − x i + ( ω − ∆ x i + ∆ x i + ∆ x i x − x i x i + − x i ! (B.15)˜ c [ i , i + = ω x − x i x i + − x i + (1 − ω ) x − x i x i + − x i , ! (B.16)and ω [ i , i + i : = ω ( f [ i , i + = f i ) = − f i + − f i f i − f i − ∆ x i ∆ x i + (B.17) ω [ i , i + i + : = ω ( f [ i , i + = f i + ) = + − f i − f i − f i + − f i ∆ x i + ∆ x i . (B.18)The corresponding B´ezier curves for di ff erent ω -parameters( ω = ω = ∆ x i / ( ∆ x i + ∆ x i + ) for continuousquadratic interpolations) are also shown in Fig. B.2. We notethat the B´ezier interpolation gives a continuous function in thecomplete interval [ x i − , x i + ] only for those ω -values that definethe parabola connecting all three data points. Appendix C: 2D B ´ezier interpolation
To interpolate upwind and downwind quantities, a 2D interpola-tion scheme is required. Fig. C.1 displays the geometry for a 2Drectangular area, with grid points indicated by the black dots.With this setup, we perform a 2D B´ezier interpolation by ap-plying three 1D B´ezier interpolations along the x -axis on each y -level at ( j − j ), ( j + y at the desired x -coordinate. Within the cyanshaded interval, we obtain with the 1D B´ezier interpolation givenby Eqs. (B.13)-(B.16): f ( x , y ) = ˜ a y ˜ a ( j − x f i − , j − + ˜ a y ˜ b ( j − x f i , j − + ˜ a y ˜ c ( j − x f i + , j − + ˜ b y ˜ a ( j ) x f i − , j + ˜ b y ˜ b ( j ) x f i , j + ˜ b y ˜ c ( j ) x f i + , j + ˜ c y ˜ a ( j + x f i − , j + + ˜ c y ˜ b ( j + x f i , j + + ˜ c y ˜ c ( j + x f i + , j + , (C.1)where the subscripts of the interpolation coe ffi cients indicate thecoordinate used for each 1D interpolation. We note that all up-wind and downwind interpolations are performed in the upperright interval of a given surface, in order to obtain a simple rep-resentation of the Λ -matrix elements. PSfrag replacements i − , j − i , j − i + , j − i − , j i , j i + , ji − , j + i , j + i + , j + x , j + x , jx , j − x , y ) Fig. C.1.
2D interpolation for upwind or downwind quantitiesrequired in the cyan shaded area. The 2D B´ezier interpolationconsists of three 1D interpolations to obtain the values at thedesired x -coordinate (indicated by red dots), followed by a 1Dinterpolation along the y-coordinate using the obtained values atthe red dots. Appendix D: ALO coefficients
In this section, we derive the Λ -matrix coe ffi cients used to con-struct the approximate Λ -operator. We note that the obtained ma-trix elements can also be used for any other (2nd or lower order)interpolation scheme using the same geometry, with di ff erent in-terpolation coe ffi cients though.For a source function set to unity at grid point ( i jk ) andzero everywhere else, we consider all 27 points ranging from( i − α, j − β, k − γ ) to ( i + α, j + β, k + γ ). The correspondingmatrix coe ffi cients are derived from Eq. (29), using the dis-cretized equation of radiative transfer, Eqs. (12) / (20), with up-wind and downwind interpolations defined by Eqs. (22) and (23).We further consider only the Λ Ω ν -operator, since the integrationover frequency and / or solid angle is straightforward. Each Λ -matrix element then corresponds to the intensity (resulting from S i jk =
1) at a considered grid point p (not necessarily identical to( i jk )), and consists of an emission term (defined by the interpo-lated source functions and optical-depth steps at the correspond-ing upwind, current, and downwind points), and the irradiationfrom the upwind point (defined by the upwind intensity and up-wind optical-depth step). The upper and lower panel of Fig. D.1show an example in 2D considering the points ( i − α, j − β ) and( i , j − β ) for a source function S i j = Ω , ν . The m , n -th Λ -element is written as Λ nm , with matrix indices n , m cal-culated from Eq. (30). While n corresponds to the 3D indices ofthe local grid point ( S i jk = S n = m represents the neighbour-ing point, ( i − α, j − β, k − γ ). Applying Eq. (29) to the specificintensity at point m , we obtain: Λ nm = I m ( S = e n , Φ B = = Λ i jki − α, j − β, k − γ = I i − α, j − β, k − γ (cid:16) S i jk = δ ˜ i , i δ ˜ j , j δ ˜ k , k (cid:17) = a i − α, j − β, k − γ S ( i − α, j − β, k − γ )u (cid:16) S i jk = δ ˜ i , i δ ˜ j , j δ ˜ k , k (cid:17) + b i − α, j − β, k − γ S ( i − α, j − β, k − γ )p (cid:16) S i jk = δ ˜ i , i δ ˜ j , j δ ˜ k , k (cid:17) + c i − α, j − β, k − γ S ( i − α, j − β, k − γ )d (cid:16) S i jk = δ ˜ i , i δ ˜ j , j δ ˜ k , k (cid:17) + d i − α, j − β, k − γ I ( i − α, j − β, k − γ )u (cid:16) S i jk = δ ˜ i , i δ ˜ j , j δ ˜ k , k (cid:17) ,
21. Hennicker et al.: 3D short-characteristics method in the winds from OB stars
PSfrag replacements I i − α, j − β = I i − α, j − β = I i − α, j − β u = S i − α, j − β u = Λ iji − α, j − β = I i − α, j − β S i − α, j − β p = S i − α, j − β = S i − α, j − β d , S i , j = I i , j − β u , S i , j − β u = Λ iji , j − β = I i , j − β S i , j − β p = S i , j − β = S i − α, j − β d , S i , j = I i − α, j − β = Λ iji − α, j − β I i − α, j − β = n (see upperpanel) PSfrag replacements I i − α, j − β = I i − α, j − β = I i − α, j − β u = S i − α, j − β u = Λ iji − α, j − β = I i − α, j − β S i − α, j − β p = S i − α, j − β = S i − α, j − β d , S i , j = I i , j − β u , S i , j − β u = Λ iji , j − β = I i , j − β S i , j − β p = S i , j − β = S i − α, j − β d , S i , j = I i − α, j − β = Λ iji − α, j − β I i − α, j − β = n (see upperpanel) Fig. D.1.
2D example for calculating the Λ Ω ,ν -matrix elementsat a grid point ( i − α, j − β ) (upper panel) and ( i , j − β ) (lowerpanel). The matrix elements correspond to the intensity at theconsidered grid points calculated for a source function S i j = i − α, j − β ).with boundary contribution Φ B , n -th unit vector e n , and δ ˜ i , i , δ ˜ j , j , δ ˜ k , k the Kronecker- δ for all possible x ˜ i , y ˜ j , and z ˜ k coordinates,respectively. S u and S d are the upwind and downwind sourcefunctions corresponding to a considered short characteristic atgrid point p ↔ ( i − α, j − β, k − γ ), S p is the source functionat the grid point , I u is the upwind intensity, and a , b , c , d arethe integration coe ffi cients for this particular short characteris-tic. All upwind and downwind quantities are to be interpolatedfrom neighbouring grid points. We use the notation w , ˆ w , ˜ w , toidentify di ff erent interpolation coe ffi cients corresponding to theupwind source function, upwind intensity, and downwind source S p , ↔ ( i jk ). Then, S ( ijk )d =
0, and S ( ijk )u , function, respectively. Using Eqs. (22) and (23) to interpolateupwind and downwind quantities, we find: Λ i jki − α, j − β, k − γ = a i − α, j − β, k − γ · (cid:2) w A S i − α, j − β, k − γ + w B S i − α, j − β, k − γ + w C S i − α, j − β, k − γ + w D S i − α, j − β, k − γ + w E S i − α, j − β, k − γ + w F S i − α, j − β, k − γ + w G S i − α, j − β, k − γ + w H S i − α, j − β, k − γ + w I S i − α, j − β, k − γ + w J S i − α, j − β, k − γ + w K S i − α, j − β, k − γ + w L S i − α, j − β, k − γ + w M S i − α, j − β, k − γ + w N S i − α, j − β, k − γ + w O S i − α, j − β, k − γ + w P S i − α, j − β, k − γ + w Q S i − α, j − β, k − γ + w R S i − α, j − β, k − γ + w S S i − α, j − β, k − γ + w i − α, j − β, k − γ S i − α, j − β, k − γ (cid:3) + b i − α, j − β, k − γ S i − α, j − β, k − γ + c i − α, j − β, k − γ · (cid:2) ˜ w A S i − α, j , k − γ + ˜ w B S i − α, j , k − γ + ˜ w C S i , j , k − γ + ˜ w D S i − α, j , k − γ + ˜ w E S i − α, j , k − γ + ˜ w F S i , j , k − γ + ˜ w G S i − α, j , k + ˜ w H S i − α, j , k + ˜ w I S i , j , k + ˜ w J S i − α, j − β, k + ˜ w K S i − α, j − β, k + ˜ w L S i , j − β, k + ˜ w M S i − α, j − β, k + ˜ w N S i − α, j − β, k + ˜ w O S i , j − β, k + ˜ w P S i , j − β, k − γ + ˜ w Q S i , j − β, k − γ + ˜ w R S i , j − β, k − γ + ˜ w S S i , j − β, k − γ (cid:3) + d i − α, j − β, k − γ (cid:2) ˆ w A I i − α, j − β, k − γ ( S i jk = + · · · + ˆ w S I i − α, j − β, k − γ ( S i jk = (cid:3) , with the upwind intensity interpolated from the same points asthe upwind source function, and a compact notation for the in-terpolation coe ffi cients (with skipped superscripts). Since only S i jk = Λ i jki − α, j − β, k − γ = c i − α, j − β, k − γ ˜ w ( i − α, j − β, k − γ )I . (D.1)The matrix element for a point ( i − α, j − β, k − γ ) with a non-vanishing source function at point ( i jk ) is thus solely given bythe integration coe ffi cient c i − α, j − β, k − γ from the discretized equa-tion of radiative transfer multiplied with the interpolation co-e ffi cient for the downwind source function of point I (corre-sponding to grid point ( i jk ), see Fig. 1). The other neighboursare obtained analogously, without vanishing incident intensities,however. Accounting also for the interpolation of upwind sourcefunctions and intensities when necessary, we find: Λ i jki , j − β, k − γ = c i , j − β, k − γ ˜ w i , j − β, k − γ H + d i , j − β, k − γ ˆ w i , j − β, k − γ S Λ i jki − α, j − β, k − γ (D.2) Λ i jki + α, j − β, k − γ = c i + α, j − β, k − γ ˜ w i + α, j − β, k − γ G + d i + α, j − β, k − γ ˆ w i + α, j − β, k − γ S Λ i jki , j − β, k − γ (D.3) Λ i jki − α, j , k − γ = c i − α, j , k − γ ˜ w i − α, j , k − γ O + d i − α, j , k − γ ˆ w i − α, j , k − γ I Λ i jki − α, j − β, k − γ (D.4)
22. Hennicker et al.: 3D short-characteristics method in the winds from OB stars Λ i jki , j , k − γ = c i , j , k − γ ˜ w i , j , k − γ N + d i , j , k − γ · h ˆ w i , j , k − γ H Λ i jki − α, j − β, k − γ + ˆ w i , j , k − γ I Λ i jki , j − β, k − γ + ˆ w i , j , k − γ S Λ i jki − α, j , k − γ i (D.5) Λ i jki + α, j , k − γ = c i + α, j , k − γ ˜ w i + α, j , k − γ M + d i + α, j , k − γ · h ˆ w i + α, j , k − γ D Λ i jki − α, j − β, k − γ + ˆ w i + α, j , k − γ H Λ i jki , j − β, k − γ + ˆ w i + α, j , k − γ I Λ i jki + α, j − β, k − γ + ˆ w i + α, j , k − γ S Λ i jki , j , k − γ i (D.6) Λ i jki − α, j + β, k − γ = c i − α, j + β, k − γ ˜ w i − α, j + β, k − γ L + d i − α, j + β, k − γ ˆ w i − α, j + β, k − γ I Λ i jki − α, j , k − γ (D.7) Λ i jki , j + β, k − γ = c i , j + β, k − γ ˜ w i , j + β, k − γ K + d i , j + β, k − γ · h ˆ w i , j + β, k − γ R Λ i jki − α, j − β, k − γ + ˆ w i , j + β, k − γ H Λ i jki − α, j , k − γ + ˆ w i , j + β, k − γ I Λ i jki , j , k − γ + ˆ w i , j + β, k − γ S Λ i jki − α, j + β, k − γ i (D.8) Λ i jki + α, j + β, k − γ = c i + α, j + β, k − γ ˜ w i + α, j + β, k − γ J + d i + α, j + β, k − γ · h ˆ w i + α, j + β, k − γ R Λ i jki , j − β, k − γ + ˆ w i + α, j + β, k − γ G Λ i jki − α, j , k − γ + ˆ w i + α, j + β, k − γ H Λ i jki , j , k − γ + ˆ w i + α, j + β, k − γ I Λ i jki + α, j , k − γ + ˆ w i + α, j + β, k − γ S Λ i jki , j + β, k − γ i (D.9) Λ i jki − α, j − β, k = c i − α, j − β, k ˜ w i − α, j − β, k F + d i − α, j − β, k ˆ w i − α, j − β, k O Λ i jki − α, j − β, k − γ (D.10) Λ i jki , j − β, k = c i , j − β, k ˜ w i , j − β, k E + d i , j − β, k · h ˆ w i , j − β, k N Λ i jki − α, j − β, k − γ + ˆ w i , j − β, k O Λ i jki , j − β, k − γ + ˆ w i , j − β, k S Λ i jki − α, j − β, k i (D.11) Λ i jki + α, j − β, k = c i + α, j − β, k ˜ w i + α, j − β, k D + d i + α, j − β, k · h ˆ w i + α, j − β, k M Λ i jki − α, j − β, k − γ + ˆ w i + α, j − β, k N Λ i jki , j − β, k − γ + ˆ w i + α, j − β, k O Λ i jki + α, j − β, k − γ + ˆ w i + α, j − β, k S Λ i jki , j − β, k i (D.12) Λ i jki − α, j , k = c i − α, j , k ˜ w i − α, j , k S + d i − α, j , k · h ˆ w i − α, j , k F Λ i jki − α, j − β, k − γ + ˆ w i − α, j , k O Λ i jki − α, j , k − γ + ˆ w i − α, j , k I Λ i jki − α, j − β, k i (D.13) Λ i jki jk = a i jk w i jk + b i jk + d i jk · h ˆ w i jk E Λ i jki − α, j − β, k − γ + ˆ w i jk F Λ i jki , j − β, k − γ + ˆ w i jk N Λ i jki − α, j , k − γ + ˆ w i jk O Λ i jki , j , k − γ + ˆ w i jk H Λ i jki − α, j − β, k + ˆ w i jk I Λ i jki , j − β, k + ˆ w i jk S Λ i jki − α, j , k i (D.14) Λ i jki + α, j , k = a i + α, j , k w i + α, j , k S + d i + α, j , k · h ˆ w i + α, j , k D Λ i jki − α, j − β, k − γ + ˆ w i + α, j , k E Λ i jki , j − β, k − γ + ˆ w i + α, j , k F Λ i jki + α, j − β, k − γ + ˆ w i + α, j , k M Λ i jki − α, j , k − γ + ˆ w i + α, j , k N Λ i jki , j , k − γ + ˆ w i + α, j , k O Λ i jki + α, j , k − γ + ˆ w i + α, j , k G Λ i jki − α, j − β, k + ˆ w i + α, j , k H Λ i jki , j − β, k + ˆ w i + α, j , k I Λ i jki + α, j − β, k + ˆ w i + α, j , k S Λ i jki jk i (D.15) Λ i jki − α, j + β, k = c i − α, j + β, k ˜ w i − α, j + β, k R + d i − α, j + β, k · h ˆ w i − α, j + β, k L Λ i jki − α, j − β, k − γ + ˆ w i − α, j + β, k F Λ i jki − α, j , k − γ + ˆ w i − α, j + β, k O Λ i jki − α, j + β, k − γ + ˆ w i − α, j + β, k I Λ i jki − α, j , k i (D.16) Λ i jki , j + β, k = a i , j + β, k w i , j + β, k I + d i , j + β, k · h ˆ w i , j + β, k E Λ i jki − α, j , k − γ + ˆ w i , j + β, k F Λ i jki , j , k − γ + ˆ w i , j + β, k H Λ i jki − α, j , k + ˆ w i , j + β, k I Λ i jki jk + ˆ w i , j + β, k K Λ i jki − α, j − β, k − γ + ˆ w i , j + β, k L Λ i jki , j − β, k − γ + ˆ w i , j + β, k N Λ i jki − α, j + β, k − γ + ˆ w i , j + β, k O Λ i jki , j + β, k − γ + ˆ w i , j + β, k R Λ i jki − α, j − β, k + ˆ w i , j + β, k S Λ i jki − α, j + β, k + i (D.17) Λ i jki + α, j + β, k = a i + α, j + β, k w i + α, j + β, k H + d i + α, j + β, k · h ˆ w i + α, j + β, k D Λ i jki − α, j , k − γ + ˆ w i + α, j + β, k E Λ i jki , j , k − γ + ˆ w i + α, j + β, k F Λ i jki + α, j , k − γ + ˆ w i + α, j + β, k G Λ i jki − α, j , k + ˆ w i + α, j + β, k H Λ i jki jk + ˆ w i + α, j + β, k I Λ i jki + alpha , j , k + ˆ w i + α, j + β, k J Λ i jki − α, j − β, k − γ + ˆ w i + α, j + β, k K Λ i jki , j − β, k − γ + ˆ w i + α, j + β, k L Λ i jki + α, j − β, k − γ + ˆ w i + α, j + β, k M Λ i jki − α, j + β, k − γ + ˆ w i + α, j + β, k O Λ i jki + α, j + β, k − γ + ˆ w i + α, j + β, k R Λ i jki , j − β, k + ˆ w i + α, j + β, k S Λ i jki , j + β, k i (D.18) Λ i jki − α, j − β, k + γ = c i − α, j − β, k + γ ˜ w i − α, j − β, k + γ C + d i − α, j − β, k + γ ˆ w i − α, j − β, k + γ O Λ i jki − α, j − β, k (D.19) Λ i jki , j − β, k + γ = c i , j − β, k + γ ˜ w i , j − β, k + γ B + d i , j − β, k + γ · h ˆ w i , j − β, k + γ N Λ i jki − α, j − β, k + ˆ w i , j − β, k + γ O Λ i jki , j − β, k + ˆ w i , j − β, k + γ Q Λ i jki − α, j − β, k − γ + ˆ w i , j − β, k + γ S Λ i jki − α, j − β, k + γ i (D.20) Λ i jki + α, j − β, k + γ = c i + α, j − β, k + γ ˜ w i + α, j − β, k + γ A + d i + α, j − β, k + γ · h ˆ w i + α, j − β, k + γ M Λ i jki − α, j − β, k + ˆ w i + α, j − β, k + γ N Λ i jki , j − β, k + ˆ w i + α, j − β, k + γ O Λ i jki + α, j − β, k + ˆ w i + α, j − β, k + γ Q Λ i jki , j − β, k − γ + ˆ w i + α, j − β, k + γ S Λ i jki , j − β, k + γ i (D.21) Λ i jki − α, j , k + γ = c i − α, j , k + γ ˜ w i − α, j , k + γ Q + d i − α, j , k + γ · h ˆ w i − α, j , k + γ C Λ i jki − α, j − β, k − γ + ˆ w i − α, j , k + γ F Λ i jki − α, j − β, k + ˆ w i − α, j , k + γ I Λ i jki − α, j − β, k + γ + ˆ w i − α, j , k + γ O Λ i jki − α, j , k i (D.22) Λ i jki , j , k + γ = a i , j , k + γ w i , j , k + γ O + d i , j , k + γ · h ˆ w i , j , k + γ B Λ i jki − α, j − β, k − γ + ˆ w i , j , k + γ C Λ i jki , j − β, k − γ + ˆ w i , j , k + γ E Λ i jki − α, j − β, k + ˆ w i , j , k + γ F Λ i jki , j − β, k + ˆ w i , j , k + γ H Λ i jki − α, j − β, k + γ + ˆ w i , j , k + γ I Λ i jki , j − β, k + γ + ˆ w i , j , k + γ N Λ i jki − α, j , k + ˆ w i , j , k + γ O Λ i jki jk + ˆ w i , j , k + γ Q Λ i jki − α, j , k − γ + ˆ w i , j , k + γ S Λ i jki − α, j , k + γ i (D.23)
23. Hennicker et al.: 3D short-characteristics method in the winds from OB stars Λ i jki + α, j , k + γ = a i + α, j , k + γ w i + α, j , k + γ N + d i + α, j , k + γ · h ˆ w i + α, j , k + γ A Λ i jki − α, j − β, k − γ + ˆ w i + α, j , k + γ B Λ i jki , j − β, k − γ + ˆ w i + α, j , k + γ C Λ i jki + α, j − β, k − γ + ˆ w i + α, j , k + γ D Λ i jki − α, j − β, k + ˆ w i + α, j , k + γ E Λ i jki , j − β, k + ˆ w i + α, j , k + γ F Λ i jki + alpha , j , k + ˆ w i + α, j , k + γ G Λ i jki − α, j − β, k + γ + ˆ w i + α, j , k + γ H Λ i jki , j − β, k + γ + ˆ w i + α, j , k + γ I Λ i jki + α, j − β, k + γ + ˆ w i + α, j , k + γ M Λ i jki − α, j , k + ˆ w i + α, j , k + γ N Λ i jki jk + ˆ w i + α, j , k + γ O Λ i jki + alpha , j , k + ˆ w i + α, j , k + γ Q Λ i jki , j , k − γ + ˆ w i + α, j , k + γ S Λ i jki , j , k + γ i (D.24) Λ i jki − α, j + β, k + γ = c i − α, j + β, k + γ ˜ w i − α, j + β, k + γ P + d i − α, j + β, k + γ · h ˆ w i − α, j + β, k + γ C Λ i jki − α, j , k − γ + ˆ w i − α, j + β, k + γ F Λ i jki − α, j , k + ˆ w i − α, j + β, k + γ I Λ i jki − α, j , k + γ + ˆ w i − α, j + β, k + γ L Λ i jki − α, j − β, k + ˆ w i − α, j + β, k + γ O Λ i jki − α, j + β, k i (D.25) Λ i jki , j + β, k + γ = a i , j + β, k + γ w i , j + β, k + γ F + d i , j + β, k + γ · h ˆ w i , j + β, k + γ B Λ i jki − α, j , k − γ + ˆ w i , j + β, k + γ C Λ i jki , j , k − γ + ˆ w i , j + β, k + γ E Λ i jki − α, j , k + ˆ w i , j + β, k + γ F Λ i jki jk + ˆ w i , j + β, k + γ H Λ i jki − α, j , k + γ + ˆ w i , j + β, k + γ I Λ i jki , j , k + γ + ˆ w i , j + β, k + γ K Λ i jki − α, j − β, k + ˆ w i , j + β, k + γ L Λ i jki , j − β, k + ˆ w i , j + β, k + γ N Λ i jki − α, j + β, k + ˆ w i , j + β, k + γ O Λ i jki , j + β, k + ˆ w i , j + β, k + γ P Λ i jki − α, j − β, k − γ + ˆ w i , j + β, k + γ Q Λ i jki − α, j + β, k − γ + ˆ w i , j + β, k + γ R Λ i jki − α, j − β, k + γ + ˆ w i , j + β, k + γ S Λ i jki − α, j + β, k − γ i (D.26) Λ i jki + α, j + β, k + γ = a i + α, j + β, k + γ w i + α, j + β, k + γ E + d i + α, j + β, k + γ · h ˆ w i + α, j + β, k + γ A Λ i jki − α, j , k − γ + ˆ w i + α, j + β, k + γ B Λ i jki , j , k − γ + ˆ w i + α, j + β, k + γ C Λ i jki + α, j , k − γ + ˆ w i + α, j + β, k + γ D Λ i jki − α, j , k + ˆ w i + α, j + β, k + γ E Λ i jki jk + ˆ w i + α, j + β, k + γ F Λ i jki + alpha , j , k + ˆ w i + α, j + β, k + γ G Λ i jki − α, j , k + γ + ˆ w i + α, j + β, k + γ H Λ i jki , j , k + γ + ˆ w i + α, j + β, k + γ I Λ i jki + α, j , k + γ + ˆ w i + α, j + β, k + γ J Λ i jki − α, j − β, k + ˆ w i + α, j + β, k + γ K Λ i jki , j − β, k + ˆ w i + α, j + β, k + γ L Λ i jki + α, j − β, k + ˆ w i + α, j + β, k + γ M Λ i jki − α, j + β, k + ˆ w i + α, j + β, k + γ N Λ i jki , j + β, k + ˆ w i + α, j + β, k + γ O Λ i jki + α, j + β, k + ˆ w i + α, j + β, k + γ P Λ i jki , j − β, k − γ + ˆ w i + α, j + β, k + γ Q Λ i jki , j + β, k − γ + ˆ w i + α, j + β, k + γ R Λ i jki , j − β, k + γ + ˆ w i + α, j + β, k + γ S Λ i jki , j + β, k + γ i (D.27)Since the integration and interpolation coe ffi cients need to becalculated only once at each considered grid point (here denotedby ( u , v, w ), to avoid confusion), we obtain the Λ -matrix coe ffi -cients by substituting indices. For Eq. (D.1), we find: Λ i jki − α, j − β, k − γ = Λ u + α,v + β,w + γ u vw = c u vw ˜ w u vw I , (D.28) and proceed analogously for all other elements in Eqs. (D.2)-(D.27). Thus, the ALO can be calculated in parallel to the formalsolution scheme.(D.28) and proceed analogously for all other elements in Eqs. (D.2)-(D.27). Thus, the ALO can be calculated in parallel to the formalsolution scheme.