Issue 
A&A
Volume 670, February 2023



Article Number  A41  
Number of page(s)  22  
Section  Stellar atmospheres  
DOI  https://doi.org/10.1051/00046361/202243997  
Published online  03 February 2023 
Progress towards a 3D Monte Carlo radiative transfer code for outflow wind modelling
^{1}
Masaryk university, Faculty of Science,
Kotlářská 2,
Brno, Czech Republic
email: fisak@physics.muni.cz
^{2}
Astronomical Institute of the Czech Academy of Sciences,
Fričova 298,
251 65
Ondřejov, Czech Republic
^{3}
Heidelberger Institut für Theoretische Studien,
SchlossWolfsbrunnenweg 35,
69118
Heidelberg, Germany
Received:
10
May
2022
Accepted:
24
October
2022
Context. Radiative transfer modelling of expanding stellar envelopes is an important task in their analysis. To account for inhomogeneities and deviations from spherical symmetry, it is necessary to develop a 3 D approach to radiative transfer modelling.
Aims. We present a 3 D Monte Carlo code for radiative transfer modelling, which is aimed to calculate the plasma ionisation and excitation state with the statistical equilibrium equations, moreover, to implement photonmatter coupling. As a first step, we present our Monte Carlo radiation transfer routines developed and tested from scratch.
Methods. The background model atmosphere (the temperature, density, and velocity structure) can use an arbitrary grid referred to as the modGrid. The radiative transfer was solved using the Monte Carlo method in a Cartesian grid, referred to as the propGrid. This Cartesian grid was created based on the structure of the modGrid; correspondence between these two grids was set at the beginning of the calculations and then kept fixed. The propGrid can be either regular or adaptive; two modes of adaptive grids were tested. The accuracy and calculation speed for different propGrids was analysed. Photon interaction with matter was handled using the Lucy’s macroatom approach. Test calculations using our code were compared with the results obtained by a different Monte Carlo radiative transfer code.
Results. Our method and the related code for the 3 D radiative transfer using the Monte Carlo and macroatom methods offer an accurate and reliable solution for the radiative transfer problem, and are especially promising for the inclusion and treatment of 3 D inhomogeneities.
Key words: stars: atmospheres / stars: winds, outflows / radiative transfer / methods: numerical
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
The stellar wind is a kind of an outflow from the star that is stable on a longterm scale. The stellar radiation field is often strong enough to accelerate matter and create a stellar wind. Intensive radiationdriven stellar winds are observed, for example, around Otype stars and WolfRayet stars. Stellar wind affects many astrophysical processes, including stellar evolution, the chemical composition of galaxies, and the dynamics of interstellar matter. A star can lose a significant amount of matter per unit of time; this amount is characterized with a physical quantity named the ‘massloss rate’, which is important for the stellar evolution and for the galaxy chemical composition. The accelerating mechanism of atoms by line radiation was suggested by Milne (1926); linedriven winds were later described in detail by Lucy & Solomon (1970). The first hydrodynamical solution of the linedriven wind was provided by Castor et al. (1975, also known as the CAK model). This model assumes stationary, homogeneous, and spherically symmetric wind with a uniform outflow.
There is a wide range of existing codes for the stellar wind modelling. Stellar winds of hot massive stars are routinely modelled assuming spherical symmetry (i.e. 1 D), and taking into account the possibility of departures from the local thermodynamic equilibrium (LTE). The latter approach is usually referred to as NLTE. The codes solve the radiative transfer equation (hereafter RTE) and the kinetic equilibrium equations (the NLTE line formation problem), typically supplemented by the equation for temperature (either radiative equilibrium or thermal balance). There exist three such widely used codes, namely CMFGEN (e.g. Hillier 1987, 1990, 2012; Hillier & Miller 1998, and references therein), FASTWIND (e.g. SantolayaRey et al. 1997; Puls et al. 2005, 2020), and PoWR (e.g. Hamann & Gräfener 2003, 2004; Sander et al. 2017). In addition to these codes used for analysis of observed spectra, other codes construct spherically symmetric 1 D wind models using the solution of hydrodynamic equations (Pauldrach et al. 2001; KrtiCka & Kubát 2017; Sander 2017; Sundqvist et al. 2019).
The assumptions of spherical symmetry and time independence were later shown not to be correct. The spherical symmetry can be broken by the rotation of the star (Puls et al. 1993; Owocki & Cramner 1994; Petrenz & Puls 2000), by the magnetic field (udDoula 2015), or by the accretion onto the compact object (Blondin et al. 1990; Feldmeier et al. 1999).
These processes are not the only ones which cause the wind asymmetry. Stability analysis of stellar winds done by Carlberg (1980); Owocki & Rybicki (1984); Owocki & Puls (2002); Feldmeier et al. (1997), and many others showed that density and velocity are not monotonically changing functions. The corresponding instability was found also by radiative hydrodynamical simulations in 1D (Feldmeier et al. 1995; Owocki et al. 1988; Runacres & Owocki 2002) or in 2D (e.g. Dessart & Owocki 2003, 2005; Sundqvist et al. 2018; Driessen et al. 2021). The inhomogeneities are called clumps. They can be introduced even into simplified 1D models using fiducial factors such as the clumping factor, or a volume filling factor (e.g. Abbott et al. 1981), or an effective opacity formalism (Sundqvist & Puls 2018). However, there are possibilities to solve the radiative transfer in clumped media precisely.
Consequently, there is a growing need for full 3D NLTE modelling of massive star winds. Although there is a number of 3 D radiationhydrodynamics codes for LTE atmospheric modelling (e.g. Nordlund & Stein 2009; Ludwig & Steffen 2016; Freytag et al. 2019), which were successful especially in the modelling of convection in cool stars, related NLTE calculations are usually limited to the solution of the NLTE line formation problem (see Asplund & Lind 2010, for a review). Nevertheless, first steps towards full NLTE wind models have already been taken; however, due to the complexity of the problem caused mainly by the necessity to include a general velocity field, only simplified problems were solved. The key problem is the way how the radiative transfer problem is treated. Several different numerical methods have been applied to solve this problem. Adam (1990), Lobel & Blomme (2008), and Hennicker et al. (2018) used the finite volume method, while Papkalla (1995), Korčáková & Kubát (2005), Georgiev et al. (2006), Zsargó et al. (2006), Leenaarts & Carlsson (2009), Ibgui et al. (2013), Štěpán & Trujillo Bueno (2013), and Hennicker et al. (2020) used the short characteristics method. All these codes require enormous computing time to solve even the simplest line formation problems.
An alternative to the solution of the multidimensional radiative transfer equation by means of these methods which determine the intensity of radiation is to use Monte Carlo (hereafter MC) radiative transfer (see Noebauer & Sim 2019, for a review). This method is very useful for the optically thin media such as circumstellar environments, and planetary nebulae. The main advantage of the MC method is that the solution is calculated naturally in a 3 D space. This allows us to fully include some 3 D phenomena such as wind inhomogeneities, that is clumping.
There were some intermediate steps towards the 3 D models. Oskinova et al. (2004, 2006) and Sundqvist et al. (2010) calculated pseudo2 D models, later on followed by a pseudo3 D model by Sundqvist et al. (2011). There are several codes for the MC 3D RTE calculation: HDUST (Carciofi et al. 2004, 2017; Carciofi & Bjorkman 2006), a 3 D NLTE radiative transfer code for Be stars, hot star winds, etc. The PYTHON code (Long & Knigge 2002, extended by Higginbottom et al. 2013) is capable of including analytical wind models, as well as disc, and supernova models. It uses cylindrical or a spherical grid and it also can read in arbitrary geometries. Šurlan et al. (2012, 2013) calculated wind clumping with the MC code. They solved radiative transfer in doublet resonance lines.
Also the backward ray tracing method is being used – this method can calculate spectra based on the observer’s position. This is useful for the calculation of emergent radiation from nonsymmetric objects. The standard Monte Carlo method deals with photon propagation through a medium. Usually we observe a star via detectors occupying a very small area. A very large amount of photons must be propagated to create a spectrum without a significant noise. The ray tracing method evades this disadvantage and solves RTE along specific emergent rays and then we get spectra for the given direction.
There are also several MC codes calculating the radiative transfer through the supernova ejecta. These codes calculate timedependent spectra to take the rapid variations of circumstellar ejecta into account. As an example, the code Tardis (Kerzendorf & Sim 2014; Vogl et al. 2019; Kerzendorf et al. 2022) enables the fast calculation of supernova spectra. More sophisticated codes (Lucy 2005; Kromer 2009; Kromer & Sim 2009) also simulated radioactive reactions in the ejecta, which provide another source of radiation. There are codes calculating the radiation of the star with both the supernova ejecta and the circumstellar disc.
We present the 3 D radiative transfer MC code which we developed from the scratch; the code is suitable for outflow wind modelling. In Sect. 2 we describe the NLTE MC radiative transfer method which is implemented in the code. The general code structure is present in Sect. 3, while in Sect. 4, we present the results of testing the propagation grid. Physical processes in the outflow, which are implemented in the code, are described in Sect. 5. Computational details are presented in Sect. 6, while comparison of our code with the Tardis code is presented in Sect. 7. In Sect. 8, we summarize and outline future work.
Fig. 1 Scheme of the macroatom interactions. 
2 NLTE Monte Carlo radiative transfer
Our radiative transfer code is based on Monte Carlo method of Lucy (2002, 2003). This method traces propagation of indivisible energy quanta (energy packets) through matter. The energy packets can have form of radiation energy packets (rpackets), internal energy packets (ipackets), and kinetic energy packets (kpackets). The packets may change their form in radiationmatter interactions.
The matter is quantized into macroatoms. Detailed properties of energy packets and macroatoms are described in Lucy (2002).
2.1 Radiative packets
Radiation in this method is quantized to radiation energy packets (rpackets), which behave such as basic photonic quanta, but their energy is quite different and does not correspond to a single photon, it is an ensemble of photons. The propagation and changes of properties of these packets represent the radiation energy transport through the atmosphere and the interaction of radiation with matter. The radiation packets can be changed during radiationmatter interaction (see Fig. 1) to packets of internal energy (ipackets) or kinetic (thermal) energy (kpackets).
Although the method is general and can be applied to many astrophysical radiation sources, here we restrict ourselves to stars and circumstellar regions. Since the Sobolev approximation is often used in supernovae modelling, consequently, we adopted this approximation to enable easy comparison with results obtained by supernova codes. Relaxation of this approximation in our code will be discussed in a forthcoming paper.
2.1.1 Creation of radiative packets
Creation of rpackets follows the lower boundary condition at the stellar surface. Packet position, frequency, and direction have to be randomly chosen. Starting positions of packets are chosen to be equally distributed across the stellar surface. If we denote the stellar radius as R_{*}, the initial packet position in the spherical coordinates is (1a) (1b) (1c)
Here τ and σ are randomly generated numbers in the interval [0, 1] (see also Noebauer & Sim 2019, Eqs. (26) and (27)).
The initial propagation direction of the rpacket is randomly chosen in the outward direction from the stellar surface: another random unit vector (2a) (2b) (2c)
Here τ′ and σ′ are randomly generated numbers in the interval [0, 1] (see also Noebauer & Sim 2019, Eq. (47)). The propagation direction in the stellar coordinate system is calculated as a double rotation (two rotationmatrix multiplication), then the initial packet direction d_{pack} is in the form (3)
where r_{i} represent the Cartesian coordinates of the position vector (Eq. (1a)). The packet initial frequency follows the radiation distribution emerging from the stellar surface, which may be the Planck distribution or, more precisely, the calculated synthetic spectrum of the stellar photosphere. Packets are created at the beginning of every iteration. As the packet is not a photon, its energy must be established to ensure the flux energy conservation. The packet energy E (emitted per unit time) is equal to (4)
and is the same for all energy packets regardless their frequency. Here N is the total number of generated packets and L_{star} is the stellar luminosity.
2.1.2 Propagation and interaction of radiative packets
The radiation energy packets (rpackets) are propagated through the wind. We denote physical quantities expressed in the rest frame (the frame connected with the centre of the star, hereafter RF) using the index #, and physical quantities expressed in the comoving frame (hereafter CMF) using the index ↝. The interaction of matter and radiation can occur only at a point where the comoving frame photon frequency v^{↝}equals the line transition frequency (v^{↝} = v_{line}).
Consider an rpacket with a rest frame frequency v^{#} at a position r with density p(r) and temperature T(r). We express the CMF frequency using the approximate Doppler shift formula (5)
where υ is the velocity vector of matter (in the RF) at a given point r, and n is the direction unit vector of a light ray. As the packet propagates through the medium with a velocity gradient, its CMF frequency is changing.
The continuum optical depth is calculated using the equation (6)
where r is the packet position, r_{R} is the position of the closest possible line interaction (resonance point) of the rpacket with matter and χ_{cont} (r′) is the continuum opacity, (7)
The first term corresponds to Thomson scattering, the second one to freefree processes, and the third one to boundfree processes. N_{I} is a total number of all atomic ions and N_{L} is a total number of energy levels of all ions.
The line optical depth for a transition between a lower level l and an upper level u integrated along a path s from 0 to s_{0} is given by (see Kromer 2009, Eq. (4.15)) (8)
where ϕ is the line profile, n_{l}, n_{u} are number densities of atoms at levels l (lower) and u (upper), respectively, ɡ_{l} ɡ_{u} are corresponding statistical weights, and B_{1u} is the Einstein coefficient for a transition from a level l to a level u. In the Sobolev approximation this equation simplifies to (9)
where v_{lu} is the transition frequency. We note that the derivative (ds/dv^{↝})v_{lu} is negative. If we assume spherically symmetric velocity fields, this expression simplifies (Castor 1974, see also Noebauer & Sim 2015, Eq. (22)) to (10)
where µ is angle cosine between n and υ(r), and (11)
according to Noebauer & Sim (2015, Eq. (7)). Here f_{lu} is the oscillator strength for a transition from a level l to a level u.
In Monte Carlo radiative transfer, each rpacket is allowed to travel a random optical depth τ_{rand} and then it interacts. This optical depth is (see, e.g. Whitney 2011) (12)
where z is a random number. We can decide which type of the process (line or continuum) happens by comparison of three optical depths, random optical depth τ_{rand} (Eq. (12)), line optical depth τ_{line} (Eq. (10)), and continuum optical depth τ_{cont} (Eq. (6)). The line and continuum optical depths are calculated by integration along the photon path. The continuum optical depth increases continuously, while the line optical depth in the Sobolev approximation changes by jumps at particular resonance points. If (13)
the photon continues its path, otherwise it interacts. The jump in line optical depth caused by the Sobolev approximation enables a simplification of determination whether a line or continuum transition happens. The type of interaction (line or continuum) is chosen according to a simple rule (as in Kromer 2009), it corresponds to the first transition which caused invalidity of Eq. (13).
If a line transition happens, then the specific line is chosen according to the scheme mentioned above, namely the one which caused invalidity of Eq. (13) is chosen.
If a continuum process is chosen, it has to be decided (using random numbers) which specific process occurs. Currently, boundfree and freefree transitions, and scattering on free electrons (Thomson scattering) are taken into account. The Thomson scattering causes reemission of a new packet in a random direction (here we assume isotropic scattering on free electrons) with the same CMF frequency as of the original packet. The CMF quantities are recalculated to corresponding RF values. A freefree process can transform an rpacket into thermal kinetic energy (a kpacket). Boundfree processes contribute to both kinetic and internal energy (k and ipackets, respectively). Following Lucy (2003, Eq. (27)), the packet converts to an ipacket with a probability w_{i} = v_{i}/v^{↝}, where v_{i} is the ionization edge frequency. The packet is transformed to the kpacket otherwise (w_{k} = 1 − v_{i}/v^{↝}).
Freefree processes transform an rpacket into a kpacket and vice versa. Boundfree processes contribute to both internal energy and kinetic energy. For each rpacket we decide (again using random numbers) whether it changes to a kpacket or to an ipacket. The line event changes the rpacket to an ipacket (macroatom activation).
2.2 Macroatom statistical equilibrium equations
Let us write the kinetic (statistical) equilibrium equation for the level i (see Hubeny & Mihalas 2015, Chapter 14) (14)
where R_{ij} represents the radiative rates from the atomic state i to the atomic state j, C_{ij} the collisional rates, and NL denotes the total number of levels for the given atom (including all atomic ions). Let us assume that the levels are ordered with increasing excitation energy, ∀_{i} : ɛ_{i} < ɛ_{i+1} starting with the ground level of the neutral atom and ending with the fully ionized atom (with only one energy level). Energy of a level is equal to the sum of ionization energy of an ion and excitation energy of the state. According to Lucy (2002), the sums in Eq. (14) can be split into two parts: transitions to levels with excitation energies either lower or higher than the excitation energy of the level i, (15)
We define the total rate for the transition i → l as Ƥ_{il} = n_{i}(R_{il} + C_{il}), after a simple manipulation we get (see Lucy 2002, Eq. (4)): (16)
It is possible to rewrite this equation in the form of energy absorbed and emitted by processes connected with the level i (see Lucy 2002, Eq. (5)). To this end, we define the unit volume energy rates: (radiative energy emission rates) and (collisional energy emission rates), and the absorption energy rates (radiative) and (collisional). These terms express energy absorbed or emitted per unit time and unit volume. The radiative energy rates can be expressed as (17a)
the collisional energy rates as (17b)
Let us express the sum of these rates using total transition rates Ƥ, (18a) (18b)
Multiplying Eq. (16) with ɛ_{i} and substituting for , from Eq. (18) we get (19)
which is equivalent to Eq. (16). Terms with Ƥ_{ij} in Eq. (19) describe energy flows between individual states. The Eq. (19) deals with macroscopic energy flow rates in a volume element. Flows can be quantized into indivisible energy packets and atoms in the volume element can be taken as a macroatom with discrete energy.
2.3 Internal energy packets
Absorption of the radiation energy packet (the rpacket) via excitation or ionization activates the macroatom. Internal energy packets (ipackets) represent excitation energy in excited atoms. Every ipacket must contain information about atom(s) it represents: an atomic number, an ion number, and a level index. The ion number and the level index may change during the ipacket processing. To describe particular processes changing excitation and ionization state, we use following index notation: i – the current level, l – a lower level, u – an upper level, m – a level in a lower ionization state than the level i, and p an upper level ionization state. If necessary, we specify the level with a lower index (E – element, I – ion, L – level) and its lower index classifies its concrete value, for example represents the basic level of the ion k (the state is denoted as i).
Let us denote the total energy loss (left hand side of Eq. (19)) caused by transitions from the energy level i as (see Kromer 2009, Eq. (4.38)) (20)
where and denote radiative and collisional deactivation (introduced in Eq.(17)), respectively, and and mean internal upward and downward jumps, respectively.
The macroatom processing is probabilistic. Every possible process occurs with a probability p_{if}, where i denotes the initial (level, ion, element) and f the final atomic state (again: level, ion, element). Here f stands for l, u, m, or p, as introduced above. The internal processes occur until the macroatom is deactivated by one of the noninternal (deactivation) processes, namely^{1} radiative deexcitation from the state i to the state l – the probability of this process is equal to: (21a)
where is a rate in the Sobolev approximation, denoted with an upper index S, radiative recombination (considered only from an ion ground level): (21b)
R_{im} is the rate for recombination (Eq. (B.11)), the index i_{Lo} means the lowest level of the ion, collisional deexcitation – in this case: (21c)
C_{il} is the collisional deexcitation rate (see Appendix B.2), and collisional recombination (considered only from an ion ground level): (21d)
where C_{im} is the collisional recombination rate (see Appendix B.2).
The internal processes (which do not deactivate the macroatom) taken into account include internal downward jump within the current ionization state: (22a)
internal downward jump to the lower ionization state: (22b)
internal upward jump within the current ionization state: (22c)
where is the rate of transition to the upper level, and, finally, internal upward jump to a higher ionization state: (22d)
Equations (21) and (22) contain rates which are not corrected for stimulated emission. If we want to use the rates corrected for the stimulated emission, we replace rates with Eqs. (B.8) and (B.9). By analogy, recombination can be treated as negative photoionization via the Eq. (B.15). There are no downward internal jumps to the ground level of the lowest ionization stage (the neutral atom), since the energy ɛ_{l} of this level is zero.
The choice of a particular process. The basic algorithm is shown in the Fig. 2). Let us denote the total rate of all transitions from the ith level as (23)
where type has a value from the set rad deexc., rad. recomb., col deexc, col recomb, internal: upward and downward jump within or to lower/upper ionization state. Let us number the variable type as integers running from 1 to the maximal number of processes included. Based on this definition of L_{i} we define the cumulative rate (24)
used in the upper panel of Fig. 2. Clearly ℒ_{0} = 0. We can also define a cumulative probability for a transition (25)
Here x is an index of the last transition taken into account. We generate a random number z ϵ [0, 1] (with the uniform distribution). The random number z is then multiplied with the total rate L_{tot} = ℒ_{N} where N is the total number of types of processes (N ≤ 8, cf. Eqs. (21) and (22)). Then we check which interval (ℒ_{q−1}, ℒ_{q}] includes the number z·L_{tot}: the condition ℒ_{q−1} ≤ z L_{tot} < ℒ_{q} (where q ϵ {1, 2,…, N} and ℒ_{0} = 0) is tested whether it is satisfied. For this selected index q we search in a similar way for an index x satisfying the relation ℒ_{q, x−1} < z L_{tot} < ℒ_{q, x} and these indexes q and x are denoted as Q and X, respectively. This twostep process is faster than a onestep process which goes through all transitions once and determines directly a specific transition.
When the radiative deexcitation in line (Eq. (21a)) or the radiative recombination occurs, the ipacket transits into the rpacket. We have to calculate its new frequency and direction. The direction unit vector is sampled before the calculation of rates using the equation (26)
where Θ and Φ are determined using Eqs. (1c) and (1b), respectively.
The rates depend on the Sobolev optical depth via Eqs. (B.4) and (B.7), see also Klein & Castor (1978). For spherical nonhomological velocity fields µ does not vanish and Eq. (10) depends on µ.
If the radiative deexcitation (RD) occurs a new packet frequency must be calculated in the following way: We find all possible downward transitions from the state i. The probability of choosing lth line is given by (27)
The next step is a choice of a process similar to the lower part of the Fig. 2. The chosen line sets up the CMF frequency of a new created rpacket. Currently a δfunction is assumed as a line profile. The rpacket direction is sampled randomly according to Eq. (26). The sampled frequency is connected to the CMF and it must be transformed to the RF.
The probability of choosing a specific radiative recombination transition from those introduced by Eq. (21b) is defined similarly to the case of radiative deexcitation, namely (28)
The only difference between the treatment of radiative deexcitation and radiative recombination is the sampling of a new packet frequency. In the case of radiative recombination the frequency is sampled from the equation (29)
where z ϵ [0, 1] is chosen randomly, the emissivity is given by Eq. (B.22), and v_{i} is the ionization frequency. The method of sampling of a frequency is following: first the integral on the right hand side is calculated. Second, the left hand side integral is repeatedly evaluated as a function of its lower boundary until the integral value is approximatelly the same as the right hand side. The lower integration boundary value v^{↝} is a new rpacket frequency.
In the case of collisional deexcitation (Eq. (21c)) and recombination (Eq. (21d)), the ipacket transits into a kpacket. No determination of packet frequency and direction is needed in this case.
In the case of internal jumps only new actual state of macroatom is set. New transition probabilities (Eqs. (21) and (22)) are calculated in the current new state and the process is repeated.
Fig. 2 Scheme of choosing the corresponding process. There is saved information for concrete transitions (lower rectangle represents transition of a specific type with a definition (Eq. (25))). Number of these transitions can be very high. For saving a computer time, first of all it is determined which process occurs. All total transition probabilities for specific processes are determined with Eq. (23). After the process determination the programme goes through all possible transitions and chooses the right one which satisfies the condition described in the text. In the example here eight processes are possible, out of them the fourth process is chosen and it contains eleven transitions. The fifth transition is then selected. 
2.4 Kinetic energy packets
Kinetic energy packets represent free electrons in the medium. These packets do not contain any information about their status and they are assumed not to move in space. The only thing that has to be done is a choice of their interaction. The total collisional cooling rate is equal to (Kromer 2009) (30)
where 𝒩_{I} is the total number of ions (of element k), and finally 𝒩_{L} is the total number of levels (of ion j of element k). Individual processes leading to a change of a kpacket include collisional excitation (31)
where C_{iu} is given by Eq. (B.17), freefree emission (32)
where C_{0} = 1.426 × 10^{−27} in CGS units, q_{i} is the charge of ion i, and N_{i} is the ion i concentration, freebound transition (radiative recombination) (33)
the cross section can be found in the Eq. (B.21). The first term represents thermal and ionization energy converted to radiant energy, and the second term the ionization energy spontaneously converted to radiant energy in a given fb transition. Finally, collisional ionization (34)
where c_{im} is written in the Eq. (B.20). A process is chosen randomly, the same procedure is used as in the case of the ipackets: the rates are calculated, then a random number is generated and the same machinery as described in Fig. 2 is processed.
When the collisional excitation (Eq. (31)) or ionization (Eq. (34)) is chosen, the kpacket transforms to an ipacket. A corresponding transition (or ion) must be chosen, because ipacket has to contain information, which element (ion, level) will be set during the initialization.
Both freefree (Eq. (32)) and freebound (Eq. (33)) transitions change kpackets to rpackets. A new frequency of this packet is in the freefree case determined from the equation (35)
z ϵ [0, 1] is chosen randomly. In the case of freebound transition the equation is as the same as the Eq. (29).
3 The general code structure
The current version of the code calculates radiative transfer through a given medium. The code can handle a broad variety of media: stellar atmosphere, circumstellar disc, expanding atmosphere of supernova, etc. The main advantage of using the MC code is that we solve the problem in 3 D (or how many dimensions we want). Thus we are not restricted by the requirement of the model symmetry and radiative transfer through an atmosphere with irregular mass distribution can be calculated. On the other hand, the MC method is difficult to apply to optically very thick media. The calculation is extremely time expensive because the number of interactions of every photon (in this case packet) is very huge. This is the case of deep layers of the stellar photosphere. In our stellar wind case this problem does not appear.
The general code structure is in the flow chart Fig. 3. The main input of the code is a precalculated model of the structure of the flow. This part of the input includes density, temperature, and velocity profiles. The next part of the input includes information about the star: effective temperature, luminosity, stellar radius, and a lower boundary condition, which can be a photospheric spectrum calculated by another stellar atmosphere code. The main output of the MonteCarlo code is the emergent flux.
The code currently assumes the local thermodynamical equilibrium, but it is possible to use precalculated populations of energy levels from another code (e.g. PoWR, Hamann & Gräfener 2003, 2004, Sander et al. 2015).
There are two types of numerical grids used in our code: the model grid (hereafter modGrid) and the propagation grid (hereafter propGrid), see Appendix A. The modGrid represents the hydrodynamical model (input) in isolated points. Direct radiative transfer calculations in modGrid can be difficult and depend on modGrid geometry. To avoid this property, the propGrid is created. The propGrid is a Cartesian grid, which allows radiation transfer calculations regardless the input model geometry. Every propGrid cell is connected to a modGrid cell. We assume every physical quantity in the propGrid cell to be constant. Creating the propGrid is done only at the beginning of the calculations. The connection of modGrid and propGrid is described in the following subsection. All technical information about modGrid and propGrid is in Appendix A.
The basic propGrid is a regular shaped grid. The modGrid can be very irregular: with a larger number of grid points in some areas and less grid points otherwise. To avoid unnecessary dense propGrid we developed an adaptive grid. There are two types of grids: octgrid, see Fig. A.3, and onelevel subgrid, see Fig. A.4. Both grids are more deeply described in the Appendix A.
Fig. 3 Scheme describing the actual code routine. 
3.1 A connection between the model and the propagation grid
The basic propGrid is based on the modGrid size and the density of model grid points. The dimensions of the modGrid are equal to the Cartesian maximal dimensions of the propGrid, thus the modGrid fits to the rectangular parallelepiped shape of the propGrid.
During calculation of propagation of an energy packet through the propGrid, we need actual values of density, temperature, etc., which are stored in the modGrid. Since both grids are independent, the definition of their connection is necessary. Every propGrid cell (hereafter PC) must ‘know’ which modGrid point (modGrid line or modGrid surface for 2 D and 1 D models, respectively) it belongs to. The corresponding modGrid point (line/surface) is chosen using following steps:
Choose a propagation cell represented by its centre with coordinates (x_{0}, y_{0}, z_{0}).

Calculate distance l between this cell (i.e. its centre) and all possible modGrid points after the expression following from the metrics theory (a distance of two sets, see Munkres 2000) (36)
where (x, y, z) is a set of position vectors defining the whole set of modGrid points. In the 1D spherically symmetric case it is a surface of a sphere , where R_{i} is the sphere radius; in the 2D axially symmetric case the set is defined by circles parallel to the xy plane , where R_{i} is a circle radius and Z_{i} is the zcoordinate of the circle. In the 3D case the set of position vectors is a single point.
Select the modGrid point (line / surface) with the lowest distance.
In the 1 D case test whether the propGrid point is in the wind area, i. e. its centre radius r_{C} satisfies the condition
The physical variables are constant in the propGrid cell volume except of the velocity field, which must be continuous. Analytical velocity fields are used or an interpolation methods must be implemented in the case of discrete velocity fields.
3.2 Subtleties of the propagation of a packet through the adaptive grid
Indexes of all six neighbouring cells are stored for each propGrid cell, which simplifies treatment of energy packets moving through the propGrid. This concept is simple for a regular grid, however, it becomes more complicated in the case of an adaptive grid. If a packet flies from the current cell to another one, a problem can occur: the cell division is not regular, thus the packet can move from one grid level to another. This problem must be solved more sophistically.
Every cell ‘knows’ information about its connection to the lowerlevel cell (a cell which is divided by the current subgrid) it belongs to, and the first upperlevel cell (the first cell of a subgrid dividing the current cell). These indexes are zero if the upper or lower level grids do not exist. For the purpose of packet propagation it ‘knows’ also indexes of neighbouring cells of the same level appertaining to the same lowerlevel cell. In numbers, if the index of a neighbouring cell is larger than zero, then the cell exists and it is at the same level. Otherwise it is equal to zero or −99 in the case of the basic grid, which indicates the end of the computational domain; in the case of neighbour cell it indicates that the next cell must be found in a different tree branch, see Fig. 4. The index zero means that there is no other cell on the same level, the index −99 is the edge of the whole propagation grid. The packet propagation through the PC is calculated in the highest level cells. It cannot continue outside the subgrid but it has to find the neighbouring cell even in another subgrid.
An example of propagation is shown in Fig. 4. The packet starts in the highest level in the cell 15, it moves to the cell 16 because it is the same subgrid. But there is no cell beyond the current cell and the programme has to look one level lower – to the cell 12 which is divided by this subgrid. But the packet is on the border of this subgrid too, thus once more it goes to the lower level (cell 7). Here it finds an existing neighbouring cell 8. Furthermore, it determines an existence of a subcell – a higher level subgrid. And it moves to the cell 19, which is calculated using the Eqs. (A.1) and (A.2). This cell is the highest one and the propagation can continue. This process can be generalized also for subgrids with a variable number of subcells.
One can see that there are several options of propagation. The most simple case is the transition from one cell to another within the same subgrid. Then the Eq. (A.1) can be used. We know the actual cell, the index of the first cell of the current subgrid.
If the packet reaches the outer boundary of the current subcell there is no possibility to propagate in the same level. The programme must skip to the lower level and continue propagation. Only if the current cell has its subgrid we have to skip again upwards. The Eq. (A.1) is used for the index calculation. N_{0} is saved – it is a saved index in the lower level cell. For example the x+ border: we know, that the cells on this border have N_{x} = 1, the indexes N_{y}, N_{z} are calculated with the equation
where R is the coordinate of the x y z corner of the given subgrid.
Fig. 4 Scheme of the packet propagation through the grid. The bottom scheme shows a propagation grid. The upper scheme shows the same grid, but expressed in a tree structure. The green line is an example of packet’s path. It is simple to show this path on the lower graph. The upper graph explicitly shows the adaptive propagation grid level structure including the illustration of cell indexing. The numbers in the figure describe the logic of cell indexing in the code. 
4 Testing the propagation grids
To check the effectivity of propGrid refinements (construction of the adaptive propGrid) described in Appendix A.3, we constructed an artificial wind model for the purpose of testing. Included atomic data were imported from the Opacity project (Delahaye et al. 2016). The atomic levels were split according to their principal and orbital quantum numbers (see Table 1). We assumed a 1 D spherically symmetric wind with a homologous velocity field, (37)
where V_{∞} is the terminal wind velocity and R_{∞} is the radius of the outer model boundary, r is the radial distance. Within the wind we created a clump represented by a spherical shell with either one or two density peaks. Then the massdensity of the wind is described by (38)
This equation is a power law density distribution in the interval starting at R_{*} and ending at R_{∞}. The clump is represented by an additional function ϱ_{cl}(r). For a singledensitypeak spherical shell located between R_{1} and R_{2} this function is described by the Gaussian distribution (39)
The singledensitypeak clump with a Gaussian density distribution has the width R_{2} − R_{1} and the density maximum is at (R_{1} + R_{2})/2 (the shell centre). For a doubledensitypeak clump, the massdensity of the wind is described by (40)
with the halfwidth σ, and r_{01} and r_{02} are positions of densitypeak centres fulfilling the condition R_{1} < r_{01} < r_{02} < R_{2}. The function ϱ_{cl}(r) is equal to zero outside the region (R_{1}, R_{2}) both for single and doubledensitypeak profiles. The quantities ϱ_{0} and ϱ_{cl} are density parameters. Numerical values of these parameters for Eqs. (39) and (40) for our testcase are listed in the Table 2. The distance between modGrid points is 0.5R_{*} outside the clump and 0.01R_{*} inside the clump.
Let us introduce a set of variables which describe how well the particular propGrid (or adaptive propGrid) corresponds to the modGrid. There exists an association between modGrid and propGrid (see Sect. 3.1). For each modGrid cell I we calculate a ratio of the total volume of all propGrid cells associated with a given modGrid cell and the volume of that particular model cell VI. This ratio is expressed as (41)
where V_{PC} represents the propGrid cell volume belonging to the modGrid cell I (see Eq. (36), PC → I indexes propGrid cells belonging to the Ith model cell). Ideally, χ_{I}.= 1 However, this condition is not always satisfied since we try to fit modGrid cells of a general shape by a set of rectangular propGrid cells. To compare the global quality of connection between the propGrid and modGrid, we calculated the standard deviation as (42)
where N_{modGrid} is the total number of included model cells connected with at least one propagation cell. For a good correspondence between propGrid and modGrid, this number should be as close to zero as possible. The relative covering is a percentage of model cells with at least one propGrid cell connected to.
We compared three types of grids, namely a regular one (Appendix A.2), the octgrid (Appendix A.3.1), and a onelevel subgrid (Appendix A.3.2). Every grid type has special characteristics and behaves differently based on the number of cells.
Numbers of atomic levels included.
4.1 Regular propagation grid test
First we tested the regular grid. We calculated several spectra with different division of the basic grid (N^{B}, N^{B} = 5^{3}, 10^{3}, 20^{3}, 50^{3}, 100^{3}, 150^{3}) and plotted the calculated spectra in the region around the hydrogen Lya line (see the left panel of Fig. 5 and Table 3). Figure 5 shows a PCygni type Lya line with an additional absorption component caused by the spherical density clump. The position of this absorption peak at 1128 Å corresponds to the velocity of the clump (0.072c). It is nicely seen how the grid spatial resolution affects the resulting line profiles. For the coarsest resolution with N^{B} = 5^{3} the Lya line is very weak, and the clump absorption is almost not seen. For N^{B} = 10^{3} the PCygni line appears, however, the clump absorption peak is rather weak. It is stronger for N^{B} = 20^{3}. For finer grids we obtain better resolution. This is consistent with the parameter describing the relative covering, which is below 100% for low number of propGrid points (see the last column in Table 3). Finally, the profiles for N^{B} = 100^{3} and N^{B} = 150^{3}, are almost identical, which indicates that the grid with N^{B} = 100^{3} is sufficient to resolve of the spherical clump. The center and right panel in Fig. 5 include several blueshifted absorption lines (hydrogen Banner series and neutral helium lines), which originate from absorption in the spherical clump. The pure absorption lines with an exception of a weak Ha line are not present in the case with no clump (pink line). The Lya and Ha lines have PCygni profiles originating in the wind.
The calculated 〈χ^{2}〉 for given selected cases is shown in Table 3. The basic (regular) propGrids (with 10^{3} and 50^{3} cells in Table 3) do not correspond to the modGrid well. Increasing the number of propGrid cells to 100^{3} or 150^{3} fits the modGrid better (see Table 3). The space diagonal of the cell is equal 0.87(R_{2}−R_{1}) or 0.56(R_{2}−R_{1}) (R_{1}, R_{2} are introduced in Eq. (39)) in the case of the propGrid 100^{3} or 150^{3}, respectively. However, this improvement is done for the price of increasing the number of grid points also in the parts of the model, where it is not necessary.
Parameters of the regular grids for the singlepeak tests.
Fig. 5 Comparison of spectral line profiles calculated with different basic propGrids in the case of the single clump model. Left panel: hydrogen Lya line with a PCyg profile and the clump absorption component. Central and right panel: Hydrogen Balmer lines and neutral helium absorption lines. Line labels correspond to their vacuum wavelengths. 
4.2 Adaptive propagation grids tests
A significantly better and computationally cheaper fit to the model grid can be obtained using adaptive grids. We tested both methods described in Appendix A.3 for subgrid creation.
For the test of the octgrid we selected basic grids with 10^{3} and 50^{3} cells, both of them having insufficient relative covering. We created four propGrids with a different number of subcells. As a tool for subcell generation we used virtual points (see Appendix A.3). The more virtual points the more subcells are generated. We generated 10^{3}, 10^{4}, 10^{5}, and 10^{6} virtual points. Corresponding numbers of generated subcells are listed in Table 4. The numbers in last two columns of the Table characterize the quality of fit of the propGrid to the modGrid. With the exception of the 10^{3} basic grid with 10^{3} virtual points all propGrid show 100% covering. For basic 10^{3} propGrid and 10^{5} virtual points we obtain a propGrid with a standard deviation comparable to the 150^{3} regular propGrid, which need less than half of the CPU time for packet propagation. With 10^{6} virtual points we obtain significantly better standard deviation (0.06) with still lower consumption of the CPU time than for the 150^{3} basic propGrid. The results of calculations with a double densitypeak clump are shown in the Figs. 6 and 7. One can see that the double profile appears with a better resolution for propGrid with lower standard deviation. The double minima correspond to the double shaped density wind structure. We note that insufficient covering of the propGrid causes erroneous absorption profile with only one peak. Corresponding adaptive propGrids are plotted in Fig. 8.
Tests have shown us some advantages and disadvantages of the selected propagation grid types. The regular grid provides a good approximation, however, a sufficient resolution requires a dense propGrid in the whole volume. Much better way is to use the octgrid, which works well for less dense basic grids: the time necessary to create the octgrid is acceptable since it is done only once, at the beginning of the calculations, and the time of the packet propagation is not worse than the regular grid 150 × 150 × 150, but with a better resolution with a parameter equal to 0.06 and compared to the clump size, the space diagonal (of a rectangular prism) of the smallest propGrid cell is equal to 0.03σ, where σ s defined in Eq. (40). The octgrid with denser regular grid is not as efficient as the previous case, grid creating is more time consuming with the similar efficiency. The onelevel subgrid type consumes the most time for the grid creation and a large amount of memory as well.
The propagation cell input parameters must be chosen wisely. This can be seen in the Fig. 8. The right figure shows the onelevel subgrid based on the 10 × 10 × 10 regular propagation grid. The high density of the subgrid is also located in the large area outside the spherical clump. For the onelevel subgrid the choice of denser basic grid is more appropriate. One can see that in both cases the basic grid cells are divided even outside the clump, which is caused by the virtual points distribution. The virtual points are not distributed in the clump only but also among the other model cells.
The adaptive grid provides a more suitable approach for the modGrid representation. It does not slow the packet propagation, however its creation is really time consuming in some cases. This could be prevented by parallelization of the grid creation, which is a very difficult task, or enable reading already created propagation grid, if the old one is sufficient enough.
5 Physical processes in the outflow
The state of the matter (such as temperature, electron density, etc) must be determined before the packets propagation calculation. In this section we introduce the implemented approximations for the plasma state calculation and the approximation included in the radiative part of the code.
Computational statistic of the propagation grids in several cases.
Ionization and excitation state of plasma
The ionization equilibrium calculation plays a crucial role in radiative transfer. One of the best ways is to solve the statistical equilibrium equations. The ionization state and the occupation numbers for the atomic energy levels depend on the radiation field and the temperature. The implementation is not easy, thus we temporarily use a simpler approximation of LTE.
The LTE approximation presupposes the thermodynamic equilibrium for the ionization fraction and energy levels populations is satisfied locally (on small distance scales). Thus all occupancy numbers depend only on local values of temperature and electron density, via the Boltzmann Eq. (C.1) and the Saha Eq. (C.2). These equations do not provide reasonable approximation of excitation and ionization balance in hot star winds. The NLTE approximation is more appropriate.
The partition function is also needed to determine plasma state. Eq. (C.4) is used in the code as an initial approximation. Later on more appropriate approximation will be used.
6 Computational details
The energy flow is in the MC method represented via the packet transport. A packet is a MC energy quantum. The initial packet position is the stellar surface (the inner boundary condition). A packet propagates through the ejecta until: it reaches the inner boundary, it reaches outer boundary, or it interacts n_{max} times (a user defined integer). Spectrum is calculated only if the packet reaches the outer boundary.
Since the propagation of packets is independent it can be parallelized quite easily. Every processor propagates n packets and the spectrum is calculated using informations of all propagated packets. The code is parallelized with the MPI (mpich) library.
6.1 rpackets
Every time the rpacket is propagated the optical depth must be calculated. The optical depth is composed of line and continuum parts. The line optical depth is limited in frequency, and, consequently, in differentially expanding media also in space, by a narrow spectral profile. The continuum optical depth is nonzero for any frequency.
It is advantageous to create a list of all possible atomic transitions. Every line has its transition frequency and we arrange the list according to this frequency. The packet remembers then the last line it visited. The next line can be only the one closest to the last line. This saves the computational time.
Opacities in continua are more computationally expensive. The photoionization cross section depends on the CMF frequency. Since the cross section is not expressed as an analytical formula but stored as a table for specific values of frequency, interpolation is necessary. This process consumes plenty of time. User can choose an approximate way: choose number of included cross sections for each element. The number can be larger than zero (number of levels with included photoionization cross sections, other levels’ photoionization cross section is zero), zero (no photoionization) or –1, all possible data included.
Fig. 6 Comparison of spectra calculated with different propGrids in the case of the double clump density profile. Upper panels: with different number of cells in the basic propGrid. Middle panels: the octgrid: with different number of the virtual points, basic grid size is 10 × 10 × 10. Lower panels: the octgrid: with different number of virtual points, basic grid size is 50 × 50 × 50. 
Fig. 7 Comparison of spectra calculated with different propGrids in the case of the double peak density profile. Upper panels: the onelevel type of the propagation grid with the number of the virtual points as a parameter, 10 × 10 × 10 basic grid. Middle panels: onelevel grid type, 50 × 50 × 50 basic grid. Lower panels: onelevel grid type, 100 × 100 × 100 basic grid. 
6.2 ipackets
The ipacket processing – calculation of its rates and the inner transitions can be very time consuming. We have implemented a tool to make the computations faster. The most time consuming parts of calculation of the ipacket rates are the recombination parts, see Eq. (B.11): the integrals must be calculated every single step. Integrals contain exponentials which slow calculations, and, consequently, also the packet propagation down. We prefer to use the precalculated tables. These tables are created during the initialization processes. The function values are saved for every included ion with the photoionization cross section data are dependent on temperature. The intervals are split up linearly, thus the corresponding interval for the given temperature can be calculated with a simple interpolation.
Fig. 8 Adaptive propagation grids calculated for the doubleclump spherically symmetrical model. We have plotted only one quadrant in the zero plane. The red “circle” is the central star, blue circles represent the model grid shells and the grey and the blue Cartesian grid represent propagation grid. It is clear that in the area of the clump the density of the propGrid is much larger. Left: Octgrid, basic grid: 50^{3}, 10^{6} virtual points, right: onelevel subgrid, basic grid: 10^{3}, 10^{6} virtual points. 
6.3 kpackets
The main challenge of the kpackets implementation is in the calculation optimalization. Every time a packet changes into the kinetic energy all possible rates must be calculated, which takes some time, especially in the case with many ion levels included.
The first method to save some computation time is to order the rates from the largest to the smallest one. The time needed for the calculation which process was chosen gets lower – in the code the corresponding rates are saved in an array. The code generates a random number, which is multiplied with the total rate. Then it goes through the whole array and looks for an interval corresponding to the randomly calculated rate. If the first rates are the largest one, the most of the finding will end in the beginning and do not have to go through all indexes, which clearly saves the time.
The second method consists of saving the rates from one last model cell. If the packet returns to the previous model cell it is not necessary to calculate rates again.
6.4 propGrid resolution and number of packets
Our code contains many changeable parameters which are important for proper spectra calculation. This implies that inappropriately chosen parameters may cause incorrect output. We tested spectral sensitivity on the propGrid parameters and number of packets. We chose four propGrid sizes N^{B}: 10 × 10 × 10, 50 × 50 × 50, 100 × 100 × 100, and 150 × 150 × 150, and calculated emergent spectra with a different number of packets. Emergent spectra are plotted in Fig. 9. The case of the 10^{3} propGrid shows insufficient accuracy with too shallow spectral lines, while the other three grids yield sufficient accuracy. Different tests of the accuracy of the adaptive propGrids described in the Sect. 4 (Figs. 5–7) support the fact that the 10^{3} propGrid does not produce reliable emergent spectrum. What minimum size of propGrid is good enough depends also on the modGrid resolution. Spectra calculated in Fig. 9 are well converged in the case of propGrid N^{B} = 50^{3}, whilst the presence of a spherical clump in Sect. 4 needed much better resolution or the adaptive propGrid. If small structures are present (clumps, etc.) the propGrid cells must be small enough to describe them.
On the other hand, the number of packets affects the noise in the spectra, which is clearly visible in the plot in Fig. 9. Furthermore if a bad propGrid is chosen, increasing the number of packets does not cause better fit of spectra. Both propGrid and number of packets must be chosen properly. The relative difference among spectra are shown in Fig. 10. One can see that relative difference between two spectra decreases with increasing number of packets used for spectra calculation.
7 A comparison with the Tardis code
To check basic performance of the code and validate it before further steps, we performed comparison for a specific case using another code, which is somewhat similar. To this end we chose the wind model as simple as possible, consisting only from hydrogen and helium because of a low number of spectral lines. The wind was assumed to be spherically symmetric, radiation flux at the lower wind boundary was simplified as the Planck function, atomic level populations and ionization balance were calculated using the LTE approximation. The adopted velocity field was homologous (Eq. (37)). Input model parameters are summarized in Table 5. For comparison calculations we used the code Tardis^{2}. Our code reads the same input wind structure (chemical composition, mass density, velocity structure, temperature) and the atomic data as the Tardis code.
7.1 A brief description of the Tardis code
The Tardis code (Kerzendorf & Sim 2014) is a Monte Carlo code for radiative transfer modelling of a supernova explosion. This code is written in Python and C languages. It allows one to calculate a supernova spectrum for several parameters describing the physical state of outflow. It can calculate atomic level populations in both LTE and NLTE approximation, and it can also iterate the temperature structure. The final spectrum is denoised using the method of virtual packets.
Numerical grid. The Tardis model grid is, similarly to the wind model (see Appendix A.1), spherically symmetric and divided to several shells. Everything is calculated in the spherical coordinates and the code does not use a propagation grid.
Fig. 9 A comparison of calculated spectra for different regular propGrid sizes (parameters are written above every figure). Each plot contains a comparison of spectra generated from a different number of packets (see the legend of upper right plot). 
Fig. 10 Number of packets is one of the crucial factors for the magnitude of ratio: signal to noise in calculated spectra. Above: a comparison of calculated spectra for different number of packets. Below: relative changes between two spectra for each line. Spectra calculated for the grid N^{B}=50^{3}. 
7.2 Input model for comparison calculations
There are several parameters which have to be defined. We assumed elemental abundances to be homogeneous in the whole wind. The mass density dependence in Tardi s is defined by the formula (43)
where ϱ_{0} is density in the lowest boundary shell (where υ = υ_{0}) at an initial time t_{0}. In this equation t_{expl} is a time parameter and the exponent κ has to be chosen by the Tardi s user. The Eq. (43) is timedependent, but both Tardi s and our code calculate a single time frame. We selected t_{0} = 1 day, t_{expl} = 13 days, and κ = −2. The Tardi s code assumes a homologous velocity field (37). The excitation and ionization equilibrium is calculated using the LTE approximation.
Atomic data. The Tardis atomic data for hydrogen and helium are sourced from the Kurucz database. They differ from the atomic data in the Sect. 4. This database contains hydrogen atom with 24 energy levels, neutral helium: 48 energy levels, and once ionized helium: 24 energy levels. We used these data for our tests.
Fig. 11 Comparison between emergent spectra calculated by our MC code and by Tardis. Upper panels: spectra comparison, Lower panels: relative difference between the spectral fluxes. Left: Model with only resonance scattering included. Every time the packet interacts with an atom, it is directly emitted in the random direction with the same CMF frequency. Right: Downbranch mode: resonance scattering and fluorescence allowed. The deexcitation is possible also to an arbitrary lower atomic level. 
A set of parameters used for the comparison with the Tardis code.
7.3 Properties of the calculated models
The size of computational domain is defined by the velocity field, because of homologous approximation υ ∝ r. We calculated models with only line interactions included. Line profiles are approximated as a δ function, (44)
where i denotes the i–th line and v_{i} is its transition frequency. The line interaction in Tardis can be treated in three modes (see Kerzendorf & Sim 2014, Table 1), namely ‘scatter’, where all transitions are handled as resonance line scattering, ‘downbranch’, where internal transitions in macroatoms are not considered, and ‘macroatom’, where the full macroatom scheme is used.
The main output is the emergent spectrum averaged over all angles. We calculated spectra for all three Tardis line treatment modes (scatter, downbranch, macroatom) and compared the spectra between the Tardi s code and our code. The results for the resonance and downbranch cases are in Fig. 11 and for the macroatom case in Fig. 12. Strong Lyman and weaker Balmer lines are present in all figures. While differences in the continuum flux are negligible, differences in lines are more pronounced. It is seen especially in plots of relative differences (lower panels of Figs. 11 and 12).
There is a difference close to the wavelength corresponding to the wind boundary velocity V_{∞} in the blue wing of the hydrogen Lyα line. It is caused by a velocity shift between Tardi s and our models and would disappear if the spectra were corrected for it. The velocity shift causes also a difference near the centre of the Hα lines for all three test cases. Generally, the differences in the Hα line increase with more complex line treatment approximation. The largest difference appears in the macroatom case. There are several possibilities to explain these differences. The Tardis code and our code are not identical. The differences in emergent spectra can be caused by differences between grids used in both codes. Our code creates a Cartesian propGrid, which approximates the modGrid and ‘pixelizes’ the physical model. The Tardis code calculates only with one spherical grid corresponding to our modGrid.
To understand these differences better, we did several additional tests to see the spectra sensitivity on several parameters change. The most sensitive is a change depending on the model velocity structure. The change of V_{∞} and R_{∞} caused large shifts of the Lyα and the other Lyman lines. The Balmer lines similar to Ha were not that sensitive and their shifts were not as significant as those in the Lyman series. These changes reflect the optical depth of the outflow in particular line frequencies. We tested the lower boundary condition: the change of temperature changed the whole continuum and did not provide any better agreement with the Test case. The stellar radius change affects the effective temperature as well and did not explain any difference. The tests including the input supernova model change did not bring a significant change.
To check the sensitivity of our calculations to a grid resolution, we tested models calculated with our code for several different grid resolutions. Spectra do not differ significantly in the cases of regular grid N^{B} > 50.
We found differences also between the plasma states, the electron densities differ by 0.7% in maximum, and ionization fractions are also different by about 1%. This affects the opacity structure and, consequently, the Monte Carlo rates.
Fig. 12 Comparison between emergent spectra calculated by our MC code and by Tardis. Line interaction – the macroatom mode: the resonance scattering, fluorescence, and the internal jumps (in the same ionization state) are included in the model. Upper panel: spectra comparison. Lower panel: relative difference between the spectral fluxes. 
8 Summary and future work
We presented the first version of our new Monte Carlo code for stellar wind modelling, which currently performs the calculation of synthetic spectra of objects with an outflow. The MC method was based on Lucy’s Macroatom model (Lucy 2002, 2003). The Macroatom approach is a NLTE radiative transfer treatment. However, as we concentrated on the radiative transfer problem, the ionization balance and the level population calculation were evaluated using the Saha and the Boltzmann equation, because it is much more simple than the calculation using the statistical equilibrium equations.
The code starts with reading the input model and creates the model grid (modGrid), which contains physical properties of the circumstellar outflow. Now we use a precalculated hydrodynamical (or pure analytical) model of stellar outflow. Then it creates a computation domain where the propagation of radiative energy packets (which describe radiation) is processed. This domain is called propGrid and is Cartesian. Every packet represents an energy which propagates and changes its form until it reaches the outer boundary. The inner and outer boundaries are set up to be consistent with the modGrid.
The propGrid is a Cartesian grid composed of cells of a block shape in general. The propGrid may be regular or adaptive. The parameters of this grid can affect the output, thus we tested the spectra calculation based on the propGrid parameters. The test of the adaptive grid possibilities was done using an artificial ‘single clump’ and ‘double clump’ density model. The model without such clump can be calculated efficiently using a regular propGrid. On the other hand, the model with an added clump requires an enhanced grid cell density in the clump and a sparser (i.e. unchanged) grid cell density elsewhere. In the case of the regular propGrid, the grid density is the same in the whole computation domain. It has to be as dense as the modGrid in the region of its highest resolution. This construction can be very memory consuming. Based on this fact we created an adaptive propGrid, which has a high grid density only in regions with highest modGrid resolution and is sparse elsewhere.
To test the results of radiative transfer, we compared the synthetic spectra calculated with our code and with the Tardis code (Kerzendorf & Sim 2014) for simple test cases with limited photonatom interactions. We adopted three kinds of approximations of line interaction treatment. There are differences between resulting emergent spectra, increasing with the considered complexity of line interaction. The reason of this difference is most probably the difference in treatment of Monte Carlo radiative transfer between our code and Tardis.
The code can solve the radiative transfer through the circumstellar outflow with an analytically described velocity field. As a next step we shall develop NLTE ionization and excitation balance in full 3 D models with a velocity field defined only in discrete points. This will be reported in future papers.
Acknowledgements
This research has made use of the NASA’s Astrophysics Data System Abstract Service. J.F. was partly supported by the Erasmus+ mobility project. The Astronomical Institute Ondřejov is supported by the project RVO:67985815. This work was partly supported by a grant GA ČR 2234467S. Computational resources were supplied by the project “eInfrastruktura CZ” (eINFRA CZ LM2018140) supported by the Ministry of Education, Youth and Sports of the Czech Republic.
Appendix A Numerical grids
Appendix A.1 Model grid
The model grid (hereafter also modGrid) is a set of points in ndimensional space (n ∈ {1, 2, 3}), where all physical characteristics of a given circumstellar environment are defined. Selected physical quantities such as temperature are recalculated during model iterations.
The type of the modGrid corresponds to the assumed problem symmetry. For the most general 3 D case the modGrid is simply a set of evenly distributed isolated grid points. For a 2 D model the modGrid is represented by evenly distributed grid lines and for a 1 D model the modGrid is represented by evenly distributed grid surfaces. For example, spherically symmetric models are 1 D and the grid points are characterized only by radial distance from the grid centre. Consequently, a grid point is a surface of a sphere with a given radius.
The initial physical model may be calculated by another (e.g. hydrodynamic) code, but also input from previous run of the current code is possible. As the initial model does not need to follow the rules of our modGrid (since it may come from a different numerical code), a transformation of the input model to our modGrid definition is necessary in such case.
Appendix A.2 Basic propagation grid
The energy packet propagation itself is not calculated in the modGrid directly. Following Kromer & Sim (2009), a Cartesian grid called the ‘propagation grid’ (hereafter propGrid) is created and all energy packet propagation calculations are processed in this grid. It is advantageous to locate the star to the central part of the propGrid, then the stellar centre is located at the grid origin.
The main advantage of the Cartesian grid construction is that we employ the relative simplicity of performing calculations of models with variable symmetries in a Cartesian grid (which is rectangular) compared to using more general coordinates (spherical, cylindrical, or even more general ones) to lower the amount of necessary calculations. In addition, if we want to include a model with a specific general geometry, we have to modify only the part of the code which reads the input model and the part which connects the model and propGrid. The energy packet propagation itself remains unchanged.
The basic propGrid has to be large enough to include the whole modGrid and consists of cells, where represents a number of cells for the coordinate i (i ∈ {x, y, z}), see Fig. A.1. All cells have widths , , and , respectively. The widths may be different for each coordinate. The cell orientation (the plus and minus direction) is in the Fig. A.2. The corresponding cell’s index N(n_{x}(x), n_{y}(y), n_{z}(z)) is defined as (A.1)
(see numbering in Fig. A.1), N_{0} is the index of the first cell, in the case of the basic propagation grid N_{0} = 1. This equation can be also used for the determination of neighbouring cells. For example, a neighbour of the cell N in the ydirection placed at (n_{x}, n_{y} + 1, n_{z}) has an index . The indexes of neighbouring cells are calculated and stored during the propGrid creation.
Coordinates of a cell corresponding to an energy packet at a position described by the radius vector r = (r_{x}, r_{y}, r_{z}) can be calculated as (A.2)
where ⌊⌋ denotes the floor of a number, w_{i} is the cell width, and N_{i} is the number of cells for the coordinate i, i ∈ {x, y, z}. To prevent possible numerical problems in cases when r points exactly to the cell boundary, we add a small number ε, which can be the smallest number allowed by machine accuracy. This helps to overcome the situation, when the sum should be exactly an integer, while its numerical representation is slightly lower than the expected value. Then the expression gives the correct value of the cell coordinate n_{i}.
Fig. A.1 Scheme of the basic propGrid and its coordinates. For simplicity, only two dimensions, coordinates x and y, are plotted here. Each cell has widths and . Numbers below the figure and on the left denote and for particular cells, respectively. The cell indexes following from Eq. (A.1) are written in the centers of selected cells. Generalization to three dimensions is straightforward. 
Fig. A.2 Propagation cell and its orientation to the coordinate system. The boundaries are signed as follows: the blue x+, the red x−, the green y+, the yellow y−, the front z−, and the back z+. The colour cube was created using the code from https://tex.stackexchange.com. 
Appendix A.3 Adaptive propagation grid
The regular propagation grid is not the most efficient one for every case. As an example, let us consider a circumstellar disc without a wind. Such disk has usually a conical shape. To calculate a reliable model we need to choose the propGrid with small cells, which in the case of a regular grid would be also outside the disc where very little matter is present. To use memory more efficiently we need to introduce a propagation grid in which the almost empty space outside the disk is not considered in such detail.
Every propagation cell represents a homogeneous area of the outflow numerical model. We want the propGrid to correspond to the modGrid as accurately as possible. The best and the simplest case is a modGrid with regularly distributed grid points (every point’s location can be represented as r = n_{x}e_{x} + n_{y}e_{y} +n_{z}e_{z}, where n_{x,y,Z} were defined in Eq. (A.l) and e_{x,y,z} is the Cartesian orthonormal basis). In this case, the propGrid can be created as a regular net with the model points located in the center of each propagation cell.
In a general case, the modGrid points need not to be distributed equally. In some regions the number of modGrid points in a unit volume can be much larger than in the rest of the space. A natural solution handling this problem is to create a denser propGrid only in the region with more modGrid points, and not in the whole computational domain. The concept of the adaptive propagation grid is described, for example, in Lunttila & Juvela (2012).
The adaptive propagation grid (hereafter adaptive propGrid) is created in the following way: We start from the basic propGrid described in the Section A.2. Every grid cell is tested if a subgrid has to be created. For a 3 D modGrid, modGrid grid points are counted in every basic propGrid cell. If their number is larger than one, then a subgrid is created. Depending on the subgrid type (described later), this process can be recursively repeated for subgrid cells, if necessary.
Due to an inherent symmetry, modGrid points in 1D and 2 D cases form grid surfaces and grid lines, respectively. Consequently, subgrids can not be created such as in the 3 D case, a method for their creation has to be developed first. We define socalled virtual points for this purpose. The virtual point is a point with coordinates defined by modGrid surface or line (2 or 1 coordinates, respectively), supplemented by randomly chosen remaining coordinates. For example, in the 1D case the radius is defined but the rest two coordinates, the azimuthal and longitudinal ones, must be calculated. In total, N virtual points are generated. The number of virtual points 𝒩_{i} belonging to the given model cell is proportional to the spherical shell area (A.3)
where N is the total number of virtual points, r_{i} represents the radial distance of the ith shell and N_{MC} the total number of model cells. Its azimuthal and inclination angle are chosen randomly. Virtual points are counted for every basic propGrid cell. Similarly to the 3D case, if this number is larger than one, then a subgrid is created, again with the possibility of recursively generated subgrids.
Fig. A.3 Adaptive propGrid ‘octgrid’. Every cell can be recursively divided into four (in a 2D case) or eight (in the 3D case) cells until the predefined minimum cell limit is reached. 
Fig. A.4 Adaptive propGrid onelevel subgrid. In every basic cell (red colour) a subcell with number of cells n_{x},n_{y},n_{z} and variable widths w_{x},w_{y}, and w_{z} is created. 
The main advantage is that the method of virtual points can be used for a wide range of geometries without complicated calculations. Virtual points are created before setting the propGrid up and are not used after the propGrid is created, thus they do not occupy the computer memory during the rest of the calculation.
Adaptive propGrid indexing. The created subgrid is referred to as a one level higher grid than the basic one. Every cell of the subgrid has saved the index of its parent cell (lowerlevel cell) and the divided cell (the parent cell) in the lower level has saved the first cell index of the subgrid (upperlevel cell). If there is no upper/lower cell, the corresponding index is equal to zero.
The indexing is illustrated in Fig. 4. First the basic grid is numbered as described in the Appendix A.2. Then the programme goes through the grid from the cell with the lowest index and if a cell contains a subgrid, the subgrid is indexed. In the next step it recursively checkes the subgrid and looks for its existing subgrids, again starting with the cells with lowest indexes. The process continues in a loop until it reaches the end of the given subgrid level. Then it goes back to the lowerlevel cell and continues to the cell with the index higher by one. In following subsections implemented methods of subgrid creation are described.
Appendix A.3.1 The octgrid
The first subgrid scheme is based on the octree grid. It is called ‘octgrid’ (e.g. Jonsson 2006; Acreman et al. 2010; Robitaille 2011). It is created based on the following rules:
A cell is divided into 2^{3} subcells, which are in a geometrical sense similar with the divided cell. The subcell dimensions are half the dimensions of the original cell.
Subcells can be created recursively, the total number of nested cells is limited by a minimal cell width.
A 2 D example is shown in the Fig. A.3. This subgrid scheme is useful for models with nonregularly distributed grid points as well as for models with inhomogeneities of model physical parameters or their derivatives. Such models require larger number of propGrid cells in specific regions and subgrid nesting allows smooth transition between smallest and largest cells.
This type of grid has been already described in Barnes & Hut (1986) and Kurosawa & Hillier (2001). However, while Kurosawa & Hillier used calculation of an emissivity at random points, we created virtual points.
Appendix A.3.2 Onelevel subgrid
The second type of grid we call the ‘onelevel subgrid’ (see Fig. A.4). In this case the basic cell is divided into a regular subgrid, generally with any number of equal subcells. The number of subcells may differ from the number of the basic grid cells. The grid has the following properties:
A cell is divided into a subgrid where is a number of subcells in the coordinate i.
Further division of the subgrid is forbidden.
Appendix B Radiative and collisional rates used in the code
Here we specify detailed expressions for rates used in our calculations.
Appendix B.1 Radiative rates
Appendix B.1.1 Boundbound rates
The radiative rate for an upward boundbound transition from the level i to the level u (excitation) is equal to (B.1)
The inverse rate (deexcitation from the level i to the level l) is (B.2)
In these equations, A_{il}, B_{il}, and B_{il} represent the Einstein coefficients, ϕ_{iu}(v) is the line profile of a transition between levels i and u, and J(v) the mean radiation field intensity.
For the case of the Sobolev approximation we approximate the line profile by a delta function, (B.3)
where v_{iu} denotes the line centre frequency. Then the absorption rate for transitions from the level i to the level u simplifies to (for derivation see Klein & Castor 1978) (B.4)
where the upper index S stands for Sobolev, (B.5)
denotes a probability that a scattered photon escapes. The quantity is the Sobolev optical depth for the transition between levels i and u. Similarly, the emission rate from the level i to the level l is (B.7)
where β_{il} is defined similarly to the Eq. (B.6).
If we treat stimulated emission as negative absorption the absorption and emission rates become (B.8)
respectively.
Appendix B.1.2 Boundfree and freebound rates
The radiative ionization rate from the level i to the level p (here the ground level of the next higher ion) (B.10)
and the radiative recombination rate (i is the ground level of an ion here, the corresponding population is denoted as ), m corresponds to a level in the next lower ion: (B.11)
where α_{mi}(v) is the photoionization cross section from the level m to the level i. The equation consists of two terms: spontaneous and stimulated recombination rate, (B.12)
is the SahaBoltzmann factor.
If we assume that the stimulated recombination is negative photoionization, the photoionization rate is given by (B.15)
Appendix B.2 Collisional rates
For collisional rates from the level i to the level f (which may be l, u, m, or p, see section 2) the expressions of C_{if} for boundbound and boundfree transitions are given below. The inverse rates are calculated using (see Hubeny & Mihalas 2015, Eq. 9.52).
Appendix B.2.1 Boundbound transitions
The boundbound collisional rates are calculated using the van Regemorter (1962) approximation (B.17)
(see also Hubeny & Mihalas 2015, Eq. 9.58), here c_{0} = 5.46510. 10^{−11} is a constant, I_{H} = 13.6eV is ionization potential of hydrogen and hv_{iu} = ϵ_{u} – ϵ_{i} is the transition energy between levels i and u. The function Γ is defined as (B.18)
is given by (here n, n′’ denote principal quantum number, l, l′ the orbital quantum number) (B.19)
and E_{1} (x) is the exponential integral function.
Appendix B.2.2 Boundfree transitions
The collisional boundfree rates are calculated using an approximate formula (see Hubeny & Mihalas 2015, Eq. 9.60) (B.20)
Appendix B.2.3 Freebound transitions
The cross section for the spontaneous freebound transition from the level i to the level m is (B.21)
where the SahaBoltzmann factor is defined in (B.14).
Appendix B.3 Emissivities
The freebound emissivity (B.22)
where α_{ip}(v) is the photoionization cross section for a transition from the level i to the level p.
The freefree emissivity (B.23)
N_{j} represents the population of the ion j, the freefree cross section, see (Kromer 2009, Eq. (3.46)) (B.24)
If we put the Eq. (B.23) into Eq. (35) we can find easily an analytic formula (B.25)
where z ∈ (0, 1) is a random number.
Appendix C Ionisation and excitation equilibrium
The occupation numbers are in the LTE approximation calculated via the Boltzmann excitation formula, which we use in the form (C.1)
where n denotes the atomic level number densities, ɡ is the statistical weight, ε is the level energy, and the indexes l and u denote lower and upper level of the transition, respectively. The ionization balance follows the Saha equation, here in the form (C.2)
where C is the Saha constant (C.3)
N denotes ion number densities, U denotes the partition function, ε denotes energies, and the indexes p and j denote atomic states. The partition function of the ion i is (C.4)
where N_{L} is the total number of included levels for ion i. The sum runs only over the levels included in the calculation.
References
 Abbott, D. C., Bieging, J. H., & Churchwell, E. 1981, ApJ, 250, 645 [Google Scholar]
 Acreman, D. M., Harries, T. J., & Rundle, D. A. 2010, MNRAS, 403, 1143 [NASA ADS] [CrossRef] [Google Scholar]
 Adam, J. 1990, A&A, 240, 541 [NASA ADS] [Google Scholar]
 Asplund, M., & Lind, K. 2010, in IAU Symposium, 268, Light Elements in the Universe, eds. C. Charbonnel, M. Tosi, F. Primas, & C. Chiappini, 191 [Google Scholar]
 Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., Kallman, T. R., Fryxell, B. A., & Taam, R. E. 1990, ApJ, 356, 591 [Google Scholar]
 Carciofi, A. C., & Bjorkman, J. E. 2006, ApJ, 639, 1081 [Google Scholar]
 Carciofi, A. C., Bjorkman, J. E., & Magalhães, A. M. 2004, ApJ, 604, 238 [NASA ADS] [CrossRef] [Google Scholar]
 Carciofi, A. C., Bjorkman, J. E., & Zsargó, J. 2017, in The Lives and Death Throes of Massive Stars, 329, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, 390 [NASA ADS] [Google Scholar]
 Carlberg, R. G. 1980, ApJ, 241, 1131 [Google Scholar]
 Castor, J. I. 1974, ApJ, 189, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
 Delahaye, F., Zwölf, C. M., Zeippen, C. J., & Mendoza, C. 2016, J. Quant. Spec. Radiat. Transf., 171, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2003, A&A, 406, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Owocki, S. P. 2005, A&A, 437, 657 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Driessen, F. A., Kee, N. D., & Sundqvist, J. O. 2021, A&A, 656, A131 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Feldmeier, A., Puls, J., Reile, C., et al. 1995, Ap&SS, 233, 293 [Google Scholar]
 Feldmeier, A., Puls, J., & Pauldrach, A. W. A. 1997, A&A, 1997, 878 [Google Scholar]
 Feldmeier, A., Shlosman, I., & Vitello, P. 1999, ApJ, 526, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Höfner, S., & Liljegren, S. 2019, IAU Symp., 343, 9 [NASA ADS] [Google Scholar]
 Georgiev, L. N., Hillier, D. J., & Zsargó, J. 2006, A&A, 458, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R., & Gräfener, G. 2003, A&A, 410, 993 [CrossRef] [EDP Sciences] [Google Scholar]
 Hamann, W.R., & Gräfener, G. 2004, A&A, 427, 697 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2018, A&A, 616, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennicker, L., Puls, J., Kee, N. D., & Sundqvist, J. O. 2020, A&A, 633, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Higginbottom, N., Knigge, C., Long, K. S., Sim, S. A., & Matthews, J. H. 2013, MNRAS, 436, 1390 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J. 1987, ApJS, 63, 969 [NASA ADS] [Google Scholar]
 Hillier, D. J. 1990, A&A, 231, 116 [NASA ADS] [Google Scholar]
 Hillier, D. J. 2012, in IAU Symposium, 282, From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards, & I. Hubeny, 229 [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2015, Theory of Stellar Atmospheres (Princeton, N.J.: Princeton University Press) [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jonsson, P. 2006, MNRAS, 372, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Kerzendorf, W. E., & Sim, S. A. 2014, MNRAS, 440, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Kerzendorf, W., Sim, S., Vogl, C., et al. 2022, 1S.5281/zenodo.6662839 [Google Scholar]
 Klein, R. I., & Castor, J. I. 1978, ApJ, 220, 902 [NASA ADS] [CrossRef] [Google Scholar]
 Korcáková, D., & Kubát, J. 2005, A&A, 440, 715 [CrossRef] [EDP Sciences] [Google Scholar]
 Kromer, M. 2009, PhD thesis, Technische Universität München, Germany [Google Scholar]
 Kromer, M., & Sim, S. A. 2009, MNRAS, 398, 1809 [Google Scholar]
 KrtiCka, J., & Kubát, J. 2017, A&A, 606, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kurosawa, R., & Hillier, D. J. 2001, A&A, 379, 336 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leenaarts, J., & Carlsson, M. 2009, in ASP Conf. Ser., 415, The Second Hinode Science Meeting: Beyond DiscoveryToward Understanding, eds. B. Lites, M. Cheung, T. Magara, J. Mariska, & K. Reeves, 87 [Google Scholar]
 Lobel, A., & Blomme, R. 2008, ApJ, 678, 408 [NASA ADS] [CrossRef] [Google Scholar]
 Long, K. S., & Knigge, C. 2002, ApJ, 579, 725 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 2002, A&A, 384, 725 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2003, A&A, 403, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2005, A&A, 429, 19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B., & Solomon, P. M. 1970, ApJ, 1970, 879 [NASA ADS] [CrossRef] [Google Scholar]
 Ludwig, H. G., & Steffen, M. 2016, Astron. Nachr., 337, 844 [NASA ADS] [CrossRef] [Google Scholar]
 Lunttila, T., & Juvela, M. 2012, A&A, 544, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Milne, E. A. 1926, MNRAS, 86, 459 [NASA ADS] [Google Scholar]
 Munkres, J. R. 2000, Topology, 2nd edn. (Upper Saddle River, NJ: Prentice Hall) [Google Scholar]
 Noebauer, U. M., & Sim, S. A. 2015, MNRAS, 453, 3120 [Google Scholar]
 Noebauer, U. M., & Sim, S. A. 2019, Living Rev. Comput. Astrophys., 5, 1 [CrossRef] [Google Scholar]
 Nordlund, Å., & Stein, R. F. 2009, in AIP Conf. Ser., 1171, Recent Directions in Astrophysical Quantitative Spectroscopy and Radiation Hydrodynamics, eds. I. Hubeny, J. M. Stone, K. MacGregor, & K. Werner, 242 [Google Scholar]
 Oskinova, L. M., Feldmeier, A., & Hamann, W. R. 2004, A&A, 422, 675 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Oskinova, L. M., Feldmeier, A., & Hamann, W. R. 2006, MNRAS, 372, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Cramner, S. R. 1994, ApJ, 1994, 887 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Puls, J. 2002, ApJ, 568, 965 [NASA ADS] [CrossRef] [Google Scholar]
 Owocki, S. P., & Rybicki, G. B. 1984, ApJ, 284, 337 [Google Scholar]
 Owocki, S. P., Castor, J. I., & Rybicki, G. B. 1988, ApJ, 335, 914 [NASA ADS] [CrossRef] [Google Scholar]
 Papkalla, R. 1995, A&A, 295, 551 [NASA ADS] [Google Scholar]
 Pauldrach, A. W. A., Hoffmann, T. L., & Lennon, M. 2001, A&A, 375, 161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Petrenz, P., & Puls, J. 2000, A&A, 2000, 956 [Google Scholar]
 Puls, J., Owocki, S. P., & Fullerton, A. W. 1993, A&A, 279, 457 [NASA ADS] [Google Scholar]
 Puls, J., Urbaneja, M. A., Venero, R., et al. 2005, A&A, 435, 669 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Puls, J., Najarro, F., Sundqvist, J. O., & Sen, K. 2020, A&A, 642, A172 [EDP Sciences] [Google Scholar]
 Robitaille, T. P. 2011, A&A, 536, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Runacres, M. C., & Owocki, S. P. 2002, A&A, 381, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A. A. C. 2017, in The Lives and DeathThroes of Massive Stars, 329, eds. J. J. Eldridge, J. C. Bray, L. A. S. McClelland, & L. Xiao, 215 [NASA ADS] [Google Scholar]
 Sander, A., Shenar, T., Hainich, R., et al. 2015, A&A, 577, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sander, A. A. C., Hamann, W. R., Todt, H., Hainich, R., & Shenar, T. 2017, A&A, 603, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 SantolayaRey, A. E., Puls, J., & Herrero, A. 1997, A&A, 323, 488 [NASA ADS] [Google Scholar]
 Štepán, J., & Trujillo Bueno, J. 2013, A&A, 557, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., & Puls, J. 2018, A&A, 619, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., & Feldmeier, A. 2010, A&A, 510, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Puls, J., Feldmeier, A., & Owocki, S. P. 2011, A&A, 528, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Owocki, S. P., & Puls, J. 2018, A&A, 611, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sundqvist, J. O., Björklund, R., Puls, J., & Najarro, F. 2019, A&A, 632, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Šurlan, B., Hamann, W. R., Kubát, J., Oskinova, L. M., & Feldmeier, A. 2012, A&A, 541, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Šurlan, B., Hamann, W. R., Aret, A., et al. 2013, A&A, 559, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 UdDoula, A. 2015, IAU Symp., 307, 321 [NASA ADS] [Google Scholar]
 van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Vogl, C., Sim, S. A., Noebauer, U. M., Kerzendorf, W. E., & Hillebrandt, W. 2019, A&A, 621, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Whitney, B. A. 2011, Bull. Astron. Soc. India, 39, 101 [NASA ADS] [Google Scholar]
 Zsargó, J., Hillier, D. J., & Georgiev, L. N. 2006, A&A, 447, 1093 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Detailed expressions for radiative (R) and collisional (C) rates can be found in the Appendix B.
All Tables
All Figures
Fig. 1 Scheme of the macroatom interactions. 

In the text 
Fig. 2 Scheme of choosing the corresponding process. There is saved information for concrete transitions (lower rectangle represents transition of a specific type with a definition (Eq. (25))). Number of these transitions can be very high. For saving a computer time, first of all it is determined which process occurs. All total transition probabilities for specific processes are determined with Eq. (23). After the process determination the programme goes through all possible transitions and chooses the right one which satisfies the condition described in the text. In the example here eight processes are possible, out of them the fourth process is chosen and it contains eleven transitions. The fifth transition is then selected. 

In the text 
Fig. 3 Scheme describing the actual code routine. 

In the text 
Fig. 4 Scheme of the packet propagation through the grid. The bottom scheme shows a propagation grid. The upper scheme shows the same grid, but expressed in a tree structure. The green line is an example of packet’s path. It is simple to show this path on the lower graph. The upper graph explicitly shows the adaptive propagation grid level structure including the illustration of cell indexing. The numbers in the figure describe the logic of cell indexing in the code. 

In the text 
Fig. 5 Comparison of spectral line profiles calculated with different basic propGrids in the case of the single clump model. Left panel: hydrogen Lya line with a PCyg profile and the clump absorption component. Central and right panel: Hydrogen Balmer lines and neutral helium absorption lines. Line labels correspond to their vacuum wavelengths. 

In the text 
Fig. 6 Comparison of spectra calculated with different propGrids in the case of the double clump density profile. Upper panels: with different number of cells in the basic propGrid. Middle panels: the octgrid: with different number of the virtual points, basic grid size is 10 × 10 × 10. Lower panels: the octgrid: with different number of virtual points, basic grid size is 50 × 50 × 50. 

In the text 
Fig. 7 Comparison of spectra calculated with different propGrids in the case of the double peak density profile. Upper panels: the onelevel type of the propagation grid with the number of the virtual points as a parameter, 10 × 10 × 10 basic grid. Middle panels: onelevel grid type, 50 × 50 × 50 basic grid. Lower panels: onelevel grid type, 100 × 100 × 100 basic grid. 

In the text 
Fig. 8 Adaptive propagation grids calculated for the doubleclump spherically symmetrical model. We have plotted only one quadrant in the zero plane. The red “circle” is the central star, blue circles represent the model grid shells and the grey and the blue Cartesian grid represent propagation grid. It is clear that in the area of the clump the density of the propGrid is much larger. Left: Octgrid, basic grid: 50^{3}, 10^{6} virtual points, right: onelevel subgrid, basic grid: 10^{3}, 10^{6} virtual points. 

In the text 
Fig. 9 A comparison of calculated spectra for different regular propGrid sizes (parameters are written above every figure). Each plot contains a comparison of spectra generated from a different number of packets (see the legend of upper right plot). 

In the text 
Fig. 10 Number of packets is one of the crucial factors for the magnitude of ratio: signal to noise in calculated spectra. Above: a comparison of calculated spectra for different number of packets. Below: relative changes between two spectra for each line. Spectra calculated for the grid N^{B}=50^{3}. 

In the text 
Fig. 11 Comparison between emergent spectra calculated by our MC code and by Tardis. Upper panels: spectra comparison, Lower panels: relative difference between the spectral fluxes. Left: Model with only resonance scattering included. Every time the packet interacts with an atom, it is directly emitted in the random direction with the same CMF frequency. Right: Downbranch mode: resonance scattering and fluorescence allowed. The deexcitation is possible also to an arbitrary lower atomic level. 

In the text 
Fig. 12 Comparison between emergent spectra calculated by our MC code and by Tardis. Line interaction – the macroatom mode: the resonance scattering, fluorescence, and the internal jumps (in the same ionization state) are included in the model. Upper panel: spectra comparison. Lower panel: relative difference between the spectral fluxes. 

In the text 
Fig. A.1 Scheme of the basic propGrid and its coordinates. For simplicity, only two dimensions, coordinates x and y, are plotted here. Each cell has widths and . Numbers below the figure and on the left denote and for particular cells, respectively. The cell indexes following from Eq. (A.1) are written in the centers of selected cells. Generalization to three dimensions is straightforward. 

In the text 
Fig. A.2 Propagation cell and its orientation to the coordinate system. The boundaries are signed as follows: the blue x+, the red x−, the green y+, the yellow y−, the front z−, and the back z+. The colour cube was created using the code from https://tex.stackexchange.com. 

In the text 
Fig. A.3 Adaptive propGrid ‘octgrid’. Every cell can be recursively divided into four (in a 2D case) or eight (in the 3D case) cells until the predefined minimum cell limit is reached. 

In the text 
Fig. A.4 Adaptive propGrid onelevel subgrid. In every basic cell (red colour) a subcell with number of cells n_{x},n_{y},n_{z} and variable widths w_{x},w_{y}, and w_{z} is created. 

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.