© Ch. Rab et al. 2022

1 Introduction

Understanding the dispersal and evolution of protoplanetary disks, the birthplaces of planets, is crucial to understanding the planet-formation process. Disk winds, either thermally or magnetically driven, likely play a crucial role in the disk mass loss and – in the case of magnetically driven winds – angular momentum extraction. Although several theoretical models exist to understand the physics of disk winds and their impact on disk evolution (see reviews of Gorti et al. 2016; Ercolano & Pascucci 2017; Pascucci et al. 2022; Lesur et al. 2022), questions remain as to: whether magnetohydrodynamic (MHD) or thermally driven photoevaporative winds dominate in disks; the nature of their interplay; and whether or not the outflows evolve from being MHD-dominated in young objects to thermally dominated at later ages (e.g. Ercolano & Pascucci 2017). Observationally constrained theoretical models are necessary to tackle these questions. Specifically, wind dynamics models are necessary as an input to astrochemical models (especially for molecular tracers) in order to obtain synthetic observables.

The most commonly observed tracers of outflows and winds in the protoplanetary disk stage are forbidden lines of atoms and low-ionisation species such as [OI], [SII], [NII], [FeII], and [NeII], with the most prominent being the [OI] 0.63 µm (e.g. Baldovin-Saavedra et al. 2012; Rigliaco et al. 2013; Simon et al. 2016; Fang et al. 2018; Banzatti et al. 2019). High-spectral-resolution observations (∆v ≈ 4–7 km s⁻¹) of the [OI] 0.63 µm commonly show complex line profiles that can be decomposed into several Gaussian components. The two main components are a high-velocity component (HVC), showing a shift of the peak flux location of v_p < −30 km s⁻¹ (relative to the systemic velocity) and a low-velocity component (LVC) showing blueshifted emission with v_p > −30 km s⁻¹. The HVC is often attributed to jets and strong outflows, whereas the LVC is supposed to trace slower disk winds. The LVC is also often further decomposed into two Gaussian components, a broad LVC (BLVC; i.e. full width at half maximum > 30 km s⁻¹) that dominates the line wings and a narrow LVC (NLVC) with a blueshifted peak of a few km s⁻¹, but single Gaussian fits to LVCs seem to be as common (see e.g. Banzatti et al. 2019).

Alternatives to emission lines are absorption lines, such as CII λ1335 Å, which was used by Xu et al. (2021) to study winds and outflows in T Tauri stars. These authors find that their findings can be best explained by a combination of magnetically and thermally driven flows. Molecular tracers of slow disk winds observed at high spectral resolution are less common and are mostly restricted to CO (e.g. Pontoppidan et al. 2011; Klaassen et al. 2013; Banzatti et al. 2022) and H₂ (e.g. Beck et al. 2008; Beck & Bary 2019). However, molecular tracers are of particular interest as they can trace different physical regions of the wind than the atomic tracers and provide further constraints for measuring mass-outflow rates.

A consistent physical interpretation of all the available observational data is challenging and often the different tracers are interpreted individually, in particular as simultaneous observations of multiple tracers come with significant difficulty. However, Gangi et al. (2020) presented new data for the [OI] 0.63 µm and the ortho–H₂ (hereafter o–H₂) 2.12 µm spectral lines for a sample of 36 young disk-bearing stars. The data set of these latter authors is unique as both lines were observed simultaneously, and for seven targets, where both lines are detected, kinematic signatures of winds for both tracers were identified. Gangi et al. (2020) interpreted their data as a potential indication of magnetically driven winds, as their analysis indicates that both spectral lines originate from similar regions of the disk wind. Also, several of the previously mentioned observational studies postulate that MHD-driven winds (e.g. Fang et al. 2018; Whelan et al. 2021) or a combination of MHD- and thermally driven winds (e.g. Xu et al. 2021) are common in T Tauri stars. The wealth of observational data provides strong constraints for theoretical models. However, a direct comparison of the observational data to the model is challenging, as it requires detailed models for the disk dynamics, evolution, thermal structure, radiation environment, and chemistry, as well as proper modelling of synthetic observables.

Several disk-wind modelling approaches exist. Ercolano et al. (2008a) and Gorti et al. (2009) presented static photoevaporative disk-wind models with a focus on radiative transfer and ionisation physics, with their main aim being to interpret atomic wind tracers. The first models that coupled X-ray and extreme ultraviolet (EUV) radiation physics to hydrodynamic models for photoevaporative winds were performed by Owen et al. (2010) and more recently by Picogna et al. (2019, 2021). By applying a post-processing step, those models can be used to predict atomic line emission (e.g. Ercolano & Owen 2016) or for direct comparison to observations (e.g. Weber et al. 2020). In addition to the photoevaporative models, Weber et al. (2020) used semi-analytic MHD wind models and suggested that the multi-component [OI] spectral line profiles could be explained by a combination of thermally driven (LVC) and MHD-driven (HVC) wind components.

More consistent wind models, coupling the dynamics, radiation physics, and thermochemistry for photoevaporative winds, are presented by Wang & Goodman (2017) and Nakatani et al. (2018a,b), who show that molecular hydrogen can survive in several of their models (i.e. depending on the included heating sources). The first theoretical models for MHD outflows and winds that include chemistry were presented by Panoglou et al. (2012). These authors follow the chemical evolution along individual wind stream lines, and one of their main findings is that, even in the T Tauri stage, molecular hydrogen can survive in the wind regions. Full MHD disk-wind models including thermochemistry were built by Wang et al. (2019) and Gressel et al. (2020). These latter models indicate that both magnetic and thermal effects (i.e. radiation) can play an important role in slow disk winds. Such self-consistent disk-wind models are computationally expensive and therefore a direct comparison to observations is difficult; only Gressel et al. (2020) presented synthetic observables for their models, but these authors focused on far-infrared (FIR) and sub-millimetre (submm) line emission.

In this work, we aim to interpret observations of spectral lines o–H₂1–0 S(1) at 2.12 µm and the [OI] ¹D₂–³P₂ at 0.63 µm in the context of photoevaporative disk-wind models. Following the approach of Weber et al. (2020) for modelling the atomic line tracers, we use existing radiation-hydrodynamic photoevaporative disk-wind models and post-process them with a radiation thermo-chemical disk code that includes molecular chemistry. This allows us to model the molecular and atomic line emission simultaneously for [OI] 0.63 µm and o–H₂ 2.12 µm. With this efficient approach, we aim to interpret the observations of Gangi et al. (2020), focusing on the line kinematics and the physical origin of the line emission.

We first explain our modelling approach in Sect. 2. In Sect. 3, we present our results and compare the models to the observational data with a focus on line kinematics. In Sect. 4, we discuss our interpretation of the line kinematics in the context of our models, and simpler approaches and possible next steps for theoretical disk-wind modelling. In Sect. 5, we summarise our results and main conclusions.

2 Method

2.1 Radiation thermochemical modelling

For the physical disk-wind density structure and the velocity field we use the 2D EUV/X-ray (XUV) photoevaporative disk-wind models presented in Weber et al. (2020). For these models, we use the PLUTO (Mignone et al. 2007) hydrodynamic code, extended with a parameterised treatment for the gas temperature (Picogna et al. 2019). The parameterised treatment is based on detailed gas photoionisation and radiative-transfer disk models built with MOCASSIN (e.g. Ercolano et al. 2003, 2005, 2008b). As shown in Weber et al. (2020), those models are in good agreement with observed kinematic signatures of the [OI] 0.63 µm spectral line; however, they do not include the chemical modelling of molecules such as molecular hydrogen. Therefore, in addition to these models, we use the radiation thermochemical disk code PRODIMO (PROtoplanetary DIsk MOdel¹, Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Woitke et al. 2016). PRODIMO uses wavelength-dependent continuum radiative transfer to calculate the dust temperature and the radiation field in the disk. Furthermore, PRODIMO solves consistently for the gas temperature and the chemical abundances and has a module to produce synthetic observables such as spectral line profiles.

For this work, we added an interface to PRODIMO that allows us to use the 2D gas density and velocity structure from hydrodynamic models as input. Figure 1 shows the gas density and velocity structure for the disk and wind component as used in PRODIMO. This structure remains the same for all the models presented in this work (i.e. as the X-ray luminosity is always the same). For the dust structure, we adopt the same dust properties as Weber et al. (2020) and Ercolano et al. (2008a); an Mathis, Rumpl, Nordsieck (MRN) dust size distribution (Mathis et al. 1977) and no settling (see also Sect. 2.4 and Fig. A.1). Besides this interface, it was not necessary to add any other new physical or chemical processes to the PRODIMO model. To validate the approach, we compared the results from PRODIMO to the results of Weber et al. (2020), who used MOCASSIN for postprocessing, and find reasonable agreement for both the thermal structure and the synthetic [OI] 0.63 µm spectral line predictions. The details of this comparison are described in Appendix B. We note that we did not attempt to benchmark PRODIMO and MOCASSIN, but simply compared the results of both codes for the same physical disk+wind structure. Although there are certain differences in the results, the general agreement is quite remarkable, considering that both codes were developed independently and use different methods for the radiative transfer and thermal-balance modelling for example.

We note that we neglect dynamics when modelling the chemistry, which means that we likely underestimate the amount of molecular hydrogen that might be transported from the disk into the wind region. However, as we model a T Tauri star with a significant far-ultraviolet (FUV) radiation field due to ongoing accretion, it is unlikely that the molecular hydrogen survives for long periods in the disk wind (see also Panoglou et al. 2012). Nevertheless, because of this simplifying assumption, the results for the o–H₂ 2.12 µm emission likely represent lower limits, as only the molecular hydrogen that is formed and survives in the wind region is modelled (for more details see Sect. 3.1). In future work, we will consider the dynamics and thermochemical processes consistently. However, such models are complex and computational expensive. The approach used in the present work is therefore a first step towards more self-consistent models and allows efficient direct comparison of existing hydrodynamic models to observations (see also Sect. 4.3).

Fig. 1 Two-dimensional disk and wind gas density structure and velocity field (black arrows). This structure is the same in all the presented models. The red solid contour shows where the radial visual extinction Av,rad is equal to unity; here, we refer to this layer as the disk surface. The white dashed contours correspond to the values shown in the colour bar.

2.2 Chemical modelling

For the chemical modelling, we use a chemical network with 100 gas and ice phase species and 1288 chemical reactions, including gas-phase chemistry (i.e. photo-reactions); freeze-out of atoms and molecules; photo-, thermal-, and cosmic-ray desorption of ices; X-ray chemistry (see Aresu et al. 2011; Meijerink et al. 2012); H₂ formation on grains; and excited H₂ chemistry. This chemical network is described in detail in Kamp et al. (2017), and was used for example to model multi-wavelength observational data of protoplanetary disks (DIsc ANAlysis project, DIANA²). Also, to the network of Kamp et al. (2017), we add collisional ionisation of H by electrons using the data from Janev & Smith (1993).

For this work, we focus on molecular hydrogen and therefore summarise the main chemical processes for the H₂ chemistry, in particular the formation and destruction pathways used in our model. We note that PRODIMO is a flexible code and other options and chemical rates can be used. Here, we only describe the options adopted in this work.

For H₂ formation on dust grains, we follow the analytical model of Cazaux & Tielens (2002, 2004, 2010); furthermore, we include gas-phase formation of H₂ via H⁺ and H⁻ (for details see Woitke et al. 2009; Thi et al. 2020). The main H₂ destruction reaction relevant in the wind region is photo-dissociation of H₂. The rate is calculated using the frequency-dependent FUV radiation field from the radiative transfer modelling and detailed photo-dissociation cross-sections (Allison & Dalgarno 1969; van Dishoeck 1988) and by correcting for self-shielding using the approximation of Draine & Bertoldi (1996) (for details see Woitke et al. 2009). We also tested the H₂ self-shielding function of Wolcott-Green & Haiman (2011), but having found no significant impact on our results (i.e. line profiles), we decided to remain with the prescription of Draine & Bertoldi (1996). Furthermore, we include excited molecular hydrogen chemistry as described in Kamp et al. (2017). All other gas-phase chemistry reactions involving H₂ are taken from the UMIST 2012 Database (McElroy et al. 2013).

2.3 Line excitation and synthetic observables

For the calculation of the non-local thermodynamic equilibrium (NLTE) level populations and the line cooling for O and H₂, the escape probability method, including UV pumping, is used in PRODIMO (Woitke et al. 2009). Also included are fluorescent UV pumping and chemical pumping by OH photo-dissociation for O (Acke et al. 2005), and pumping by its formation on dust grains for H₂. The details of those pumping mechanisms are described in detail in Appendix A of Woitke et al. (2011).

The collisional and radiative data for O are from the National Institute of Standards and Technology (NIST) atomic spectroscopic database (Ralchenko 2009), the Leiden Atomic and Molecular Database (LAMDA Schöier et al. 2005), Krems et al. (2006), and Störzer & Hollenbach (2000); including collisions with o–H₂, p–H₂, H, H⁺, and e⁻. For H₂, the data are from Wolniewicz et al. (1998), Wrathmall et al. (2007), and Le Bourlot et al. (1999), including collisions with o–H₂, p–H₂, H, and He. The H₂ ortho/para ratio is assumed to be at thermal equilibrium according to the local gas temperature; for further details, see Woitke et al. (2009) and Woitke et al. (2011).

For a proper comparison to observations, we use the line radiative transfer module of PRODIMO (Woitke et al. 2011, Appendix A) to calculate the spectral line profiles for the [OI] 0.63 µm and o–H₂ 2.12 µm lines at various disk inclinations i. For the local line width, we assume , where v_therm is the thermal broadening and v_turb the turbulent broadening, which is fixed to 0.15km s⁻¹ (see also Thi et al. 2019). The modelled line profiles are then convolved to the spectral resolution of the observations (R ≈ 50 000). Furthermore, we use the routine for the line profile decomposition (i.e. fitting Gaussian functions to the line profiles) of Weber et al. (2020), which follows the fitting approach used for the observational data (Banzatti et al. 2019; Gangi et al. 2020). Including the Gaussian fitting for the synthetic observables allows a direct comparison between the derived quantities, such as the full width at half maximum (FWHM), and the velocity peak location v_p of the Gaussian component(s) (see Sect. 3.3).

Fig. 2 Stellar input spectra for the three different accretion luminosities Lacc.

2.4 Model series

Apart from the underlying physical disk structure, we use the same stellar properties as Weber et al. (2020). The disk is irradiated by a M_* = 0.7 M_⊙ star using the XUV (EUV and X-rays) spectrum presented in Ercolano et al. (2008a, 2009). The X-ray luminosity is fixed to L_X = 2 × 10³⁰ ergs⁻¹, whereas for the UV radiation field we consider three different accretion luminosities of L_acc = 2.6 × 10⁻², 0.31, and 1 L_⊙ assuming a black-body spectrum with T_eff = 12 000 K (see Fig. 2). We note that we neglect any screening of the XUV emission from the star by possible accretion flows or inner MHD winds, which might be present in case of strong accretion (i.e. relevant for the PE-3 model). For screening columns N_H ≳ 10²²cm⁻², XUV photoevaporation might be completely suppressed, assuming that the screening covers all solid angles at all times (Ercolano et al. 2009). However, this inconsistency in our models does not affect our main results significantly. In any case, and as shown below, the high-accretion models are not expected to produce strong wind signatures for H₂.

The amount of dust in the wind is relevant for H₂ formation and for shielding of H₂ from photo-dissociation (additionally to the self-shielding), but it is not well constrained from observations. To simulate dust entrainment in the wind, we adopt a very rough approximation and simply reduce the amount of dust in the wind region by a factor of 10 (i.e. gas-to-dust mass ratio g/d = 1000). This is based on the dust entrainment model for photoevaporative winds of Franz et al. (2020, 2022) (see also Booth & Clarke 2021; Hutchison & Clarke 2021; Rodenkirch & Dullemond 2022). We want to emphasise that our approach here is not entirely realistic and should only give a first indication as to the relevance of the amount of dust in the wind for the interpretation of the o–H₂ 2.12 µm spectral line. For this new wind dust structure, we run models again for all three accretion luminosities. The model series parameters and names are summarised in Table 1.

As an example, we show physical quantities in Fig. 3 as calculated by PRODIMO for the PE-2 and the PE-2 gd models. Shown are the gas temperature T_gas, dust temperature T_dust, and the FUV radiation field χ (in units of the Draine field Draine & Bertoldi 1996; Woitke et al. 2009), which is most relevant for H₂ photo-dissociation.

3 Results

3.1 H₂ in the wind and the line-emitting regions

Figure 4 shows the main line-emitting regions for [OI] 0.63 µm and o–H₂ 2.12 µm for a face-on inclination. Appendix C (Table C.1) also summarises the average physical properties within those line-emitting regions.

Figure 4 shows that the o–H₂ 2.12 µm line emits from close to the disk surface (the disk-wind interface), but, depending on the accretion luminosity, a significant contribution of the line flux can also come directly from the wind region. This implies that in those models a significant amount of molecular hydrogen can survive in the photoevaporative wind region. We emphasise that our model represents a lower limit to the H₂ abundance in the wind as it neglects dynamical effects, and hence no H₂ is transported from the disk into the wind. Although H₂ would be quickly dissociated and might not survive for a long time (i.e. the photo-dissociation timescale is shorter than the flow timescale), the dynamics can only enhance the amount of H₂ in the wind (see Panoglou et al. 2012) and therefore also the contribution to the total line flux from the wind region.

For the models with less dust in the wind, the situation is similar. However, due to the lower amount of dust, FUV radiation can penetrate deeper into the wind and disk and therefore photo-dissociation of H₂ becomes more efficient. Consequently, the main emitting region of o–H₂ 2.12 µm moves towards the disk surface and emission from the wind region becomes less significant. However, for the PE-1 gd model (lowest accretion luminosity), significant emission is still coming from the wind regions (Fig. 4, bottom left panel).

Compared to [OI] 0.63 µm, o–H₂ 2.12 µm is emitted from deeper layers in the disk and from farther out (up to r ≈ 30 au) in the disk. This is because H₂ can only survive in regions where it is efficiently shielded from photo-dissociation. Although there can be some overlap in the line-emitting regions of the two lines in the radial direction, the two lines rather trace distinct regions of the disk. Furthermore, the atomic oxygen line is excited at higher temperatures (several thousand K; see Table C.1 and Weber et al. 2020) than o–H₂ 2.12µm, which mainly comes from regions with T_gas ≲ 1000 K. Whether or not those distinct emitting regions can be traced with the spatially unresolved observations is discussed in Sects. 3.3 and 4.

Fig. 3 Two-dimensional gas temperature Tgas, dust temperature Tdust, and FUV radiation field χ (in units of the Draine field Draine & Bertoldi 1996; Woitke et al. 2009) for the PE-2 (top row) and PE-2 gd (bottom row) models. The red solid contour shows where the radial visual extinction AV,rad is equal to unity, which roughly separates the wind region and the rotating disk. The white dashed contours correspond to the values shown in the colour bars.

Fig. 4 Main emitting region for the [OI] 0.63 µm (orange boxes) and o–H2 2.12 µm (red boxes) spectral lines. The columns show the three models with varying Lacc (increasing from left to right). The top row shows the models with g/d = 100 and the bottom row with g/d = 1000 in the wind. The coloured contours show the molecular hydrogen number density . The borders of the line emission boxes indicate where the total emission reaches 15% and 75% in the radial (integrated inside-out) and vertical (integrated from top to bottom) directions. The solid white line shows where the radial visual extinction AV,rad reaches unity, and the white dotted line Tgas = 1000 K. The grey arrows indicate the 2D velocity field (vertical and radial components of the wind).

3.2 Line luminosities

In Fig. 5, we compare the modelled line luminosities for [OI] 0.63 µm and o–H₂ 2.12 µm to the observations. We mainly use the observational data from Gangi et al. (2020), but also discuss further data from the literature for targets where line flux measurements for both lines are available (see Appendix D).

In addition to the wind models presented here, we include modelling results from the DIANA (DIsc ANAlysis) project. For this project, multi-wavelength observational data (continuum and line fluxes) for several individual disks were modelled using PRODIMO and the results are published in Woitke et al. (2019). Those disk models do not include any wind component, but can still be used to decipher the approximate impact of stellar properties and disk structure on the line fluxes. The [OI] 0.63 µm and o–H₂ 2.12 µm were not included in the fitting procedure at that time, but predictions for the line fluxes were included in the models and those published values are directly used here.

From Fig. 5, we see that the models populate the o–H₂ 2.12 µm versus [OI] 0.63 µm line luminosity space reasonably well, but also that both the wind and the DIANA models tend to underpredict the o–H₂ 2.12 µm line luminosities. We discuss this in more detail in Sect. 4.2. For [OI] 0.63 µm, all models are within the range of the observations, except for one outlier in the DIANA models. This is simply because this object, ET Cha, is a very small (r < 10 au) disk (Woitke et al. 2011; Ginski et al. 2020).

The agreement of the DIANA models with the observations for three out of four targets for which observational data are available is within a factor of two to three for both lines; only for one target is the o–H₂ 2.12 µm luminosity significantly lower (by about two orders of magnitude). overall, this is an encouraging outcome considering that the models were not optimised to fit the o–H₂ 2.12 µm and [OI] 0.63 µm line luminosities. The three DIANA models with the highest [OI] 0.63 µm and o–H₂ 2.12 µm luminosities are Herbig Ae/Be stars with L_* > 10 L_⊙ and these therefore also have a high luminosity in the FUV range. observational data are not available for these targets. The Gangi et al. (2020) sample includes only one Herbig Ae/Be star, MWC 480, for which a DIANA model also exits. For MWC 480, the DIANA model also agrees within a factor of three with the [OI] 0.63 µm luminosity, but for o–H₂ 2.12 µm only an observational upper limit exists, which is consistent with the model prediction.

There is also a tendency towards slightly too low o–H₂ 2.12µmluminosities for the photoevaporative wind models with respect to the observational data. However, at least for the high-accretion luminosity and/or the reduced dust-to-gas mass ratios in the wind, the predicted o–H₂ 2.12µmluminosities are similar to the line luminosities observed by Gangi et al. (2020). If the additional literature values are also taken into account, the wind models populate the o–H₂ 2.12 µm versus [OI] 0.63 µm line luminosities remarkably well. We note that unlike the DIANA models, the disk-wind models use only one fiducial disk-wind structure and only the accretion luminosity is varied. The disk structure, but also the stellar properties (e.g. higher L_acc leads to higher line luminosities due to additional heating), can have a significant impact on the line luminosities, which is indicated by the significant scatter of the DIANA models seen in Fig. 5.

Our results are also consistent with the thermo-chemical disk models of Nomura & Millar (2005) and Nomura et al. (2007). Nomura et al. (2007) studied the impact of X-ray and FUV radiation and dust size distribution on H₂ line fluxes and find line luminosities in the range of 0.001–0.2 × 10⁻⁵ L_⊙ for o–H₂ 2.12 µm. This is consistent with our results, which show a range of 0.003–0.3 × 10⁻⁵ L_⊙ (including both wind and DIANA models; see Fig. 5). We note that Nomura et al. (2007) also studied the impact of H₂ pumping by X-ray radiation, which is not included in our model. These authors find that as long as FUV radiation is present, X-ray pumping is not a dominating factor. However, they argue that for disks around stars with strong X-rays and weak UV, it might be possible to observe H₂ emission from cooler regions excited by X-ray pumping, but more detailed modelling of this process is required.

The comparison of the models to the data shows that the former can reproduce reasonably well the observed line luminosities for both o–H₂ 2.12 µm and [OI] 0.63 µm. However, we also note that our fiducial photoevaporative disk-wind model is not a good representation of certain targets. one example is DG Tau, as already shown by Weber et al. (2020) for [OI] 0.63 µm (see also Cabrit et al. 1999). DG Tau shows [OI] 0.63 µm blueshifted emission (with respect to the systemic velocity) at velocities < −30 km s⁻¹, which cannot be produced by the photoevaporative disk-wind models presented here. Additionally, o–H₂ 2.12 µm shows strongly blueshifted emission at velocities of 10–20 km s⁻¹ (Gangi et al. 2020) and DG Tau shows one of the highest o–H₂ 2.12 µm luminosities (see Fig. 5). Also, other observational data show indications of both MHD-driven and photoevaporative flows, although shocks might also play a role (e.g. Agra-Amboage et al. 2014; Güdel et al. 2018). Furthermore, DG Tau is strongly variable in the near- and mid-infrared (e.g. Varga et al. 2017; Gangi et al. 2020). This complex environment might be the reason for the strong o–H₂ 2.12 µm luminosities, which cannot be matched by our models orbythe models of Nomura et al. (2007). Such ascenario might also be responsible for high o–H₂ 2.12 µm line fluxes observed in other targets.

Nevertheless, both types of models, the DIANA and the disk-wind models, predict [OI] 0.63 µmand o–H₂ 2.12 µm line luminosities within the observed range. This indicates that the physical (i.e. heating and cooling), chemical, and line excitation model used in PRODIMO is well suited to modelling the o–H₂ 2.12 µm and [OI] 0.63 µm spectral lines, including the line kinematics.

Fig. 5 o–H2 2.12 µm versus [OI] 0.63 µm line luminosities. The dark grey points with error bars show the observations from Gangi et al. (2020, 2021), and the light-grey points show observational data collected from the literature (see Appendix D for details). The orange symbols show the photoevaporative disk-wind models (inclination i = 40°) with varying accretion luminosities (lowest(highest) line luminosities are for the lowest(highest) accretion luminosity). The square symbols are for the models with reduced dust-to-gas mass ratio in the wind. The brown diamonds are the results from the DIANA project (PRODIMO models without a wind component). We mark DG Tau as an example of a target that likely cannot be modelled with a pure photoevaporative wind. We also mark the two independent observations for GM Aur (see Appendix D).

Fig. 6 Example of the line-fitting results for the PE-1 model. The top figure shows the results for a spectral resolution R = 50 000 (i.e. similar to the observations). The bottom figure shows the same model, but the line profiles were convolved to a spectral resolution of R = 100 000 (see Sects. 4.1.1 and 4.1.4). Each individual panel shows the line profile with the model resolution (blue), convolved to the target spectral resolution (black) and the Gaussian fit (the NLVC) to the convolved profile (red solid line). In case of multiple Gaussians (dashed and dotted lines), the orange dashed line indicates the one chosen as the NLVC. This component is used to determine vp and the FWHM (also given in each panel). The fitting results for all other models are shown in Appendix F.

3.3 Line kinematics

In this section, we focus on the comparison of the kinematic properties derived directly from the observed line profiles to the model results. Those kinematic quantities are the shift of the peak location v_p and the FWHM or half-width at half-maximum (HWHM), as presented in Gangi et al. (2020). To derive v_p and the FWHM for the modelled line profiles, we follow the approach of Gangi et al. (2020) and fit the modelled line profiles with Gaussian components (see Sect. 2.3), focusing on the blueshifted NLVC (|v_p| < 30 km s⁻¹). Figure 6 shows the results of this fitting procedure for the PE-1 model as an example. Shown are the results for two spectral resolutions R = 50000 and R = 100 000. The first one is similar to the observations and those results for v_p and the FWHM are used for the comparison of the models to the observations. The high-resolution results are used to investigate the impact of the limited spectral resolution of the observations on the derived line kinematics and are discussed in Sects. 4.1.1 and 4.1.4, but we show them here for ease of comparison.

As already shown in Weber et al. (2020), photoevaporative wind models are unlikely to be the origin of the observed HVCs seen for [OI] 0.63 µm. For o–H₂ 2.12 µm, Gangi et al. (2020) only considered one component, the NLVC, as none of the targets show a HVC or a clear BLVC. For comparison with the models, we only use targets from the observational sample that have a detection and derived kinematic quantities for the NLVC for both lines. From those 15 targets, 8 show v_p > 0 km s⁻¹ and only 7 have a blueshifted NLVC in both lines. For the [OI] 0.63 µmline profiles of these latter 7 targets, two cases show a clear HVC, 2 show a BLVC, and 1 case shows no other component. For the remaining 2 cases, the high-velocity part of the profiles could not be fitted by a Gaussian (for details see Appendix E).

As seen in Fig. 6 (see also Appendix F), we find blueshifted components in both lines for almost all of our modelled profiles. However, especially for o–H₂ 2.12 µm, the derived v_p is often close to zero (see also Sect. 3.3.1). Figure 6 also shows that at intermediate inclinations there are several cases where the profiles are not well fitted by a Gaussian (or multiple Gaussians). The reason is that, at those inclinations, the shape of the line profiles is strongly affected by the complex wind velocity field (see also Weber et al. 2020), the absorption of the redshifted emission by the dust, and, in the case of o–H₂ 2.12 µm, also by some self-absorption. We note that this self-absorption effect is limited to a narrow range of inclinations of about i ≈ 60° ± 5° in all our models. However, the inclination at which this self-absorption effect might be seen depends on the actual disk and wind structure. The issue of non-Gaussian components in the line profiles becomes especially apparent in the models with high-spectral resolution (see Fig. 6). This is discussed further in Sects. 4.1.1 and 4.1.4.

3.3.1 Peak velocity

In Fig. 7, we compare the model results for the peak velocity v_p to the observations. The figure shows that the models are largely consistent with the data. One clear exception is the target DO Tau, which shows v_p < −10 km s⁻¹ for both spectral lines, whereas the largest shifts predicted in the models are v_p ≈ −6 km s⁻¹. This is expected, as photoevaporative disk-wind models do not predict such high velocities for [OI] 0.63 µm(e.g. Weber et al. 2020). This implies that photo-evaporation is not the main physical mechanism driving the wind or outflow in DO Tau. However, DO Tau is likely a special case as there is indication for a recent stellar encounter (Winter et al. 2018) and recent observations clearly show a complex environment (Huang et al. 2022).

Another discrepancy is that all models show v_p < 0 km s⁻¹ for both spectral lines, whereas for the majority of the observed targets v_p > 0 km s⁻¹. However, this does not necessarily mean that there is no observed signature of a wind in those targets. In particular, the observed line profiles for [OI] 0.63 µm still show blueshifted components at high velocities, which are likely signatures of winds or outflows. However, for o–H₂ 2.12 µm, it is often not clear whether or not there is a blueshifted component at all, which can rather indicate that there is no significant amount of H₂ in the wind or that the data quality is not sufficient to clearly identify the NLVC; the latter is likely also true for [OI] 0.63 µm. In the pure photoevaporative disk–wind models it is likely easier to identify the NLVC for the [OI] 0.63 µm compared to the observations, as the models do not include a high-velocity component. For o–H₂ 2.12 µm, we suspect that limitations of the data are the main cause for the group of targets with v_p > 0 km s⁻¹ as inspection of the individual targets does not clearly reveal redshifted emission. For example, the o–H₂ 2.12 µm line profile for BP Tau, which is the target with the highest v_p, shows a high level of noise. For most of the other targets with observed v_p > 0 km s⁻¹ the error bars indicate that the observations are also consistent with v_p ≈ 0 km s⁻¹, that is, no observed wind component or a slow wind. Considering those arguments, in what follows, we focus on the targets with v_p < 0 km s⁻¹ for our comparisons to the observational data.

For the models, 7 out of 30 show a v_p < −3 km s⁻¹ for o–H₂ 2.12 µm (we assume σ(v_p) ≈ 3 km s⁻¹ as a typical error in the observations). As seen from Fig. 7, this is the case for models with lower L_acc in particular (i.e. the PE-1 model). This is because photo-dissociation of H₂ is less efficient because of the lower flux of FUV photons. Also, all those models have an inclination of i ≲ 60°, indicating that the detection of disk winds is difficult for highly inclined disks, at least with spatially unresolved observations.

As reported in Gangi et al. 2020, and as clearly seen in Fig. 7, the observational data show that v_p(H₂) ≈ v_p([OI]). This is also true for the models, although there might be a slight trend towards v_p([OI]) < v_p(H₂) for the models with higher accretion luminosities (i.e. PE-2 and PE-3). In the framework of photoevaporative winds, such a trend would indicate that [OI] 0.63 µm traces higher wind velocities than o–H₂ 2.12 µm, and that o–H₂ 2.12 µm is emitted from smaller heights (i.e. closer to the disk surface) than [OI] 0.63 µm, as the wind velocity increases with height. This is indeed the case for our photoevaporative wind models, as seen in Fig. 4 and discussed in Sect. 3.1.

According to Gangi et al. (2020), there is no correlation between v_p and disk inclination i. As we present only one physical wind–disk structure, the model series is not well suited to investigating this non-correlation. Nevertheless, it is still interesting to check the existing models for a correlation with i. As seen in Fig. 8, there is a slight trend of increasing v_p with inclination within each group of models with the same L_acc, that is, if the highest inclination is excluded (those models always show v_p ≈ 0 km s⁻¹), where the trend is stronger for [OI] 0.63 µm than for o–H₂ 2.12 µm. Considering the expected errors in real observations, such a correlation would be hard to identify (see also Gangi et al. 2020). Furthermore, Fig. 8 also shows that both L_acc and the gas-to-dust mass ratio have a stronger (or at least similar) impact on v_p than the inclination itself. Although we present only one physical structure, our model series indicates that it is challenging to observe a clear correlation of v_p with inclination and that other properties of the targets, such as the accretion luminosity, also have to be considered.

Fig. 7 Velocity peak location vp of o–H2 2.12 µm vs vp [OI] 0.63 µm. Filled circles with error bars show the observed values; the black colour is for targets with vp < 0 km s−1 (i.e. blue-shifted wind component) for both the o–H2 2.12 µm and [OI] 0.63 µm; light-grey marks targets with vp > 0 km s−1 (i.e. no clear indication of a disk wind). The coloured symbols are for the models. Shown are models for the three different accretion luminosities and with a reduced dust-to-gas mass ratio in the wind (coloured diamonds). The size of the coloured symbols scales with the inclination; smallest symbols are for i = 0° largest are for i = 80°. The thick dashed line indicates vpH2 = vp[OI] and the light grey stripe indicates vpH2 = vp[OI] ± 3 km s−1. Top panel: all data points; bottom panel: zoom-in (brown box) showing all the models and the observational data for targets with clear wind signatures, excluding DO Tau.

Fig. 8 Velocity peak location vp versus disk inclination for the [OI] 0.63 µm (top panel) and o–H2 2.12 µm (bottom panel). Filled circles with error bars show the observed values; the black colour is for targets with vp < 0 km s−1 (i.e. blueshifted wind component) for both the o–H2 2.12 µm and [OI] 0.63 µm; light-grey marks targets with vp > 0 km s−1 (i.e. no clear indication of a disk wind). The coloured symbols are for the models. Shown are the same models as in Fig. 7.

3.3.2 Full width at half maximum

The FWHM or HWHM of a spectral line can be used to determine the maximum radius of the emitting region of the line. Gangi et al. (2020) used the relation (1)

where G is the gravitational constant and M_* is the stellar mass. This relation assumes a thin disk in Keplerian rotation. This approximation does not fully represent the physical conditions for the spectral lines studied here, as the lines are not emitted in the midplane or in a thin layer, and because the spectral line profiles are strongly affected by the wind velocity (see Weber et al. 2020). Nevertheless, it is a rough indicator of the maximum emitting radius R_K and it can be directly applied to observational data. To allow for a direct comparison between the models and the data of Gangi et al. (2020), we also use Eq. (1) and measure the HWHM of the modelled spectral lines from the fitted Gaussian of the respective NLVC.

In Fig. 9, we compare the derived HWHM scaled by the stellar mass as a function of disk inclination to the results from Gangi et al. (2020). The model results are in good agreement with the data for both [OI] 0.63 µm and o–H₂ 2.12 µm. The measured HWHM from the modelled [OI] 0.63 µm profiles show a peak at medium inclination (i = 40°), but the scatter is also significant, especially if the models with lower dust content in the wind are also considered. Such a trend is not really visible in the data, but the error bars are large and only very few measurements for disks close to face-on or edge-on exist. For o–H₂ 2.12 µm, the HWHM does not show a dependence on disk inclination either in the models or in the available observational data.

Figure 9 also indicates that, for the models, the HWHM for [OI] 0.63 µm is systematically larger than for o–H₂ 2.12 µm in the range of sin(i) = 0.2–0.9. For those 16 models, the ratio HWHM_[OI] / ranges from ≈ 0.9 to ≈ 2.4, except for one outlier with a ratio of four and two models with a ratio of smaller than unity. The situation is similar for the subset of the observational data with v_p < 0 km s⁻¹ in both lines (all sin(i) are within 0.2 and 0.9). For those targets, the ratio HWHM_[oi]/ is in the range of ≈ 0.6 to ≈ 2.5, and 2 out of 7 have a ratio smaller than unity. We further discuss this in Sect. 4.1.

4 Discussion

4.1 Interpretation of observed line kinematics

As discussed in Gangi et al. (2020), the observational data show that the line kinematics for the two spectral lines [OI] 0.63 µm and o–H₂ 2.12 µm are similar. Gangi et al. (2020) argue that this might indicate that the lines trace similar physical regions of a disk wind, and that such a scenario is more consistent with a centrifugal MHD-driven wind as studied in Panoglou et al. (2012). Panoglou et al. (2012) follow the thermochemical evolution (including FUV and X-ray irradiation) along wind streamlines derived from MHD models. These authors found that the region where H₂ can exist in the wind evolves with the evolutionary stage (i.e. Class 0/I/II) of the target; in particular, wind temperatures become higher and the shielding of H₂ less efficient because of the stronger FUV radiation field impinging on the disk at later stages. We note that Panoglou et al. (2012) did not produce synthetic observables for their models, and therefore a direct comparison to observational data is not possible. Nevertheless, their results for the Class II/T Tauri stage are in general agreement with our results, as we also find that molecular hydrogen only survives close to the disk surface.

In this work, we show that observed line kinematics for [OI] 0.63 µm and the o–H₂ 2.12 µm are consistent with a photoevaporative wind. This indicates that both scenarios, that is, MHD- and photoevaporative driven winds, might be consistent with the kinematics derived from currently available o–H₂ 2.12 µm observations. The main argument from Gangi et al. (2020) for MHD winds is that the observational data seem to indicate that the [OI] 0.63 µm and o–H₂ 2.12 µm are tracing similar regions of the wind. This is in contradiction to our modelling results, which indicate that [OI] 0.63 µm comes from regions closer to the star and higher up in the wind compared to the emitting region of o–H₂ 2.12 µm (see Fig. 4). Nevertheless, the kinematic quantities derived from our models are roughly consistent with the data. There may be several reasons for this: the spectral resolution of the observational data is not high enough to discriminate the two emitting regions; the thin-disk approximation to determine the line emitting region is perhaps too simple, and the physical properties (i.e. the wind velocity) are similar in both line-emitting regions (i.e. because of different heights). In the following sections, we discuss these possibilities in more detail, and separately for the peak velocity and FWHM.

Fig. 9 HWHM (left panel: [OI], right panel H2) scaled by the square root of the stellar mass as a function of inclination. The coloured symbols show the models; the gd models (diamonds) are slightly shifted along the inclination axis for clarity. The black (grey) symbols show the observations, where the black coloured symbols mark the targets that show blueshifted peaks in the respective line profile. The solid lines show the maximum emission radius RK derived for the given HWHM and inclination assuming Keplerian rotation for a thin disk (see Eq. (1)).

4.1.1 Impact of spectral resolution on peak velocity

To determine the shift in the peak velocity v_p, the spectral resolution is most relevant. As seen from the data, many targets are around v_p ≈ 0 km s⁻¹ with error bars as big as ±3 km s⁻¹, and therefore it is unclear whether there is a wind signature or not.

Using the models, we have the possibility to study the impact of the spectral resolution. Figure 10 shows the same observational data as in Fig. 7 but for the models we now assume a spectral resolution of R = 100000, which is two times higher than for the data. The general picture for the models does not change significantly. As can be seen in Fig. 10, some models now show v_p ≲ −3 km s⁻¹ for one line but v_p ≈ 0 km s⁻¹ for the other, which increases the scatter in the v_pH₂ versus v_p[OI] plane. This is a consequence of the fitting procedure. For those models, the line profiles are now fitted by two (or more) Gaussian components, instead of one in the low-spectral-resolution case. Especially for the PE-1 model, we find that out of the ten line profiles, four are fitted by multiple Gaussians for R = 100000 (see Fig. 6). For example, for i = 0°, the PE-1 model shows a narrow peak and a significantly blueshifted component due to the wind (first column in Fig. 6). This profile was still fitted by a single Gaussian in the low-resolution case, but two Gaussians were identified in the high-resolution case. This makes the identification of the NLVC in the model ambiguous. A similar issue arises at higher inclinations for o–H₂ 2.12 µm because of self-absorption (see Sect. 3.3).

We note that for the models with R = 50 000, the identification of the NLVC is straightforward as all modelled profiles are fitted with a single Gaussian except for one case. In the PE-1 gd model, the line profile for the o–H₂ 2.12 µm at i = 60° (see Fig. F.2) requires a two-Gaussian fit, but this is due to the self-absorption effect. For the high-resolution models, 9 out of the 60 modelled line profiles resulted in a multi-Gaussian fit. Most of these are for the PE-1 and PE-1 gd models (7 out of 20), where we see the strongest wind signatures in the line profiles. Although an unambiguous identification of the NLVC is possible for the majority of our models by applying the Gaussian fitting procedure, our results also indicate that a clear identification of the NLVC might not be feasible for higher spectral resolution observations. For high-resolution observational data, the identification of the NLVC will likely be even more ambiguous, in particular for the cases where broad and high-velocity components are also present, as they might not be well represented by Gaussian components at all. On the other hand, such high-resolution observations might allow more detailed studies of the line shape without Gaussian fitting; for example, the blueshift of the profile could be determined simply by measuring the velocity at the actual peak of the line profile.

For this experiment, we did not adapt the fitting procedure, and the fitting mechanism selects the component with v_p closer to zero. This is rather an issue of the definition of the NLVC and the Gaussian-fitting procedure, which might need to be revised for higher spectral resolution data. However, it does not affect the conclusion that, in the models, both lines still show signatures of blueshifted emission. Nevertheless, this exercise shows that an increase in the spectral resolution by a factor of two does not change the general picture that, for both the models and the data, the derived v_p is similar for both lines in many cases.

Fig. 10 Same as the bottom panel of Fig. 7, but for the modelled line profiles a spectral resolution of R = 100 000 is used.

Fig. 11 Ratio of the FWHM of [OI] 0.63 µm and o–H2 2.12 µm. Shown are the same models as used for Fig. 9. Additionally, we show the model series called Keplerian (crosses), where we assume that the velocity field in the disk and wind is purely Keplerian (i.e. thin-disk approximation). The black data points with error bars mark the targets with blueshifted vp in both lines.

4.1.2 Interpreting the observed line width

Figure 11 shows the ratio of the FWHM of the two studied lines to compare the measured line widths. The model results show that the FWHM for [OI] 0.63 µm is larger than that for o–H₂ 2.12 µm for 20° < i < 60°. This is also the case in the data but is not as evident there. In particular, at sin(i) ≈ 0.6, there are two targets with a ratio close to one: XZ Tau and DG Tau. XZ Tau has a ramp-like [OI] 0.63 µm line profile and therefore has some uncertainty in the fit (Gangi et al. 2020). DG Tau is likely not well represented by our model, as already discussed in Sect. 3.2. The one target with a ratio < 1 at sin(i) ≈ 0.8 is GM Aur, which is a transitional disk with a large inner hole and rings (see e.g. Huang et al. 2020), which might be the reason for the observed narrow and similar line widths for the [OI] 0.63 µm and o–H₂ 2.12 µm (see also discussion in Gangi et al. 2020); however, at those high inclinations, our model also predicts ratios close to unity.

As already mentioned, a direct comparison of individual targets to our models is not feasible, as here we only present a fiducial model. Nevertheless, Fig. 11 clearly shows that the photoevaporative wind models are consistent with the observationally derived kinematics, even though in our models the [OI] 0.63 µm and o–H₂ 2.12 µm lines are not emitted from the same region in the wind or disk. In the following, we want to focus on this aspect and use our models to get a better understanding of the limitations that the simple Keplerian thin-disk approximation and the limited spectral resolution put on our interpretation of the line width.

4.1.3 Impact of the wind velocity on the line width

The presence of a wind can affect the shape of the line profile and therefore also affects the measured FWHM in the disk. This was already discussed in detail by Weber et al. (2020) for [OI] 0.63 µm and other atomic tracers. An important conclusion of Weber et al. (2020) is that in most cases the width of the emission lines is not a good indicator of the physical region from where the lines were emitted. At low inclinations, the line width is not determined by Keplerian rotation, and at higher inclinations, different velocity components are projected into the same observed velocity bin, which makes the localisation of the physical emitting region difficult.

For the wind tracers studied in this work, the most interesting case is the face-on view (i = 0°), for which Keplerian broadening does not play a role (e.g. TW Hya). As can be seen from Fig. 11, in the model PE-1 (lowest L_acc), the measured FWHM of o–H₂ 2.12 µm is larger than for [OI] 0.63 µm. This is caused by the traced wind velocities. In the PE-1 model, H₂ survives higher up in the wind region (see Fig. 4) than in the other models, and consequently the profile is dominated by the wind velocities, causing the line broadening (see top-panel of Fig. F.1). The average wind velocity in the o–H₂ 2.12 µm -emitting region is ≈ 2.5 times higher than in the region emitting [OI] 0.63 µm (see Table C.1). The opposite is true for the PE-3 model (highest L_acc), where H₂ is more efficiently dissociated and the line emission origin moves towards the disk surface where the wind velocities are approaching zero. In the models with reduced dust mass in the wind, the FWHMs of the two lines are almost identical. In those models, both line-emitting regions are closer to the disk surface and trace generally lower wind velocities, and the impact on the line broadening is limited. Nevertheless, the models predict that the ratio of the FWHM also depends on the stellar properties (at least at low inclinations) and show the potential of using atomic and molecular tracers for studying the wind velocity structure in different regions. However, a comprehensive study of this effect requires more observations of targets with varying stellar properties and low-inclination disks.

For the two lines studied here, the impact of the wind velocities on the FWHM should be similar, and hence it is unlikely that this is why we measure a similar FWHM for the two lines using the Keplerian thin-disk approximation.

Fig. 12 Same as Fig. 11 but a spectral resolution of R=100 000 is used for the modelled line profiles.

4.1.4 Impact of the spectral resolution on line width

The measured FWHM of the line profiles also depends on the spectral resolution of the observations. Fig. 12 shows the FWHM ratio of the two lines again but the modelled line profiles are now convolved to a spectral resolution of R = 100 000 (i.e. the same models as discussed in Sect. 4.1.1; see also Fig. 6 for an example). A comparison with the results shown in Fig. 11 (low resolution) shows that, for 14 models the ratio increases, for 6 it decreases, and for 10 it remains similar (within 5%). The biggest changes are for sin(i) ≈ 0.85 where we now see clear differences in the FWHM of [OI] 0.63 µm and o–H₂ 2.12 µm, but for the models with reduced dust content the ratios at that inclination are still similar. Nevertheless, for the modelled photoevaporative wind scenario, the lower spectral resolution of R ≈ 50 000 tends to make the FWHM of the two lines more similar, and higher spectral resolution observations would be helpful in order to better constrain the line origin. However, the general picture for the models remains the same, in particular the impact of the wind velocity on the FWHM, as discussed in Sect. 4.1.3.

4.1.5 The thin-disk approximation

In addition to the wind velocity field, the Keplerian thin-disk approximation used for Eq. (1) also neglects the vertical position of the emission region in the disk, with the latter actually affecting the rotation velocity (e.g. Rosenfeld et al. 2013). Looking at the physical emitting regions in the model (see Fig. 4), one can see that [OI] 0.63 µm and o–H₂ 2.12 µm are emitted from different heights. This means that even if the two lines emit from the same radius, the measured velocities will be different.

Figure 11 also shows the FWHM ratios for PE-2 models, where we replace the actual velocity field with pure Keplerian rotation, also neglecting the height of the disk (models with suffix Keplerian). This should simulate the thin-disk approximation, although we note that, in the model, the emission is still vertically extended. To derive the FWHM for those pure Keplerian line profiles, we followed the same Gaussian fitting procedure as for the wind models (see Fig. F.3).

Figure 11 shows that for pure Keplerian models, the FWHM ratio for the PE-3 model series remains similar to that in the wind models. The reason is that the PE-3 models show the weakest wind features in the line profiles. For both lines, the emission is coming from closer to the disk surface, and in case of o–H₂ 2.12 µm there is less H₂ in the wind due to the strong FUV radiation field (see Fig. 4). Therefore, the line profiles are mostly dominated by the Keplerian velocity field of the disk in any case.

For the PE-1 and PE-2 models, the FWHM ratio increases for i ≈ 20°–60° by factors of up to ≈ 1.7 and ≈ 1.4, respectively. The PE-1 models show the strongest wind features in the line profiles, especially for o–H₂ 2.12µm, as more H₂ can survive in the wind causing broader line profiles in general. In the Keplerian models, this broadening due to the wind is neglected and therefore the FWHM ratio increases. This implies that considering a more realistic velocity field (such as that in the models) brings the measured FWHMs for the lines closer to each other and therefore makes their emitting radii appear more similar, although they are not.

Nevertheless, the general picture for the measured FWHM in the models remains similar, and is also consistent with the observational data for the FWHM. This means that the trend in the FWHM ratio of the models is mainly driven by inclination. This is a consequence of the different physical emitting regions. The [OI] 0.63 µm always comes from radii ≲2 au tracing higher rotation velocities, whereas the o–H₂ 2.12 µm flux is dominated by the emission from a larger area and only relatively little emission is coming from the inner region. For increasing inclinations, the measured projected velocities increase and hence the line gets broader. Although this happens for both lines, for [OI] 0.63 µm this high-velocity emission from the inner disk has a stronger impact on the spatially integrated line profile as the total emitting area remains small.

4.2 Theo–H₂ 2.12 µm luminosity

As mentioned above, our models tend to systematically under-predict the o–H₂ 2.12 µm line luminosities (see Sect. 3.2). Although we do not expect this to impact our conclusions on the line kinematics, we feel it necessary to discuss the potential reasons for the overly low o–H₂ 2.12 µm luminosities and possibilities for improvements of the models. In Fig. 13, we present the luminosity ratio of o–H₂ 2.12 µm to [Ol] 0.63 µm () as a function of L_[OI], as this more clearly shows the overly low with respect to L_[OI] in the models. It is apparent from Fig. 13 that especially for L_[OI] around 10⁻⁵ L_⊙ our models under-predict by about an order of magnitude.

One reason might be that, for such cases, the layer emitting the majority of o–H₂ 2.12 µm is colder than the [OI] 0.63 µm – emitting regions, pointing towards an inaccurate vertical temperature gradient in our models. Further investigation of this potential issue would require a detailed study of the heating and cooling mechanism for models with varying stellar and disk properties that cover this observational parameter space. Furthermore, the detailed line excitation mechanisms for o–H₂ 2.12 µm mightneed to be re-evaluated (see e.g. Nomura & Millar 2005; Nomura et al. 2007), but such investigations are out of the scope of this paper.

Another possibility is that for targets with high , the o–H₂ 2.12 µm line luminosity also includes contributions from outflows and jets (e.g. DG Tau). The spatially resolved o–H₂ 2.12 µm observations of Beck & Bary (2019) indicate a complex spatial distribution of the o–H₂ 2.12 µm line emission for some targets, indicating shock-excited emission. Although such emission should usually show higher velocities, a potential contribution to the NLVC cannot be excluded. We note that, for GM Aur, the observations of Beck & Bary (2019) indicate compact and centralised emission, pointing towards a disk origin. For this target, the predicted o–H₂ 2.12 µm flux of the GM Aur DIANA model is in reasonable agreement with the observations (within a factor of three), indicating that for such targets our model works quite well.

As the o–H₂ 2.12 µm traces non-Keplerian velocity fields in several cases, dynamical effects might also be highly relevant. For the wind scenario, additional molecular hydrogen could be transported into the wind regions and might survive long enough to contribute significantly to the line flux. However, even without a disk wind, dynamical effects such as advection of gas through the ionisation front at the H⁺/H/H₂ transition, as studied by Maillard et al. (2021) for photo-dissociation regions, might also have significant effects for the disk scenario. However, to study such effects, more sophisticated models are required (see discussion in Sect. 4.3).

Another important point is that both the observational sample and the presented models are biased. From Fig. 5, it is clear that the observations of Gangi et al. (2020) are only sensitive down to of ≈10⁻⁶ L_⊙ and therefore might miss a large sample of disks with low (i.e. such as TW Hya). The bias in the model stems from the facts that we only present a fiducial disk-wind model and the DIANA models were not selected to cover the physical properties of the Gangi et al. (2020) sample. To better understand the origin of the o–H₂ 2.12 µm emission and the shortcomings of the models presented here, both increasingly sensitive observations and a detailed modelling of selected targets covering the observed range of physical properties are necessary.

Fig. 13 Line luminosity ratio of o–H2 2.12 µm and [Ol] 0.63 µm as a function of the [Ol] 0.63 µm luminosity. The observational data points and models shown are the same as in Fig. 5. We omit the one DIANA model with the lowest line luminosities for clarity. The upper limits for the line ratios (downward arrows) are a consequence of the o–H2 2.12 µm line luminosity upper limits.

4.3 Model limitations and future aspects

In this work, we focus on observables of disk winds from molecular and atomic species by post-processing radiation–hydrodynamical simulations with a thermochemical code. Although the radiation–hydrodynamic models used consider thermochemical processes for atoms (i.e. photo-ionisation), molecules are neglected and are only included in the postprocessing step. However, thermochemical processes can also influence the physical disk-wind properties; for example molecular cooling might affect the thermal structure and hence the wind-launching region. There are a few models that include thermochemical processes in the (magneto) hydrodynamic modelling, compensating for the high computational cost with drastic approximations in the radiative transfer (e.g. Wang & Goodman 2017; Nakatani et al. 2018a,b; Wang et al. 2019; Gressel et al. 2020). The majority of these models are not suitable for producing synthetic observations, given the uncertainties on the derived temperature structures, and this was only attempted in Gressel et al. (2020), although without a direct comparison to existing observations.

While molecular cooling is not expected to affect the wind structure and mass-loss rates significantly (Owen et al. 2010; Sellek et al. 2022), a computationally efficient self-consistent approach with an appropriate treatment of the radiative transfer would be highly desirable. To that end, we are in the process of coupling the existing hydrodynamic photoevaporative wind models with the thermochemical code PRIZMO (Grassi et al. 2020), especially designed for this purpose. We will make use of the existing models such as MOCASSIN and PRODIMO to verify the new code but also to identify the most relevant thermochemical processes that need to be considered in a self-consistent model, the aim being to build a computationally efficient model that allows for direct comparison to observations.

Furthermore, modelling of various potential wind tracers is required. To the best of our knowledge, this study is the first attempt at comparing disk-wind model predictions for both atomic and molecular wind tracers directly to observations. However, to fully understand disk winds and their origin, multi-wavelength studies and observations are necessary. For example, the models of Gressel et al. (2020) indicate that the CI line at 492 GHz (609 µm) is a tracer that could potentially be used to tell apart a magnetocentrifugal from a photoevaporative wind. More recently, Xu et al. (2021) presented observations of the CII resonance line at 1335 Å of a sample of protoplanetary disks. These authors interpret the observed absorption line profiles with a semi-analytic thermal-magnetic wind model, indicating that both wind mechanisms might be important. Furthermore, CO ro-vibrational lines in the mid-infrared also show signatures of winds and outflows (Banzatti et al. 2022 and references therein). However, modelling these various disk wind tracers with self-consistent models could be challenging or even unfeasible due to computational limitations.

A similar approach as used in this work might be preferable, and we plan to use existing dynamic wind-model grids (e.g. Picogna et al. 2021) to explore the impact of disk structure (e.g. transitional disks) and stellar properties on atomic and molecular observables in various wavelength regimes. This will also include predictions for spectrally and spatially resolved observations (including spectro-astrometry, e.g. Whelan et al. 2021) of disk wind tracers. Currently, this is only feasible for molecular hydrogen to some extent (see Beck & Bary 2019), but upcoming facilities such as the Extremely Large Telescope (ELT) and potentially GRAVITY+³ at the Very Large Telescope will open up new possibilities and will increase the sample size of spatially resolved observations. Studying models that produce suitable synthetic observables for those instruments will be especially useful for identifying the wind tracers that allow us to distinguish wind-driving mechanisms from observations.

5 Summary and conclusions

In this work, we present photoevaporative disk-wind models to interpret atomic ([OI] ¹D₂ – ³P₂ at 0.63 µm) and molecular (o–H₂ 1 – 0 S(1) at 2.12 µm) line emission from protoplanetary disks that show signatures of winds. The modelling approach consists in the post-processing of hydrodynamic X-ray photoevaporative disk-wind models, using the radiation-thermochemical model PRODIMO to produce synthetic spectral line emission.

We used this framework to model the observations of [OI] 0.63 µm and o–H₂ 2.12 µm from Gangi et al. (2020), including the line fluxes and the kinematic signatures. In this first attempt, we neglect the dynamics for the chemical calculations. Therefore, the models might under-predict the wind signatures (i.e. blueshift of the peak emission) of molecular hydrogen for cases with strong and dense winds. Nevertheless, we find that X-ray-driven photoevaporative wind models are generally consistent with the observed wind signatures for both the [OI] 0.63 µm and o–H₂ 2.12 µm spectral lines in the currently available observational data.

We find blueshifted peaks in the o – H₂ 2.12 µm synthetic line profiles in the range of v_p = 0 to ≈ −6 km s⁻¹. Also, the measured FWHMs of the modelled lines are consistent with the data. However, we find that it is not required that both lines be emitted from the same regions close to the star, as suggested by Gangi et al. (2020) to explain the observationally derived kinematic wind signatures. In our models, the outer radius of the physical emitting regions for the o–H₂ 2.12 µm is always significantly larger than for [OI] 0.63 µm, that is, by factors of a few to ≈ 10. The complex velocity field structure, the different vertical emitting layers, and the limited spectral resolution make the measured FWHM appear similar. Thus, we conclude that for currently available observational data, the approach of using the FWHM of the [OI] 0.63 µm and o–H₂ 2.12 µm spectral line profiles and assuming only Keplerian broadening is not well suited to deriving the physical emitting radii of those two lines. Therefore, this method does not allow us to discriminate between magnetically driven and photo-evaporative disk winds. Further modelling of atomic and molecular line emission for both photoevaporative and magnetically driven winds is required in order to determine whether or not it is at all possible to determine the wind-driving mechanism from such spatially unresolved line emission.

The observational data are still scarce and most targets in the sample of Gangi et al. (2020) do not show clear signatures of blueshifted emission in molecular hydrogen, which might be a consequence of the limited spectral resolution of the data. On the other hand, our models indicate that simply not enough molecular hydrogen can survive in low-density winds (such as photo-evaporative winds) in particular if the target has a high accretion luminosity, which results in efficient photo-dissociation of molecular hydrogen by FUV radiation. Consequently, any o–H₂ 2.12 µm emission becomes dominated by that from the high-density regions at the disk surface where the kinematics becomes dominated by Keplerian rotation.

Further improvements to the models are required. In particular, an efficient coupling of the dynamics, chemistry, and thermal balance will allow more accurate predictions of molecular line emission to be produced. However, those models are computationally very expensive, meaning that a modelling approach such as the one presented here will still be required for a comprehensive exploration of the vast parameter space (e.g. disk structure, varying stellar properties, chemical networks). In particular, such models will also be able to produce spatially resolved synthetic observables, which will be useful for exploring the capability of upcoming facilities such as the ELT.

Appendix A Dust density structures

In Fig. A.1 we show the dust density structure for our fiducial model and for the model where we lowered the dust density in the wind by a factor of 10 (gd models), resulting in a gas-to-dust mass ratio of 1000 (see Sect. 2.4).

Fig. A.1 Dust density structure for the models with a constant gas-to-dust mass ratio of g/d = 100 in the disk and wind (top panel) and the model with a factor 10 lower dust density g/d = 1000 in the wind region (bottom panel). The red solid contour shows where the radial visual extinction AV,rad equals unity. The white dashed contours correspond to the values shown in the colour bar.

Appendix B Comparison to MOCASSIN

Here we compare the results of the PRODIMO model to the results of the MOCASSIN model of Weber et al. (2020). We do not attempt to do a full benchmark test. Rather, we are interested in how the two codes compare if the same physical input properties are used, such as the disk–wind density structure, stellar properties, and dust opacities. The two codes are conceptually different, in particular for the chemistry (i.e. PRODIMO includes molecules) and the heating and cooling processes. It is therefore especially interesting to see the differences in the code for e.g. the temperature structure.

Fig. B.1 shows a comparison for the electron number density ne (ionisation degree) and the gas temperature T_g for the 2D structure. Fig. B.2 shows the comparison for vertical cuts through the disk structure. These figures show that the general features in n_e and T_g are very similar. However, there are also differences. The regions with T_g > 10000 K are more extended in the MoCASSIN model. This is likely caused by differences in the X-ray radiative transfer. MoCASSIN places the EUV/X-ray source above and below the star (≈ 10 R_*) whereas in PRODIMO the EUV/X-ray source is at the position of the star (for technical reasons). Placing the emitting source higher above the disk midplane allows the radiation to penetrate deeper into the disk/wind. This will also affect ne in the wind due to direct ionisation of hydrogen. However, those effects should not be significant and they only affect low-density regions (Ercolano & Glassgold 2013). Otherwise, the gas temperatures agree within a factor of two in the disk wind region (i.e. high up in the disk; see also Fig. B.2). At the transition from the wind region to the disk the differences in T_g become significant, which is expected as MoCASSIN does not include molecules, which can become important coolants in that region.

For the electron density (abundance), the differences are more significant (about a factor of 10), in particular in the outer disk (r ≳ 3 au). This is somewhat expected, because PRODIMO includes molecules, which affect the ionisation chemistry. Furthermore, differences in wavelength-dependent ionisation cross-sections likely also play a role in addition to the differences in the radiative transfer methods.

Considering that PRODIMO and MOCASSIN are conceptually different codes, the results compare reasonably well. We therefore conclude that using PRODIMO in combination with radiation-hydrodynamic photoevaporative wind models is a suitable approach that also allows molecular line emission to be modelled, in particular molecular hydrogen.

Fig. B.1 Comparison of the MOCASSIN and PRODIMO model. The first row shows the MOCASSIN model, the bottom row the PRODIMO model. From left to right the 2D structure of the hydrogen density (as reference, identical in both models), electron density and the gas temperature for the inner 10 au are shown.

Fig. B.2 Comparison of the PRODIMO (orange colour) and the MOCASSIN (blue colour) model for vertical cuts through the disk at r ≈ 0.3, 1, 3 and 10 au (columns from left to right). The rows (top to bottom) show the comparison for the hydrogen number density (for reference, identical in the two models), gas temperature, and electron abundance (relative to the total hydrogen density), respectively.

Appendix C Physical properties in the line-emitting regions

Table C.1 summarises certain physical properties of the line-emitting regions shown in Fig. 4. To calculate those quantities, we simply average the respective quantity over the whole emitting region. From an observational perspective, those values are only strictly valid for a disk viewed face-on, as not all emission is seen depending on the disk inclination. Nevertheless, those values are representative of the physical conditions for line excitation. We note that a mass-weighted average does affect the absolute numbers, but the general picture and trends remain the same.

Appendix D Line luminosities from the literature

In addition to the data from Gangi et al. (2020), we also collected data from the literature for targets with observations for both the [OI] 0.63 µm and o–H₂ 2.12 µm lines. The collected line luminosities and the adapted distances to calculate them are listed in Table D.1. We note that for this work we only use the luminosities from the literature data as the kinematic data for the o–H₂ 2.12 µm were often not available (i.e. due to low spectral resolution).

From the collected targets, only GM Aur is also included in Gangi et al. (2020), and these authors report 0.18 ± 0.02 × 10⁻⁵ L_⊙ and 0.4 ± 0.1 × 10⁻⁵ L_⊙, for the o–H₂ 2.12 µm and the [OI] 0.6 µm line, respectively. Those values are a factor 4.6 higher for o–H₂ 2.12µm and a factor 4.2 lower for [OI] 0.6 µm, compared to the values listed in Table D.1. The origin of those differences are unclear, but might be related to variability (although this would be rather strong, see Gangi et al. 2020), different definitions of the NLVC for [OI] 0.63 µm or unknown systematic errors (i.e. observations were done with different instruments).

Appendix E Selected targets for comparison to the models

Table E.1 summarises the kinematic properties of the targets used for the comparison of the measured line kinematics to the models. All those targets have a blueshifted NLVC in both [OI] 0.63 µm and o–H₂ 2.12 µm. All data listed in Table E.1 are from Gangi et al. (2020).

Appendix F Gaussian fitting of line profiles

Figs. F.1 and F.2 show the results of the Gaussian fitting process (see Sect. 2.3) for all photoevaporative models presented in this work using the observational spectral resolution of R ≈ 50000.

Fig. F.3 shows the fitting results of the modelled profiles assuming a purely Keplerian velocity field. As opposed to the models with the wind velocity field, these latter profiles are always symmetric and hence the fitting routine always gives v_p = 0 km s⁻¹. The impact of Keplerian assumption concerning the FWHM is discussed in Sect. 4.1.5.

Fig. F.1 Fitting of the modelled line profiles. From top to bottom, the models PE-1, PE-2, and PE-3 are shown. Please note the different scaling on the y axis (flux) for each panel. For each model, we show the [OI] 0.63 µm (first row) and the o–H2 2.12 µm (second row) for five different inclinations (columns from left to right). Each panel shows the line profile with the model resolution (blue) convolved to a resolution of R = 50000 similar to the observations (black), and the fit (the NLVC) to the convolved profile (red solid line).

Fig. F.2 Same as Fig. F.1 but for the models with a factor of ten lower dust density in the wind. The red and orange dotted lines show the individual components of the fit (not seen in the case of a single Gaussian). In the case of multiple Gaussians, the orange dotted line indicates the one chosen as the NLVC used to determine vp and the FWHM.

Fig. F.3 Same as Fig. F.1 for the PE-1, PE-2, and PE-3 models, but assuming a purely Keplerian velocity field (see Sect. 4.1.5 for details).

Acknowledgements

The authors are grateful to the referee for a constructive and positive report that improved the paper. The authors want to thank M. Gangi and B. Nisini for useful discussions on their data and for providing some of it in digitized form. We acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Research Unit “Transition discs” – 325594231. This research was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. CHR is grateful for support from the Max Planck Society. This research has made use of NASA’s Astrophysics Data System. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), matplotlib (Hunter 2007) and scipy (Virtanen et al. 2020).

References

All Tables

All Figures

	Fig. 1 Two-dimensional disk and wind gas density structure and velocity field (black arrows). This structure is the same in all the presented models. The red solid contour shows where the radial visual extinction Av,rad is equal to unity; here, we refer to this layer as the disk surface. The white dashed contours correspond to the values shown in the colour bar.
In the text

	Fig. 2 Stellar input spectra for the three different accretion luminosities Lacc.
In the text

Fig. 3 Two-dimensional gas temperature Tgas, dust temperature Tdust, and FUV radiation field χ (in units of the Draine field Draine & Bertoldi 1996; Woitke et al. 2009) for the PE-2 (top row) and PE-2 gd (bottom row) models. The red solid contour shows where the radial visual extinction AV,rad is equal to unity, which roughly separates the wind region and the rotating disk. The white dashed contours correspond to the values shown in the colour bars.

In the text

Fig. 4 Main emitting region for the [OI] 0.63 µm (orange boxes) and o–H2 2.12 µm (red boxes) spectral lines. The columns show the three models with varying Lacc (increasing from left to right). The top row shows the models with g/d = 100 and the bottom row with g/d = 1000 in the wind. The coloured contours show the molecular hydrogen number density . The borders of the line emission boxes indicate where the total emission reaches 15% and 75% in the radial (integrated inside-out) and vertical (integrated from top to bottom) directions. The solid white line shows where the radial visual extinction AV,rad reaches unity, and the white dotted line Tgas = 1000 K. The grey arrows indicate the 2D velocity field (vertical and radial components of the wind).

In the text

Fig. 5 o–H2 2.12 µm versus [OI] 0.63 µm line luminosities. The dark grey points with error bars show the observations from Gangi et al. (2020, 2021), and the light-grey points show observational data collected from the literature (see Appendix D for details). The orange symbols show the photoevaporative disk-wind models (inclination i = 40°) with varying accretion luminosities (lowest(highest) line luminosities are for the lowest(highest) accretion luminosity). The square symbols are for the models with reduced dust-to-gas mass ratio in the wind. The brown diamonds are the results from the DIANA project (PRODIMO models without a wind component). We mark DG Tau as an example of a target that likely cannot be modelled with a pure photoevaporative wind. We also mark the two independent observations for GM Aur (see Appendix D).

In the text

Fig. 6 Example of the line-fitting results for the PE-1 model. The top figure shows the results for a spectral resolution R = 50 000 (i.e. similar to the observations). The bottom figure shows the same model, but the line profiles were convolved to a spectral resolution of R = 100 000 (see Sects. 4.1.1 and 4.1.4). Each individual panel shows the line profile with the model resolution (blue), convolved to the target spectral resolution (black) and the Gaussian fit (the NLVC) to the convolved profile (red solid line). In case of multiple Gaussians (dashed and dotted lines), the orange dashed line indicates the one chosen as the NLVC. This component is used to determine vp and the FWHM (also given in each panel). The fitting results for all other models are shown in Appendix F.

In the text

Fig. 7 Velocity peak location vp of o–H2 2.12 µm vs vp [OI] 0.63 µm. Filled circles with error bars show the observed values; the black colour is for targets with vp < 0 km s−1 (i.e. blue-shifted wind component) for both the o–H2 2.12 µm and [OI] 0.63 µm; light-grey marks targets with vp > 0 km s−1 (i.e. no clear indication of a disk wind). The coloured symbols are for the models. Shown are models for the three different accretion luminosities and with a reduced dust-to-gas mass ratio in the wind (coloured diamonds). The size of the coloured symbols scales with the inclination; smallest symbols are for i = 0° largest are for i = 80°. The thick dashed line indicates vpH2 = vp[OI] and the light grey stripe indicates vpH2 = vp[OI] ± 3 km s−1. Top panel: all data points; bottom panel: zoom-in (brown box) showing all the models and the observational data for targets with clear wind signatures, excluding DO Tau.

In the text

Fig. 8 Velocity peak location vp versus disk inclination for the [OI] 0.63 µm (top panel) and o–H2 2.12 µm (bottom panel). Filled circles with error bars show the observed values; the black colour is for targets with vp < 0 km s−1 (i.e. blueshifted wind component) for both the o–H2 2.12 µm and [OI] 0.63 µm; light-grey marks targets with vp > 0 km s−1 (i.e. no clear indication of a disk wind). The coloured symbols are for the models. Shown are the same models as in Fig. 7.

In the text

Fig. 9 HWHM (left panel: [OI], right panel H2) scaled by the square root of the stellar mass as a function of inclination. The coloured symbols show the models; the gd models (diamonds) are slightly shifted along the inclination axis for clarity. The black (grey) symbols show the observations, where the black coloured symbols mark the targets that show blueshifted peaks in the respective line profile. The solid lines show the maximum emission radius RK derived for the given HWHM and inclination assuming Keplerian rotation for a thin disk (see Eq. (1)).

In the text

	Fig. 10 Same as the bottom panel of Fig. 7, but for the modelled line profiles a spectral resolution of R = 100 000 is used.
In the text

	Fig. 11 Ratio of the FWHM of [OI] 0.63 µm and o–H2 2.12 µm. Shown are the same models as used for Fig. 9. Additionally, we show the model series called Keplerian (crosses), where we assume that the velocity field in the disk and wind is purely Keplerian (i.e. thin-disk approximation). The black data points with error bars mark the targets with blueshifted vp in both lines.
In the text

	Fig. 12 Same as Fig. 11 but a spectral resolution of R=100 000 is used for the modelled line profiles.
In the text

	Fig. 13 Line luminosity ratio of o–H2 2.12 µm and [Ol] 0.63 µm as a function of the [Ol] 0.63 µm luminosity. The observational data points and models shown are the same as in Fig. 5. We omit the one DIANA model with the lowest line luminosities for clarity. The upper limits for the line ratios (downward arrows) are a consequence of the o–H2 2.12 µm line luminosity upper limits.
In the text

	Fig. A.1 Dust density structure for the models with a constant gas-to-dust mass ratio of g/d = 100 in the disk and wind (top panel) and the model with a factor 10 lower dust density g/d = 1000 in the wind region (bottom panel). The red solid contour shows where the radial visual extinction AV,rad equals unity. The white dashed contours correspond to the values shown in the colour bar.
In the text

	Fig. B.1 Comparison of the MOCASSIN and PRODIMO model. The first row shows the MOCASSIN model, the bottom row the PRODIMO model. From left to right the 2D structure of the hydrogen density (as reference, identical in both models), electron density and the gas temperature for the inner 10 au are shown.
In the text

	Fig. B.2 Comparison of the PRODIMO (orange colour) and the MOCASSIN (blue colour) model for vertical cuts through the disk at r ≈ 0.3, 1, 3 and 10 au (columns from left to right). The rows (top to bottom) show the comparison for the hydrogen number density (for reference, identical in the two models), gas temperature, and electron abundance (relative to the total hydrogen density), respectively.
In the text

Fig. F.1 Fitting of the modelled line profiles. From top to bottom, the models PE-1, PE-2, and PE-3 are shown. Please note the different scaling on the y axis (flux) for each panel. For each model, we show the [OI] 0.63 µm (first row) and the o–H2 2.12 µm (second row) for five different inclinations (columns from left to right). Each panel shows the line profile with the model resolution (blue) convolved to a resolution of R = 50000 similar to the observations (black), and the fit (the NLVC) to the convolved profile (red solid line).

In the text

	Fig. F.2 Same as Fig. F.1 but for the models with a factor of ten lower dust density in the wind. The red and orange dotted lines show the individual components of the fit (not seen in the case of a single Gaussian). In the case of multiple Gaussians, the orange dotted line indicates the one chosen as the NLVC used to determine vp and the FWHM.
In the text

	Fig. F.3 Same as Fig. F.1 for the PE-1, PE-2, and PE-3 models, but assuming a purely Keplerian velocity field (see Sect. 4.1.5 for details).
In the text