Issue 
A&A
Volume 616, August 2018



Article Number  A140  
Number of page(s)  21  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201731858  
Published online  11 September 2018 
3D radiative transfer: Continuum and line scattering in nonspherical winds from OB stars
^{1}
LMU München, Universitätssternwarte, Scheinerstr. 1, 81679 München, Germany
email: levin@usm.lmu.de, levin@usm.unimuenchen.de
^{2}
Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
^{3}
Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
Received:
30
August
2017
Accepted:
29
April
2018
Context. State of the art quantitative spectroscopy utilizes synthetic spectra to extract information from observations. For hot, massive stars, these synthetic spectra are calculated by means of 1D, spherically symmetric, NLTE atmosphere and spectrumsynthesis codes. Certain stellar atmospheres, however, show strong deviations from spherical symmetry, and need to be treated in 3D.
Aims. We present and test a newly developed 3D radiative transfer code, tailored to the solution of the radiation field in rapidly expanding stellar atmospheres. We apply our code to the continuum transfer in windablation models, and to the UV resonance line formation in magnetic winds.
Methods. We have used a 3D finitevolume method for the solution of the timeindependent equation of radiative transfer, to study continuum and linescattering problems, currently approximated by a twolevelatom. Convergence has been accelerated by coupling the formal solver to a nonlocal approximate Λiteration scheme. Particular emphasis has been put on careful tests, by comparing with alternative solutions for 1D, spherically symmetric model atmospheres. These tests allowed us to understand certain shortcomings of the methods, and to estimate limiting cases that can actually be calculated.
Results. The typical errors of the converged source functions, when compared to 1D solutions, are of the order of 10–20%, and rapidly increase for optically thick (τ ≳ 10) continua, mostly due to the order of accuracy of our solution scheme. In circumstellar discs, the radiation temperatures in the (optically thin) transition region from wind to disc are quite similar to corresponding values in the wind. For MHD simulations of dynamical magnetospheres, the line profiles, calculated with our new 3D code, agree well with previous solutions using a 3DSEI method. When compared with profiles resulting from the socalled analytic dynamical magnetosphere (ADM) model, however, significant differences become apparent.
Conclusions. Due to similar radiation temperatures in the wind and the transition region to the disc, the same linestrength distribution can be applied within radiation hydrodynamic calculations for optically thick circumstellar discs in “accreting highmass stars”. To properly describe the UV line formation in dynamical magnetospheres, the ADM model needs to be further developed, at least in a large part of the outer wind.
Key words: radiative transfer / methods: numerical / circumstellar matter / stars: magnetic field / stars: winds, outflows
© ESO 2018
1. Introduction
Hot, massive stars are key tools to interpret the Universe. For example, they both enrich the interstellar medium with highly processed material due to their strong stellar winds, as well as end their lives as supernovae, enriching the surroundings with metals even further. Additionally, the latter produce strong shocks within the interstellar medium, thus triggering star formation.
Massive stars are thought to be the progenitors of massive black holes (M_{BH} ≳ 25 M_{⊙}) when their final fate is a direct collapse into a black hole instead of a supernova explosion. The collapse into a black hole, however, is only possible when the progenitor stars exhibit a weak wind, or, more generally, do not lose too much mass during their lifetime. Belczynski et al. (2016) proposed a lowmetallicity environment to explain the massive black hole merger GW150914 (M_{1} ≈ 36 M_{⊙}, M_{2} ≈ 29 M_{⊙}, Abbott et al. 2016) observed at the Laser Interferometric GravitationalWave Observatory (LIGO), whereas Petit et al. (2017) showed that magnetic surface fields could also lower the mass loss to form such progenitor stars.
The physical properties of massive stars are generally derived by means of quantitative spectroscopy, that is by modelling their atmospheres numerically and comparing the calculated synthetic spectra with observations. When using onedimensional (1D), spherically symmetric codes (e.g. CMFGEN Hillier & Miller 1998, PHOENIX^{1} Hauschildt 1992, POWR Gräfener et al. 2002, WMBASIC Pauldrach et al. 2001, FASTWIND Puls et al. 2005; Rivero González et al. 2012, as a nonexhaustive list for massive star atmospheric models), this has meanwhile become a routine job and is widely applied. Within the last decade, however, it became more and more evident, from both the theoretical and observational side, that many stars differ from spherical symmetry, rendering the results from 1D codes questionable for such objects. The deviations can be of different origin and shape:
(i) Magnetic winds: Wade et al. (2012) showed within the Magnetism in Massive Stars (MiMeS) survey, that about 7% of all (Galactic) OBstars have detectable magnetic fields. Magnetohydrodynamic (MHD) calculations from udDoula & Owocki (2002) and udDoula et al. (2008) revealed, that large scale magnetic fields, possibly of fossil origin (Alecian et al. 2013), have a direct impact on the stellar wind by channeling the wind outflow along magneticfield lines, often producing disclike structures around the magnetic equator. Depending on the competition between wind energy and magneticfield energy (characterized by the so called Alfvénradius, R_{A}), and the rotational velocities (characterized by the Keplerradius, R_{K}), a dynamical (R_{A} < R_{K}) or centrifugal (R_{A} > R_{K}) magnetosphere may form (see also Petit et al. 2013). In the former case, confined material in the equatorial disc falls back onto the stellar surface, whereas in the latter case, material is supported by centrifugal forces, forming a quite stable and strongly confined disc structure. Townsend & Owocki (2005) were able to explain the observed Balmerline variability in σ Ori E, a magnetic Bp star, by applying the obliquerotator model. Recently, Owocki et al. (2016) developed a simplified model, the “analytic dynamical magnetosphere” (ADM), in order to provide a framework for the analysis of magnetic winds, and were able to reproduce the observed H_{α}line variations of HD 191612. However, a test of the ADM model for UV resonance lines is still missing, and will be addressed in this paper.
(ii) Accretion and decretion discs: Both premainsequence and a significant fraction of massive mainsequence stars are observed to have circumstellar discs, from which material is either accreted onto the star, or decreted by different mechanisms, respectively. For classical Be stars, Kee et al. (2016) proposed that the destruction of the circumstellar disc by windablation, that is by radiative forces along the disc surface, is the major mechanism, which possibly plays also a significant role for stars that already have ignited hydrogen burning, however are still accreting material through a disc. These objects will be named “accreting highmass stars” within this paper. Although Kee (2015) applied reasonable approximations in order to derive the incident intensity at the disc’s surface layers, the local ionization stages are not treated explicitly. Since the lineforce depends on the linestrength distribution function (e.g. Castor et al. 1975; Puls et al. 2000), which crucially depends on the local ionization stages of the ions and therefore also on the mean intensity, a consistent treatment of the (pseudo) continuum radiative transfer (RT) is needed.
(iii) Fast rotation: Stellar rotation is a natural consequence of the classical starformation scenario, that means of the collapse of an initially rotating molecular cloud. Within the VLTFLAMES Tarantula survey (VFTS), RamírezAgudelo et al. (2013) derived the distribution of projected rotational velocities of (LMC) Otype stars, and found that it peaks at relatively low values (v sin i ≈ 50–100 km s^{−1}), and has an extended tail towards higher ones, thought to be due to binary interaction. For example, the very rapidly rotating Ostar, VFTS102, rotates at a nearly critical rotation rate (Dufton et al. 2011). Nearly critically rotating stars are highly distorted, with a ratio of equatorial to polar radius approaching a value of three half. Additionally, the emergent flux at the stellar surface becomes latitudedependent due to gravity darkening (see von Zeipel 1924 for uniform rotation, and Maeder 1999; Maeder & Meynet 2000 for shellular rotation).
(iv) Binary interaction: By a detailed investigation of binarity within the VFTS, Sana et al. (2013) found that about 50% of the observed Ostar population in the Tarantula nebula have binary companions, that will interact with the primary star during its lifetime via Rochelobeoverflow, or even merging. Roughly 30%–40% of these stars have periods less than six days (see also Sana & Evans 2011), which influences already the mainsequence evolution. Due to the mass transfer, the donor star will be stripped and spins down, whereas the companion will spin up due to the transport of angular momentum (e.g. de Mink et al. 2013). In case of two massive stars being part of the binary system, their winds will collide and lead to strong shocks and Xray emission (e.g. Prilutskii & Usov 1976; Cherepashchuk 1976; Stevens et al. 1992; Pittard 2009), in addition to phaselocked variations of recombination lines (see Sana et al. 2001).
Since all these effects break the symmetry of many observed objects, a 3D treatment of the RT has to be applied to derive the (possibly dynamical) spectral energy distributions. This, however, is computationally expensive due to the dependence of the radiation field on frequency and direction, in addition to the three spatial dimensions. Nevertheless, solution schemes have been developed already from the early 70s on, starting with the longcharacteristics (LC) method (see Jones & Skumanich 1973 and Jones 1973), over the shortcharacteristics (SC) method (see Kunasz & Auer 1988), to finitevolume methods (FVM, see Adam 1990). Within the LC method, the equation of radiative transfer (EQRT) is solved for a given direction along a ray intersecting a considered grid point. Thus, for a 3D grid with N^{3} grid points, and assuming on average N/2 grid points along the ray, N^{4}/2 operations have to be performed to obtain the intensities at all grid points for a single direction. Each operation consists of interpolating the physical properties from the 3D grid onto the ray, and solving the formal integral along all points on the ray. The SC method solves the EQRT only across each grid cell, with incident intensities on the cell surfaces obtained by interpolation from the already known grid points. The number of operations is thus reduced by a factor of N/2. For both methods, the accuracy of the applied interpolation scheme is crucial for the final performance of the obtained solution.
Instead of considering discrete rays, as in the above methods, the FVM solves the EQRT for a single direction, by considering the complete set of rays with the given direction entering the volume of a grid cell, and assuming that the physical properties at the corresponding grid point are representative averages for the complete cell. Though this still gives N^{3} operations for the solution for a single direction, all interpolations are avoided. Thus, when compared to the LC or SC method, the amount of work needed for each operation is reduced considerably. Interestingly, the development of methods underwent an evolution from computationally most expensive methods (LC) to cheaper methods (SC and FVM).
The above mentioned tools are called formal solvers, since they calculate the radiation field for given source and sinkterms. However, the most important problem arises in scattering dominated atmospheres, where the sources (and sinks) depend on the radiation field, and vice versa. Such atmospheres need to be calculated by applying, for instance, the socalled Λiteration and associated acceleration techniques (e.g. Cannon 1973).
Meanwhile, various 2D and 3D codes exist, with widely differing assumptions on the underlying geometry and designated applications. In this respect, one important point is the assumption of local thermodynamic equilibrium (LTE), which cannot be applied in expanding atmospheres, when the radiation field dominates the level populations. For the latter problem, the coupling between the nonLTE (NLTE) rate equations and radiation field is mostly performed using the aforementioned acceleration techniques (e.g. the accelerated Λiteration, ALI). We briefly mention specific codes designated to the multiD RT in stellar atmospheres (for other applications, there are many more such codes; e.g. for multiD codes related to the ionization balance in the interstellar medium, see Weber et al. 2013 and references therein):
WIND3D (Lobel & Blomme 2008) is a 3D FVM code, which has been developed to calculate the line transport in scattering dominated environments using Cartesian coordinates, and with level populations approximated by a twolevelatom (TLA). It adopts the classical Λiteration scheme, which has poor convergence properties for optically thick, scatteringdominated lines. WIND3D is thus restricted to the treatment of weak lines. For those, however, Lobel & Blomme (2008) were able to model the time variations of discrete absorption components (DACs), as observed for the Si IV 1400 Å doublet of the B0.5Ib supergiant HD 64760. Together with the hydrodynamic code ZEUS3D, they were able to reproduce these DACs, supporting the suggestion of Mullan (1984), that they arise from corotating interaction regions (see also Cranmer & Owocki 1996).
ASTAROTH (Georgiev et al. 2006; Zsargó et al. 2006) is a 2D SC code, which is capable of solving the RT in parallel with the NLTE rate equations for axisymmetric problems with nonmonotonic velocity fields. It uses spherical coordinates and includes a (local) ALI scheme. Many tests have been performed by a comparison to 1D spherically symmetric models from CMFGEN, giving errors of a few percent only. Zsargó et al. (2008) applied the code to investigate the H and He ionization stages in the envelopes of B[e] supergiants, and modelled, among others, the H_{α} line.
IRIS (Ibgui et al. 2013a) is a 3D SC code, which, to our knowledge, has been only applied for studying laboratorygenerated radiative shocks (Ibgui et al. 2013b) thus far. Ibgui et al. (2013a) performed test cases (searchlight beam test and 1D planeparallel models), which show an astonishing accuracy for the solution of intensity, mean intensity, radiative flux, and radiation pressure tensor, on nonuniform Cartesian coordinates and including velocity field gradients. By now, the code assumes LTE. However, the inclusion of scattering terms is planned for the future.
MULTI3D (Leenaarts & Carlsson 2009) is a 3D SC code developed to accurately solve the NLTE rate equations in cool FGKtype stars. Amarsi et al. (2016) used this code to predict, among other lines, the O I 777 nm line for a grid of 3D hydrodynamical models, by setting up a 23level model atom. This code, however, has been developed (and optimized) for the application to cool stars (with subsonic velocity fields, at most), and cannot be used in our context.
PHOENIX/3D (Hauschildt & Baron 2006 and other papers in this series) is a 3D LC solver, which is capable of solving the RT together with the multilevel NLTE rate equations. They use spherical, cylindrical or Cartesian coordinates, and implemented a nonlocal ALI scheme. With the extension from Seelmann et al. (2010), arbitrary velocity fields can be included as well. To our knowledge, PHOENIX/3D has not been applied, thus far, to real, multidimensional, problems.
In this paper, we develop a 3D solution scheme to calculate continuum and line source functions in expanding atmospheres. As a first step, we implement a FVM to determine the formal solution, since this method has already been applied in the context of 3D structures (e.g. Adam 1990; Stenholm et al. 1991; Lobel & Blomme 2008), and is supposed to have the lowest computational cost, when compared to the SC and LC methods. To calculate optically thick continua and strong lines, an ALI technique is required, and developed for the first time within the 3D FVM framework by means of a nonlocal approximate Λoperator (ALO, see Sect. 3). In Sect. 4, we present, also for the first time, an error analysis of the 3D FVM for the case of 1D spherical atmospheres, both for optically thin and optically thick environments. This error analysis allows us to understand certain shortcomings of the methods, and to estimate limiting cases that can actually be calculated. As first applications for our 3D solver, we study specific aspects related to the circumstellar discs of accreting highmass stars (Sect. 5), and to the UV line formation in magnetic winds (Sect. 6).
2. Basic assumptions
The problems considered within this paper assume an (almost) stationary atmospheric structure, meaning that the densities, velocities and boundary conditions are assumed to be constant in time, or are at least changing much slower than the radiation field. Thus, we use the timeindependent EQRT (e.g. Mihalas 1978):(1)
In the following, we skip the explicit notation for the spatial (r) and directional (n) dependencies, and write the frequency dependence as subscript, or, when appropriate, skip also this notation. The source function, S ν, and the opacity, χ_{ν}, consist of the sum of continuum and line processes. By splitting the continuum emissivity and opacity into true^{2} and scattering processes, the corresponding source function can be parameterized by a single parameter, ϵ_{c}, the ratio of true absorption to total opacity, which is called the thermalization parameter:(2)
where J_{ν} is the mean intensity, and B_{ν} the Planck function. As a first approach, we treat the line transfer similarly by considering a TLA only. This is well suited to describe resonance lines, but for future applications the complete rate equations need to be taken into account, of course. The line source function in the TLA approach reads(3) (4)
with C_{ul} and A_{ul} being the collisional rate (from the upper to the lower level) and the Einsteincoeffcient for spontaneous emission, respectively. J̄ = 1/4π ∫ I_{ν}Φ_{x} dx dΩ is the profileweighted and frequencyintegrated, angleaveraged intensity, the socalled scattering integral. The profile function, Φ_{x}, has been approximated by a Doppler profile. Further, we did not calculate the profile in frequency space, but rather in the variable x_{CMF}, describing the comoving frame (CMF) frequency shift from line centre, in units of a fiducial Doppler width, :(5)
where ν_{0} and are the linecentre transition frequency, and the fiducial thermal velocity, respectively. The fiducial width is required to enable a depthindependent frequency grid. x_{CMF} is related to the corresponding observer’s frame quantity via x_{CMF} = x_{OBS} − n · V, with the local velocity vector in units of . In most cases, we do not label x explicitly to distinguish between comoving and observer’s frame (nor from the spatial xcoordinate), since the meaning should be clear wherever it occurs.
The profile function (in frequencyspace) is (6)
is the ratio between the local thermal velocity (including microturbulence, see Sect. 3.1) and the fiducial velocity, T is the local temperature, and m_{A} is the mass of the considered ion. With the profile function, normalized in xspace,(8)
the line opacity, which needs to be described in frequencyspace (because the RT and Planck function are formulated w.r.t. frequency), is given by (9) (10)
with (g f ) the g f value of the considered transition, and n_{l}, n_{u}, g_{l}, g_{u} the occupation numbers and statistical weights of the lower and upper level, respectively. Since n_{l} ≫ n_{u} for resonance lines, we have neglected stimulated emission on the righthand side. χ̄_{0} is the frequency integrated opacity^{3}, and is related to χ̄_{L} by(11)
We finally parameterize the continuum and line opacities in terms of the Thomsonopacity, χ_{Th} = n_{e}σ_{e},(12) (13)
where we use k_{L} as a depthindependent parameter, since the ratio n_{l}/n_{e} remains almost constant in the atmosphere for resonance lines (almost frozenin ionization).
The radiation field depends on the source functions via the EQRT, and the source functions depend on the radiation field via Eqs. (2) and (3). Although this coupled problem could be solved directly, at least in principle, we note already here that this would require the inversion of a very large matrix, which is computationally prohibitive in 3D calculations (see Sect. 3.2). Thus, we apply a Λiteration, for which a formal solution (FS) of the EQRT is obtained, given the continuum and line source functions. These are subsequently updated according to Eqs. (2) and (3), given the formal solution of the previous iterate. For large optical depths and low ϵ_{c} (ϵ_{L}), however, the strong nonlocal coupling between the source functions and the radiation field results in the well known convergenceproblem of the classical Λiteration. This problem is based on the fact that, due to scattering processes, photons can travel over many meanfreepaths before being destroyed or escaping from the atmosphere. On the other hand, information is propagated in each iteration step only over roughly one meanfreepath. Therefore, a large number of iteration steps would be required, until a consistent solution was obtained (if at all). Thus, acceleration techniques are urgently needed, and developed in this work in terms of an ALI scheme using operatorsplitting techniques (Cannon 1973).
3. Numerical methods
The ALI scheme can be split into two parts, namely the calculation of the FS, and the construction of a new iterate, by means of an approximate Λoperator. There are various solution schemes to obtain the FS, differing in required computing power, and using different assumptions on the geometry and physical properties of the considered problems, as outlined in the introduction. Among those, the FVM (Adam 1990) is expected to be the fastest and most simple one, and is used here as a first step to tackle 3D radiative transfer problems.
Since we aim at modelling nonmonotonic velocity fields, for which a comovingframe formulation, if at all, is very complicated to implement, we solve the EQRT in the observer’s frame. We also use a Cartesian coordinate system, which has the advantage of constant angular directions w.r.t. the spatial grid. Anyhow, a description of the atmospheric structure in spherical coordinates loses its advantages for the problems considered in this (and future) work.
In the first part of this section, we describe the FVM, and the angular and frequency integration methods. Main focus is put on the development of the ALI scheme tailored to the FVM (Sect. 3.2).
3.1. Finitevolume method
The main ideas of the FVM originate from heat transfer (Patankar 1980), and have been already applied to radiative transfer problems in accretion discs by Adam (1990), as well as to the formation of discrete absorption components in hot star winds (Lobel & Blomme 2008). Although the derivation of the discretized EQRT can be found, for instance, in Adam (1990), we outline the basic ideas in more detail in Appendix A. The final solution scheme reads:(14)
with α, β, γ set to ±1 for directionvector components n_{x}, n_{y}, n_{z} ≷ 0. All quantities except χ_{c}, χ̄_{L}, and the source terms, depend on the considered direction n as well as on frequency. Equation (14) represents a pure upwind scheme, with projected ∆τsteps^{4} calculated from a centraldifferencing approach (following Patankar 1980, see also Appendix A). Due to the upwind scheme, and because the coeffcients a_{i jk} … e_{i jk} ∈ [0, 1], the solution method is unconditionally stable (Adam 1990).
Contrasted to our centraldifferencing approach, Adam (1990) and Lobel & Blomme (2008) used backward differences for the calculation of the ∆τsteps. We have tested both methods, and found superior results when using central differences.
Instead of using an equidistant grid, a gridconstruction procedure is applied in order to resolve the important regions of the atmosphere. In the case of nearly spherical winds, this basically consists of setting up a 3D spherical grid in the first octant, with uniformly distributed latitude and azimuthangles, and with radial coordinates such that ∆τ_{c} and ∆v_{r} are (nearly) equidistant. ∆τ_{c} and ∆v_{r} are the stepsizes in continuum opticaldepth and radial velocity, respectively. The final (Cartesian) x (y, z)coordinates are chosen such that they correspond to the distribution of the x_{s} (y_{s}, z_{s})coordinates from the spherical grid. If not indicated explicitly, we use N_{x} = N_{y} = N_{z} = 133 grid points for the continuum, and N_{x} = N_{y} = N_{z} = 93 grid points for the line formation, distributed over the complete range, [−R_{max}, R_{max}]. We allow for a higher grid resolution within the continuum calculations, since only one frequency point needs to be considered, and we thus have lower computation times anyhow. As a comparison, Lobel & Blomme (2008) used only N_{x} = N_{y} = N_{z} = 71 grid points for their (optically thin) models. A typical grid in the x–zplane is shown in the top panel of Fig. 3. The applied grid resolution is required to properly resolve the resonance zones (where the CMF line opacity and corresponding profile function is nonnegligible), and to obtain reasonable opticaldepth steps within our (firstorder^{5}) method. We emphasize that a large number of grid points would be required to properly describe the radiation field in the optically thick limit, since a firstorder scheme is generally unable to reproduce the (secondorder) diffusion equation. This problem, however, arises only at large optical depths (e.g. when considering the winds of Wolf–Rayet stars), and is only of minor importance for the winds of OB stars, where the total continuum optical depth is of the order unity. To consider only the regions where we actually know the structure of the atmosphere, we calculate the RT within a predefined “calculation volume”, defined by the constraint r(x_{i}, yj, z_{k}) ∈ [R_{∗}, 1.1R_{max}].
Angular integration. The mean intensity at each grid point is obtained from the solutions of the EQRT for many directions, where the distribution of these directions over the unit sphere depends on the used quadrature formula. The resulting intensities are numerically integrated via(15)
with θ_{m}, φ_{l} being the colatitude (measured from the zaxis) and the azimuth (measured from the xaxis), respectively, and w_{m}, w_{l} the corresponding integration weights. The projection factor sin(θ_{m}), and the normalization have been included into w_{m}, w_{l}. Several different techniques and directional distributions have been tested, including the standard trapezoidal rule with equidistant θ, φgrids, as well as using the approach of Lobel & Blomme (2008), who basically use an ansatz dΩ = const, and obtain different φgrids on each θlevel. To obtain a fair angular resolution of the unit sphere, that means without preferring certain directions, we split the integration into a sum over all octants, with the same nodes and weights within each octant. The mean relative errors^{6} of the mean intensity for a zeroopacity model, using N_{x} = N_{y} = N_{z} = 133 grid points and different angular grids, are summarized in Table 1. The corresponding theoretical solution has been calculated from J^{(theo)} = WI_{c}, with the dilution factor, and I_{c} the intensity emitted from the stellar core.
Mean relative errors of the mean intensity for a zeroopacity model, using either Gauss–Legendre quadrature, or the trapezoidal rule with integration nodes from Lobel & Blomme (2008).
The best results have been obtained using a Gauss–Legendre quadrature, for which the integration nodes are fixed, and cannot be chosen arbitrarily. Using this quadrature, already N_{θ}·N_{ϕ} = 512 grid points would be suffcient for an accurate integration, at least in principle. We note, however, that the solution in the outer part of the atmosphere oscillates around the theoretical value, because the star is not resolved by the angular grid in these regions. This effect becomes even worse for the trapezoidal rule. The corresponding errors, however, contribute only weakly to the mean error, because most grid points are located near to the star. Additionally, the dependence of the profile function on the local projected velocity (see Eq. (8)) requires a much finer angular resolution in the line case. Thus, to safely avoid artefacts from the angular integration scheme, we use N_{θ} · N_{ϕ} = 2048 throughout this paper. This is a factor of three lower than the number of angular integration points used by Lobel & Blomme (2008), (N_{θ}·N_{ϕ})^{Wind3D} = 6400. We finally note that the calculation of integration nodes for the Gauss–Legendre quadrature requires the inversion of a large matrix, which directly shows the advantage of a Cartesian grid, for which the angular nodes and weights are the same at all grid points, and can be calculated prior to the complete solution scheme. Using a spherical grid, Gaussian quadrature would be computationally very timeconsuming, because the involved angles w.r.t. the coordinate system depend on position in the atmosphere.
Frequency integration. The integration over frequency is only performed when solving for J̄ in the line case. Due to its simplicity, we apply the trapezoidal rule with equidistant steps,(16)
The integration is performed over the complete frequency range, , with maximum velocity, υ_{∞}, corresponding thermal velocity, υ_{th}(R_{max}), and a rapidly vanishing Doppler profile for x_{CMF}/δ > 3^{7}. To resolve the profile function everywhere in the atmosphere, we require ∆x_{CMF}/δ ≲ 1/3 at each position, which is achieved by defining , and choosing the number of frequency points, N_{F}, such that is the minimum thermal velocity (including microturbulence) of the atmosphere. However, since the comoving frame frequency depends on the projected velocity, the corresponding integration nodes may not be centred around the profile maximum anymore. This issue is only of minor importance and has been checked by a comparison to model calculations that use ∆x_{OBS} ≲ 1/6. In atmospheres with very low microturbulent velocities, a number of N_{F} ≈ 1200 frequency points would be required, for typical values of (if no microturbulent velocities were present), and typical wind terminal speeds, υ_{∞} = 2000 km s^{−1}. When using a rather high microturbulent velocity, υ_{turb} = 100 km s^{−1}, N_{F} is reduced considerably, by a factor of ten. Fortunately, such high values are not untypical in the winds of hot stars (see next paragraph). Again, when compared to Lobel & Blomme (2008), who used N_{F} = 100 frequency points for a thermal velocity of υ_{th} = 30 km s^{−1}, we have a much finer resolved frequency grid, with ≈15 frequency points distributed over the complete line profile, whereas they resolved the profile function with only three frequency points.
Microturbulence. To correctly treat the radiative transport in the line, we need to resolve the resonance zones, and demand that ∆(n · v/υ_{th}) ≈ 1/3. Again, for low microturbulent velocities, we would require a resolution of at least 1200 grid points per spatial dimension (for a spherical wind accounting for both hemispheres). On the other hand, including a large microturbulent velocity (υ_{turb} = 100 km s^{−1}), as applied here, results in a much lower required resolution, both in space and frequency. Due to the linear (cubic) scaling of computation time with the number of frequency (spatial) grid points, this results in a reduction of computation time by a factor of 10^{4}, when compared to models without large microturbulent velocities. Putting it the other way round, and already mentioning here that typical modelcalculations take about 50 min of wallclock time per iteration (and using 16 processors, see next paragraph), a large microturbulent velocity is needed thus far to keep the computation time on a reasonable scale.
Hamann (1981) showed that such microturbulent velocities can indeed be used to correctly model the black absorption troughs observed in the P Cygni profiles of hot star winds. From a theoretical point of view, a large velocity dispersion mimicks the effects of multiply nonmonotonic velocity fields, as originating from the linedriven instability (see, e.g. Lucy 1983; Puls et al. 1993).
Parallelization and timing. In the line case, a number of N_{θ} × N_{ϕ} × N_{F} formal solutions of the EQRT are required to calculate the scattering integrals. In order to obtain accurate angular and frequency integrals, the minimum number of integration nodes is basically fixed, and the computation time can only be reduced further by parallelization. In our case, the most simple and straightforward procedure is a parallelization in frequency, which was done here using OPENMP. The computation time of a typical model with N_{x} = N_{y} = N_{z} = 93, N_{θ} · N_{ϕ} = 2048, N_{F} = 139 spatial, angular and frequency grid points, respectively, is about 50 min of wallclock time per iteration on a 16 CPU INTEL XEON X5650 (2.67 GHZ) machine. As a reference, WIND3D (Lobel & Blomme 2008) requires about 30 min per iteration for their grid parameters (). To obtain a more meaningful comparison, we scale the computation time by the number of applied angular and frequency grid points, and perform a test calculation, which uses the same spatial grid as Lobel & Blomme (2008). Although our CPUs are more modern and faster, our code performs only slightly better than WIND3D, with computation times , per iteration, per angular and frequency point, and per CPU. We note, however, that our algorithm requires at least a factor of two more operations per formal solution, since we additionally compute the (nonlocal) approximate Λoperator in parallel (see next section).
3.2. Λiteration
By now, we are able to construct a formal solution for a given source function. In this section, we discuss the iteration procedure. The following discussion considers the line case alone, with a frequencyindependent background continuum (i.e. constant continuum opacity and source function), assumed to be known (either in form of an optically thin continuum, or from previous calculations in the absence of the line). For convenience, we summarize the basic ideas and corresponding acceleration techniques via the ALI from first principles.
Matrix equation. To show that the Λ formalism can also be applied to our 3D formulation, we derive a matrix equation for the scattering integral, and show that this equation is consistent with an affne representation of the Λ operator. Since, however, the final ALO will be calculated differently, a detailed description of the involved matrices is skipped in the following.
All 3D quantities are expressed as vectors of length N_{x} × N_{y} × N_{z}, by introducing a unique ordering of the (i, j, k)triple:(17)
where i, j, k ∈ [1, N_{x}_{, y, z}], and m ∈ [1, N_{x}N_{y}N_{z}]. Replacing the (i, j, k)indices in Eq. (14) with this unique index, m, and after collecting all intensity terms on the lefthand side, we obtain:(18)
This equation is written in matrix form,(19)
where I_{inc} refers to the boundary condition. Since we use a TLA with given opacity (i.e. the opacity does not change during the iteration), and consider timeindependent boundary conditions, we combine the constant vectors from Eq. (19), Ĩ_{inc} := Q · S^{(c)} + I_{inc}. Inverting the matrix T and integrating over all angles and frequencies gives:(20)
The diagonal matrix Φ_{x} and the vector Φ_{B} describe the local profile function, and the contribution of the boundary conditions and the background continuum to the scattering integral, respectively. As shown below, an explicit calculation of these quantities is not required to obtain the finally used ALO. We note that also for the continuum case, which is calculated close (w.r.t. frequency) to the line, Eq. (20) is applicable, with a different Λmatrix and boundary contribution though. A comparison of Eq. (20) with the common Λoperator formalism (i.e. formally writing J̄= Λ[S L]) directly shows that the Λoperator is an affne operator, that means a linear operator given by the Λmatrix plus a constant displacement vector Φ_{B} (see also Puls 1991), also for our 3D method.
From Eqs. (3) and (20), we could formulate an explicit solution of the radiation field already now. This, however, would require the calculation, storage and inversion of the complete Λmatrix, which is computationally prohibitive:
Firstly, the Λmatrix is a full matrix with N_{x}N_{y}N_{z} × N_{x}N_{y}N_{z} elements, which would require, for typical grid sizes of N := N_{x} = N_{y} = N_{z} = 93, N^{6} ≈ 6.5 × 10^{11} numbers, equivalent to 5.2 TB data to be stored in memory, when doubleprecision numbers are used. Secondly, the Λmatrix elements can be obtained, at least in principle, by inversion (see also Puls 1991),(21)
with e_{n} being the nth unit vector. Thus, N_{x} × N_{y} × N_{z} formal solutions would be needed to calculate the complete Λmatrix, which again is computationally prohibitive on reasonably wellresolved grids.
An iterative solution is therefore the only possibility to solve problems of this kind. Due to the well known convergence problems of the classical Λiteration, we directly focus on the ALI. For completeness, let us mention here a similar approach, the “nonlinear multigrid method” (see Fabiani Bendicho et al. 1997 and references therein), which has even better convergence properties, and, in contrast to the ALI method, does not depend on the spatial resolution of the grid. For simplicity, however, we only implement an ALI scheme. The ALI is an operatorsplitting technique, which splits the original Λoperator into the combination
with an appropriately chosen ALO, Λ^{∗} (see Cannon 1973). Appropriately chosen means that Λ^{∗} should be easily calculated (preferentially in parallel with the formal solution), and easily inverted. Moreover, the ALO should reflect the basic physical properties of the original Λmatrix, in order to significantly boost the convergence.
ALI. Using Eq. (22), where now the first term acts on the current iterate of the source function, and the second one on the previous iterate, we obtain in combination with Eq. (3) (23) where the approximate relation becomes an exact one for the converged solution. Here we have used the notation and definitions of the diagonal matrix, ζ := 1 − ϵ _{L}, and the thermal contribution vector, Ψ := ϵ _{L} · B_{ν}(T), from Puls & Herrero (1988).
Again, all quantities are ordered according to Eq. (17). After some algebra, we obtain(24)
Since also the ALO is an affne operator (analogous to the original Λoperator), the displacement vector Φ_{B} only cancels if it remains constant over subsequent iteration steps, that is, when the background and the boundary conditions remain constant over the iteration. (In realistic, multilevel NLTE calculations, this means in practice that the ALI cycle shows a fast convergence only when the continuum is close to convergence).
Given a previous iterate of the source function, S^{(n−1)}, and
the corresponding formal solution, Eq. (24) is used to calculate the next iterate, S^{(n)}. The detailed choice of the ALO is the crucial point, and finally determines the convergence behaviour. We could, for instance, choose Λ^{∗} = Λ, which would result in the direct solution via inversion, and is computationally not feasible, as discussed above. On the other hand, choosing Λ^{∗} = 0 would result in the classical Λiteration, with the known convergence problems. Olson et al. (1986) showed that an ALO containing only the diagonal of the exact Λmatrix is very effcient, because the matrix (1 − ζ · Λ^{∗}) becomes diagonal, and Eq. (24) could be solved by a simple scalar division. Furthermore, the diagonal, that means the local part, already contributes most to the radiative transfer (at least in the critical optically thick case), and thus, is quite a good approximation for the original Λoperator. Such an ALO corresponds to the well known Jacobiiteration (see also Trujillo Bueno & Fabiani Bendicho 1995 for a thoughtful discussion, also about a Gauss–Seidel method with successive overrelaxation in the context of the ALI).
In 3D calculations, however, a diagonal ALO will not converge fast enough (see Sect. 4.1). To achieve faster convergence rates, a multiband ALO is favourable, as already shown by Olson & Kunasz (1987) for 1D cases, and extended to a 3D, longcharacteristics solver by Hauschildt & Baron (2006). For such ALOs, the matrix (1 − ζ · Λ^{∗}) is sparse, whereas its inverse is a full matrix, and cannot be stored due to the N^{6} scaling of required memory. Therefore, we have already formulated Eq. (24) as a fixpoint iteration, A · S^{(n)} = b, which can be solved for the new iterate by applying Jacobi or Gauss–Seidel methods. We found that a Jacobiiteration, coupled with the storage of the iterationmatrix in coordinateformat (COO)^{8}, is particularly fast and easy, because its computationally most expensive term is a matrixvector multiplication, which reduces to N_{NZ} operations only, where N_{NZ} is the number of nonzero elements (e.g. Tessem 2013).
Constructing the ALO. To construct a multiband ALO as aimed at above, we need to to calculate the corresponding elements of the exact Λmatrix. This could be done, in principle, by using Eq. (20), which would require the inversion of T. Due to the upwind scheme, however, we can simply use Eq. (21), in combination with Eqs. (14) and (16). Since the (m, n)th element of the Λmatrix describes the impact of a nonvanishing source term at point n ↔ (i′ j′ k′) onto a point m ↔ (i jk), the local contribution is given by n = m, whereas the coupling with directly neighbouring points is found from:
n(i − 1, j, k) = m − 1 for n_{x} > 0
n(i, j − 1, k) = m − N_{x} for n_{y} > 0
n(i, j, k − 1) = m − N_{x}N_{y} for n_{z} > 0
n(i + 1, j, k) = m + 1 for n_{x} < 0
n(i, j + 1, k) = m + N_{x} for n_{y} < 0
n(i, j, k + 1) = m + N_{x}N_{y} for n_{z} < 0.
One big advantage of our method is that the exact elements of local and neighbouring terms can be easily calculated from Eq. (14),
We call this ALO “direct neighbour” (DN)ALO, to discriminate from the “nearest neighbour” ALO from Hauschildt & Baron (2006), who use all 26 surrounding grid points and the local term, whereas we are using the local term and the contribution from the six direct neighbours only^{9}. Although it would be possible to include also the other neighbouring terms in our calculations, we note that the calculation of the ALO elements in parallel to the formal solution requires already 50% of the calculationtime in our case, which would increase rapidly when including even more terms for the ALO. On the other hand, the inversion of the ALO, that means the calculation of the new iterate via Jacobiiterations, requires only about 0.5% of the calculation time needed for the complete FS. We emphasize that the ALO actually needs to be calculated only once, because the opacity of the (simplified) TLA remains constant over subsequent iteration cycles. When considering multilevel atoms (as planned in the future), the situation changes, and the opacity depends on the occupation numbers, and thus, also on the radiation field. We therefore implemented the calculation of the ALO in parallel to the formal solution at each iteration step already at the current stage of our code.
To accelerate the iteration scheme further, we implemented the extrapolation technique from Ng (1974, see also Olson et al. 1986). In order to use independent extrapolations, the Ngacceleration is applied in every fifth iteration step. We finally note that the convergence behaviour depends on the grid resolution, which determines the opticaldepth steps, ∆τ_{x, y, z}, and affects the coeffcients c_{i jk}, d_{i jk}, e_{i jk}. The finer the grid, the poorer the convergence behaviour (Kunasz & Olson 1988).
Input parameters for the spherically symmetric models used for our test calculations.
To summarize, we have implemented a 3D FVM in Cartesian coordinates, to determine the FS for a given source function, and to calculate, in parallel, a DNALO applied within an ALI scheme.
4. Spherically symmetric models
Before applying our method to first nonspherical test problems, we have checked its reliability by investigating the convergence properties, and by comparing our 3D solution for a spherical wind with an accurate 1D solution. The 1D solution for the line case has been found by a raybyray solution scheme in p–zgeometry, which has been formulated in the comoving frame and accelerated by an appropriate ALO (Puls 1991). The 1D solution for the continuum transport has been found by applying the Rybickialgorithm (e.g. Mihalas 1978).
Three major error sources of our 3D code will be discussed in the following: Firstly, errors occurring due to false convergence (Sect. 4.1), secondly, errors originating from the incorporation of the boundary conditions and from numerical diffusion (Sect. 4.2), and finally, errors occurring for optically thick media, due to the order of accuracy of our method (Sect. 4.3).
The spherically symmetric models to be compared with have been calculated from a prescribed βvelocity law, the equation of continuity, and a temperature structure from the 1D code (which plays almost no role in our test cases).
For opacities and source functions, see Sect. 2, and the required electron density, n_{e}, has been derived from a completely ionized H and He plasma, with N_{He}/N_{H} = 0.1. For all following model calculations, we used a fixed set of prototypical input parameters, summarized in Table 2. These parameters roughly correspond to the wind from ζ Pup, when assuming an unclumped wind. Within the line transfer, we considered a generic UV resonance transition, with different linestrengths following Eq. (12). To calculate the thermal width, we used m_{A} = 12 m_{p}, m_{p} being the proton mass. The total width of the line profile, however, is mainly controlled by the (large) turbulent velocities. Different optical depths,
(for the model considered in Table 2), and scattering properties of the model atmosphere were simulated by varying the scaling factors, k_{c}, k_{L}, and the thermalization parameters, ϵ_{c}, ϵ_{L}, in the continuum and line case, respectively.
Fig. 1. Convergence behaviour of the continuum transfer with ϵ_{c} = 10^{−6} (left), and of the line transfer with ϵ_{L} = 10^{−6} (right), for different acceleration techniques applied to the RT, and for a spherical model atmosphere. The optical depth scale varies according to k_{c} = 10, k_{L} = 10 (top) and k_{c} = 100, k_{L} = 10^{5} (bottom). See text. 

Open with DEXTER 
4.1. Convergence behaviour
To test the convergence behaviour of our ALI implementation with corresponding ALO, we only considered scattering dominated atmospheres by setting ϵ_{c} and ϵ_{L} to 10^{−6}. We discuss the continuum transfer in the absence of a line, and the line transfer assuming an optically thin continuum. In both cases, we applied the classical Λiteration, and compare to ALI schemes using a diagonal or DNALO, with Ngextrapolation switched on or off.
Pure continuum. The left panel of Fig. 1 shows the maximum relative corrections of the mean intensity after each iteration step, for an intermediate grid resolution of N_{x} = N_{y} = N_{z} = 93 spatial points and N_{θ}N_{ϕ} = 968 angular points. Two different continua are shown, referring to an optically thick (k_{c} = 100), and marginally optically thick (k_{c} = 10) case. Obviously, the classical Λiteration converges only very slowly (if at all) in the high opticaldepth regime, requiring iterations until the convergence criterion of maximum relative corrections being less than 10^{−3} is fulfilled. We emphasize that a (steep) gradient of the “convergence curve” is required to achieve actual convergence, rather than this (arbitrarily chosen) number.
For the optically thick problem, even the diagonal and also the DNALO have severe convergence problems, with and , respectively. To ensure fast convergence, the Ngacceleration is urgently needed, since it boosts the relative corrections significantly, reducing the number of required iterations for the optically thick problem to (lower left panel of Fig. 1).
Line case. With the same setup as above, we display the maximum relative corrections of the scattering integral after each iteration step in the right panel of Fig. 1. We applied two different linestrength parameters, k_{L} = 10, 10^{5}, which correspond to a weak and strong line, respectively. Since the line transport is restricted to the finite widths of the resonance zones, and is therefore intrinsically much more local than the continuum, the convergence behaviour is accelerated significantly already by the diagonal and DNALO. For the strong line, the required number of iterations until convergence is reduced from to , and , for the diagonal and DNALO, excluding or including the Ngacceleration, respectively (lower right panel of Fig. 1).
4.2. Boundary conditions and zeroopacity models
Two important points still to be explained refer to the incorporation of boundary conditions and to numerical diffusion. The latter can be best understood by considering continuum models with an opacity set to zero.
Fig. 2. Boundary conditions for two different points, r_{1}, r_{2}, and different directions, n_{1}, n_{2}. Green: a ray originating from the stellar photosphere. To calculate the intensity at r_{1} in direction n_{1}, the intensities at the neighbouring grid points, and , need to be known. A boundary condition is required for grid point , while the intensity at results within the “normal” RTscheme. Red: a ray originating from outside the photosphere. For the grid point r_{2}, a boundary condition has to be specified at the corresponding phantom point, . 

Open with DEXTER 
Boundary conditions. At the outer boundary, the intensities coming from outside are set to zero, whereas those coming from inside are calculated within the RT scheme, and need not to be specified. Close to the star, the situation is more complicated: The inner grid points are generally located close to the boundary, but not directly on it, due to the difference of Cartesian and spherical grids. The intensities at those grid points are calculated by the standard FVMRT (Eq. (14)), but using different grid cells, defined by the intersection(s) of the original grid cell with the (spherical) photosphere (see Fig. 2). The intensity for certain directions needs to be specified at those socalled phantom points. Figure 2 displays phantom points corresponding to two distinct grid points, r_{1}, r_{2}, and different ray directions, n_{1}, n_{2}. Radiation originating from the stellar surface (direction n_{1} in Fig. 2) is set to B_{ν}(T_{eff})^{10} at the corresponding phantom point, . For some grid points and ray directions, however (e.g. direction n_{2} in Fig. 2), even intensities incident onto the photosphere need to be specified at the corresponding phantom point.
Within the corehalo approximation used here, inwards directed intensities should, in principle, be set to zero at the lower boundary, that is I_{phantompoint} = 0. If, on the other hand, the lower boundary is located at significant optical depths (as for the majority of the cases considered here), a specification within a simplified diffusion approximation, I_{phantompoint} ≈ Bν(Teff), is more appropriate. Although the diffusion approximation is no longer justified when concentrating on purely scattering atmospheres, backscattering of photons in such environments mimicks a similar effect, at least if the optical depths are not too low.
Mean and maximum relative errors of our 3D solution scheme, when compared to an “accurate” 1D solution.
We have tested this issue by considering an optically thin model as the most extreme testbed, using both alternative descriptions for those critical directions. The mean relative errors for both alternatives are of the same order, and do not significantly differ from those arising under more physical conditions discussed later (with larger optical depths at the lower boundary, see Sect. 4.3 and Table 3). Thus, we apply B_{ν}(T_{eff}) as the inner boundary condition for those critical rays as well, also because this procedure is less timeconsuming, since it avoids conditional clauses in the innermost loop of the code.
Zeroopacity models. Numerical diffusion is a major error source within the FVM (it also occurs in the SC method, though to a lesser extent, see Kunasz & Auer 1988 and Ibgui et al. 2013a for an analysis of this effect in the context of multiD SC methods). Quite generally, rays propagating parallel to the grid axes are nearly undisturbed by this effect, whereas those propagating at different angles are effectively widened, due to the diffusion of intensity into neighbouring cells. This diffusion results from the finite size of the grid cells, and the competition between the ∆x/∆y, ∆x/∆z and ∆y/∆z terms in the discretized EQRT (Eq. (14)), for any given direction, and can only be minimized by increasing the grid resolution.
To obtain a qualitative measure of numerical diffusion errors, we performed a searchlight beam test with ray direction n= (1, 0, 1), by setting the opacity to zero. For consistency, the critical boundary conditions have been set to zero.
The top panel of Fig. 3 shows the contours of the specific intensity in the x–zplane, and the middle panel displays the intensity through an aperture perpendicular to the ray direction. In our projection, the aperture appears as a straight line in the x–zplane, specified by its distance from the origin and the impact parameter. Finally, the bottom panel of Fig. 3 shows the intensity along the considered direction, at the centre of the beam.
Evidently, numerical diffusion plays a crucial role. We expect corresponding errors to be most severe when photons propagate over large distances, for instance, for optically thin continua, or, in the line case, before they hit the resonance zones. Since the resonance zones are mostly quite narrow, while the pathlength of freely propagating line photons is usually quite large (at least if the continuum is comparatively weak, as in most realistic conditions), numerical diffusion errors are of major significance for the line transfer.
To address the impact of numerical diffusion on the final solution, we consider the mean intensity for zeroopacity models. Figure 4 shows the mean intensity (scaled by its theoretical value obtained from the dilution factor) for such models on spherical surfaces at two distinct radii, r = 1.1R_{∗} and r = 3R_{∗}. We find a clear pattern of the shape of the mean intensities: Close to the star, the mean intensities are reasonably accurate on the axes, in contrast to the regions away from the axes, where they become overestimated. Far from the star, the situation reverses, with reasonable results away from the axes, and an overestimate on the axes. This behaviour is explained in the following.
Fig. 3. Top: specific intensities for direction n= (1, 0, 1) in the x–zplane, for a zeroopacity model. Overplotted is the spatial grid with a typical size of N_{x} = N_{y} = N_{z} = 133 grid points. The thick lines indicate the theoretical boundary of the beam. Middle: projected intensity profile through the area indicated by a straight line in the top panel, covering all rays with direction n and impact parameter p = ±2R_{∗}, and theoretical intensity profile denoted by dotted lines. The central distance from the origin is 3R_{∗}. Bottom: specific intensity along n, compared to the theoretical value indicated by the dotted line. 

Open with DEXTER 
Fig. 4. Contours of the mean intensity for a zeroopacity model, normalized by its theoretical value, on spherical surfaces, r = 1.1R_{∗} (left) and r = 3R_{∗} (right). The line of sight is along the vector n= (1, 1, 1). N_{x} = N_{y} = N_{z} = 133 grid points have been used. 

Open with DEXTER 
For a given grid point on or close to a coordinate axis, and far from the stellar surface, corerays, that means those originating from the stellar surface, remain nearly undisturbed by numerical diffusion. Without diffusion, only such corerays would contribute to the mean intensity. Due to numerical diffusion, however, also noncore rays contribute, that means those which form a certain angle w.r.t. the considered grid axis, since they have been fed with intensity by corresponding corerays propagating in the same direction (widening of the effective aperture, see above). Consequently, the mean intensity becomes overestimated.
For grid points far from the star, and away from the major axes, corerays and noncore rays are both affected by numerical diffusion, resulting in an under and overestimation of the intensity, respectively. Consequently, there is a significant cancellation of both effects, and the mean intensity remains close to its expected value.
At points located on the grid axes close to the star, numerical diffusion plays only a minor role, mostly because the contributing noncore rays are propagating almost perpendicular to the considered axis, which means parallel to one of the other axes, with negligible diffusion errors.
Away from the axes, and close to the star, contributing noncore rays are significantly inclined w.r.t. the coordinate axes, and thus strongly fed by diffusion effects. Thus, the mean intensities become overestimated.
Due to the different effects for different ray directions and for different regions in the atmosphere, any symmetry will be broken. This error cannot be avoided, and is minimized only for higher grid resolutions. As the important part of the radiative transfer is mainly located near to the star (where the densities are largest), and the numerical diffusion errors are not too large in this regime (up to a radius of r ≲ 3–4R_{∗}), the 3D solution scheme should deliver at least qualitatively correct results.
4.3. Variation of optical depth
We finally compare the results from our 3D solutions for a spherically symmetric wind with corresponding ones from 1D solutions using spherical coordinates, as a function of optical depth. Within the continuum transport, the optical depth has been varied by increasing the linear scaling factor, k_{c} = (10^{0}, 10^{1}, 10^{2}), which defines a radial opticaldepth scale, τ_{r} = (0.17, 1.7, 17.0), at the lowermost point. The solutions for the mean intensities, together with the 1D solution, are shown in Fig. 5, top panel. The line transitions have been calculated for the same model, assuming an optically thin continuum in order to extract the error from the line transport alone. The adopted linestrength parameters, k_{L} = (10^{0}, 10^{3}, 10^{5}), describe a weak, intermediate and strong line, respectively. The corresponding solutions are shown in the bottom panel of Fig. 5. To ensure convergence, we have used the DNALO together with the Ngacceleration (see Sect. 4.1) for all test calculations. Thus, the differences between the 1D and 3D solution originate from the formal solution scheme alone, and not from a false convergence. The mean and maximum relative errors are summarized in Table 3.
Fig. 5. Top: mean intensities (scaled by the emitted intensity from the stellar core, I_{c}) for the continuum transport, with ϵ_{c} = 10^{−6}, and k_{c} = (10^{0}, 10^{1}, 10^{2}) corresponding to τ_{r} = (0.17, 1.7, 17.0), from left to right. Bottom: line source functions (scaled by I_{c}) assuming an optically thin continuum, and with ϵ_{L} = 10^{−6}, and k_{L} = (10^{0}, 10^{3}, 10^{5}) from left to right. The red line corresponds to the accurate 1D solution, the blue line is the solution along the xaxis, and the dots (mainly located in between the red and the blue line) correspond to the solution of (arbitrarily) selected grid points (only a subset of all grid points is shown to compress the images). 

Open with DEXTER 
Fig. 6. Emergent line profiles for spherically symmetric models with different linestrengths, k_{L} = (10^{0}, 10^{3}, 10^{5}), from left to right. The red line corresponds to the accurate 1D solution, and the black line is the solution corresponding to our 3D calculations. 

Open with DEXTER 
The mean errors are increasing together with the optical thickness, because the FVM becomes a firstorder scheme for increasing opticaldepth steps^{11}. The errors in the line case do not exceed roughly 25%, and are generally lower than those for the continuum, because the RT is much more local (see above)^{12},
and thus, the error is not being propagated through the complete grid. On the other hand, numerical diffusion plays a larger role (see Sect. 4.2), resulting in minimum errors of ≳10% even for very weak lines.
For optically thick continua (with τ_{r} ≳ 20), the continuum transfer breaks down, giving mean errors larger than 100%. Figure 5 shows that the errors for the continuum and line transport are mostly due to an overestimation of the mean intensities and scattering integrals, respectively. With respect to the local distribution of the errors, we find a similar behaviour as for the zeroopacity models (on vs. away from the axes for different radii). The maximum errors are mostly found at the same locations as in Fig. 4.
Finally, we conclude that the overall error is a combination of errors introduced by numerical diffusion, and the firstorder scheme for optically thick environments. The line transport is generally more reliable than the continuum transport, which should be treated with higher spatial resolution.
4.4. Emergent flux profile
To calculate the line profiles, we implemented a postprocessing LC method. Within this method, the converged source function is used to derive the formal solution along characteristics in a cylindrical coordinate system, with the zaxis aligned with the line of sight. The emergent intensity is finally integrated over the projected stellar disc, which gives the flux (see, e.g. Lamers et al. 1987; Busche & Hillier 2005; Sundqvist et al. 2012, for more details). To obtain the source functions at each position of a given characteristics, a 3D trilinear interpolation is performed in logspace. Though this could be applied to the velocity components and opacities as well, we calculate those properties from analytic expressions whenever possible. This way, we avoid errors introduced by the interpolation scheme, which have been found to be of the order of only few percent for a βvelocity law, but have rather strong influence for the ADM models, where the interpolation of shearvelocities is not a simple task.
The emergent profiles for the same models as above are shown in Fig. 6, with differences between the 1D and 3D models originating from the line source function alone. Generally, when compared to the 1D solution, the line profiles from our 3D code overestimate the emission part due to the larger source functions. The absorption edge is slightly redshifted from its theoretical value at x_{obs} ≈ 1+v_{th}/v_{∞}, because the calculation volume extends only up to 13.2R_{∗} (where v(13.2R_{∗}) = 0.92 v_{∞}). This issue, however, could be improved by enlarging the size of the calculation volume.
Overall, despite the slightly enhanced emission peak, we find that the line profile can be reproduced by our FVM, in combination with an accurate postprocessing routine, for a wide range of linestrength parameters. Thus, our method allows for a qualitative interpretation of line profiles even for the most extreme test cases, that is for strong scattering lines.
5. Windablation
Using 3D radiationhydrodynamic simulations, Kee (2015) and Kee et al. (2016) modelled the ablation of circumstellar discs around massive stars, due to radiative line driving. They showed that a significant lineforce arises due to the coupling of nonradially streaming photons to the nonradial velocity field of circumstellar discs (see also Kee et al. 2016, their Fig. 1). The lineforce has been calculated within a Sobolevapproach, by means of linestrength distribution functions. Contrasted to the original formulation by Castor et al. (1975), they followed the parameterization by Gayley (1995). The full 3D line acceleration can then be written as^{13}
Stellar and wind parameters for the windablation model.
where n and v describe the direction of the considered ray and the velocityvector, respectively. ρ is the density, Q̄, Q_{0}, and α describe the linestrength distribution, and were taken from the calibration of Puls et al. (2000, their Table 2), for the considered T_{eff}.
For optically thin continua (e.g. in classical Be stars), the incident intensity I(n) can be directly replaced by the intensity originating from the stellar core, I_{c}, and Eq. (35) can be solved by quadrature, for a given density and velocity structure. For accreting highmass stars (see Sect. 1), that means for massive objects in their late formation phases, however, the circumstellar discs are optically thick, and at least two major problems arise:
Firstly, due to absorption and scattering processes, the incident intensity at a considered point needs to be calculated by a global solution of the radiation field, which is very time consuming in hydrodynamic simulations. Kee (2015) developed an effcient method to delimit the contribution by either calculating the absorption part alone (giving a lower limit of the incident intensity), or assuming the disc to be optically thin (giving an upper limit). A comparison between the irradiation obtained from their method to the irradiation obtained from our 3D code (including scattering of photons) will be presented in a forthcoming paper, and shall not be discussed here.
In this paper, we only consider the second problem of optically thick environments: Since the disc partly blocks the irradiation from the star, the radiation field might become considerably reduced. However, it is this (ionizing) radiation field, which mainly determines the ionization stages of the considered plasma, and consequently might influence the linestrength distribution function. To address this issue, we proceed as follows: First, the radiation field for a specific hydrodynamic structure (see below) is computed by our 3D code, and the resulting mean intensity is translated to a corresponding radiation temperature (using J_{ν} =: W · B_{ν}(T_{rad}) with dilution factor W; the derived radiation temperature would correspond to T_{eff} if the star had an optically thin, spherically symmetric atmosphere). Corresponding linestrength parameters could again be obtained from Puls et al. (2000). In regions where the local radiation temperature is similar to the effective one, one can safely assume that the parameters of the linestrength distribution remain at their original input values, and indeed can be used to calculate the line force throughout all following hydrotimesteps. If, on the other hand, T_{rad} differed significantly from T_{eff}, this would imply that these parameters would need to be consistently adapted within the hydrodynamical evolution.
For our analysis, we used a wind+Kepleriandisc model similar to the initial conditions for the accreting O7star system as considered by Kee et al. (2016). This model describes such objects as already defined in the introduction as accreting highmass stars (see Hosokawa et al. 2010; Kuiper et al. 2016). The wind and stellar parameters are given in Table 4, following Kee et al. (2016). A radial optical depth of the disc, τ_{disc} = 1400, has been adopted. We approximated the continuum by pure Thomsonscattering, ϵ_{c} = 0, to ensure frequency independence. This is a fair assumption for the 500–2000 Å range, where the majority of linedriving happens (e.g. Puls et al. 2000). Of course, we would expect thermalization in the disc’s deeper layers. Due to the dominating ρ^{−α}dependence of the line force (Eq. (35)) and the large densities inside the disc, however, most of the windablation occurs at the surface layers, and we do not need to care about the details in the inner parts. This fact is even more important, since it allows us to apply our 3D FVM method, although being aware of the large errors of the continuum transfer for optically thick media. To ensure that the transition region from the wind to the disc is not subject to (larger) numerical uncertainties, we have performed a test calculation with doubled grid resolution . Although we found, as expected, differences in the inner part of the disc, our results for the outer part and the wind region are (almost) identical. We can, therefore, safely assume that the obtained solutions, at least in the aforementioned regions, are only mildly affected by numerical artefacts.
Fig. 7. Density structure for the windablation model with τ_{disc} = 1400, in the x–zplane. 

Open with DEXTER 
Fig. 8. As Fig. 7, but for the radiation temperature. Additionally, the density contours corresponding to a decrease in line force by f = 10 (solid) and f = 100 (dashed) are displayed. Both contours indicate the transition region from wind to disc (see text). The grey colour corresponds to radiation temperatures less than the colourcoded minimum value. 

Open with DEXTER 
The density structure and radiation temperature (the latter computed by our code) are shown in Figs. 7 and 8, respectively. The radiation temperature in the wind (here: along the zaxis) exceeds the effective temperature by a factor of roughly 1.25. In order to ensure that this is not a numerical effect, we have checked this issue by calculating the same wind model, however applying an optically thin disc with τ_{disc} = 1.4 × 10^{−3}. For such a model, T_{rad} and T_{eff} turned out to be fairly identical. We thus conclude that the enhancement of radiation temperature in our original model is due to additional irradiation of the wind from the disc, by scattering off photons from the disc. Most likely, this effect will induce latitudinal lineforce components (also to be addressed in a forthcoming paper).
Windablation dominates in the transition region between wind and disc. We define this region by calculating the decrease in lineforce by a certain factor, f, due to density effects alone, that means assuming the same ionization stages and the same velocity structure. Such a reduction of lineforce or lineacceleration is easily cast into an enhancement of density via Eq. (35),
where the radiusdependent quantities from the wind can be measured along the zaxis. In this picture, f should be chosen such that the corresponding decrease in lineforce represents the border from the wind region to the region where the lineforce is negligible (i.e. inside the disc). As a first guess, we adopted f = 10, and display the corresponding density contour in Fig. 8. Since a factor f = 10 seems to be somewhat artificial, we additionally display the density contour corresponding to f = 100.
From our simulations, we then find that both contours are located within a range of T_{rad} between roughly 31 and 33 kK, which is of the same order as the effective temperature, T_{eff} = 36 kK. We thus conclude that the ionization stages at the disc surface are not changing too much, when compared to the ionization stages in the wind, and that the linestrength parameterization of the wind can also be used to calculate the lineforce at the surface of such optically thick circumstellar discs. Due to significant scattering of photons off the disc, a multiD description of the radiative transfer might need to be incorporated into the hydrodynamic simulations, to account for all forcecomponents.
6. Dynamical magnetospheres: HD 191612
As a first application to line transitions, we modelled UV resonance lines in dynamical magnetospheres, that means atmospheres which form in slowly rotating magnetic OBstars (in contrast to the socalled centrifugal magnetospheres, which form in fast rotating magnetic OBstars). As a prototypical case, we considered the Of?p star HD 191612, which has a negligible equatorial rotation speed of v_{rot} ≈ 1.4 km s^{−1} (Howarth et al. 2007; Sundqvist et al. 2012). Marcolino et al. (2013) already calculated corresponding resonance lines for this star, by extending the 3D formal solver developed by Sundqvist et al. (2012) to a “3D Sobolev with exact integration” method (SEI, Lamers et al. 1987), and applying this method to a set of 100 twodimensional MHDsimulation snapshots, equidistantly distributed over the azimuthangle to enable a 3D description of the atmosphere. At least for the H_{α} line (where the source function is taken from prototypical 1D NLTEcalculations), such a patchingtechnique produces quite similar results as full 3D MHD simulations (see udDoula et al. 2013). In Sect. 6.1, we use the same simulations as a benchmark for our 3D code, and compare the obtained line profiles to those from Marcolino et al. (2013). In Sect. 6.2, we calculate analogous line profiles for the ADM model developed by Owocki et al. (2016), to investigate in how far their simplified description of the magnetosphere can be used as a reasonable substitute for elaborate MHD simulations. We already note here, that such a simplified approach would be favourable to MHD simulations, because it provides (within the applied approximations) an average, steady state solution for the magnetospheric structure, and avoids timeconsuming hydrodynamic simulations.
6.1. MHD models
To understand the behaviour of the line profiles presented below, we first explain the basic characteristics of (nonrotating) magnetic winds. For a more detailed discussion, we refer the reader to the seminal work by udDoula & Owocki (2002) and udDoula et al. (2008). These authors introduced a magnetic dipole field as an initial condition, and evolved the (initially spherical) stellar wind according to the MHD equations. Within ideal MHD^{14}, the material follows the (closed) magneticfield lines in regions where the magnetic energy exceeds the kinetic energy of the wind (close to the star), whereas, in the opposite case, the field lines follow the (almost radial) mass flow (far from the star). The border of both regions can be roughly described by the AlfvénRadius, , with windconfinement parameter and B_{p} the polar magneticfield strength evaluated at the stellar radius (see udDoula & Owocki 2002).
Within closedfield regions, material originating from opposite footpoints shocks (and accumulates) in the equatorial plane. Due to the 1/ρ^{α}dependence of the lineforce (see Eq. (35)), the netforce becomes dominated by gravity, and produces an inflow along the magnetic field lines in a “snakelike” pattern.
In the openfield regions, the presence of the magnetic field, together with a frozenin mass flow, results in a density decrease when compared with spherically symmetric models, due to the fasterthanradial expansion of the flowtube area (see Fig. 7 in udDoula & Owocki 2002). Consequently, the lineforce becomes increased, resulting in higher terminal velocities than in 1D nonmagnetic models. A single snapshot and an azimuthal average of the applied MHD simulations are shown in Fig. 9.
Based on such MHD simulations, Sundqvist et al. (2012) calculated corresponding H_{α}line profiles, while Marcolino et al. (2013) investigated the UV resonance line formation. To remain consistent with the calculations by Marcolino et al. (2013), we apply ϵ_{L} = 0, and use their description of the linestrength parameter, κ_{0}, originally introduced by Hamann (1980). κ_{0} is related to the linestrength parameter from Eq. (12) by
with I_{He} = 2 and Y_{He} = N_{He}/N_{H} = 0.1, the number of free electrons per helium atom, and helium abundance by number, respectively.
Although we use a microturbulent velocity of v_{turb} = 100 km s^{−1} for the determination of the source function, we calculate the final line profile (somewhat inconsistently) for v_{turb} = 50km s^{−1}, as done by Marcolino et al. (2013). The line profiles obtained from our 3D code and the SEI line profiles from Marcolino et al. (2013) are displayed in Fig. 10, for two different linestrength parameters, κ_{0} = 0.1 and 1.0, respectively.
Fig. 9. Upper panel: density structure for an example snapshot from the MHD simulations for HD 191612, as performed by Sundqvist et al. (2012). Lower panel: azimuthal average of the MHD simulations. In both figures, the density has been normalized by a typical downflow density, , with Ṁ_{B}= 0 from Table 5, and v_{esc} ≈ 800 km s^{−1} the photospheric escape velocity. The velocity field is displayed by arrows, with the length of the velocity vectors limited to 0.5 v_{esc}. We additionally show the dipole magnetic field of the ADM models used in Sect. 6.2 (solid lines, and thick solid line for R_{A} = 2.7R_{∗}). The corresponding magnetic axis is aligned with the zaxis. The grey colour corresponds to densities outside the range indicated on the right. 

Open with DEXTER 
The agreement between the two methods is excellent. The minor differences in the emission part are related to two effects: Firstly, the methods for determining the source functions (SEI implying very narrow resonance lines vs. FVM accounting for much broader ones, due to v_{turb} = 100 km s^{−1}) are quite different, and a certain deviation must be present. Secondly, the (general) overestimation of the scattering integrals and thus source functions due to the FVM might play a role as well. Also the absorption parts of the line profiles observed equatoron (lower panel of Fig. 10) are not perfectly matched. This (small) effect is most likely simply due to different formulations of the numerical solvers.
A comparison of these line profiles with those from corresponding spherically symmetric models has already been performed by Marcolino et al. (2013), and we summarize only the most important characteristics: (1) The absorption trough for poleon and equatoron observers extends beyond the 1D terminal velocity, as expected from the MHD atmospheric structure. We note that such a large extension has not been observed for HD 191612. (2) The emission for equatoron observers is reduced (compared to the 1D case), due to the lower densities in the emission plane (e.g. the polar plane for linecentre frequencies with x_{OBS} = 0). (3) The particular form of the line profiles is determined by the different mapping of projected velocities for different observer directions.
Fig. 10. UV resonanceline profiles for the MHD models, as obtained from our 3Dcode (black) and from the SEI method by Marcolino et al. (2013; red). Two different linestrength parameters, κ_{0} = 0.1 and 1.0, have been used. For convenience, the line profiles for κ_{0} = 0.1 have been shifted vertically by 1.5. The upper and lower panels show the synthetic line profiles for poleon and equatoron observers, respectively. The abscissa has been scaled to v_{∞} = 2700 km s^{−1}, the “observed” 1D value applied by Marcolino et al. (2013). 

Open with DEXTER 
Given the overall agreement of the two different methods, we conclude that the SEI and our 3D FVM solutions are consistent. Additionally, we are highly confident that the line formation is described correctly (at least qualitatively), since both methods are completely independent. Under this assumption we are able to study the UV line formation within the ADM model.
6.2. ADM models
Owocki et al. (2016) developed an analytic description of dynamical magnetospheres, in order to set a framework similar to the βvelocityfield prescription for spherically symmetric winds. This ADM formalism provides a timeindependent, steadystate solution for dynamical magnetospheres, which is comparable to the average of several MHDsimulation snapshots, and has been corroborated by a comparison of synthetic H_{α} lines with observations. The formation of resonance lines within the ADM framework, however, has not been analysed yet, and is the focus of this section. For that purpose, we aim at modelling the MHD atmospheric structure from above with the ADM method, and compare the resulting line profiles.
Within the ADM method, Owocki et al. (2016) divide the atmosphere into two major zones. The border between both regions is given by the condition r_{m} = R_{A}, where the apexradius, r_{m}, is defined as the distance between the origin and the intersection of magnetic equator and closed dipole magneticfield line attached to a considered point (see left panel of Fig. B.1 for clarification). In the following, we call these two regions the “closedfield region” (r_{m} < R_{A}) and the “outer wind” (r_{m} > R_{A}).
The closedfield region consists of three different components:
windupflow component: The magnetic loops are fed with material ejected from the stellar surface. The matter flow follows the dipole magneticfield lines, with absolute velocities calculated from a βvelocity law, using β = 1.
postshock component: The collision of outflows following the Bfield lines from opposite foot points leads to a shock at the magnetic equator, resulting in a hot and dense postshock region. The extent of this region is controlled by a (dimensionless) cooling parameter, χ_{∞}, where 1/χ_{∞} describes the effciency of radiative cooling by Xray emission (see udDoula et al. 2014 for details). In test calculations, however, this component turned out to have only very small influence on the UV line formation. Thus, and to keep the model as simple as possible, we neglect the postshock component in this work.
cooleddownflow component: As the postshock gas cools, its density increases, and the lineforce decreases. Thus, the cooled and compressed gas is pulled back onto the stellar surface by gravity, resulting in a downflow starting at the magnetic equator. The gas is accelerated from zero velocity along the Bfield lines to the escape speed at the stellar surface.
For their H_{α} analysis, Owocki et al. (2016) only considered the cooled downflow, because of the mostly larger densities of this component. Since the infall occurs in episodic infall events, the closedfield region is actually highly structured, and the authors found rather large clumping factors ⟨ρ^{2}⟩/⟨ρ⟩ (of the order of several tens). Under the assumption of clumps that are optically thin then, this clumping factor can be used to translate the actual (structured) densitydistribution to the mean opacities and emissivities of recombination lines (ρ^{2}processes, see e.g. Puls et al. 2008). For the UV resonance line formation (linear in ρ), microclumping has no direct impact on the mean opacities. Therefore, and because of the different densities and velocities within the upflow and downflow components, an explicit description of the structured medium is required when considering UV resonance lines. As it is not a priori clear how to treat the combination of the above mentioned components, we consider four different approaches, and model the closedfield region by:

(i)
Applying only the cooleddownflow component.

(ii)
Introducing a statistical treatment, where the probabilities of using either the windupflow or the cooled downflowcomponent when calculating the radiative transfer are here defined as
This approach preferentially chooses the component with higher density and lower velocity^{15}, in other words that component with the larger timescale for the matter flow.

(iii)
Introducing fluxtubes that alternating consist of the downflow and upflow component.

(iv)
Applying only the windupflow component.
The models are ordered such that the contribution of the windupflow component is increasing from model (i) to (iv).
As a zerothorder approximation, Owocki et al. (2016) model the outer wind (at r_{m} > R_{A}) by the windupflow component, that means by a flow following closed magneticfield lines even in that region. This is a fair assumption for modelling the polar regions, since it accounts for the faster than radial decline of the density (see also Owocki & udDoula 2004). Moreover, the magnetic field lines are nearly radial in these regions, thus resulting in a nearly radial outflow similar to the MHD simulations. On the other hand, the velocity vectors near the equatorial regions are modelled with a large latitudinal component, whereas they are radially directed within the (more realistic) MHD simulations. Thus, a match of the ADM and MHD magnetospheric structure in the equatorial region cannot be achieved within the standard formulation. With respect to UV line formation, this is the major drawback of the ADM formalism, and will influence the line formation (see below).
To set the base density, Owocki et al. (2016) introduced the massloss rate of the star if it had no magnetic field, Ṁ_{B}= 0, which determines the loopfeeding rate. With the input parameters from Table 5, the dynamical magnetosphere can be modelled according to the recipe from Owocki et al. (2016)^{16}. We used the values of Ṁ_{B}= 0 and in Table 5, right panel, to adapt the ADM model to the MHD simulations. For our model parameters, the AlfvénRadius, R_{A} = 2.7R_{∗}, has been calculated from the massloss rate, terminal velocity and magneticfield strength. We stress that the adopted massloss rate is not necessarily the “true” one, nor the massloss rate the star would have if no magnetic field was present. For the applied ADM models, the resulting density stratification, magneticfield lines and velocity vectors in the xzplane are shown in the left panel of Fig. B.1. Here and in the following, the equatorial plane coincides with the plane of the magnetic equator, since we assume the magnetic axis to be aligned with the zaxis.
Compared to the MHD structure (see Fig. 9), the ADM densities in the closedfield region are best represented by model (ii) and (iii), that is by a combination of downflow and upflow component. In the outer wind near the equator, the densities are underestimated due to the aforementioned different description of the velocity field (see Fig. B.1, left panel).
Fig. 11. UV resonanceline profiles obtained from our 3D code, for the MHD simulation (black, as Fig. 10), and the different ADM models, (i)–(iv) (see text). The linestrength parameter has been set to κ_{0} = 1. The upper and lower panels show the synthetic profiles for poleon and equatoron view, respectively. As in Fig. 10, the abscissa has been scaled to v_{∞} = 2700 km s^{−1}, the “observed” 1D value applied by Marcolino et al. (2013). 

Open with DEXTER 
Once again we apply ϵ_{L} = 0, and compare the corresponding line profiles with linestrength parameter, κ_{0} = 1, for equatoron and poleon observers, with the line profiles obtained from the MHD simulations (see Fig. 11). For clarification, Fig. B.1 additionally displays all line profiles with their emission and absorption parts. The differences between the profilesets can be explained as follows.
For poleon observers. With increasing contribution from the upflow component, the emission peak becomes broader, because the emitting volume at intermediate to high velocities increases. Simultaneously, the cooled downflow component with only low absolute velocities decreases, resulting in a lower emission peak near the line centre. An exception is model (i), for which the emission peak at low frequency shifts is lowest, because the emission of the upflow component near the star (with high densities and low velocities) is missing. When compared to the line profiles from the MHD simulations, the best result is obtained for model (ii), that is for the statistical description of upflow and downflow component in the closedfield region. Even for this model, however, we only get a relatively poor match with the MHD profiles. Since our ADM models cover a large range of combinations of upflow and downflow component (including the most extreme cases of a pure upflow and a pure downflow), this finding suggests that the outer wind region is inadequately modelled. Indeed, the major differences of the line profiles can be explained (at least qualitatively) by the different description of the outer wind: (1) In the ADM models, the emission peaks close to line centre (i.e. at x_{OBS} ≈ 0, with corresponding resonance zones at projected velocities n · u ≈ 0) are underestimated compared to the MHD model, since within all ADM models also the mass flow in the outer wind is adopted to follow closed magnetic field lines. This assumption becomes problematic in equatorial regions, since here the ADM wind flows almost perpendicular to the plane, whereas it is almost radial in the MHD case. Consequently, when viewed pole on, only large projected velocities are present in the corresponding area of the equatorial plane, where the latter creates a large part of lowvelocity emission in the MHD model. This part is now missing in the ADM models, and the emitting area is almost limited to the downflow component (with generally low projected velocities). Thus, the profile becomes shallower than in the MHD case. (2) Within the blue absorption trough, the absorption column in front of the star is slightly decreased, because the velocity vectors are once again following the magnetic field lines, and do not perfectly match the MHD simulations. The differences of the line profiles can thus be explained by the different description of the outer wind region alone.
For equatoron observers. The emission peak of the ADM models becomes stronger and shifted to the blue side with increasing contribution of the upflow component. Additionally, the absorption part on the blue side increases, while it decreases on the red side. This behaviour is readily explained: as the upflow contribution in front of the star (with projected velocities directed towards the observer, thus affecting the blue side of the profile) grows, the downwind contribution (with projected velocities directed away from the observer, and affecting the red side of the profile) is diminished. Consequently, the absorption in front of the star increases on the blue side, and decreases on the red one. Again, when compared to the MHD model, none of the obtained profiles provides a good agreement. While model (iv) reproduces the absorption part on the blue side relatively well, the redsided absorption part is highly underestimated. On the other hand, a better model for the redsided absorption (e.g. model iii) underestimates the absorption part on the blue side. In fact, it is not possible to simultaneously model the blue and red absorption by only tuning the composition of the closedfield region, suggesting that (at least) the outer wind needs to be treated differently. For instance, assuming a radial outflow in the equatorial plane of the outer wind region would increase the bluesided absorption (and emission), while preserving the rather good behaviour of model (iii) on the red side.
Taking all this evidence together, we conclude that the (present) ADM model needs to be improved for the modelling of UV resonance lines, at least in the outer wind. Such a reformulation then needs to include a consistent description of the actual velocity and density stratification, accounting for the delicate interplay between Bfield and wind.
7. Summary and conclusions
In this paper, we introduced a newly developed 3D FVM code, which solves the equation of radiative transfer for continuumand linescattering problems (the latter approximated by a twolevelatom in the present version). An observer’sframe formulation allows us to consider arbitrary velocity fields (and density structures).
Within the ALI scheme, the code iterates the source functions to convergence, using a nonlocal approximateΛoperator, and extrapolating subsequent iterates by the Ngformalism. For the most challenging problems of optically thick, scattering dominated atmospheres, we obtained a satisfying convergence behaviour, with relative corrections between subsequent iterates of less than 10^{−3} within 20 iteration steps. Due to this convergence behaviour, we were able to analyse the performance of the 3D FVM for such optically thick, scattering dominated atmospheres.
A comparison of spherically symmetric problems calculated with our 3D code and an accurate 1D solver shows that the FVM requires a relatively high spatial resolution. Continuum errors are of the order of ≲20% for marginally optically thick atmospheres, that is for typical Ostar conditions, however increase rapidly for larger optical depths, due to the firstorder scheme. These errors can only be reduced by applying either a higher grid resolution, or when using more accurate solution schemes (e.g. the SC method with appropriate interpolations). To analyse rapidly expanding stellar winds with large continuum optical depths, the development of a 3D SC method is planned for the future in our group.
The line transfer, on the other hand, performs much better, with relative errors less than 25% even for strong lines. The resulting profiles as obtained from a postprocessing LC method are reasonably accurate.
Due to significant numerical diffusion, intrinsic to the FVM, we found a minimum error (for optically thin continua and weak lines) of roughly 10%. Additionally, any symmetry of a considered problem is broken, due to the distinct behaviour of numerical diffusion for different ray directions and in different regions of the atmosphere. Numerical diffusion errors, however, could be minimized by increasing the grid resolution in the outer parts of the atmosphere. Accounting for the sound reproduction of line profiles for spherically symmetric models, we conclude that our code is ready to be used also for arbitrary 3D atmospheric structures, at least if the continuum displays an only moderate optical depth.
As a first application to continuumscattering problems, we estimated the radiation temperatures of windablation models, focusing on the transition region between a linedriven wind and the optically thick circumstellar disc (as present during the late phases of massive star formation in accreting high mass stars). We found a reduction of radiation temperatures by only few percent, which indicates that the ionization stages in this region are (almost) the same as in the wind. Thus, a linedistribution formalism with the same set of linestrength parameters as used in the wind can be applied to obtain the line acceleration that finally ablates the disc. Because of the numerical inaccuracies of the FVM, our findings must be taken with caution, and possibly rechecked with more elaborate methods or a much higher grid resolution. To analyse the complete evolution of optically thick circumstellar discs, the impact of continuum scattering on latitudinal forces still needs to be investigated, and is left to future studies.
As a benchmark for our code regarding the line transfer in nonspherical models, we considered the same MHD simulations of dynamical magnetospheres as used by Marcolino et al. (2013), and compared the resulting UV resonanceline profiles to those obtained from their 3DSEI analysis. The profiles as viewed both polaron and equatoron are in excellent agreement, indicating that our 3D FVM performs well also for such atmospheric models. We additionally applied the analytic dynamical magnetosphere framework (ADM, Owocki et al. 2016), and modelled the corresponding atmospheric structure by adopting four different descriptions of the closedfield regions. A comparison between the obtained line profiles and those for the MHD simulations from above showed significant differences. These were explained by the (somewhat insuffcient) description of the outer wind region within the (present) ADM formulation, primarily in the equatorial plane. An improvement of the underlying assumptions is planned for future work.
The projected ∆τsteps represent the optical depth of the cell for a given direction, and are easily obtained from Eq. (14) and the definition of the coeffcients c_{ijk} … e_{ijk}, e.g. c_{ijk} =: 1/(1 + ∆τ_{x}).
We emphasize that for the line transfer, the ratio of the photon destruction probability, ϵ_{L}, to the photon escape probability is less than unity for all our models, indicating that the linetransfer is, in principal, nonlocal. However, for the considered spherically symmetric problems, no multiple resonances arise, and the line is formed within a single, welllocalized resonance region.
Strictly speaking, Eq. (35) holds only when the strongest line is optically thick. See Kee (2015), Chap. 2, for a complete derivation and discussion.
Acknowledgments
We thank our anonymous referee for valuable comments and suggestions. Many thanks to V. Petit and S. Owocki for fruitful discussions about the ADM. LH gratefully acknowledges support from the German Research Foundation, DFG, under grant PU 117/91. NDK acknowledges support from the German DFG under grant KU 2849/31, which funds the Emmy Noether research group on “Accretion Flows and Feedback in Realistic Models of Massive Star Formation”. JOS acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie SklodowskaCurie grant agreement no. 656725.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Adam, J. 1990, A&A, 240, 541 [NASA ADS] [Google Scholar]
 Alecian, E., Wade, G. A., Catala, C., et al. 2013, MNRAS, 429, 1001 [NASA ADS] [CrossRef] [Google Scholar]
 Amarsi, A. M., Asplund, M., Collet, R., & Leenaarts, J. 2016, MNRAS, 455, 3735 [NASA ADS] [CrossRef] [Google Scholar]
 Belczynski, K., Holz, D. E., Bulik, T., & O’Shaughnessy, R. 2016, Nature, 534, 512 [NASA ADS] [CrossRef] [Google Scholar]
 Busche, J. R., & Hillier, D. J. 2005, AJ, 129, 454 [NASA ADS] [CrossRef] [Google Scholar]
 Cannon, C. J. 1973, ApJ, 185, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [NASA ADS] [CrossRef] [Google Scholar]
 Cherepashchuk, A. M. 1976, Sov. Astron. Lett., 2, 138 [NASA ADS] [Google Scholar]
 Cranmer, S. R., & Owocki, S. P. 1996, ApJ, 462, 469 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., Langer, N., Izzard, R. G., Sana, H., & de Koter, A. 2013, ApJ, 764, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Dufton, P. L., Dunstall, P. R., Evans, C. J., et al. 2011, ApJ, 743, L22 [NASA ADS] [CrossRef] [Google Scholar]
 Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. 1997, A&A, 324, 161 [Google Scholar]
 Gayley, K. G. 1995, ApJ, 454, 410 [NASA ADS] [CrossRef] [Google Scholar]
 Georgiev, L. N., Hillier, D. J., & Zsargó, J. 2006, A&A, 458, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gräfener, G., Koesterke, L., & Hamann, W.R. 2002, A&A, 387, 244 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R. 1980, A&A, 84, 342 [NASA ADS] [Google Scholar]
 Hamann, W.R. 1981, A&A, 93, 353 [NASA ADS] [Google Scholar]
 Hauschildt, P. H. 1992, J. Quant. Spectr. Rad. Transf., 47, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Howarth, I. D., Walborn, N. R., Lennon, D. J., et al. 2007, MNRAS, 381, 433 [NASA ADS] [CrossRef] [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013a, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., et al. 2013b, ASP Conf. Ser., 474, 66 [NASA ADS] [Google Scholar]
 Jones, H. P. 1973, ApJ, 185, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, H. P., & Skumanich, A. 1973, ApJ, 185, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Kee, N. D. 2015, PhD Thesis, University of Delaware, USA [Google Scholar]
 Kee, N. D., Owocki, S., & Sundqvist, J. O. 2016, MNRAS, 458, 2323 [NASA ADS] [CrossRef] [Google Scholar]
 Kuiper, R., Turner, N. J., & Yorke, H. W. 2016, ApJ, 832, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectr. Rad. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Kunasz, P. B., & Olson, G. L. 1988, J. Quant. Spectr. Rad. Transf., 39, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lamers, H. J. G. L. M., CerrutiSola, M., & Perinotto, M. 1987, ApJ, 314, 726 [NASA ADS] [CrossRef] [Google Scholar]
 Leenaarts, J., & Carlsson, M. 2009, ASP Conf. Ser., 415, 87 [NASA ADS] [Google Scholar]
 Lobel, A., & Blomme, R. 2008, ApJ, 678, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1983, ApJ, 274, 372 [NASA ADS] [CrossRef] [Google Scholar]
 Maeder, A. 1999, A&A, 347, 185 [NASA ADS] [Google Scholar]
 Maeder, A., & Meynet, G. 2000, A&A, 361, 159 [NASA ADS] [Google Scholar]
 Marcolino, W. L. F., Bouret, J.C., Sundqvist, J. O., et al. 2013, MNRAS, 431, 2253 [NASA ADS] [CrossRef] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edition (San Francisco: W. H. Freeman and Co) [Google Scholar]
 Mullan, D. J. 1984, ApJ, 283, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Ng, K.C. 1974, J. Chem. Phys., 61, 2680 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, G. L., & Kunasz, P. B. 1987, J. Quant. Spectr. Rad. Transf., 38, 325 [NASA ADS] [CrossRef] [Google Scholar]
 Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spectr. Rad. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & udDoula, A. 2004, ApJ, 600, 1004 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., udDoula, A., Sundqvist, J. O., et al. 2016, MNRAS, 462, 3830 [NASA ADS] [CrossRef] [Google Scholar]
 Patankar, S. V. 1980, Numerical Heat Transfer and Fluid Flow (Hemisphere Publishing Corporation), Series on Computational Methods in Mechanics and Thermal Science [CrossRef] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petit, V., Owocki, S. P., Wade, G. A., et al. 2013, MNRAS, 429, 398 [NASA ADS] [CrossRef] [Google Scholar]
 Petit, V., Keszthelyi, Z., MacInnis, R., et al. 2017, MNRAS, 466, 1052 [NASA ADS] [CrossRef] [Google Scholar]
 Pittard, J. M. 2009, MNRAS, 396, 1743 [NASA ADS] [CrossRef] [Google Scholar]
 Prilutskii, O. F., & Usov, V. V. 1976, Sov. Astron., 20, 2 [NASA ADS] [Google Scholar]
 Puls, J. 1991, A&A, 248, 581 [NASA ADS] [Google Scholar]
 Puls, J., & Herrero, A. 1988, A&A, 204, 219 [NASA ADS] [Google Scholar]
 Puls, J., Owocki, S. P., & Fullerton, A. W. 1993, A&A, 279, 457 [NASA ADS] [Google Scholar]
 Puls, J., Springmann, U., & Lennon, M. 2000, A&AS, 141, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Vink, J. S., & Najarro, F. 2008, A&ARv, 16, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RamírezAgudelo, O. H., SimónDíaz, S., Sana, H., et al. 2013, A&A, 560, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rivero González, J. G., Puls, J., Najarro, F., & Brott, I. 2012, A&A, 537, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sana, H., & Evans, C. J. 2011, IAU Symp., 272, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Sana, H., Rauw, G., & Gosset, E. 2001, A&A, 370, 121 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seelmann, A. M., Hauschildt, P. H., & Baron, E. 2010, A&A, 522, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stenholm, L. G., Stoerzer, H., & Wehrse, R. 1991, J. Quant. Spectr. Rad. Transf., 45, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Stevens, I. R., Blondin, J. M., & Pollock, A. M. T. 1992, ApJ, 386, 265 [NASA ADS] [CrossRef] [Google Scholar]
 Sundqvist, J. O., udDoula, A., Owocki, S. P., et al. 2012, MNRAS, 423, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Tessem, T. 2013, Master’s Thesis, University of Bergen, Norway [Google Scholar]
 Townsend, R. H. D., & Owocki, S. P. 2005, MNRAS, 357, 251 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Trujillo Bueno, J., & Fabiani Bendicho, P. 1995, ApJ, 455, 646 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., & Owocki, S. P. 2002, ApJ, 576, 413 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., Owocki, S. P., & Townsend, R. H. D. 2008, MNRAS, 385, 97 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., Sundqvist, J. O., Owocki, S. P., Petit, V., & Townsend, R. H. D. 2013, MNRAS, 428, 2723 [NASA ADS] [CrossRef] [Google Scholar]
 udDoula, A., Owocki, S., Townsend, R., Petit, V., & Cohen, D. 2014, MNRAS, 441, 3600 [NASA ADS] [CrossRef] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Wade, G. A., Howarth, I. D., Townsend, R. H. D., et al. 2011, MNRAS, 416, 3160 [NASA ADS] [CrossRef] [Google Scholar]
 Wade, G. A., et al. (MiMeS Collaboration) 2012, ASP Conf. Ser., 464, 405 [NASA ADS] [Google Scholar]
 Weber, J. A., Pauldrach, A. W. A., Knogl, J. S., & Hoffmann, T. L. 2013, A&A, 555, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsargó, J., Hillier, D. J., & Georgiev, L. N. 2006, A&A, 447, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zsargó, J., Hillier, D. J., & Georgiev, L. N. 2008, A&A, 478, 543 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: The discretized EQRT within the FVM
Following the ideas of Patankar (1980) and Adam (1990), we here describe the discretization scheme of the timeindependent equation of radiative transfer, Eq. (1). At each grid point, Eq. (1) is integrated over a finite control volume (see Fig. A.1). Applying Gauss’s theorem, we obtain:
Here and in the following, we omit the notation for the explicit frequency dependence. The lefthand side of Eq. (A.1) describes the intensity propagating into and out of the controlvolume surfaces, and the righthand side corresponds to the gridcell contribution from sources and sinks. Assuming that the variables at the grid points are appropriate mean values within the corresponding control volume, the righthand side is easily integrated, and gives for grid point (i, j, k),
Since we are using Cartesian coordinates, the integral on the lefthand side can be readily calculated:
Again, assuming that the intensities at the midpoints of the controlvolume surfaces are representative averages of the corresponding surfaces, we obtain:
Since the controlvolume coordinates are positioned at the midpoints of the grid coordinates, we substitute:
Finally, we use the upwind approximation to replace the (unknown) intensities at the controlvolume surfaces:
Fig. A.1. Geometry used within the controlvolume approach: The discretized 3D spatial grid is shown in blue. The dashed lines indicate the control volume, corresponding to a gridpoint (i, j, k). The controlvolume surfaces are located at the centre between the gridpoint coordinates. 

Open with DEXTER 
Combining Eqs. (A.1), (A.2), (A.4), (A.5)–(A.7), and the definitions of α, β, γ, the discretized EQRT finally reads:
where we already have separated the continuum and line contribution of the opacity and source function. Collecting terms, and solving for I_{i}jk leads to:
Appendix B: UV line profiles for different ADM models
Fig. B.1. Left panel: as Fig. 9, but for the corresponding ADM model structures. To clarify the definition of the apex radius, r_{m}, we have displayed a specific value, r_{m} = 2R_{∗}, as a red arrow, where this value corresponds to all points located on the red magneticfield line. In the closedfield region (inside r_{m} = R_{A}, displayed by a thick line), the models contain, from top to bottom: (i) the cooleddownflow component alone. (ii) A statistical approach for the downflow and upflow component. (iii) Alternating flux tubes with cooleddownflow and windupflow component. (iv) The windupflow component alone. Middle panel: as Fig. 11, for poleon observers, and for the different ADM models (i)–(iv). The dashed and dotted lines display the emission part and the absorption part of the line profiles, respectively. Right panel: as middle panel, but for equatoron observers. 

Open with DEXTER 
All Tables
Mean relative errors of the mean intensity for a zeroopacity model, using either Gauss–Legendre quadrature, or the trapezoidal rule with integration nodes from Lobel & Blomme (2008).
Input parameters for the spherically symmetric models used for our test calculations.
Mean and maximum relative errors of our 3D solution scheme, when compared to an “accurate” 1D solution.
All Figures
Fig. 1. Convergence behaviour of the continuum transfer with ϵ_{c} = 10^{−6} (left), and of the line transfer with ϵ_{L} = 10^{−6} (right), for different acceleration techniques applied to the RT, and for a spherical model atmosphere. The optical depth scale varies according to k_{c} = 10, k_{L} = 10 (top) and k_{c} = 100, k_{L} = 10^{5} (bottom). See text. 

Open with DEXTER  
In the text 
Fig. 2. Boundary conditions for two different points, r_{1}, r_{2}, and different directions, n_{1}, n_{2}. Green: a ray originating from the stellar photosphere. To calculate the intensity at r_{1} in direction n_{1}, the intensities at the neighbouring grid points, and , need to be known. A boundary condition is required for grid point , while the intensity at results within the “normal” RTscheme. Red: a ray originating from outside the photosphere. For the grid point r_{2}, a boundary condition has to be specified at the corresponding phantom point, . 

Open with DEXTER  
In the text 
Fig. 3. Top: specific intensities for direction n= (1, 0, 1) in the x–zplane, for a zeroopacity model. Overplotted is the spatial grid with a typical size of N_{x} = N_{y} = N_{z} = 133 grid points. The thick lines indicate the theoretical boundary of the beam. Middle: projected intensity profile through the area indicated by a straight line in the top panel, covering all rays with direction n and impact parameter p = ±2R_{∗}, and theoretical intensity profile denoted by dotted lines. The central distance from the origin is 3R_{∗}. Bottom: specific intensity along n, compared to the theoretical value indicated by the dotted line. 

Open with DEXTER  
In the text 
Fig. 4. Contours of the mean intensity for a zeroopacity model, normalized by its theoretical value, on spherical surfaces, r = 1.1R_{∗} (left) and r = 3R_{∗} (right). The line of sight is along the vector n= (1, 1, 1). N_{x} = N_{y} = N_{z} = 133 grid points have been used. 

Open with DEXTER  
In the text 
Fig. 5. Top: mean intensities (scaled by the emitted intensity from the stellar core, I_{c}) for the continuum transport, with ϵ_{c} = 10^{−6}, and k_{c} = (10^{0}, 10^{1}, 10^{2}) corresponding to τ_{r} = (0.17, 1.7, 17.0), from left to right. Bottom: line source functions (scaled by I_{c}) assuming an optically thin continuum, and with ϵ_{L} = 10^{−6}, and k_{L} = (10^{0}, 10^{3}, 10^{5}) from left to right. The red line corresponds to the accurate 1D solution, the blue line is the solution along the xaxis, and the dots (mainly located in between the red and the blue line) correspond to the solution of (arbitrarily) selected grid points (only a subset of all grid points is shown to compress the images). 

Open with DEXTER  
In the text 
Fig. 6. Emergent line profiles for spherically symmetric models with different linestrengths, k_{L} = (10^{0}, 10^{3}, 10^{5}), from left to right. The red line corresponds to the accurate 1D solution, and the black line is the solution corresponding to our 3D calculations. 

Open with DEXTER  
In the text 
Fig. 7. Density structure for the windablation model with τ_{disc} = 1400, in the x–zplane. 

Open with DEXTER  
In the text 
Fig. 8. As Fig. 7, but for the radiation temperature. Additionally, the density contours corresponding to a decrease in line force by f = 10 (solid) and f = 100 (dashed) are displayed. Both contours indicate the transition region from wind to disc (see text). The grey colour corresponds to radiation temperatures less than the colourcoded minimum value. 

Open with DEXTER  
In the text 
Fig. 9. Upper panel: density structure for an example snapshot from the MHD simulations for HD 191612, as performed by Sundqvist et al. (2012). Lower panel: azimuthal average of the MHD simulations. In both figures, the density has been normalized by a typical downflow density, , with Ṁ_{B}= 0 from Table 5, and v_{esc} ≈ 800 km s^{−1} the photospheric escape velocity. The velocity field is displayed by arrows, with the length of the velocity vectors limited to 0.5 v_{esc}. We additionally show the dipole magnetic field of the ADM models used in Sect. 6.2 (solid lines, and thick solid line for R_{A} = 2.7R_{∗}). The corresponding magnetic axis is aligned with the zaxis. The grey colour corresponds to densities outside the range indicated on the right. 

Open with DEXTER  
In the text 
Fig. 10. UV resonanceline profiles for the MHD models, as obtained from our 3Dcode (black) and from the SEI method by Marcolino et al. (2013; red). Two different linestrength parameters, κ_{0} = 0.1 and 1.0, have been used. For convenience, the line profiles for κ_{0} = 0.1 have been shifted vertically by 1.5. The upper and lower panels show the synthetic line profiles for poleon and equatoron observers, respectively. The abscissa has been scaled to v_{∞} = 2700 km s^{−1}, the “observed” 1D value applied by Marcolino et al. (2013). 

Open with DEXTER  
In the text 
Fig. 11. UV resonanceline profiles obtained from our 3D code, for the MHD simulation (black, as Fig. 10), and the different ADM models, (i)–(iv) (see text). The linestrength parameter has been set to κ_{0} = 1. The upper and lower panels show the synthetic profiles for poleon and equatoron view, respectively. As in Fig. 10, the abscissa has been scaled to v_{∞} = 2700 km s^{−1}, the “observed” 1D value applied by Marcolino et al. (2013). 

Open with DEXTER  
In the text 
Fig. A.1. Geometry used within the controlvolume approach: The discretized 3D spatial grid is shown in blue. The dashed lines indicate the control volume, corresponding to a gridpoint (i, j, k). The controlvolume surfaces are located at the centre between the gridpoint coordinates. 

Open with DEXTER  
In the text 
Fig. B.1. Left panel: as Fig. 9, but for the corresponding ADM model structures. To clarify the definition of the apex radius, r_{m}, we have displayed a specific value, r_{m} = 2R_{∗}, as a red arrow, where this value corresponds to all points located on the red magneticfield line. In the closedfield region (inside r_{m} = R_{A}, displayed by a thick line), the models contain, from top to bottom: (i) the cooleddownflow component alone. (ii) A statistical approach for the downflow and upflow component. (iii) Alternating flux tubes with cooleddownflow and windupflow component. (iv) The windupflow component alone. Middle panel: as Fig. 11, for poleon observers, and for the different ADM models (i)–(iv). The dashed and dotted lines display the emission part and the absorption part of the line profiles, respectively. Right panel: as middle panel, but for equatoron observers. 

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