Self-consistent modelling of line-driven hot-star winds with Monte Carlo radiation hydrodynamics

Radiative pressure exerted by line interactions is a prominent driver of outflows in astrophysical systems, being at work in the outflows emerging from hot stars or from the accretion discs of cataclysmic variables, massive young stars and active galactic nuclei. In this work, a new radiation hydrodynamical approach to model line-driven hot-star winds is presented. By coupling a Monte Carlo radiative transfer scheme with a finite-volume fluid dynamical method, line-driven mass outflows may be modelled self-consistently, benefiting from the advantages of Monte Carlo techniques in treating multi-line effects, such as multiple scatterings, and in dealing with arbitrary multidimensional configurations. In this work, we introduce our approach in detail by highlighting the key numerical techniques and verifying their operation in a number of simplified applications, specifically in a series of self-consistent, one-dimensional, Sobolev-type, hot-star wind calculations. The utility and accuracy of our approach is demonstrated by comparing the obtained results with the predictions of various formulations of the so-called CAK theory and by confronting the calculations with modern sophisticated techniques of predicting the wind structure. Using these calculations, we also point out some useful diagnostic capabilities our approach provides. Finally we discuss some of the current limitations of our method, some possible extensions and potential future applications.


INTRODUCTION
Predicting massive star evolution is dramatically complicated by the presence of powerful winds, constituting an important massloss mechanism (see, for example, overview in Kudritzki & Puls 2000). The evolutionary timescales of the star and its luminosity are, for example, significantly affected by this continuous mass loss. In addition to influencing the stars themselves, these winds also affect their environments by injecting energy and momentum into the interstellar medium and supplying chemically enriched material. For a detailed study of all these processes and effects a firm understanding of the mechanisms driving stellar winds is required.
In the case of hot O and B stars, the main driver for the mass outflow has been identified as the momentum transfer mediated by a multitude of atomic line interactions of the radiation field in the wind material (Lucy & Solomon 1970;Castor et al. 1975). Over the years, this picture has been significantly refined and a "standard model" (see, for example, review by Puls et al. 2008) for line-driven winds has emerged, with important contributions from Castor et al. (1975); Pauldrach et al. (1986); Kudritzki et al. (1989). With the inclusion of multi-line effects, in particular by ⋆ unoebauer@mpa-garching.mpg.de Vink et al. (2000), this model has become very successful in describing hot-star winds and predicting their properties. However, observational diagnostics have identified a number of puzzles that challenge this "standard model" (see again Puls et al. 2008). Among these, the so-called clumping problem takes a prominent role. Contrary to the assumptions of the "standard model", hotstar winds seem to exhibit strong temporal and spatial variability. Intense effort has been invested into understanding the properties and the nature of these clumpy outflows, both from an observational and theoretical viewpoint. In light of these efforts, we have developed a new radiation hydrodynamical approach to numerically model line-driven mass outflows self-consistently. Due to its reliance on Monte Carlo radiative transfer techniques, this approach should be well-suited for time-dependent and multidimensional self-consistent studies of hot-star winds. In previous works, a Monte Carlo radiative transfer scheme has already been successfully coupled with a fluid dynamical method, resulting in the radiation hydrodynamical approach MCRH (Noebauer et al. 2012;Noebauer 2014). In this study we present extensions to MCRH and demonstrate the utility of this method for the investigation of linedriven winds. This work has been partially carried out during the PhD project of Noebauer (2014).
The presentation of our approach is structured as follows: in c 0000 RAS Section 2, we briefly recapitulate some basic concepts of hot-star winds and outline different solution strategies to address this problem. We also highlight the advantages of Monte Carlo radiative transfer techniques for the study of line-driven mass outflows. We continue with a detailed description of the numerical techniques relevant for our approach in Section 3. In Section 4, we present a series of simple wind simulations to test our approach and in Section 5, we confront our method with an alternative, modern and sophisticated technique of predicting the wind structure. The obtained results and potential future improvements of our method are discussed in detail in Section 6.

Solving the stationary hot-star wind problem
Numerous techniques and approaches have been used to solve for the structure of hot-star line-driven winds and derive their principle properties. Many of these share the common feature of addressing the problem in a "radiation-hydrostatics" fashion. By assuming that the mass outflow is smooth, spherically-symmetric and in a stationary state, the wind structure may be derived from global energy conservation considerations (e.g. Abbott & Lucy 1985, Lucy & Abbott 1993, Sim 2004 or by solving the momentum equation alone (e.g. Castor et al. 1975, Pauldrach et al. 1986, Kudritzki et al. 1989). If included, inhomogeneities in the mass outflow, i.e. clumps, are typically incorporated in a parameterised fashion, for example, by introducing a clumping factor (e.g. Schmutz 1997;Hamann & Koesterke 1998). In the simplified, stationary situation, the momentum equation reduces to u du dr + 1 ρ dP dr + GM⋆ r 2 = g rad , after introducing the wind velocity u, its density ρ, the thermodynamic pressure P , the stellar mass M⋆ and the radiative acceleration g rad . Solving this equation is still very challenging, since the bulk of the radiative acceleration derives from the momentum transfer mediated by interactions of the radiation field emitted by the central star with a large number of atomic line transitions. The contribution due to Thomson scattering may be taken into account by reducing the stellar mass in the gravity term by a factor (1−Γe), using the Eddington factor with respect to electron scattering, Γe. Various approaches have been used to calculate the linedriving force. We briefly sketch the basic principle of the so-called Castor-Abbott-Klein (CAK) theory (after Castor et al. 1975; see in particular Abbott 1982, Pauldrach et al. 1986, Kudritzki et al. 1989 for improvements of the original approach) since it will be used for comparison in this work. Castor et al. (1975) found an approximate analytic expression for the line-driving force. In particular, the linedriving force is expressed as a multiple of the electron-scattering contribution and this so-called force-multiplier is approximated by a power-law parameterisation. Given this form of the radiative acceleration, an analytic solution to the momentum equation may be identified. By performing a critical point analysis, the wind velocity is found to follow a β-type law with β = 0.5 in the original CAK case of Castor et al. (1975), in which the central star is assumed to radiate like a point source and with β ≈ 0.8 if the finite extent of the star is taken into account (see Pauldrach et al. 1986; Kudritzki et al. 1989). Throughout this work, we reference this improved formulation of the CAK approach as MCAK ("Modified-CAK"). The mass-loss rate of the wind scales asṀ In addition to the stellar parameters (luminosity L⋆, photospheric radius R⋆ and central mass M⋆), the parameters of the power-law fit, k and α, appear. The functional form of the coefficient C(α) differs slightly between the CAK and MCAK approaches. More details about C(α) and the CAK/MCAK methods, as used in this work, are provided in Appendix A. A significant advantage of these radiation-hydrostatic techniques lies in the analytic expressions for the wind structure they either provide or in the possibility to efficiently solve the problem numerically (see, e.g., Abbott & Lucy 1985;Pauldrach et al. 1994;Vink et al. 1999;Puls et al. 2005), even if more sophisticated treatments of the line-driving force are included. In both cases, wind models may be constructed for entire grids of different stellar parameters (e.g. Abbott 1982;Vink et al. 2000;Puls et al. 2005;Muijres et al. 2012). A drawback lies in the omission of a selfconsistent treatment of temporal and spatial variations in the wind flow.

A dynamical approach
To account for, follow and solve the temporal and spatial variability of stellar winds, the full radiation hydrodynamical problem the radiatively-driven mass outflow constitutes has to be addressed (e.g. Owocki et al. 1988;Feldmeier 1995;Feldmeier et al. 1997;Dessart & Owocki 2003;Dessart 2004;Dessart & Owocki 2005). In such an approach, a numerical solution to the full radiation hydrodynamical equations describing the conservation of mass, momentum and energy has to be formulated. In contrast to a pure hydrodynamical problem, the energy and momentum transfer between the radiation field and the fluid material are included. These terms are found by solving the co-evolution of the radiation field in parallel. However, due to the complexity of the full radiation hydrodynamical problem, a simplified treatment of the radiation field and radiative transfer is typically adopted (c.f. Owocki et al. 1988;Feldmeier 1995;Dessart & Owocki 2003). The radiation hydrodynamical equations are given in Appendix B in the form in which they are used in this study.
In this work we propose a new numerical approach to perform such dynamical self-consistent simulations of hot-star winds. The distinguishing feature of our method is the incorporation of Monte Carlo radiative transfer techniques. In particular, we build upon the general-purpose Monte Carlo radiation hydrodynamics scheme MCRH, developed by Noebauer et al. (2012);Noebauer (2014) and adapt it to the line-driven wind environment. Other Monte Carlo-based radiation hydrodynamical techniques have been developed by Nayakshin et al. (2009), Acreman et al. (2010, Haworth & Harries (2012) and most recently by Roth & Kasen (2015), who place a particular emphasis on the incorporation of implicit Monte Carlo methods (see Fleck & Cummings 1971), and by Harries (2015).

Advantages of Monte Carlo techniques
Monte Carlo techniques have become established as a competitive and successful numerical approach to solve radiative transfer problems due to the increasing availability of computational resources and the more and more levels of sophistication added to the formalism (see specifically, Abbott & Lucy 1985, Lucy 1999a, Lucy 2002, Lucy 2003, Carciofi & Bjorkman 2006. The appeal of using Monte Carlo techniques derives from their advantages in treating interaction physics (see, for example, Pozdnyakov et al. 1983 for Monte Carlo calculations of comptonisation or Kasen et al. 2006 andSim 2009, for presentations of fully-fledged three-dimensional Monte Carlo codes to solve multi-frequency radiative transfer in Type Ia supernova ejecta) and from the ease with which arbitrary geometrical configurations are treated (see, for example, Camps et al. 2013, describing a Monte Carlo radiative transfer scheme on a Voronoi mesh). These advantages all derive from the microscopic viewpoint that Monte Carlo techniques take by solving the propagation history of a large number of representative photons probabilistically. They are very relevant for the study of hot-star winds, in which a multitude of atomic line interactions are the main driver, multiple scatterings frequently occur and a clumpy irregular structure of the mass outflow is expected. Thus, Monte Carlo techniques have been frequently used to address the radiative transfer problem in hot-star winds, for example by Abbott & Lucy (1985), Lucy & Abbott (1993), Schmutz (1997) Vink et al. (1999, Vink et al. (2000), Sim (2004), Müller & Vink (2008), Muijres et al. (2012), Šurlan et al. (2012).
In this work, we aim at exploiting the advantages of Monte Carlo techniques for solving the line-driving problem in a dynamical and self-consistent approach by using a Monte Carlo radiation hydrodynamical technique. This strategy distinguishes itself form previous Monte Carlo-based investigations by its treatment of the coupled co-evolution of the radiation-wind state. Moreover, the radiative acceleration is determined locally, which sets us apart from some earlier Monte Carlo-based studies, such as Abbott & Lucy (1985), Lucy & Abbott (1993), Vink et al. (1999) and Sim (2004), in which the momentum transfer onto the wind material was derived from global considerations without insisting on local consistency (however, see Müller & Vink 2008, for developments towards a consistent local force balance in Monte Carlo-based approaches).

NUMERICAL METHODS
In order to adequately address the radiation-matter coupling problem in the stellar wind environment, we have extended our numerical framework whose underlying methodology and first implementation has been described in Noebauer et al. (2012). In the following, the basic principles of our approach are briefly recalled and the new extensions described in detail.

Basic operator-split approach
In our approach, a simple operator-splitting scheme is employed, coupling a time-dependent Monte Carlo radiative transfer simulation with a finite-volume fluid dynamical calculation. In the latter, the piecewise parabolic method (PPM, Colella & Woodward 1984) is used to construct a series of Riemann problems at the interfaces of the computational domain, which are solved with a standard Riemann solver (HLLC, see Toro et al. 1994;Batten et al. 1997) to determine the mass, energy and momentum flux through the interfaces and eventually the evolution of these conserved quantities in the cells. In the following radiative transfer step, a Monte Carlo simulation is performed to determine the evolution of the radiation field and reconstruct the energy and momentum transfer between the radiation field and the fluid material. In the final step of the splitting approach, the momentum and energy of the fluid are altered in accordance with these transfer terms.
This basic outline of the numerical scheme remains valid for the stellar wind applications. Alterations and extensions of the scheme only concern the inclusion of additional physical effects or adopted simplifications that are relevant for an adequate treatment of the stellar wind environment. In particular, the gravitational field originating from the central star is taken into account, an isothermal treatment of fluid dynamics is adopted and the radiative transfer scheme is generalised to incorporate frequency-dependent resonant line interactions. These alterations are systematically introduced and presented in detail in the following sections.

Isothermal hydrodynamics and gravity
In this work, we reduce the complexity of the hydrodynamical calculations in our modelling approach of hot-star winds by adopting an isothermal approximation. This simplification, which has been used in numerous investigations of line-driven hot-star winds (e.g. Abbott 1982;Abbott & Lucy 1985;Owocki et al. 1988;Vink et al. 1999Vink et al. , 2000Sim 2004), is justified by the continuous irradiation of the wind material by the central star, which together with a characteristic cooling time scale that is faster than the flow time scale sustains the wind material roughly at the effective temperature of the star (Klein & Castor 1978). In the isothermal approximation, only the mass and momentum conservation equations have to be addressed in the fluid dynamical calculation. With the equation of state of an isothermal flow directly connecting the fluid density and pressure a solution of the energy conservation equation is not required. In the case of an ideal gas with constant temperature Tiso and mean molecular weight µ, the isothermal sound speed aiso is aiso = kBTiso µ .
We modify the fluid dynamical solution procedure in the isothermal version of MCRH by reconstructing only the primitive variables density and the velocity. Using the isothermal equation of state, the pressure may be directly determined and the Riemann problems at the cell interfaces solved. Here, we use a particular version of the HLLC Riemann solver, adopted from the ATHENA code 1 (Stone et al. 2008), which is tailored to the particular form of the simple waves in isothermal flows.
An important contribution to the momentum balance in hotstar winds derives from the gravitational pull exerted by the central object Following Colella & Woodward (1984), we incorporate the effect of the star's gravitational pull into the isothermal fluid dynamical calculation by accounting for the additional momentum density contribution ∆p n = 1 2 ∆t(ρ n g n + ρ n+1 g n+1 ), in the final determination of the new fluid dynamical state at the end of the splitting step. Here, the superscripts n and n + 1 mark quantities at the beginning and the end of a time step with length ∆t = t n+1 − t n .

Monte Carlo principles
At the heart of the Monte Carlo radiative transfer machinery lies the discretization of the radiation field into a large number of quanta, so called packets, which are launched into the computational domain. In a series of random number experiments, the packets describe radiation-matter interactions and the temporal evolution of the radiation field. The basic layout of our Monte Carlo method has already been presented in detail in Noebauer et al. (2012). Here, we only recall some important aspects and dedicate the bulk of the presentation to the additions related to the frequency-dependent calculation.

Discretization, initialization and propagation
We adopt the indivisible energy packet scheme of Abbott & Lucy (1985) and Lucy (1999a). Thus, Monte Carlo packets are interpreted primarily as parcels of radiative energy. To model the continuous illumination of the outflowing material by the central star, a number of packets are injected into the computational domain through the lower boundary in each radiative transfer step, Assuming that the star radiates as a black body with effective temperature T eff , new packets, each with energy ε, are launched from the inner boundary (located at R⋆) during a time step of length ∆t. Neglecting any limb-darkening effect, the initial propagation direction of these packets follows Here and in the following, µ denotes the cosine of the directional angle between the trajectory and the symmetry axis and z represents a random number, uniformly drawn from the interval [0, 1]. The frequency assignment of the packets is chosen to reflect the Planck function 2 This frequency assignment is adopted for this current study for simplicity and convenience. Sampling a realistic atmosphere model, however, does not pose a conceptual or implementation challenge (see, e.g., Abbott & Lucy 1985). After determining the initial properties of the packets, they are launched and initiate the propagation process. As detailed in Noebauer et al. (2012), we follow a mixed-frame approach and perform the packet propagation procedure in the lab frame (LF) in which the computational grid is defined. However, all interaction processes are treated in the local co-moving frame (CMF), in which the fluid is at rest and where the material functions take convenient forms. Transformation rules relate quantities between these two frames. These are derived in detail, for example, in Mihalas & Mihalas (1984) and are listed partially in Noebauer et al. (2012).
During their propagation, packets interact whenever they have covered a pathlength equivalent to the optical depth which is determined for every packet in a random number experiment. The procedure with which the nature of the interaction event is determined is presented in the next section. Once the details of the interaction process are established, the packet properties are updated accordingly. We highlight the case of line interactions, which we treat as resonant scatterings in the Sobolev approximation (after Sobolev 1960; see also, for example, Lamers & Cassinelli 1999 for a summary of this approximation and Pauldrach et al. 1986 for a discussion of its applicability to hot-star winds). In this case, the packet is assigned a new propagation direction in the CMF according to the Sobolev escape probability Here, τs denotes the Sobolev optical depth, which will be explicitly introduced in Section 3.7 and which depends on the direction of the photon trajectory. The changes in the LF frequency and energy of the packet follow from the Doppler effect and from energy conservation in the CMF. Assuming that the process of sampling the Sobolev escape probability returned the CMF propagation direction µ a 0 and accounting for the appropriate transformation laws, the resonant line scattering process results in the new packet quantities: Here and in the following, a subscribed "0" denotes CMF quantities and the superscripts b and a distinguish packet properties before and after the interaction process. Moreover, the standard notation of special relativity is adopted, introducing β = u/c and γ = 1/ 1 − β 2 . Above and in the following we keep full account of the relativistic terms and only in the final expression reduce the accuracy to O(u/c), which is appropriate for modelling hot-star winds.
In addition to the physical radiation-matter interactions, packets also experience so-called numerical events due to the spatial and temporal discretization of the computational domain. When a packet intercepts a grid cell boundary, some packet properties, such as the optical depth distance to the next interaction, have to be redetermined to be compatible with the physical conditions in the target cell. At the end of the time step, the trajectories of all remaining packets are interrupted and their properties stored to allow them to resume their propagation during the next simulation cycle.

Optical depth summation
The sole purpose of the packet propagation process lies in the determination of packet interaction histories from which the radiation field characteristics can be reconstructed. For the stellar wind application, this concerns primarily the line-driving force. Crucial for the determination of the packet trajectories is the decision about where and how the packets interact. As already mentioned, packets experience physical radiation-matter interactions after having accumulated an integrated optical depth equal to the random number experiment outcome (11). For the current work, the wealth of possible interaction processes has been restricted to include only frequency-dependent bound-bound processes, which we treat as resonant scatterings in the Sobolev limit (see previous section). All continuum processes are neglected, apart from Thomson scattering, which is incorporated approximately by reducing the mass of the central star as outlined in Section 2.1.
Once a packet interacts, the location of this event in optical depth space has to be translated into a physical position within the computational domain. This procedure is trivial in the grey case, since optical depth and pathlength only differ by the opacity, a constant multiplicative factor in the CMF, but complicated in the presence of frequency-dependent processes. We adopt a simplified version of the technique of Mazzali & Lucy (1993) to locate the line-interaction events packets perform. On its trajectory, a packet propagates freely to the Sobolev point of the next line with which it comes into resonance. Each time such a resonance point is reached, the optical depth is incremented by the full line optical depth of the corresponding transition. The packet undergoes an interaction once the value drawn in (11) is surpassed by the optical depth accumulated. If this occurs during the instantaneous increases at one of the Sobolev points, the packet undergoes a resonant line interaction, otherwise it may leave the current grid cell uninterrupted. This procedure may be easily extended to include additional interaction types, in particular frequency-independent processes, such as Thomson scatterings (see Mazzali & Lucy 1993), but for the current work we have omitted to do so.

Monte Carlo estimators
In reconstructing the radiation field characteristics from the ensemble of packet interaction histories, we follow the volume-averaged estimator approach proposed by Lucy (1999a) and refined by Lucy (2003Lucy ( , 2005. This formalism aims at reducing the statistical fluctuations inherent to the Monte Carlo approach by increasing the number of contributions to the packet census. For the case of frequency-independent processes being the only interaction channel, adequate estimators have already been derived and presented by Noebauer et al. (2012) and similarly by Roth & Kasen (2015).
To determine the radiative acceleration due to spectral line interactions, we consider the momentum transferred in such an event.
Assuming that these interactions occur as resonant scatterings, a packet transfers momentum onto the material. Estimators for the radiation force can be obtained by summing over the transfer terms of all interacting packets. To reduce the statistical fluctuations in these estimators we follow the suggestion of Lucy (1999b) and include all packets that come into resonance with a line and weight their contributions with the corresponding interaction probability given by (1 − e −τs ).
Taking the forward-backward symmetry of the re-emission into account, thus cancelling all terms that are of odd power in µ a 0 , the following estimator for the radiation force [c.f. Equation (B5)] due to line interactions are obtained: Here, the volume of a grid cell, ∆V appears. The superscript b has been dropped and only terms of O(u/c) have been retained.
Using this estimator, the radiation force is calculated and the fluid momentum updated in the final splitting step.

Ionization and level population
The strength of the different line transitions, as encoded in the corresponding values for the Sobolev optical depth, τs, depends on the excitation and ionisation balance in the wind. Thus, the Monte Carlo radiative transfer routine requires a separate calculation step which determines the ionisation and the level population in each cell of the computational grid. In this first study, we follow Abbott & Lucy (1985), Vink et al. (1999) and Sim (2004) and adopt an approximate non-LTE treatment. We stress, however, that nothing in our radiative transfer or hydrodynamics formalism requires these approximations. A full, non-LTE scheme for calculating ionisation and excitation states could be incorporated following the methods described by Lucy (2003, see Section 6.3). Following our simplified strategy, we determine the ionisation balance by applying the modified nebular approximation (see, e.g. Mazzali & Lucy 1993), This expression relates the number densities of two successive ion stages, Nj, Nj+1 with the electron number density Ne. Compared to a pure LTE calculation based on the Saha equation, whose results are denoted by the asterisks, modifications due to the dilution of the radiation field and due to recombination effects are taken into account. As a consequence, the dilution factor W , the ratio of the electron and radiation temperatures, Te and TR, and the fraction ζj of recombination processes going directly to the ground state appear in the above formulation. In the determination of the level population, we again follow Abbott & Lucy (1985) and identify levels with no permitted electromagnetic dipole transitions to lower energy levels as metastable; these levels are assumed to obey Boltzmann statistics with "1" denoting the ground state. For all other levels, the effect of radiative processes on the level population is crudely taken into account by including a dilution factor (see Vink et al. 1999;Sim 2004) For the current work, we do not solve Equation (18) iteratively to determine the ionisation state. Instead, we use a predefined characteristic electron density Ne/W and calculate the ionisation accordingly (see, e.g. Abbott 1982, for a comparable strategy) in all our wind simulations. Also, the radiation temperature is treated in a simple fashion according to the prescription rather than relying on predictions of realistic atmosphere models.
With the excitation and ionisation balance accessible through Equations (18), (19) and (20)  packet trajectory along the direction µ may be calculated according to Here, ν0 denotes the rest frame frequency of the transition and f l its oscillator strength. Also, the statistical weights g l,u associated with the lower and upper energy levels connected by the transition and the electron charge and mass, e and me, appear. Finally, we stress again that the above expressions only depend on the physical conditions at the Sobolev point, an essential feature which is highlighted by the subscribed rs.

Velocity interpolation
As photons propagate through the wind material, their CMF frequency is continuously Doppler-shifted by the varying wind velocity. To reproduce this behaviour in the Monte Carlo simulation, an interpolation scheme is required to reconstruct the evolution of the wind velocity within the computational grid cells. For the current work, we do not expect the formation of shocks but the presence of a smoothly varying velocity field. Thus, in constructing the algorithm, an artificial introduction of discontinuous jumps in the interpolated wind velocity should be avoided: otherwise, packets may skip frequency regions which may be populated by line transitions. We prevent this by devising an interpolation scheme that insists on continuity in the reconstructed wind velocity at the interfaces between grid cells. First, the velocity at the interfaces are determined by linear interpolation between the cell-averaged values provided by the finite-volume fluid dynamical calculation. The interface values are then used in a second linear interpolation step to determine the evolution of the wind velocity within the cell. Figure 1 graphically illustrates this procedure.
Similarly to the spatial considerations, care has be taken when packets reach the end of the time step. In order to ensure consistency in the CMF frequency between successive time steps, the velocity is also interpolated linearly in time between the results of two preceding fluid dynamical splitting steps. Otherwise the changes of the fluid velocity due to fluid dynamics and the radiation-matter coupling could introduce discontinuous jumps in the CMF frequency.

Boundary conditions
In the hot-star wind calculations, the inner boundary of the computational domain is set to the stellar surface, R⋆. So-called ghost cells, which are required for the fluid dynamical reconstruction step at the domain boundaries, are inserted below this radius. The fluid state in these inner ghost cells is found by using a temporally constant hydrostatic density stratification and by extrapolating the velocity field from the innermost region of the domain into the ghost cells using a low-order polynomial. For most wind calculations presented in this work, the hydrostatic profile has been normalised (arbitrarily) to ρ0 = 10 −11 g cm −3 . However, no significant changes in the calculated wind structure was found if this value was modestly varied. In the future, predictions from atmospheric models may be used for this process instead. With the wind density fixed but the possibility of the velocity in the ghost cells to float, the mass flow through the boundary may quickly adjust itself to the conditions in the wind. This situation closely resembles the boundary conditions presented by Owocki et al. (1988) for their time-dependent studies.
As detailed in Section 3.4, the inner edge of the domain constitutes an inflow boundary for the radiation field. Any Monte Carlo packet that is back-scattered through the inner boundary is discarded from the ensuing propagation process.
The outer edge of the computational domain is transparent to the radiation field. All escaping packets are recorded to construct spectra or other diagnostics during the postprocessing steps. Regarding the fluid state, the outer boundary simulates the outflow of wind material. To this end, the values for the fluid variables in the ghost cells are found by extrapolating the wind flow. To minimise the effect of stochastic fluctuations introduced by the Monte Carlo radiative transfer step, a linear regression scheme is used.
This outline of the boundary treatment concludes the description of the numerical techniques which are relevant for this study and which have been added to MCRH. The performance of all these extensions has been verified in a series of test calculations. A detailed description of this procedure and the different test problems used during this process may be found in Noebauer (2014).

HOT-STAR WIND TEST CALCULATIONS
We now aim at demonstrating the capability of our approach to self-consistently solve the line-driving problem in hot-star winds by considering an idealised representation of the problem. Specifically, we test the basic accuracy of our method by exploring whether it produces results compatible with the CAK and MCAK theory when comparable assumptions about the physical effects at work are adopted (smooth flow, Sobolev approximation, etc.). At the same time, we investigate how these simple calculations compare with simulations in which the full details as outlined in Section 3 are incorporated. Consequently, this series of wind simulations, constitutes the basic testing and validation step before we confront our method with the technique of Müller & Vink (2008) in Section 5. Finally, in Section 6, we discuss how some of the simplifications and approximations adopted in the current work may be relaxed and  Table 2. Properties of ζ-Puppis according to Puls et al. (1996). We adopt these stellar parameters in our wind test calculations.
Parameter Value M⋆ 52.5 M ⊙ L⋆ 10 6 L ⊙ T eff 4.2 × 10 4 K a more realistic description of the line-driving problem be achieved within our numerical framework in the future.

Atomic data
In all wind calculations, we assume that the material has solar composition. We include the elements listed in Table 1 and adopt the corresponding abundances of Asplund et al. (2009). Table 1 also identifies the ionisation stages that are taken into account for each element. Information about the atomic structure and the line transitions in these ions are taken from the same database used in the study of colliding winds by Parkin & Sim (2013), which is based on the Kurucz & Bell (1995) atomic data set. To reduce the computational effort in the Monte Carlo radiative transfer steps, we do not account for all line transitions recorded in the data set, but only those with log gf > −6. We have explicitly verified that the inclusion of more weak lines does not affect the outcome of our simulations. When required, the recombination fractions for the included ions, ζi, are adopted from the PYTHON code (Long & Knigge 2002).

Stellar parameters
We perform all wind calculations for fixed sets of stellar parameters. For the simulations series constituting the basic testing process, we consider a system that is similar to the well-studied O-star ζ-Puppis. The basic stellar parameters for these calculations are listed in Table 2 and are adopted from Puls et al. (1996). These values imply a stellar radius of R⋆ = 1.317 × 10 12 cm. Additionally, we use Ne/W = 10 15 cm −3 in the ζ-Puppis calculations. Performing an ionisation and excitation calculation as outlined in the previous section, we find a mean electron scattering cross section of σe = 0.34 cm 2 g −1 throughout the wind, which corresponds to the Eddington factor Γe = 0.5. In all MCRH wind simulations presented in the following, we account for the effect of electron scattering by reducing the stellar mass by the factor (1 − Γe), as described in Sections 2.1 and 3.5.  (24). In the transition region, between the optically thick and thin regimes, a power-law fit according to Equation (A4) has been performed (dashed blue). The parameters of this fit are included in Table 3 and the predicted wind state (in terms of t), according to the CAK theory and the obtained k and α values, is marked by a green cross. For comparison, the force multiplier values of Abbott (1982), determined for a star with T eff = 4 × 10 4 K and Ne/W = 1.8 × 10 8 , 1.8 × 10 11 , 1.8 × 10 14 cm −3 are included (grey open symbols).

CAK fitting
In the CAK theory, the wind structure may be readily determined from the stellar properties and the power-law parameters k and α of the force multiplier. However, the wind characteristics, massloss rate and terminal wind speed, are sensitive to the exact values of these parameters. Given our aim at assessing the utility of our approach for solving the line-driving problem, we refrain from consulting previous studies, such as Abbott (1982) or Pauldrach (1987), which provide these parameters for a wide range of stellar conditions. Instead, we use the Monte Carlo routine of MCRH to determine the values for k and α. This way, we ensure that the MCRH-CAK/MCAK comparisons are not obscured by differences in atomic data sets, the ionisation and excitation treatments or the stellar parameters. In particular, we determine the values for k and α by carrying out the force-multiplier summation (c.f. Abbott 1982)  Figure 2 together with the fitting curve. As a reference, the results of Abbott (1982), obtained for comparable physical conditions are included. In general, both calculations agree, but noticeable differences are present, which are most likely a consequence of our simplified treatment of ionisation and excitation, differences in the atomic data sets and of the use of realistic stellar atmosphere models instead of a simple black body by Abbott (1982). Given the purpose of our calculations, these differences are irrelevant but highlight the need for consistency be- Table 3. Overview of the main CAK parameters determined in the fitting process and used throughout this work. In addition, the wind properties, as predicted by the CAK theory according to these parameters are listed in the bottom rows. The adopted values for Ne/W , ρ 0 and Γe are included for completeness.  Table 3. Based on the force multiplier parameters, the CAK predictions for the massloss rate and the terminal wind speed are calculated and also provided in the table. Including the finite-cone effect according to the MCAK approach increases the terminal wind speed by a factor of ∼ 2.68 to about u∞ ≈ 2360 km s −1 and reduces the mass-loss rate toṀ ≈ 2 × 10 −5 M⊙ yr −1 . With these modifications, the terminal wind speed is compatible with previous investigations of the wind of ζ-Puppis (e.g. Puls et al. 1996). However, the density in our wind is too high. Specifically, the mass-loss rate is about a factor of 3.5 larger than established by Puls et al. (1996). Again, our simplified treatment of the wind ionisation and excitation conditions is most likely the cause for these discrepancies. We stress once more, however, that as long as we use the simplified description of the wind conditions consistently in the CAK/MCAK calculations and the MCRH simulations, the conclusions drawn from comparing the corresponding results are unaffected.

General parameters
In the MCRH simulations, the evolution of the wind material is followed until a steady-state structure emerges. All calculations are carried out on a non-uniform spherical mesh with 512 cells, which span the region between 1 R⋆ and 10 R⋆. The cell size increases outwards from 1.76 × 10 −3 R⋆ at the stellar surface to 1.73 × 10 −1 R⋆ at the outer edge of the domain. During each time step, the incident radiation field is discretised by 5 × 10 4 packets, which sample the black-body spectrum in the wavelength range between λmin = 228 Å and λmax = 22800 Å. The lower edge corresponds to the ionisation edge of He II. Any radiation more energetic than this threshold is assumed to be rapidly removed by bound-free absorptions by helium and consequently not included in our consideration (see Sim 2004, for a similar strategy). For the initial wind state in the MCRH calculations, we adopt a structure that is similar to the CAK/MCAK solution but scaled up or down. By experimenting with other initial wind configurations, we have explicitly verified that the final stationary wind solution, which is found in the calculations, is insensitive to the details of the initial state. Furthermore, we have also ensured that enough Monte Carlo packets are used in the simulations. Increasing the number does not change the overall shape of the radiative acceleration but only reduces the strength of the stochastic fluctuations (see Figure 4).

Simulation series layout
Using the general parameters outlined in the previous sections, we perform a series of MCRH simulations of a ζ-Puppis-like wind. In this series we successively increase the level of physical detail (see Section 4.7). For the first stage, the point-source approximation and an unattenuated radiation field will be used, similar to Castor et al. (1975). In the second set of calculations, the finite extent of the star will be taken into account (similar to Friend & Abbott 1986, Pauldrach et al. 1986and Kudritzki et al. 1989. During these two first steps of the series, the ionisation and excitation calculation is further simplified by setting all recombination fractions ζi = 1 and by assuming a Boltzmann excitation formula [i.e. W = 1 in Equation (20)]. Finally, in the last stage of the series, the assumption of an unattenuated radiation field will be dropped. In this calculation, we then use the full ionisation and excitation description as outlined in Section 3.7, including recombination fractions and the geometric dilution factor in the excitation balance. This leads to a mild variation of the degree of ionisation in the wind. When applicable, the results of the MCRH calculations will be confronted with the predictions of the CAK and MCAK theory. These comparisons are performed on the basis of the stationary wind structure found in all calculations. With this in mind, we further reduce the computational complexity in all MCRH simulations by following all Monte Carlo packets until they escape the wind during each time step. This is equivalent to assuming a light propagation speed much larger than the fluid flow and constitutes a common strategy in Monte Carlo radiative transfer approaches, for example in PYTHON (Long & Knigge 2002) or TARDIS (Kerzendorf & Sim 2014). Adopting this technique absolves us of the need to invest computational power into the initial build-up of the radiation field. As long as we are solely interested in the final stationary wind state, this simplification is justified. For future calculations, which aim at investigating the dynamical behaviour of the wind structure, this modification of the Monte Carlo radiative transfer scheme will be dropped.

MCRH CAK/MCAK module
To better judge the outcome of the MCRH-CAK/MCAK comparison, we include an additional set of numerical calculations performed with our approach. Instead of relying on the Monte Carlo radiative transfer procedure to determine the line-driving force, we incorporate the analytic power-law expressions of the CAK/MCAK approach. In each simulation cycle, after the fluid dynamical and central gravity substeps, the instantaneous value of the dimensionless optical depth parameter t is determined in each cell and the force multiplier M (t) calculated accordingly. In this procedure, the central star may be treated either as a point source [see Equation (A4)], or its finite extent may be taken into account [see Equations (A8) and (A9)]. To determine the required instantaneous velocity gradient in each cell, we rely on the algorithm by Fornberg (1988) to construct finite differences on arbitrarily spaced grids. We will refer to all numerical radiation hydrodynamical calculations performed with this module as CAK-RH and MCAK-RH respectively.

Calculation in the point-source limit
In the first set of MCRH calculations, the central star is approximated by a point source and the radiation field is assumed to be unattenuated. Within the Monte Carlo approach, both simplifications are realised by injecting the packets only on radial trajectories, i.e. replacing the random number experiment (9) by µ = 1, and by not changing the propagation direction in interaction events, i.e. µ a 0 = µ b 0 . In Figure 3, the final stationary wind state, which establishes in the MCRH calculation, is shown in terms of the density, velocity and mass-loss rate. The comparison with the results of the CAK-RH simulation and the analytic predictions according to the CAK theory shows a very good agreement. This positive finding is a first indication for the accuracy and utility of our Monte Carlo-based scheme to address the line-driving problem in hot-star winds selfconsistently.

Including the finite-cone effect
During the second stage of the wind simulation series, the pointsource approximation is dropped and the finite extent of the central star is taken into account, analogously to the MCAK approach of Pauldrach et al. (1986) and Kudritzki et al. (1989). In the MCRH calculations during this stage, we allow Monte Carlo packets to propagate on non-radial rays as well by sampling the initial direction according to Equation (9).
As in the point-source calculations, the stationary wind structure obtained with MCRH is compared with the analytic predictions, now according to the MCAK theory as outlined in Appendix A. Again, we include the numerical results calculated with the MCAK-RH version of our scheme. When accounting for non-radial photon propagation paths, the Monte Carlo-based results agree again very well with the analytic predictions and the MCAK-RH calculation, as illustrated by Figure 3.
As expected from numerous previous studies, most notably from Friend & Abbott (1986), Pauldrach et al. (1986) and Kudritzki et al. (1989), the inclusion of the finite-cone effect leads to higher wind velocities, in our case by a factor of 2.5 in terms of the velocity at r = 10 R⋆ compared to the point-source case. At the same time, the mass-loss rate drops by a factor of 0.5. Figure  4 illustrates the difference in the radial dependence of the radiative acceleration, which is responsible for the change in the structure of the wind flow.

Full inclusion of scatterings
After having established the basic applicability of our approach to the line-driving problem in the first stages of the simulation series, we take another step towards a realistic description of the radiation field in hot-star winds by dropping the unattenuation approximation and accounting for the full scattering process in interactions. Now the entire procedure described in Section 3.4 is performed: in particular, the emergent propagation directions in atomic line interactions are drawn according to Equation (12). Accounting fully for the multiple scattering phenomenon, these calculations reach beyond the capabilities of the basic CAK and MCAK approaches (as outlined in Appendix A). We emphasise, that the key consequence of multiple scattering lies in the capability of line interactions to lengthen the photon propagation trajectory. By this process, photons may potentially exert a stronger acceleration onto the wind material, compared to the unattenuated case in which they may also interact multiple times but their trajectories are always straight lines. Once the unattenuation of the radiation is dropped in the MCRH simulations, some packets may potentially be back-scattered onto the stellar disc. To counteract the luminosity loss induced by this process, we follow Lucy & Abbott (1993) and scale the packet energies by a constant factor, which ensures that in each time step, net radiative energy amounting to the luminosity L⋆ is streaming into the wind. This process may be interpreted as a colour correction of the stellar spectral energy distribution (see Lucy 1999b, for a similar strategy in the context of calculating synthetic observables for supernovae). Accounting both for the finite-cone effect and the full scattering procedure, we again determine the stationary wind structure with MCRH. The resulting wind velocity, the density stratification, and the mass-loss rate are included in the summary plot of Figure  3. Compared to the unattenuated calculations with the finite-cone effect, we find a slightly slower wind and, as shown in Figure 4, a weaker radiative acceleration. At first glance, this seems to contradict the statement about the effect of multiple scattering in the introductory part of this section. But one has to bear in mind that in the CAK/MCAK description of the line-driving problem, in each interaction the full photon/packet momentum is transferred onto the wind material. This is comparable to a purely absorbing medium, with the important difference that the photon trajectory is not terminated in the CAK/MCAK description. Instead, the same photon may still contribute to the acceleration in regions of the wind located at larger radii. In the full scattering case, in which the re-emission of the line-interaction is taken into account, no momentum is transferred onto the wind on a CAK/MCAK like trajectory since it involves forward-scatterings only. Thus, the momentum transfer in the CAK/MCAK case may generally be overestimated and only the photons that are on trajectories that have been significantly lengthened due to many successive scatterings may contribute comparably to the CAK/MCAK description. This phenomenon is illustrated in Figure 5. A general reduction of the line-driving force once the unattenuation assumption is dropped has already been found in previous studies, e.g. Abbott & Lucy (1985). Red lines corresponds to calculations, in which the finite-cone effect is included and green to those which also include the full scattering procedure. All MCRH results are presented by solid lines. Where applicable, the analytic predictions according to the CAK and the MCAK theory are included as dotted lines. In these cases, also the results obtained with CAK-RH/MCAK-RH are given as an additional reference (dashed lines).

Additional diagnostics
By recording the details of all interaction events performed by the Monte Carlo packets in the final calculation of the simulation series, the origin of the radiative acceleration can be studied. The contributions of the different elements and ionisation stages, averaged over the entire computational domain, are illustrated in Figure  6. This highlights that the line-driving force mainly derives from interactions with lines of iron, nickel some intermediate mass elements and the CNO group. The relative importance of the different contributions, however, changes throughout the wind, as shown in the right panel of the same figure. In our simulation, lines of iron group elements mostly contribute in the inner part of the wind, close to the photosphere. Further out, the intermediate mass elements grow in importance. This finding is compatible with the investigation of Vink et al. (1999), in which the importance of iron for the radiative acceleration in the lower parts of the wind has been highlighted in the context of the bi-stability jump. We stress, however, that due to the simplified treatment of ionisation and excitation in our simulations, the results presented in Figure 6 should not be over-interpreted, but viewed as an illustration of the diagnostic possibilities of our approach.
By recording the interaction histories of all packets we can also investigate the importance of multiple scattering in our simulation. Figure 7 shows the number of interactions performed by the packets that escaped through the outer edge of the computational domain. In our particular simulation, most packets never explicitly interacted on their way out. They may have still contributed to the line-driving force as long as they came into resonance with at least one line transition (see discussion in Section 3.6). The majority of those that interacted, however, performed multiple scattering events.

COMPARISON WITH DETAILED CALCULATIONS
Having established the basic performance of our method under simplified conditions, we now compare our approach with modern, sophisticated techniques for determining the structure of hot-star winds, in particular with the approach developed by Müller & Vink (2008, MV08 hereafter). This method constitutes an advancement of the original approach of Vink et al. (1999) and Vink et al. (2000). A Monte Carlo radiative transfer calculation is used to determine the local line-driving force for a given static wind structure.
Here a velocity structure that is more general and flexible than the β-type law is used. According to the reconstructed line acceleration, the mass-loss rate and the parameters of the wind velocity law are iteratively updated until a converged wind structure has been found.

Parameters
MV08 test their approach by predicting the wind structure of an O5-V main sequence star. For direct comparison, we adopt parameters in our simulation to match their calculation. In Table 4, all quantities which have been changed with respect to the calculations presented in Section 4 are listed. These choices imply an Eddington factor of Γe = 0.210, which is very close to the value quoted by MV08.
The comparison with the Müller & Vink (2008) technique will be based on MCRH calculations that incorporate all techniques described in Section 3. Since no CAK or MCAK-like simulations are performed in this context, the CAK fitting procedure of Section 4.3 does not have to be repeated for the current stellar parameters.

Results
With the stellar parameters listed in Table 4 a full scattering MCRH simulation, similar to the calculation presented in Section 4.7.3, is performed to determine the structure of the wind of the O5-V main sequence star. The result is shown in Figure 8 in terms of the stationary velocity and mass-loss rate and compared to the structure found by MV08. To obtain the comparison MV08 wind velocity, we use the approximate expression, Equation 39 in MV08, which is only strictly valid in the supersonic wind regime.
As seen in Figure 8, both approaches predict winds that are quite similar in shape. However, the velocity structure found by MCRH accelerates quite quickly towards the terminal wind speed, whereas the wind velocity law predicted by MV08 approaches its final value in a gentler fashion. The two winds deviate slightly in absolute values: the MCRH wind reaches a velocity of u = 3065 km s −1 at r = 10 R⋆, while the MV08 solution lies at u = 2719 km s −1 . The mass-loss rates of the winds agree on a similar, namely a ten percent level, withṀ = 7.97×10 −7 M⊙ yr −1 being obtained in the MCRH calculation andṀ = 8.99×10 −7 M⊙ yr −1 found by MV08.
These minor discrepancies may partly be related to the different solution strategies followed in the two approaches. In our method, a radiation hydrodynamical calculation is relaxed to a steady state solution without any major restrictions on the shape of the line acceleration and velocity structure. The MV08 approach relies on a parameterisation of the line driving force and in turn of the velocity law, thus allowing only wind structures of a certain family. Differences between the two calculations will also arise because of the very simplistic treatment of ionisation and excitation in our method. Overall, however, the level of agreement between the two calculations is good and lends credence to their complementary use in the study of stellar winds.

DISCUSSION
The results of the wind calculations and the diagnostic capabilities presented in the previous sections, in particular the good agreement between the MCRH results and the CAK/MCAK predictions and the comparison to the MV08 technique, clearly demonstrate the utility and accuracy of our Monte Carlo-based radiation hydrodynamical approach for solving the line-driving problem in hot-star winds. In the following, we briefly discuss some particular features to our approach and also comment on some of its current limitations. We also sketch possible alterations of our approach aimed at alleviating these shortcomings and discuss potential future applications.

Numerical performance
In general, Monte Carlo radiative transfer techniques are rather computationally demanding, since many Monte Carlo quanta have to be processed to reach the desired statistical fidelity. This drawback is typically balanced by the very favourable parallelisation properties of such calculations (see discussion in section 6.4). In our case, the Monte Carlo-based radiation hydrodynamical simulations are indeed costly. In the full scattering simulation, the final stationary state established after about 4800 simulation cycles. On of these takes roughly 14 s on a Intel Xeon E5520 processor if 5 × 10 4 packets describe the incident radiation field. This increases to about 210 s once 10 6 packets are used. In these calculations about 10 GFlops (for 5 × 10 4 packets) and 160 GFlops (10 6 packets) were executed per cycle. By comparison, only 0.015 s are required and 1.3 MFlops executed per cycle if the CAK-RH version is used instead. We emphasise, however, that all calculations presented in this work were carried out serially on a single processing core. We stress, that the performance of the Monte Carlo scheme may be significantly improved by devising a parallelisation scheme and distributing the workload onto many cores. For all calculations, the executable was produced with the gcc compiler, version 4.7.3, using the optimisation level -O3.

Influence of Monte Carlo noise
The probabilistic nature of Monte Carlo techniques leads to the inevitable introduction of stochastic fluctuations (see, for example, Carter & Cashwell (1975) for a discussion of Monte Carlo errors). However, Noebauer et al. (2012) and Roth & Kasen (2015) demonstrated that, with appropriate reconstruction techniques, the fluctuations can be controlled and Monte Carlo-based radiation hydrodynamical simulations can indeed be performed. The simulations performed in this work and presented in the previous section illustrate that the Monte Carlo noise, which is clearly visible in the reconstructed radiative acceleration (c.f. Figure 4), does not prohibit the use of Monte Carlo methods to study the line-driving problem. In fact, the obtained velocity and density structure (shown in Figure 3) is remarkably smooth. Due to their pure stochastic nature the fluctuations cancel in an average sense (but see also remarks in Section 6.4).

Full non-LTE treatment
The calculations presented here assumed either LTE or relied on the approximate non-LTE prescription presented in Section 3.7. Non-LTE effects, however, play an important role in stellar winds (see, e.g., Pauldrach 1987;Pauldrach et al. 1994;Puls et al. 2005, for non-LTE calculations of hot-star winds) and should be included in future calculations.
In the general non-LTE situation, the level populations of excited states are influenced by the radiation field. Lucy (2002) and Lucy (2003) describes a simple and elegant procedure to reconstruct these radiative rates from a Monte Carlo simulation. In a first step of refining the physical detail in the interaction treatment of our approach, this procedure could be adopted to determine a non-LTE excitation and ionisation balance during the Monte Carlo radiative transfer step. This strategy has already been followed in previous Monte Carlo-based studies, for example by Sim et al. (2005), using the PYTHON code (Long & Knigge 2002) and in non-LTE radiative transfer calculations in Keplerian discs by Carciofi & Bjorkman (2006).
To go beyond the pure resonant scattering approximation in the line-interaction treatment, either the simple branching scheme of Lucy (1999b) (see Mazzali 2000, for Monte Carlo-based radiative transfer calculations in Type Ia ejecta using this scheme) may be incorporated or the so-called Macro-Atom formalism introduced by Lucy (2002Lucy ( , 2003 may be used. The latter fully accounts for all processes that may excite or de-excite line-transitions, including collisional processes. We emphasise, that the utility and accuracy of these techniques has been carefully examined and established (Lucy 2002(Lucy , 2003(Lucy , 2005. Moreover, they have already been incorporated into a number of Monte Carlo radiative transfer frameworks, including ARTIS (Kromer & Sim 2009), PYTHON (Sim et al. 2005) and TARDIS (Kerzendorf & Sim 2014).

Multidimensionality and non-Sobolev treatment
Observations and also theoretical investigations indicate that linedriven winds are neither stationary nor smooth, but feature a variable and inhomogeneous structure (see, e.g., overview by Puls et al. 2008). These findings advocate a dynamical and fully multidimensional treatment of the problem. Monte Carlo techniques are suited for this task since the formalism readily generalises (conceptually) to arbitrary geometrical configurations and also exhibits excellent scaling properties on parallel processors -a desirable feature when facing the increased numerical workload of multidimensional calculations (see, for example, comments by Kasen et al. 2006 andBaes et al. 2011, as well as the scaling tests by Roth &Kasen 2015 andHarries 2015). With this in mind an obvious extension of our method comprises a generalisation to multidimensional geometries to address inhomogeneous line-driven outflows. In this context, some work should be invested into devising an efficient interfacing of the standard parallelisation strategies of Monte Carlo and fluid dynamical techniques (see e.g. discussion in Baes et al. 2011 anddevelopments by Harries 2015.).
When investigating multidimensional line-driven mass outflows, a generalisation of the line interaction treatment that goes beyond the Sobolev approximation used here should also be considered: this would be necessary to study the line-driving instability (Lucy & Solomon 1970;Owocki 1994) which should occur in linedriven outflows and may play a part in understanding the clumping mechanism (see discussion in Puls et al. 2008 andVink 2015). Non-Sobolev Monte Carlo radiative transfer schemes have already been developed and used, for example by Knigge et al. (1995) and Kusterer et al. (2014) in the context of winds of cataclysmic variables. However, the computational costs of such schemes are significantly higher than their Sobolev counterparts. Moreover, the inherent Monte Carlo noise could potentially be more problematic than described in 6.2, since small perturbations in the line-driven wind may, in principle, self amplify in non-Sobolev calculations. Further investigation is required to asses the relevance of this potential caveat.

CONCLUSIONS
In this work we have introduced a new approach to solve linedriven stellar winds self-consistently by using a Monte Carlo-based radiation hydrodynamical approach. The key feature of this technique lies in the reliance on a Monte Carlo radiative scheme. This technique offers a number of advantages when dealing with complex interaction physics, in particular when multi-line effects play a role. Moreover, this Monte Carlo-based technique readily generalises to multidimensional geometries, a very advantageous feature for potential studies of inhomogeneous outflows.
Establishing the utility and accuracy of the introduced approach in solving the line-driving problem was the main focus of this work. Consequently, we designed the calculations to capture the essence of the line-driving problem, but adopted a number of simplifications. These do not interfere with the general applicability of our method to solve the local line-driving problem, but reduce the computational complexity. Most importantly, we adopted the Sobolev approximation and a simple and approximate non-LTE treatment to determine the excitation and ionisation balance.
Using our scheme, we successfully solved for the stationary structure of one-dimensional, spherically symmetric hot-star winds achieving good agreement with the predictions of the CAK and MCAK theory. We demonstrated that our method can also go beyond the capabilities of CAK and MCAK by dealing with an attenuated radiation field and the effects of multiple scattering. For the particular physical conditions investigated in Section 4, these effects lead to a reduction of the line-driving force and thus a slower wind velocity structure compared to the MCAK predictions. Finally, we compared results obtained with our approach to those computed by MV08 for an O5-V main sequence star and found good agreement, given the difference in approach and simplifications made here.
The successful outcome of our wind simulations demonstrates the possibility to use a Monte Carlo-based radiation hydrodynamical approach to model line-driven mass outflows. Thus, this approach holds promise for detailed multidimensional and selfconsistent studies of inhomogeneous winds, including multi-line effects. With a fully multidimensional version of our approach, line-driven outflows in systems other than hot-star winds may be addressed as well. In particular, the winds emanating from accretion discs in cataclysmic variables (e.g. Proga et al. 1998;Noebauer et al. 2010) or active galactic nuclei (e.g. Proga et al. 2000;Proga & Kallman 2004;Higginbottom et al. 2013) may be investigated self-consistently and in great detail with such a Monte Carlo-based approach.
if the central star is assumed to radiate as a point-source (again following Castor et al. 1975). It should be noted that, throughout this work, we neglect any influence of changes in the ionisation, which is equivalent to setting δ = 0 in the force multiplier formulation proposed by Abbott (1982).
In the point-source limit, the solution to the momentum equation (1), with the CAK line-driving force, results in the wind velocity law with the terminal speed being a multiple of the local escape speed from the photosphere, uesc: The constant mass-loss rate of the wind is given bẏ Once the finite extent of the star is taken into account, the force multiplier is modified by a finite-cone correction factor We follow Kudritzki et al. (1989) and predict the wind structure according to their approximate analytic solution technique. In particular, we adopt their proposed approach for the case of a "frozen-in" ionisation state (δ = 0, c.f. Kudritzki et al. 1989, section 4.1). The wind velocity is assumed to be very close to a β-type law and the mass-loss rate decreases with respect to the original CAK case according toṀ The wind velocity in this approach follows from performing the integration with Z(v, α, β) = fN (v, α, β) In these expressions, the inverse of the radial distance relative to the photosphere, v = R⋆/r, is used. Throughout this work, the MCAK equations are solved for β = 0.8 (see Pauldrach et al. 1986;Kudritzki et al. 1989, for a motivation of this value).