An Outbursting Protostar: The Environment of L1251 VLA 6

Young protostars that undergo episodic accretion can provide insight into the impact on their circumstellar environments while matter is accreted from the disk onto the protostar. IRAS 22343+7501 is a four-component protostar system with one of those being a fading outbursting protostar, which was obseved with the Karl G. Jansky Very Large Array (VLA) and is referred to as L1251 VLA 6. Given the rarity of young stellar objects undergoing this type of accretion, L1251 VLA 6 can elucidate the fading phase of the post-outburst process. Here, we examine structure in the disk around L1251 VLA 6 at frequencies of 10 and 33 GHz with the VLA. We model the disk structure using Markov Chain Monte Carlo (MCMC). This method is then combined with a parametric ray-tracing code to generate synthetic model images of an axisymmetric disk, allowing us to characterize the radial distribution of dust in the system. The results of our MCMC fit show that the most probable values for the mass and radius are consistent with values typical of Class I objects. We find that the total mass of the disk is 0.070−0.2+0.031M☉ and investigate the conditions that could cause the accretion outburst. We conclude that the eruption is not caused by gravitational instability and consider alternative explanations and trigger mechanisms.


INTRODUCTION
The star formation process begins in the dense regions of interstellar molecular clouds in which a gravitational collapse causes gas and dust to heat up and begin forming young stellar objects (YSOs).The rotation of the collapsing cloud core produces an accretion disk, which can add mass to the protostar.Mass buildup can be accelerated by episodic accretion in YSOs.Theoretical considerations (e.g., Vorobyov & Basu 2015) suggest that Sun-like stars build up their mass during periods of significantly enhanced accretion (i.e., episodic accretion).The most substantial outbursts manifest during the embedded phase, wherein the accretion disk is supplied with gas originating from the collapsing envelope.Episodic accretion can play a crucial role in the evolution of circumstellar disks, and it is considered to have an impact on planet formation.Outbursts have been shown to potentially change the chemistry and mineralogy of the surrounding circumstellar disk ( Ábrahám et al. 2009;Rab et al. 2017); could spur the growth of small solids through, e.g., evaporation and recondensation from a rapid evolution of the ice line (Cieza et al. 2018); and may lead to the large spread of observed protostellar luminosities (Fischer et al. 2023).Furthermore, in the common model of star formation, YSOs surrounded by circumstellar and/or circumbinary disks regularly generate strong, episodic jets that enable interpretation of the evolutionary state and activity of those YSOs (Arce et al. 2007;Pech et al. 2010;Vorobyov et al. 2018).
While it is predicted that all Sun-like stars undergo several outbursts during their pre-main sequence evolution (Audard et al. 2014), there are very few detected young stars that currently experience accretion outbursts (Fischer et al. 2023).This is in large part due to the phenomenon becoming observable only as the circumstellar envelope thins (Audard et al. 2014).Given that there are only a few dozen stars of this sort that are known, it is crucial for our understanding to examine each new item (Semkov et al. 2021).Among those first discovered were the original prototypes FU Orionis and EX Lupi (Herbig 1989).Subsequent to Audard et al. (2014), even more outbursting sources have been found, as the Gaia satellite has been remarkably successful in discovering new, outbursting sources (e.g., Hillenbrand et al. 2018;Szegedi-Elek et al. 2020;Hodapp et al. 2019;Hillenbrand et al. 2019;Hodapp et al. 2020;Cruz-Sáenz de Miera et al. 2022;Nagy et al. 2023;Siwak et al. 2023), making this an important phenomenon to study more broadly.These discoveries have expanded our understanding of this phenomenon, and it is now known that outbursts occur in various YSOs, encompassing different classes and intermediate/ambiguous properties.Mechanisms of mass accumulation and disk heating, leading to accretion bursts, are discussed in several theoretical papers (e.g., Bell & Lin 1994;Bell et al. 1995;Vorobyov & Basu 2006;Zhu et al. 2009;Bae et al. 2014).The existence of an envelope is presupposed in most scenarios, therefore FUors may signify the change from embedded to optically visible young stars.In one of the proposed mechanisms, episodic outbursts are caused by matter accumulation in the inner disk region as a result of angular momentum rearrangements in the outer disk caused by gravitational instabilities (GI), as proposed in Zhu et al. (2009).The inner disk heating causes rising ionization, which in turn initiates instabilities, including magnetorotational and thermal, which propel the gas to accrete from the disk to the star (Bell & Lin 1994).The interactions of the circumstellar disk with a planet or nearby stellar partner on an eccentric orbit might be another potential triggering mechanism (Lodato & Clarke 2004;Reipurth et al. 2004;Dong et al. 2022).
FU Orionis type stars, which we will refer to as FUors, can exhibit a rapid increase in magnitude (∆V > 5) lasting 10s of years (Audard et al. 2014;Fischer et al. 2023).EX Lupi type stars, which we will refer to as EXors, have smaller and shorter lived outbursts.The increased optical and near infrared magnitudes of these objects are due to an increase in the accretion rate.The increased accretion heats a substantial part of the circumstellar disk, and in general the whole circumstellar matter, evidenced by the brightening of large-scale reflection nebulae, e.g.HBC 722 in Miller et al. (2011).Then, when observed at multiple wavelengths, the brightening is seen in a very broad wavelength range from the optical to the far-infrared, e.g.V1647 Ori in Muzerolle et al. (2005), meaning a significant increase in the total bolometric luminosity as well.In order to build an accurate profile of outbursting protostars, detailed observations with broad spectral coverage, spatial scales, and covering the epochs of before, during, and after outburst are imperative.
The question still remains on whether all young Sun-like stars undergo eruptive phases or not.To answer this, it is important to study the properties of outbursting disks.For instance, the disk mass (M disk ) is a particularly relevant parameter regarding whether or not the disk can be gravitationally unstable.The M disk for eruptive stars have been attempted to be measured, such as in Liu et al. (2018); Cieza et al. (2018); Kóspál et al. (2021).An ALMA survey presented in Kóspál et al. (2021) found that FUor disks are more compact and more massive than disks of so-called "regular" (non-eruptive) Class I or Class II sources.For reliable mass measurements, there needs to be long wavelength data/observations because of optical depth effects as disks may be optically thick even at millimeter wavelengths.VLA, operating at cm wavelengths, is the perfect instrument for this.The observed emission from X band data is likely the mixture of thermal and free-free emission, and observations in multiple bands and/or high spatial resolution are necessary to disentangle the different contributions.
Here, we present VLA observations of 2MASS J22352345+7517076, also known as L1251 VLA 6 and hereafter referred to as VLA 6, an embedded eruptive YSO belonging to the five-component protostellar system IRAS 22343+7501.Embedded in the Lynds 1251 molecular cloud, the source's distance is 350 +46 −38 pc as presented in Kun et al. (2019).Rosvick & Davidge (1995) identified a cluster of five near-infrared sources associated with IRAS 22343+7501 (RD95 A, B, C, D, and E).The cluster is a source of several molecular outflows (Sato & Fukui 1989;Nikolić et al. 2003;Kim et al. 2015), the Herbig-Haro jet HH 149 (Balazs et al. 1992), and radio continuum jet sources (VLA 6, 7, and 10, Reipurth et al. 2004).VLA 6 exists within a four component system with stars RD95A referred to as VLA 7, RD95B and RD95C referred to as VLA 10 with the young outbursting protostar referred to as VLA 6 and is identical with RD95D.We also present observations of VLA 6 that were conducted using the Faint Object infraRed CAmera (FORCAST) with the Stratospheric Observatory for Infrared Astronomy (SOFIA) for greater insight into whether the source is declining at wavelengths of 25.3 and 31.5 µm.Onozato et al. (2015) reported the dramatic brightening of VLA 6 at mid-IR wavelengths by comparing IRAS, Akari, and WISE data.The history of the source's luminosity evolution is clearly outlined in Kun et al. (2019), and especially evident in their Figure 6.Briefly, over the course of a 35-year period, the flux tripled in the first 6 years    2019) utilized an SED, corrected for interstellar extinction, and assumed the bolometric luminosity in the low-state is dominated by the central star to determine a mass range of 1.6 − 2.0 M ⊙ for VLA 6 based on a bolometric luminosity that increased from 32 L ⊙ to 165 L ⊙ between 1983 and 2010 and its proximity to the 10 5 -year isochrone.The peak accretion rate, also from Kun et al. (2019), was estimated from their accretion disk model to be slightly above 10 −4 M ⊙ yr −1 , typical of an FUor-type outburst.The exact cause of the outburst, and indeed the underlying mechanism of outbursts in general, is not well constrained.The IRAS 22343+7501 cluster presents a valuable opportunity to observe a young four-component star system and the underlying effects these stars potentially have on each other, one of which is experiencing an accretion outburst.Given the rarity of YSOs undergoing this type of accretion, VLA 6 can provide insight into the fading post-outburst process.
In this paper, we describe the new observational data in Section 2, look at results in Section 3, examine disk parameters and constraints in Section 4, discuss our results in Section 5, and present our summary and conclusion in Section 6.

OBSERVATIONS
The IRAS 22343+7501 cluster was observed with the VLA.We collected data at 33 GHz data on 2019 November 08 (ID 19B-088, PI White) and at 10 GHz on 2020 December 10 (ID 20B-096, PI White).An overview of the observations is given in Table 1.

Ka Band Data
The 33 GHz observations were centered on VLA 6 using J2000 coordinates RA = 22 h 35 min 23.46 s and δ = +75 • 17 ′ 7. ′′ 60.The data were taken using 27 antennas in the D configuration with baselines ranging from 0.035 − 1.03 km.The total time on-source was 831 s.
The observations used a Ka Band tuning setup with 4 × 2.048 GHz basebands and rest frequency centers of 28.976 GHz,31.024 GHz,34.976 GHz,and 37.024 GHz.The data yield an effective frequency of 33 GHz (9.1 mm).The raw data was reduced and calibrated using the CASA 6.1.2-7pipeline (CASA Team et al. 2022), which included bandpass, flux, and phase calibrations.Quasar J0019+7327 was used for phase and bandpass calibrations and quasar 0542+498=3C147 was used as a flux calibration source.The uncertainty of the absolute flux calibration of the VLA at these frequencies is typically ∼ 5%.The pipeline calibrated data was inspected in CASA and additional flagging was performed to account for radio frequency interference (RFI).All the significantly impacted channels and time bins were flagged and the pipeline was ran again.These included data from specific antennas such as ea06, ea14, ea26, and spectral windows 34-49 for antenna ea17.All of the RFI was individually removed using the CASA task flagdata.
The final calibrated dataset was imaged with a natural weighting and cleaned using CASA's CLEAN algorithm down to a threshold of 1.5 × σ RM S , presented in Figure 1.The observations achieve a sensitivity of 19 µJy beam −1 as measured in the CLEAN ed image.The size of the resulting synthesized beam is 3. ′′ 94 × 2. ′′ 07 at a position angle of −68.3 • .

X Band Data
The 10 GHz observations were centered on VLA 6 using RA = 22 h 35 min 24.34 s and δ = +75 • 17 ′ 9. ′′ 20.The data were taken using 24 antennas in the A configuration with baselines ranging from 0.68 − 36.4 km.The total time on-source was 1792 s.The observations used rest frequency centers of 8.999 GHz and 10.999 GHz.The data yield an effective frequency of 10 GHz (30 mm).
The raw data was reduced and calibrated using the CASA 6.1.2-7pipeline (CASA Team et al. 2022), which included bandpass, flux, and phase calibrations.The pipeline calibrated data was inspected in CASA and additional flagging was performed to account for RFI.All the significantly impacted channels and time bins were flagged and the pipeline was run again.The continuum image is presented in Figure 2. The observations achieve a sensitivity of 100 µJy beam −1 as measured in the CLEAN ed image.The size of the resulting synthesized beam is 0. ′′ 54 × 0. ′′ 26 at a position angle of −60.23 • .

SOFIA FORCAST Data
Observations of VLA 6 were conducted using FORCAST with the SOFIA telescope on August 31, 2018 at wavelengths of 25.3 and 31.5 µm (PLANID: 06 0165, PI Huard).Data reduction and calibration were performed by the SOFIA team using the SOFIA data-reduction pipeline (Herter et al. 2018).
Level 4 images, which were artifact-corrected and flux-calibrated in units of Jy/pixel, were downloaded from the SOFIA archive.Aperture photometry was performed on the images using a 12-pixel aperture and a sky annulus between 20-25 pixels.The aperture includes the neighboring source VLA 7, and to correct for its contribution to the total flux measured in the aperture, the peak fluxes of both sources were compared and the flux ratio was used to correct the composite flux measured in the large aperture.These observations can be used to further refine our understanding of the varying nature of fading at different wavelengths of VLA 6.

RESULTS
Figure 1 and Figure 2 show naturally weighted images of the combined Ka and X band data sets generated using the CASA task tclean.In the Ka band, three sources that correspond to sources VLA 7, VLA 10, and VLA 6 are detected, while in the X band, only VLA 6 is detected.We used CASA's imfit tool to derive basic parameters such as coordinates, sizes, inclinations, position angles, and fluxes.This enabled us to find best-fit 2D Gaussian models for each source component, which are presented in Table 2.It was not possible to determine a precise size for VLA 10 as it is measured to be nearly a point source.The peak flux of VLA 6 is detected with the Ka band at a 120 sigma level and with the X band at a 44 sigma level.The radii were estimated from the deconvolved best-fitted 2D Gaussian, we took the FWHM of the major axis and divided it by 2(2ln2) to obtain σ, which we use as a proxy for the radius of the disk emission at the specific wavelength based on their primary beam corrected images.At millimeter wavelengths, protoplanetary disks are usually optically thin with the exception of the core regions of particularly massive systems (Andrews et al. 2009(Andrews et al. , 2010)).This means that the majority of dust grains are responsible for the observed emission, and the total flux is closely related to the total mass of small grains.Under this assumption, the following formula can be used to determine the dust masses using millimeter fluxes:

Disk Masses
where d is the distance to the target, F ν is the flux density, κ ν is the dust opacity, and B ν (T ) is the Planck function evaluated at a temperature of T (e.g., Andrews et al. 2013).We incorporate a distance of 350 +46 −38 pc based on the average parallax of the 15 known optically visible members of the L1251 star forming region included in Gaia DR2 as presented in Kun et al. (2019), set the frequency to the 33 GHz measurement, and make standard assumptions about the dust opacity (κ ν = 10(ν/1000 GHz) cm 2 g −1 ; Beckwith et al. 1990) and beta = 1.For the temperature, we used different values to check how its variation affects our mass estimate.In order to obtain the total disk mass, we assume a 100:1 gas to dust ratio where we multiply our M dust by 100.We obtained total disk masses of 0.092, 0.055, 0.040, and 0.031 M ⊙ using dust temperatures of 30, 50, 70, and 90 K, respectively.
Following from Equation 1, we can instead look at disk mass as we vary the internal density, ρ, instead of κ ν .We know that if 1/κ ν = 4 3 sρ, and Ω s = π s 2 d 2 , then we get the following alternate equation: where ρ is the internal density, Ω s is the solid angle of a single grain, and s is the size as presented in White et al. (2016), where here s is set to 9.08 mm.Similarly with Equation 1, we obtain a disk mass by assuming a 100:1 gas-todust ratio where we multiply our M dust by 100.For grain size, it is most common to assume the size is equal to the wavelength of the observations as blackbody grains emit most efficiently at wavelengths that are about the same as their size.We use a range of different values for ρ including 1000, 1500, 2000, 2500, and 3000 kg/m 3 .Table 3 presents our results.
By comparing the results from Equations 1 and 2, we are able to determine how varying the density and temperature, in addition to other parameters, can affect the potential mass of the disk.It is apparent that varying temperature and density can result in M disk values varying by a factor of 2. It is evident that for the Equation 1 value results, the M disk decreases as the temperature increases.The same scenario similarly appears in Table 3 in addition to the M disk increasing as the density increases.Kim et al. (2015) found the total circumstellar mass to be 0.033 M ⊙ (adjusting for distance to that of 350 pc assumed here) based on SMA 1.3 mm dust continuum observations, which is at the lower end of our 95% confidence interval.Now that we have been able to determine M disk , we are able to begin generating disk models in the RADMC-3D radiative transfer code, discussed in 4.1, that can be compared to the images.

Spectral Index
We measured the spectral index α of our target, VLA 6, by constructing comparisons from the two continuum spectral windows separately with ν 1 and ν 2 , and measuring the targets' fluxes (F) at these two frequencies.Then, we calculated α using: and its uncertainties σ with: We calculated indices using data measured with the JCMT/SCUBA-2 in 2014, (Pattle et al. 2017), and SMA on 2007 Oct 17 (Kim et al. 2015).The indices are also calculated for our 33 GHz Ka band and 10 GHz X band observations.We found the spectral index from the 666 GHz micron JCMT data point to 231 GHz is 5.670 ± 0.018, 231 GHz to 33 GHz is 2.018 ± 0.053, and 33 GHz to 10 GHz is −0.65 ± 0.23.We are able to see how the spectral index using data at 33 GHz or above is positive, while the 33-10 GHz slope is negative.This could be the result of jet emissions that contaminates the flux measurement / contributes to the measured flux at lower wavelengths.With the flux measured at 33 GHz, and assuming that this is purely disk thermal emission, we extrapolate the flux to 10 GHz using an adopted spectral slope (either −3 for blackbody or −4.7 for ISM-like grains).This gives a predicted flux for 10 GHz that can be compared to the noise level of the 10 GHz image.In this case, we see that the 10 GHz flux in the VLA image is higher than what the SED predicts, which we will discuss in Section 5.4.

ANALYSIS
In this section, we analyze the emission of VLA 6 to characterize in detail the spatial distribution of dust in the system.We adopt a modeling approach for the 33 GHz data specifically taking into account the protoplanetary disk density profile equations laid out in Andrews et al. (2009) and then follow the same technique in White et al. (2020) that uses MCMC to constrain the most probable values.This model is then combined with a parametric ray-tracing code to generate synthetic model images of an axisymmetric disk with an MCMC fitting algorithm, allowing us to characterize the radial distribution of dust in the system.

Radiative Transfer Modeling
We use radiative transfer model fitting, as in Dullemond et al. (2012), in order to constrain the disk parameters of VLA 6.This technique gives us a more accurate model than fitting a simple Gaussian to the image shown in Figures 1 and 2. The disk structure model and radiative transfer calculations are defined on a spatial grid in spherical coordinates {r, ϕ}, as described in Andrews et al. (2009).To obtain the physical parameters of the spatially resolved disk, we use the RADMC-3D code (Dullemond et al. 2012) to build radiative transfer models with the Python interface radmc3dPy1 .
The density profile of a trial protoplanetary disk model is given by: ρ(r, ϕ) = Σ(r, ϕ) where Σ is the surface density profile, H p is the pressure scale height, and z is the height above the disk midplane.The surface density profile includes a power-law inner disk and an exponential out tapering (Andrews et al. 2009): where R c is the characteristic radius of the disk, γ is the power-law exponent of the radial surface density profile, and Σ c is the surface density normalization at the inner radius.
The pressure scale height is defined as: where h c is the ratio of the pressure scale height over radius at R c and ψ is the degree of flaring for the disk.These steps follow White et al. (2020).
To generate realistic dust absorption and scattering properties, we followed the same approach as White et al. ( 2020) and generated the RADMC-3D input opacity files with the OpacityTool2 program (Toon & Ackerman 1981;Woitke et al. 2016).This program calculates dust opacities by using a volume mixture of 60% amorphous silicates (e.g., Dorschner et al. 1995), 15% amorphous carbon (e.g.Zubko et al. 1996), and a 25% porosity.Bruggeman mixing is used to calculate an effective refractory index, and a distribution of hollow spheres with a maximum hollow ratio of 0.8 (Min et al. 2005) is included to avoid Mie theory scattering artifacts.The disk is populated by 0.1 − 15000 µm grains following a power-law size distribution of s −3.5 .Finally, the total disk mass is calculated assuming a gas-to-dust mass ratio of 100:1.

Modeling Formalism
To converge on the best fit model, we use Metropolis-Hastings MCMC model fitting, similar to that used in White et al. (2020).We varied eight free parameters: stellar radius (R star ), the total disk mass (M disk ), degree of disk flaring (Ψ), characteristic radius (R c ), power law exponent of the surface density profile (γ), and scale height ratio (h c ).The motivation for varying the stellar radius is that we wanted to allow the model to explore potential disk properties that also consider a range of stellar radii.Setting a fixed parameter of 9.0, as hypothesized by Kun et al. (2019), leads to parameter values that do not converge.A model is created for each individual non-point source, which for this system are VLA 6 and source VLA 7, given their locations in right ascension and declination.
Given the detailed calculation in RADMC-3D, a continuum image, which we will call the "trial sky model", is then produced and projected to a trial inclination and position angle for each source.Note-a Rstar is a fixed parameter for VLA 7, so there are no 95% confidence regions.
The model is first run on VLA 6 with a trial sky model being the result.The trial sky model of the individual source is then combined with the sky models of other sources to create a combined sky model.The point sources for protostars VLA 10 are included at this time, and a trial model for VLA 7 with given stand-in fixed parameters instead of free parameters is created.
The combined sky model is then attenuated by the primary beam and convolved with the synthetic beam for the observational setup in order to later compare the convolved model image with the data image.The convolution consists of a series of Fast Fourier Transform calculations.
To assess the likelihood of a given trial convolved model, a χ 2 is calculated as where σ is the observed σ rms for a given observation multiplied by the synthetic beam size in pixels (see Booth et al. 2016).
A model is accepted if a random number drawn from a uniform distribution [0,1] is less than α, where This process is then repeated for VLA 7 where the model runs with the free parameters to create a new trial sky model.The trial sky model of this individual source is then combined with the sky models of other sources.The point sources for protostars VLA 10 are included at this time, and the trial model for VLA 6 with the previously determined best-fit parameters instead of free parameters is created.
For both VLA 6 and VLA 7, we ran 100 chains with 1,000 links each beyond the burn-in period of 100 links.The most probable values are summarized in Table 5.The 95% confidence regions are also presented, indicated from the posterior distributions in Figure 8.These confidence intervals provide uncertainties for the MCMC fit and correspond roughly to 2σ uncertainties.Figure 4 shows the data image (left) compared with the synthesized best-fit model (center right) and the residuals (right).The lack of residuals after the best-fit model is subtracted from the data verifies a good fit a posteriori.

DISCUSSION
We found that the resulting total disk mass of 0.070 M ⊙ from the MCMC fitting in comparison to the dust calculations in Section 3.1 is accurately reflected.The disk mass value is similar to the resulting value when temperature is set to 50 K in Equation 1.The value is also similarly consistent when temperature is 50 K and density is 3,000 kg/m 3 in Equation 2. This value also closely resembles to when the temperature is 30 K and density is 2,000 kg/m 3 .

Classification
In Figure 5, we show disk masses as a function of characteristic disk radii for FUors and non-eruptive Class I and Class II sources, reproduced from Kóspál et al. (2021).The results of our MCMC fit show that the most probable values for the mass and radius are more consistent with Class I as seen in Figure 5.In comparison with VLA 7, the eruptive star has less massive, smaller disk with flatter surface density profile, higher scale height, and more flaring based on the parameters in Table 5.However, the confidence intervals largely overlap, implying no significant difference between these two disks.
The results presented in Kóspál et al. (2021) suggest that FUor disks differ from Class I/II disks in the distribution of their masses and radii.The FUors presented have a median FUor disk radius of 27 au and median FUor M disk of 0.08 M ⊙ .The Class I disks that are plotted in blue and were sampled from Sheehan & Eisner (2014, 2017) have a median disk radius of 127 au and median M disk of 0.02 M ⊙ .Class II comparisons were drawn from Andrews et al. (2010), and show a median disk radius of 40 au and median M disk of 0.03 M ⊙ .
From these trends, it is evident that FUor disks tend to be more massive and compact than regular Class I/II disks.However, we find that with a characteristic disk radius of 260 au, VLA 6 is larger for its disk mass of 0.07 M ⊙ , coinciding more with regular Class I disks, rather than FUor disks.We also find that the mass-radius ratio of VLA 6 and VLA 7 are similar.It is interesting to note that despite these similarities, VLA 7 is not an outbursting source like VLA 6.Thus, it appears that the mass-radius ratio is not a property that distinguishes FUor-like protostars from non-outbursting YSOs in all cases.This result is in contrast with Kóspál et al. (2021) that raised the possibility that FUor eruptions were specific to a Class of YSO's with a unique compact disk structure that undergoes a different evolutionary path.The larger radius of VLA 6, instead, suggests that there may be a population of outbursting stars that deviate from these assumed FUor properties.We propose that these deviations may be dependent on the cause of outburst (as discussed in Section 5.2).
All FUor and Class II sources presented in Figure 5 were modeled similarly with RADMC-3D and their parameters were constrained with MCMC fitting as our own modeling for VLA 6.The data in Kóspál et al. (2021) used 1.3 mm observations while we present here 9 mm data, therefore, there will be a potential difference in the obtained disk mass or disk radius.The Class I disks were also modeled with RADMC-3D though did not constrain the parameters with the same approach, using a power law surface density profile without exponential taper.Though the modeling formalism differs from that of the FUor samples, Kóspál et al. (2021) found that setting the characteristic radius as the outer disk radius, R disk from Sheehan & Eisner (2017), produced similar spectral profiles and nearly identical total disk masses.Given the lack of studies on Class I disks with the same model as presented in this paper, these coinciding results make the parameters from Sheehan & Eisner (2017) most comparable to that of the FUor samples.
Historically, young outbursting protostars were classified by well defined categories.Observed diversity, as described in Fischer et al. (2023), shows that this is no longer necessarily the case.Therefore, it is possible there are many young outbursting stars that are similar to VLA 6 but are yet to be discovered and labeled as such.According to the classification scheme set by Connelley & Reipurth (2018), VLA 6 is classified as a "peculiar" object that was similar but distinct from classical FUors; in the following discussion we explore several additional characteristics of VLA 6 in the context of other FUors for reference, yet the exact designation should not impede the interpretation of this particular source.

Gravitational Instability / Causes of Outburst
The larger and small-scale structure of FUor disks can reveal the source of the episodic accretion that causes the outbursts.While the exact mechanisms that drive increased disk material accretion onto the protostar are unknown, it could be a result of the disk becoming gravitationally unstable, possibly in conjunction with a magnetorotational instability.GI is one possible mechanism for outbursts (Vorobyov & Basu 2006;Zhu et al. 2009;Kuffmeier et al. 2018).Kóspál et al. (2021) checked the criterion for GI in a sample of FUor disks and found that about 2/3 of them may be gravitationally unstable.There are multiple metrics to calculate the GI, the more detailed of which is Toomre's Q parameter.GI occurs if the value of Toomre's Q parameter is Q = c s Ω/πGΣ ≤ 1.4,where c s is the sound speed, Ω is the epicyclic frequency for a Keplerian disk, and Σ is the column density at a given radius in the disk.Based on our best-fitting radiative transfer model, we calculated the Toomre Q parameter and the results are presented in Figure 6.
We found that VLA 6 is gravitationally stable, consistent with the other 1/3 of the FUor disks that were examined in Kóspál et al. (2021).This contrasts with predictions in Kun et al. (2019) that suggested GI was the cause of a 10 −3 M ⊙ clump accretion onto the star within 10 au.Appearances of GI earlier on in the history of the outbursting process may be considered.Present results of there being no evidence of GI implies that there are other underlying causes for the enhanced accretion.
According to Hartmann & Kenyon (1996), mass is added onto the central star so rapidly that the magnetosphere which may launch an X wind is likely crushed out of existence, implying that the strong outflows that still occur during these outbursts must be caused by an extended disk wind.Several molecular outflows originate from the cluster (Sato & Fukui 1989;Nikolić et al. 2003;Kim et al. 2015); Balázs et al. (2004) reports on the Herbig-Haro jet HH 149, and radio continuum jet sources have been observed in VLA 6, 7, and 10 by Reipurth et al. (2004).Close encounters in multiple systems and flybys have been shown to trigger outbursts (Borchert et al. 2022).As predicted by models, interactions between the stars and their disks for close binaries cause accretion pulses to occur (Artymowicz & Lubow 1996).Recent findings indicate that for less massive disks (M disk ≤ 0.1M ⊙ ), close stellar encounters and thermal instability are more likely causes of outbursts, whereas disk fragmentation, GI, and magnetorotational instability are more mass-dependent mechanisms (Cieza et al. 2018).The low disk mass of 0.07 M ⊙ may explain why the mass-radius ratio of VLA 6 deviates from that of the FUors sampled in Kóspál et al. (2021) as seen in Figure 5.Further investigation with higher angular resolution will have to be conducted in order to determine the cause of the outburst in this system.

Disk Flaring
Figure 7 shows a 2D slice of the model disk's density structure for VLA 6 as a function of radius.The disk is flared due to the outburst heating up the disk which increases the vertical mixing and therefore making the disk thicker.The black contour lines correspond to the disk temperature.These move outward from the outburst as the central effective temperature increases.

Possible Jet Emission
We can confidently determine that the X band detection coincides with the outbursting VLA 6 source, looking at the contours in Figure 1.The X band data can give insight into emission observed from episodic jets.At this frequency, free-free emission is most likely from an ionized jet (Galván-Madrid et al. 2015).Other possibilities for emission include non-thermal emission that may come from the disk, disk wind, or stellar wind.Given the distance, 1 mJy in the X band may be too bright to be dust emission as the dust mass is too high.The X band flux is higher than what the SED modeling in Fig 3 would imply at that frequency, and it is more likely to consider that free-free emission is obscured by dust at high-frequency.We use our best-fit radiative transfer model to produce the spectral energy distribution presented in Figure 3 (blue line) from our MCMC results.The flux value predicted by this model for a disk with 0.070 M ⊙ mass would produce a thermal emission of 60 µJy at 10 GHz, whereas what we observe is about 1200 µJy, producing a factor that can go up to 60 higher.
Based on our radiative transfer modeling of the disk emission, the thermal disk flux is negligible in the 8-12 GHz range.At 10 GHz, the model predicts a thermal flux of 60 µJy, but we measure 1158 µJy, so the thermal flux is probably only about 5% of the measured total flux at 10 GHz.Therefore, the whole range of free-free emission is closely constrained by the measured X band flux density at any given moment (Liu et al. 2018).The flux densities at higher frequencies that cannot be explained by free-free emission must be assigned to dust thermal emission, as these are well fitted by our radiative transfer model in Figure 3.

Figure 1 .
Figure 1.Continuum image of the Ka band resolved with VLA 7, VLA 10, and VLA 6 visible.White contours represent X band emission at a 4 σ level, which clearly coincides with the source VLA 6.The white ellipse in the lower left corner represents the size and orientation of the synthesized beam for the Ka band image.

Figure 2 .
Figure 2. Continuum image of the X band marginally resolved with white contours around emission at a 4 σ level.The white ellipse in the lower left corner represents the size and orientation of the synthesized beam.

Figure 3 .
Figure 3. Spectral energy distribution displaying values presented in Table 4.The blue line represents the spectrum as performed by the radmc3dPy package from best-fit model.

Figure 4 .
Figure 4. (Left) Naturally weighted Ka band VLA image of the IRAS 22343+7501 system.(Center) Convolved model image showing the most-probable model.(Right) Residual image after subtracting the convolved model from the data.

Figure 5 .
Figure 5.Total M disk (in M⊙) as a function of characteristic radii.The orange points are the FUors presented in Pérez et al. (2010); Pérez et al. (2020); Cieza et al. (2018); Kóspál et al. (2021).The black points are T Tauri disks from Andrews et al. (2010).The blue points are Class I disks from Sheehan & Eisner (2014, 2017).FU Ori, VLA 6, and VLA 7 are labeled.The linear fits are presented through the solid lines, where M disk ∝ R 1.33 c represents the FUors, M disk ∝ R 0.99 c represents Class I and M disk ∝ R 1.40 c represents Class II.

Figure 6 .
Figure 6.GI plot with the upper limit of Toomre's Q parameter of 1.4 for VLA 6 represented by the black line.The calculations for VLA 6 are represented by the blue line.

Figure 7 .
Figure 7. 2D slice of the model disk's density structure for VLA 6.

Figure 8 .
Figure8.Posterior distributions for the MCMC chains for VLA 6 after removing burn-in.Histograms for each parameter are plotted at the top of the corresponding column, while the plots for the rest of the grid show the distribution of walkers across slices through parameter space for the corresponding pair of parameters.The best-fit value for each parameter is plotted with a blue line.