© The Authors 2023

1 Introduction

The stellar wind is a kind of an outflow from the star that is stable on a long-term scale. The stellar radiation field is often strong enough to accelerate matter and create a stellar wind. Intensive radiation-driven stellar winds are observed, for example, around O-type stars and Wolf-Rayet 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 ‘mass-loss 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); line-driven winds were later described in detail by Lucy & Solomon (1970). The first hydrodynamical solution of the line-driven 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 thermody-namic 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. Santolaya-Rey 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 (ud-Doula 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 hydro-dynamical 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 radiation-hydrodynamics 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 multi-dimensional 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 pseudo-2 D models, later on followed by a pseudo-3 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 non-symmetric 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 time-dependent spectra to take the rapid variations of circum-stellar 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 macro-atom 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 (r-packets), internal energy packets (i-packets), and kinetic energy packets (k-packets). The packets may change their form in radiation-matter interactions.

The matter is quantized into macro-atoms. Detailed properties of energy packets and macro-atoms are described in Lucy (2002).

2.1 Radiative packets

Radiation in this method is quantized to radiation energy packets (r-packets), 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 radiation-matter interaction (see Fig. 1) to packets of internal energy (i-packets) or kinetic (thermal) energy (k-packets).

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 r-packets 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 $$r_{\text{init}}=\left(\right(R_{\ast }\sin \left(\right(\text{Θ}\left)\right)\cdot \cos \left(\right(\text{Φ}\left)\right),R_{\ast }\sin \left(\right(\text{Θ}\left)\right)\cdot \sin \left(\right(\text{Φ}\left)\right),R_{\ast }\cos \left(\right(\text{Θ}\left)\right)\left)\right),$$(1a) $$\text{Φ}=2\pi \tau ,$$(1b) $$\cos \left(\right(\text{Θ}\left)\right)=2\sigma -1.$$(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 r-packet is randomly chosen in the outward direction from the stellar surface: another random unit vector $$n_{\text{init.}}=\left(\right(\sin \left(\right(t\left)\right)\cdot \cos \left(\right(p\left)\right),\sin \left(\right(t\left)\right)\cdot \sin \left(\right(p\left)\right),\cos \left(\right(t\left)\right)\left)\right)$$(2a) $$p=2\pi {\tau }^{\prime },$$(2b) $${\cos }^2\left(\right(t\left)\right)={\sigma }^{\prime }.$$(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 rotation-matrix multiplication), then the initial packet direction d_pack is in the form $$\begin{array}{l}{d}_{\text{pack}}=\left(\right(r_3\cdot \sin \left(\right(t\left)\right)\cdot \cos \left(\right(p\left)\right)+r_1\cdot \cos \left(\right(t\left)\right)\cdot \cos \left(\right(p\left)\right)-r_2\cdot \sin \left(\right(p\left)\right),)\hfill \\ r_3\cdot \sin \left(\right(t\left)\right)\cdot \sin \left(\right(p\left)\right)+r_1\cdot \cos \left(\right(t\left)\right)\cdot \sin \left(\right(p\left)\right)+r_2\cdot \cos \left(\right(p\left)\right),\hfill \\ r_3\cdot \cos \left(\right(t\left)\right)-r_1\cdot \sin \left(\right(t\left)\right)\left)\right),\hfill \end{array}$$(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 $$E=\frac{L_{\text{star}}}{N},$$(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 (r-packets) 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 co-moving frame (hereafter CMF) using the index ↝. The interaction of matter and radiation can occur only at a point where the co-moving frame photon frequency v^↝equals the line transition frequency (v^↝ = v_line).

Consider an r-packet 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 $$v^{\rightsquigarrow }=v^{\#}\left(\right(1-\frac{n\cdot \upsilon }{c}\left)\right),$$(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 $${\tau }_{\text{cont}}={\displaystyle \underset{r}{\overset{r_{\text{R}}}{\int }}\text{d}{r}^{\prime }}\cdot {\chi }_{\text{cont}}\left(\right(r^{\prime }\left)\right),$$(6)

where r is the packet position, r_R is the position of the closest possible line interaction (resonance point) of the r-packet with matter and χ_cont (r′) is the continuum opacity, $${\chi }_{\text{cont}}\left(\right(r,v^{\rightsquigarrow }\left)\right)={\chi }^{\text{Th}}\left(\right(r\left)\right)+{\displaystyle \sum _{i=1}^{N_{\text{I}}}{\chi }_i^{\text{ff}}\left(\right(r,v^{\rightsquigarrow }\left)\right)+}{\displaystyle \sum _{i=1}^{N_{\text{L}}}{\chi }_i^{\text{bf}}\left(\right(r,v^{\rightsquigarrow }\left)\right)}.$$(7)

The first term corresponds to Thomson scattering, the second one to free-free processes, and the third one to bound-free 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₀ is given by (see Kromer 2009, Eq. (4.15)) $${\tau }_{\text{line}}=n_l\frac{B_{lu}hv}{4\pi }\left(\right(1-\frac{n_u{g}_l}{n_l{g}_u}\left)\right){\displaystyle \underset{v^{\rightsquigarrow }\left(\right(s=0\left)\right)}{\overset{v^{\rightsquigarrow }\left(\right(s=s_0\left)\right)}{\int }}\text{d}{v}^{\rightsquigarrow }\varphi \left(\right(v^{\rightsquigarrow }\left)\right)\frac{\text{d}s}{\text{d}{v}^{\rightsquigarrow }},}$$(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 $$\begin{array}{l}{\tau }_{\text{line}}=n_l\frac{B_{lu}hv}{4\pi }\left(\right(1-\frac{n_u{g}_l}{n_l{g}_u}\left)\right)\hfill \\ \times {\left(\right(\frac{\text{d}s}{\text{d}{v}^{\rightsquigarrow }}\left)\right)}_{v_{lu}}\times \left\{\right\{\begin{array}{ll}0\hfill & v_{lu}\notin \left[\right[v^{\rightsquigarrow }\left(\right(0\left)\right),v^{\rightsquigarrow }\left(\right(s_0\left)\right)\left]\right],\hfill \\ -1\hfill & v_{lu}\in \left[\right[v^{\rightsquigarrow }\left(\right(0\left)\right),v^{\rightsquigarrow }\left(\right(s_0\left)\right)\left]\right],\hfill \end{array}\hfill \end{array}$$(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 $${\tau }_{\text{line}}=\frac{\kappa \rho c}{v_{lu}}{\left[\right[{\mu }^2\frac{\text{d}\upsilon }{\text{d}r}+\left(\right(1-{\mu }^2\left)\right)\frac{\upsilon }{r}\left]\right]}^{-1},$$(10)

where µ is angle cosine between n and υ(r), and $$\kappa \rho =\frac{\pi e^2}{m_{\text{e}}c}{f}_{lu}{n}_l\left(\right(1-\frac{n_u}{n_l}\frac{g_l}{g_u}\left)\right),$$(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 r-packet is allowed to travel a random optical depth τ_rand and then it interacts. This optical depth is (see, e.g. Whitney 2011) $$\begin{array}{cc}{\tau }_{\text{rand}}=-\ln \left(\right(1-z\left)\right)=-\ln \left(\right(z^{\prime }\left)\right),& z,z^{\prime }\in \left[\right[0,1\left]\right],\end{array}$$(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 $${\tau }_{\text{rand}}\ge {\tau }_{\text{cont}}+{\tau }_{\text{line}},$$(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, bound-free and free-free 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 free-free process can transform an r-packet into thermal kinetic energy (a k-packet). Bound-free processes contribute to both kinetic and internal energy (k- and i-packets, respectively). Following Lucy (2003, Eq. (27)), the packet converts to an i-packet with a probability w_i = v_i/v^↝, where v_i is the ionization edge frequency. The packet is transformed to the k-packet otherwise (w_k = 1 − v_i/v^↝).

Free-free processes transform an r-packet into a k-packet and vice versa. Bound-free processes contribute to both internal energy and kinetic energy. For each r-packet we decide (again using random numbers) whether it changes to a k-packet or to an i-packet. The line event changes the r-packet to an i-packet (macro-atom activation).

2.2 Macro-atom statistical equilibrium equations

Let us write the kinetic (statistical) equilibrium equation for the level i (see Hubeny & Mihalas 2015, Chapter 14) $${\displaystyle \sum _{j\ne i}^{\text{NL}}{n}_i\left(\right(R_{ij}+C_{ij}\left)\right)=}{\displaystyle \sum _{j\ne i}^{\text{NL}}{n}_j\left(\right(R_{ji}+C_{ji}\left)\right)},$$(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, $$\begin{array}{l}{\displaystyle \sum _{l=1}^{i-1}{n}_i\left(\right(R_{il}+C_{il}\left)\right)+}{\displaystyle \sum _{u=i+1}^{\text{NL}}{n}_i\left(\right(R_{iu}+C_{iu}\left)\right)}\hfill \\ ={\displaystyle \sum _{l=1}^{i-1}{n}_l\left(\right(R_{li}+C_{li}\left)\right)+}{\displaystyle \sum _{u=i+1}^{\text{NL}}{n}_u\left(\right(R_{ui}+C_{ui}\left)\right).}\hfill \end{array}$$(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)): $${\displaystyle \sum _{l=1}^{i-1}\left(\right({\mathcal{P}}_{li}-{\mathcal{P}}_{il}\left)\right)+}{\displaystyle \sum _{u=i+1}^{\text{NL}}\left(\right({\mathcal{P}}_{ui}-{\mathcal{P}}_{iu}\left)\right)}=0.$$(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: ${\dot{E}}_i^{\text{R}}$ (radiative energy emission rates) and ${\dot{E}}_i^{\text{C}}$ (col-lisional energy emission rates), and the absorption energy rates ${\dot{A}}_i^{\text{R}}$ (radiative) and ${\dot{A}}_i^{\text{C}}$ (collisional). These terms express energy absorbed or emitted per unit time and unit volume. The radiative energy rates can be expressed as $$\begin{array}{cc}{\dot{A}}_i^{\text{R}}={\displaystyle \sum _{l=1}^{i-1}{n}_l{R}_{li}\left(\right({\epsilon }_i-{\epsilon }_l\left)\right),}& {\dot{E}}_i^{\text{R}}={\displaystyle \sum _{l=1}^{i-1}{n}_i{R}_{il}\left(\right({\epsilon }_i-{\epsilon }_l\left)\right),}\end{array}$$(17a)

the collisional energy rates as $$\begin{array}{cc}{\dot{A}}_i^{\text{C}}={\displaystyle \sum _{l=1}^{i-1}{n}_l{C}_{li}\left(\right({\epsilon }_i-{\epsilon }_l\left)\right),}& {\dot{E}}_i^{\text{C}}={\displaystyle \sum _{l=1}^{i-1}{n}_i{C}_{il}\left(\right({\epsilon }_i-{\epsilon }_l\left)\right).}\end{array}$$(17b)

Let us express the sum of these rates using total transition rates Ƥ, $${\dot{A}}_i^{\text{R}}+{\dot{A}}_i^{\text{C}}={\displaystyle \sum _{l=1}^{i-1}{\mathcal{P}}_{li}{\epsilon }_i}-{\displaystyle \sum _{l=1}^{i-1}{\mathcal{P}}_{li}{\epsilon }_l}$$(18a) $${\dot{E}}_i^{\text{R}}+{\dot{E}}_i^{\text{C}}={\displaystyle \sum _{l=1}^{i-1}{\mathcal{P}}_{il}{\epsilon }_i}-{\displaystyle \sum _{l=1}^{i-1}{\mathcal{P}}_{il}{\epsilon }_l.}$$(18b)

Multiplying Eq. (16) with ɛ_i and substituting for ${\displaystyle \sum {}_{l=1}^{i-1}}\left({\mathcal{P}}_{li}-{\mathcal{P}}_{il}\right){\epsilon }_i$, from Eq. (18) we get $$\begin{array}{l}{\dot{E}}_i^{\text{R}}+{\dot{E}}_i^{\text{C}}+{\displaystyle \sum _{u=i+1}^{\text{NL}}{\mathcal{P}}_{iu}{\epsilon }_i}+{\displaystyle \sum _{l=1}^{i-1}{\mathcal{P}}_{il}{\epsilon }_l}\hfill \\ ={\dot{A}}_i^{\text{R}}+{\dot{A}}_i^{\text{C}}+{\displaystyle \sum _{l=1}^{i-1}{\mathcal{P}}_{li}{\epsilon }_i}+{\displaystyle \sum _{u=i+1}^{\text{NL}}{\mathcal{P}}_{ui}{\epsilon }_i,}\hfill \end{array}$$(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 macro-atom with discrete energy.

2.3 Internal energy packets

Absorption of the radiation energy packet (the r-packet) via excitation or ionization activates the macro-atom. Internal energy packets (i-packets) represent excitation energy in excited atoms. Every i-packet 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 i-packet 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 $i_{I_k{L}_0}$ 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)) $${\dot{E}}_i^{\text{tot}}={\dot{E}}_i^R+{\dot{E}}_i^C+{\dot{E}}_i^{\mathrm{int},down}+{\dot{E}}_i^{\mathrm{int},up},$$(20)

where ${\dot{E}}_i^{\text{R}}$ and ${\dot{E}}_i^{\text{C}}$ denote radiative and collisional deactivation (introduced in Eq.(17)), respectively, and ${\dot{E}}_i^{\text{int,\hspace{0.17em}up}}$ and ${\dot{E}}_i^{\text{int,\hspace{0.17em}down}}$ mean internal upward and downward jumps, respectively.

The macro-atom 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 macro-atom is deactivated by one of the non-internal (deactivation) processes, namely¹ radiative deexcitation from the state i to the state l – the probability of this process is equal to: $$p_{il}^{\text{R}}=n_i{R}_{il}^{\text{S}}\left({\epsilon }_i-{\epsilon }_l\right)/{\dot{E}}_i^{\text{tot}},$$(21a)

where $R_{il}^{\text{S}}$ is a rate in the Sobolev approximation, denoted with an upper index S, radiative recombination (considered only from an ion ground level): $$p_{im}^{\text{R}}=n_i{R}_{im}\left({\epsilon }_{i_{{\text{L}}_0}}-{\epsilon }_m\right)/{\dot{E}}_i^{\text{tot}},$$(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: $$p_{il}^{\text{C}}=n_i{C}_{il}\left({\epsilon }_i-{\epsilon }_l\right)/{\dot{E}}_i^{\text{tot}},$$(21c)

C_il is the collisional deexcitation rate (see Appendix B.2), and collisional recombination (considered only from an ion ground level): $$p_{im}^{\text{C}}=n_i{C}_{im}\left({\epsilon }_{i_{{\text{L}}_0}}-{\epsilon }_m\right)/{\dot{E}}_i^{\text{tot}},$$(21d)

where C_im is the collisional recombination rate (see Appendix B.2).

The internal processes (which do not deactivate the macro-atom) taken into account include internal downward jump within the current ionization state: $$p_{il}^{\mathrm{int},\text{down}}=n_i\left(R_{il}^{\text{S}}-C_{il}\right){\epsilon }_l/{\dot{E}}_i^{\text{tot}},$$(22a)

internal downward jump to the lower ionization state: $$p_{im}^{\mathrm{int},\text{down}}=n_i\left(R_{im}+C_{im}\right){\epsilon }_m/{\dot{E}}_i^{\text{tot}},$$(22b)

internal upward jump within the current ionization state: $$p_{iu}^{\mathrm{int},up}=n_i\left(R_{iu}^{\text{S}}+C_{iu}\right){\epsilon }_i/{\dot{E}}_i^{\text{tot}},$$(22c)

where $R_{iu}^{\text{S}}$ is the rate of transition to the upper level, and, finally, internal upward jump to a higher ionization state: $$p_{ip}^{\mathrm{int},up}=n_i\left(R_{ip}+C_{ip}\right){\epsilon }_i/{\dot{E}}_i^{\text{tot}}.$$(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 i-th level as $$L_i^{\text{type}}={\displaystyle \sum _{\text{j}}{p}_{ij}^{\text{type}},}$$(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 $${ℒ}_k={\displaystyle \sum _{t\text{ype}=1}^k{L}_i^{\text{type}}}$$(24)

used in the upper panel of Fig. 2. Clearly ℒ₀ = 0. We can also define a cumulative probability for a transition $${ℒ}_{i,x}={ℒ}_{i-1}+{\displaystyle \sum _{y=1}^x{p}_{iy}^{\text{type}}.}$$(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) 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 two-step process is faster than a one-step 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 i-packet transits into the r-packet. We have to calculate its new frequency and direction. The direction unit vector is sampled before the calculation of rates using the equation $$n=\left(\text{sin}\text{Θ})\right)\cdot \text{cos}\left(\text{Φ}\right),\cdot \text{sin}\left(\text{Θ}\right)\cdot \text{sin}\left(\text{Φ}\right),\text{cos}\left(\text{Θ})\right),$$(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 non-homological 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 l-th line is given by $$w_{il}^{\text{RD}}=\frac{p_{il}^{\text{R}}}{{\displaystyle {\sum }_{l^{\prime }}{p}_{i{l}^{\prime }}^{\text{R}}}}.$$(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 r-packet. Currently a δ-function is assumed as a line profile. The r-packet 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 $w_{im}^{\text{RR}}$ from those introduced by Eq. (21b) is defined similarly to the case of radiative deexcitation, namely $$w_{im}^{\text{RR}}=\frac{p_{im}^{\text{R}}}{{\displaystyle {\sum }_{m^{\prime }}{p}_{i{m}^{\prime }}^{\text{R}}}}.$$(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 $${\displaystyle \underset{v^{\rightsquigarrow }}{\overset{\infty }{\int }}\text{d}v{\left({\eta }_{im}^{\text{fb}}\left(v\right)\right)}^{\rightsquigarrow }=z{\displaystyle \underset{v_i}{\overset{\infty }{\int }}\text{d}v{\left({\eta }_{im}^{\text{fb}}\left(v\right)\right)}^{\rightsquigarrow },}}$$(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 r-packet frequency.

In the case of collisional deexcitation (Eq. (21c)) and recombination (Eq. (21d)), the i-packet transits into a k-packet. No determination of packet frequency and direction is needed in this case.

In the case of internal jumps only new actual state of macro-atom 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) $$c={\displaystyle \sum _{i=1}^{{𝒩}_{\text{I}}}{n}_i{C}_i^{\text{ff}}+}{\displaystyle \sum _{i=1}^{{𝒩}_{\text{L}}}\left(n_i{C}_{im}^{\text{fb,sp}}+n_i{C}_{ip}^{\text{ion}}+{\displaystyle \sum _{i^{\prime }=i=1}^{{𝒩}_{\text{L}}}{n}_i{C}_{i{i}^{\prime }}^{\text{exc}}}\right),}$$(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 k-packet include collisional excitation $$C_{iu}^{\text{exc}}=n_i{C}_{iu}\left({\epsilon }_u-{\epsilon }_i\right),$$(31)

where C_iu is given by Eq. (B.17), free-free emission $$C_i^{\text{ff}}=C_0{q}_i^2{T}_e^{1/2}{N}_i{n}_e,$$(32)

where C₀ = 1.426 × 10⁻²⁷ in CGS units, q_i is the charge of ion i, and N_i is the ion i concentration, free-bound transition (radiative recombination) $$C_{mi}^{\text{fb,sp}}=N_m{n}_{\text{e}},\left({\alpha }_i^{\text{E,spont}}-{\alpha }_i^{\text{spont}}\right)\left({\epsilon }_{m_{L_0}}-{\epsilon }_i\right),$$(33)

the cross section ${\alpha }_i^{\text{E,spont}}$ 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 f-b transition. Finally, collisional ionization $$C_{im}^{\text{ion}}=c_{im}\left({\epsilon }_{m_{L_0}}-{\epsilon }_i\right),$$(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 i-packets: 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 k-packet transforms to an i-packet. A corresponding transition (or ion) must be chosen, because i-packet has to contain information, which element (ion, level) will be set during the initialization.

Both free-free (Eq. (32)) and free-bound (Eq. (33)) transitions change k-packets to r-packets. A new frequency of this packet is in the free-free case determined from the equation $${\displaystyle \underset{v^{\rightsquigarrow }}{\overset{\infty }{\int }}\text{d}v{\left({\eta }_{j_{{\text{I}}_j}}^{\text{ff}}\left(v\right)\right)}^{\rightsquigarrow }=z{\displaystyle \underset{0}{\overset{\infty }{\int }}\text{d}v{\left({\eta }_{j_{{\text{I}}_j}}^{\text{ff}}\left(v\right)\right)}^{\rightsquigarrow },}}$$(35)

z ϵ [0, 1] is chosen randomly. In the case of free-bound 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 photo-spheric spectrum calculated by another stellar atmosphere code. The main output of the Monte-Carlo 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 (here-after 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 one-level 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:

1. Choose a propagation cell represented by its centre with coordinates (x₀, y₀, z₀).

2. 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) $$l=\min \left(\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2+{\left(z-z_0\right)}^2}\right),$$(36)

   where (x, y, z) is a set of position vectors defining the whole set of modGrid points. In the 1-D spherically symmetric case it is a surface of a sphere $\left\{x^2+y^2+z^2=R_i^2\right\}$, where R_i is the sphere radius; in the 2-D axially symmetric case the set is defined by circles parallel to the xy plane $\left\{x^2+y^2=R_i^2\wedge z=Z_i\right\}$, where R_i is a circle radius and Z_i is the z-coordinate of the circle. In the 3-D case the set of position vectors is a single point.

3. Select the modGrid point (line / surface) with the lowest distance.

4. 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

$$R_{\ast }<r_C<R_{\infty }.$$

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 prop-Grid 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 lower-level cell (a cell which is divided by the current subgrid) it belongs to, and the first upper-level 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 lower-level 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 sub-cell 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₀ 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 $$N_i=\lfloor \frac{r_i-R_i}{{\omega }_i}\rfloor +1$$

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, $$\upsilon \left(r\right)=V_{\infty }\frac{r}{R_{\infty }},$$(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 mass-density of the wind is described by $$\varrho \left(\right(r\left)\right)=\left\{\right\{\begin{array}{ll}{\varrho }_0{\left(\right(\frac{r}{R_{\ast }}\left)\right)}^{-2}\hfill & R_1\ge r>R_{\ast },\hfill \\ \begin{array}{cc}{\varrho }_0{\left(\right(\frac{r}{R_{\ast }}\left)\right)}^{-2}& {\varrho }_{\text{cl}}\left(\right(r\left)\right)\end{array}\hfill & R_2\ge r>R_1,\hfill \\ {\varrho }_0{\left(\right(\frac{r}{R_{\ast }}\left)\right)}^{-2}\hfill & R_{\infty }\ge r>R_2,\hfill \\ 0\hfill & r\ge R_{\infty }.\hfill \end{array}$$(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 single-density-peak spherical shell located between R₁ and R₂ this function is described by the Gaussian distribution $${\varrho }_{\text{cl}}\left(\right(r\left)\right)={\varrho }_{\text{cl}}\exp \left[\right[-\frac{{\left(\right(r-\frac{R_1+R_2}{2}\left)\right)}^2}{{\left(\right(R_2-R_1\left)\right)}^2}\left]\right].$$(39)

The single-density-peak clump with a Gaussian density distribution has the width R₂ − R₁ and the density maximum is at (R₁ + R₂)/2 (the shell centre). For a double-density-peak clump, the mass-density of the wind is described by $${\varrho }_{\text{cl}}\left(\right(r\left)\right)={\varrho }_{\text{cl}}\left[\right[\exp \left(\right(-\frac{{\left(\right(r-r_{01}\left)\right)}^2}{{\sigma }^2}\left)\right)+\exp \left(\right(-\frac{{\left(\right(r-r_{02}\left)\right)}^2}{{\sigma }^2}\left)\right)\left]\right]$$(40)

with the half-width σ, and r₀₁ and r₀₂ are positions of density-peak centres fulfilling the condition R₁ < r₀₁ < r₀₂ < R₂. The function ϱ_cl(r) is equal to zero outside the region (R₁, R₂) both for single- and double-density-peak profiles. The quantities ϱ₀ 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 $${\chi }_I=\frac{{\displaystyle \sum _{\text{PC}\to I}{V}_{\text{PC}}}}{V_{\text{I}}},$$(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 I-th 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 $$\sqrt{\langle {\chi }^2\rangle }=\sqrt{\frac{1}{N_{\mathrm{modG}\text{rid}}}{\displaystyle \sum _{I=1}^{N_{\text{modGrid}}}{\left({\chi }_I-1\right)}^2},}$$(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 one-level subgrid (Appendix A.3.2). Every grid type has special characteristics and behaves differently based on the number of cells.

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³, 10³, 20³, 50³, 100³, 150³) 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 P-Cygni 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³ the Lya line is very weak, and the clump absorption is almost not seen. For N^B = 10³ the P-Cygni line appears, however, the clump absorption peak is rather weak. It is stronger for N^B = 20³. 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³ and N^B = 150³, are almost identical, which indicates that the grid with N^B = 100³ is sufficient to resolve of the spherical clump. The center and right panel in Fig. 5 include several blue-shifted 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 P-Cygni profiles originating in the wind.

The calculated 〈χ²〉 for given selected cases is shown in Table 3. The basic (regular) propGrids (with 10³ and 50³ cells in Table 3) do not correspond to the modGrid well. Increasing the number of propGrid cells to 100³ or 150³ fits the modGrid better (see Table 3). The space diagonal of the cell is equal 0.87(R₂−R₁) or 0.56(R₂−R₁) (R₁, R₂ are introduced in Eq. (39)) in the case of the propGrid 100³ or 150³, 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.

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 P-Cyg 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³ and 50³ 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 sub-cells are generated. We generated 10³, 10⁴, 10⁵, and 10⁶ 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³ basic grid with 10³ virtual points all propGrid show 100% covering. For basic 10³ propGrid and 10⁵ virtual points we obtain a propGrid with a standard deviation comparable to the 150³ regular propGrid, which need less than half of the CPU time for packet propagation. With 10⁶ virtual points we obtain significantly better standard deviation (0.06) with still lower consumption of the CPU time than for the 150³ basic propGrid. The results of calculations with a double density-peak 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 $\sqrt{\langle {\chi }^2\rangle }$ equal to 0.06 and compared to the clump size, the space diagonal (of a rectangular prism) of the smallest prop-Grid 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 one-level 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 one-level 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 one-level 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.

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 r-packets

Every time the r-packet 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 one-level type of the propagation grid with the number of the virtual points as a parameter, 10 × 10 × 10 basic grid. Middle panels: one-level grid type, 50 × 50 × 50 basic grid. Lower panels: one-level grid type, 100 × 100 × 100 basic grid.

6.2 i-packets

The i-packet 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 i-packet 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 pre-calculated 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 double-clump 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: 503, 106 virtual points, right: one-level subgrid, basic grid: 103, 106 virtual points.

6.3 k-packets

The main challenge of the k-packets 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³ prop-Grid 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³ 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³, whilst the presence of a spherical clump in Sect. 4 needed much better resolution or the adaptive prop-Grid. 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². 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 NB=503.

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 $$\varrho \left(\upsilon ,t_0,t_{\text{expt}}\right)={\varrho }_0\left(\frac{t_0}{t_{\text{expt}}}\right){\left(\frac{\upsilon }{{\upsilon }_0}\right)}^K,$$(43)

where ϱ₀ is density in the lowest boundary shell (where υ = υ₀) at an initial time t₀. 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 time-dependent, but both Tardi s and our code calculate a single time frame. We selected t₀ = 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.

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, $${\varphi }_i\left(v\right)\propto \delta \left(v-v_i\right),$$(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, ‘down-branch’, 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 macro-atom 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 Macro-atom model (Lucy 2002, 2003). The Macro-atom 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 hydro-dynamical (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 photon-atom 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 cir-cumstellar 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 22-34467S. Computational resources were supplied by the project “e-Infrastruktura CZ” (e-INFRA 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 n-dimensional 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 mod-Grid 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 $N_x^{\text{B}}\times N_y^{\text{B}}\times N_z^{\text{B}}$ cells, where $N_i^B$ represents a number of cells for the coordinate i (i ∈ {x, y, z}), see Fig. A.1. All cells have widths $w_{\text{x}}^{\text{B}}$, $w_{\text{y}}^{\text{B}}$, and $w_{\text{z}}^{\text{B}}$, 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 $$\begin{array}{l}{N}^B(n_{\text{x}}\left(x\right),n_{\text{y}}\left(y\right),n_{\text{z}}\left(z\right))=\hfill \\ \left(N_0-1\right)+n_x+\left(n_y-1\right)\cdot N_x^{\text{B}}+\left(n_z-1\right)\cdot N_x^{\text{B}}\cdot N_y^{\text{B}},\hfill \end{array}$$(A.1)

$N^B\left(n_{\text{x}}\left(x\right),n_{\text{y}}\left(y\right),n_{\text{z}}\left(z\right)\right)\in \left\{1,\dots ,N_x^{\text{B}},N_y^{\text{B}},N_z^{\text{B}}\right\}$ (see numbering in Fig. A.1), N₀ is the index of the first cell, in the case of the basic propagation grid N₀ = 1. This equation can be also used for the determination of neighbouring cells. For example, a neighbour of the cell N in the y-direction placed at (n_x, n_y + 1, n_z) has an index $N+N_x^{\text{B}}$. 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 $$n_i=\lfloor \frac{r_i}{w_i}+\frac{N_i^{\text{B}}}{2}\rfloor +1$$(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 $r_i/w_i+N_i^{\text{B}}/2$ should be exactly an integer, while its numerical representation is slightly lower than the expected value. Then the expression $\lfloor r_i/w_i+N_i^{\text{B}}/2+\epsilon \rfloor$ 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 $w_x^{\text{B}}$ and $w_y^{\text{B}}$. Numbers below the figure and on the left denote $n_x^{\text{B}}$ and $n_y^{\text{B}}$ 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_xe_x + n_ye_y +n_ze_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 so-called 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 1-D 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 $${𝒩}_i=N_{\text{points}}\lceil \frac{r_i^2}{{\sum }_{l=1}^{N_{\text{MC}}}{r}_i^2}\rceil ,$$(A.3)

where N is the total number of virtual points, r_i represents the radial distance of the i-th 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 3-D 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 2-D case) or eight (in the 3-D case) cells until the predefined minimum cell limit is reached.

Fig. A.4 Adaptive propGrid one-level subgrid. In every basic cell (red colour) a subcell with number of cells nx,ny,nz and variable widths wx,wy, and wz 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 (lower-level cell) and the divided cell (the parent cell) in the lower level has saved the first cell index of the subgrid (upper-level 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 lower-level 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:

1. A cell is divided into 2³ subcells, which are in a geometrical sense similar with the divided cell. The subcell dimensions are half the dimensions of the original cell.

2. 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 non-regularly 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 One-level subgrid

The second type of grid we call the ‘one-level subgrid’ (see Fig. A.4). In this case the basic cell is divided into a regular sub-grid, 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:

1. A cell is divided into a subgrid $n_{\text{x}}^{\text{S}}\times n_{\text{y}}^{\text{S}}\times n_{\text{z}}^{\text{S}}$ where $n_i^{\text{S}}$ is a number of subcells in the coordinate i.

2. 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 Bound-bound rates

The radiative rate for an upward bound-bound transition from the level i to the level u (excitation) is equal to $$R_{iu}=B_{iu}{\displaystyle \int \text{d}v{\varphi }_{iu}\left(v\right)J\left(v\right).}$$(B.1)

The inverse rate (deexcitation from the level i to the level l) is $$R_{il}=A_{il}+B_{il}{\displaystyle \int \text{d}v{\varphi }_{il}\left(v\right)J\left(v\right).}$$(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, $$\varphi \left(v\right)=\delta \left(v-v_{iu}\right)$$(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) $$R_{iu}^{\text{s}}={\beta }_{iu}{\beta }_{iu}{J}_{iu},$$(B.4)

where the upper index S stands for Sobolev, $$J_{iu}={\displaystyle \int \text{d}v\cdot \varphi \left(v\right)J\left(v\right)={\displaystyle \int \text{d}v\cdot \delta \left(v-v_{iu}\right)J\left(v\right)=J\left(v_{iu}\right),}}$$(B.5)

and $${\beta }_{iu}=\frac{1}{{\tau }_{iu}^{\text{s}}}\left[1-\exp \left(-{\tau }_{iu}^{\text{s}}\right)\right]$$(B.6)

denotes a probability that a scattered photon escapes. The quantity ${\tau }_{iu}^{\text{s}}$ 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 $$R_{il}^{\text{S}}={\beta }_{il}\left(A_{il}+B_{il}{J}_{il}\right).$$(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 $${\tilde{R}}_{iu}^{\text{S}}=\left(B_{iu}-\frac{n_u}{n_i}{B}_{ui}\right){\beta }_{iu}{J}_{iu},$$(B.8)

and $${\tilde{R}}_{il}^{\text{S}}={\beta }_{il}{A}_{il},$$(B.9)

respectively.

Appendix B.1.2 Bound-free and free-bound rates

The radiative ionization rate from the level i to the level p (here the ground level of the next higher ion) $$R_{ip}=4\pi {\displaystyle \underset{v_i}{\overset{\infty }{\int }}\text{d}v\frac{{\alpha }_{ip}\left(v\right)}{hv}J\left(v\right),}$$(B.10)

and the radiative recombination rate (i is the ground level of an ion here, the corresponding population is denoted as $n_{i_{l_0}}$), m corresponds to a level in the next lower ion: $${\begin{array}{l}{R}_{im}={\left(\frac{n_m}{n_{i_{l_0}}}\right)}^{\ast }\hfill \\ \times 4\pi {\displaystyle \underset{v_m}{\overset{\infty }{\int }}\text{d}v\frac{{\alpha }_{mi}\left(v\right)}{hv}\left(\frac{2h{v}^3}{c^2}+J\left(v\right)\right)\exp \left(-\frac{hv}{k_{\text{B}}{T}_{\text{e}}}\right),}\hfill \end{array}}^{\ast }$$(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 $\left(R_{im}^{\text{spont}}\right)$ and stimulated $\left(R_{im}^{\text{stim}}\right)$ recombination rate, $$R_{im}^{\text{spont}}=\frac{n_{i_{l_0}}}{n_i}{n}_{\text{e}}{\varphi }_i\left(T_{\text{e}}\right)4\pi {\displaystyle \underset{v_i}{\overset{\infty }{\int }}\text{d}v\frac{{\alpha }_{mi}\left(v\right)}{hv}\frac{2h{v}^3}{c^2}\exp \left(-\frac{hv}{k_{\text{B}}{T}_{\text{e}}}\right),}$$(B.12)

and $$R_{im}^{\text{stim}}=\frac{n_{m_{{\text{L}}_0}}}{n_i}{n}_{\text{e}}{\varphi }_i\left(T_{\text{e}}\right)4\pi {\displaystyle \underset{v_i}{\overset{\infty }{\int }}\text{d}v\frac{{\alpha }_{mi}\left(v\right)}{hv}J\left(v\right)\exp \left(-\frac{hv}{k_{\text{B}}{T}_{\text{e}}}\right),}$$(B.13)

Here $${\varphi }_i\left(T_{\text{e}}\right)=\frac{n_i}{n_{m\left(0\right)}^{\ast }{n}_{\text{e}}},$$(B.14)

is the Saha-Boltzmann factor.

If we assume that the stimulated recombination is negative photoionization, the photoionization rate is given by $${\tilde{R}}_{ip}={\gamma }_i-\frac{n_{p_{{\text{L}}_0}}}{n_i}{n}_{\text{e}}{\alpha }_i^{\text{stim}},$$(B.15)

where $${\gamma }_i-\frac{n_{p_{{\text{L}}_0}}{n}_{\text{e}}}{n_i}{\alpha }_i^{\text{stim}}=4\pi {\displaystyle \underset{v_i}{\overset{\infty }{\int }}\text{d}v\left(1-\frac{n_{p_{{\text{L}}_0}}}{n_i}\frac{n_i^{\ast }}{n_{p_{l_0}}^{\ast }}\exp \left(-\frac{hv}{k_{\text{B}}{T}_e}\right)\right)}J\left(v\right).$$(B.16)

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 bound-bound and bound-free transitions are given below. The inverse rates are calculated using $n_i^{*}{C}_{if}=n_f^{*}{C}_{fi}$ (see Hubeny & Mihalas 2015, Eq. 9.52).

Appendix B.2.1 Bound-bound transitions

The bound-bound collisional rates are calculated using the van Regemorter (1962) approximation $$C_{iu}=14.5\cdot n_{\text{e}}{c}_0{T}^{\frac{1}{2}}\left(\frac{I_{{}^H}}{h{v}_{iu}}\right)f_{iu}\frac{h{v}_{iu}}{k_{\text{B}}{T}_{\text{e}}}\exp \left(-\frac{h{v}_{iu}}{k_{\text{B}}{T}_{\text{e}}}\right)\Gamma \left(\frac{h{v}_{iu}}{k_{\text{B}}{T}_{\text{e}}}\right),$$(B.17)

(see also Hubeny & Mihalas 2015, Eq. 9.58), here c₀ = 5.46510. 10⁻¹¹ 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 $$\Gamma \left(x\right)=\max \left(\overline{ɡ}\text{,0}.276\cdot \text{exp}\left(x\right)E_1\left(x\right)\right),$$(B.18)

$\overline{g}$ is given by (here n, n′’ denote principal quantum number, l, l′ the orbital quantum number) $$\widetilde{ɡ}\sim \{\begin{array}{lll}0.7,\hfill & \text{for\hspace{0.17em}transations}\hfill & <n,l>\to <n,l^{\prime }>,\hfill \\ 0.2,\hfill & \text{for\hspace{0.17em}transations}\hfill & <n,l>\to <n^{\prime },l^{\prime }>,\hfill \end{array}$$(B.19)

and E₁ (x) is the exponential integral function.

Appendix B.2.2 Bound-free transitions

The collisional bound-free rates are calculated using an approximate formula (see Hubeny & Mihalas 2015, Eq. 9.60) $$C_{ip}=n_{\text{e}}\frac{1.55\cdot {10}^{13}}{T_{\text{e}}^{\frac{1}{2}}}{\overline{ɡ}}_i\frac{{\alpha }_{ip}}{\frac{h{v}_i}{k_{\text{B}}{T}_{\text{e}}}}\exp \left(-\frac{h{v}_i}{k_{\text{B}}{T}_{\text{e}}}\right),$$(B.20)

where ${\overline{g}}_i$ is equal to $${\overline{ɡ}}_i\sim \{\begin{array}{ll}1\hfill & \text{Charge of\hspace{0.17em}ionization\hspace{0.17em}state\hspace{0.17em}of\hspace{0.17em}lavel\hspace{0.17em}\hspace{0.17em}}i\text{\hspace{0.17em}\hspace{0.17em}is\hspace{0.17em}\hspace{0.17em}}\mathrm{0,}\hfill \\ 2\hfill & \text{Charge of\hspace{0.17em}ionization\hspace{0.17em}state\hspace{0.17em}of\hspace{0.17em}lavel\hspace{0.17em}\hspace{0.17em}}i\text{\hspace{0.17em}\hspace{0.17em}is\hspace{0.17em}\hspace{0.17em}}\mathrm{1,}\hfill \\ 3\hfill & \text{Charge of\hspace{0.17em}ionization\hspace{0.17em}state\hspace{0.17em}of\hspace{0.17em}lavel\hspace{0.17em}\hspace{0.17em}}i\text{\hspace{0.17em}\hspace{0.17em}is\hspace{0.17em}\hspace{0.17em}}\ge 2.\hfill \end{array}$$

Appendix B.2.3 Free-bound transitions

The cross section for the spontaneous free-bound transition from the level i to the level m is $${\alpha }_i^{\text{E,spont}}={\varphi }_i\left(T_{\text{e}}\right)4\pi {\displaystyle \underset{v_i}{\overset{\infty }{\int }}\text{d}v}\frac{{\alpha }_{im}\left(v\right)}{h{v}_i}\frac{2h{v}^3}{c^2}\exp \left(-\frac{hv}{k_{\text{B}}{T}_{\text{e}}}\right)$$(B.21)

where the Saha-Boltzmann factor is defined in (B.14).

Appendix B.3 Emissivities

The free-bound emissivity $${\eta }_{ip}^{\text{bf}}\left(v\right)=n_i^{\ast }{\alpha }_{ip}\left(v\right)\frac{2h{v}^3}{c^2}\exp \left(-\frac{hv}{k_{\text{B}}{T}_{\text{e}}}\right).$$(B.22)

where α_ip(v) is the photoionization cross section for a transition from the level i to the level p.

The free-free emissivity $${\eta }_i^{\text{ff}}{n}_e{N}_j{\alpha }_i^{\text{ff}}\left(v,T_{\text{e}}\right)\frac{2h{v}^3}{c^2}\exp \left(-\frac{hv}{k_{\text{B}}{T}_{\text{e}}}\right),$$(B.23)

N_j represents the population of the ion j, the free-free cross section, see (Kromer 2009, Eq. (3.46)) $$\begin{array}{l}{\alpha }_{i_{{\text{I}}_k}}^{\text{ff}}\left(v,T\right)=\frac{4e^6}{3hc}{\left(\frac{2\pi }{3m_{\text{e}}^2{k}_{\text{B}}}\right)}^{1/2}\frac{{\left(k-1\right)}^2}{T^{1/2}}\frac{{ɡ}_{\text{ff}}\left(j\right)}{v^3}\hfill \\ =3.69255\times {10}^8\cdot {\left(k-1\right)}^2{ɡ}_{\text{ff}}\left(j\right)\cdot T^{-1/2}\times v^{-3}.\hfill \end{array}$$(B.24)

If we put the Eq. (B.23) into Eq. (35) we can find easily an analytic formula $$v=\frac{k_{\text{B}}T}{h}\log \left(z\right),$$(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 $$\frac{n_u}{n_l}=\frac{{ɡ}_u}{{ɡ}_l}\exp \left(-\frac{{\epsilon }_u-{\epsilon }_l}{k_{\text{B}}{T}_{\text{e}}}\right),$$(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 $$\frac{N_j}{N_p{n}_{\text{e}}}=\frac{U_j}{U_p{T}_{\text{e}}^{3/2}}\exp \left(\frac{{\epsilon }_{p_{\text{L}}{}_{{}_0}}-{\epsilon }_{j_{{\text{L}}_0}}}{k_{\text{B}}{T}_{\text{e}}}\right),$$(C.2)

where C is the Saha constant $$C=\frac{1}{2}{\left(\frac{2\pi {ħ}^2}{m_{\text{e}}{k}_{\text{B}}}\right)}^{\frac{3}{2}},$$(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 $$U_i\left(T_{\text{e}}\right)={\displaystyle \sum _{j=1}^{{\text{N}}_{\text{L}}}{ɡ}_j\exp \left(-\frac{{\epsilon }_j}{k_{\text{B}}{T}_{\text{e}}}\right)},$$(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.