Constraining Protoplanetary Disk Winds from Forbidden Line Profiles with Simulation-based Inference

Abstract

Protoplanetary disks are the sites of vigorous hydrodynamic processes, such as accretion and outflows, and ultimately establish the conditions for the formation of planets. The properties of disk outflows are often inferred through analysis of forbidden emission lines. These lines contain multiple overlapping components, tracing different emission regions with different processes that excite them: a high-velocity component (tracing a jet), a broad low-velocity component (tracing inner disk wind), and a narrow low-velocity component (tracing outer disk wind). They are also heavily contaminated by background spectral features. All of these challenges call into question the traditional approach of fitting Gaussian components to the line profiles, and cloud the physical interpretation of those components. We introduce a novel statistical technique to analyze emission lines in protoplanetary disks. Simulation-Based Inference is a computationally efficient machine learning technique that produces posterior distributions of the parameters (e.g. magnetic field, radiation sources, geometry) of a representative wind model when given a spectrum, without any prior assumption about line shapes (e.g. symmetry). In this pathfinder study, we demonstrate that this technique indeed accurately recovers the parameters from simulated spectra without noise and background. A following work will deal with the analysis of observed spectra.

tablenum

1 Introduction

Protoplanetary disks (PPDs) are composed of molecular gas and dust and encircle most low-mass pre-main-sequence stars, referred to as protostars. These disks play a pivotal role in the formation of planets, making them a subject of immense interest in research. Over time, PPDs disperse, typically within a few million years, likely due to a combination of stellar accretion, planet formation, and outflows (Alexander et al., 2014; Ercolano & Pascucci, 2017; Manara et al., 2023). To investigate the characteristics and evolution of PPDs, researchers utilize tracers of disk winds that provide insights into various physical conditions and regions (Güdel et al., 2014). Notably, emission lines from oxygen, sulfur, and neon have emerged as prime candidates for tracing disk winds, enabling the study of their physical properties through accurate modeling of line ratios. The detection of the $[\textrm{Ne}\textsc{ii}]12.8\micron$ fine structure line has confirmed the existence of warm ionized disk winds, often exhibiting a blueshift relative to the stellar velocity, indicating an outflow moving toward the observer, while the emission from the receding outflow in the other half of the disk is obstructed by the dust (Font et al., 2004; Alexander, 2008; Gorti & Hollenbach, 2009; Pascucci & Sterzik, 2009).

Detailed spectroscopic observations have revealed the presence of both fast outflows reaching speeds of hundreds ${\rm\,km\,s^{-1}}$, as well as slower outflows with velocities $\lesssim 10{\rm\,km\,s^{-1}}$ (Hartigan et al., 1995). In certain luminous Class I and II systems, the fast outflows develop into highly collimated jets that extend for parsecs and manifest as Herbig-Haro objects (Bally, 2016; Pascucci et al., 2023, and references therein). The rate at which mass is lost through these outflows is closely linked to the rate of accretion (Hartigan et al., 1995; Ellerbroek et al., 2013; Rigliaco et al., 2013; Natta et al., 2014). The acceleration of these jets is believed to involve magnetohydrodynamic (MHD) mechanisms, although the exact launch locations—whether in close proximity to the star or more broadly from the PPD—remain uncertain (Frank et al., 2014). The lower-velocity component of the outflows is likely associated with the disk itself, supported by evidence such as inclination-dependent line widths, lower blue shifts, and higher estimated mass-outflow rates (Fang et al., 2018; Simon et al., 2016; Banzatti et al., 2019).

Despite the valuable information extracted from these observations, limitations arise from the scarcity of high-resolution infrared spectra, hindering a comprehensive understanding of PPDs, their evolutionary trends, and detailed line profiles. However, optical forbidden lines, including $[\textsc{Oi}]\,6300\textrm{\AA}$ and $[\textsc{Sii}]\,4068\textrm{\AA}$, offer an alternative avenue with abundant high-resolution optical spectra (Hartigan et al., 1995). These lines exhibit distinguishable low-velocity components (LVCs), emission blueshifted by less than 30${\rm\,km\,s^{-1}}$, and high-velocity components (HVCs, Hartigan et al., 1995). Recent studies have successfully modeled LVCs with magnetothermal wind models (Nemer & Goodman, 2024), although recent observations have shown that LVCs can be further deconstructed into broader components (BCs), with full width at half maximum (FWHM) of at least 40 ${\rm\,km\,s^{-1}}$, and narrower components (NCs, Rigliaco et al., 2013; Simon et al., 2016; Fang et al., 2023). This decomposition into different components is usually done with a multi-Gaussian fitting, where the choice of the number and shapes of the Gaussians is chosen by minimizing the residual after fitting (Fang et al., 2018).

The exploration of the relationship between BC and NC FWHMs and the inclination of observed systems suggests that BC likely originates from a wind launched at radii of $0.05$-$0.5{\textsc{au}}$ ($0.5$-$5{\textsc{au}}$ for the NC), assuming Keplerian rotation. From these conclusions about emission regions, it was suggested that the NC traces an outer photoevaporative wind, while the BC has an MHD wind signature (Simon et al., 2016). However, it is more likely that the emission originates from a range of radii spanning a range of temperatures and densities, and the observed components are clearly overlapped with no explicit separation into distinct components tracing distinct regions. Furthermore, the observed spectra are commonly contaminated with background spectral features that have to be processed before fitting the lines with any model which introduces further difficulty in analyzing the optical spectra (Banzatti et al., 2019). These issues put the Gaussian fitting procedure into question as it was suggested that analysis of the emission lines should be done in more robust way that is not based on the symmetry of the lines (Ballabio et al., 2020).

In this paper, we introduce a Simulation-based Inference (SBI) approach to extract spectrum information on the basis of state-of-the-art wind models. With an innovative combination of convolutional neural networks and neural density estimation techniques we can interpret forbidden emission line spectra of PPDs, without the limitations of Gaussian profile fitting. We test the efficacy of our approach with simulated spectra in this paper and will employ the same methods to analyze a sample of observed spectra in a subsequent paper.

2 Methods

We use the analytic wind models of Nemer & Goodman (2024) to generate simulated spectra, then compress them with a neural spectrum encoder, and train a neural density estimator to carry out SBI as the mapping from model parameters to (compressed) spectra. Afterwards, we test the trained SBI network with a sample of simulated spectra (with known parameters) that it has not seen during training. We summarize the different steps below.

2.1 The wind model

In accordance with the methodology outlined by (Bai et al., 2016, hereafter BYGY), we undertake the solution of the steady-state, axisymmetric ideal magnetohydrodynamic (MHD) equations governing fluid flow along each poloidal magnetic field line, under the condition of a specified density $\rho_{0}$ and poloidal Alfvén speed $V_{\textsc{a}0}=B_{\textsc{p}0}$ at the base of the wind. These winds possess a finite sound speed $c_{\rm s}$, which remains spatially invariant along each magnetic field line with radial variation across lines. Termed “magnetothermal”, these winds are characterized by their acceleration resulting from a confluence of magnetic forces (represented by $V_{\textsc{a}}$) and gas pressure (represented by $c_{\rm s}$). The configuration of the poloidal magnetic field is delineated by straight field lines that intersect the surface of the accretion disk at an inclination angle $\theta^{\prime}$ from the disk midplane, with the specific intersection points denoted by cylindrical coordinates $(R_{0},z_{0})$. Anticipating these intersections to lie several scale heights above the midplane, where FUV radiation from the star can penetrate and ionize the gas sufficiently to couple it with the field, we set $z_{0}=0.15R_{0}$ to approximate this condition. The dynamics of the flow along each magnetic field line are governed by four dimensionless input parameters $(\theta^{\prime},z_{0}/R_{0},q,V_{\textsc{a}0}/V_{\textsc{k}0},c_{\rm s}/V_{\textsc{k}0})$ and three dimensional scaling factors $(R_{0},\rho_{0},V_{\textsc{k}0})$, with the latter representing the Keplerian orbital velocity at the base of the wind, $V_{\textsc{k}0}=\sqrt{GM_{*}/R_{0}}$. Notably, for $q<1$ the scaling of the poloidal magnetic field strength with radius smoothly transitions from $B\propto R^{-1}$ to $B\propto R^{-2}$ at $R\gtrsim R_{0}/q$.

2.2 Postprocessing with Cloudy

We utilize the radiative transfer software Cloudy (Ferland et al., 2017) to analyze the density and velocity distributions within our models of stellar winds. This analysis, coupled with appropriate energetic photon luminosities, allows us to compute the emission of the [Oi] 6300Å line, following the methodology outlined in Nemer & Goodman (2024). Our focus lies predominantly on the inner-to-intermediate regions of the disk, spanning distances from (0.04-4 au), as these regions contribute the majority of the emission observed in the line. The computations are conducted on a spherical grid, with the radial cell size being variable and automatically adjusted by Cloudy to optimize the accuracy of the calculations. The grid comprises 361 points uniformly distributed in the $\sin\theta\in[0,1]$ coordinate space within the range [0,1]. Under the assumption of the disk midplane being optically thick, we exclude the lower half of the disk from our analysis. To determine the luminosities of the emission lines, we integrate their emissivities along radial paths. Each such path is associated with its unique density profile and is integrated over one hemisphere $(0\leq\theta\leq\pi/2)$, a condition pertinent to optically thin emission lines.

We construct a three-dimensional axisymmetric volume and conduct computations to derive the spectral profiles of our lines under various lines of sight, corresponding to different inclinations. The resultant line profile is intimately linked to the local gas temperature and flow velocity within the volume. The emitted radiation from each cell undergoes broadening based on the local gas temperature, and its distribution across velocity bins is determined by the projections of the flow onto the observer’s line of sight at a specified inclination angle. Aggregating the contributions from all velocity bins yields the composite line profile, which is subsequently adjusted to conform to the calculated emission luminosity. To address instrumental effects, we incorporate convolution with a Gaussian function characterized by a standard deviation of $\sigma=8{\rm\,km\,s^{-1}}$, mirroring the resolution capabilities of observational instruments.

In total, we generate 27,860 simulated [Oi] 6300Å spectra from sampling 9 input parameters (as described in the following section and shown in Fig. 1): the Alfvén velocity $v_{A}$, the sound speed $c_{s}$, the geometry of the poloidal field lines $\theta_{B}$, the inner radius of the disk $R_{\rm in}$, the luminosity of various radiation sources $L_{\rm stell},L_{\rm FUV},L_{\rm EUV},L_{\rm Xray}$, and inclination of the system $inc$.

2.3 Simulation-Based Inference using Normalizing Flows

Our goal is to infer the posterior distribution of the 9 model parameters (PPD attributes), $\theta=\{v_{A},c_{s},\theta_{B},R_{\rm in},L_{\rm stell},L_{\rm FUV},L_{\rm EUV},L_{\rm Xray},inc\}$, given a spectrum of forbidden emission lines, $x$: $p(\theta\,|\,x)$. With the conventional Bayesian probabilistic inference approach, the posterior would be evaluated using Bayes’ rule: $p(\theta\,|\,x)\propto p(\theta)\,p(x\,|\,\theta)$. $p(\theta)$ is the prior distribution and $p(x\,|\,\theta)$ is the likelihood, which is often evaluated assuming a Gaussian functional form. Then, the posterior is inferred by sampling it using Markov Chain Monte Carlo (MCMC) methods. Sampling a 9-dimensional posterior distribution with MCMC requires many ($\gg 10,000$) evaluations of the posterior. Each evaluation in our case involves running a full forward model that requires 400-500 CPU hours. This makes the conventional inference approach computationally infeasible.

SBI offers an alternative approach. The latest SBI methods use neural density estimators (NDEs) to directly approximate the posterior using a training dataset, $\{(x_{i},\theta_{i})\}$ (e.g., Alsing et al., 2019; Wong et al., 2020; Dax et al., 2021; Zhang et al., 2021; Hahn et al., 2022). With NDEs, we can efficiently and accurately estimate high-dimensional posteriors with just a few thousands of training data (e.g., Hahn et al., 2022). This approach has additional benefits. SBI does not make any assumptions about the functional form of the likelihood. Furthermore, once the neural posterior estimator is trained, generating the posterior for a new spectrum requires no additional model evaluations and training.

Below we describe our SBI approach in detail. We follow the SBI approach in Hahn & Melchior (2022) with an added data compression/feature extraction step that reduces the dimensionality of the spectrum and improves the overall performance of the SBI model.

2.3.1 Prior Distribution

In SBI, the prior distribution is set by the distribution of the parameters in our training data. In our case, the prior distribution is not uniform as there are physical limitations that prohibit certain combinations of input parameters. The choice of the prior is based on the exploratory work by Nemer & Goodman (2024) done with this wind model:

• For the radiation sources, $L_{i}$, the central value of the prior distribution was found to reproduce the average [Oi] 6300Å observed line luminosity, ratio to [Oi] 5577Å, and the shape of the LVC. The width of the prior distribution for these sources is chosen to mimic the spread in the quantities as reported by observations (Rigliaco et al., 2013).

• The sound speed, $c_{s}$, chosen to reflect a temperature $\sim 5,000-10,000K$ in the line emitting region (EO16) since the temperature is prescribed to be constant along the field lines in the wind solution.

• The inclination, $inc$, is chosen to be uniform since there is no specific preferred inclination for PPDs.

• The inner radius of the model, $R_{\rm in}$, is chosen to be bimodal to reflect both classical PPD systems, with a primordial disk, and disks with inner cavities where there is little to no emission from inner radii.

• Disk magnetic field, on the other hand, has no solid measurements, so the magnitude and the geometry are somewhat ambiguous. Hence, $v_{A}$, which is a proxy for the field magnitude, is chosen within a range that would produce reasonable line width in accordance with observations; we did sample outlier field values that would produce extreme line widths as this has also been seen in observations (Banzatti et al., 2019).

• Poloidal field lines parallel, or close to parallel, to the rotation axis or the midplane can only accelerate the gas in extended regions ($>5{\textsc{au}}$) which is outside the expected emission region (Fang et al., 2023). So, the choice of field angles $\theta_{B}$ is constrained between $20^{o}$ and $70^{o}$

We present the distribution of parameters in the training data, i.e., our prior distribution (blue), in Fig. 1.

Figure 1: The distribution of model parameters, $\theta=\{v_{A},c_{s},\theta_{B},R_{\rm in},L_{\rm stell},L_{\rm FUV},L_{\rm EUV},L_{\rm Xray},inc\}$. The blue contours show the distribution of physical parameters used to generate the training spectra. The red contours instead show the latent space after compressing the spectra with our trained encoder Spender. The agreement between the two distributions illustrates that the encoder successfully extracts the essential information about the parameters from the full spectrum.

2.3.2 Data Compression

The simulated spectrum has 4,000 spectral elements centered at the central wavelength. However, the success of past works in decomposing the spectra into Gaussian components strongly suggest that it can be efficiently expressed in a much lower dimensional space, without the loss of significant information. In this work, rather than using Gaussian decomposition, we use the convolutional encoder from Spender (Melchior et al., 2023; Liang et al., 2023b, a), an autoencoder neural network architecture designed to extract compact representation of spectra. With the Spender encoder, we can efficiently extract relevant features of the spectra while also relaxing any assumptions on the shape of the different components.

Spender employs a convolutional encoder $g(x)$ with an attention layer to identify the most prominent patterns or “features” in the spectrum $x$. The input spectrum of the encoder first goes through three convolutional layers, which account for correlations among the dominant features and segment the spectrum in wavelength. The attention mechanism then identifies the most influential spectral features, which are finally passed through a fully connected multilayer perceptron to produce the compressed spectrum in latent parameter space. In Melchior et al. (2023), the compressed spectrum is passed through a decoder and the Spender autoencoder is trained by minimizing the reconstruction loss between the input spectrum and the output of the decoder.

In this work, we only use the Spender encoder without invoking the decoder. We train it to predict the input model parameter values used to generate the simulated spectrum. More specifically, we minimize the L2 norm between the input and predicted model parameters. As a result, our compression reduces the spectrum from 4,000 spectral elements to $N=9$ dimensions, equivalent to the number of input parameters. In other words, the encoder acts as a data compression and initial parameter regression step. By training our encoder to predict the model parameters, we design the encoder to extract the features with the most constraining power on the model parameters. A similar data compression approach was used in Chen et al. (2023); Lemos et al. (2023).

To train the encoder, we first reserve 15% of the data for testing. Then, we split the rest of the data into an 80% training and 20% validation sets. We train the encoder for 100 epochs with the Adam optimizer (Kingma & Ba, 2014), the 1Cycle schedule (Smith & Topin, 2017) with a maximum learning rate of $5\times 10^{-3}$, and a batch size of 128. On an NVIDIA V100 GPU, the training takes 5 minutes approximately. Our final trained encoder has comparable training and validation losses. This suggests that the network is not overfitting and we conclude that the training procedure is overall stable. As shown in Fig. 1, the predicted distribution of latent parameter/compressed spectrum for the test sample (red) is similar to the input parameters of the training sample (blue). This demonstrates that the encoder successfully compresses the spectrum into a lower dimensional latent space that encodes the important features in the spectrum.

2.3.3 Simulation Based Inference

From the compressed spectrum $s=g(x)$, we obtained the parameter posteriors $p(\theta\,|\,s)$ using SBI based on “normalizing flows” (e.g., Tabak & Vanden-Eijnden, 2010; Tabak & Turner, 2013; Kobyzev et al., 2019). In short, normalizing flow models utilize an invertible bijective transformation $f$ to map a complex target distribution to a simpler base distribution. This transformation needs to be invertible and have a tractable Jacobian for evaluation of the target distribution from the base distribution. A neural network, with parameters $\phi$, is trained to represent $f$. Among various normalizing flow models (e.g., Germain et al., 2015; Durkan et al., 2019), we use Masked Autoregressive Flow (MAF; Papamakarios et al., 2017) models implemented in Python package sbi (Tejero-Cantero et al., 2020), similar to Hahn & Melchior (2022).

Our goal is to train a MAF model, $q$, with hyperparameters $\phi$, that accurately estimates the posterior, i.e. the conditional probability distribution of the parameters given the (compressed) data: $p(\theta\,|\,s)\approx q_{\phi}(\theta\,|\,s)$. We do this by minimizing the KL divergence between $p(\theta\,|\,s)$ and $q_{\phi}(\theta\,|\,s)$, which is equivalent to maximizing the total log likelihood, $\sum_{i}\log q_{\phi}(\theta_{i}\,|\,s_{i})$ over the training set. In practice, we first split the training data into a training and validation set with a 90/10 split. Then we use the Adam optimizer with a learning rate of $5\times 10^{-4}$ to maximize the log likelihood. We use early stopping to prevent overfitting: we evaluate the total log likelihood on the validation data at every training epoch and stop the training when the validation log likelihood fails to increase after 20 epochs. We determine the architecture of our MAF (e.g., number of blocks, hidden layers) through experimentation and selecting the model that produces the highest validation log likelihood. Our final MAF model has 3 Masked Autoencoder for Distribution Estimation (MADE; Germain et al., 2015) blocks, each with 2 hidden layers and 80 hidden units.

To apply SBI directly to observations, we require simulated spectra with realistic noise and observational systematics. Available PPD observations are usually processed for background spectral features (including photospheric and telluric absorption lines) utilizing stellar templates of known systems and precise measurements of the system’s radial velocity, as outlined by Banzatti et al. (2019); Fang et al. (2018). At present, however, we lack sufficient information about these processing steps and their uncertainties to include their effects on our mock spectra. Consequently, in this paper, we do not add noise and background to the simulated spectra and focus on demonstrating the overall feasibility of our SBI framework for analyzing forbidden line profiles. In the subsequent paper, when we apply our framework to observations, we will supplement our model with a realistic background spectral features model that is tailored to the specific observations.

3 Results

Figure 2: The posterior estimate $q_{\phi}(\theta\,|\,x)$ for a randomly selected test spectrum. The inferred posterior using SBI is consistent with the truth parameter values (blue).

To demonstrate how the posterior of a spectrum is estimated using SBI, we show in Fig. 2 the posterior estimate $q_{\phi}(\theta\,|\,x_{i})$ for a randomly selected test spectrum. We overlay the true input parameters of the test spectrum (vertical and horizontal blue lines) for reference. The estimated $q_{\phi}(\theta\,|\,x)$ is in excellent agreement with the truth. The spread in the distribution is a consequence of the loss of information, either in the mapping from parameters to spectra or in the encoding from spectra to initial parameters. To our knowledge, this is the first time a method is presented to robustly connect wind model parameters to observed line shapes without the need for parametric fitting with Gaussians. Each off-diagonal panel presents the marginalized 2D posterior of different parameter pairs, where the contours represent the 68, 95, and 99 percentiles. The diagonal panels present the marginalized 1D posterior of each parameter with the 16 and 84 percentiles (dashed) marked.

To further test whether $q_{\phi}(\theta\,|\,x)$ robustly approximates the true posterior $p(\theta\,|\,x)$, we take the maximum a posteriori (MAP) of the estimated posterior and perform forward modeling on it to predict the corresponding spectrum. We then examine how closely the MAP simulated spectrum reproduces the original spectrum for a given sample. This suggests that the parameters, in which SBI places the highest confidence, adequately represent the original spectrum. To estimate the MAP, we rank the posterior samples by their log probability and select the parameters with the highest value. In Fig. 3, we present the comparison between the line profiles of the simulated spectra (red) and the MAP spectra (green). We find good agreement between the spectra despite their complex, non-Gaussian shape. It is noteworthy that three of the modelled spectra in Fig. 3 show a slight double peak, while one shows an extreme double-peaked profile; a similar feature was reported in previous modelling work (Weber et al., 2020). The double-peaked profile is due to the rotational motion of the gas, and is expected from any system that exhibits Keplerian motion, but is rarely observed in PPDs. The slightly-double-peaked profiles are expected to be smoothed out when noise and background features are added. Nevertheless, even if the predicted double-peaked emission does not represent any real situation, it helps to have these out-of-sample predictions to avoid overfitting the neural networks with biases from our prior knowledge of the observations.

Figure 3: Comparison between the simulated spectrum (red) and the forward-modeled spectrum of the MAP estimate from $q_{\phi}$ (green). They are in overall good agreement with each other.

Next, we validate the accuracy of $q_{\phi}$ using simulation-based calibration (SBC, Talts et al., 2018). SBC examines the distribution of the rank statistics of the true parameter values within the marginalized posteriors. In Fig. 4, we present the SBC of each model parameter for the normalizing flow posteriors (blue) using 4179 test samples. The shape of this distribution can tell us about the quality of predicted posteriors compared with the true ones, including whether it is biased and whether it is over- or under-confident. Any asymmetry in the distribution implies that the posterior estimates are biased. For example, a U-shaped distribution where the true parameter values are more often at the lowest and highest ranks means that the posterior estimates are narrower than the true posteriors. We do not find these features for any of the model parameters in Fig. 4. Hence, SBI provides unbiased and accurate estimates of the true posteriors (black dashed) for all model parameters.

Figure 4: SBC plot of the model posteriors for 4179 simulated spectra. The histogram in each panel represents the distribution of the rank statistic of the true value within the marginalized model posterior (blue) for each parameter. For the true posterior, the rank statistics will have a uniform distribution (black dashed). The rank statistic distribution of the model is nearly uniform for all of the parameters.

4 Summary and Discussion

The motivation behind the analysis done in this paper is to provide a more capable framework for analyzing optical observations of PPDs. The line profiles and ratios provide key information about the emitting region which helps understand the physical properties of the system at different evolutionary stages. Conventionally, observations of optical forbidden emission lines were analyzed using Gaussian fitting (Fang et al., 2018; Banzatti et al., 2019). Notably, observations exhibit multiple components that are assumed to trace distinct emission regions.

Current models of PPDs indeed confirm the existence of multiple emission regions, and as a result, multiple components (Weber et al., 2020). However, our models show that the lines are typically asymmetric, broad, and trace a range of physical properties (temperature, density, radiation). Moreover, the observed spectra are usually contaminated with background spectral features (photospheric and telluric absorption lines) that need to be corrected for before analyzing the optical line profiles. This step involves accurately determining the location of these absorption lines, by correctly measuring the system’s radial velocity, and then using stellar templates of known systems to subtract the unwanted features. This step is crucial in the analysis because the recovered line shape and shift need to be as precise as possible for correct analysis. In addition, the the different components of the line usually overlap which introduces further uncertainty when disentangling them (Banzatti et al., 2019).

We present a novel method using machine learning to interpret the observations, which eliminates the subjectivity inherent in the current methods including Gaussian fitting. The asymmetry and overlapping components in the line shape introduce a degeneracy in the Gaussian fitting, and potentially valuable information could be lost with the Gaussian assumption. The asymmetry is an indication that emission comes from wide spread regions besides the regions where bulk emission comes from, and that part needs to be included in the models for a sound analysis of the observations. Also, the extraction of distinct components from overlapping emission could be tricky and subjective, especially if there are multiple asymmetric constituents.

One significant benefit of employing machine learning in this analysis lies in mitigating the distortion caused by background spectral features and noise. Addressing this point in detail will be postponed to a future paper, given that publicly available observations typically have these contaminants corrected for, and reintroducing them to the spectra requires meticulous handling. However, our approach promises to better quantify the relationship between uncertainties in the wind parameters and uncertainties in the background features such as the systemic velocity, stellar and atmospheric spectral templates, etc.

The fidelity of the analysis provided here is largely sensitive to the validity of the model used to describe PPD systems. The radiation sources used in the radiative transfer simulations are generic, but PPD systems in reality span a wide range of stellar, far ultraviolet, and X-ray sources depending on the type of the system. Using generic radiation sources proved useful in interpreting observations in the past (Fang et al., 2023), and further refinement of the sources to specific systems would provide better constraints on the observed systems. We also do not include jets in our model of disk winds, and those are suspected to be responsible for one of the overlapping emission components (the HVC). Yet, we can apply the model to disks that do not display the high velocity component and those represent about half of the available sample of PPDs (Fang et al., 2018; Banzatti et al., 2019). Nonetheless, our SBI scheme is sufficiently modular that other wind models could be used if such models became available.

We show in this paper that modern machine learning-based techniques are effective in analyzing PPD emission spectra. They offer robust analysis of complex line shapes beyond the simple Gaussian assumption. With the emergence of sophisticated models for emission spectra, SBI becomes increasingly compelling to extract parameter constraints in a computationally efficient way. We show the successful application of SBI to the analysis of simulated PPD spectra. This work shall be extended to further analyze actual observations and other complicated systems.