© ESO 2021

1 Introduction

As pointed out already by Bohr (1920), in the limit of large quantum numbers n, all atomic spectra appear ‘hydrogenic’ with the ionisation energy of a single electron far away from the ionic core well described by the classical Rydberg formula, E_n ∝n⁻². Such highly excited atoms occur frequently in different astronomical objects where conditions are favourable; for instance, at low densities and in intense photo-ionising fields. The theory of Rydberg transitions in the cosmic nebular environment was first developed by Menzel and collaborators in a series of papers starting in 1937 (e.g. Menzel 1937). The possibility to detect such cosmic lines at radio wavelengths was first discussed by van de Hulst (1945) and clarified quantitatively by Kardashev (1959), specifically in the tenuous nebular medium which is exposed to intense UV light, where highly excited Rydberg states are populated by radiative recombination of the parent ion and a free electron; hence the name, (radio) ‘recombination’ lines.

The first detections of cosmic radio recombination lines from hydrogen were made by Dravskikh & Dravskikh (1964) and Höglund & Mezger (1965) towards H II regions. Subsequently, recombination lines from helium were detected towards H II regions, and carbon lines were detected from a cooler medium surrounding H II regions (Lilley et al. 1966; Pankonin et al. 1977). Chaisson et al. (1972) found a weak feature on the blueshifted side of the C158α line, and identified it as due to one or a combination of the heavy element species Mg, Si, S, and Fe. Safe identification of the carrier is often hampered by the fact that the line separations – only due to the differences in reduced mass – of elements heavier than carbon are smaller than the line widths. However, in the cold ρ Oph dark cloud, Chaisson (1975) managed to identify a sulphur recombination line. Recombination lines have now been detected also towards other ionising sources like planetary nebulae (Roelfsema et al. 1991; Bachiller et al. 1992; Sánchez Contreras et al. 2017), luminous blue variables like η Car (Cox et al. 1995), and young stellar objects (Martin-Pintado et al. 1989; Thum et al. 2013). Comprehensive reviews of recombination lines are given in Brown et al. (1978) and Gordon & Sorochenko (2002). Also the Sun has been detected in Rydberg transitions at submillimetre (H and possibly Mg Clark et al. 2000a,b) and mid-infrared (Mg and Al Murcray et al. 1981; Brault & Noyes 1983; Chang & Noyes 1983) wavelengths. Here the conditions are different, and a different approach is required to describe the excitation (Carlsson & Rutten 1992); therefore, using the nomenclature of Rydberg transition lines (RTLs) rather than recombination lines is more appropriate.

Here wepresent the detection of spectral lines at about 232 and 354 GHz towards the moderate-temperature, low-mass star HD 101584 using the Atacama Large Millimeter/submillimeter Array (ALMA). We argue that these lines are RTLs of neutral atoms of elements heavier than carbon. A simple model in strict local thermal equilibrium (LTE) is used to estimate the viability of this proposal, and to estimate the physical characteristics of the line-emitting region and identify the most likely carrier (or carriers) of the lines.

2 HD 101584

HD 101584 is a binary system with a significant infrared excess. For the primary star, we adopt a spectral type of A6Ia and effective temperature of 8500 K; see Kipper (2005) who determined the temperature from spectral synthesis of high-resolution optical spectra. It has been difficult to establish whether its current state represents evolution beyond the asymptotic giant branch (AGB) or the red giant branch (RGB). It was originally classified as a post-AGB star (Parthasarathy & Pottasch 1986; Bakker et al. 1996), but Olofsson et al. (2019) classified it as a post-RGB star, the consequence of premature termination of the RGB caused by common-envelope (CE) evolution that ended before the two stars merged.

Olofsson et al. (2019) studied the circumstellar environment of HD 101584 in detail using ALMA. This object is morphologically and kinematically complex, consisting of a bipolar high-velocity outflow (seen almost pole-on) surrounded by an hourglass structure, an equatorial density enhancement (seen almost face-on), and a central compact source (whose geometry and orientation are unknown). A total of 12 different molecular species (not counting isotopologues) have been identified in the circumstellar medium. Here we recapitulate our present knowledge of the near environment of HD 101584. For more information on this source, see Olofsson et al. (2019) and references therein.

Many of the molecular line emissions observed using ALMA show a compact component that coincides in space with the position of HD 101584, the central compact source (CCS) in the nomenclature of Olofsson et al. (2019); detected molecular species are CO, SiO, SiS, CS, SO, SO₂, OCS, and H₂S. The emission from this component is partially resolved, and its estimated mass is ≈0.03 [D/1 kpc]² M_⊙, which corresponds to ≈7% of the total circumstellar gas mass. The angular size varies with the observed molecular transition, but lies in the range 0.′′ 1–0.′′15 (FWHM of Gaussian fit to the brightness distribution deconvolved with the synthesised beam), that is, within ≈75 [D/1 kpc] au of HD 101584. The morphology of this component has not been determined, except for the innermost structure observed with the VLTI/PIONIER instrument. Kluska et al. (2020) detected an almost circular structure interpreted to be the sublimation rim of a circumbinary dust disc seen close to face on. The radius of this ring is about 2 [D/1 kpc] au. The lines discussed in this paper most likely come from a region within a radius of about 5 [D/1 kpc] au from the binary. The binary separation is estimated to be about 0.5 au.

The systemic velocity of HD 101584 is estimated to be 41.5 ± 0.15 km s⁻¹ (in LSR scale). We have derived the average velocity from Gaussian fits to six ALMA lines that show strong emission from the CCS: CO(2–1), C¹⁸ O(2–1), SO(5₅−4₄), SO(5₆−4₅), SO₂ (16_3,13−16_2,14), and H₂S(2₂₀−2₁₁); see Olofsson et al. (2019). The systemic velocity is crucial for deriving the rest frequencies of the lines discussed in this paper.

The parallax of HD 101584 in the second Gaia data release suggests a distance as large as 2 kpc (Gaia Collaboration 2016, 2018). The third Gaia data release maintains this distance, but the excess noise required in the fitting process has increased (Gaia Collaboration 2021). However, as argued by Olofsson et al. (2019) there are reasons to suspect that this distance is not correct, the primary one being that the estimated angular distance between HD 101584 and its companion is comparable to the parallax. Consequently, we follow their adopted procedure of presenting all results for a nominal distance of D = 1.0 kpc. From a fit to the spectral energy distribution, Olofsson et al. (2019) derived a luminosity and radius of the star of 1600 L_⊙ and 19 R_⊙, respectively,at 1 kpc, based on the assumption that the radiation of the infrared excess is not isotropic.

3 Observations

3.1 ALMA data

The ALMA data were obtained during cycles 1 (May 2014, TA1), 3 (October 2015, TA2; September 2016, TA3), and 6 (November 2018, TA4) with 35 to 46 antennas of the 12 m main array in one frequency setting in each of Band 6 (TA1–TA3; 232 GHz) and Band 7 (TA4; 354 GHz). In both settings, the data set contains four 1.875 GHz spectral windows with 1920 channels each, except for one spectral window in the cycle 3 data that has 960 channels covering 1.875 GHz. The final velocity resolutions are 3.0 and 1.0 km s⁻¹ for the Band 6 and 7 data, respectively. The baselines range from 13 to 16 196 m for band 6 and 15 to 1398 m for Band 7. This means a highest angular resolution of about 0.′′03 and 0.′′ 2 for the two bands, respectively. Bandpass calibration was performed on J1107-4449 (TA1, TA2, TA3) and J1037-2934 and J0635-7516 (TA4), and gain calibration on J1131-5818 (TA1), J1132-5606 (TA2, TA3, TA4). Flux calibration was done using Ceres and Titan (TA1), J1131-5818 (TA2), J1150-5416 (TA3), and J0635-7516 and J1037-2934 (TA4). Based on the calibrator fluxes, we estimate the absolute flux calibration to be accurate to within 5%.

The data were reduced using different versions of the Common Astronomy Software Applications package (CASA; McMullin et al. 2007). After corrections for the time and frequency dependence of the system temperatures, and rapid atmospheric variations at each antenna using water vapour radiometer data, bandpass and gain calibration were performed. For the Band 6 setting, data obtained in three different configurations were combined. Subsequently, self-calibration was performed on the strong continuum for the Band 6 data. Imaging was carried out using the CASA tclean algorithm after a continuum subtraction was performed on the emission line data. The final line images were created using Briggs robust weighting (0.5). This resulted in beam sizes of 0.′′ 028 × 0.′′026 (38°) and 0.′′ 22 × 0.′′17 (–48°) for the Bands 6 and 7 data, respectively. Typical rms noises per channel are ≈0.6 and ≈1.3 mJy beam⁻¹ for the Bands 6 and 7 data at 3 and 1 km s⁻¹ resolution, respectively.

3.2 Archival optical spectra

The archive of the European Southern Observatory (ESO) contains some excellent high-resolution spectra of HD 101584. These spectra were obtained with the ESO VLT/UVES instrument through the observing program 266.D-5655(A) and were retrieved by us as pipeline-processed data products. All 19 exposures were made on 2003-02-04 between 09:00 and 09:30 UT (JD 2 452 674.9). The spectra in the visual and red region (wavelengths 4727–10 426 Å) have a resolving power of R = 74 450, corresponding to a velocity resolution of 4.0 km s⁻¹. The blue spectra (3760–4990 Å) have R = 65 030 or a resolution of 4.6 km s⁻¹. We performed a detailed analysis of small wavelength intervals to search for emission lines that might be associated with the region that emits the observed mm-wave emission lines.

Fig. 1 X30α (upper spectrum; scaled by a factor of three and biased by 25 mJy, 3 km s−1 resolution) and X26α (lower spectrum; 1 km s−1 resolution) lines observed towards HD 101584. Continuum emissions have been subtracted.

4 Line identification

We detect two lines, at rest frequencies of 232.0245 and 353.8100 GHz (adopting a systemic velocity of 41.5 km s⁻¹) in our ALMA observations towards HD 101584; see Fig 1. Both have Gaussian line shapes with full widths at half maximum (FWHMs) of about 9 km s⁻¹. Both line emissions are spatially unresolved and centred on the continuum peak, which is believed to be the position of the star. The best upper limit to the line-emitting source size is obtained for the 232.0245 GHz line, with a FWHM ≲ 0.′′02 assuming a Gaussian source (after deconvolving with the beam).

The source is modestly rich in molecular lines with a line density of less than three lines per GHz (Olofsson et al. 2019), but we have found no identification in terms of relevant molecular lines¹. Furthermore, the line widths differ from those of the molecular lines obtained within an aperture of about 0.′′ 2, Sect. 6.4. On the other hand, these rest frequencies are consistent with an identification in terms of the X30α and X26α RTLs from an element X heavier than carbon. As an example, the S30α and S26α lines lie at 232.0233 and 353.8093 GHz. In Sect. 5 we present a simple model that makes the identification in terms of atomic RTLs viable. Nearby Xnβ lines lie outside the frequency coverage of our observations.

Direct identification of the carrier(s) is more difficult because the RTLs from heavy elements are close together in rest frequency, for instance, the ²³Na and ⁵⁶Fe 30α and 26α lines lie at 232.0217 and 232.0250 GHz and 353.8069 and 353.8112 GHz, respectively. That is, they lie on each side of our determined rest frequencies, and well within the observed line widths that correspond to ≈7 and 11 MHz for the X30α and X26α lines, respectively. Therefore, theoretical modelling is required to make further progress.

A useful formula for calculating rest frequencies of RTLs (ignoring the fine structure of hydrogen and the electronic structure of the core for non-hydrogenic atoms), n + Δn →n, is given by the Rydberg formula of Towle et al. (1996), $$\begin{array}{lll}{\nu }_{n+\Delta n,n}\hfill & =\hfill & 3.289842\times {10}^6(1-5.48593\times {10}^{-4}/m_{\text{a}})\hfill \\ \hfill & \hfill & \times \left[\frac{1}{n^2}-\frac{1}{{(n+\Delta n)}^2}\right][\text{GHz}],\hfill \end{array}$$(1)

where m_a is the atomic mass of the species in atomic mass units (in the usual nomenclature Δn = 1 is an α line, Δn = 2 a β line, etc.).

The continuum flux from the line-emitting source is difficult to estimate due to contamination with dust continuum emission. The Band 6 data have the highest resolution and using different uv-tapering we estimate flux densities of 7.5, 10.1, and 13.2 mJy for beam sizes of about 0.′′025, 0.′′ 083, and 0.′′ 15, respectively. This gives a flux density that increases only slowly with beam area, which suggests that the 7.5 mJy within an estimated source size of 0.′′ 011 (the FWHM of a Gaussian source after deconvolving with the beam) is relatively free from dust emission at this frequency². The Band 7 continuum source size is 0.′′19, as measured by an 0.′′23 beam, and the flux density is 61 mJy. Even if the dust contribution increases only slowly with beam area, the contribution from dust continuum is likely to be significant in this case.

The line characteristics are summarised in Table 1. The line strengths are 7.3 and 25 mJy and the line-to-free-free-continuum flux density ratios are >1.0 and >0.4 for the 30α and 26α lines, respectively. Upper limits to the corresponding H, He, and C lines are 1 and 2 mJy in the two bands, respectively. That is, for the 26α line the X/H line intensity ratio is ≳10. The lines are perfectly Gaussian, within the noise levels, and there is no evidence of low-intensity, broad line wings.

5 A simple model in strict LTE

The atomic lines reported here in the mm-wave spectrum of HD 101584 arise in the immediate environment of the binary system, within a fewastronomical units, as shown by the observations. A complete analysis of the data requires a non-LTE model that takes into account chemistry, ionisation, excitation, and radiative transfer in a source of arbitrary geometry. However, this is far beyond the scope of this paper. Instead, we have devised a model in strict LTE that will explain the presence of heavy-element RTLs and a free-free continuum, at their observed intensities, and the simultaneous absence of the corresponding C and H RTLs. With such a model it is possible to construct models of the RTL emission, from an assumed geometry of the gas, with the smallest number of free parameters. The model can also be used to easily make predictions for future observations. Because of the complexity of the problem, it is difficult to argue in detail whether or not such an approach is justified, but here we present some arguments in favour of the LTE model.

In our opinion, the primary star with its modest luminosity (≈1600 [D/1 kpc]² L_⊙) and moderate temperature (≈8500 K) is not able to produce an ionised nebula that can explain the observed RTL and continuum flux densities and, at the same time, the absence of the corresponding C and H RTLs. First, the high density that is required to explain the observed line and continuum flux densities (see below) will lead to efficient absorption of the ionising radiation already close to the star (the density is much higher here than in the traditional nebular case). The alternative, that the RTL region is located at a distance from the star with an absorption-free cavity in between, appears less likely. Second, such a nebula, produced by a star with the characteristics of HD 101584, would produce H and C RTLs that are stronger than the corresponding heavy-element lines, in stark contrast to the observed result.

The question remains as to whether or not strict LTE is well justified. In order for the ionisation, chemical abundances, internal excitation, continuum brightness, and line intensities to be well characterised by equilibrium conditions, the rates of all microscopic processes must be in detailed balance at the same value of the temperature. The densities of the best-fitting LTE models (see below), n_H ≈ 10¹¹–10¹² cm⁻³, are higher by orders of magnitude than typical densities in the interstellar medium where non-equilibrium prevails. On the other hand, this density is still much lower than the density in a typical stellar photosphere at the same temperature, where LTE is a useful starting approximation, but where subtle non-LTE effects can nevertheless be found in ionisation and line formation of some elements. Simple estimates of chemical, photodissociation, and ionisation timescales, as well as comparisons of collisional and radiative strengths, suggest that the LTE assumption is, if not perfect, at least acceptable. In addition, the line flux density ratio, 26α/30α = 3.4 ± 0.4, is not too far from the blackbody ratio (in the Rayleigh-Jeans limit) for emission at 232 and 354 GHz, namely 2.3, and we therefore suspect that the lines are formed in thermal equilibrium. The slightly higher value for the observed flux density ratio than the blackbody ratio can be attributed to an optical depth effect, that is, the lines are at least partially optically thin and the optical depth of the 26α transition is higher than that of the 30α transition.

Consequently, we assume that the abundances of molecular species, the ionisation balance, and the level populations are all governed by the processes that tend to drive the gas toward the equilibrium state at the kinetic temperature prevailing in the near circumstellar environment, that is, we assume that strict local LTE applies. This assumption is also the reason behind our choice to deliberately refer to the spectroscopic features as RTLs rather than radio recombination lines throughout this paper in order to steer the reader away from any prejudice that the lines arise mainly through the process of radiative recombination or that they should behave like the well-known features in nebular spectra. Further, the equilibrium temperature must be low enough that ionisation of hydrogen, helium, and carbon does not occur. We compare the observational constraints (Table 1³) to a simple model of a uniform, isothermal parcel of gas in strict LTE. This parcel of gas is characterised by two parameters: a number density of hydrogen in all forms, n_H, and a temperature, T. To this should be added the parameters defining the geometry and kinematics. Once the true geometry and kinematics of the source have been established, a more sophisticated modelling can be attempted.

5.1 Geometry

We have chosen a geometry in light of the VLTI/PIONIER results of Kluska et al. (2020), who infer a face-on (as seen from Earth), circumbinary dust disc with a sublimation rim at about 2 [D/1 kpc] au from the central object. The region in which the RTLs are formed could therefore consist of dust-free gas surrounding the binary system and extending along the upper and lower surface layer of this dust disc (the estimated source size is larger than the diameter of the dust rim, see below).

Considering this, we have chosen a simple geometry in terms of a cylinder seen along its axis, where R is the radius of the source in the plane of the sky and Δz is the depthof the cylinder along the line of sight. The apparent angular diameter of the emitting region is θ = 2R∕D, and we specify the geometry of the emitting region by the parameter ρ = Δz∕R. The case ρ = 1 corresponds to a cylinder ofradius R and height Δz = R viewed face-on at a distance D. When ρ ≪ 1, the modelled flux refers to that of a very thin disc, while the ρ = 2 case is close to that of a spheroid. We adopt ρ = 0.1, that is, a geometrically thin disc, as the most likely geometry for describing the RTL region. This also minimises the volume of the relatively warm gas required to explain the observations.

In strict LTE, it is possible to derive a lower limit to the size of the line-emitting region. The flux densities of the lines cannot exceed that of a blackbody at a given temperature. For a temperature of 3000 K (it is shown below that this is a reasonable temperaturein this context), both lines give about 10 mas as the lower limit to the source size. This agrees very well with the measured Gaussian size of the continuum source at 232 GHz, namely 11 mas. This is our best estimate for the source size because the 232 GHz data have a much higher angular resolution than those at 354 GHz, and the continuum data have a much higher S/N (≈70) than the line data (≈13). We make the reasonable assumption that the line and continuum emissions come from the same region. As shown by van Hoof (2000) the relation between the true source size and the size of a Gaussian fit to its brightness distribution depends on the geometry of the source; for instance, a uniform brightness disc is about 1.5 times larger than the FWHM of the Gaussian fit to its brightness distribution. We will nevertheless adopt 11 mas as the size of the RTL region, that is, the assumed thin disc extends to R = 5.5 [D/1 kpc] au. The reason is twofold: first, and most importantly, the resulting density and temperature depend only weakly on the source size (see Sect. 6.2.3 and Appendix D), and secondly, in our case the true source geometry is unknown.

Of course, other geometries are not excluded, keeping in mind that this is at the heart of the complex circumstellar morphology as shown by our ALMA data (Olofsson et al. 2019). In all likelihood, an accretion disc, an associated jet, and the innermost part of an hourglass structure are located in the region (Olofsson et al. 2019). However, as shown in Appendix D, results such as mass, density, and temperature of the gas depend little on the assumed geometry.

5.2 Chemistry

We consider an ideal gas composed of hydrogen and the 14 other most abundant elements of comparable or smaller first ionisation potential. Helium is included, although it remains fully neutral at the temperatures of interest here. Two different sets of abundances are compared in the analysis, as summarised in Table A.1. The reference abundances are those of the solar photosphere (Y ; Set 1), and Set 2 consists of abundances depleted by a factor of d.

In addition to the elemental composition we take into account the formation of some simpler molecules, as this may have important consequences for the atomic abundances, the atomic ionisation balance, and for the dominant sources of continuous opacity in the gas, which depend on collisions between free electrons and ions and the principal neutral species H, H₂, and He. We therefore calculated the molecular abundances (at each grid point) for the species given in Table B.1⁴. We used the most recent values for partition functions and dissociation energies in the calculation of the equilibrium abundances.The molecular partition function is a sum over all low-lying electronic states and enough of their vibrational and rotational levels to yield a converged value at temperatures of interest. The full list of included molecules, their dissociation energies, and abundances in the reference model is presented in Table B.1.

It should be noted here that recombination lines from molecules can be ignored because the recombination of a molecular ion and an electron preferentially leads to a dissociation of the molecule.

5.3 Ionisation balance

All elements X_i are taken to exist in two ionisation states, the neutral atom X${}_i^0$ and first ion X${}_i^{+1}$, except hydrogen, which also forms a potentially important anion, H⁻. Each internal energy state in the gas is labelled with the indices i for the element, j = 0,±1 for the charge state, and k for the quantum state of energy E_ijk. In strict LTE, the population of state ijk, p_ijk, can be expressed as a number density given by the Boltzmann formula, $$\frac{p_{ijk}}{p_{ij}}=\frac{g_{ijk}\exp (-E_{ijk}/k_{\text{B}}T)}{Q_{ij}(T)},$$(2)

where $$Q_{ij}(T)={\sum }_k{g}_{ijk}\exp (-E_{ijk}/k_{\text{B}}T)$$(3)

is the partition function and p_ij the total density of the ion X${}_i^j$ summed over all its states, p_ij = ∑_kp_ijk. For each element i, the sum over all states of all its ions is related to the total number density of hydrogen in all forms, n_H, by $${\sum }_j{p}_{ij}\equiv {\sum }_j{\sum }_k{p}_{ijk}=n_{\text{H}}{10}^{(Y_i-Y_{\text{H}})}.$$(4)

The ionisation balance is governed by the Saha equation $$\frac{p_{i1}{p}_e}{p_{i0}}=\frac{Q_{i1}{Q}_e}{Q_{i0}}\exp (-I_i/k_{\text{B}}T),$$(5)

where $$Q_e={\left(\frac{2\pi m_e{k}_{\text{B}}T}{h^2}\right)}^{3/2}$$(6)

is the partition function of the continuum of free-electron states (translational energies), I_i is the first ionisation potential of element i, k_B is Boltzmann’s constant, h is Planck’s constant, and m_e is the electron mass. Overall charge neutrality is assumed, so that the density of free electrons is given by $$p_e={\sum }_i{p}_{i1}-p_{\text{H},-1}.$$(7)

Given a set of element abundances Y_i and the two parameters T and n_H, Eqs. (2)–(7) are solved simultaneously by iteration. The result is a complete set of densities, p_ijk. We computed the partition functions of neutral atoms and first ions by direct summation over tables of energy levels derived from the NIST Atomic Spectra Database⁵.

5.4 Radiative transfer

Once the populations have been determined, the continuum opacity can be evaluated. At radio, mm, submm, and far-infrared wavelengths, the important processes in the weakly ionised gas are electron–ion bremsstrahlung (free-free), electron–neutral bremsstrahlung, and Thomson scattering. The electron–ion free-free absorption coefficients are computed using the accurate Gaunt factors of van Hoof et al. (2014). We compute the electron–neutral (e⁻ + H) free–free absorption coefficients based on the analytical approximation of Stallcop (1974), which agrees within 5% with accurate calculations (Bell & Berrington 1987; John 1994d) at wavelengths λ > 1 μm and temperatures T = 1400 to 10 080 K. We also include e⁻ + He (John 1994a,b) and e⁻ + H₂ (John 1994c) in the electron–neutral free–free opacity. In principle, electrons on polar, neutral molecules like CO and H₂O produce very large free–free opacities per molecule (John 1975b,a); however, they contribute relatively little to the total compared with the more abundant H,H₂, and He in our models. In order to describe the opacity at mid-infrared and shorter wavelengths (λ < 50 μm), we include non-relativistic Thomson scattering by free electrons and H⁻ bound-free absorption (cross-sections from McLaughlin et al. 2017b,a).

These absorption coefficients are conventionally expressed as an inverse length per unit density of the atom or ion and per unit electron pressure, κ, where the electron pressure is P_e = p_ek_BT. The optical depth in the continuum through a uniform medium of length Δz is thus $${\tau }_{\text{c}}=P_{\text{e}}(p_{\text{H}}{\kappa }_{\text{e},\text{neutral}}+p_{\text{e}}{\kappa }_{\text{e},\text{ion}})\Delta z.$$(8)

Under the assumption of a Gaussian line profile with FWHM Δυ, the corresponding optical depth at the centre of a line of element i in ion stage j is given by $${\tau }_{n^{\prime }{n}^{″}}=3.738\times {10}^{-7}{A}_{n^{\prime }{n}^{″}}{p}_{ij{n}^{″}}\left(\frac{g_{n^{\prime }}[1-\exp (-h\nu /k_{\text{B}}T)]}{g_{n^{″}}{(\nu /c)}^3\Delta \upsilon }\right)\Delta z,$$(9)

where $g_{n^{\prime }}$ and $g_{n^{″}}$ are the statistical weights of the respective upper and lower states of the transition n′ → n″, $A_{n^{\prime }{n}^{″}}$ is the spontaneous transition probability (s⁻¹), c is the speed of light [cm s⁻¹], and the line width has dimensions (km s⁻¹). The product $p_{ij{n}^{″}}\Delta z$ is the lower-state column density withdimensions [cm⁻²]. With the assumption of LTE, the excitation temperature of any transition is equal to T, the source function is equal to the Planck function B_ν(T), and the observable flux density of either the line peak or continuum is $$S_{\nu }=B_{\nu }(T)[1-\exp (-\tau )]\pi {\left(\frac{R}{D}\right)}^2,$$(10)

that is, the line and continuum emissions are treated separately. This requires, for instance, optically thin emission, a situation that applies reasonably well to our models. To summarise the model, Eqs. (2) through (10) show how the observable flux densities of line and continuum depend on the parameters n_H, T, θ (or R), ρ, and D for a specified set of elemental abundances, where n_H = $n_{{\text{H}}^0}+n_{{\text{H}}^{+}}+n_{{\text{H}}^{-}}+2(n_{{\text{H}}_2}+n_{{\text{H}}_2^{+}})+3n_{{\text{H}}_3^{+}}$.

5.5 Line broadening: Doppler motions and blending of emitters

The shapes and widths of the line profiles contain valuable information. The observed lines are fit very well by single Gaussian functions, with a FWHM of Δυ ≈ 8.5 km s⁻¹ (taken from the X26α line that is observed with the highest spectral resolution). It is significant that the 26α and 30α lines show the same value of line width, and that this line width is substantially larger than the thermal line widths at the relevant temperatures (e.g. for neutral Mg it is only 2.4 km s⁻¹ at T = 3000 K).

Typically this would be explained by a combination of thermal Doppler broadening and micro-turbulence, with only minor contributions from any larger scale kinematical motions. For simplicity, we adopt this description in the determination of the best-fit equilibrium model, the reference model, and separate the two components into thermal Doppler motions (Δυ_th) and micro-turbulent motions (Δυ_turb). We further assume that both have Gaussian broadening functions along the line of sight through the emitting region. The Gaussian widths add in quadrature Δυ² = $\Delta {\upsilon }_{\text{th}}^2$ + $\Delta {\upsilon }_{\text{turb}}^2$.

Broadening of the lines results also from the fact that each line is most likely a blend of lines due to a mixture of different species. The slight differences in reduced mass of the species lead to slight differences in line frequencies; see Eq. (1). These velocity differences are independent of n for any α transition.

In the modelling, we obtained solutions by linear combination of optical depths of overlapping lines of the different atoms in the Rydberg approximation. The best-fitting models have line optical depths that are less than or of the order of unity, whichjustifies simply summing the optical depth at each point on a densely sampled frequency grid. In the constructionof a theoretical profile, each of the line components is weighted by the relative abundance of the atomic species.

Rydberg atoms are very sensitive to perturbations by collisions with other atoms and electrons and by magnetic forces. Line profiles broadened by such perturbations provide limits on gas density and magnetic field strength as discussed in Appendix C and Sect. 6.2.2 below.

6 Results from modelling mm-wave RTLs

The model described in the preceding section can be used to gain insight into the conditions that make mm-wave RTLs observable towards stars that are only moderately warm. We have several observable properties and a model with four free parameters. We start by explaining the adopted procedure for finding the best-fit model.

6.1 Procedure for selecting the best-fit model

To find the model that best fits the observations is not trivial, and therefore we describe the procedure employed:

• Adopt a distance D and a set of elemental abundances.

• Select a value of the geometrical parameter ρ.

• For one line, specify the peak flux density S_ν and observed line width Δυ.

• Specify a grid of values of the parameters n_H and T.

• For each set of density and temperature, compute the peak line flux densities for a list of lines, recalling that the observed line might be an unresolved blend of several emitters.

• Iterate over values of the source radius R until the model matches the line flux.

• Examine the resulting grid of solutions to see which of them also satisfies the upper limit on the hydrogen line flux, the flux of the adjacent continuum, and the observed constraints on angular size of the source.

• For the same choice of D, abundances, and ρ, repeat the above steps for a second transition.

• Determine whether there is some parameter space that satisfies both transitions and their continuum fluxes simultaneously.

• Repeat the entire procedure for different choices of elemental abundances and ρ.

The line list for the RTLs includes two lines (H and C) and one blend (all abundant heavier elements that contribute to the blend, Na to Ni). The densities of the neutral species have been computed explicitly in strict LTE. The mathematical model is non-linear in several regards. The temperature enters in the arguments of exponential functions in the Saha and Boltzmann equations describing the ionisation balance and excitation. The source path length Δz = ρR is a term in each optical depth and thus enters in the arguments of exponential functions. The coupled parameter for the radius enters as R² in the flux densities. Thus, the solutions must be found by successive iterations. Remarkably, in the present case, the solutions appear to converge on a relatively small volume in parameter space.

Fig. 2 Left: LTE solutions with ρ = 0.1 that yield a line flux density of 7.3 mJy in the 30α blend for solar abundances. Blue contours display equal values of the angular diameter θ in mas. Dark red contours show the 232 GHz continuum flux density in mJy. Green contours show the flux density of the corresponding hydrogen line in mJy. The polygon covers the space of primary solutions that allow for ±10% changes on the continuum flux density and the angular diameter for illustrative purposes. The shaded rectangles indicate how the primary solution shifts when the value of ρ is changed by a factor of two (0.05 (cyan) and 0.20 (green)), and when a source size of 17 mas and ρ = 0.1 (violet) and 2 (blue) are adopted. Right: same as left panel but for the 26α blend and a line flux density of 25 mJy, and the 354 GHz continuum flux density in mJy. The cyan solid circle indicates the location of the reference model.

6.2 The reference model

The behaviour of the models is illustrated in two contour diagrams shown in Fig. 2. These diagrams of the (log n_H,T)-plane display contours of equal angular diameter (θ in mas, blue), continuum flux density (f_c in mJy, brown), and hydrogen line flux (f_H in mJy, green). These models were computed for ρ = 0.1 (a cylinder with a depth equal to a tenth of its radius), undepleted solar abundances, a nominal distance of D = 1.0 kpc, and the line flux densities given in Table 1, that is, all points in the diagrams give the correct observed line intensities. However, the observational results will limit the area of acceptable solutions, for instance, the upper limits on the H30α and H26α lines imply that the temperature must be T≲ 4000 K, for the adopted solar abundances, and the measured source size is about 11 mas.

The X30α transition provides the most restrictive constraints, because of the higher angular resolution of the observations. As shown in Fig. 2: left, there are two regions of the (n_H, T) plane that satisfy all of the observed constraints of the 30α transition and adjacent continuum. We pick the midpoint of the polygonic area in the left panel of Fig. 2 as the reference model. The physical characteristics of this region are T = 2800 K and n_H = 8 × 10¹¹ cm⁻³, and its size is R = 5.5 au. The resulting gas mass is ≈7 × 10⁻⁴ M_⊙, that is, only about 2% of the estimated mass of the CCS component. The ionisation fraction is low, ≈ 3 × 10⁻⁵, and CO dominates over H₂ with [CO]/[H₂] ≈ 6 and [CO]/[H] ≈ 3 × 10⁻⁴ (a significant fraction of the C is bound in CO). The corresponding H and C recombination lines are more than a million times weaker than the heavy-element blend. The properties of the reference model are given inTables 2 and 3. The corresponding solutions for the 26α transition are displayed in the right panel of Fig. 2. Remarkably, the solutions for this transition cluster in nearly the same part of parameter space as for the 30α transition.

The contours of the 232 GHz continuum flux density do not close in the upper part of the diagram. As a result, there is a secondary solution at higher density (2.2 × 10¹³ cm⁻³) and lower temperature (2050 K). However, this solution is implausible for two reasons. A very large continuum optical depth is required in order for this model to produce the observed line flux within the observed angular size. This would imply that the τ_c ~1 surface of the continuum source is much more extended than the line emission, in contradiction to the 232 GHz line and continuum maps. Secondly, the secondary solution implies much higher molecular abundances than the primary model because it is both denser and colder. The result is that the predicted LTE emission line fluxes of lines of CO and SiO in the observed bands are much larger than those observed; see Sect. 6.4.

There are no reasonable solutions for the X30α line with depleted interstellar abundances. Depleted abundances of the heavier elements (especially Mg, Si, and Fe) cause the corresponding hydrogen line to be stronger, relative to the lines of the heavier elements, than in the undepleted case. The fluxes of H30α and/or the 232 GHz continuum are too high over the entire parameter space explored (T = 1500 to 10 000 K).

Fig. 3 Computed model profiles of the 30α (left) and 26α (right) transitions, excluding collisional broadening and radiative damping, are shown in black (truncated at a value of0.5, the peak value is 1.0). The blended profiles of the eight atoms and isotopes are shown in colours to illustrate their relative contributions and the mass-shifts in frequency. All of the profiles are slightly asymmetrical owing to the fine structure as described in Appendix C. A turbulent velocity Δυturb = 7.5 km s−1 is adopted so that the fitted FWHMs of the model profiles match the observed width, 8.5 km s−1.

6.2.1 Identification of the line emitters

Our LTE models further indicate that the observed lines are blends of the corresponding RTLs of Na, Mg (including its three stable isotopes), Al, and Ca, with smaller contributions from Si and Fe if these elements are present in the gas at nearly solar abundances (even though Si has a greater abundance than any of these elements, its contribution to the line emission in strict LTE is smaller, because its ionisation potential is significantly higher); see Table 2. Mg is the largest contributor, providing 74% of the total intensity.

We computed theoretical line profiles for the 26α and 30α transitions for the blend of these eight atomic species (including isotopic variants) in the reference model. The computed model profile represents the weighted sum of the line profiles of the species. Each line component is broadened by a Gaussian of total Doppler width Δυ. The thermal width is calculated separately for each atom, because it depends upon the atomic mass, while a single value of the micro-turbulent width was assumed for all (Δυ_turb = 7.5 km s⁻¹ is adopted so that the fitted FWHMs of the model profiles match the observed width, 8.5 km s⁻¹). The computed profiles are displayed in Fig. 3. The profile shapes of the two transitions are very similar, although the 30α profile is slightly narrower in velocity units than the 26α profile (zero velocity corresponds to the rest frequency of ²⁴Mg in the simple Rydberg formula). The small shift of the blended profile relative to Mg is quantified in Fig. 4 through a least-squares fit of single Gaussian functions to the computed profiles. The fits determine a central rest frequency for the model blends, which are 232 022.635 ± 0.002 MHz and 353 808.691 ± 0.004 MHz for the 30α and 26α transitions, respectively. These values agree very well with the rest frequencies of the observed lines estimated assuming a systemic velocity of 41.5 km s⁻¹; see Table 1. The fits also show that the deviations from a Gaussian shape are small, less than 3.0 and 4.2% for the 30α and 26α transitions, respectively.

6.2.2 Line broadening

As outlined in Appendix C we can expect line broadening due to fine structure, collisional broadening, radiative damping, and Zeeman splitting in the presence of an external magnetic field, in addition to that of Doppler shifts due to thermal and turbulent motions.

Unfortunately, there is some uncertainty about the collisional broadening coefficients: adopting the nominal values, the total damping (collisional plus radiative) corresponds to a Lorentzian line shape with FWHMs of Γ = 2.63 and 3.28 MHz for 26α and 30α, respectively, in the reference model. With this choice, the convolved profile has a Doppler core that can be fit rather well by a Gaussian of Δυ = 8.5 km s⁻¹, as in the observed profile, provided that the turbulent width is Δυ_turb = 5.5 km s⁻¹. The deviation from a Gaussian profile is modest, reaching at most 7% of the peak intensity in the blue wing. According to the present theory, the electron-impact collisional broadening has a strong n-dependence at higher values of n. If this applies also to lower values of n, the 30α line is expected to be markedly broader than the 26α line if the density is high enough. This can be used to estimate an upper limit to the density because the observed lines are well fitted with Gaussians of the same width and no indications of Lorentzian wings. Using this (see Appendix C) it appears that we rule out densities higher than n_H ≳ 3 × 10¹² cm⁻³, a result independent of the excitation models.

Finally, models with Zeeman splitting have been computed assuming a uniform magnetic field throughout the RTL region. A strong upper limit to the magnetic field of 2 G can be set because otherwise the 30α line would have a larger line width than the observed one even in the absence of turbulence and damping (and would show a noticeable Zeeman triplet structure). Allowing a combination of turbulence and damping that smooths out the structure of the Zeeman triplet decreases the upper limit on the magnetic field to about 1 G.

Fig. 4 Computed model profiles of the 30α (left) and 26α (right) transitions in Fig. 3 are shown in black, together with the best-fitting single Gaussians in cyan. The blue solid line is the difference between them.

6.2.3 Dependence on geometry and distance

An analysis of the sensitivity of the model results to the adopted geometry and distance is presented in Appendix D and summarised in Tables D.1 and D.3. The primary model solution is not very sensitive to the thickness of the cylinder. As an example of the dependence on geometry, the primary models for ρ between 0.05 and 0.20 areshown in the contour diagram (left panel of Fig. 2), and they both fall within the hexagon of preferred solutions for ρ = 0.1. Even for ρ = 2, the temperature and density of the primary solution are not substantially different, T = 2650 K, n_H = 2 × 10¹¹ cm⁻³, and are certainly well within the uncertainties of our adopted model.

The primary model is not very sensitive to the source size either. If we adopt a source size of 17 mas, that is, the proper size of a uniform brightness disc for which the Gaussian size is 11 mas, and ρ = 0.1, the resulting temperature and density, 2725 K and 4 × 10¹¹ cm⁻³, are only marginally different from those of the reference model, and are well within the uncertainties of our adopted model. The effect of changing both the source size (17 mas) and the geometry (ρ = 2) is illustrated by the blue rectangle in the left panel of Fig. 2. We note that a larger source results in lower optical depths for both lines and continuum.

Similarly, we recomputed the models for different values of the distance, D = 0.5 and 2 kpc. The primary solutions are very similar to those for the adopted distance of 1 kpc because the radius R simply scales to match the observed angular size. Only the gaseous mass depends significantly on the distance. The primary solutions for ρ = 0.1 have gas masses of 1.2 × 10⁻⁴ and 4.2 × 10⁻³ M_⊙ for the closer and more distant alternatives, respectively, that is, the gas mass scales roughly as D^2.5.

6.3 Constraints from optical emission lines

HD 101584 has an optical spectrum that is rich in emission and absorption lines. Bakker et al. (1996) made a thorough analysis of their spectra obtained with the ESO Coudé Auxiliary Telescope, but we note here that they assumed an edge-on geometry as opposed to the face-on geometry advocated by Olofsson et al. (2019) and Kluska et al. (2020). Here we only select a few emission lines of low optical depth in order to perform a sanity check of the model derived from the radio data.

Numerous emission lines in the visible spectrum can be identified as low-excitation transitions in neutral atoms and a few singly ionised ions. While many of these lines may arise in the chromosphere and inner wind of the star, it is possible that some lines originate in the more extended RTL region that is responsible for the mm-wave line emission. We focused on lines that are most likely to be optically thin, and extracted VLT/UVES data from the ESO archive; see Sect. 3.2. A selection of such lines are listed in Table E.1. The amplitude, centre, and width of each line was determined from a simultaneous least-squares fit of a low-order polynomial continuum and several Gaussian functions to a small wavelength interval, typically 50 Å. The radial velocities, near 41.5 km s⁻¹ with respect to the LSR, are in excellent agreement with the systemic velocity of the star and with the derived radial velocities of the X26α and X30α lines at mm wavelengths. However, the observed line widths are significantly larger than those of the mm wave lines, in the range of 15 to 20 km s⁻¹. These types of lines were also discussed by Bakker et al. (1996).

It is interesting to consider whether the RTL-region contributes to any of the observed emission lines in the visible spectrum. Although a full analysis of the visible spectrum is beyond the scope of this paper, a few quantitative comments can be made. We tried to ensure that the description of the principal continuum opacities applies both in the long-wavelength limit at mm wavelengths and throughout the visible spectrum. At visible wavelengths, both H⁻ bound-free opacity and Thomson scattering become important compared with the bremsstrahlung that dominates at mm wavelengths. It is re-assuring that the reference model of the RTLs is not in conflict with the intensities of the visible emission lines; see Appendix E.

6.4 Predictions for radio lines

As can be seen in Table 2, the ionisation degree is low (≈10⁻⁴) in the RTL region and the warm gas consists of predominantly neutral hydrogen, while the molecular fraction is dominated by CO at an abundance of ≈3 × 10⁻⁴ with respectto hydrogen. The question is whether other atomic or molecular lines can be detected from the RTL region. The results presented inTable B.1 show that in the case of molecules, only CO and possibly SiO are viable candidates. Molecular nitrogen, N₂, has a high abundance, but its lines are due to weak electric quadrupole transitions. Predictions are made for the LTE fluxes of some interesting atomic and molecular lines in Table 4.

In order to compare with our ALMA data, we re-analysed the molecular line data for the CCS component (see Sect. 2) using only the data obtained with the most extended configuration in Band 6, which therefore has the highest angular resolution, ≈ 0.′′03, which is three times the size of the RTL region. The S/Ns are high enough to determine the source sizes from the global line profiles, but not high enough to determine their brightness distributions as a function of velocity. We present the results for the CO(2–1) and SiO(5–4) lines because these species have high abundances in the RTL region, and for the H₂S(2₂₀−2₁₁) and SO₂ (16_3,13−16_2,14) lines because they provide additional kinematic information. The SO₂ line comes from a transition with an upper state energy corresponding to 148 K, that is, a transition at high energy compared to those of the other molecular lines. We also added the upper limit obtained for the CO(υ = 1, J = 3–2) line at an angular resolution of ≈ 0.′′2. The results are summarised in Table 5.

All brightness distributions are larger than the RTL region by about a factor of ten or more (all sizes given here are FWHM of a Gaussian fitted to the brightness distribution and deconvolved with the beam). Except for the SiO line, the lines are reasonably fitted by a Gaussian line profile. For CO, the line width is substantially larger than that of the RTLs, while those of H₂S and SO₂ are substantially smaller. Only the SiO line has a line width comparable to the RTL line width, but this line is clearly non-Gaussian; it is double-peaked, a result corroborated by the profile of the ²⁹ SiO(5–4) line. Therefore, the source sizes and the line profiles clearly show that the RTL region contributes only marginally, if at all, to the observed molecular lines. Only, the vibrationally excited CO rotational line deviates from this conclusion, but whether or not LTE is a good approximation for the vibrational excitation of CO is uncertain.

We also estimated the molecular flux densities within an aperture of 0.′′ 03 at the centre position. For CO and SiO, these come out close to those estimated for the reference model (Table 4). This suggests that the RTL region may contribute significantly to these observed intensities. In contrast, the H₂S and SO₂ lines must come from a region well outside the RTL region because the predicted abundances of these molecules are very low in the RTL region (Table B.1) and the lines are very narrow.

With the present data it is not easy to paint a consistent picture of the geometry and kinematics of the different line brightness distributions of the CCS component, including the RTL region. It would seem reasonable that the density and temperature of the gas declines with distance from the centre. This would explain why the more easily excited CO line emission is more extended than those of the other species. Chemistry may also play a role, restricting the three-atom molecules to the higher densities. The nature of the velocity field remains an open question because the observed lines (from mm wave to optical) show a number of different line profiles and line widths.

7 Discussion

7.1 The origin of the RTL-region

HD 101584 is believed to be in the aftermath of a CE process that ended about 800 yr ago without a stellar merger (Olofsson et al. 2019). This resulted in the termination of the RGB evolution of the primary and only its core remains, while the low-mass companion lies at a distance of about 0.5 au from HD 101584. The complexity of this process, for which little observational information exists, and the absence of models for stars in this stage make it difficult to speculate about the origin of the warm gas. On the other hand, it does not appear unlikely that about 10⁻³ M_⊙ of warm gas can leave the star at the end of the CE evolution. The gas is most likely stationary – although small effects of expansion and/or rotation cannot be excluded – but must be characterised by substantial micro-turbulence (especially if the optical lines discussed above also originate from this region). The exact process responsible for keeping the gas heated and turbulent remains unknown.

7.2 A new tool for studying cooler stars?

The detection of RTLs from low-ionisation-potential elements towards HD 101584, combined with the high sensitivity of ALMA, opens up the possibility of also detecting such lines towards other objects of cool to moderate temperature. In particular, red giants on the RGB or AGB may be detectable if they have enough gas in their immediate surroundings with temperatures around 3000 K and particle densities of the order 10¹² cm⁻³. The studies of radio continuum emission from low-ionisation gas in the extended atmospheres of AGB stars indicate that this may very well be the case (Reid & Menten 2007; Matthews et al. 2018; Vlemmings et al. 2019). To this category we can also add red supergiants. At temperatures in excess of about 4000 K we can expect the hydrogen lines to also become detectable. This is, for instance, the case for pre-PNe (without ionising radiation) and objects like yellow supergiants (e.g. Betelgeuse, O’Gorman et al. 2017) and hypergiants. The H30α line has been reported in ALMA observations of the Fried Egg Nebula (Wallström et al. 2017), which is associated with the yellow hypergiant star Henize 3–1379 (IRAS 17163–3907).

Close binaries – like HD 101584 – of various types are also likely to be of interest in this context. In this category are the symbiotic stars. Several such stars, notably RR Tel, show detectable free–free continuum emission at cm wavelengths. An interesting example in the high-mass regime is VV Cep, which is an eclipsing binary star system consisting of a red supergiant that fills its Roche lobe when closest to the companion, a Be star that is most likely on the main sequence (Hagen Bauer et al. 2008). Here substantial amounts of warm gas can be produced through shocks and the hot radiation field.

The lines due to Rydberg transitions may therefore provide independent information on the warm gas in the vicinity of moderate-temperature stellar objects, and hence provide constraints on the origin and heating of such regions. Rydberg transitions of heavy elements are also susceptible to strong Zeeman splitting, and consequently can be used to estimate stellar magnetic fields, a property that is notoriously difficult to measure for evolved stars (e.g. Vlemmings 2014).

8 Conclusions

We identified two lines towards the post-giant star HD 101584 that can be attributed to the Rydberg transitions 30α and 26α of neutral atoms of elements heavier than carbon. We used a model in strict LTE and a source geometry in the form of a geometrically thin disc to interpret the data and to identify the line carrier(s). A solution that fits all the observational constraints (i.e. the measured line flux densities, upper limits to the continuum flux densities, and the estimated size of the emitting region) includes a gas mass of about 10⁻³ M_⊙, a density of about 10¹² cm⁻³, a temperature of about 2800 K, a size of ≈10 au, and a turbulent velocity field in the range 5.5–7.5 km s⁻¹, assuming a distance of 1 kpc and solar abundances of the elements. The element Mg contributes about 75% of the line fluxes. The contributions from Na, Al, and Ca lie in the range of 5–10% each. The corresponding H and C lines are much weaker. Using the Zeeman effect, we set an upper limit to the average magnetic field in the line-emitting region of 1 G.

Further observations are required to determine in more detail the geometry and kinematics of the region where the Rydberg transition lines are formed. The gas appears to be stationary, and the process responsible for keeping the gas heated and turbulent over a region of about 10 au remains to be identified; it may be connected to processes active during the CE phase.

Although HD 101584 is not a standard object, we speculate that the types of lines discussed here should be detectable towards other moderately warm and evolved objects, in particular taking the sensitivity of ALMA into account, and serve as a new tool for their study.

Acknowledgements

We are grateful to the referee (Albert Zijlstra) for providing insightful and constructive comments on the paper. H.O. and W.V. acknowledge support from the Swedish Research Council. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2012.1.00248.S, #2015.1.00078.S and #2018.1.00753.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. H.O. acknowledges support from the Nordic ALMA Regional Centre (ARC) node based at Onsala Space Observatory. The Nordic ARC node is funded through Swedish Research Council grant No 2017-00648. This research has made use of the services of the ESO Science Archive Facility. Based on observations collected at the European Southern Observatory under ESO programme 266.D-5655(A).

Appendix A Elemental abundances

The abundances are expressed on the standard logarithmic scale [H] = 12.0, so that the relative number densities are log (n_X∕n_H) = Y (X) − Y (H) + logd(X); see Table A.1. For Set 1 we adopt solar abundances and all d(X) ≡ 1.0). The depletion factors for Set 2 are typical of well-mixed, diffuse interstellar gas where refractory elements have already been condensed into dust. Only elements up to Ni with ionisation potentials lower than that of hydrogen are considered.

Appendix B Chemistry

Table B.1 lists the molecules that have been included in the chemistry and ionisation balance in the LTE models. For each molecule X, the fractional abundance is given with respect to H in the reference model. The adopted dissociation energy D₀ is expressed in wavenumber units of cm⁻¹ (one electron volt is eV = 8065.5439 cm⁻¹), and the source of the partition function is listed. For most of the included molecules, tables of partition functions were taken from the ExoMol Database (Tennyson et al. 2016)⁶ and used via interpolation. In a few cases, notably H₂, CO, and N₂, the partition functions were calculated directly by summing over states listed in extensive spectroscopic tables compiled for the Leiden Atomic and Molecular Database for Astronomy (LAMDA, Schöier et al. 2005; van der Tak et al. 2020)⁷.

Although 30 diatomic and triatomic molecules have been included in the equilibrium computation of ionisation, only a few of them have significant abundances. Notice that the primary solutions reflect conditions in which a small fraction of the hydrogen is in molecular form, while a significant fraction of the total carbon and oxygen is consumed in CO. The relatively high abundances of CO, N₂, and SiO reflect their high dissociation energies. At T ≈ 2800 K and n_H ≈ 10¹² cm⁻³ in a gas of solar abundances, the presence of these most abundant molecules does have a small effect on the overall ionisation balance and thus the intensity of the free-free continuum emission. In addition, these models can be used to predict the intensity of molecularline emission that should accompany the atomic lines under the adopted assumptions about strict LTE and solar abundances.

Appendix C Additional line broadening

Because Rydberg atoms interact strongly with electrons, other atoms, radiation, and magnetic fields, it is worthwhile exploring whether or not the line profiles contain additional information about the conditions in the emitting regions. This can be done through computation of theoretical line profiles based on the conditions in the reference model described above. In order to do this properly, the fine structure of the Rydberg transitions will have to be added.

C.1 Fine structure

Rydberg states of principal quantum number n possess slightly non-degenerate fine structure in their orbital substates (quantum number ℓ = 0 to n − 1). Hence, a Rydberg transition n′ →n″ is composed of many fine-structure components (n′, ℓ′) →(n″, ℓ″) that obey the selection rules ℓ′ = ℓ″ ± 1. Although the fine-structure splitting is negligible at n = 25 to 30 in hydrogen (Towle et al. 1996), atoms like Mg, Al, and Ca have parental ionic cores with rather high values of the static dipole polarisability, which gives rise to larger fine-structure splittings than in hydrogen. Thus, the 26α transition of each atomic emitter has 51 fine-structure components and the model profile for this transition is a blend of 408 lines for the atomic species considered in the reference model; see Table C.1. The same figures for the 30α transition are 59 and 472. In the construction of a theoretical profile, each of the fine-structure components is weighted by the relative abundance of the atomic species and by an intensity factor, which is proportional to the upper-state statistical weight multiplied by the spontaneous transition probability.

We computed the energy of each state of each atom, labelled (n, ℓ), from the solution of the Dirac equation for a hydrogenic atom with n ≫ 1 as presented in Eqs. (2) through (4) of Towle et al. (1996). An additional term, − α_dP(n, ℓ), is included to account for the polarisation of the ionic core, where P(n, ℓ) is a function of n and ℓ, and α_d is the electric dipole polarisability (Chang & Noyes 1983; Chang 1987). In H itself, with a bare proton core, the polarisability vanishes and this term is zero. In Mg, on the other hand, the Mg⁺ core has α_d = 33.05±0.03 $a_0^3$ (atomic units, where a₀ is the Bohr radius (Mitroy et al. 2010). We computed an ℓ-resolved line list for each transition n′ →n″. The number of fine-structure components (distinct lines) in an α transition is 2n″–1, where double-prime indicates the lower state. The spontaneous transition probabilities (Einstein A-coefficients) for these hydrogenic transitions were computed through use of the computer program of Storey & Hummer (1991). We assume that pure L-S coupling applies. For Rydberg states in the limit of high-n, all atoms are assumed to behave like degenerate hydrogenic spin-doublets of total angular momentum j = ℓ ± 1/2. The degeneracy is lifted in the presence of external fields, as we explain below in the section on Zeeman splitting. In the absence of external fields, the blending of lines of different atoms and the overlap of their fine-structure components broaden the line profile and make it slightly asymmetrical. The fine-structure splittings in the 26α and 30α transitions correspond to only a few tenths of a km s⁻¹, and hence will not contribute substantially to the line widths measured here.

C.2 Collisional broadening

Atoms in Rydberg states interact with electrons and neutral H atoms through elastic and inelastic collisions that perturb their energies and thus broaden their line profiles. This collisional broadening depends on the densities of the collisionpartners and on the kinetic temperature. In ionised nebulae, recombination lines suffer collisional broadening mainly with electrons and ions. In contrast, our reference model is only weakly ionised, with n_e ∕n_H = 3.4 × 10⁻⁵; as a result, neutral hydrogen and electrons are the most important collision partners.

The electron-impact broadening of nebular recombination lines has been discussed in the literature; unfortunately, most of that work is not applicable here. Watson (2006) presented a formula for the collisional broadening based on the classic work of Griem (1967). As pointed out by Alexander & Gulyaev (2016), that formula contains errors. Beyond that, Griem’s approximation was originally intended for cm- and dm-wave transitions with n > 100; the approximation breaks down when n < 70 and is thus not appropriate for the transitions of interest here. Brocklehurst & Leeman (1971) computed the electron–hydrogen cross-sections and rates numerically, and evaluated the line-broadening without the restrictive assumptions of semi-classical theory. Their results could be approximated by a simple fitting formula for the collisional line width of α transitions, $$\frac{\Delta {\nu }_{\text{e}}}{{\nu }_0}=1.43\times {10}^{-5}{\left(\frac{n}{100}\right)}^{7.4}{\left(\frac{{10}^4\text{K}}{T}\right)}^{0.1}\left(\frac{n_{\text{e}}}{{10}^4{\text{cm}}^{-3}}\right),$$(C.1)

where T and n_e are the temperature and number density of electrons, respectively (Brocklehurst & Seaton 1972). Because the original computations did not extend to n < 109 and because the dependence on n is so steep, it is not certain that this formula is accurate at n = 26 and 30. If we simply apply the formula for T = 2805 K and n_e = 2.8 × 10⁷ cm⁻³, as in the reference model, then the electron-impact widths (FWHM) are Δυ_e = 0.63 and 1.83 km s⁻¹ for n = 26 and 30, respectively. We note that the collisional broadening produces a Lorentzian line shape rather than a Gaussian. Moreover, if electron-impact broadening were to dominate, the width of the 30α line would be almost three times larger than that of 26α in velocity units according to the formula of Brocklehurst & Seaton (1972).

Collisional broadening by neutral atoms has also been considered, particularly in connection with laboratory experimentson Rydberg atoms in gases of noble or alkali atoms (Omont 1977; Kaulakys 1984). Neutral–atom collisions are also important in line-broadening in stellar atmospheres. Specifically, Hoang-Binh & van Regemorter (1995) applied the impulse approximation to compute the broadening of mid-infrared lines of Mg I by atomic hydrogen in the photosphere of the Sun, where these lines are observed in emission at wavelengths near 12 μm. These are also Rydberg transitions, but at the relatively low values of the principal quantum number involved, n = 5, 6, and 7, the orbital and spin-orbit fine structures are well resolved. If we adapt the theoretical widths calculated by Omont (1977) to the case of H on Mg, we estimate Δν_H ≈ 1.8 MHz at both n = 26 and 30, corresponding to velocity widths Δυ_H = 1.5 and 2.2 km s⁻¹, respectively,in the reference model. Even though the neutral–neutral interactions are much weaker than electron–neutral interactions involving Rydberg states, the estimated contributions to the line broadening are comparable in our case because of the low electron fraction.

In both the cases of electron and neutral collisions, the collisional broadening should be re-computed numerically with modern methods, because both elastic and inelastic collisions are important and earlier approximations at n > 100 are known to break down when n ~20 to 30, but this is beyond the scope of this paper. The above estimates suggest that the collisional broadening might become noticeable ifthe line profile can be observed with sufficiently high S/N to reveal Lorentzian wings. However, the existing data show line profiles that are fit very well by single-Gaussian functions, meaning that, even with the uncertainties in the theoretical line widths, the current data on HD 101584 suggest that the neutral density in the emitting region cannot be as large as 10¹³ cm⁻³, otherwise the Lorentzian wings would be apparent. This kind of limit on density is not model dependent; moreover, it could be refined with improved theory.

C.3 Radiation damping

Atoms in Rydberg states also interact with electromagnetic radiation. The natural line broadening (radiation damping) via spontaneous transitions can be further enhanced through stimulated emission and absorption in the ambient radiation field. Considering a transition from upper state u to lower state l, the energies of these states are diffuse because both states have finite lifetimes. The Lorentzian FWHM of the radiation damping is Δν_rad = γ∕(2π) where $$\begin{array}{lll}{\gamma }_{\text{ul}}\hfill & =\hfill & {\sum }_{\text{i}<\text{u}}\left(A_{\text{ui}}+{\overline{J}}_{\text{ui}}{B}_{\text{ui}}\right)+{\sum }_{\text{k}>\text{u}}{\overline{J}}_{\text{ku}}{B}_{\text{ku}}+{\sum }_{{\text{i}}^{\prime }<\text{l}}\left(A_{{\text{li}}^{\prime }}+{\overline{J}}_{{\text{li}}^{\prime }}{B}_{{\text{li}}^{\prime }}\right)\hfill \\ \hfill & \hfill & +{\sum }_{{\text{k}}^{\prime }>\text{l}}{\overline{J}}_{{\text{kl}}^{\prime }}{B}_{{\text{kl}}^{\prime }},\hfill \end{array}$$(C.2)

and A_ik is the spontaneous transition probability for transition i → k, and ${\overline{J}}_{\text{ik}}$ the angle-averaged intensity of radiation at the frequency ν_ik of the transition. The Einstein coefficients for stimulated emission, B_ik = A_ik/(2hν³/c²), and absorption, B_ki = g_iB_ik/g_k are related to A_ik and to the statistical weights of the upper and lower states g_i and g_k. When the radiation can be represented by a blackbody diluted by a geometrical factor w ≤ 1, then $${\overline{J}}_{\text{ik}}=w{B}_{{\nu }_{\text{ik}}}(T_{\text{rad}})=\frac{w2h{\nu }_{\text{ik}}^3/c^2}{\exp \left(h{\nu }_{\text{ik}}/k_{\text{B}}{T}_{\text{rad}}\right)-1},$$(C.3)

in terms of the Planck function B_ν at temperature T_rad. Thus, the damping coefficient can be re-written $$\begin{array}{lll}{\gamma }_{\text{ul}}\hfill & =\hfill & {\sum }_{\text{i}<\text{u}}{A}_{\text{ui}}\left(1+b({\nu }_{\text{ui}})\right)+{\sum }_{\text{k}>\text{u}}{A}_{\text{ku}}b({\nu }_{\text{ku}})\frac{g_{\text{k}}}{g_{\text{u}}}\hfill \\ \hfill & \hfill & +{\sum }_{{\text{i}}^{\prime }<\text{l}}{A}_{{\text{li}}^{\prime }}\left(1+b({\nu }_{{\text{li}}^{\prime }})\right)+{\sum }_{{\text{k}}^{\prime }>\text{l}}{A}_{{\text{kl}}^{\prime }}b({\nu }_{{\text{kl}}^{\prime }})\frac{g_{{\text{k}}^{\prime }}}{g_{\text{l}}},\hfill \end{array}$$(C.4)

where b(ν) = w/$\left(\exp (h\nu /k_{\text{B}}{T}_{\text{rad}})-1\right)$. Our model of the line-emitting region assumes strict LTE at the kinetic temperature T. This means that w = 1 and T_rad = T should be a good approximation at frequencies less than that of the 26α transition, owing to the inverse frequency dependence of the predominant opacity source. However, some of the stimulated transitions occur at much higher frequencies where the emitting region may be transparent to radiation from the much hotter visible photosphere of HD 101584, for which w < 1 and T_rad ~ 8500 K.

The computation of Δν_rad is straightforward for hydrogenic transition frequencies and transition probabilities and is done here by summing all rates for states up to n = 1000. The radiation damping of the 26α and 30α lines at T_rad = 2805 K and w = 1 gives Δν_rad corresponding to 0.072 and 0.083 km s⁻¹, respectively. Even in the case of the hot photosphere T_rad = 8500 K shining through the emission region at w = 1/2, the damping widths would be increased by only about 50%. In either case, the radiation damping is negligible in comparison with other broadening mechanisms. In general, this broadening effect on Rydberg transitions might be most important at very low radio frequencies (n > 200) in gas exposed to intense non-thermal continuum radiation, as in the narrow-line regions of quasars (cf. Wadiak et al. 1983).

Both the collisional line broadening and the radiation damping produce Lorentzian profiles. Two Lorentz functions combine to form a new Lorentz function whose width is simply the sum of the first two Lorentz widths. The total profile is obtained by a convolution of the final Lorentz function with the profile for the blend.

C.4 Magnetic fields and Zeeman broadening

Magnetic fields can cause a line broadening in radio and mm-wave RTLs through the Zeeman effect. The blending of lines of different emitting atoms, each with the orbital fine structure described above, is further complicated in the presence of external fields. A one-electron, hydrogenic state with spin s = 1/2, labelled by the quantum numbers n, ℓ, and j in L-S coupling, has a pair of degenerate levels for each j = ℓ ± 1/2. The degeneracy, that is, the number of unresolved sublevels, is 2j+1. A magnetic field of strength B lifts this degeneracy, and yields splittings in the weak-field limit of the Zeeman effect given by $$E(n,\mathcal{l},j,m_j)=E(n,l)+g{\mu }_{\text{B}}B{m}_j,$$(C.5)

where m_j = −j, − j + 1, …, + j − 1, + j. In L-S coupling the Landé g-factor is $$g=1+\frac{j(j+1)-\mathcal{l}(\mathcal{l}+1)+s(s+1)}{2j(j+1)},$$(C.6)

and g ≈ 1 for the levels of high j and ℓ that contribute most to the line profile. We note that μ_B B∕h = 1.3996 B MHz (B is in the unit of gauss), so that a field strength of the order of B ≈ 1 G can producea splitting of the order of 1 MHz, which may be significant in the lines of interest here.

C.5 The blended line profile

We have computed theoretical line profiles for the 26α and 30α transitions taking into account all of the effects of blending and broadening. Eight atomic species are included, as listed in Table C.1. The computed model profile represents the weighted sum of the scaled optical depths of the various (n, ℓ) line components. Each line component is broadened by a Gaussian of total Doppler width Δυ. The thermal width is calculated separately for each atom, because it depends upon the atomic mass, while a single value of the micro-turbulent width was assumed for all.

Further, we computed a series of model profiles with several values of the damping width, Γ, the FWHM of the Lorentzian function that describes the combined effect of collisional broadening and radiation damping. The total profile is a convolution of the blend with a Lorentzian function of width Γ. The collisional broadening term scales directly with density. As discussed above, there is some uncertainty about the collisional broadening coefficients: if we adopt the nominal values, then the total damping is Γ = 2.63 and 3.28 MHz for 26α and 30α, respectively,for the reference model. With this choice, the convolved profile has a Doppler core that can be fit rather well by a Gaussian of Δυ = 8.5 km s⁻¹, like the observed profile, provided that the turbulent width is Δυ_turb = 5.5 km s⁻¹ (in contrastto the best model with no damping wings that requires Δυ_turb = 7.5 km s⁻¹). The deviation from a Gaussian profile is modest, and reaches about 7% of the peak intensity in the blue wing. It therefore appears that the collisional broadening effect will not be conspicuous in the analysis of only one line.

The situation becomes more interesting if the electron-impact collisional broadening does indeed have the strong n-dependence predicted by the existing theory at higher values of n. In that case, the best-fitting Gaussian width of the 30α transition is expected to be somewhat different from that of the 26α transition. Figure C.1 shows the corresponding 30α profiles for the conditions that yield effective line widths of 8.5 km s⁻¹ in the 26α line. Here, the nominal collisional broadening of the 30α line gives a profile with a best-fitting Gaussian width that is clearly larger, Δυ = 10.2 km s⁻¹. The collisional broadening width would need to be a factor of two times smaller in order for the effective widths of the two transitions to be approximately the same, if both lines are truly formed at the same density and temperature and turbulent width. This factor of two is within the uncertainty of the extrapolated theory of collisional broadening. On the other hand, collisional broadening two times larger than the nominal case would yield noticeable Lorentzian wings and a 50% difference in fitted FWHM between the two transitions. These model profiles were computed with reference to an LTE excitation model with density n_H ≈ 8 × 10¹¹ cm⁻³. Thus it appears that the observed line profiles can already rule out densities higher than n_H ≳ 3 × 10¹² cm⁻³, independent of excitation models.

Finally, we have computed models with fully expressed Zeeman splitting for several combinations of magnetic field strength B (assumed uniform through the source), turbulent width Δυ_turb, and damping width Γ. A few examples of Zeeman-broadened profiles are shown in Fig. C.2. The combination of turbulence and damping tends to smooth out the structure of the Zeeman triplet when B < 2 G, so that the 26α profiles canall be made to approach Gaussian shapes with Δυ = 8.5 km s⁻¹ when the damping is small and B ≤ 1.5 G). However, the corresponding 30α profiles have Gaussian-fitted widths Δυ ≥ 9.0 km s⁻¹ when B ≥ 1.0 G. As the observed line widths of the two transitions are approximately equal, this suggests an upper limit to the field strength of 1 G when the damping is small. The observed width of the 30α transition sets a strong upper limit of B < 2 G, because the fitted width of the computed profile is 9.2 km s⁻¹ even when Δυ_turb = 0 and Γ = 0 at B = 2 G. Again, the strong limit on magnetic field is independent of any excitation model, and is insensitive to uncertainties in the collisional broadening theory, but it depends on the field geometry.

Fig. C.1 Series of model profiles of the 30α line for different values of the Lorentzian width Γ and turbulent widths that give Δυ = 8.5 km s−1 in the 26α line. The fitted Gaussian widths are indicated: they exceed the observed width when Γ > 1.7 MHz. The deviation from a single Gaussian is shown for the nominal case as a dashed green curve.

Fig. C.2 Series of model profiles including Zeeman splitting. Examples with small damping (Γ = 0.34 MHz in 26α) are shown in black. The structure of the 26α profile is smoothed by turbulence and damping to approximate a Gaussian shape with B in the range 0.5–1.6 G. The corresponding 30α profile is measurably broader, >9.0 km s−1 when B >1 G. A strong upper limit of B < 2 G is implied bythe fact that the 30α width is too large even when turbulence and damping are turned off, which then leaves the Zeeman triplet resolved (green curve).

Appendix D Dependences on geometry and distance

Additional models were computed by us in order to explore the sensitivity of the solutions to the geometrical parameter ρ and the distance D. We considered values of ρ from 0.05 to2.0 and distance ranges from 0.5 to 2.0 kpc. All models were constrained to fit the observed results at 232 GHz: θ = 11.0 mas, S_ν (X30α) = 7.3 mJy, and S_c (232) = 7.5 mJy. Details of additional models are summarised in Tables D.1 and D.3.

The models all behave similarly and the differences in best-fitting temperature (within 100 K) and total density (within a factor of 3) are relatively small. The primary solutions all have continuum optical depths of the order of unity at 232 GHz, and the X30α line optical depths for each contributing element are small enough that simply combining them in the optically thin limit is a good approximation.

Appendix E Optical lines

Table E.1 shows a comparison with integrated line fluxes predicted by the reference model, including a correction for extinction. The simple equilibrium model was not tuned in any way to match the optical spectrum. The model nicely reproduces the integrated fluxes of the two Ca II forbidden lines. It under-predicts the 6300 and 6363 Å lines of O I by a factor of six and the two C I lines by a much larger factor. However, the amounts of atomic O and C in the gas are very sensitive to the abundance of CO, which in turn depends quadratically on density and exponentially on temperature in equilibrium. In the model, all of the tabulated lines have a peak optical depth τ < 0.7. Weak and forbidden lines were selected for a preliminary analysis because lines of higher optical depth would require a much more sophisticated treatment of radiative transfer through an extended atmosphere. Such analysis is far beyond the scope of the present work. Even so, it would be interesting for future work to consider that the RTL region can explain the fluxes of prominent forbidden lines in the visible spectrum within an order of magnitude.

References

All Tables

All Figures

	Fig. 1 X30α (upper spectrum; scaled by a factor of three and biased by 25 mJy, 3 km s−1 resolution) and X26α (lower spectrum; 1 km s−1 resolution) lines observed towards HD 101584. Continuum emissions have been subtracted.
In the text

Fig. 2 Left: LTE solutions with ρ = 0.1 that yield a line flux density of 7.3 mJy in the 30α blend for solar abundances. Blue contours display equal values of the angular diameter θ in mas. Dark red contours show the 232 GHz continuum flux density in mJy. Green contours show the flux density of the corresponding hydrogen line in mJy. The polygon covers the space of primary solutions that allow for ±10% changes on the continuum flux density and the angular diameter for illustrative purposes. The shaded rectangles indicate how the primary solution shifts when the value of ρ is changed by a factor of two (0.05 (cyan) and 0.20 (green)), and when a source size of 17 mas and ρ = 0.1 (violet) and 2 (blue) are adopted. Right: same as left panel but for the 26α blend and a line flux density of 25 mJy, and the 354 GHz continuum flux density in mJy. The cyan solid circle indicates the location of the reference model.

In the text

Fig. 3 Computed model profiles of the 30α (left) and 26α (right) transitions, excluding collisional broadening and radiative damping, are shown in black (truncated at a value of0.5, the peak value is 1.0). The blended profiles of the eight atoms and isotopes are shown in colours to illustrate their relative contributions and the mass-shifts in frequency. All of the profiles are slightly asymmetrical owing to the fine structure as described in Appendix C. A turbulent velocity Δυturb = 7.5 km s−1 is adopted so that the fitted FWHMs of the model profiles match the observed width, 8.5 km s−1.

In the text

	Fig. 4 Computed model profiles of the 30α (left) and 26α (right) transitions in Fig. 3 are shown in black, together with the best-fitting single Gaussians in cyan. The blue solid line is the difference between them.
In the text

	Fig. C.1 Series of model profiles of the 30α line for different values of the Lorentzian width Γ and turbulent widths that give Δυ = 8.5 km s−1 in the 26α line. The fitted Gaussian widths are indicated: they exceed the observed width when Γ > 1.7 MHz. The deviation from a single Gaussian is shown for the nominal case as a dashed green curve.
In the text

Fig. C.2 Series of model profiles including Zeeman splitting. Examples with small damping (Γ = 0.34 MHz in 26α) are shown in black. The structure of the 26α profile is smoothed by turbulence and damping to approximate a Gaussian shape with B in the range 0.5–1.6 G. The corresponding 30α profile is measurably broader, >9.0 km s−1 when B >1 G. A strong upper limit of B < 2 G is implied bythe fact that the 30α width is too large even when turbulence and damping are turned off, which then leaves the Zeeman triplet resolved (green curve).

In the text