Line Profiles of Forbidden Emission Lines and What They Can Tell Us About Protoplanetary Disk Winds

Emission in forbidden lines of oxygen, neon, and other species are commonly used to trace winds from protoplanetary disks. Using Cloudy, we calculate such emission for parametrized wind models of the magnetothermal type, following Bai et al. These models share characteristics with both photoevaporative and magnetocentrifugal winds, which can be regarded as end members, and are favored by recent theoretical research. Both broad and narrow low-velocity components of the lines can be produced with plausible wind parameters, something that traditional wind models have difficulty with. Line luminosities, blueshifts, and widths, as well as trends of these with accretion luminosity and disk inclination, are in general accordance with observations.


INTRODUCTION
Protoplanetary disks (PPDs) of molecular gas and dust surround most low-mass pre-main-sequence stars (protostars).As the term suggests, PPDs are understood to be the sites and sources of planet formation, so their properties and evolution are of great interest.Observations show that PPDs typically disperse within a few Myr, presumably through some combination of accretion onto the star, sequestration in planets, and arXiv:2311.01802v1[astro-ph.SR] 3 Nov 2023 outflows (winds) (Ercolano & Pascucci 2017;Manara et al. 2023;Alexander et al. 2014).
Spectroscopic observations show that PPDs are often accompanied by outflows at hundreds of km s −1 as well slower outflows at ≲ 10 km s −1 (Hartigan et al. 1995).In some of the more luminous Class I&II systems, the fast outflows develop into resolved, narrowly collimated jets that may extend for parsecs and terminate in Herbig-Haro objects (Bally 2016;Pascucci et al. 2023, and references therein).Estimated mass-loss rates for both types of outflow correlate with accretion rates (Hartigan et al. 1995;Ellerbroek et al. 2013;Rigliaco et al. 2013;Natta et al. 2014).It is generally believed that magnetohydrodynamic (MHD) mechanisms are needed to accelerate the jets, but it is not known whether these are launched very close to the star or more broadly from the PPD (Frank et al. 2014).A disk origin seems more secure for the lower-velocity component because of the dependence of line widths on inclination, the lower blue shifts themselves, and the higher estimated mass-outflow rates (Fang et al. 2018;Simon et al. 2016;Banzatti et al. 2019).
Whether fast or slow, magnetized disk winds naturally relate outflow to accretion by extracting orbital angular momentum from the disk.Accretion may occur via a turbulent disk viscosity independently of any wind, however, or of a magnetic connection between the disk and its wind (Ercolano & Pascucci 2017).In most accretion disks, the main source of turbulence is probably magnetorotational instability (MRI) (Lesur et al. 2022); if magnetic fields are needed to explain the disk's effective viscosity, then it would seem natural that any outflow from the disk should be magnetized as well.In PPDs, however, theoretical work suggests that MRI may be suppressed in regions likely associated with the low-velocity optical line emission (∼1-20 AU) (Gammie 1996;Bai & Stone 2013;Turner et al. 2014).Hydrodynamic mechanisms not requiring magnetic fields exist in principle, but it is not clear whether they can produce the level and type of turbulence needed to drive accretion in PPDs at observed rates (Turner et al. 2014;Fromang & Lesur 2019).
In any case, low-velocity disk winds almost certainly exist, so if they are not magnetized, then they must be launched by thermal pressure.We use the term "photoevaporative" for unmagnetized disk winds because in all models of such winds of which we are aware, the main source of heat is high-energy photons (FUV, EUV, X-ray) impinging on the surface of the disk and raising the sound speed to a significant fraction of the local escape velocity (Hollenbach et al. 1994;Gorti & Hollenbach 2009;Gorti et al. 2016, and references therein).There is little doubt that photoevaporation must occur, because protostars are observed to produce these hard photons-except perhaps in the EUV, i.e. 100eV > hν > 13.6eV, which is difficult to observe because it is heavily absorbed by the ISM.Photoevaporation can also be caused by irradiation of the outer parts of the PPD by other stars or protostars in the natal stellar association (Adams et al. 2004;Haworth et al. 2016, e.g.).Photoevaporative wind models have been successful in explaining many aspects of the observations, including line profiles and ratios (Ercolano & Owen 2016).Pure photoevaporation cannot drive accretion directly because, absent a magnetic connection, the outflow exerts no torque on the disk.
Nevertheless, a correlation between accretion and outflow could be explained if the wind were driven primarily by FUV (for which there is evidence, e.g.Simon et al. 2016), since the FUV probably comes mainly from the accretion boundary layer.Photoevaporative winds probably cannot explain the broader component (FWHM> 40 km s −1 ) of the low-velocity emission (blueshift< 30 km s −1 ).Simon et al. (2016) found, from the relation between the lines' FWHM and disk inclination, that the broad component is broadened by the rotation of the disk at disk radii ≤ 0.5 AU.At those implied radii, gas at temperatures compatible with the observed lines would not be able to escape the gravitational potential, or at least not at rates sufficient to explain the line luminosity (Simon et al. 2016;Ercolano & Pascucci 2017;Weber et al. 2020).Whelan et al. (2021) and Fang et al. (2023) performed spatially-resolved measurements of the velocity profile of optical forbidden emission lines and found the emission to be inconsistent with a thermal wind and supportive of an MHD disk wind origin of these lines.
In addition to providing a dynamical linkage between accretion and outflow, magnetized winds have the additional advantage in principle of being able to accelerate the outflow to a given velocity without heating it as much as would be required by a purely photoevaporative wind, hence allowing the highervelocity gas to remain largely neutral, as suggested by some blue-shifted [OI] lines, and even molecular lines.However, many early magnetic models for PPD winds overdid this, explicitly or implicitly following the paradigm created by Blandford & Payne (1982, hereafter BP82): the BP82 solutions neglect thermal pressure altogether; more importantly perhaps, their fiducial solution has a relatively large Alfvén radius, so that the specific angular momentum of the outflow along a given field line is asymptotically 30 times larger than the orbital angular momentum where that line meets the disk (the footpoint or wind base).Since the angular momentum carried off by the wind balances what is lost by the disk, models based on this "magnetocentrifugal" paradigm tend to underpredict the ratio of the wind mass-loss rate to the accretion rate (Konigl & Pudritz 2000;Wardle & Königl 1993;Weber et al. 2020).An additional issue is that if the magnetic field is in near equipartion with the gas and well coupled to it even at the midplane, as in some of these models, then the surface density of the disk corresponding to observed accretion rates would be far too small to form planets; this however might allow such solutions to explain some transition disks, where substantial accretion is observed despite a drastic reduction of the surface density within a few to a few tens of AU (Wang et al. 2019).
Recent theoretical models and numerical simulations of MHD disk winds invoke weaker magnetic fields that are launched from upper layers of the disk where exposure to hard photons (FUV, EUV, & X-rays) provides sufficient ionization, and where the gas pressure and density are much lower than at the midplane (Bai & Stone 2013;Gressel et al. 2015;Bai et al. 2016;Bai 2017;Béthune et al. 2017;Wang et al. 2019;Gressel et al. 2020).Mass-loss rates tend to be high-comparable to accretion rates-and field lines heavily loaded, so that Alfvén radii (R A ) are only modestly larger than footpoint/launching radii (R 0 ).Magnetic and thermal pressure are at least as important as centrifugal force in launching such winds (except in some of those in Gressel et al. 2020), whence they have been dubbed "magnetothermal" by Bai et al. (2016).
Magnetothermal winds can be regarded as part of a continuum of which photoevaporative and cold magnetocentrifugal winds are end members.Unlike the former, they directly connect outflow to accretion; unlike the latter, they can reconcile observed accretion rates with disk surface densities that are high enough to form planets.To date, there has been very little systematic work to explore the emission-line luminosities and profiles predicted by magnetothermal winds. 1   The availability of disk-wind tracers that probe different physical conditions and regions, along with the modelling of these tracers, makes it possible to understand the different types of disk winds and their role in PPD evolution.Forbidden emission lines from oxygen (together with other atomic species) and low ionization stages of elements such as sulfur and neon are the prime candidates for tracing disk winds, and their line ratios, if modelled correctly, can provide information about the physical properties of the emitting region.The presence of a warm ionized disk wind was first confirmed through the detection of the [NeII]12.8µmfine structure line (Pascucci & Sterzik 2009;Alexander 2008;Gorti & Hollenbach 2009;Font et al. 2004).The line was often found to be blueshifted by a few km s −1 relative to the stellar velocity, indicating an origin in an outflow that is approaching the observer, while the emission from the receding outflow in the other half of the disk is obstructed by the dust.The lack of high-resolution IR spectra stands in the way of collecting a large sample of PPDs to identify evolutionary trends, or the details (or components) of the line profile.The correlation observed between the neon line (luminosity, width and blueshift) and the system variables (inclination, accretion luminosity, etc..) that were observed was successfully produced by thermal winds which led the community to accept that the [NeII] line, in fact, traces a photoevaporative wind.However, there is a degeneracy amongst photoevaporative wind models in predicting the [NeII] line depending on the type of radiation used to drive these winds; both X-ray and EUV driven winds can successfully predict the line, although with vastly different mass loss rates (∼ 2 orders of magnitude), and the latter fails to explain the emission of other forbidden lines (Ercolano & Owen 2016).
More progress could be made with optical forbidden lines, such as those of [OI] 6300 Å and [SII] 4068 Å, since high-resolution spectra are abundant in the optical.Optical forbidden lines, e.g.[OI] 6300 Å, have long been known to possess a low velocity component (LVC), emission blueshifted by less than 30 km s −1 , in addition to a high-velocity component (HVC) tracing fast (100 km s −1 ) collimated jets (e.g., Hartigan 1 at least for the optical forbidden lines, which are the main focus of this paper; Gressel et al. (2020) however have made interesting predictions for the [OI] 63.2µm line and for millimeter-wave lines of [CI] and HCN based on their MHD simulations.et al. 1995).EO16 simulated forbidden emission line profiles from a thermal wind that is driven by X-ray photoevaporation and found line luminosities that are in good agreement with observed LVCs.However, observational studies of high-resolution spectra by Rigliaco et al. (2013) and Simon et al. (2016) showed that the LVC can also be decomposed into a broad component (BC) with FWHM of at least 40 km s −1 and a narrower component (NC).We note that similarities between observed spectra of [OI] 6300 Å (composite BC and NC) and the modelled data from Ballabio et al. (2020) (that is not fitted with 2 gaussian components) led them to conclude that gaussian fitting is not the best approach and more complex, physically-motivated model fitting (such as wind solutions) may yield greater insight than a gaussian decomposition alone.
From the relationship between the full width at half maximum (FWHM) of the BC and NC and the inclinations of the sample of systems observed, Simon et al. (2016) found that the broad component likely forms in a wind that is launched at radii of 0.05-0.5AU(0.5-5AU for the NC) if the width reflects keplerian rotation.Launching a photoevaporative wind from within 0.5AU requires high temperatures to achieve escape velocity.Such high temperatures will tend to ionize the gas, leaving insufficient OI to explain the observed line emission.More evidence for physically distinct emission regions for the broad and narrow components is that the [OI] 5577/6300 line ratio is typically higher in the BC than the NC (Simon et al. 2016).Pascucci et al. (2020) detected the NC of the OI lines in spectra of both full disks (presumed to be younger PPDs) and disks depleted of dust in their inner regions (presumed to be older), albeit the emission was weaker for the latter group, which we shall refer to in this paper as "transition disks" (TDs). 2 The [NeII] line showed an opposite relationship where the [NeII] line is enhanced in transition disks and almost absent in full disks.This led Pascucci et al. (2020) to conclude that the Ne line traces an outer photoevaporative wind driven by X-ray emission that can penetrate the inner regions in depleted disks while the OI BC emission comes from an inner MHD wind.The dependence of the forbidden line emission on the accretion luminosity (FUV) and the X-ray luminosity, and the inference that [OI] line emission comes from distinct regions created a controversy around line excitation mechanisms.Collisions with electrons and neutral hydrogen, FUV pumping, and OH dissociation have all been proposed to explain the emission (Nemer et al. 2020).Weber et al. (2020) have recently computed optical forbidden-line emission for a magnetized wind based on the self-similar model of BP82.They were able to model the [OI] lines, confirming the previous finding that the broad low-velocity component traces an inner disk MHD wind.As explained above, however, the BP82 wind is colder and more strongly magnetized, with a larger Alfvén radius, than is currently thought appropriate for the low-velocity component of PPD winds; this causes Weber et al. (2020)'s LVC to be broader, more double-peaked, and more blue-shifted than is usually observed, as they themselves point out.Weber et al. (2020) and others (Pascucci et al. 2020;Simon et al. 2016) postulate that the observed emission of these lines likely traces multiple regions of the disk with different excitation and launching mechanisms, and that to reproduce these observed lines, one has to use a hybrid model that includes a combination of the above-mentioned processes.
In this paper we focus on the line profiles of oxygen, sulfur and neon to confirm that they can indeed trace a PPD wind, and whether those lines are able to distinguish between the different types of winds.We construct idealized wind solutions of the magnetothermal type following Bai et al. (2016), while varying parameters to include photoevaporative and magnetocentrifugal winds as limiting cases.We then postprocess these solutions with a radiation transport code (Ferland et al. 2017) to predict line luminosities and profiles.

The wind model
Following (Bai et al. 2016, hereafter BYGY), our wind solutions use a prescribed poloidal magnetic geometry: the field lines are straight, at angle θ ′ to the midplane, with footpoints at cylindrical coordinates (R 0 , z 0 ) on the "surface" of the disk.By default, z 0 = 0.15R 0 , to approximate our expectation that the base of the wind lies several scale heights above the midplane where FUV from the star is able to penetrate and provide sufficient ionization to couple the gas to the field.The strength of the poloidal field varies as R −1 near the wind base, and transitions smoothly to R −2 at R ∼ q −1 R 0 along the line, q being an adjustable parameter (see BYGY).We solve the usual steady, axisymmetric, ideal-MHD equations for the flow along each field line given a prescribed density ρ 0 and poloidal Alfvén speed V A0 = B P0 at the base.The wind also has a finite sound speed c s , which is taken to be constant along each field line.Thus the flow along each line is determined by the five3 dimensionless input variables (θ ′ , z 0 /R 0 , q, V A0 /V K0 , c s /V K0 ) plus three dimensionful scalings (R 0 , ρ 0 , V K0 ), the last of these being the keplerian orbital velocity at the wind base, V K0 = GM * /R 0 .After one has solved for the flow, the output variables include the positions of the critical points R slow magnetosonic < R Alfvénic < R fast magnetosonic and the mass and angular momentum fluxes along the line.Some more discussion of the fluxes is given in the Appendix.
Like Bai et al. (2016), we call these winds "magnetothermal" because they are accelerated by a combination of magnetic forces (represented by B P and V A ) and gas pressure (represented by c s ), whereas the magnetocentrifugal winds of Blandford & Payne (1982) are completely cold and are launched by magnetic forces alone.(In that case, it is necessary that z 0 = 0 and θ ′ < 60 • .)In the opposite limit that V A /c s = 0, the wind is launched entirely by the thermal pressure and sound speed of the gas, and therefore imitates a photoevaporative wind.We represent these unmagnetized winds by the same general scheme, with ghostly poloidal field lines that exert no force on the flow but coincide with the streamlines.There is then only one critical point on each line-the sonic point.While we postprocess these winds with Cloudy, as described below, the sound speed assumed in the dynamical model is not necessarily consistent with the temperature and molecular weight found by Cloudy.

Postprocessing with Cloudy
We have used the two-dimensional grids of emissivities from our radiation transport calculations and our grids of the wind structure to construct a three-dimensional axisymmetric volume and calculate the spectral profiles of our lines when viewed along different lines of sight (different inclinations).Based on the assumption that the disk midplane is optically thick, we omit the bottom half of the disk.The luminosities of the lines are calculated by integrating their emissivities over volume, as if each emitted line photon escapes the system.This approximation is important to our calculations, as our implementation of Cloudy is capable of radiative transfer in only one spatial dimension (the radial), so we comment here on its justification.The optical depths of the forbidden lines are small because of their small A coefficients, so line scattering is neglected.Extinction by dust grains is included but is relatively important.For the fiducial model, the visual extinction along a radial ray at the base of the wind (its densest part) is a few percent.
The hydrogen column density there is ∼ 10 22 cm −2 , the adopted grain size 0.1µm, and the gas-to-dust ratio ∼ 0.01% by mass (i.e., depleted by a factor of a 100 compared to the ISM).
Our idealized wind solution does not include the whole disk, but we expect little emission to come from the regions not covered.We assume that the midplane is opaque to all lines of interest, including the [NeII] 12.81µm line, and consider emission from one hemisphere only; this assumption is open to question, especially for the transition disks.Also, the wind region is confined by the poloidal field lines, so for any choice of field line angle other than 90 • , there will be an unmagnetized region for which the field lines would be anchored outside of the simulated region.Both of these regions are shown in blue in Fig. 1 representing a very low density region.
We apply Cloudy to the density and velocity fields of our wind models, together with appropriate luminosities in energetic photons, to calculate the [OI] 6300 Å and 5577 Å, [SII] 4068 Å, and [NeII] 12.81µm line emission.Our approach follows Nemer et al. (2020): we use the same radiation sources, chemical composition, grains, and stellar parameters.Namely, we represent the radiation source by four blackbodies with various effective temperatures.We used a blackbody with T eff = 4250K, M = 0.7M ⊙ , and R * = 2.5R ⊙ to represent the photospheric luminosity (L * = 1.83L ⊙ ).Another with T eff = 12000K was used to simulate the accretion radiation with a luminosity equal to the stellar luminosity (similar to the observed accretion luminosity (Rigliaco et al. 2013)).Finally we included two more blackbodies to mimic EO16's rather complex X-ray SED: one with T eff = 10 6 K, and another with T eff = 10 7 K, dividing EO16's X-ray luminosity of 2 × 10 30 ergs/s equally between the two.We used the same elemental abundances as reported by EO16.
We also include dust with uniform grain size 0.1µm (D'Alessio et al. 1998).We use dust-to-gas ratio 10 −4 by mass in the wind (Brown et al. 2013;D'Alessio et al. 2006;Gorti & Hollenbach 2008) However, we feed Cloudy different distributions of gas density and velocity appropriate to both photoevaporative and MHD disk winds.We concentrate on the inner-to-intermediate regions of the disk (0.04-4 AU), since most of the emission in the lines listed above comes from those regions.The calculations were performed on a spherical grid with variable cell size in the radial direction (adjusted automatically by Cloudy to optimize the calculation), and 361 points uniformly spaced in sin θ ∈ [0, 1].We integrate the line emissions over radius and over one hemisphere in (0 ≤ θ ≤ π/2), supposing that the disk obscures the emission from the other hemisphere.Cloudy is capable only of spherical (or slab) geometries, so we use Cloudy along radial rays through the (non-spherical) wind, each such ray having its own density profile.
Cloudy produces a discrete spectrum of all the simulated lines, which we endow with line profiles suited to the local gas temperature and flow velocity.Each cell on our poloidal grid represents a ring centered on the symmetry axis; we broaden the emission from the cell according to the gas temperature (we do not allow for turbulent broadening within cells) and assign its emission to velocity bins according to the projections of the flow onto the observer's line of sight at a specified inclination.The sum of the contributions from all such rings gives the line profile, and to account for instrumental broadening, we then convolve with a gaussian dispersion σ = 8 km s −1 to represent the resolution of current observations.The profiles simulated in this way are usually asymmetric and skewed to the blue side, as has been found by Ballabio et al. (2020).Observed spectra have a varying continuum level contributed by the stellar photosphere, as well as contamination from telluric lines.Uncertainties in the continuum level may influence the inferred line widths and (especially for asymmetric lines) blueshifts.
Information about the wind dynamics can be deduced from the line shape as can be seen in Fig 1-7.The blueshift is higher at the lowest inclinations (face-on) and decreases with increasing inclination for some of our models, but in other models, the blueshift increases with inclinations up to ∼ 20 o − 40 o , but then decreases at still higher inclination.A similar trend is reported from observations for the LVC (Banzatti et al. 2019) as expected for a conical wind.We fit the line profiles with a single gaussian, and then we obtain the blueshift and the FWHM from the first and second moment of the fit, respectively.We fit the gaussian to the top 99% of the line to avoid numerical errors while fitting down to the continuum level.Those fits are all listed in Table.2. We chose to fit with a single gaussian, contrary to the multi-gaussian fitting that is done for observations and then assigning each of the components to a distinct emission region, which shall be discussed further below.

Line data and comparison with observations
Table 1 shows the input parameters and line luminosities for representative models of each of the three types considered in this work: magnetothermal, photoevaporative, and magnetocentrifugal.The models are intended to be representative of the winds proposed by Bai et al. (2016), Ercolano & Owen (2016), and Weber et al. (2020), respectively.The inner radius of the disk is set to 0.04AU, and the streamlines, which coincide with magnetic field lines for the magnetized winds, thread the disk at an angle θ ′ = 60 • from the midplane.
1 Table 1.Salient properties of representative wind models.The sound speed (c s ) and the Alfvén speed (V A ) are normalized by the orbital speed except for the photoevaporative wind model, where c s is constant, as described in appendix A. Z 0 is the height of the wind base above the midplane in units of R 0 , the cylindrical radius of the wind base.L 1 is the [OI] 6300 Å luminosity, L 2 [OI] 5577 Å luminosity, L 3 [SII] 4068 Å, and The observed forbidden line luminosities range between ∼ 10 −7 − 10 −3 L ⊙ with a median around 10 −5 L ⊙ for the OI lines and 10 −6 L ⊙ for the [NeII] line, while the [OI] 6300 Å/[OI] 5577 Å ratio ranges from 1 to 15 with a median around 7. The range of blueshifts from observations is between 0 and 300 km s −1 for all lines; any line with blueshift more than 30 km s −1 is considered an HVC and below that it is an LVC.The range of FWHM is 10 to 200 km s −1 , with lines wider than 40 km s −1 considered to be BC (Simon et al. 2016;Fang et al. 2018;Banzatti et al. 2019).and density for our fiducial model (magnetothermal 1 in Table 1)

Magnetothermal wind
We choose the input parameters for our fiducial model ("magnetothermal 1 " in Table 1) to include the effects of both magnetic and thermal forces, and we optimize this input to achieve realistic mass loss and accretion rates.The details of our choices are summarized in the Appendix.Here we look at the spatial distribution of the line emission in our fiducial model.As can be seen from Fig. 1, the emission in the optical forbidden lines is concentrated in the inner regions (< 1AU) and low latitudes (10 − 30 • ), whereas the [NeII] emission is more extended and mostly from higher latitudes.There are multiple reasons for this disparity.We find that each line is most luminous in regions where the density is close to its critical density, density for our fiducial model with a different scaling to the base density, (magnetothermal 2 in Table 1).
as found previously by Ballabio et al. (2020).The line with the highest critical density ([OI] 5577 Å) comes from the innermost regions of our fiducial model, while [NeII] 12.81µm, which has the lowest critical density, comes from more extended regions.Gas temperature also has an effect.The relatively small excitation energy of the fine-structure [NeII] line (≈ 0.1eV) prefers lower temperatures (∼ 1000K) than the optical lines (∼ 8000K).As can be seen from Fig. 1, the inner regions are hotter than the outer regions.A third factor is the penetration of the ∼ 1keV X-rays needed to form [NeII], as discussed below.
The simulated line luminosities shown in Table 1 are generally within the range observed.We find that our predicted line luminosities are sensitive to the FUV continuum illuminating the wind, as is also observed ( §1).The FUV is the main heat source for our model winds, and collisional excitation of the forbidden lines depends on the gas temperature, as discussed above.The FUV heating comes mainly from the grain photoelectric effect, gas-grain collisions, H-photoionization, and photoionization of excited hydrogenic species.In addition, the oxygen lines can be excited by FUV pumping (Nemer et al. 2020).The line luminosities for the fiducial magnetothermal wind model are broadly consistent with the observations (Pascucci et al. 2020), and as expected, the line emission positively correlates with the FUV luminosity for the optical lines.This is expected, because the optical lines are excited by collisions, to which FUV contributes by heating the gas, and by FUV pumping.On the other hand, [NeII] line shows a negative correlation with the FUV flux in our models.This has not been reported by observers.The [NeII] line has a lower excitation energy, so the increase in temperature due to the increase in FUV flux tends decrease emission in this line .(Wan et al. 2019).As is well known (e.g.Draine 2011), the rate of excitation of ions by electron collisions tends to scale as T −1/2 , whereas the excitation of neutrals is exponentially suppressed if kT e is less than the energy of the transition, which is ∼ 2eV for the OI lines.
Line ratios can be used to constrain physical conditions in the emitting regions, provided that one understands how the upper states of the line transitions are populated.It was reported in EO16 that the observed [OI] 6300 Å/([OI] 5577 Å) ratios are mostly in the range of 1-15 with a median of about 5. Some authors (EO16, and others) assume that the lines are purely collisionally excited, so that the line ratio should depend only on the density and the temperature of the emitting region.On the other hand, Gorti et al. (2011) argue that to produce the observed line ratios by collisional excitation requires either very high densities or high temperatures (T e > 30, 000K).They suggest that the line ratios could more naturally be produced by OH dissociation, but Fang et al. (2018) argue against this by comparison with [SII] 4068 Å. Nemer et al. (2020) point out the importance of FUV pumping for the oxygen lines, an effect that is included in our present calculations.In short, caution is needed when using these line ratios to diagnose physical conditions in the gas.
The line profiles in Fig. 1 show a single peak with a modest blueshift and both broad and narrow width, hence would be classified as NC or BC depending on the veiwing angle.The [OI] 6300 Å comes mostly from lower latitudes and intermediate radii, as suggested by the line shape.On the other hand, the [NeII] line's emission is more extended and samples higher latitudes in the wind where the poloidal velocities are larger.As a result, the [NeII] line has a clear blue tail and higher mean blueshift; if more emission comes from those layers, the blue tail could turn into the BC.The blueshifts are 0−3 km s −1 for all the lines.

Alternate scaling for the wind-base density
Here we present a variant of the magnetothermal wind model in which the density at the base scales as R −5/2 , more steeply than in the fiducial model (∼ R −1.33 ).As noted in the Appendix, the steeper density scaling implies a mass-loss rate d Ṁwind /d ln R 0 ∝ R −1 0 , and an accretion rate through the disk with the same scaling if the wind torque is the main driver of accretion.Such a disk and wind could not be in steady state but would develop a central hole.This might be an mechanism for producing transition disks, but our main purpose here is to explore how a more radially concentrated wind affects the line profiles.
We can see from Fig. 2 that the density has a different structure; notably, the gas is more concentrated toward higher latitudes, where the field lines emanate from footpoints at small radii.The line profiles in Fig. 2 display distinct peaks that could be interpreted as broad and narrow components.To the best of our knowledge, this is the first physical wind model to display both a BC and a NC.As shown below, these components come from different regions in the wind.The blueshifts are 7 − 30 km s −1 for the optical lines and 10 − 80 km s −1 for the [NeII] line; the latter is at the upper limit of the LVC classification in observations (Pascucci et al. 2020).This range of blueshifts in the [NeII] is achieved by a single wind model (magnetothermal 2 ) at various inclinations and is subject to the gaussian fitting of the line in Fig. 2.
By varying the input parameters for these magnetothermal models, we can obtain higher or lower blueshifts, depending on the balance between magnetic and thermal forces.We find that the density at the wind base and the FUV and X-ray radiation sources correlate positively with the blueshift, but the effect is small: a few km s −1 for an order-of-magnitude increase in the radiation source or doubling of the base density.
The blueshift is most sensitive to the assumed Alfven velocity at the wind base.When increasing the Alfven velocity by an order of magnitude from the fiducial value (as an extreme following BYGY), we find that the blueshift increases by an order of magnitude from ∼ 15 − 93 km s −1 which is a significant change and can affect the classification of the line.The blueshifts and line profiles are less sensitive to the sound speed specified in our wind solutions, as long as the latter is reasonably consistent with the temperatures at which the lines of interest are strong: 3-10 km s −1 .(We choose our wind parameters with this in mind, but we have not attempted an iterative reconciliation of the dynamical sound speed with the temperatures computed by Cloudy.This might be achievable but would complicate our wind models considerably, as they currently take the sound speed to be constant along each field line.) To investigate the origins of the two components seen in [OI] 6300 Å, we show in Fig. 3 the contributions to the line profile from restricted ranges of radii within the wind.Evidently, the broader peak is produced mainly within 0.5AU.As we include larger radii, we start to see a second narrower and less blueshifted peak emerging.These findings are consistent with Simon et al. (2016), who found that the BC traces a gas from 0.05 -0.5 AU, while the NC was found to trace the wind beyond 0.5 AU.
Furthermore, we fit multiple gaussians to the full profile at 20 o inclination, for which we find that three gaussian components give the best fit.In Fig. 4 we show these three components plotted against the emission profile from the regions described in Fig. 3 (i.e.Fig. 3b, 3c, and 3d).One might hope that the decomposition of the line profile into gaussian components would represent distinct emission regions in a one-to-one fashion.But as can be seen from Fig. 4, this is not the case.In fact, the issue that stands out the most is that the wing of any of the regional profiles is not well captured by any single component of the the 3 gaussian fit.

Photoevaporative wind
We construct unmagnetized wind models as described in §2.1 for comparison with the photoevaporative winds of Ercolano & Owen (2016); Ballabio et al. (2020);Picogna et (2019).As described in the Appendix, we fix the sound speed at 5.8km s −1 for all radii.Note that the parameters of this model are different from those of the other four, in order to better resemble those of the above-cited works: the streamlines make angle 75 • (vs.60 • ) with the midplane, the wind base lies at z 0 = 0.25R 0 (vs.0.15R 0 ), and the radial extent is 17 AU (vs. 4 AU).The sonic point for the streamlines are extended beyond the emission region of the [OI] 6300 Å and are indicated for a few streamlines The optical line luminosities are on the same order  1).
of magnitude as the fiducial model.The [NeII] line emission, on the other hand, is an order of magnitude higher than in our fiducial model which is more consistent with observations than the fiducial model.The oxygen line ratio is consistent with the observations.
As can be seen from Fig. 5, the line profiles are narrower than in our fiducial model and have negligible blueshift.Lacking magnetic forces, the gas is accelerated to lower velocities (i.e.smaller blueshifts) that exceed the local escape velocity only at larger radii, where rotational broadening is less.Actually, gas escapes from every radius, but the sonic point lies far from the wind base along streamlines launched from the inner disk, so that the density of the escaping gas is exponentially small and contributes negligibly to line emission.The [NeII] line is narrower than the oxygen lines at high inclination, suggesting that it  1).
originates at larger radii, and shows a distinct blueshift.In these photoevaporative models the [OI] 6300 Å emission comes mostly from an almost hydrostatic atmosphere above the disk rather than a full-fledged wind, i.e. from below the sonic point.

Magnetocentrifugal wind
Since the input parameters for such an MHD wind are not well constrained, we have chosen values inspired by the fiducial model of BP82, as did Weber et al. (2020) [model "magnetocentrifugal 1 " in Table 1].
The sound speed is half of that of the fiducial model (BP82 neglected thermal pressure altogether), the Alfven speed is large in order to produce a large Alfvén radius (R A /R 0 = √ 30 in BP82), and the wind is anchored at the mid plane (as for a cold and thin disk); We omit the emission from below the "surface" of the disk as done by Weber et al. (2020).
These choices produce a faster and (in post-processing) somewhat cooler wind than our fiducial magnetothermal model.The line luminosities are on the same order of magnitude as in the fiducial model.
Comparison of the density profile between the magnetothermal and the magnetocentrifugal models suggests that the densities are similar.Actually, we should have decreased the density at the wind base in inverse proportion to (R A /R 0 ) 2 so that the fiducial and magnetocentrifugal winds exerted the same torque on the disk, and hence made for the same accretion rate.This would give a lower wind density and fainter lines (by a factor of 4).b But we chose the current density normalization (and slightly higher Alfvén speed) to produce similar accretion and mass loss rates as Weber et al. (2020) for better comparison.
We find that the oxygen emission comes from the surface of the disk where the gas is not significantly accelerated, similar to our findings for the magnetothermal wind model (Fig. 6).Here the OI6300 is centered at the stellar velocity and represents the NC, while the [NeII] comes from higher layers and exhibits a blue tail.

Alternate density scaling
As with the magnetothermal models, here we explore the effects of as steeper wind-base density (scaling as R −5/2 rather than as R −1.33 as for constant Σ FUV ).With this steeper density scaling (model "magnetocentrifugal 2 " in Table 1), the strong magnetization and large Alfvén radii of the magnetocentrifugal model produce excessive blueshifts and strongly asymmetric line profiles from the higher latitudes, as can be seen from Fig. 7.As a result, those lines have multiple gaussian components, which would likely be interpreted as HVC +BC, although the HVC is often attributed to a jet rather than a disk wind.The BC component could be qualitatively consistent with the observations of Fang et al. (2018), who concluded that this represents a slow disk wind.The highly blueshifted lines of the present model (magnetocentrifugal 2 ) are in agreement with Weber et al. (2020), but the luminosities are an order of magnitude higher than their models, as discussed above.

X-rays and Ne II
The [Ne II] 12.8µm line has been postulated to be a unique signature of X-ray-irradiated protoplanetary disks (Glassgold et al. 2007;Alexander 2008), and subsequent observations of this line seem to favor an X-ray-driven photoevaporative wind Pascucci et al. (2020).In this work, we find that the NeII 12.8µm line traces an outer (up to 20-15AU) disk wind and is strongly correlated with the X-ray radiation source assumed.In particular, the [NeII] line traces an upper wind layer that has sufficient temperatures (∼ 1000K) to excite the fine-structure level.
As explained in Glassgold et al. (2007), most of the singly ionized Neon is formed through X-ray K and L-shell ionization by photons with energies ∼ 1keV, followed by charge exchange with neutral hydrogen.
There are a few calculations done for the rate coefficient of collisional fine-structure excitation of the NeII line; we have consulted the most recent data from (Wan et al. 2019), where we see that the collision strength decreases with increasing temperature, though rather slowly.We created a planar-slab model using Cloudy to investigate the emission process for [NeII] line.We found that the emission is proportional to the X-ray and FUV sources and density (although it saturates for higher densities) while it decreases for temperatures above 1000K.This confirms that the [NeII] line emission is collisionally excited and highly dependant on the ionization of the local environment.

DISCUSSION
This work studies the emission details of the forbidden lines based on the magnetized wind models of Bai et al. (2016).Our aim is to demonstrate that the LVC (both BC and NC) can be produced by magnetized wind models.Current direct numerical simulations of magnetized winds rarely cover a wide enough radial range to exhibit both the BC and the NC.Weber et al. (2020) attempted to model forbidden lines using the self-similar magnetocentrifugal wind of Blandford & Payne (1982); they were able to reproduce the line width and blueshift of the BC but not the line luminosities, nor the NC.We believe that this is because the magnetocentrifugal BP82 solutions, although pioneering and conceptually fundamental, are qualitatively different from the more weakly magnetized and heavily loaded magnetothermal winds that are characteristic of recent theoretical PPD work.
Photoevaporative wind models were shown to be successful in reproducing the observed NC of the LVC, that the line prevails more strongly in transition disks (TD), which (as the name suggests) may represent a later phase of PPD evolution (but see Owen 2016).They believe that this correlation is mainly due to the key role that the X-ray radiation plays in creating enough ionized neon to observe this line emission; the dependence of the [NeII] line on the X-ray source in PPDs was also suggested in the theoretical work of Glassgold et al. (2007) before such a large sample of observations was available.Hence it is believed that the X-ray radiation is less absorbed in the inner regions for disks with inner holes allowing it to reach the upper layer of the wind to excite the NeII line (where the temperature is appropriate for emission).
We wanted to compare line profiles from our fiducial model with other magnetized wind models.Weber on the parameters of the radiation sources, the sound speed of the wind, but are most sensitive to the assumed magnetic field and the density at the wind base.We confirm that the poloidal velocity component is mainly responsible for the blueshift, whereas the width of the line is determined by the local keplerian velocity, which depends on the distance of the emission region from the star.This is clear from the magnetized models, where we see broader, more blueshifted emission that comes from a wind launched in the inner regions accelerated predominantly by the magnetic field.Conversely, the photoevaporative model displays narrower, less blueshifted lines because the emission comes from larger radii where the rotational broadening is less.
Solar mass T-Tauri stars have a wide range of observed accretion luminosities, but the majority of these systems have L acc ∼ 0.01 − 0.1 L sol (Manara et al. 2023;Fang et al. 2018).Usually, the accretion luminosity is not measured directly, because the source is mostly in the UV and is efficiently absorbed before reaching the observer, but inferred from tight correlations with H lines (Fang et al. 2018).France et al.
(2023) directly measures the UV H 2 , Ly α , C IV 1600 Å bump, and FUV continuum between 912 − 1760 Å (the major constituents of the accretion FUV continuum) using HST observations, and confirm the results of previous studies about the typical luminosity of this source.We approximate the shape of the FUV source with a blackbody with T ef f = 12,000K following the previous theoretical work (Ercolano & Owen 2016;Wang et al. 2019;Gorti & Hollenbach 2009).Our choice of L acc = 1.8 L sol , although within observational constraints, is an extreme case according to the observed sample.To better compare with observations , we show in Fig. 8 the relationship between L acc and L OI for a wide range of L acc .We find that for log(L acc /L sun ) ≥ -2, the simulated OI luminosities are -6.2 ≤ log(L OI /L sun ) ≤ -4 which is within the observed range for these lines (Banzatti et al. 2019;Rigliaco et al. 2013;Fang et al. 2018).Though, it seems that the relationship plateaus for log(L acc /L sun ) ≤ -2 as was reported by Ercolano & Owen (2016) in their independent theoretical work perhaps because there is a process that suppresses the OI excitation for lower temperatures which is not included in current models.Moreover, Fang et al. (2018)  It is commonly accepted that the [OI] lines can be broken into broad and narrow components (BC, NC), and that these components trace distinct regions.This conclusion relies on the multi-gaussian fitting procedure used in analyzing the observations.We can see from our models that the lines are highly asymmetric and fitting a gaussian (or multiple gaussians) can be misleading.Moreover, we saw that a simple change in the level of the continuum assumed could change the shape of the fit to the line.The inference of distinct regions based on a decomposition of the observed lines into gaussians has been challenged by many authors who point out that the gaussian fitting does not clearly separate the flux of the overlapping component in the observed spectra (Banzatti et al. 2019), and that the line is inherently asymmetric (Ballabio et al. 2020) and not purely keplerian (McGinnis et al. 2018).The above-mentioned observational limitations could introduce subjectivity to the fitting procedure.In future, we hope to explore whether fitting, both the modelled and observed spectra, with more robust procedures such as machine-learning techniques, could substantially enhance the analysis.
Our models have their limitations and could be made more self-consistent.The sound speed assumed in launching the flow is prescribed and constant along each field line (or streamline), and no attempt is made to reconcile it with the temperature computed by Cloudy in post-processing.As in Bai et al. (2016), the geometry of the magnetic field is prescribed (with two adjustable parameters, θ ′ and q); the flow along each line is computed in steady-state, but force balance across lines is ignored.Ambipolar diffusion and other non-ideal MHD effects have been neglected, though these could be important for launching the wind and for heating it after it is launched.Our models already have several adjustable parameters, however, so that while it will some day be important to make more complete and realistic models for comparison with data, the present state of both observations and theory may not yet justify the effort involved.

CONCLUSION
In this work, we process the analytic wind solution from Bai et al. (2016) through the radiation-transport code Cloudy to obtain simulated emission spectra of PPDs.We then modify the analytic wind solution to produce models very similar to photoevaporative and magnetocentrifugal wind solutions to compare with our fiducial magnetothermal model and discuss the discrepancies.Our aim is to provide insight into current studies of wind tracers, so we focus on common observables such as the blueshift and FWHM in relation to physical properties like the magnetic field and sound speed.
Our simulations suggest that magnetic forces are needed to launch a wind close to the star, where the BC of the oxygen lines forms.Hence, magnetothermal wind models are needed to explain the observed complex spectra.Photoevaporative wind models can explain only the NC, because at temperatures compatible with the forbidden oxygen lines, the sound speed is insufficient to launch a vigorous wind from small radii.
Furthermore, at least in our photoevaporative models, much of the oxygen emission comes from a quasihydrostatic atmosphere above the disk, as reflected by its small blueshifts, whereas the [NeII] better samples the disk wind.
In this work, we have imitated standard observational procedure by fitting gaussian profiles to each of our model lines and declaring the blueshift and FWHM based on these fits.However, as our computed line shapes are often far from gaussian and have extended tails that might be difficult to discern above the continuum in a real observation, the use of gaussian fits may not make the best use of the available data for diagnosing the dynamics and physical conditions of disk winds.

Figure 4 .
Figure 4. Comparing individual components of the multi-gaussian fitting for the full profile in Fig.3against the emission profiles of distinct regions for the magnetothermal 2 model for i=20 (i.e.Fig.3b,3c,3d)

Figure 5 .
Figure 5. Polar plots of emissivity for [OI] 6300 Å and [NeII] 12.81 µm along with their line profiles, temperature and density for photoevaporative wind model.Note that the velocity scale is smaller than for the magnetized models.The sonic point (from the wind base) for a streamline anchored at 0.5AU = 5.7AU, at 5AU = 1.6AU, and at 10AU = 1.4AU well beyond the [OI] 6300 Å emission region but captures some of the [NeII] emission as seen in the line profile.

Figure 7 .
Figure 7. Polar plots of emissivity for [OI] 6300 Å and [NeII] 12.81µm along with their line profiles, temperature and density for the magnetocentrifugal wind model with a different base density scaling (magnetocentrifugal 2 in Table1).

leading
photoevaporative wind models fail to reproduce the BC of the LVC, however, and have some discrepancies with the observed line.Weber et al. (2020) models suffer from much lower luminosities than observed, and double-peaked emission for the BC.Ballabio et al. (2020) models cannot reproduce the observed line widths and blueshifts simultaneously, and they could only model NLVC but not the BLVC.
et al. (2020)  performed the most recent of these simulations, where they relied on the analytic description of an MHD wind based on the work of(Blandford & Payne 1982) axisymmetric self-similar solutions for a magnetocentrifugally driven wind.The problem with these "cold" magnetocentrifugal wind models is that they are too strongly magnetized and too lightly mass-loaded, with the result that the velocities are too large, leading to excessive blueshifts and highly skewed profiles.Our magnetothermal model based on the work ofBai et al. (2016) is successful in reproducing both the BC and NC of LVC of the [OI] lines, and the LVC of the [NeII] line; the latter usually lacks sufficient observational resolution to be decomposed into BC and NC(Pascucci et al. 2020).Our fiducial model represents only one of many plausible choices for the input parameters.The predicted line shapes depend finds the correlation to follow log L OI ∼ 0.6 log L acc − 4.07 similar to our model with a relation log L OI ∼ 0.77 log L acc − 4.74 for log(L F U V /L sun ) ≥ -2.It is evident that the FUV source has a strong connection to the OI emission, and careful treatment of this source, informed by recent observational work, should be included in future theoretical work.

Figure 8 .
Figure 8. [OI] 6300 luminosity vs accretion luminosity relation based on the simulation results of our fiducial model (orange).We overlay a fit to the data points (dashed) which is explained further in the text.

Table 2 .
A table that lists the blueshifts and FWHM of the [OI] 6300 Å line for different inclinations from a single gaussian fitting .HM 0 F W HM 20 F W HM 40 F W HM 60 F W HM 80