© ESO, 2013

1. Introduction

Although massive stars (with masses higher than about 8 M_⊙) play a key role in the evolution of the Universe, their late stages of evolution as well as the mechanisms that lead to their formation and early evolution are still poorly known (Zinnecker & Yorke 2007). Massive stars ionize their surroundings at the early stages of their evolution and also at some evolved stages (such as supergiant B[e] stars) forming Ultra-Compact HII regions (UC-HII). Statistical analysis of the number of UC-HII regions has shown that these objects seem to be confined for longer times than those derived from their crossing times, if it is assumed that they constantly expand at the speed of sound (Churchwell et al. 1990). This raises the question of why high-pressure HII regions do not quickly expand into their surrounding environments of lower pressure.

Hollenbach et al. (1994) proposed that the long lifetime of UC-HII regions can be explained if they are formed by photoevaporation of the neutral material from the rotating accretion disk left after star formation. This hypothesis would be consistent with two observational results: i) the substantial fraction of UC HII regions (~30%) characterized by broad hydrogen-radio recombination lines (with linewidths of 70−200 km s^-1), which are interpreted as arising from low velocity ionized outflows associated to neutral disks (Jaffe & Martín-Pintado 1999); and ii) the high spatial resolution observations of massive protostars, which have shown that massive star formation seems to proceed, at least in some cases, through accretion of material from a circumstellar disk (Nielbock et al. 2007; Jiménez-Serra et al. 2009; Kraus et al. 2010; Davies et al. 2010; Preibisch et al. 2011). Evidence of circumstellar disks have also been found toward evolved massive stars such as supergiant B[e] stars (Miroshnichenko et al. 2005). However, the kinematics of circumstellar disks and ionized envelopes around pre- and post-main sequence massive stars are still poorly understood.

Since its discovery, the UC-HII region of MWC349A has turned out to be a key object for the study of the physical and kinematical properties of these regions thanks to its peculiar characteristics. The star MWC349A is the most intense radio-continuum source at centimeter wavelengths (Braes et al. 1972) and also one of the most intense sources in the NIR and mid-IR (Tafoya et al. 2004). The position of this source in the HR diagram clearly corresponds to a massive star (Kraus 2009). According to its spectral features, it is classified as a B[e] star, which are characterized by their forbidden emission lines ([Fe II] and [O I]) in the optical spectrum and its NIR excesses (Lamers et al. 1998). Nevertheless, its evolutionary stage, as often happens with B[e] stars, is under dispute because of the difficulty caused by the overlap of the pre- and post-main sequence B[e] evolutionary tracks in terms of effective temperature and luminosity (Marston & McCollum 2008) and the lack of optical absorption lines (Hamann & Simon 1988; Andrillat et al. 1996). However, it is one of the few massive stars with a well-established circumstellar disk (White & Becker 1985; Danchi et al. 2001) and an ionized outflow expanding at nearly constant velocity (Olnon 1975; Tafoya et al. 2004).

The most peculiar characteristic of MWC349A that allows us to study in detail the kinematics of its ionized regions is that, until the detection toward η Carinae (Cox et al. 1995), Cepheus A HW2 (Jiménez-Serra et al. 2011) and MonR2-IRS2 (Jiménez-Serra et al. 2013), it was the only known celestial source with hydrogen radio-recombination line (RRL) masers at millimeter wavelengths (Martín-Pintado et al. 1989a). The strong maser amplification of mm and submm RRLs and the strong radio-continuum emission of MWC349A offer a unique possibility to obtain spectra and images with a high spectral and spatial resolution, down to 5 mas (Planesas et al. 1992; Martín-Pintado et al. 2011). The comparison of the H30α radio-recombination maser line profile, the H30α centroid map and the free-free radio-continuum emission with the predictions from a non-local thermodynamic equilibrium (non-LTE) 3D radiative transfer MOdel for REcombination LInes (MORELI), has provided stringent constraints on the disk and outflow kinematics of MWC349A as well as its electron temperature and electron density distribution (Martín-Pintado et al. 2011). The observations were explained by assuming two kinematical components, a thin virtually edge-on Keplerian ionized disk around a massive star (~30−60 M_⊙), responsible for the maser spikes observed in the H30α, and a rotating and expanding ionized wind with a terminal velocity of about 70 km s^-1. It also led to the remarkable conclusion that the outflow seems to be launched from the Keplerian disk at a radius of ≤20 AU from the star. Magnetic wind models (both X-wind models, Shu et al. 1994; and disk wind models, Blandford & Payne 1982) seem to be the most suitable to explain the region where the wind emerges and accelerates up to its terminal velocity. These results are extremely important since they provide clues on the origin of ionized winds from massive stars with circumstellar disks at very small spatial scales.

This paper gathers all the centimeter, millimeter, submillimeter, and far-infrared observations of the continuum and RRLs toward MWC349A to establish, for the first time, a model that reproduces the bulk of the observed features found over this wide frequency range. We carried out a thorough analysis of all the data by using a 3D non-LTE radiative transfer model to constrain the structure, physical conditions, and the kinematics. In Sect. 2, we discuss the basic physics incorporated in the MORELI radiative transfer model required for the interpretation of both the continuum and RRL emission. In particular, we describe the line formation under non-LTE conditions and the effects of saturation of masers. In Sect. 3, we describe the different geometries, physical and kinematical structure included in MORELI. Then, in Sect. 4, we report the first model for MWC349A that consistently explain most of the observational data available to date. Also, in Sect. 4, we show how the main parameters of the ionized gas were constrained by comparing the observations with the predictions of MORELI. While in Sect. 4 we focus on the H30α centroid map and its line profile, as well as the integrated-line intensities of RRLs from cm to the mid-IR wavelengths, in Sect. 5 we show the fit of our predictions for the rest of RRLs observed so far. Finally, conclusions are drawn in Sect. 6.

2. Non-LTE radiative transfer model

2.1. Radiative transfer integration

The resulting spectrum for every line-of-sight is obtained by discretizing the sources in regular 3D grid cells with sizes of dx, dy, and dz, where the z axis refers to the axis along the line-of-sight (with the observer toward negative z values), x the revolution symmetry axis projected into the plane of the sky, and y the axis orthogonal to the x − z plane. Thus, this discretization observed in the plane of the sky, (x,y), is equivalent to a mesh of points. We solve the radiative transfer equation for every cell along the line-of-sight by assuming constant properties in the whole cell volume (1)where I_ν is the specific radiation intensity at the frequency ν, S_ν is the source function, and τ_ν is the optical depth (defined as dτ_ν = κ_νdl, with κ_ν the absorption coefficient of the cell with physical length l).

We consider, as an approximation, that the electron density of the source is large enough to assume that the diffuse radiation is locally absorbed leading to a subsequent photoionization, the so-called on-the-spot approximation (Osterbrock 1989). The outgoing radiation of cell i for a given frequency, ν, is then given by (2)where I_ν,i − 1 is the background emission (which is given by the cosmic microwave background radiation at the rear side of the source where i = 1), B_ν(T_e) the source function (given by Planck’s law for a black body with electron temperature T_e), and τ_ν,i the optical depth from the cell i.

2.1.1. Continuum emission

The optical depth for the free-free continuum radiation is calculated by using the absorption coefficient (3)where m_e and e are the electron mass and charge, h and k the Boltzmann and Planck constants, c the speed of light, ϵ₀ the vacuum permittivity, N_e and N_i the electron and ionic densities and g_ff,ν and g_bf,ν the Gaunt factors for free-free and bound-free transitions, respectively (see Appendix A).

2.1.2. LTE RRL emission

In the formation of RRLs, one has to simultaneously consider the contribution of both the continuum and line emission. According to the quantum theory of the radiation, spontaneous emission, stimulated emission, and absorption are the three basic processes through which photons and atoms interact (Einstein 1916). The probability per unit time for an electron to spontaneously decay from the upper level m to the lower level n is given by the Einstein coefficient of spontaneous emission, A_mn (Towle et al. 1999). On the other hand, we use the spectral energy density, u_ν, to define the Einstein B-coefficients. The transition probability per unit time for absorption and stimulated emission are P_a = u_νB_nm and P_st = u_νB_mn, respectively, where u_ν is related to I_ν(4)and where Ω is the solid angle.

Using the Einstein coefficients described previously and treating the stimulated emission as a negative absorption, the line absorption coefficient depends on the electron population of the involved levels, N_m and N_n, according to the expression (5)where B_nm is the Einstein coefficient for absorption, g_n and g_m are the statistical weights of the electronic levels n and m respectively (with g_k = 2k² for k = n,m), and Φ_ν the line profile (see Appendix B).

The values of N_m and N_n are determined by solving the equation system of statistical equilibrium considering all the processes involved in the excitation of electronic levels. Under LTE conditions, collisions dominate over the radiative processes and the electronic levels are populated according to the Saha-Boltzmann distribution. In this case, Eq. (5) yields: (6)where B_nm is estimated according to the Dirac theory described by Towle et al. (1999).

2.1.3. Non-LTE RRL emission

The absorption coefficient discussed in the previous section must be modified for the case of non-LTE since the electron population departs from that given by the Saha-Boltzmann distribution as a result of radiative transitions dominating over collisional transitions. To account for this one introduces for each electronic level, n, a correction factor known as the departure coefficient, b_n, that relates the electron population in the non-LTE case, N_n, with the LTE case, , (7)These departure coefficients depend on both collisional and radiative processes. They have been computed independently by several authors taking into consideration different ranges of electron density and electron temperature. To date, the two most extensive b_n coefficient tables are those reported by Storey & Hummer (1995) and Walmsley (1990). In Fig. 1 we show their values for an ionized gas of 10 000 K for different electron density ranges. The MORELI code incorporates the parametrization of the b_n values obtained in both studies for the case of an optically thin HII-region for all RRLs except for the Lyman lines. We remark that this procedure will give us approximate results since the departure coefficients would actually depend on the radiation field for our particular assumed geometry and physical structure, which will differ from those assumed in the models used in Storey & Hummer (1995) and Walmsley (1990).

As described by Dupree & Goldberg (1970), the RRL absorption coefficient can be then expressed as in the non-LTE case, where β_mn is defined (8)For certain combinations of n, T_e, and N_e found in some UC-HII regions (i.e. in MWC349A and Cep A HW2), the level populations can be inverted and, subsequently, the second term of the latter equation can be larger than unity. It implies that β_mn < 0 and, therefore, the incident radiation is amplified. We can see in Fig. 2 that this is the case for an ionized gas with electron density of 10⁷ cm^-3 and electron temperature of 10 000 K as found in MWC349A.

Using the above mentioned departure coefficients, it is possible to derive the non-LTE source function. And from the latter, by integrating the radiative transfer equation, one obtains the brightness temperature at the center of the line (Dupree & Goldberg 1970) (9)where T_bg is the background temperature respectively. This equation reduces to Eq. (2) for the LTE case (b_n = 1).

2.2. Stimulated and maser emission

Soon after the discovery of RRLs, it was noticed that for the typical physical conditions found in ionized regions, the RRLs should be emitted under non-LTE conditions (Goldberg 1966), mainly stimulated emission. This type of emission is due to the amplification of the background continuum and line emission when the optical depth is negative but larger than − 1, that is, − 1 < τ_ν,c + τ_ν,l < 0. Nevertheless, if the pumping mechanisms producing the inversion of population in the ionized gas are intense enough, the amplification of the radiation can become very intense, τ_ν,c + τ_ν,l ≪ −1. This is referred to as maser amplification. In this case, the line intensity (Eq. 9) would have an approximately exponential dependence on optical depth and, therefore, on physical conditions (electron density and electron temperature) (10)Nevertheless, the exponential growth of the intensity along a line of sight with the required physical conditions for maser amplification cannot continue indefinitely. This is because the efficiency for stimulated emission reaches an upper limit when each pump event results in the release of a maser photon (Strelnitski et al. 1996b). In this case the intensity increases linearly with the optical depth and it is referred to as a saturated maser.

2.3. Saturation effects

The full treatment of saturation effects should take into account the changes in the electron populations (the b_n coefficients) induced by the maser radiation from all directions. This full treatment is clearly out of the scope of this paper. As a first approximation, we can estimate the saturation effects by considering only the changes in b_n coefficients caused by the incident radiation in the line-of-sight. From the resolution of the equations of statistical equilibrium under a particular set of approximations, we obtain the absorption coefficient for the saturated case (Strelnitski et al. 1996b) (11)where κ_ν is the absorption coefficient for the unsaturated regime (J_ν ≪ J_ν,sat), J_ν is the intensity averaged over the whole solid angle (), and J_sat is the saturation intensity, which is defined as (12)with C_t the collisional coefficient from the maser levels to other levels and C_mn the collisional coefficient between the two involved maser levels. The sum C_t + C_mn is derived from Table B.1 (C_t + C_mn = Δν_l/⁽2π⁾).

For strong maser emission (J_ν ≫ J_ν,sat), the absorption coefficient given in Eq. (11) becomes κ_ν,sat ≈ κ_νJ_ν,sat/J_ν. Substituting this into the radiative transfer Eq. (1) results in an equation whose integration gives a linear dependence of the outgoing intensity with the optical depth (13)where τ_ν,l, sat is the line optical depth taking into account the saturation effects (see Eq. (11)) and T_sat the saturation temperature.

To estimate the saturation effects, we use MORELI to compute the saturation degree, J_ν,sat/J_ν, assuming a solid angle for the maser beam, Ω_m.

3. Physical structure and kinematic considerations for the source model

3.1. Geometry

The MORELI code offers the possibility to consider different source geometries for the wind, such as a sphere, a cylinder (specified by its radius and length), a double-cone (specified by the semi-opening angle, θ_w, and length) or a hyperboloid (specified by its minimum width, semi-opening angle and length). In all cases MORELI incorporates the possibility that the strong radiative pressure coming from the central star has swept up all the ionized gas in the region within a radius r_min. For all the different geometries except for the sphere, one also needs to specify the inclination angle, θ_i, of the revolution axis with respect to the plane of the sky.

Although physically an unbounded ionized wind could be extended virtually to infinity, with no well definite boundaries, our model limits the size to a maximum radius, r_max, for all the geometries. This radius is one of the critical parameters in our model to adequately sample most of the flux in the RRLs tracing the outer regions of the ionized outflow. The MORELI code allows this radius to be introduced by the user or to be estimated from the effective radius, R_eff, defined to contain a specific percentage of the free-free radio-emission produced for an isothermal spherical ionized wind isotropically expanding at constant velocity with electron density distribution given by (see Sect. 3.4). Since we do not consider only this particular wind structure for the UC-HII regions, we essentially trace all the emission by considering that the source has a size n_R times that of the effective radius, that is (14)For a partially optically thick wind, the derived effective radius containing half of the free-free radio-continuum flux is given (Panagia & Felli 1975) (15)where μ is the main atomic weight per electron and Ṁ the mass loss rate in units of M_⊙/year. Instead, for an optically thick wind, the effective radius containing the percentage p of the free-free radio-continuum flux results (see Appendix C) (16)While the cylinder geometry is adequate to model sources with very collimated radio-jets (i.e. CRL618; Martín-Pintado et al. 1988), the double-cone and the hyperboloid geometries are adequate to model sources with wide ionized bipolar outflows.

3.2. Modelling of the kinematics

3.2.1. Disk kinematics

For the particularly interesting case of a source with a neutral disk and an expanding ionized outflow represented by a double-cone geometry, the model allows us to consider a kinematic component consisting in an ionized Keplerian rotating layer located next to the neutral disk. The ionized layer will be specified by its opening angle relative to the conical surface of the double-cone, θ_d (see Fig. 3), and the radius up to which the Keplerian rotation extends, r_d. Hereafter, we will refer to this layer as the ionized Keplerian rotating disk. The velocity rotation component is added to the expansion, being its projection along the line-of-sight (see Appendix D.1) given by (17)where V_Kepler is the Keplerian velocity at a radius equal to unity (), M the central mass of the star, and G the universal gravitational constant.

3.2.2. Outflow kinematics

The model also considers different possibilities for the kinematics of the ionized gas depending on the geometry and on the features of the source to model. For all geometries it is possible to consider the kinematics given by the Hollenbach model (Hollenbach et al. 1994). It is included in a simplified way considering two discontinuous kinematic components, a Keplerian rotating disk confined in its gravitationally bounded ionized atmosphere, and a rotating and expanding wind at the terminal velocity, v₀, outside. The radius at which the discontinuity occurs (known as gravitational radius, R_g) is defined as the distance at which the most likely thermal velocity equals the orbital velocity of the Keplerian disk due to the star’s gravitational potential. Subsequently, this radius is given by the expression (Hollenbach et al. 1994) (18)The required value for R_g to explain the observations can be also considered as a free input parameter in the model. The procedure applied is an approximation since there is not a specific radius from which the ionized gas can escape, but a continuous range of radii where the different individual ions with different velocities reach the escape velocity.

Another possibility incorporated in MORELI is to consider a constant acceleration for the wind, reaching its terminal velocity at a radius r_a. This is the case of the modelling of Cep A HW2 (Jiménez-Serra et al. 2011).

The projected velocity along the line-of-sight for those radii where the wind has reached its terminal velocity (r > r_a) or for the case when it is assumed that the wind reaches the terminal velocity suddenly at r_a = 0 is given by (19)While for those radii where the wind is still being accelerated (r < r_a) the expression results (see Appendix D.2) (20)%Finally, we note that in addition to the expansion, we can also consider another kinematic component for the ionized outflow: its rotation following a Keplerian law similar to that of the disk.

3.3. Electron temperature and electron density distribution

The MORELI code provides the possibility of considering electron temperature, T_o, and electron density, N_e, gradients specified by power-law with indices, b_t and b_d respectively. The values of T_o and N_e are specified in the model at a particular radius. Moreover, for the case of the double-cone geometry, we also incorporate a possible dependence of the electron density on the angle with respect to the revolution axis of the cone, θ (see Fig. 3), since physically we expect the higher density close to the neutral disk if it is photoevaporating. Thus, the electron density distribution assumed in the model is given by (21)where θ_a ≡ θ_w + θ_d is the semi-opening angle of the whole ionized gas, θ_a − θ is the angle between an ionized gas cell and the conical surface of the neutral disk (see Fig. 3) and θ₀ is a free dimensionless input parameter used to model the angular electron density dependence (θ₀ = ∞ if one assumes there is no dependence of N_e on the angle).

Although Eq. (21) refers to both the ionized outflow and ionized circumstellar disk located next to the neutral disk, the model also considers the possibility that this disk has a different electron density, N_e,d, and/or electron temperature, T_d, than the ionized wind to account for the possible difference in cooling due to its higher electron density.

The geometry of the source, and the electron density and electron temperature distributions are constrained by fitting the radio-continuum images and the observed spectral energy distribution (SED).

4. The case of MWC349A

In the following sections we describe in detail how the different observations have been used to constrain the values of the input parameters of the model for the ionized gas in MWC349A. Table 1 shows the parameters that fit best the overall data set available to date at the radio and millimeter wavelengths. The results presented in this paper have been derived by assuming a double-cone geometry with two kinematical components (ionized outflow and ionized circumstellar disk) and with sudden acceleration of the ionized wind. Since the outflow is partially optically thick, we have assumed the effective radius derived by Panagia & Felli 1975 (see Sect. 3.1). Although the distance for MWC349A is under debate (Meyer et al. 2002), we have considered the typically assumed distance of 1.2 kpc (Cohen et al. 1985).

The best fit was obtained by using the b_n coefficients reported by Storey & Hummer 1995 (see Sect. 4.2).

4.1. Constraining the physical structure and geometry of MWC349A on the basis of the SED and radio-continuum images

The radio-continuum maps obtained with the best angular resolution (at 7 mm, see Fig. 4) clearly show a bipolar morphology; there is a prominent dark lane in the East-West direction with a position angle of 8° (measured counterclockwise from East). The lack of free-free emission is interpreted as being caused by a neutral circumstellar disk. This interpretation is also supported by the dust emission from the disk in the near-infrared (Danchi et al. 2001) and its polarization perpendicular to the disk (Aitken et al. 1990; Yudin 1996). The radio-continuum emission shows an ionized bipolar structure arising above and below the neutral disk. The geometry of the ionized outflow has been modelled by considering a double-cone. The radio-continuum images are used to constrain the semi-opening angle of the double cone, θ_a (Fig. 3), to ~57°. This opening angle is consistent with that found with a semi-analytical model of ionized winds around sgB[e] stars (Kraus & Lamers 2003). The radio-continuum maps also show that for a given radius, the intensity is larger close to the neutral disk. In the isothermal case, this fact indicates that the electron density decreases from the edge of the ionized disk toward the revolution axis (0 < θ₀ < ∞).

Next we have constrained the physical structure of MWC349A by considering the whole set of continuum measurements available to date for the spectral range covering from radio-frequencies to the mid-IR range, where there is a clear excess of emission with respect to that expected from the free-free emission of the wind (ν ~ 6000−10   000 GHz as shown in Fig. 5). For simplicity we have assumed an isothermal wind expanding at nearly constant velocity in order to explain its radio-continuum spectral index of 0.62 (Fig. 5).

The best fits to the radio-continuum images and SED (see Figs. 4 and 5, respectively) were obtained by considering the values of n_R, b_d, T_e, N_e⁽r = 1,θ = θ_a⁾, and θ₀ shown in Table 1. The value derived for θ₀ is consistent with that obtained by White & Becker (1985). Based on the SED, which is rising until at least 5.8 THz (Fig. 5), we derive an inner radius of the ionized wind of r_min = 3 AU or smaller, since larger r_min would imply optically thin free-free emission at ν < 5.8 THz. We note that the quality of the fit is not very sensitive to the consideration of a thin layer of ionized material located next to the neutral disk with different electron temperature, T_d (e.g. 9450 K obtained by the fit of the RRL profiles). Although the set of input parameters that best fit the SED is not unique, the fit to the radio-continuum maps breaks this degeneracy. Figure 4 shows that our radio-continuum predictions are in very good agreement with the maps reported by Tafoya et al. (2004). We also note that the constraints on the parameters that are affected by some uncertainties were improved by using the RRL observations as described below.

4.2. Constraining the orientation, kinematics, and physical structure of MWC349A on the basis of the H30α centroid map and profile

Martín-Pintado et al. (2011) explained the H30α RRL emission obtained with the Plateau de Bure Interferometer by considering a model with two kinematic components made up of an expanding outflow and an ionized Keplerian rotating disk (as represented in Fig. 3). Besides the explanation of the velocity peak separation of the maser spikes, the biggest achievement of the quoted model was to explain the behaviour of the H30α RRL centroid map. Its interpretation is not straightforward since β_mn < 0 holds for a wide range of electron density and electron temperature values (Strelnitski et al. 1996b) and, therefore, the stimulated amplification is not produced in a small region with very specific physical conditions, but in an extended region (see bottom panels of Fig. 7). Thus, in the case of RRL maser emission, the centroid map must be interpreted as the averaged position where the emission for every radial velocity is produced.

Martín-Pintado et al. (2011) have shown that the centroid map of the H30α (see Fig. 6) for radial velocities comprised between the two maser spikes (at − 16 km s^-1 and 32 km s^-1, respectively) is well fitted to a straight line. This behaviour suggests that the emission between these radial velocities mainly arises from the Keplerian ionized layer. This is also supported by the spatial distribution of the emission predicted by our model (see Fig. 7).

The physical reason for the emission from the maser spikes to originate mainly in the ionized Keplerian disk is because the electron densities in this region are close to that optimum for maser amplification of the H30α RRL. For example, an electron density of 4.0 × 10⁷ cm^-3 (the optimum density for an electron temperature of 10⁴ K; Strelnitski et al. 1996b) is reached at the border between the ionized and the neutral disk at a radius ~50 AU according to the results of our model.

On the other hand, the emission at larger radial velocities mainly arises from a region farther away from the star and located outside of the ionized circumstellar disk. Although the overall H30α centroid map clearly constrains the kinematic features of the outflow (Martín-Pintado et al. 2011), it leaves open some uncertainties due to the high number of input parameters required for the modelling. However, the RRL profiles lead to strong constraints on their values.

In particular, the H30α RRL profile shows two different spectral features (see Fig. 8), a narrow double-peak maser emission at radial velocities of ~− 16 and 32 km s^-1 arising from the ionized Keplerian disk, and two broader features at slightly larger radial velocities, ~− 45 and ~71 km s^-1, showing a clear asymmetry (hereafter wing humps). Although the two components could be confused, the model shows that while the features of the narrow maser spikes depend on the kinematics of the Keplerian layer, the wing humps depend on the characteristics of the outflow. As described in Sect. 4.2.1, these two maser components arising from two different regions with different kinematics were used to constrain some of the input parameters of the model.

Figure 8 shows that while the computed profile obtained by using the set of b_n coefficients provided by Walmsley (1990) reproduces the main features of the observational profile, it clearly fails to predict its intensity. Neither set of parameters (with an electron temperature larger than 4000 K) can explain the intensities of the observed spikes. On the contrary, the Storey & Hummer (1995) departure coefficients with the input parameters given in Table 1 are able to reproduce much better both the features of the profile and its peak intensity. Thus we have assumed the Storey & Hummer (1995) departure coefficients in the model in order to significantly improve the fit of the maser spike intensities. Unlike the profile, the modelled H30α centroid map remains similar for both set of coefficients.

We note that the observational profile of both the H30α and the rest of double-peaked lines (see Sect. 5) show significant differences between the intensities of the blue- and red-shifted peaks, while our model predicts roughly the same intensities for both spikes. However, as we can see in Fig. 8, the maser peak intensities are strongly dependent on the set of departure coefficients assumed by MORELI. On the other hand, since most of the emission of the maser spikes comes from very small regions (see Fig. 7) and it is strongly dependent on optical depth, clumps with slightly different electron temperatures or electron densities would produce significant changes in the peak intensities. Thus, differences in the electron temperature and electron density of clumps located in the blue- and red-shifted regions of MWC349A, together with the uncertainties in the departure coefficient values (see Sect. 2.1.3 and Fig. 2), could explain the observed asymmetries between the two maser peaks. Likewise, the presence of different clumps with slightly different physical conditions over time would also explain the rapid time variations of maser RRLs (Martín-Pintado et al. 1989b). This variability is seen in time scales as short as 30 days (Thum et al. 1992).

To explain the H30α maser peak-to-peak separation (see Table 4), we have considered a central mass for MWC349A of 38 M_⊙, unlike the ~60 M_⊙ obtained if we use the Walmsley (1990)b_n coefficients (Martín-Pintado et al. 2011). As thoroughly discussed by Kraus (2009), the location of MWC349A in the HR diagram is subjected to a high uncertainty, particularly because of poor determination of its luminosity (log T_eff = 4.37 ± 0.07, log (L/L_⊙) = 5.7 ± 1.0). This makes its location in the HR diagram consistent with the evolutionary tracks of a wide range of central masses of rotating and non-rotating massive stars (Schaller et al. 1992; Meynet et al. 1994). Nevertheless, if we take into account just the estimated central values of the effective temperature and luminosity (log T_eff = 4.37, ), we derive a mass for the star of ~38 M_⊙. This value is consistent with our results. These results seem to suggest again that the Storey & Hummer departure coefficients provide a better description of the stimulated amplification of the radiation in MWC349A.

For a 38 M_⊙ central mass, the best fit for both the H30α centroid map (Fig. 6) and the RRL profiles was obtained by considering the following input parameter values: i) outflow terminal and turbulent velocities of ~60 and 15 km s^-1, respectively; ii) Keplerian rotation velocity component for the the outflow (see Fig. 6) in the same sense as the ionized disk; iii) a Keplerian rotating ionized disk layer with an opening angle of 6.5° located next to the neutral disk and extended from 0.05 AU to 130 AU; iv) an electron temperature of 9450 K for the ionized Keplerian disk; and v) inclination angle of the disk axis with respect to the plane of the sky of 8°, where the front-side is tipped up. Next, we discuss how changes in the value of these parameters affect the fitting of the observational data.

4.2.1. Dependence of the modelled H30α centroid map and RRL profiles on the input parameters

4.2.1.1. Inclination angle, θ_i

The inclination angle of the disk axis is the key parameter that makes the emission asymmetric with respect to the east-west plane. Therefore, it is responsible for the north-south loops of the H30α centroid map shown in Fig. 6. The larger the inclination angle, the larger the asymmetry, with larger heights for the loop peaks at high radial velocities (~50 km s^-1), which correspond to the emission of the outflow. Thus, the inclination angle is clearly constrained between 4.5° and 15°, with the front side of the disk tipped up with respect to the line of sight. Thus, we explain that the red-shifted and blue-shifted loops (see Fig. 6) occur north and south of the disk respectively.

4.2.1.2. Terminal velocity, v₀

The model also shows that the radial velocity range where the maser emission of the outflow is produced depends on its terminal velocity. Thus, we can rule out terminal velocities larger than 100 km s^-1 since then the peaks of the wing humps would be located at larger radial velocities, clearly distinguished as high-velocity shoulders (see Fig. 9a). In contrast, a too low terminal velocity would imply that the wing hump emission would appear at nearly the same radial velocities as the maser spikes and, therefore, we would not be able to distinguish the contribution of the two different kinematic components (ionized disk and outflow). Thus, we conclude that the terminal velocity must be comprised between 40 and 100 km s^-1. This is also supported by the way in which the predicted H30α centroid map changes with this velocity. However, the fit to the linewidths of the H41α and H76α RRLs (cf. Sect. 5.1) provided the strongest constraint on the outflow terminal velocity. Thus, finally we considered a value of 60 km s^-1.

4.2.1.3. Turbulent velocity, v_tu

Turbulent motions in the ionized wind are needed in order to explain the lack of a double-peaked profile observed for the H66α line profile (see Sect. 5). We note that the two observed peaks would be due to stimulated emission from the ionized outflow and not from the disk, opposite to that found in mm RRLs. Thus, we need to consider a turbulent velocity. We conclude that the best choice for the turbulent velocity to explain the line profiles and the centroid map is v_tu ~ 15 km s^-1, which is comparable to the sound speed for 12 000 K. Larger turbulent velocities would increase the predicted linewidths.

4.2.1.4. Opening angle of the Keplerian ionized layer, θ_d

Given a semi-opening angle for the ionized wind (θ_a ~ 57° as constrained from radio-continuum maps, Sect. 4.2), the contribution of the emission arising within the ionized Keplerian ionized disk is determined by the opening angle of this disk, θ_d = θ_a − θ_w (see Sect. 3.2). The larger the opening angle of the disk, the lower the flux arising from the outflow and, therefore, the weaker the wing hump intensities. Thus, too low opening angles for the disk would significantly increase the wing hump intensities compared to observed levels. In addition, in this case the asymmetry between the red- and blue-shifted wing humps would also increase. We conclude that our model is very dependent on θ_d for angles smaller than 6°. In particular, we rule out that θ_d could be smaller than 4.5° (see Fig. 9b). Likewise, if one considers a larger angle (i.e. 15°), the wing hump intensities would decrease and the maser peak intensities would increase above the observed values (see Fig. 9b). Taking into account the aforementioned considerations, we conclude that the opening angle of the Keplerian ionized layer must be between 4.5° and 15°.

4.2.1.5. Electron temperatures

Regarding the electron temperature of the ionized outflow and of the ionized Keplerian disk, our model shows that the centroid map is less sensitive to changes in the electron temperature than the line profiles. We consider an electron temperature of 12 000 K for the outflow to adequately fit the radio-continuum data (Sect. 2) and the intensities of the wing humps. If we considered lower values, the maser effect would increase and, consequently, also the wing hump intensities as shown in Fig. 9c.

On the other hand, we have considered a different electron temperature for the Keplerian ionized disk, T_d, from that of the ionized outflow to increase the intensity of the maser spikes. Since the maser intensities depend strongly on the temperature (as shown in Fig. 9d), we constrain the electron temperature in the Keplerian ionized disk between 8500 and 11 000 K. We adopt T_d = 9450 K as the final value for the Keplerian ionized layer.

4.2.1.6. Size of the disk, r_d

In the equatorial region located beyond the circumstellar disk, the observed radio-continuum emission at cm wavelengths (Tafoya et al. 2004) must be due to the ionization of the gas by the dust-scattered ionizing radiation. This is because we observe emission in the equatorial region for radii larger than that of the beam for the 2, 3.6, 6, and 20 cm radio-continuum images (White & Becker 1985; Rodríguez & Bastian 1994; Cohen et al. 1985; Tafoya et al. 2004). From the beam sizes and the constrained semi-opening angle for the ionized gas, θ_a, we derive an upper limit to the radius of ~130 AU. This upper limit is consistent with the derived size of 75 AU from the 3.8 μm image (Danchi et al. 2001). On the other hand, since the Keplerian rotating ionized disk is formed by the photoevaporation of the neutral disk, we assume that it reaches at least the same radius as the neutral disk. Thus, we have assumed a radius for the Keplerian ionized disk of 130 AU. We have to stress that the results of the model are only very sensitive to the disk size for a radius smaller than 60 AU (50 mas) since the bulk of the mm-wavelengths RRL emission (for which the maser emission is mainly produced in the ionized disk) comes from a region within this radius (i.e. as happens in the H30α line seen in Fig. 7).

4.3. Fit to the line-integrated intensities

Using the structure and physical properties derived in previous sections, we carry out the modelling of the integrated intensity of RRLs to study the frequency range where the maser emission is taking place. For that purpose, we predict the Hnα integrated line intensities both for the non-LTE and LTE cases for the quantum numbers from n = 4 to n = 80 as shown in Fig. 10. For the non-LTE case, we show the predictions using both the Walmsley (1990) and the Storey & Hummer (1995) departure coefficients.

Firstly, we see that the two integrated line flux predictions obtained under the non-LTE assumption for the two sets of departure coefficients converge for Hnα lines at low frequencies (with n ≳ 41). The non-LTE predictions agree extremely well with the observations and show a constant offset with respect to the LTE prediction. This is a clear sign that RRL emission is stimulated at low-frequencies as claimed by Martín-Pintado et al. (1993) and as supported by the analysis of RRLs profiles described in Sect. 5.

For RRLs with quantum numbers smaller than 41, we can see that the predicted integrated intensities for the non-LTE case become much larger than the LTE prediction, especially for RRLs with n < 30. This is what one expects since the lower the electronic level, n, the larger the importance of the mechanisms causing the electron population inversion, as in agreement with the departure coefficients (Strelnitski et al. 1996b). Only the non-LTE case is able to explain the large increase of the observed integrated line intensities. The large dependence of the integrated line intensity on the quantum numbers for 21 ≤ n ≤ 41 is proof that the dominant process of amplification of the emission at such RRLs is maser.

Likewise, it is remarkable that the predictions obtained with the Storey & Hummer (1995) coefficients are different from those derived using the Walmsley (1990) coefficients. This allows us to discriminate which set of departure coefficients is the best to describe the maser amplification under the physical conditions of MWC349A. For RRLs with n ≥ 30, we clearly see that the Storey & Hummer (1995) departure coefficients are more appropiate since their integrated line intensities agree extremely well with the published data for 30 ≤ n < 41.

Thus, we conclude that the maser emission at such frequencies is consistent with the physical conditions and the kinematics that we have deduced based on the radio-continuum and RRLs profiles. This supports the proposed kinematic model and the estimated physical parameters. Nevertheless, the integrated line intensity predictions increase very strongly with the frequency compared to that observed for n < 30, as one expects from an exponential increase of the brightness temperature with the optical depth (Sect. 2.3). It indicates that saturation effects in the masers play an important role for n < 30. Assuming an equivalent maser beam solid angle, 4π/Ω_m, of 60 (Thum et al. 1994a), we have checked that the degree of saturation, J_ν,sat/J_ν, strongly increases from the H30α to the H26α, and especially for the H21α. Table 2 displays the radius, electron density, temperature of saturation, brightness temperature and degree of saturation derived from MORELI in those pixels where the upper limit for the degree of saturation is maximum. We clearly see that the degree of saturation for the H26α and, above all, for the H21α are large enough to significantly contribute to the decrease of the maser amplification of the radiation in at least some of the pixels.

Finally, the model shows that the non-LTE effects dominate the emission of RRLs down to n = 7, while RRLs with smaller quantum numbers seems to be formed under LTE conditions as discussed by Thum et al. (1998).

5. Analysis for RRL profiles toward MWC349A

In the following we will use the inferred kinematic properties of the ionized wind and disk for MWC349A to show that the predicted results of our model are also consistent with the observation of RRLs other than H30α. First of all, we will focus in the RRLs with quantum numbers n > 30. Since these RRLs trace outer regions characterized by gas with smaller electron density than those traced by the H30α RRL, their maser emission is expected to be weaker (Martín-Pintado et al. 1989a). Thus, the modelling of their profiles is less sensitive to the set of departure coefficients considered and to small variations of the input parameters. Furthermore we will also study the profiles of RRLs with n < 30, which trace the kinematics of the innermost regions of the disk.

The predicted profiles should explain the change from the double-peaked profile observed for the H35α (as for the H30α) into an apparently simple component for the H41α and also the large change in the peak intensities of these RRLs. This behaviour and the linewidths are perfectly predicted by the model as shown in Fig. 11 and in Tables 3 and 4. The observed change in the RRL profiles is expected since RRLs at lower frequencies have larger continuum optical depths (see Eq. (3)) and trace outer regions where the electron density in the Keplerian ionized layer is lower than the optimum values for maser amplification (i.e. 6.3 × 10⁶ cm^-3 for H40α and 1.4 × 10⁷ cm^-3 for H35α; Table 1 in Strelnitski et al. 1996b). For this reason, the H41α emission is mainly dominated by the outflow with little contribution from the disk and, therefore, its profile shows just a single peak.

We stress that under LTE conditions, we clearly underestimate the peak intensities for all RRLs by factors of two for the H76α and H66α lines, and by larger factors for the rest of the RRLs for which the stimulated emission is stronger. Thus, we rule out the possibility that even apparently single-peaked profiles like that of the H41α line could be formed in LTE.

5.1. Low-frequency radio recombination lines. The outflow

Quantitatively, we see that the predicted peak intensities under non-LTE conditions for apparently single-peaked RRL profiles agree very well with those reported in the literature (see Table 3), with the exception of the only spatially-resolved cm-wavelength RRL to date, the H66α. We predict its zero-intensity linewidth (Δv ~ 100 km s^-1), although we clearly underestimate the peak intensity reported by Loinard & Rodríguez (2010) even for non-LTE emission (i.e. peak intensity of 9.0 mJy for the LTE case versus the 16.3 mJy predicted in the non-LTE case, and the 24.9 mJy observed by Loinard & Rodríguez 2010). Thus, we conclude that even at this low frequency, the stimulated amplification of the radiation is still significant as was also claimed by Martín-Pintado et al. (1993). This is also supported by the available RRL profile with the lowest frequency, the H76α line. Its asymmetric profile is a clear indication that stimulated emission plays a key role for this RRL (Escalante et al. 1989). In this case, our non-LTE prediction matches quite well with the reported RRL profile, while the LTE prediction greatly underestimates its peak intensity and it is not able to explain the asymmetry in the profile (see Fig. 12). The larger intensity of the blue-shifted peak is due to the existence of a larger background continuum emission to be amplified by the foreground ionized material approaching us. Thus, we conclude that even at 2 cm, the consideration of non-thermal emission is essential to explain the observations.

5.2. High-frequency radio-recombination lines. The inner disk: a distorted Keplerian-rotation

For RRLs with quantum numbers n ≥ 30, our model accurately predicts both the peak intensities (see Table 4) and the maser peak-to-peak velocity separations, (v_red − v_blue)_obs, which does not depend so strongly on the physical conditions but on the kinematics of the region traced by the RRLs. Table 4 shows that the predicted change of (v_red − v_blue) from H30α to H39α is consistent with the observations when one assumes the Storey & Hummer departure coefficients. Since the maser peaks mainly trace the Keplerian-rotating ionized disk, we conclude that the derived central mass for the source is consistent with the reported (v_red − v_blue)_obs.

The kinematics of the innermost region of the ionized circumstellar disk can only be studied by using the RRLs at higher frequencies since the continuum becomes optically thinner with increasing frequency, tracing the innermost parts. At present there are velocity-resolved observations for three transitions with frequencies higher than that of the H30α. These are the H27α, H26α, and H21α. Table 4 shows the predicted peak intensity for the three maser lines are clearly overestimated in the unsaturated case, especially for the H21α which is roughly 10¹⁴ times more intense than that observed (Thum et al. 1994b). Thus, as discussed, saturation effects could prevent the maser intensity from increasing with an exponential dependence.

Table 4 shows that our modelling is not able to explain the observed maser peak-to-peak velocity separations, (v_red − v_blue). The predicted peak-to-peak velocity separation increases with the frequency of the RRL, as expected from the assumption that the disk follows a Keplerian rotation, while observations show that it remains approximately constant from the H30α to the H21α line.

This fact could be explained if the submm RRL maser emission occured in the same regions as for the H30α. This would be the case if there was a lack of ionized material inside the region where the bulk of the H30α emission arises or an inner rarefication of the disk as claimed by Kraus et al. (2000). Nevertheless, in such a case, it would not be possible to explain the radio-continuum spectral index of 0.62 up to the frequency of the H21α line. On the other hand, saturation effects of interaction between masing transitions could explain that the maser lines are produced in the same region (Strelnitski et al. 1996b). Hengel & Kegel (2000) studied the effects of the radiation field on the electronic level population by using a radiative transfer model for a spherically symmetric source. They derived the ranges of electron density and quantum numbers where maser emission arises. While the optimum electron density values for unsaturated maser amplification are very sensitive to the quantum numbers of the RRLs (Fig. 6 in Strelnitski et al. 1996b), the saturation of the maser lines counteract this behaviour such that the optimum electron density of the bulk of maser lines are contained in a small range of electron densities and quantum numbers (Fig. 4 in Hengel & Kegel 2000). Thus, saturated maser lines would arise from regions with similar electron densities, i.e. from regions of similar radial distance and thus similar radial velocities. In such a case, the peak-to-peak velocity separation for saturated lines would not increase with the RRL frequency but would remain constant as observed for n < 30. This provides good evidence that RRLs with n < 30 are saturated and, therefore, our model is not able to reproduce these lines.

A third possibility to explain the constant peak-to-peak velocity separation from the H30α to submm RRLs results from the inner regions of the disk not following the Keplerian rotation of the outer parts. The kinematic discontinuity could likely be because of disturbances in the kinematics of the ionized circumstellar disk in the region where the outflow is launched.

5.3. Fit of the Hnβ radio-recombination lines

Finally, we have used the Hnβ transitions to further constrain the physical parameters derived from our modelling. The predicted Hnβ line intensities for LTE and non-LTE assumptions are shown in Table 5. The physical structure and the departure coefficients assumed by MORELI predict a significant amplification for Hnβ lines, but much weaker than for the Hnα lines due to their lower | β_mn |. However, contrary to what was found in the Hnα transitions, the Hnβ transitions can be best explained by assuming LTE emission. The difference between the observed and predicted LTE intensities are within ~4 times the rms noise of the data. The good fit obtained for the observed Hnβ transitions could represent another proof that supports the quantitative parameter values obtained by our modelling, but only if the Hnβ lines are not amplified.

However, the amplification predicted by our model for the Hnβ lines is unclear. A possible explanation could be related to the fact that RRLs sample the physical conditions where the continuum emission is optically thin. Thus, pairs of Hnα and Hn′β transitions with similar frequencies arise from regions with similar electron densities. For a certain electron density, the b_n coefficients might represent roughly the electronic population for a range of quantum numbers (Hnα lines), but not for higher quantum numbers (Hn′β lines) for which the electronic levels are closer to thermalization. As an example, we can consider the H30α and H38β pair at similar frequencies mainly sample regions with electron densities of ~10⁷ cm^-3. However, we can see in Fig. 13 that the β_mn coefficients for H30α and H38β significantly differ. In particular, the β_mn coefficients for H38β are close to zero while for H30α we find β_mn = −4. As we can also see in Fig. 2, the numerical error in the estimation of the β_mn coefficients increases when β_mn gets close to thermalization. If the electronic levels actually thermalize faster when there is high electron density than when the model predictions are considered, then the β_mn coefficients for the Hnβ lines would be closer to zero than our model assumes and the amplification would be negligible. In this case the H38β line would appear to be emitted under LTE conditions as observed. Thus, the uncertainty for the β_mn when the electronic population approaches thermalization might explain why the Hnβ are not amplified as expected from our model.

6. Conclusions

This paper describes the MORELI 3D radiative transfer model to predict RRL emission under LTE and non-LTE conditions in a variety of geometries, physical structures, and kinematics. Previous versions of this model were used in order to study different UC-HII regions: CRL618 (Martín-Pintado et al. 1988), MWC349A (Martín-Pintado et al. 1989a, 2011), Cepheus A HW2 (Jiménez-Serra et al. 2011), and MonR2-IRS2 (Jiménez-Serra et al. 2013). In this paper, we used MORELI to explain the continuum and RRL data observed in MWC349A in a coherent way. This is the first complete modelling of MWC349A, which not only constrains the model parameters from a subset of RRL data but also for the whole set of observed RRL. Thus, our new results for the geometry, struture, and kinematics of MWC349A provide the best description of the characteristics of the ionized wind and outflow in this source, and of the central mass of the source.

Mainly on the basis of the explanation of the complex behaviour of the H30α centroid map and its RRL profile, we have strongly constrained the geometry, disk inclination, physical conditions, and kinematics of MWC349A. These results were also derived with a systematic study of the whole set of radio-continuum and RRL observations. We found that the set of departure coefficients b_n provided by Storey & Hummer (1995) describes much better the intensities of mm and submm masers than those given by Walmsley (1990) for the physical structure and geometry derived for MWC349A. With these coefficients we are able to explain for the first time intensities of the maser spikes as well as the peak-to-peak separation between the spikes of mm RRLs.

Thus, our modelling is supported by the fact that they explain the bulk of the different observational data available (SED, radio-continuum emission maps, RRL profiles, and velocity-integrated line intensities for RRLs) at very different wavelengths (from mm to cm wavelengths). In addition, our model provides a general view of the dominant processes involved in the amplification of RRLs for the different frequency ranges. While maser amplification seems to dominate for Hnα RRLs with n < 41, stimulated amplification is still important even up to frequencies as low as that of the H76α, 15 GHz. However, our model does not reproduce the peak intensities, peak-to-peak velocity separation, and velocity-integrated line intensities for submm and infrared RRLs with n < 30. Thus, one of the open issues to be addressed in the future is taking into account the saturation effects. High spatial resolution and velocity-resolved studies by using high-frequency RRLs (with n < 30) are also required to establish the kinematics of inner regions.

The main difficulty found for our non-LTE modelling is to predict the observed peak intensities for the Hnβ line. While we expected the Hnβ lines to be out of LTE and with significant amplification, their observational peak intensities agree very well with the predictions for LTE emission. High spatial resolution observations of Hnβ should be perfomed in order to explain this puzzling behaviour. However, precise calculations of the b_n coefficients for the geometry and physical conditions in MWC349A might explain these discrepancies.

In summary, we have used MORELI to better constrain the characteristics of both ionized components on MWC349A, disk and outflow, as proved by the large amount of observational data explained by our modelling. Further progress depends on observations which combine high angular and high velocity resolution, like those possible with ALMA.

---

Acknowledgments

This work has been partially funded by MICINN grants AYA2010-21697-C05-01, FIS2012-39162-C06-01, and Astro-Madrid (CAM S2009/ESP-1496). Alejandro Báez-Rubio acknowledges support from grant JAE predoct (2009), CSIC, Spain. We would also like to thank Josefina Torres for providing Fig. 3 and L.F. Rodríguez for providing the observational radio-continuum images of MWC349A (Fig. 4) and the H76α RRL profile (Fig. 12). Finally, we are also grateful to the anonymous referee and Malcom Wamsley for very valuable comments.

References

Appendix A: Gaunt factors

In this appendix we show the equations used by MORELI to estimate the Gaunt factors. In the literature we find different approximations for the free-free Gaunt factor for different frequency and temperature regimes. MORELI uses the analytic approximation given by Gronenschild & Mewe (1978), which reproduces the computational results of Karzas & Latter (1961) with an accuracy of 10% in the frequency and electron temperature ranges where 10^-2 < hν/⁽kT_e⁾ < 10³: (A.1)where K₀ is the modified Bessel function of the second kind, Z the effective charge of the ionized gas, and ν₀ the hydrogen ionization frequency. This approximation is equivalent to that of Leitherer & Robert (1991) in the Rayleigh-Jeans regime (hν ≪ kT_e) provided that the frequency exceeds the plasma resonance frequency, ν ≫ ν_p (where ν_p = N_ee²/⁽πm_e^)1/2). Thus, the Gronenschild & Mewe expression can be used in a broader range of frequencies than the previous one.

In the absence of dust, free-free continuum emission is the dominant process for wavelengths λ ≳ 10 μm, while at NIR and optical wavelengths for which g_bf,ν ≃ 1, the bound-free processes become significant. This is shown in Fig. A.1, where we plot the relative contribution of the bound-free and free-free Gaunt factors. The bound-free Gaunt factor is calculated by using the expression given by Brussaard & Van de Hulst (1962)(A.2)where g_n(ν) is approximated to 1 for all the frequencies. This approximation yields results with an accuracy of ~10−20%.

Appendix B: Radio-recombination line profiles

The RRL profiles of UC-HII regions are due to the combination of different broadening mechanisms whose convolution gives the observed profile. In Sect. 3.2 we have shown the large-scale motions that give as a result the broadening of the profile. However, even in every cell of the model (Sect. 2), there is a profile broadening caused by the microscopic motions of the ionized gas, such as turbulent and thermal motions.

Regarding the thermal motions, one can consider the approximation that the ionized gas is well described by the Maxwell velocity distribution. It subsequently produces a Gaussian profile because of the Doppler shifting. Likewise, the turbulent motions can be also described by a Gaussian distribution governed by a parameter known as turbulent velocity, v_tu. Thus, the velocity half-width at half-maximum height, Δv_g, of the Gaussian profile due to both contributions can be estimated using the expression (Mezger & Hoglund 1967) (B.1)Additionally, there is a contribution due to the broadening of the electronic levels in regions with high electron density. This is known as pressure broadening and it produces a Lorentzian profile. This broadening can be specified by the half-width at half-maximum height, Δv_l, which has been computed numerically using analytic expressions that approximate the results derived from classical and semi-classical theory (Gee et al. 1976) for the different quantum number ranges of the RRLs (see Table B.1).

The convolution of the aforementioned Gaussian and Lorentzian distributions results in a Voigt profile, which can be expressed as a function of Δν_l and Δν_g (where ) by (B.2)This integral is numerically computed in MORELI following the formulation of Kielkopf (1973) to describe the RRL emission in every cell of the source. Besides the intrinsic line profile, the projection of the large-scale motions (modelled in Sect. 3.2) along the line-of-sight gives rise to a dynamical broadening of the different RRL profiles.

Appendix C: Effective radius for an optically thin wind

The total free-free continuum flux for a spherical outflow in the circular region in the plane of the sky of radius R_eff is given by (C.1)where ρ and R_eff are in units of the radius of the central hole without ionized gas, r_min, with and R_eff > 1.

The optical depth along the line-of-sight located at a radius ρ from the center of the source, τ_ν(ρ), for the spherical wind isotropically expanding at constant velocity, is given by (C.2)By using the optical depths given above, the integration of Eq. (C.2) in a circumference of radius R_eff for the optically thin case, 1−e^{− τ_ν(ρ)} ≈ τ_ν(ρ), results as (C.3)If we define R_eff as the radius containing a given percentage, p, of the total continuum flux, it results from Eq. (C.3) the following effective radius in units of r_min: (C.4)

Appendix D: Kinematics

Appendix D.1: Disk kinematics

In this appendix, we derive Eq. (18) (Sect. 3.2.2) describing the velocity component along the z axis for a Keplerian rotating ionized disk. We assume that the only disk velocity component is the rotation around the revolution axis specified by the x_d axis. Thus, in order to describe its velocity, the cylindrical polar system of coordinates unit vectors (e_ρ,e_ϕ,e_xd) and coordinates (ρ,ϕ,x_d) is the most suitable. In this case the velocity is given by (D.1)Describing Eq. (D.1) in a Cartesian coordinate system with the same polar axis and unit vector (y_d,z_d,x_d)¹, and assuming a Keplerian rotation velocity, Eq. (D.1) results as (D.2)where .

Since in the general case, the disk could be tilted with respect to the line-of-sight with an angle θ_i (see Fig. 3), we apply a rotation of coordinates around the axis y_d to lead to the Cartesian coordinate system with unit vectors (e_y,e_z,e_x) orientated as described in Sect. 2. Such a coordinate transformation is given by the equations (D.3)In this coordinate system, Eq. (D.2) is (D.4)where , y_d = y, and z_d = cos(θ_i)z − sin(θ_i)x.

Thus, the velocity component along the line-of-sight (the z axis) yields (D.5)For a nearly edge-on disk (θ_i ≪ 1), Eq. (D.5) can be approximated as (D.6)

Appendix D.2: Outflow acceleration

In this appendix, we derive Eq. (19) (Sect. 3.2.2) describing the velocity component along the z axis for an outflow expanding radially away from the central star reaching a terminal velocity of v₀ at a radius r_a. The total radial velocity of the ionized gas at r < r_a is given by . Thus, its projection along the z axis of the Cartesian coordinate system with the unit vector (y_d,z_d,x_d) described in (D.1) is (D.7)where θ is the polar angle and ϕ the azimuthal angle.

Thus, we obtain the velocity projected to the z-axis (D.8)

All Tables

All Figures