# Spectral line fluorescence in moving envelopes of stars: I. Theoretical background, numerical code, and proof of concept

Claude Bertout

► To cite this version:

Claude Bertout. Spectral line fluorescence in moving envelopes of stars: I. Theoretical background, numerical code, and proof of concept. Astronomy & Astrophysics - A&A, 2024, 689, 10.1051/0004-6361/202449817. insu-04851542

HAL Id: insu-04851542

https://insu.hal.science/insu-04851542v1

Submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spectral line fluorescence in moving envelopes of stars

## I. Theoretical background, numerical code, and proof of concept

Claude Bertout*

Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France
e-mail: cb.work@snoqualm.org

Received 29 February 2024 / Accepted 28 May 2024

## ABSTRACT

**Context.** Fluorescence of the optical FeI lines is observed in active T Tauri stars, and is considered a defining characteristic of this class of young stellar objects.

**Aims.** The formation of optical fluorescent lines in moving media has not yet been studied in detail. This work represents a first step in investigating the fluorescence process in different types of macroscopic velocity fields: (a) accelerated outflows, (b) accelerated infalls, and (c) non-monotonic velocity fields (such as an accelerating outflow followed by a deceleration region or an accretion shock front). We aim to develop a general computer code for studying the fluorescent emission in any 2D macroscopic velocity field. As an illustration, we investigate FeI T Tauri-like fluorescent emission in these moving stellar envelopes.

**Methods.** We solved the radiative transfer equations for the lines involved in the fluorescent process, assuming spherical symmetry and a simplified atomic model. We used the framework of the generalized Sobolev theory to compute the interacting, nonlocal source functions. The emergent line fluxes were then integrated exactly.

**Results.** Because of Doppler shifts in the moving gaseous envelope, photons of the three lines involved in T Tauri star FeI fluorescence (CaII H, FeIλ3969, and H<sub>ε</sub>) interact with each other in a complex way, such that fluorescent amplification of the line flux occurs not only for FeIλ3969, but also for the other two lines, in all velocity fields that we investigated. With the assumption of local thermodynamic equilibrium populations, the line source functions of moderately optically thick lines are not strongly affected by line interactions, while they are depressed in the inner envelope for optically thick lines because of stellar photon absorption in the interaction regions. Fluorescent amplification takes place mainly in the observer's reference frame during the flux integration. We define a measure of fluorescence based on the line equivalent widths and perform a parameter study for an accretion flow with a broad range of envelope temperatures and densities while including approximate collisional de-excitation rates in the source function computations. Significant fluorescence occurs over the entire temperature range of the investigated flow, but only in the higher density range, suggesting that relatively high mass accretion rates are needed to trigger the fluorescence process.

**Conclusions.** Further comparison with observations will require solving the rate equations for the atomic populations involved, along with the radiation field computed with the method presented here. The main product of this research is the open-source computer code SLIM2 (Spectral Line Interactions in Moving Media), written in Python/NumPy, which numerically solves the fluorescence problem for arbitrary 2D velocities.

**Key words.** line: formation – radiative transfer – methods: numerical – stars: variables: T Tauri, Herbig Ae/Be

## 1. Introduction

The physical process of fluorescence, also known as photoexcitation by accidental resonance, is associated in stellar astronomy with the name Ira S. Bowen, who showed that the spectra and composition of gaseous nebulae could be explained by fluorescent emission of the ions OIII and NIII excited by neighboring HeII transitions (Bowen 1934, 1935). This solved a long-standing mystery that had even led Huggins (1898) to postulate the existence of a chemical element he called “nebulium”. Fluorescent emission was also successful in explaining features of the line spectra in long-period variables (Thackeray 1937; Willson 1972) as well as in symbiotic stars (e.g., Kastner 1991; Pereira et al. 1999; Eriksson et al. 2004; Selvelli et al. 2007; Hyung et al. 2018).

A similar mechanism was proposed by Herbig (1945) to account for the anomalous Fe I emission observed in active T Tauri stars (TTSs). Specifically, the transitions λ4063 and λ4132 of Fe I appear in strong emission in active TTSs, while other lines of the multiplet $a^3F - y^3F^0$ are in very weak emission or are not detected (Joy 1945). The two anomalous lines originate from the same $y^3F^0$ upper level together with the λ3969 line. Herbig conjectured that “the explanation lies in the near coincidence of the λ3969 line of Fe I with the red component of the H line of Ca II. Absorption of the proper frequencies by the Fe I atoms would ensure a supply of atoms to $y^3F^0$, while the other sub-levels of $J = 2$ and 4 would lack such an active source.” The λ3969.3 line is flanked by Ca II H at λ3968.5 on the blue side and by the H I $H_\epsilon$ line at λ3970.1 on the red side. Both of these lines are in strong emission in TTSs, and given the complexity of motions in the surroundings of these stars, it is conceivable that either one could be responsible for pumping atoms into the upper level of Fe I λ3969. Herbig did not consider the possible role of $H_\epsilon$ in the Fe I fluorescence, presumably

<footnote>* Retired.</footnote>
because it was believed in 1945 that TTS activity was a scaled-up version of solar-type, chromospheric activity accompanied by a strong stellar wind. However, Willson (1974, 1975) concluded from a detailed investigation of the statistical equilibrium of the lines involved that both H and H<sub>e</sub> were contributing to the fluorescence mechanism in the active TTS RW Aur, and he suggested that rotational velocities were greater than ejection velocities in the envelope of that star. While we know now that TTSs are not fast rotators, and that TTS activity is fueled mainly by accretion onto the star of disk material, it is certainly plausible that either H or H<sub>e</sub>, or both, are at work in the fluorescence process. And given the wealth of mid- and high-resolution spectroscopic data that has become available with the advent of such instruments as ESO's ESPRESSO, UVES, and X-Shooter (see, e.g., Armeni et al. 2023), it is obvious that a renewed quantitative study of fluorescence in young stars is overdue. Our aim is to present the theoretical background as well as a numerical code that enables such a study. The problem is relatively complex, and we will therefore solve it in two steps.

The first one, presented in the present paper, focuses on the formation of the λ3969 FeI line and its interactions with the CaII H and HI H<sub>e</sub> lines, using the equivalent two-level atom approach for each of the lines and the assumption of local thermodynamic equilibrium (LTE) for computing the populations of the atomic levels involved. Our aim is to create a code for solving the radiative interactions between lines that occur in the moving media that we consider and to gain some understanding of the physical conditions of the envelope that lead to amplification of the λ3969 line flux.

This first investigation will not make a direct comparison between models and observations possible, since we will not be solving the statistical equilibrium equations for the level populations of FeI self-consistently with the radiation field. In a forthcoming extension of this work, we will therefore couple the code developed here with a computation of the non-LTE rate equations for a multilevel Fe atom model, which will then allow us to compare the anomalous FeIλ4063 and λ4132 lines, as well as the nonfluorescent lines of the same multiplet, to observations.

The physical processes taking place in the most active T Tauri envelopes involve several components: circumstellar disk, magnetospheric accretion columns between disk and star, bipolar disk wind, and stellar wind. Modeling the line formation process in this dynamic environment is challenging, but considerable advances have been achieved over the past 20 years (e.g., Tambovtseva et al. 2001, 2014; Kurosawa & Romanova 2012). All of this work has been done in the framework of the generalized Sobolev approximation, which makes this problem tractable.

In a moment of remarkable insight, Sobolev (1947, English translation in 1960) realized that the Doppler shift simplifies the formation of spectral lines in moving media such as stellar envelopes and winds. The basic approximation made in what is now called the Sobolev theory is that, because of the macroscopic velocity gradient in the moving gas, the source function and optical depth of a given transition are constant over the narrow region where emitted photons are not Doppler-shifted with respect to atoms in the rest of the envelope, and equal to zero everywhere else. This is true in particular in outwardly accelerating, supersonic winds, which represent the simplest case for applying the Sobolev theory; the line formation problem then reduces to finding regions of constant radial velocity with respect to the observer, each of which corresponds to a single frequency in the emergent line profile.

In more quantitative terms, Sobolev assumed in his original derivation that the local absorption and emission profiles were constant over the width $\Delta v_{ik} = 2(u/c)v_{ik}$ and equal to zero everywhere else. Here, $v_{ik}$ is the frequency of transition *ik*, *u* is the mean thermal velocity of atoms, and *c* is the velocity of light. For a hydrogen line at $10^4$ K, the emission and absorption profile widths are on the order of 10 km s<sup>-1</sup>, so $\Delta v_{ik}/v_{ik}$ is on the order of $10^{-4}$, or 0.5 Å for a line wavelength of 5000 Å. This is to be compared to the overall line width, which will extend, for example, over 16 Å in a 500 km s<sup>-1</sup> outflow. In other words, whenever the Sobolev approximation holds, absorption and emission line profiles can be seen as delta-like functions, which is consistent with considering that the emitting regions at a given line frequency are radial velocity surfaces rather than extended volumes.

The Sobolev theory has been used extensively for modeling line formation in envelopes of stars and outflows in a variety of objects, ranging from early-type stars (e.g., Castor 1970; Lamers et al. 1987) to young stellar objects (e.g., Dmitriev et al. 2019) and from interstellar clouds (Goldreich & Kwan 1974) to bipolar outflows (e.g., Cabrit & Bertout 1986, 1990, 1992). A review of the rich literature applying Sobolev's theory in various contexts is beyond the scope of this paper, but it is fair to say that this approach has led to huge progress in our understanding of line formation in many different astrophysical situations, as Sobolev himself had foreseen in his original treatise.

Another reason for this success – besides the simplification in line formation that we have alluded to – stems from the fact that the Sobolev method was shown, by extensive comparisons with much more computationally challenging numerical solutions of the exact line transfer equations, to be sufficiently accurate for a large class of problems, especially when the approximation is used only for computing the line source function and the emergent line flux is integrated exactly (see, e.g., Bastian et al. 1980; Hamann 1981).

Originally used for the study of accelerating outflows, the Sobolev method has been generalized and extended over the years to more involved velocity fields, including decelerated outflows, accelerated inflows, rotation, and various combinations thereof. The generalization of the original Sobolev theory for treating velocity laws that lead to complex radial velocity surfaces was worked out by Rybicki & Hummer (1978) in a seminal paper. The same theoretical framework has also been useful for studying the interactions between lines that occur in some astrophysical situations, for example close doublets (cf. Olson 1982) and masers (cf. Pavlakis & Kylafis 1996). The latter authors considered the case of rectangular local emission and absorption profiles, in which case the line source function computation of the locally interacting lines becomes analytically tractable, albeit tedious.

In static media, spectral lines interact when their absorption profiles overlap locally. If macroscopic velocity gradients are present in the line-emitting region, however, it can also happen that the line of interest is Doppler-shifted onto another one that is too far in wavelength for a local interaction to occur; following Pavlakis & Kylafis (1996), we call this situation a nonlocal interaction.

Quantitatively, two lines, A and B, with rest wavelengths $\lambda_A$ and $\lambda_B$ and absorption profile widths $\Delta\lambda_A$ and $\Delta\lambda_B$ interact locally when $(\lambda_B - \lambda_A)/\lambda_A \lesssim \Delta\lambda_A + \Delta\lambda_B$. They interact non-locally whenever line-A photons emitted by matter moving at velocity $v_A$ cross
a distant region where the radial velocity, $v_r^B$, of the B-line-emitting matter is such that $(\lambda_B - \lambda_A)/\lambda_A = (v_r^B - v_r^A)/c^1$. In the example discussed earlier, the line at 5000 Å will interact to some extent with neighboring lines located within 8 Å of the transition rest wavelength.

In the present work, we calculate the line source functions and emergent line profiles of radiatively overlapping lines in moving media. While the generalized Sobolev theory à la Hummer and Rybicki remains the conceptual framework underlying the computation, we make use of the Sobolev approximation for computing the source function and, following Bertout (1984), perform an exact integration of the emerging flux, assuming Gaussian shapes for the local absorption profiles.

We make several simplifying assumptions in the subsequent discussion. The moving medium under consideration is a spherically symmetric, extended stellar envelope with maximum radius $R_{\text{env}}$ surrounding a star of mass $M_*$, radius $R_c$, and effective temperature $T_{\text{eff}}$. The envelope velocity field is, at first, assumed to be a monotonically accelerating outflow from velocity $V_c$ at the stellar core to $V_{\text{max}}$ at $R_{\text{env}}$. We then consider velocity fields that lead to more complex distant radiative coupling. We neglect continuous opacity in the envelope, which has temperature $T_{\text{env}}(r)$. The gas number density, $n_t$, is computed from the continuity equation. The emitting gas has a solar composition, and the atomic level populations and free electron number density, $n_e(r)$, follow from the assumption of LTE. In the framework of equivalent two-level atoms, $n_l$ and $n_u$, the occupation numbers of both line levels in Eq. (8), are then given by the Saha-Boltzmann equation evaluated at the local temperature and electron density.

Collisional de-excitation of atoms is an important radiation source for the lines that we investigate, so we included approximate rates, $\epsilon_k$, in the source function computation for line k. However, since most of the previous theoretical work done on line formation in moving media neglects the collisional term and since we want to be able to compare our results to previous work, we will assume $\epsilon = 0$ in the theoretical sections of this work and reserve the use of collisional terms to Sects. 3 and 4.

In the next section, we study the fluorescence process in a monotonically accelerating outflow. The Sobolev approximation for the source functions of interacting lines is discussed in Sect. 2.1, while Sect. 2.2 extends Bertout (1984) exact flux integration method to the case at hand. In Sect. 3 we generalize the analysis to the more complex radial velocity surfaces found in accreting envelopes (Sect. 3.1) and in non-monotonic velocity fields (Sect. 3.2). A discussion of the physical conditions that lead to fluorescence follows in Sect. 4. The FeI fluorescence observed in TTSs is used as an example in the various illustrations of the paper.

Finally, we briefly present in Appendix A the open-source computer code (available on GitHub under CC 4.0 license) written in Python/NumPy that solves the nonlocal radiative coupling between interacting spectral lines formed in spherical, moving envelopes of stars.

## 2. Formation of interacting lines in a monotonic, accelerated outflow

We considered an outwardly accelerating stellar outflow with properties loosely inspired by the Alfvén-wave-driven stellar wind originally proposed for TTSs by Hartmann et al. (1982). At first, we used the data they give in their table for the velocity of Model 5, interpolating between the few provided radial points, but this resulted in discontinuous derivatives of the velocity at the grid points that in turn led to spurious variations in the source function of the optically thick hydrogen line. We therefore looked for a least-square fit to the data, assuming a parameterized velocity of the form

$$v(r) = v_{\infty} \left[ 1 - (1 - v_{\text{min}}) \frac{r}{R_c} \right]^{\alpha} \quad (1)$$

The resulting parameters are given in Table 1 together with the other computation parameters, and Fig. 1 shows the velocity law together with the adopted gas density law, which follows from the continuity equation

$$n_t(r) = n_t(R_c) \frac{v(R_c)}{v(r)} \left(\frac{R_c}{r}\right)^2. \quad (2)$$

Our choice of computation parameters for the envelope properties is informed by spectroscopic data of the active TTS system SCrA (e.g., Bertout et al. 1982), which display the FeI fluorescence phenomenon discussed here. In these spectra, the intensities of the CaII H + H$_e$ blend and the CaII K line, both in strong emission, are of comparable strength. This is also the case, more generally, for many TTSs (Valenti et al. 1993). For this to occur (within the limitations of the LTE assumption), the temperature of the line formation region cannot be higher than $6-7 \times 10^3$ K. At such low temperatures, we must have a relatively high gas density in order to get significant emission line flux. In the range of values chosen here, the three lines are partially optically thick at the bottom of the envelope (Fig. 2). Clearly, near-photospheric conditions are required for producing moderately strong emission lines when using the assumptions of LTE populations and isothermal gas, and when neglecting collisions. As a consequence, the mass-flow rate determined using these assumptions and the additional condition of spherical symmetry is much higher than usually assumed for TTSs. The situation improves when collisions are taken into account, as we will see in Sect. 4. But our primary goal in this section, as already mentioned, is to present the method for solving the fluorescent line formation in the case $\epsilon = 0$, not to compare the results with observations.

<footnote n="1">¹ For the sake of simplicity, we neglected the line widths in this condition.</footnote>
<table><caption>**Table 1.** Computation parameters for the outflow of Sect. 2.</caption><thead><tr><th>R<sub>c</sub></th><th>R<sub>env</sub></th><th>M<sub>*</sub></th><th>T<sub>eff</sub></th><th>T<sub>env</sub></th><th>n<sub>t</sub>(R<sub>c</sub>)</th><th>v<sub>∞</sub></th><th>v<sub>min</sub></th><th>α</th></tr></thead><tbody><tr><td>2 R<sub>☉</sub></td><td>10 R<sub>c</sub></td><td>1 M<sub>☉</sub></td><td>4.5 × 10<sup>3</sup> K</td><td>6 × 10<sup>3</sup> K</td><td>2 × 10<sup>16</sup> cm<sup>-3</sup></td><td>250 km s<sup>-1</sup></td><td>0.02</td><td>0.95</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Velocity (Eq. (1)) and density law (Eq. (2)) of the stellar outflow.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Line absorption coefficient for the three lines considered.</figcaption></figure>

## 2.1. Sobolev source functions

We consider three spectral transitions: *B*, *C*, and *R* with rest wavelengths $\lambda_B < \lambda_C < \lambda_R$ in a wavelength interval

$$ (\lambda_R - \lambda_B) = \lambda_R \cdot v(R_{\text{env}})/c. \qquad (3) $$

For convenience, the lines are denoted *B* for blue, *C* for center, and *R* for red; these lines may stem from different chemical elements and have respective thermal widths $\Delta v_B$, $\Delta v_C$, and $\Delta v_R$.

The generalized Sobolev source functions of interacting lines can be computed by modifying the formalism developed by Rybicki & Hummer (1978) for the study of single lines emitted in velocity fields that lead to nonlocal radiative coupling. We assume in the following that readers are familiar with their work.

### 2.1.1. *B* line

We first note that in any spherical, monotonically accelerating large-scale flow, photons emitted by the bluest line (the *B* line as defined above) appear redshifted to atoms responsible for the redder *C* and *R* transitions and thus escape the medium freely. The line source of *B* is therefore unaffected by *C* and/or *R*, and its mathematical expression is the same as the one derived for a single line by Rybicki and Hummer. For an equivalent two-level atom, we then have

$$ S_B(r) = \frac{(1 - \varepsilon_B)\beta_B^c I_B^c + \varepsilon_B B_B}{\varepsilon_B + (1 - \varepsilon_B)\beta_B}, \qquad (4) $$

where $\beta_B$ is the local escape probability, $\beta_B^c$ the probability that an escaping photon will strike the stellar core, $I_B^c$ is the background continuous radiation emitted by the stellar core at line frequency, $\varepsilon_B$ the ratio of collisional to total de-excitation rates from the upper level, and $B_B$ is the Planck function at rest line frequency and local electron temperature. All these quantities are defined locally at radius $r$. With $\mu = \cos \theta$ defining the direction cosine of a ray with respect to the outward radial direction, the local escape probability is given by

$$ \beta_B = \frac{1}{2} \int_{-1}^{+1} \frac{1 - e^{-\tau_B}}{\tau_B} d\mu, \qquad (5) $$
while the probability for an escaping photon to reach the stellar core is

$$ \beta_B^c = \frac{1}{2} \int_{\mu_c}^{+1} \frac{1 - e^{-\tau_B}}{\tau_B} d\mu, \qquad (6) $$

where $\mu_c$ is the cosine of the half-angle subtended by the stellar core. The quantity $\tau_B$ in the last two equations is the B line-center optical depth in direction $\mu$. Supposing that the transition occurs between electronic levels $l$ and $u$, $\tau_B$ is given by (see, e.g., Hubeny & Mihalas 2014, Chapter 19, p.746ff)

$$ \tau_B(r, \mu) = \chi_{\text{lu}}(r)/Q(r, \mu), \qquad (7) $$

with the absorption coefficient given by

$$ \chi_{\text{lu}}(r) = (\pi e^2 / mc) f_{\text{lu}}[n_1(r) - (g_1/g_u)n_u(r)] / \Delta v_D \qquad (8) $$

and

$$ Q(r, \mu) = |dv_s/ds|, \qquad (9) $$

where $v_s$ is the local fluid velocity projected in direction $s$ with cosine angle $\mu$, $\Delta v_D = v_0 v_{th} / c$ the local thermal line width, and where the other symbols have their usual meaning. The absorption coefficients of the three considered lines, as expressed by Eq. (8), are shown in Fig. 2 for the parameters of Table 1.

### 2.1.2. C line

We next consider photons of transition C emitted in the outer parts of the envelope. They will interact with B photons emerging from deeper parts of the envelope, where the emitting matter is moving slower, whenever the line-of-sight velocity difference is equal to the rest velocity difference between the two lines. In other words, interaction between B and C photons on any given line of sight will take place whenever

$$ v_{\text{rad}}(r_B) - v_{\text{rad}}(r_C) = \Delta\lambda_{BC}/c, \qquad (10) $$

where $\Delta\lambda_{BC}$ is the wavelength difference between lines B and C. Solving Eq. (10) over the entire space defines a locus of distant points $r'$ where emitted B photons can interact with C photons emitted at $r$. This locus is called a common-point (CP) surface by Olson (1982), using the nomenclature introduced by Rybicki & Hummer (1978) to describe the line self-coupling that occurs in those velocity laws that lead to distant radiative coupling. We use the same denomination and denote the CP surface described above by $r'_B$ in the following equations.

Figure 3 shows some CP surfaces for the lines of interest in the stellar wind with properties shown in Fig. 1 and parameters given in Table 1. In that figure, the CP surfaces corresponding to a given radial point (identified by a colored dot in the figure) are represented by dashed and dotted lines sharing the same color as the radial point they are associated with. As mentioned earlier, the bluest line $\lambda3968$ is formed locally and has no associated CP surfaces. The dashed lines show the CP surfaces caused by interactions of the two redder lines C and R with their closest neighboring line (i.e., the interaction of $\lambda3969$ with $\lambda3968$ photons in the left panel, and interaction of $\lambda3970$ with $\lambda3969$ photons in the right panel of Fig. 3). Finally, the dotted lines denote the CP surface caused by interaction with the farthest spectral line (i.e., by the interaction of $\lambda3970$ with $\lambda3968$ photons; right panel). In other words, the left panel of Fig. 3, with only one set of CP surfaces, represents the situation for the C line discussed here, while the right panel represents the situation for the R line discussed in the next section$^2$.

As mentioned by Olson (1982) and also seen in our figures, the topology of CP surfaces is quite complex and varies with radial distance from the stellar core. Close to the star, the CP surface is an arc of circle surrounding the reference point and extending from the abscissa to the angle at which the line of sight intercepts the core. At larger distances it becomes a closed surface entirely surrounding the reference point, which then opens farther out as the envelope limit is reached. We see below that these topological changes are reflected in the source functions of optically thick lines.

The source function for transition C takes into account the presence of resonant CP surfaces coupling it to the B transition. It is written as (see also Olson 1982)

$$ S_C(r) = \frac{(1 - \varepsilon_C)\beta_C^c I_C^c + \varepsilon_C B_C + (1 - \varepsilon_C)S_{nl}^C(r)}{\varepsilon_C + (1 - \varepsilon_C)\beta_C}, \qquad (11) $$

where the nonlocal contribution $S_{nl}^C$ to the source is

$$ S_{C}^{nl} = \frac{1}{2} \int_{-1}^{+1} \frac{1 - e^{-\tau_C}}{\tau_C} S(r'_{B})(1 - e^{-\tau(r'_{B})}) d\mu \qquad (12) $$

<footnote>$^2$ Please note that we keep the same graphical conventions for the representation of CP surfaces in subsequent figures.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Selected CP surfaces for the spherically symmetric stellar wind, with the properties given in the text. Each surface is associated with the radial point of the same color. Note the variety of geometries for the surfaces, which give rise to discontinuities in the source function.</figcaption></figure>

and the probability that a photon strikes the core is modified to take into account the intervening B layer:

$$ \beta_C^c = \frac{1}{2} \int_{\mu_c}^{+1} \frac{1 - e^{-\tau_C}}{\tau_C} e^{-\tau(r_B')} d\mu. \qquad (13) $$

### 2.1.3. R line

Similarly, the third and most redward transition R will be affected by photons of both lines C and B, and we must locate the two CP surfaces $r''_B$ and $r'''_C$ where B and C have the same line-of-sight velocity as the gas at the considered R emission point r, using Eq. (14):

$$ \left\{ \begin{array}{l} v_{\text{rad}}(r_R) - v_{\text{rad}}(r_B) = \Delta\lambda_{BR}/c \\ v_{\text{rad}}(r_R) - v_{\text{rad}}(r_C) = \Delta\lambda_{RC}/c \end{array} \right. \qquad (14) $$

The generalization of Eqs. (10)–(13) to line R is obvious and leads to the following expression for $S_R(r)$:

$$ S_R(r) = \frac{(1 - \varepsilon_R)\beta_R^c I_R^c + \varepsilon_R B_R + (1 - \varepsilon_R)S_{nl}^R(r)}{\varepsilon_R + (1 - \varepsilon_R)\beta_R}, \qquad (15) $$

with the nonlocal integral $S_R^{nl}$ given by

$$ S_R^{nl} = \frac{1}{2} \int_{-1}^{+1} \frac{1 - e^{-\tau_R}}{\tau_R} \left[ S(r_B'') (1 - e^{-\tau(r_B')}) + S(r_C'') (1 - e^{-\tau(r_C'')}) e^{-\tau(r_B'')} \right] d\mu \qquad (16) $$

and

$$ \beta_R^c = \frac{1}{2} \int_{\mu_c}^{+1} \frac{1 - e^{-\tau_R}}{\tau_R} e^{-\tau(r_B'') - \tau(r_C'')} d\mu. \qquad (17) $$

The geometry of CP surfaces, where the various quantities entering Eqs. (12)–(17) are computed, depends only on the gas velocity law and the transitions wavelengths. It is computed numerically by integrating over (a) the entire space when computing
<figure><img src="image_5.png" /><figcaption>Fig. 4. Source functions for the three lines when they are computed independently.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Source functions for the three lines when radiative interactions are taken into account.</figcaption></figure>

the nonlocal parts of the source functions and (b) the angle subtended by the stellar core when computing the stellar continuum contribution to the source. The computation of $S_B$ is straightforward algebra in this approximation since all quantities are local, while computing $S_C$ and $S_R$ requires an iterative process to solve for the nonlocal parts $S^{nl}$ consistently.

The above discussion suggests constructing an interaction matrix $I$ to keep track of each possible line interaction, a concept that proved useful in coding the problem at hand. The matrix elements are 0 for a line that can be considered single or 1 when an interaction is present. For the current case of three lines in an accelerating outflow, with line ordering $B, C, R$, the interaction matrix is

$$I = \begin{pmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 1 & 1 & 0 \end{pmatrix}, \qquad (18)$$

meaning that the $B$ line (first row) has no interaction with itself, $C$, or $R$, while the $C$ line (second row) interacts with $B$, and the $R$ line (third row) interacts with both $B$ and $C$.

Figure 4 displays the computed source functions for the three lines under consideration when they are treated as independent, that is to say, when local equations similar to Eq. (4) are used for computing all of them. Figure 5 shows the source functions computed when considering the radiative interactions (i.e., when using Eq. (11) for the $C$ line and Eq. (15) for the $R$ line). While the source function of the $B$ line is evidently the same in both cases, discontinuities caused by the different topologies of CP surfaces are seen; but since the lines are only moderately optically thick, the induced changes remain small and do not affect the resulting source functions very much. The assumption of LTE atomic populations, which results in constant radial optical depths as the source functions are iterated, also plays a role here. Allowing for an iterative estimate of populations is thus likely to result in very different source functions. While the two-level atom assumption is clearly a limitation of the current code, it is a first step that allows us to understand the radiative couplings of the lines occurring in both the comoving reference frame and the observer's frame, and the next section deals with the latter.

## 2.2. Exact line profile integration

We refer to Bertout (1984) for a detailed discussion of the line profile integration method for single lines and focus here on the specific modifications that are needed to introduce the possibility of line interaction in the flux integration. The emergent flux $F_{x,l}$ of line $l$ at frequency displacement $x = (\nu - \nu_{0,l})/\Delta\nu_l$ is

$$F_{x,l} = 2\pi \int_{0}^{R_{\max}} I_l(p, x, z = \infty) p \, dp, \qquad (19)$$

where impact parameter $p$ and abscissa $z$ are related to the usual $r$ and $\theta$ coordinates by the relationships $r^2 = p^2 + z^2$ and $z = r \cos \theta$. The quantity $I_l(p, x, z = \infty)$ is the emergent intensity at frequency $x$ on the ray $p$, and $\Delta\nu_l$ is the width of the local absorption profile.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Regions and frequencies contributing to the line profiles of the B, C, and R lines (respectively for each column) in the wind model of Sect. 2. Shown here is a meridional cut in the half spherical stellar envelope of radius $R_{max} = 10 R_c$ surrounding a stellar core of radius $R_c$. Displayed are the constant radial velocity surfaces (CRVSs) seen by an outside observer located at infinity on the right side of the figure for the three lines considered (left-right diagonal), as well as the displaced CRVSs of neighboring lines. The scale of the color bars is in km s⁻¹, and blue (red) indicates positive (negative) line frequencies.</figcaption></figure>

The Doppler equation at $(p, z)$ is

$$x(p, z) = \frac{v(r) z}{\xi_l(r) r}, \qquad (20)$$

where $\xi_l(r)$ is the width of the line absorption profile in the same units as $v(r)$. The local absorption profile $\xi_l(r)$ results from both thermal broadening and an additional turbulent widening $v_{turb}$, chosen to be equal to 1 km s⁻¹ in the actual computations. The only constraint on the choice of the overall absorption profile width in velocity space is that its integral over velocity (or frequency) be the same for the three considered lines (see below). To provide a frequency rest frame common to these lines, we thus chose the width $\xi$ to be equal to or larger than the widest local line width expressed in velocity units. In other words,

$$\xi = (\alpha + 1)\max[\xi_l(r)], \qquad (21)$$

where $\alpha$ is typically in the range [0, 1].

In an accelerated outflow, there is a unique correspondence between line frequencies and positions on the line of sight. This is shown in the top left panel of Fig. 6, which displays the envelope’s gas radial velocities seen by an outside observer located at $z = \infty$. The color code of radial velocities is given in the bar under the figure. Impact parameter $p = 0$ contains the full range of velocities (or, equivalently, frequency displacements) present in the three emergent line profiles, ranging from $-V_{max}$ at $z = -R_{max}$ to $+V_{max}$ at $z = R_{max}$. We set up a uniform grid over this velocity range with the understanding that the velocity step $dv$, which defines the numerical resolution of the line profile computation, is small enough to allow the narrowest line absorption profile, which has width $\xi_l$, to be resolved with sufficient accuracy. If we consider an impact parameter $p_j$, it will contain velocities in the range

$$v_{\text{rad}}^{\text{min}} = V_{\text{max}} \frac{z_{\text{min}}}{R_{\text{max}}} \le v_{\text{rad}} \le V_{\text{max}} \frac{z_{\text{max}}}{R_{\text{max}}} = v_{\text{rad}}^{\text{max}} \qquad (22)$$
with

$$z_{\min} = \begin{cases} \sqrt{R_c^2 - p_j^2}, & p_j \le R_c \\ -\sqrt{R_{\max}^2 - p_j^2}, & p_j > R_c \end{cases} \quad \text{and} \quad z_{\max} = \sqrt{R_{\max}^2 - p_j^2}.$$

### 2.2.1. Line-C flux

In the observer's frame, there is no preferred flow direction, so all transitions are treated in the same way. We now consider the formation of line C in some detail, since it is our primary candidate for fluorescence, with the understanding that the two other lines are formed in a similar way. For the following discussion, the reader is invited to refer to the central column of Fig. 6, which illustrates the origins of line C.

When they reach the observer at $z = \infty$, line-C photons with apparent frequency $\tilde{x}$ traveling in the p direction originate in atomic transitions taking place in three different regions:

(a) Transition C photons are emitted from atoms in the z-range corresponding to frequencies around $\tilde{x}$ in the local line emission profile $\Delta\nu_{D_c}$. This is the situation illustrated by the central figure of Fig. 6, which displays the constant radial velocity surfaces (CRVSs) of C-emitting atoms seen by the observer. As explained above, each of the CRVSs corresponds to a frequency $\tilde{x}$ in the overall line profile.

(b) A second contribution to the C-line intensity on p stems from transition B photons that are emitted by atoms moving slower by the velocity difference between the two lines $c\Delta\nu_{BC}/\nu_C \pm \Delta\nu_{D_c}$ and are therefore located deeper in the envelope (as seen by the observer) than atoms emitting transition C frequency $\tilde{x}$. The top panel of the central column in Fig. 6 displays the frequency range of this contribution or, equivalently, the envelope parts from which it originates.

(c) Finally, there is a third contribution from transition R photons emitted by atoms that are moving faster by $c\Delta\nu_{CR}/\nu_C \pm \Delta\nu_{D_c}$ and are therefore located in front of the atoms that emit transition C at frequency $\tilde{x}$. The bottom panel of the central column in Fig. 6 displays the space-frequency distribution of this contribution.

Similar considerations hold for the other two lines and the resulting space-frequency contributions to these lines are illustrated by the first and third columns of Fig. 6.

We next considered the contributions (a) to (c) above to the optical depth and emergent C-line intensity at $\tilde{x}$ in a more quantitative way.

**Local contribution to the C-line intensity.** The z-range corresponding to the local absorption profile frequency range is $z(\tilde{x} - \xi/2) < \tilde{z} < z(\tilde{x} + \xi/2)$, and the optical depth along p is given by

$$\tau_C(p, z, \tilde{x}) = \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}-\xi/2)}^{z} \frac{\kappa_C(r')}{\Delta\nu_{D_c}(r')} \varphi(x') dz', \qquad (23)$$

where $\kappa_C$ is the line center opacity and $\varphi$ the line absorption profile (see also Eqs. (7)–(8)). Let us denote with $\tau_C^*$ the maximum value$^3$ of $\tau_C$ on the ray:

$$\tau_C^* = \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}-\xi/2)}^{z(\tilde{x}+\xi/2)} \frac{\kappa_C(r)}{\Delta\nu_{D_c}(r)} \varphi(x) dz. \qquad (24)$$

Moving in the $-z$ direction along impact parameter $p$, and remembering that the line frequency displacement decreases with $z$ until the line center is reached at $z = 0$, the line optical depth then varies as

$$\tau_C(p, z, \tilde{x}) = \begin{cases} 0, & z > z(\tilde{x} + \xi/2) \\ \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}+\xi/2)}^{z} \frac{\kappa_C(r')}{\Delta\nu_{D_c}(r')} \varphi(x') dz', & z(\tilde{x} + \xi/2) \ge z > z(\tilde{x} - \xi/2) \\ \tau_C^*, & z < z(\tilde{x} - \xi/2). \end{cases} \qquad (25)$$

In other words, the optical depth varies from 0 to $\tau_C^*$ over a short distance and then remains constant over the ray. The line intensity at $\tilde{x}$ is

$$I_C^c(p, z, \tilde{x}) = \int_0^{\tau(\tilde{x})} S_C(p', z', \tilde{x}') \exp[-\tau_C(p', z', \tilde{x}')] d\tau(p', z', \tilde{x}'), \qquad (26)$$

where the c subscript indicates that this intensity part is caused by C-line photons.

Following Bertout (1984), we rewrite Eq. (23) as

$$\tau_C(p, z, x_j) = \sum_{ij} \Delta\tau_{ij} \qquad (27)$$

<footnote>$^3$ We use a similar notation for the maximum values of the other lines' optical depths.</footnote>
with

$$ \Delta\tau_{ij} = \frac{1}{\sqrt{\pi}} \left\langle \frac{\kappa_{l,i}}{\Delta\nu_{D_C}} \left(\frac{dx}{dz}\right)_i^{-1} \right\rangle \int_{x_{j-1/2}}^{x_{j+1/2}} \varphi(x')dx'. \quad (28) $$

In Eq. (27), the summation over *i* runs over the points $z_i$ of the *z* grid defined on impact parameter *p*, while the summation over *j* runs over frequency displacements at which a contribution is expected at $z_i$ (i.e., for $x_j$ in the range $x_i - \xi/2 < x_j < x_i + \xi/2$, where $x_i = x(z_i)$). The averaged quantity in Eq. (28) is evaluated over the frequency step $\delta x$, and the integral is equal to $[\text{erf}(x_{j+1/2}) - \text{erf}(x_{j-1/2})]/2$ if the local absorption profile is Gaussian.

Equation (26) can now be integrated to obtain the emergent intensity of line *C* if it were a single line without interactions with its neighbors:

$$ I_C^c(p, z = \infty, \tilde{x}) = I_C^*(\tilde{x})e^{-\tau_C} + S_C(1 - e^{-\tau_C}), \quad (29) $$

where $I_C^*$ is the stellar continuum intensity at $\lambda_C$ and $\tau_C$ the optical depth value at $\tilde{x}$ on *p*, given by Eqs. (27) and (28).

While the similarity with the expression derived in the framework of the Sobolev theory is evident, we did not factor out the profile function, but integrated it explicitly over the region where the optical depth goes from zero to its maximum value, $\tau_C^*$, on the ray. In other words, we did not make use of the Sobolev approximation, so that Eq. (29) is equivalent to Eq. (25) in Rybicki & Hummer (1978) rather than to Eq. (26) of the same paper.

We next consider the additional contributions to the intensity on *p* caused by the neighboring lines.

*Contribution from B-line photons to line C.* As noted earlier, line *B* is displaced by $\delta v$ from line *C* such that line-B photons traveling in the *z* direction on a given impact parameter will be able to interact with atoms emitting transition *C* that are located at apparent velocities

$$ v_{\text{rad}}^{\text{min}} \le v_{\text{rad}} \le v_{\text{rad}}^{\text{max}} - \delta v \quad (30) $$

since these are Doppler-shifted to the rest frequency of line *C*. As illustrated by the top central panel of Fig. 6, line *C* will get a contribution from *B* photons that are emitted deeper in the envelope and are moving slower by an amount $\delta v$. This contribution is taken into account by additional terms in Eqs. (23) and (26). Naming $\delta x$ the frequency displacement corresponding to $\delta v$, the optical depth of line *C* along the line of sight is augmented by the quantity

$$ \tau_B(p, z, \tilde{x}) = \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}-\xi/2-\delta x)}^{z} \frac{\kappa_B(r')}{\Delta\nu_{D_B}(r')} \varphi(x')dz', \quad (31) $$

while the line intensity becomes

$$ I_C^{bc}(p, z, \tilde{x}) = \int_0^{\tau(\tilde{x})} \{S_C(p', z', \tilde{x}') \exp[-\tau_C(p', z', \tilde{x}')] + S_B(p', z', \tilde{x}') \exp[-\tau_B(p', z', \tilde{x}') - \tau_C(p', z', \tilde{x}')]\} d\tau(p', z', \tilde{x}') \quad (32) $$

since line *B* forms deeper than *C* (as seen by an outside observer). The subscript *bc* indicates that the expression for the *C* intensity now takes contributions of both *B* and *C* into account. The optical depth along impact parameter *p* for frequency $\tilde{x}$, from $z = z_{\text{max}}(p)$ inward, varies as<fnref n="4" />

$$ \tau_C^{bc}(p, z, \tilde{x}) = \begin{cases} 0, & z > z(\tilde{x} + \xi/2) \\ \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x} + \xi/2)}^{z} \frac{\kappa_C(r')}{\Delta\nu_{D_C}(r')} \varphi(x') dz', & z(\tilde{x} + \xi/2) \ge z > z(\tilde{x} - \xi/2) \\ \tau_C^*, & z(\tilde{x} - \xi/2) \ge z > z(\tilde{x} + \xi/2 - \delta x) \\ \tau_C^* + \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x} + \xi/2 - \delta x)}^{z} \frac{\kappa_B(r')}{\Delta\nu_{D_B}(r')} \varphi(x') dz', & z(\tilde{x} + \xi/2 - \delta x) > z \ge z(\tilde{x} - \xi/2 - \delta x) \\ \tau_C^* + \tau_B^*, & z < z(\tilde{x} - \xi/2 - \delta x). \end{cases} \quad (33) $$

We then integrate Eq. (32) to obtain the emergent intensity at $\tilde{x}$:

$$ I_C^{bc}(p, z = \infty, \tilde{x}) = I_C^*(\tilde{x})e^{-\tau_B - \tau_C} + S_C(1 - e^{-\tau_C}) + S_B(1 - e^{-\tau_B})e^{-\tau_C}, \quad (34) $$

where the notations are as defined in the previous paragraphs.

<footnote n="4">⁴ This assumes that *p* is close to the stellar core and $\tilde{x}$ is both positive and lies in the vicinity of line center, so that the entire frequency range of interest is present on the impact parameter.</footnote>
*Contribution from R-line photons to line C.* We next consider the contribution of the reddest of these three lines to C, which is displaced from C by a velocity amount $-\delta v^\dagger$ (or the corresponding frequency displacement $-\delta x^\dagger$), and interactions will occur with atoms moving at velocities

$$v_{\text{rad}}^{\text{min}} \le v_{\text{rad}} \le v_{\text{rad}}^{\text{max}} + \delta v^{\dagger}. \qquad (35)$$

While photons of line B contributing to C emerge from deeper layers than those of line C, we now have the opposite situation: line-R photons contributing to C originate from layers above C (as seen by an observer at $z = \infty$). The optical depth along the line of sight is increased by an amount

$$\tau_R(p, z, \tilde{x}) = \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}-\xi/2-\delta x^\dagger)}^{z} \frac{\kappa_R(r')}{\Delta v_{D_R}(r')} \varphi(x') dz', \qquad (36)$$

while the line intensity becomes

$$I_C^{\text{bcr}}(p, z, \tilde{x}) = \int_0^{\tau(\tilde{x})} \{S_B(p', z', \tilde{x}') \exp[-\tau_B(p', z', \tilde{x}') - \tau_C(p', z', \tilde{x}') - \tau_R(p', z', \tilde{x}')] \\ + S_C(p', z', \tilde{x}') \exp[-\tau_C(p', z', \tilde{x}') - \tau_R(p', z', \tilde{x}')] + S_R(p', z', \tilde{x}') \exp[-\tau_R(p', z', \tilde{x}')]\} d\tau(p', z', \tilde{x}), \qquad (37)$$

where the subscript bcr indicates that contributions of the three interacting lines to C are now taken in account. The optical depth variation along the ray is written, with the same caveat as before, as

$$\tau_C^{\text{bcr}}(p, z, \tilde{x}) = \begin{cases} 0, & z > z(\tilde{x} + \xi/2 + \delta x^\dagger) \\ \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}+\xi/2+\delta x^\dagger)}^{z} \frac{\kappa_R(r')}{\Delta v_{D_R}(r')} \varphi(x') dz', & z(\tilde{x} + \xi/2 - \delta x^\dagger) \ge z > z(\tilde{x} - \xi/2 + \delta x^\dagger) \\ \tau_R^*, & z(\tilde{x} - \xi/2 + \delta x^\dagger) \ge z > z(\tilde{x} + \xi/2) \\ \tau_R^* + \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}+\xi/2)}^{z} \frac{\kappa_C(r')}{\Delta v_{D_C}(r')} \varphi(x') dz', & z(\tilde{x} + \xi/2) > z \ge z(\tilde{x} - \xi/2) \\ \tau_R^* + \tau_C^*, & z(\tilde{x} + \xi/2) \ge z > z(\tilde{x} + \xi/2 - \delta x) \\ \tau_R^* + \tau_C^* + \frac{1}{\sqrt{\pi}} \int_{z(\tilde{x}+\xi/2-\delta x)}^{z} \frac{\kappa_B(r')}{\Delta v_{D_B}(r')} \varphi(x') dz', & z(\tilde{x} + \xi/2 - \delta x) > z \ge z(\tilde{x} - \xi/2 - \delta x) \\ \tau_R^* + \tau_C^* + \tau_B^*, & z < z(\tilde{x} - \xi/2 - \delta x). \end{cases} \qquad (38)$$

Integrating Eq. (37), we obtain the final emergent intensity of line C, which now takes all interactions into account:

$$I_C^{\text{bcr}}(p, z = \infty, \tilde{x}) = I_C^*(\tilde{x})e^{-\tau_B - \tau_C - \tau_R} + S_C(1 - e^{-\tau_C})e^{-\tau_R} + S_B(1 - e^{-\tau_B})e^{-\tau_C - \tau_R} + S_R(1 - e^{-\tau_R}). \qquad (39)$$

The line flux then follows from an integration of Eq. (39) overall impact parameters.

### 2.2.2. Lines B and R

As mentioned earlier, a reasoning similar to what has been done above for line C shows that formation of the B and R lines proceeds in exactly the same way, and the resulting expressions for the line intensities are very much the same as Eq. (39) above. The reason for this is the intrinsic symmetry of the problem, which is made very clear by Fig. 6. We can thus generalize Eq. (39) to

$$I_{B,C,R}^{\text{bcr}}(p, z = \infty, \tilde{x}) = I_{B,C,R}^*(\tilde{x})e^{-\tau_B - \tau_C - \tau_R} + S_B(1 - e^{-\tau_B})e^{-\tau_C - \tau_R} + S_C(1 - e^{-\tau_C})e^{-\tau_R} + S_R(1 - e^{-\tau_R}). \qquad (40)$$

Even though emergent intensities are basically the same for the three lines, their contributions to the emergent profiles are in fact quite different. This is because the frequencies to which each of the resonant regions contributes vary from line to line, as can also be seen in Fig. 6.

### 2.2.3. Line profiles

In order to show how the interactions affect the line profiles, we display them in two situations: in Fig. 7 the lines are computed independently, that is, without taking the interactions discussed above into account<fnref n="5" />. Unsurprisingly, the individual profiles display the characteristic features of lines formed in accelerating outflows, with absorption features reaching almost to the rest wavelength of the line and the sharp drop in intensity near zero velocity on the blue line side (see, e.g., Lamers et al. 1987). In Fig. 8, we show

<footnote n="5">⁵ All line profiles shown in the paper are computed with a numerical resolution of 1 km s⁻¹.</footnote>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Line profiles of the three lines considered and their blend for the wind model of Sect. 2 when the lines are computed independently.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Line profiles of the three lines considered and their blend for the wind model of Sect. 2 when radiative interactions between lines are taken into account.</figcaption></figure>

the emergent profiles when interactions are considered. Line C, which has the lowest optical depth of the three lines, is weak in the noninteracting case, but goes into strong emission when interactions are considered. Lines R and B have comparable optical depths, and are strongly enhanced by the interactions. The blends shown in these figures result from co-adding the individual line fluxes.

What can we conclude from this first example of line fluorescence? Clearly, line C is enhanced by photons of the two other lines in an outwardly accelerating outflow. So are B and R, and the resulting blend is a wide, triple peaked profile. The blue-displaced absorption characteristic of outflows is essentially that of the B line, somewhat filled up by emission. Given that the interacting line source functions are qualitatively quite similar to what they would be if the lines were independent (compare Figs. 4 and 5), the gist of line amplification occurs during flux integration. Actually, test computations demonstrated that line profiles obtained with source functions computed simply with the local Sobolev approximation are often qualitatively indistinguishable from those obtained with the interacting source functions used above. Again, this is likely an artifact introduced by the assumption of LTE atomic populations. But it is quite interesting to note that fluorescence caused by line amplification in the observer's frame occurs despite this limiting condition. This is of course a direct consequence of the additional emission in the resonance zones caused by Doppler-shifted neighboring lines in macroscopic velocity fields.

We next considered a more complex case, one with a different topology for the CRVSs: that of an accretion flow.

## 3. Extension to more complex velocity fields

### 3.1. Accretion flow

Free-falling accretion flows, as well as outwardly decelerating ballistic outflows or rotational flows, all lead to situations where distant parts of the emitting region can interact radiatively even when considering the case of single spectral lines; this is because the constant radial velocity surfaces discussed above are multi-valued along some rays. As we mentioned earlier, a generalization of the Sobolev theory to handle this case was worked out by Rybicki & Hummer (1978) for single lines. We recall their main results in the next paragraph before extending the same formalism to interacting lines.

#### 3.1.1. Sobolev source functions

We consider an accelerated gravitational accretion flow with velocity law

$$v(r) = -V_c \cdot (R_c/r)^\alpha, \qquad (41)$$

where $V_c$ is the free-fall velocity at $R_c$ (see Fig. 9). The computation parameters are summarized in Table 2. Since the maximum envelope gas density is 20 times lower than in the case discussed in the previous section, the line absorption coefficients of Fig. 10 are correspondingly smaller, but sufficient for producing sizable line emission in this velocity field, as seen below.

Since the envelope velocity is varying monotonically, the photons of the reddest considered line R are seen redshifted by atoms emitting bluer line photons, so that the source function of line R remains unaffected by the other transitions considered. We therefore start our discussion with transition R, to which the Rybicki & Hummer (1978) equations directly apply.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Velocity and density laws of the accretion flow discussed in Sect. 3.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Line absorption coefficient for the three lines considered.</figcaption></figure>

<table><caption>**Table 2.** Computation parameters for the accretion flow of Sect. 3.</caption><thead><tr><th>R<sub>c</sub></th><th>R<sub>env</sub></th><th>M<sub>*</sub></th><th>T<sub>eff</sub></th><th>T<sub>env</sub></th><th>n<sub>t</sub>(R<sub>c</sub>)</th><th>V<sub>c</sub></th><th>α</th></tr></thead><tbody><tr><td>2 R<sub>☉</sub></td><td>10 R<sub>c</sub></td><td>1 M<sub>☉</sub></td><td>4.5 × 10<sup>3</sup> K</td><td>6 × 10<sup>3</sup> K</td><td>1 × 10<sup>15</sup> cm<sup>-3</sup></td><td>436.76 km s<sup>-1</sup></td><td>0.5</td></tr></tbody></table>

**R transition.** The generalized Sobolev source function for the case of an equivalent two-level atom in an inwardly accelerating radial flow is given by

$$S_R(r) = \frac{(1 - \varepsilon_R)\bar{\beta}_R^c I_R^c + \varepsilon_R B_R + (1 - \varepsilon_R)S_R^{nl}}{\varepsilon_R + (1 - \varepsilon_R)\beta_R}, \quad (42)$$

with

$$S_R^{nl}(r) = \frac{1}{2} \int_{-1}^{+1} \frac{1 - e^{-\tau_R}}{\tau_R} (1 - e^{-\tau(r_R^*)}) S(r_R^*) d\mu, \quad (43)$$

where $\tau(r_R^*)$ and $S(r_R^*)$ are evaluated at the CP surface $r_R^*$ associated with $O(r)$. The third term of Eq. (42) represents the main nonlocal contribution to the source function that arises in these flows. A second nonlocal contribution appears in the term $\bar{\beta}_R^c I_R^c$, which represents the contribution of the stellar continuum flux at line wavelength, attenuated by the intervening CP layers between core and $O(r)$:

$$\bar{\beta}_R^c I_R^c = \frac{I_R^c}{2} \int_{\mu_c}^{+1} \frac{1 - e^{-\tau_R}}{\tau_R} e^{-\tau(r_R^*)} d\mu. \quad (44)$$

The left panel of Fig. 11 shows the CP surfaces associated with line $R$ (in this case the $H_\epsilon$ line), which can be compared to the CP surfaces found by Rybicki & Hummer (1978). These are self-CP surfaces in the sense that transition $R$ interacts radiatively with distant photons at the same rest wavelength.

The analogy of these equations with the case of interacting lines $B$ and $C$ (Eqs. (11)–(13)) of the previous section is obvious. In the previous case of transition $C$, the integration ran over a resonance surface corresponding to the locus of $B$ photons that were Doppler-shifted by the right amount, providing an additional source term that represents a net gain of photons for transition $C$. For a single line formed in an accretion flow, on the other hand, the local source term for $R$ is supplemented by nonlocal photons of the same transition $R$, originating from the CP surface associated with the considered emission point in the envelope.

**C and B transitions.** The expressions for $S_C(r)$ and $S_B(r)$ are similar to Eq. (42) but with the following changes in the nonlocal parts of the source function, reflecting the fact that the self-CP surface $r^*$ of the line is not the only region that must be taken into account in the nonlocal integral. Indeed, atoms emitting line $C$ will be able to absorb some blue-displaced line-$R$ photons, while atoms emitting transition $B$ will interact with blue-displaced photons from both $C$ and $R$. Figure 11, in its center and right panels, show the set of CP surfaces that are contributing to the nonlocal parts of the respective source functions. The graphical convention is the following: solid lines indicate line self-interactions, dashed lines show the CP surfaces corresponding to the nearest transition, and dotted lines indicate the CP surfaces corresponding to the third, farthest transition.
<figure><img src="image_12.png" /><figcaption>Fig. 11. Selected CP surfaces for the spherically symmetric accretion flow, with properties given in the text. Each surface is associated with the radial point of the same color.</figcaption></figure>

As previously, the condition of interaction for the C line is given by

$$v_{\text{rad}}(r_R) - v_{\text{rad}}(r_C) = \Delta\lambda_{RC}/c, \qquad (45)$$

which defines a CP surface denoted by $r'_R$ in the following expressions. With $\tau(r'_R) = \tau'_R$ and $\tau(r_C^*) = \tau_C^*$, we have

$$S_C(r) = \frac{(1 - \varepsilon_C)\bar{\beta}_C^c I_C^c + \varepsilon_C B_C + (1 - \varepsilon_C)S_C^{nl}}{\varepsilon_C + (1 - \varepsilon_C)\beta_C}, \qquad (46)$$

with

$$S_C^{nl}(r) = \frac{1}{2} \int_{-1}^{+1} \frac{1 - e^{-\tau_C}}{\tau_C} \left[ (1 - e^{-\tau_C})S(r_C^*) + (1 - e^{-\tau'_R})S(r'_R)e^{-\tau'_C} \right] d\mu \qquad (47)$$

and

$$\bar{\beta}_C^c = \frac{1}{2} \int_{\mu_c}^{+1} \frac{1 - e^{-\tau_C}}{\tau_C} e^{-(\tau_C^* + \tau'_R)} d\mu. \qquad (48)$$

For line B, two resonance conditions on the radial velocity lead to line interactions, namely

$$\begin{cases} v_{\text{rad}}(r_B) - v_{\text{rad}}(r_C) = \Delta\lambda_{BC}/c \\ v_{\text{rad}}(r_B) - v_{\text{rad}}(r_R) = \Delta\lambda_{BR}/c \end{cases} \qquad (49)$$

We must therefore locate the resonance regions corresponding to these conditions, the loci of which are noted $r_C''$ and $r_R'''$ in the following equations. In a direct extension of Eq. (15), we then have

$$S_B(r) = \frac{(1 - \varepsilon_B)\bar{\beta}_B^c I_B^c + \varepsilon_B B_B + (1 - \varepsilon_B)S_B^{nl}}{\varepsilon_B + (1 - \varepsilon_B)\beta_B}, \qquad (50)$$

with

$$S_B^{nl}(r) = \frac{1}{2} \int_{-1}^{+1} \frac{1 - e^{-\tau_B}}{\tau_B} \left[ (1 - e^{-\tau_B^*})S(r_B^*) + (1 - e^{-\tau_C''})S(r_C^*)e^{-\tau_B^*} + (1 - e^{-\tau_R''})S(r_R^*)e^{-(\tau_B^* + \tau_C'' + \tau_R''')} \right] d\mu \qquad (51)$$

and

$$\bar{\beta}_B^c = \frac{1}{2} \int_{\mu_c}^{+1} \frac{1 - e^{-\tau_B}}{\tau_B} e^{-(\tau_B^* + \tau_C'' + \tau_R''')} d\mu. \qquad (52)$$
<figure><img src="image_13.png" /><figcaption>Fig. 12. Source functions for the three lines when they are computed independently.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Source functions for the three lines when radiative interactions are taken into account.</figcaption></figure>

The interaction matrix for this case must include the self-interaction of each line, so that it becomes, with line ordering R, C, B,

$$I = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 1 & 1 & 1 \end{pmatrix}. \qquad (53)$$

Figure 12 displays the computed source functions for the three lines in consideration when they are treated as independent, that is, when Eq. (4) is used for computing all of them. Figure 13 shows the source functions computed when taking the radiative interactions into consideration, that is, when using Eq. (46) for the C line and Eq. (50) for the B line. The changes caused by core photon absorption at the CP surfaces are small for the optically thin Fe I line but lead to a strong depression of the optically thick Ca II line source function below R ~ 4. We recall that the H<sub>e</sub> line source function is the same in both cases.

### 3.1.2. Flux integration

In an accelerating infall or decelerating outflow, the unique correspondence between line frequencies and positions on the line of sight found in accelerating outflows is lost, as can be seen in the top left panel of Fig. 14, which displays the envelope’s gas radial velocities seen by an outside observer located at z = ∞. Depending on the impact parameter and envelope extent, the CRVSs are either open or closed, so that there are one or two locations on the line of sight that have the same frequency displacement. As explained in Bertout (1984), the flux integration method that we use handles this situation without modification. It finds all line frequency displacements present on a given impact parameter and orders them for the line intensity computation.

The line flux of R is made up of several contributions shown in the first column of Fig. 14. They are first the line-R photons originating in the CRVSs shown in the top panel, then the displaced C-line photons originating in deeper parts of the envelope (middle panel), and also the displaced B-line photons originating even deeper in the envelope (bottom panel). The final contribution to the red side of line R comes from C and B photons emitted by fast-moving atoms located in front of the star. As illustrated by the second column of Fig. 14, line C gets a contribution from R photons that are emitted in the outside part of the envelope (as seen by the observer) where atoms are moving slower by an amount δv, so the flux is enhanced by emission from the front part of the envelope. Line C gets an additional contribution from B photons that are coming from the deeper parts of the envelope and displaced by the right amount. We again note the additional contribution to the C-line red wing from B photons emitted by fast-moving atoms located in front of the star. Similar considerations hold for line B, and the various origins of its flux are shown in the right column of Fig. 14.

Figures 15 and 16 display the computed line profiles for two situations. In Fig. 15 the lines are computed as if they were single lines, without taking all interactions discussed above into account, while in Fig. 16 we show the profiles that result when interactions are considered. Profiles of Fig. 15 are typical accretion flow profiles with intermediate optical depths. Line C has the lowest optical depth of the three lines; it is indistinguishable from the continuum when the lines are computed independently but goes in moderate emission when interactions are taken into account. Line R has a somewhat lower optical depth than B but is enhanced by contributions of photons from the other transitions. The flux of these last two lines are similar, so the blend is in emission, with a single peak approximately centered on the strongest line R; it displays a moderately deep but wide absorption component that results from a combination of the three individual absorption features.
<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Regions contributing to the formation of the considered line profiles in the accretion flow of Sect. 4. See the caption of Fig. 6 for details.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Line profiles of the three lines considered and their blend for the accretion flow discussed here when the lines are computed independently.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Line profiles of the three lines considered and their blend for the accretion flow discussed here when radiative interactions between lines are taken into account.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Velocity and density laws of the non-monotonic accretion flow discussed in this section.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Line absorption coefficient and rates of collisional de-excitation for the three lines considered.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Selected CP surfaces for the spherically symmetric non-monotonic accretion flow, with properties given in the text. Each surface is associated with the radial point of the same color.</figcaption></figure>

## 3.2. Non-monotonic velocity fields

The method developed here readily extends to non-monotonic velocities, as found, for example, in an accreting envelope with an accretion shock close to the stellar surface. To illustrate this, we considered a qualitative TTS accretion shock toy model, the velocity and density fields of which are displayed in Fig. 17. The gas free-fall velocity increases toward the star with the usual $v(r) \propto r^{-1/2}$ radial dependence until it reaches the shock at $R_s$, assumed to be located at $1.2 R_c$. The shock is assumed to be isothermal, so the gas velocity decreases by a factor of 4 at $R_s$. This jump is followed by an (assumed) linear deceleration and reaches the stellar radius $R_c$ at $10 \text{ km s}^{-1}$. The density follows from the continuity equation, and the gas density at the shock is $n_t(R_s) = 5 \times 10^{15} \text{ cm}^{-3}$. We did not consider the complex physics involved in the shock, but simply solved the formation of our three lines with these values of velocity, temperature, and density fields.

For more realistic source functions, we introduced collisional de-excitation rates in an approximate way by using formulas from Allen's Astrophysical Quantities (Cox 2000), which actually go back to van Regemorter (1962). Clearly, more modern data should be used for a detailed comparison with observations, but at this preliminary stage a crude estimate of the collision effects is sufficient. The collisional de-excitation rates (dashed lines) and the absorption coefficients (solid lines) of the three lines are shown in Fig. 18. The collisional de-excitation rate of the calcium ion line is an order of magnitude higher than the rates of neutral hydrogen and iron lines, which are indistinguishable in this graph. The three lines are moderately optically thick in the post-shock region, where much of the fluorescent amplification occurs.

The CP surfaces for the three lines considered are shown in Fig. 19. As in previous panels showing CPs, the solid lines display the self-CP surfaces of each line, while the dashed lines show the interaction surfaces with the first neighboring line and the dotted lines display the interaction surfaces with the second neighboring line. We see that, unlike the monotonic cases, all lines are, as expected, involved in the interactions, and the CP surfaces are a mix of those found in the previous cases of monotonic flows. They introduce discontinuities in the source functions, particularly in the shock region, as can be seen in Fig. 20.

The resulting line profiles are displayed in Fig. 21. They do not resemble those found for the accretion flow of the previous section, partly because most of the flux originates in the post-shock deceleration region close to the star, and partly because of
<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Source functions of the three lines considered and their blend for the non-monotonic accretion flow discussed here.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Line profiles of the three lines considered and their blend for the non-monotonic accretion flow discussed here.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 22.** Parameter space study for the net EWs of FeIλ3969 line and its blend with H and H<sub>ε</sub>, in the case of an accretion flow with ε ≠ 0.</figcaption></figure>

the collision-term effect. Despite the fact that the optical depths of the H and CaII lines are very similar, the emergent profiles are quite different, with CaII H showing a steep central peak caused by the collision term that enhances its source function in the innermost envelope region (see Fig. 20). We note that the three lines participate in the fluorescence process, with the iron line getting contributions from both neighboring lines and thus being amplified most.

**4. Discussion and conclusion**

We have investigated a rather small parameter space in the previous sections of this paper, mainly because we required that CaII H and K have comparable strengths, as is observed in many TTSs. As discussed in Sect. 2, this requires – when LTE and isothermal envelope are assumed – that gas temperatures be restricted to a limited range around 6000K, which in turn implies, when ignoring collisional excitation, high gas densities at the envelope's bottom and resulting mass flow rates that appear unrealistically high for TTSs.
<table><caption>**Table 3.** Computation parameters for the parameter space study discussed here.</caption><thead><tr><th>R<sub>c</sub></th><th>R<sub>env</sub></th><th>M<sub>*</sub></th><th>T<sub>eff</sub></th><th>T<sub>env</sub></th><th>n<sub>t</sub>(R<sub>c</sub>)</th><th>V<sub>c</sub></th><th>α</th></tr></thead><tbody><tr><td>2 R<sub>☉</sub></td><td>10 R<sub>c</sub></td><td>0.5 M<sub>☉</sub></td><td>3.0 × 10<sup>3</sup>K</td><td>10<sup>12</sup>–10<sup>14</sup> cm<sup>-3</sup> K</td><td>6 × 10<sup>3</sup>–9 × 10<sup>3</sup> K</td><td>436.76 km s<sup>-1</sup></td><td>0.5</td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>**Fig. 23.** Source functions of the three lines considered for an accretion flow with $n_t(R_c) = 6 \times 10^{12}$ cm$^{-3}$ and $T_{env} = 8.5 \times 10^3$ K when collisional de-excitation rates are included.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Fig. 24.** Line profiles of the three lines considered for an accretion flow with the properties given in the previous figure caption.</figcaption></figure>

However, fluorescence is not limited to this narrow range of parameters, but can also take place at higher temperatures and at lower densities when we include the collisional de-excitation terms in the source functions. This is because the collision term enhances the source function, leading to sizable flux for low line optical depths. Figure 22 shows the result of a parameter space exploration that includes collisional de-excitation rates in the source function computation. There, we display the net equivalent widths (EWs) of both FeLλ3969 and the blend of the three lines, as computed in a more extended parameter space investigation. We vary $n_t(R_c)$ and $T_{env}$ in the ranges indicated in Table 3, which also gives the other computation parameters. The collisional terms are included in the approximation discussed earlier. The net EW<fnref n="6" /> of a line is defined as its computed EW when the computation includes interactions with the other two lines, from which we subtracted a baseline EW computed under the assumption that the three lines are independent. The net EW is thus a measure of the fluorescent flux.

In Fig. 22 we plot the logarithm of the net EWs, which ranges from -2 to +1.5 in the parameter space investigated. We find that a significant FeLλ3969 fluorescent flux (EW<sub>net</sub> ≳ 1 Å) is produced over a wide range of temperatures for $n_t(R_c) ≳ 4 \times 10^{12}$ cm$^{-3}$. The interacting lines at these relatively low densities are all optically thin, with pure, almost featureless emission profiles driven by the collision terms; the sources and profiles of a representative example are shown in Figs. 23 and 24. The minimum density above corresponds to a spherical mass flow rate of ≈$10^{-6}$ M<sub>☉</sub> yr$^{-1}$. This mass accretion rate value remains quite high for TTSs in general, but it may represent a plausible value for the most active members of the class (Hartigan et al. 1995). One should also recall that TTSs are highly variable, which indicates that the accretion process is not steady (cf. Bouvier et al. 2007). Fluorescence might therefore occur in periods of increased accretion but disappear during more quiescent phases, as observed recently in VW Cha (A. Armeni, priv. comm.).

As an absolute calibration of TTS mass accretion rates is notoriously difficult (see the discussion in Gullbring et al. 1998), it seems worthy to improve the atomic model and try computing the precise threshold mass accretion rate for fluorescence to set in. This will provide a useful new comparison tool for other mass flow rate determination methods.

With this work, we are now in a position to compute the formation of fluorescent lines in spherically moving media. However, as already mentioned, we are still far from being able to meaningfully compare theoretical line profiles to observed TTS ones. The reasons are threefold: (a) velocity fields in TTS envelopes are not spherically symmetric and so require a 3D analysis; (b) the atomic line level populations, which are dominated by non-LTE effects, must be computed iteratively with the radiation field; and (c) stellar envelopes are not isothermal. Point (a) can be addressed relatively easily by extending the present formalism to axisymmetric flow geometries. Point (b) requires setting up an iterative procedure that couples levels equations and the radiation field; several convergence acceleration methods have been proposed (cf. Auer 2009) and can be used to solve this problem. Point (c) is easily remedied; actually, a version of our code already makes provisions for varying envelope temperatures. Another aspect requiring

<footnote n="6">⁶ By convention, emission line components have positive EWs here.</footnote>
improvement for further progress is the approximation used for the collisional de-excitation rates, which should be replaced by more recent data. This will be a direct consequence of setting up realistic atom models for addressing point (b).

*Acknowledgements.* It is a pleasure to thank Sylvie Cabrit for useful discussions and Joli Adams for improving the English. I am grateful to an anonymous referee for useful comments and suggestions that improved the paper’s presentation.

## References

Armeni, A., Stelzer, B., Claes, R. A. B., et al. 2023, A&A, 679, A14

Auer, L. H. 2009, in *Numerical Radiative Transfer* (Cambridge, UK: Cambridge University Press), 101

Bastian, U., Bertout, C., Stenholm, L., & Wehrse, R. 1980, A&A, 86, 105

Bertout, C. 1984, ApJ, 285, 269

Bertout, C., Wolf, B., Carrasco, L., & Mundt, R. 1982, A&AS, 47, 419

Bouvier, J., Alencar, S. H. P., Harries, T. J., Johns-Krull, C. M., & Romanova, M. M. 2007, in *Protostars and Planets V*, eds. B. Reipurth, D. Jewitt, & K. Keil, 479

Bowen, I. S. 1934, PASP, 46, 146

Bowen, I. S. 1935, ApJ, 81, 1

Cabrit, S., & Bertout, C. 1986, ApJ, 307, 313

Cabrit, S., & Bertout, C. 1990, ApJ, 348, 530

Cabrit, S., & Bertout, C. 1992, A&A, 261, 274

Castor, J. I. 1970, MNRAS, 149, 111

Cox, A. N. 2000, *Allen's Astrophysical Quantities* (New York: AIP Press)

Dmitriev, D. V., Grinin, V. P., & Katysheva, N. A. 2019, Astron. Lett., 45, 371

Eriksson, M., Veenhuizen, H., Wahlgren, G. M., & Johansson, S. 2004, in *Revista Mexicana de Astronomia y Astrofísica Conference Series*, 21, eds. C. Allen, & C. Scarfe, 132

Goldreich, P., & Kwan, J. 1974, ApJ, 189, 441

Gullbring, E., Hartmann, L., Briceño, C., & Calvet, N. 1998, ApJ, 492, 323

Hamann, W. R. 1981, A&A, 93, 353

Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736

Hartmann, L., Avrett, E., & Edwards, S. 1982, ApJ, 261, 279

Herbig, G. H. 1945, PASP, 57, 166

Hubeny, I., & Mihalas, D. 2014, *Theory of Stellar Atmospheres* (Princeton: Princeton University Press)

Huggins, M. L. 1898, ApJ, 8, 54

Hyung, S., Lee, S.-J., & Lee, K. H. 2018, *J. Astron. Space Sci.*, 35, 7

Joy, A. H. 1945, ApJ, 102, 168

Kastner, S. O. 1991, Ap&SS, 185, 265

Kurosawa, R., & Romanova, M. M. 2012, MNRAS, 426, 2901

Lamers, H. J. G. L. M., Cerruti-Sola, M., & Perinotto, M. 1987, ApJ, 314, 726

Olson, G. L. 1982, ApJ, 255, 267

Pavlakis, K. G., & Kylafis, N. D. 1996, ApJ, 467, 292

Pereira, C. B., de Araújo, F. X., & Landaberry, S. J. C. 1999, MNRAS, 309, 1074

Rybicki, G. B., & Hummer, D. G. 1978, ApJ, 219, 654

Selvelli, P., Danziger, J., & Bonifacio, P. 2007, A&A, 464, 715

Sobolev, V. V. 1960, *Moving Envelopes of Stars* (Cambridge: Harvard University Press)

Tambovtseva, L. V., Grinin, V. P., Rodgers, B., & Kozlova, O. V. 2001, Astron. Rep., 45, 442

Tambovtseva, L. V., Grinin, V. P., & Weigelt, G. 2014, A&A, 562, A104

Thackeray, A. D. 1937, ApJ, 86, 499

Valenti, J. A., Basri, G., & Johns, C. M. 1993, AJ, 106, 2024

van Regemorter, H. 1962, ApJ, 136, 906

Willson, L. A. 1972, A&A, 17, 354

Willson, L. A. 1974, ApJ, 191, 143

Willson, L. A. 1975, ApJ, 197, 365
Appendix A: Numerical implementation

SLIM2 (Spectral Line Interactions in Moving Media) is a publicly available<fnref n="7" /> Python/NumPy code for solving the line formation problem in 2D moving media, for multiple lines stemming from multiple elements, in the special case of interacting lines leading to fluorescence effects. In the framework of the well-tested hybrid approach to the problem of line formation described above, it first computes the generalized Sobolev source functions and then integrates the line flux exactly.

The code makes extensive use of the NumPy ndarray((dim1, dim2), object) 3D construct, which is particularly useful for ray-
tracing procedures, although it is frowned upon by some pythonistas because if complicates array broadcasting.

Version 1.0 presented here is stable and optimized for solving the problem of forming the three-line blend around 3969Å in TTS spectra. The atomic level populations are LTE and remain constant during the iterative procedure used to compute the source functions. Plotting routines are restricted to this specific three-line case. They are included in the project so that interested users will be able to reproduce easily this publication’s figures, thus gaining experience with the code before starting on their own developments.

A more general version of the code, which is able to compute many interacting and noninteracting lines simultaneously and
allows for variable temperature in the envelope, will be made available at a later date.

The following modules are included in the GitHub TTS Fluorescence package:

* **Fluorescence_Main.py** contains a single function **main_module** that calls all other modules in the following order to perform the computation and to store computational parameters, source functions, and line profiles in a subdirectory of **results** directory marked with a date-time stamp of the start of computation.

* **import_parameters.py** imports all input parameters from **computation_parameters.conf**. For this, an ancillary module **parameter_parser.py** is needed. Several different options are possible for the computation (local or nonlocal Sobolev source, noninteracting or interacting lines, Sobolev or exact flux mode, single model or grid of models, various velocity fields, several levels of logging, etc.) so that functional parameters are numerous. There are also a few physical parameters, such as the maximum density, envelope temperature, and stellar properties, that must be defined. The meaning of each of these parameters is explained in the .conf file and default values are indicated.

* `setup_grids.py` sets up the spatial, velocity, and density grids and analyzes the velocity field structure.

* **import_line_data.py** reads the line data from **import_atomic_data.py** and computes and stores as arrays the line optical depths using **LTE_module.py**. It also computes the absorption line profiles, defines the line interaction matrices, and sets up the overall frequency grid for the flux integration. Databases for partition functions and abundances are stored in directory **input_data**.

* **source_module.py** first finds and orders all the CP surfaces discussed above on high-resolution angular grids. This step is a major CPU-time consumer, and we store the CP-surface geometries for repeated use when a string of models with the same geometrical properties are calculated. All integrals over solid angles that appear in the computation of the source function are solved by Gaussian quadratures. The local source functions are computed and iterated with the nonlocal contributions until convergence occurs. Although there are often more than two resonant surfaces on each line of sight, convergence occurs relatively rapidly for the velocity fields considered above and for the many others that we tested.

* **flux_integration.py** integrates the line intensities over the impact parameter grid. A spatial grid is defined on each impact parameter and is transformed by interpolation into a grid of line frequencies present on the given ray. Contributing frequencies to each line are then found by taking their velocity displacements into account, and resonant regions are deduced from there. Incremental line intensities are then computed, and the flux follows from integrating the emergent intensities.

* `function_library.py` contains ancillary functions for `source_module.py` and `flux_integration.py`.

* `in_out_library.py` contains functions for storing the results in a subdirectory of `results` and reading them for future use.

* Finally, **plot_library.py** reads stored results and plots various quantities, storing the plots for future reference in the same subdirectory of **results** (see above).

All these modules are heavily commented to make it easier for the user to follow the procedures used for solving the problem.

The computing time for three interacting lines with high-resolution grids (100 radial points, 2000 angles, 100 impact parameters, and 2000 frequencies across the line profiles) is typically on the order of 30 seconds on an Apple MacBook Pro with M1 Pro silicon running PyCharm.

<footnote n="7">⁷ under Creative Commons 4.0 BY-NC-SA license. The code is available at https://github.com/claude-bertout/SLIM2</footnote>