Models of Rotating Infall for the B335 Protostar

The PACS observations were obtained by the DIGIT (Dust, Ice, and Gas in Time) key program (Green et al. 2013, 2016). The SPIRE observations were obtained by the COPS (CO in Protostars) open-time program (Yang et al. 2018). The instruments, observational methods, and data reduction have been described in detail in those papers. For B335, the spectra produced by the improved processing (Yang et al. 2018) agreed with photometry on images obtained from the Herschel archive to within 35% at 500 μm, and less than 24% at shorter wavelengths (Yang et al. 2018). The PACS spectrum was used to define a photometric point at 70 μm, while PACS photometry was used for longer wavelengths.

We used Spitzer IRS data from the archive, but did a new reduction. The field is dominated by scattered light in the outflow lobes at short wavelengths (Stutz et al. 2008), making source extraction and sky subtraction difficult. Previous reductions (Yang et al. 2018) showed artifacts and were not representative of the source. The Spitzer IRS spectrum was produced using a mixed method that we describe here. We began with the basic calibrated data products from the Spitzer Heritage Archive. The Short-Low (SL2; 5–7.5 μm, SL1; 7.5–14 μm), Short-High (SH; 10–20 μm), and Long-High (LH; 20–36 μm) spectra were observed under AOR 3567360, using a raster mapping technique (Watson et al. 2004) in which six positions were observed surrounding the source position—a pair of observations at the two nominal nod positions separated by ∼1/3 of a slit length, a second pair +1/2 the slit width in the perpendicular direction, and a third pair −1/2 the slit width in the perpendicular direction. This strategy was intended to guarantee that the source position would be observed even if the observation grid was not centered on the source. In this case, the source was clearly detected in four of the six raster positions for SL2 (at +6 and −4 pixels from center for nod 1 and 2, respectively), SL1 (at +6 and −4 pixels from center), and LH (at +1 and −1 pixels from center). The source was not detected in any of the raster positions for SH, so we do not use those data. The SL observations were taken with a slit oriented about 14° west of north, so they were perpendicular to the outflow.

We attempted to use local and more distant sky to subtract excess emission in SL1 and SL2, for which the slit is sufficiently long to identify sky, but found even the off-positions to be dominated by substantial diffuse emission that reduced the spectral signal-to-noise ratio. Furthermore, as we could not remove sky background from LH, we decided instead not to remove background. The result will show the local shape of the SED well, but may need to be scaled in absolute flux density.

We describe the individual extraction methods for each module, using the SMART package (Lebouteiller et al. 2010). For SL2 and SL1, we extracted 6 pixel wide (108) "fixed-column" apertures centered on the source peak in each nod in each of the four positions in which the source was detected. We did not perform any sky subtraction. We averaged the four spectra by module, removed the order edges, and produced a single combined SL2+SL1 spectrum. The SL1 aperture was 37 by 108 (6 pixels at 18 per pixel) and SL2 was 36 by 108, while LH was 111 by 223. For LH, we extracted the full aperture (225 wide), and clearly detected the source peak centered at a position consistent with SL1/SL2 in four of the six raster positions. We averaged those four position spectra, removed the order edges, and produced a single combined LH spectrum. This spectrum was scaled down by a factor of 2 to match the photometric point at 24 μm. We used the scaled LH data to create a pseudo-photometric point at 35 μm, where no photometry was available, and we used the SL1 data to create such a point at 9.7 μm. These are definitely uncertain, but provide useful constraints.

The J2000 positions of the peak were 19:37:00.7, +7:34:11 for SL2 (λ ∼ 7 μm) and 19:37:01.0, +7:34:08 for LH (λ ∼ 29 μm). Both are consistent with the ALMA peak (Section 2.4) within the 3'' uncertainties in the Spitzer peak position.

The resulting SED is shown in Figure 1, where we plot only data from before MJD 57000. Almost all these data were taken before the luminosity increase, so we refer to this set as the quiescent SED. The continuous lines are spectrophotometric data from Spitzer IRS, with the current reduction (Section 2.2), along with the data from Herschel-PACS, and Herschel-SPIRE (Yang et al. 2018). The black points with error bars are photometric data, listed in Table 1. The integration over the SED yields an observed luminosity of 1.36 L_⊙, higher than estimated by Evans et al. (2015), due to the new distance and details of data chosen to constrain the integration. We favored the calibrated spectrophotometry over the photometric data, where the wavelengths overlapped, because of the better sampling of wavelengths. In general, the two agree well. If we use the higher values of photometry where they disagree, the observed luminosity increases to 1.76 L_⊙, but we consider that as a strong upper limit. Because of the edge-on geometry, the actual luminosity needed to match the observations is higher (Section 3.3), and the final determination of the actual luminosity depends only on matching the modeled and observed photometric points.

Zoom In Zoom Out Reset image size

Figure 2 zooms in on the spectral region from 5 to 14 μm. The spectrum is noisy but quite remarkable. The silicate feature is quite shallow for such an embedded source, while the ice features at 6.0 and 6.8 μm are very deep. The deepest absorption is around 9.0 μm, a significantly shorter wavelength than the usual maximum absorption by silicate grains. Finally, there seems to be broad absorption over the 13–14 μm region of the water libration mode. This spectrum, while uncertain, is different from all the spectra in Boogert et al. (2008). The fact that we could not do sky subtraction may affect the spectrum. In particular, the shallow silicate feature could arise in part from scattered silicate emission. Improved spectra should be available from approved observations with JWST.

Zoom In Zoom Out Reset image size

In addition to the SED, we use the radial profiles of emission at 450 and 850 μm (Shirley et al. 2000), which are sensitive to the radial variation in the density and hence the evolutionary stage. These profiles provide the strongest constraint on the age of the system.

ALMA data were obtained as part of Cycle 1 project 2012.1.00346.S (PI N. Evans) on 2014 April 27 (UT) in the C32-3 configuration of ALMA Early Science. The details of the observations and data reduction have been described by Evans et al. (2015). In that paper, the primary spectral line targets (HCN, H¹³CN, and HCO⁺ J = 4 → 3 lines and CS J = 7 → 6 line) were presented and the HCN and HCO⁺ lines were used to detect infalling motions. The clean beam for Cycle 1 was 490 × 440 mas at PA = −769.

ALMA data with higher spatial resolution were obtained with the same frequency setups in Cycle 3 (Table 3) as part of project 2015.1.00169.S and Cycle 4 2016.1.00069.S (PI. N. Evans). Six execution blocks were obtained, four in 2016, which did not meet the noise specifications, and two in 2018; when combined, the data met the noise specifications. The properties of the observations are listed in Table 2, including number of antennas, minimum and maximum baselines, and precipitable water vapor (PWV), along with notes on issues. The properties of the lines observed are given in Table 3. The spectral resolution was 122.07 kHz, except for H¹³CN, for which it was 488.28 kHz. The equivalent values of velocity resolution are listed in Table 3.

The data were calibrated and imaged using CASA (McMullin et al. 2007) version 5.4.0. Observations of strong lines were subject to a calibration error when all these data were obtained, which could artificially weaken the lines. The effect is roughly equal to the single-dish line temperature divided by the system temperature. The single-dish lines (Evans et al. 2005) have temperatures less than 3 K, versus typical system temperatures during the ALMA observations of about 180 K, so the effect of calibration error should be less than 1.7%, which is negligible given other uncertainties.

After learning of the outburst, we separated the 2016 data from the 2018 data and reduced them each separately; we now refer to them as Cycle3-16 and Cycle3-18. The delay between observations was extremely fortuitous because the lines in 2018 were much stronger, consistent with the brightening seen in the WISE data. The description of the data reduction below applies to both sets of data, except when noted.

Images were made with 0025 cells. Briggs weighting with a robust parameter of 0.5 was used and the task tclean used auto-multithresh for the clean mask, typically with 2–3 major cycles. As explained in Evans et al. (2015), removal of the continuum in the uv-plane can cause distortions in spectra with both emission and absorption. The lines and continuum were separated after cleaning with the task imcontsub, and spectra were extracted using the task specflux with the region defined to be the size and orientation of the continuum source before deconvolution. The continuum was subtracted in the image plane after careful specification of the line-free channels. A separate reduction with the continuum removed in the uv-plane produced nearly identical results, but with a poorer baseline removal. Self-calibration was used on the continuum and applied to the lines. The clean beam for the 2016 data after self-calibration was 200 × 150 mas at PA = −741. The clean beam for the 2018 data after self-calibration was 240 ×190 mas at PA = −631. The sensitivity, denoted K in kelvin-(Jy bm)⁻¹ and computed from the beam and frequency for each line, is given in Table 3 for the 2018 data.

The continuum image from the Cycle 1 observations was presented in Evans et al. (2015), showing a compact central source and faint arcs from the outflow cavity. The Cycle 3 data with better spatial resolution showed many lines of complex organic species, which are presented in a separate paper. They and the main target lines were avoided in constructing the continuum image in Figure 3.

Zoom In Zoom Out Reset image size

The J2000 peak position for continuum emission, averaged between 2016 and 2018, is 19:37:00.8967, +07:34:09.517 with a difference of only 50 mas, consistent within uncertainties with the Cycle 1 peak position. The deconvolved FWHM major and minor axes, peak intensity, and total flux based on an elliptical Gaussian fit are given in Table 4. In addition, we tabulate the total flux density within a 2'' circular area, representing the flux density arising from a region of diameter 329 au.

Table 4. ALMA Source Properties

Data set	MJD	${\theta }_{\max }$	${\theta }_{\min }$	PA	Iν	Sν (fit)	Sν (2'' )	Species
		(mas)	(mas)	(degree)	(mJy bm−1)	(mJy)	(mJy)
Cycle1	56774	491 ± 49	258 ± 59	9.6 ± 10.1	55.1 ± 2.7	93.5 ± 6.7	110 ± 12	Cont.
Cycle3-16	57586	100.1 ± 7.2	97.3 ± 8.3	43 ± 79	89.6 ± 1.0	118.5 ± 2.0	160 ± 8	Cont.
Cycle3-18	58379	186.0 ± 19	179 ± 22	175 ± 79	75.3 ± 2.0	129.8 ± 5.1	217.5 ± 8.6	Cont.
Cycle3-18	58379	789 ± 54	467 ± 34	72.1 ± 5.3	2.75 ± 0.16	24.8 ± 1.6	28.82 ± 0.63	HCO+
Cycle3-18	58379	259.3 ± 9.6	228.3 ± 34	103 ± 16	6.45 ± 0.11	14.51 ± 0.33	16.83 ± 0.86	HCN
Cycle3-18	58379	188.6 ± 10.7	156.0 ± 8.8	122 ± 17	6.04 ± 0.11	9.64 ± 0.27	8.04 ± 0.72	H13CN
Cycle3-18	58379	305 ± 26	238 ± 19	119 ± 19	1.658 ± 0.072	3.89 ± 0.23	5.31 ± 0.24	CS

Note. For the lines, the units are mJy-km s⁻¹ or mJy bm⁻¹ km s⁻¹.

Download table as:  ASCIITypeset image

The spectra of the four lines (Figure 4) from the 2016 and 2018 ALMA data, in an 02 diameter aperture before removal of the continuum in the image plane, show that the absorption of the continuum by the HCO⁺ and HCN lines is essentially complete and clearly shifted to the red of the source velocity of 8.3 km s⁻¹(Evans et al. 2005). The HCN transition has hyperfine structure (Ahrens et al. 2002; Mullins et al. 2016), but 96% of the line strength lies in three hyperfine components within 0.14 km s⁻¹ of the center frequency that we use to assign velocities. There are two very weak components (2% each) at +1.34 km s⁻¹and −1.67 km s⁻¹. H¹³CN also has hyperfine structure (Fuchs et al. 2004); again, 96% of the line strength lies within 0.11 km s⁻¹of our center frequency, with very weak (2%) features at +1.39 km s⁻¹and −1.72 km s⁻¹. The peak at velocities below the blue dip in the HCN spectrum is at 6.6 km s⁻¹, close to the offset expected for the very weak hyperfine component, but the corresponding velocity for the other weak hyperfine component (9.6 km s⁻¹) corresponds only to the shoulder to the blue of the red dip, not to the secondary peak. No dips or secondary peaks are visible in the H¹³CN spectrum, but some shoulders could plausibly be related to the weak hyperfine components.

Zoom In Zoom Out Reset image size

The CS line shows a local dip near 8.3 km s⁻¹, but not absorption below the continuum at the source velocity, while the H¹³CN peaks at the source velocity. The rise at high velocities in the H¹³CN spectrum is caused by lines of COMs, which also partly contribute to the H¹³CN line profile, along with a possible contribution from a line of SO₂ at 356755.19 MHz. The smoother line profile of the H¹³CN observations results from the lower spectral resolution chosen for that band. There are additional dips in the HCN spectrum at 6.9 and 9.9 km s⁻¹.

For comparison to models, the continuum was removed, choosing line-free channels to fit the continuum. First-order baselines were used, except for H¹³CN and HCN, for which there were too few line-free regions, so a zeroth-order baseline was used. We also removed the continuum in the uv-plane for images of the extended emission, and these images are used for some purposes.

The intensity, integrated over the line, is shown in Figure 5. The presence of other lines required tuning of the velocity interval, as indicated in the caption. The integrated intensity maps were fitted to Gaussians, and the fitted peak positions were the same for all lines to within 004, or 20% of the beam size. The average peak position of the lines agreed with that of the continuum to 0025, about 10% of the beam size. The fits were all compact, but they had different sizes, as indicated in Table 4. Except for HCO⁺, they are compact ellipses. The integrated intensity map of HCO⁺ shows extensions likely to be the edges of the outflow cavities in addition to the compact ellipse.

Zoom In Zoom Out Reset image size

In addition to the spectra toward the continuum source, the ALMA Cycle 3 data reveal some interesting features about the outflow cavity. Partial arcs can be seen in several lines/velocities, but the outflow is prominent in neither the integrated intensity images nor in channels centered on outflow velocities. The most complete image of the cavity walls is that of the CS J = 7 → 6 transition at v = 8.24 km s⁻¹, as shown in Figure 6. We interpret these arcs as the limb-brightened top and bottom of the cavity walls, because the velocity is close to the rest velocity and both blue and red lobes are clearly visible. The structure in Figure 6 is similar to the polarized continuum intensity at 0.87 mm in Figure 1 of Yen et al. (2020) but shows the west lobe more clearly. This structure is visible in all three epochs, but it is most prominent in the 2018 data.

Zoom In Zoom Out Reset image size

In the southwestern cavity wall, there is a bright knot visible in Figure 6. The HCN spectrum taken at this position (Figure 7) shows a long tail of emission extending about 17 km s⁻¹ to the red, so we refer to it as the "red blob." This feature is again seen in all three epochs. None of the other lines show any unusual line profiles at this position. The position of the red blob is 55 = 910 au SW of the protostar. An ejection traveling at 100 km s⁻¹ would take about 46 yr to travel this distance.

Zoom In Zoom Out Reset image size

The light curve of the W2 (4.5 μm) flux densities in Figure 8, along with vertical lines indicating the dates of the Spitzer, Herschel, and selected ALMA observations, implies that the luminosity of B335 increased substantially after about MJD 56000 and has since declined. The Spitzer data were obtained before the original WISE data were available, but the Herschel data were obtained close in time to the original WISE data point, clearly preceding the luminosity rise. We will assume that the Herschel data were also in what we will call the quiescent phase and produce a source model for that phase. Then, we will vary the luminosity in that model to predict the W2 data. After fitting the result, we will have a model for the luminosity as a function of time for the period covered by the W2 data.

Zoom In Zoom Out Reset image size

The continuum emission from B335 has been extensively modeled (Shirley et al. 2002; Stutz et al. 2008; Shirley et al. 2011), but those predating 2008 assumed a distance of 250 pc. Spherical models of inside-out collapse (Shu 1977) are characterized by a sound speed (c_s,eff) and an age (t_col) since the collapse began. For practical purposes inner and outer radii of the envelope are also assumed. Preliminary versions of these models for the ALMA data were discussed by Evans et al.(2015).

Models of the continuum were compared to photometry (Table 1) and radial profiles of the 450 and 850 μm emission (Shirley et al. 2002). These 1D models all had two serious problems. Very young ages (t_col < 1 × 10⁴ yr) were required to approximate the radial profiles, but these ages were incompatible with both outflow ages and ages based on modeling line emission. Also, 1D models severely underestimate the observed SED at shorter wavelengths, a well-known problem because such models have too much extinction at short wavelengths (Shirley et al. 2002). A 1D model of a spherical envelope for t_col = 5 × 10⁴ yr is compared to the observations in Figure 9. While the model can reproduce most of the far-infrared and submillimeter data, it severely underproduces the near-infrared data, and the radial profiles are far too flat to match the observations (Figure 10).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

More realistic models must consider the three-dimensional structure of the source, especially the role of the outflow (Dunham et al. 2010; Yang et al. 2017). We adopted a model that includes an infalling, rotating envelope, an embedded disk, and bipolar outflow lobes. We use the radiative transport code, HYPERION, an open-source, parallelized three-dimensional code (Robitaille 2011). The convergence criteria for HYPERION are applied to the absorption rate ($\dot{A}$) of the dust in each cell. The convergence criteria for most simulations were set so that 95% of the differences between iterations were less than a factor of 2, and the 95% value of the differences changed by less than 1.02 (about 2%).

Our modeling approach follows that of Yang et al. (2017), who explained the basic method and explored the effect of the different parameters on the observed SED. We summarize the input model construction here. All the parameters are summarized in Table 5. The envelope is a rotating, infalling envelope that was originally spherical and isothermal (Terebey et al. 1984), called the TSC model. This model is defined primarily by the effective sound speed (c_s,eff), the age since collapse began (t_col), and the initial angular velocity (Ω₀). In addition, we specify the outer radius for the envelope (r_out). The inner radius (r_in) is set by the dust evaporation temperature in an iterative way. The physical model is axisymmetric, so it is 2D, but the inclination angle adds a dimension, so the radiative transfer is done in 3D.

Because of the observational constraints on disk size (Yen et al. 2011, 2015a), the disk is not very significant in producing the observed SED, but we model it as a flared disk with parameters described by Yang et al. (2017). The disk is specified by its mass (M_disk), flaring power-law (β), and scale height at 100 au (h₁₀₀). These parameters do affect the near-infrared (3.6–8 μm) part of the SED.

We do not model the outflow itself, but only the effect of the cavity on the SED and radial profiles. The outflow lobes are characterized by their shape and the density structure inside them. The shape is defined in cylindrical coordinates, (z and ϖ), where z expresses the distance along the outflow axis and ϖ measures the perpendicular distance from that axis to the edge of the lobe. The equation of the surface is z = c₀ ϖ^1.5, where

$\begin{eqnarray}&&{c}_{0}=\displaystyle \frac{\cos {\theta }_{\mathrm{cav}}}{{r}_{\mathrm{fid}}^{0.5}{\sin }^{1.5}{\theta }_{\mathrm{cav}}},\end{eqnarray} \tag{ 1 }$

and θ_cav is the half-angle of the outflow lobe measured at r_fid = 10⁴ au, where r² = ϖ² + z². This model produces a curved shape that is qualitatively similar to the shape seen faintly in the Cycle 1 ALMA continuum data (Evans et al. 2015) and in Figure 6. The mass inside the outflow cavity is characterized by a small region of constant density (ρ_cav,0) that extends out to a radius (R_cav,0) beyond which the density declines as a power-law (ρ ∝ r⁻²). The axis of the outflow is assumed to be aligned with the rotation axis. Finally, the inclination angle of the rotation and outflow axes relative to the line of sight (θ_incl) has to be specified. A face-on rotation axis has θ_incl = 0.

In HYPERION, the luminosity source is characterized by a stellar radius (R_⋆) and temperature (T_⋆). For deeply embedded sources, these are nearly degenerate and only the stellar luminosity (L_⋆) is relevant. We kept T_⋆ fixed and varied R_⋆ to change the luminosity. The actual luminosity arises from accretion, so these parameters are not meant to be realistic. The value of T_⋆ affects the near-infrared portion of the SED somewhat.

The standard dust model is that of column 5 of Table 1 of Ossenkopf & Henning (1994) as extended to shorter wavelengths with anisotropic scattering by Young & Evans (2005). The standard Henyey–Greenstein model (Henyey & Greenstein 1941) for the angular distribution of scattering is used, unless otherwise noted.

Given the very large number of parameters in a 3D model, it is necessary to constrain as many as possible from existing information. The values of parameters that were fixed in the models are listed in Table 6 with a note that they were fixed. The distance was fixed at 164.5 pc (Watson 2020). The initial angular velocity, Ω₀, is set at 2 × 10⁻¹⁴ s⁻¹, based on the arguments in Yen et al. (2019), as noted in Section 1. The outer radius was estimated from Herschel maps by Launhardt et al. (2013) to be 4 × 10⁴ au, adjusted to the new distance. The outer radius, r_out, is set at 2 × 10⁴ au, about half the estimate from Herschel, but larger radii did not change the results. The mass enclosed in that value of r_out is 3.3 M_⊙. The opening angle (θ_cav) was set to 275, based on analysis of the CO outflow in the context of detailed source models (Shirley et al. 2011) and CO maps by Stutz et al. (2008). Most models assumed θ_incl = 87°, based on best fits to the SED by Stutz et al. (2008) using model grids and a fitter (Robitaille et al. 2006, 2007).

Table 6. Best-fit Model Parameters

Envelope parameters
tcol	40,000 yr
cs,eff	0.30 km s−1
Ω0	2.0 × 10−14 s−1	(fixed)
rout	2 × 104 au	(fixed)
Disk parameters
Mdisk	0.063	(fixed)
β	1.3
h100	4.0 au
Outflow cavity parameters
θcav	275	(fixed)
ρcav,0	4 × 10−20 g cm−3
Rcav,0	20 au
θincl	87°	(fixed)
Stellar parameters
T⋆	7000 K	(fixed)
R⋆	1.17 R⊙

Download table as:  ASCIITypeset image

The values of other parameters are the result of the optimization of the predicted SED and the radial profiles. As discussed at length in the appendix of Yang et al. (2017), the radial profiles and different regions of the SED best constrain different parameters. To compare the model predictions to the observations, a simulated SED was produced by calculating an emission image at a particular wavelength, using the ray-tracing part of HYPERION, and convolving that image with a top hat aperture of the specified size. Azimuthally averaged radial profiles were also computed from the model after simulating the beam used for the observations. Images of the model at specific wavelengths can also provide constraints. We will describe the steps used to optimize the model parameters in the order that they were taken, but the final optimization involves iteration between correlated parameters.

First, the luminosity has to be constrained. In the context of these models, a blackbody star is assumed. The stellar parameters of temperature (T_⋆) and radius (R_⋆) enter in combination ($L\propto {T}_{\star }^{4}{R}_{\star }^{2}$). To match the observed luminosity of 1.36 L_⊙ with θ_incl = 87°, the central luminosity must be 2.95 L_⊙, 2.2 times the observed luminosity. This result is a feature of models with outflow cavities, which cause the radiation to be channeled through the outflow lobes, and the effect is large in B335 because the outflow is nearly edge-on.

Second, the total mass of the envelope is best constrained by the longest wavelength photometry, given a model of dust opacities. In Shu-type models, this is set by the sound speed (c_s,eff) because the density depends on the sound speed: $\rho ={c}_{{\rm{s}},\mathrm{eff}}^{2}/(2\pi {{Gr}}^{2})$. The new, larger distance necessitated some changes in parameters used in recent modeling (Evans et al. 2015). The effective sound speed, including turbulence, c_s,eff had been taken to be 0.233 km s⁻¹ (Evans et al. 2015) in models of the source at 100 pc distance. At the larger distance, a larger sound speed is needed to match the flux into large beams at the longest wavelengths, with a value of c_s,eff =0.30 km s⁻¹ working well. The previous value was obtained by adding in quadrature the thermal sound speed for an initial T_K = 13K (Evans et al. 2005) and the microturbulent contribution inferred from line widths in single-dish beams (Evans et al. 2015). The obvious justification for increasing c_s,eff is that the Alfvén speed has not yet been included. If added in quadrature (Shu et al. 1987), the required Alfvén speed would be 0.19 km s⁻¹. The Alfvén speed is B(4π ρ)^−0.5. The large-scale magnetic field is estimated at 30.2 ± 17.7 μG (Kandori et al. 2020) but the relevant, pre-collapse density is unclear in the presence of a density gradient. The required Alfvén speed can be obtained if the relevant density is about 5 × 10⁴ cm⁻³, a very reasonable value. Zielinski et al. (2021) derived a field value of 142 ± 46 μG from observations at 214 μm and an assumed density of ρ = 5.41 × 10⁻¹⁸ g cm⁻³. This combination leads to nearly the same Alfvén speed of 0.17 km s⁻¹. The change in c_s,eff leads to changes in many other derived properties compared to previous models (Evans et al. 2015). These will be summarized in Section 5.2 after the other modeling parameters are discussed.

Third, the age of the source must be constrained. In TSC models, the age is the time elapsed since the initiation of infall (t_col). The age and sound speed determine the infall radius (${r}_{\inf }={c}_{{\rm{s}},\mathrm{eff}}{t}_{\mathrm{col}}$), which determines the shape of the density profile. In turn, the density profile is the main determinant of the normalized, azimuthally averaged, radial profile of emission at submillimeter wavelengths (i.e., the radial profile). As the source evolves, the initial ρ ∝ r⁻² density profile shifts to a flatter profile, leading to a flatter predicted radial intensity profile (Figure 25 of Yang et al. 2017). Because we have adjusted c_s,eff in the second step, we vary t_col to best approximate the observed radial intensity profiles. The observed radial profiles are annular averages, so we produce the model radial profiles in the same way, thus accounting for the effect of the outflow cavity. Models with different t_col are compared to the observations in Figure 11. The radial profile at 450 μm is best matched at t_col = 5 × 10⁴ yr, but the profile at 850 μm, which is less sensitive to the temperature distribution, favors t_col = 3 × 10⁴ yr. An age of 4 × 10⁴ yr predicts a radial profile at 850 μm that is slightly too flat and one at 450 μm that is slightly too steep, thus presenting about the best compromise. The larger distance and larger c_s,eff partially cancel out in their effects on the radial intensity profiles of the model, such that the best-fit t_col = 4 × 10⁴ yr is not very different from the age of 5 × 10⁴ yr favored by Evans et al. (2015). We adjusted M_disk to be 0.25 of the collapsed mass (M_⋆ + M_disk) because larger ratios tend to be unstable. For the age and infall rate fixed by the above considerations, M_⋆ = 0.25 M_⊙, so M_disk =0.063 M_⊙. For comparison, Evans et al. (2015) set a limit on the optically thin mass based on the Cycle 1 data; after updating to the new distance, the limit would be about 2 × 10⁻³ M_⊙. However, our modeling shows that much larger disk masses can be hidden by optically thick dust in the small disk of our models.

Zoom In Zoom Out Reset image size

Fourth, the mid-infrared part of the SED is most affected by the cavity properties, especially the density at the base of the cavity, with some effects from the disk properties. The best match was found for a radius for the region of constant density in the cavity (R_cav,0) of 20 au, and a density of the constant region (ρ_cav,0) of 4 × 10⁻²⁰ g cm⁻³. However, these parameters also affect the SED at shorter wavelengths, because the flux densities at those wavelengths are dominated by light scattered from the outflow cavity. The density that produces enough emission at 24 μm somewhat overproduces the emission in the small apertures in the IRAC bands, unless the disk and star parameters can compensate.

The near-infrared part of the SED (the IRAC bands) depends on many parameters, including the rotation rate, inclination angle, cavity properties, disk properties, and stellar properties. Deviations from the assumed simple, azimuthally symmetric structure can also have a large effect, so a good fit to the IRAC flux densities cannot necessarily be interpreted as actually representing the structure. With that caveat in mind, we describe the adjustments used to provide a decent match to the near-infrared observations. A lower T_⋆, with the larger R_⋆ needed to maintain the luminosity, produced more emission in the near-infrared. A better match was found for the (rather unrealistic) value of T_⋆ = 7000 K. We produced small grids of Ω₀, θ_incl, β and h₁₀₀ to optimize the near-infrared SED. The resulting best-fitting parameters are listed in Table 6. The common feature of these best parameters is that they minimize the emission in the IRAC bands. Some of the parameters are not very realistic. T_⋆ is quite large for a star of 0.25 M_⊙, and the disk parameters (β, h₁₀₀) produce a very flat disk. Even with these rather extreme parameters, the IRAC fluxes into small apertures are overproduced, especially at 3.6 μm. More typical stellar and disk parameters could be used if there is an extra source of extinction near the source, as will be discussed below.

The SED of the best model for the continuum is shown in Figure 12. The SED is well-matched except in the near-infrared, where the model overpredicts the emission into small beams and underpredicts the emission into large beams, and the ALMA observations, which we discuss below. At 3.6 μm, the model strongly overpredicts the emission into both large and small beams. Variations in the disk, cavity, and stellar properties were used to optimize the near-infrared, but the overprediction into small beams was a robust result. The 3D models produce a much better (though still imperfect) match to the short wavelength emission than did the 1D model with the same age. The outflow cavity has clearly allowed more radiation at mid-infrared and near-infrared wavelengths to reach the observer. The emission is essentially all scattered by the outflow cavities into the line of sight.

Zoom In Zoom Out Reset image size

Figure 13 shows a simulated image at 4.5 μm; two fan-shaped nebulosities, strongly peaked near the source and brighter along the cavity wall, are clearly seen. The shape agrees reasonably well with the image of the cavity walls indicated by the CS emission (Figure 6) and roughly with the image in Figure 9 of Stutz et al. (2008), but the model is more symmetric between the two lobes than are the infrared observations.

Zoom In Zoom Out Reset image size

The density structure of the best-fitting model is shown in Figure 14, which also defines the polar angle of the model, θ. The model has the outflow cavities along the z-axis by convention, so this image is rotated by 90° from the one in Figure 13. The radial distributions of dust density along five distinct polar angles are shown in Figure 15. The sharp increases show when the line of sight passes from the cavity into the envelope.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One glaring failure of the continuum model is the weakness of the emission predicted for the ALMA emission into a 2'' aperture. The best model predicts a flux that is a factor of 3 less than the Cycle 1 observations. Even the Cycle 1 observations may have occurred after the beginning of the rise in luminosity, and the continuum flux increases by a factor of 2 by the time of the 2018 observations. However, models with higher luminosity also fail to reproduce the ALMA data at any of the three epochs. A wide range of parameters for the disk, Ω₀, etc. were tried, but the predicted emission into a 2'' aperture was rather insensitive to these variations and was never close to the observed emission.

To compare to the ALMA data in more detail, a model image was made at the mean wavelength of the ALMA data and converted into UV-plane data using the CASA task simobserve for the Cycle 1 data. Then simanalyze was used to create an ALMA image that was fitted with an elliptical Gaussian within a 2'' circle. The emission in the model has a lower peak intensity and is much more elongated than the observed emission, being defined by the outflow cavity walls. So, neither the shape nor the flux are reproduced well. A similar issue was seen in BHR71 by Yang et al. (2020), who suggested the presence of an additional structure between the envelope and the Keplerian disk. The ALMA continuum data strongly suggest that there is a structure on the scales of ∼30 au in addition to the envelope and disk in our models. Such a structure could also provide additional extinction at short wavelengths, as suggested above in the discussion of the discrepancies in the near-infrared. The most likely candidate is something like a pseudo-disk, but a thermal radio jet could also contribute. Purser et al. (2016) found spectral indices as high as 1.6 in their study of radio jets from 5 to 23 GHz. If such a large index persisted to 345 GHz, the largest flux found by Galván-Madrid et al. (2004) could account for all the continuum that we observe. For a more typical index of 0.3, the contribution would be negligible.

In summary, 3D models with outflow cavities can reproduce observations of both the radial profile and most of the SED with TSC models of reasonable ages (t_col = 3–5 × 10⁴ yr). The outflow cavities are the critical element in resolving issues found in 1D models.

With a model for the quiescent source, we can relate the W2 photometric variation to variations in the luminosity. We ran a series of models, changing only the stellar radius (a stand-in for central luminosity), each predicting the W2 photometry. Then the relation between the model luminosity and the W2 flux density from the model was fitted with a quadratic function between about 3 L_⊙ and 25 L_⊙. The best-fitting model was

$\begin{eqnarray}&&L/{L}_{\odot }=1.837+0.845{S}_{\nu }{({\rm{W}}2)+0.001{S}_{\nu }({\rm{W}}2)}^{2},\end{eqnarray} \tag{ 2 }$

with S_ν (W2) the WISE photometry in mJy. The actual photometry as a function of time, multiplied by a factor of 0.67 to account for a slight underprediction by the models, then allows a model for the central luminosity as a function of time (Figure 16). This model predicts a luminosity increase of a factor of 5–7 over the quiescent state.

Zoom In Zoom Out Reset image size

In principle, the W1 photometry might be preferred because the W2 bandpass includes likely strong line emission from H₂ and CO (Yoon et al. 2022), perhaps balanced by absorption by CO, CO₂, and other ices. However, our model for the W2 emission matches the observations better than it does for the W1 bandpass. Time-sequence far-infrared photometry is what is really needed. There is one useful data point at 214 μm for MJD 58044, near the peak luminosity (Zielinski et al. 2021). Our luminosity model predicts L ≈ 20 L_⊙ at that date, but the 214 μm flux density is better matched in a model with L ≈ 6 L_⊙. There will be a time delay before the luminosity burst covers the full region that contributes to the far-infrared emission into relatively large beams (Francis et al. 2022), which would imply a luminosity larger than 6 L_⊙ at the time of the 214 μm observations. Consequently, the actual increase is rather uncertain, but likely a factor between 5 and 7.

To see how 3D models with outflow cavities and the small amount of rotation seen in B335 would fare, we used LIME (Brinch & Hogerheijde 2010) to model the line emission from the models computed for the dust emission with HYPERION. Because our dust model underestimates the ALMA observations, we also used a custom code (LIME-AID) to add a continuum source at the center of the model equal to the observed flux density (Yang et al. 2020). We used 5 × 10⁴ intensity grid points (pIntensity) for running models, but doubled that to 1 × 10⁵ for the final models. Using twice the number produced slightly smoother spectra, but the differences were very minor.

The spectral line models used the density and velocity fields from the TSC models, supplemented by the cavities added to model the continuum emission. In the inner regions, the gas temperature was assumed to be equal to the dust temperature from the HYPERION models, including the cavities, so the temperature is higher around the cavities. For radii larger than 2600 au, the gas temperature was increased, based on the models including gas heating by photoelectric heating by the interstellar radiation field (Figure 8 of Evans et al. 2005). Increased temperature on the outside of isolated globules is supported by the models of Herschel data (Launhardt et al. 2013). Because the adaptive gridding used by LIME was unduly distorted by the sharp discontinuity of the outflow cavity, we did not use the particle volume density inside the cavity of the continuum models, but instead set the abundance of the species being modeled to zero inside the cavity. The higher temperatures near the cavity were captured from the continuum models, but there were no changes in the gas velocity, so the outflows were not modeled.

The new variables introduced for the lines are those that characterize the microturbulence and the abundance profile. The microturbulence was taken as b = 0.12 km s⁻¹, as in previous models, and assumed to be independent of radius. The rotation was already in the TSC models, but for the line simulation, the sense of rotation matters as well. Based on a map of H¹³CO⁺ J = 1 → 0 emission, Kurono et al. (2013) found a velocity gradient from south to north, indicating that the north side is rotating away from us on large scales (r ∼ 2 × 10⁴ au). The opposite sense of rotation was inferred from the CO J = 2 → 1 maps by Stutz et al. (2008). Because the rotation in the models is clockwise as viewed from the "north pole" of the model, the left (right) side of the model is coming toward us and identified with south (north) for comparison to the data, depending on whether we follow Kurono et al. (2013) or Stutz et al. (2008). Given this uncertainty, we focus on simulating spectra directly toward the continuum source.

Chemical models were used to calculate the abundances of species in the envelope at various ages. The chemical model is described in detail in Zhao et al. (2018), so we provide here only a brief description. The network includes 21 neutral species, 31 ionic species, and 500 reactions. Freeze-out and desorption processes are included, but molecular reactions on the grain surface are not. Sulfur, and thus CS, is not included in the model. We also include all possible charge transfer reactions involving grains. The cosmic-ray ionization rate is set to 10⁻¹⁶ s⁻¹. The abundances were calculated for the temperature and density structure of the envelope at a given time. Rather than following the chemistry through the evolution, the abundances were precomputed for a range of temperatures and densities at a given chemical time to form a lookup table. The abundances are produced for a range of radii and polar angles. Because the abundances were nearly identical as a function of polar angle, we used the value for 45° for all angles and set the abundance to zero in the outflow cavity. Calculations were available for ages of 1, 3, 4, and 5 times 10⁴ yr. The abundances used in the models are shown in Figure 17. These were interpolated onto the grid used by LIME. These models were computed before the luminosity variation was discovered, so they assumed the pre-outburst luminosity.

Zoom In Zoom Out Reset image size

We focus first on the HCO⁺ emission, using the ALMA Cycle 1 data (Evans et al. 2015) and the ALMA Cycle 3 2018 data. To simulate the continuum source, we use the LIME-AID code (Yang et al. 2020) and insert a continuum source at the center matching the continuum before deconvolution, as noted above. The position angle was rotated by 90° because the model has the outflow in the north–south direction. Then we convolve the result from LIME-AID to the clean beam used for the data and extract a spectrum in the same way that we used for the data. We focus on ages of 1 × 10⁴ to 5 × 10⁴ yr, consistent with the results from the continuum modeling and past experience with the lines.

While many things affect the model line profiles, the clearest effect comes from the age of the system. Figure 18 shows line profiles for HCO⁺ and HCN J = 4 → 3 for ages from t_col = 1 × 10⁴ yr to 5 × 10⁴ yr. As the collapse proceeds, the infall velocities increase, moving the blue and red peaks apart, the lines broaden, and the self-absorption in the front part of the cloud broadens, eroding the lower-velocity emission in the red peak.

Zoom In Zoom Out Reset image size

Because of the luminosity variation, we produced models of the line emission for a range of luminosities at fixed age (t_col = 3 × 10⁴ yr). These show that the shape of the HCO⁺ line is little affected by the source luminosity, but the strength increases smoothly with luminosity (Figure 19). Because the chemical models do not include the luminosity variation, these line emission models are not self-consistent.

Zoom In Zoom Out Reset image size

To quantitatively compare models to observations, the figure of merit is the absolute value of the difference between observations and models, averaged over ±3 km s⁻¹about the v_LSR:

$\begin{eqnarray}&&\langle | \mathrm{Res}| \rangle =\langle | \mathrm{Obs}({\rm{v}})-\mathrm{Model}({\rm{v}})| \rangle .\end{eqnarray} \tag{ 3 }$

The chemical models have only the single variable of age. Ages of 1–5 × 10⁴ yr were tried for both HCO⁺ and HCN. The age of t_col = 3 × 10⁴ yr was clearly best for HCO⁺, but models with t_col = 4 × 10⁴ yr were not much worse. HCN was about equally (poorly) fitted by 3 × 10⁴ yr and 4 × 10⁴ yr. The resulting line profiles are shown in Figure 20. The model can reproduce the HCO⁺ Cycle 1 data reasonably. The centroid of the absorption dip is shifted to the red, matching the observations, but has more absorption near the systemic velocity than the observations. These failings apply to HCN as well, but in addition, the modeled HCN lines lack the broad wings seen in the observations. We are not modeling the outflow, a likely source of the broad wings in the observed HCN line.

Zoom In Zoom Out Reset image size

For the 2018 data, we clearly need to use a larger source luminosity. Our model based on the W2 photometry indicates L ≈ 20 L_⊙ when the 2018 data were taken, but this is uncertain as discussed above. We ran a series of models of increasing luminosity (controlled by increasing the stellar radius), while keeping the other variables fixed. The resulting HCO⁺ spectra in Figure 19 show that the line increases in strength as the luminosity increases. The model with L = 15 L_⊙ matches the 2018 HCO⁺ spectrum best. This luminosity worked best if t_col = 4 × 10⁴ yr as well.

The 3D model that fits the dust continuum emission best has the properties in Table 6. Both the SED and the radial profiles are reproduced reasonably well by this physical model with t_col = 4 × 10⁴ yr, though ages of 3–5 × 10⁴ yr are nearly as good. This result is a major improvement over 1D models, which could not match the whole SED and required unrealistically young ages to match the radial profiles. The same basic model, with t_col = 3 × 10⁴ yr, supplied with abundances of HCO⁺ from chemical models, gets the basics of the HCO⁺ observations correct—in particular the separation of the blue and red peaks, which is the main indicator of age. Moving to a 3D model mostly resolved the discrepancy between models of the continuum and models of the lines. The models are, however, not as effective in matching the Cycle 3 data nor the HCN line profile. Overall, models with ages t_col = 3–4 × 10⁴ yr provide the best fits to both continuum and line data.

With the revised distance and best-fitting age for the source, we need to update the basic source properties from those given in Evans et al. (2015). The best-fitting ages (t_col = 4 × 10⁴ yr from continuum modeling or t_col = 3 × 10⁴ yr from line modeling) are generally consistent with the outflow age of 2 × 10⁴ yr. We address first the properties of the source before the luminosity increase. The new distance results in an observed luminosity of 1.36 L_⊙; for the adopted inclination angle of 87°, the actual central luminosity must be 2.95 L_⊙. The greater distance also required a higher effective sound speed, c_s,eff = 0.30 km s⁻¹, leading to a higher mass inside r_out = 2 × 10⁴ au of 3.37 M_⊙. Because the mass infall rate is proportional to c_s,eff ³, it increases substantially to 6.26 × 10⁻⁶ M_⊙ yr⁻¹. With an age of t_col = 4 × 10⁴ yr (from the continuum model), the total mass that has fallen in would be 0.26 M_⊙. For t_col = 3 × 10⁴ yr (from the line models), it would be 0.19 M_⊙. The latter is more consistent with the value of 0.17 M_⊙ derived by Estalella et al. (2019) (scaled to the new distance) from analysis of the first moment of various observations (the central blue spot analysis), though either is probably consistent, given uncertainties.

For the larger mass of 0.26 M_⊙, the luminosity predicted for the infall rate is L = 43.6f_acc L_⊙. If the radiative efficiency of accretion f_acc = 0.5, the predicted luminosity is 21.8 L_⊙, close to the peak inferred from modeling the W2 variation. In an episodic accretion picture, the pre-outburst luminosity resulted from an accumulation of mass in an intermediate reservoir, such as a pseudo-disk or a Keplerian disk. The luminosity outburst would then represent a "catch-up" to match the rate at which envelope material is infalling.

While modeling in three dimensions has resolved some puzzles, newer ones are now apparent. The failure of models to reproduce the continuum emission from ALMA suggests the existence of additional components, such as a pseudo-disk or a radio jet with a large positive spectral index. Observations in the 1–10 mm wavelength region would resolve the latter possibility. The constraints based on photometry from 3 to 24 μm should be viewed skeptically in light of deep ice absorptions and strong line emission that are emerging from early observations with JWST (Yang et al. 2022). Contamination of the NEOWISE filters by likely CO rovibrational emission, H₂ emission lines, atomic lines, and ice absorptions render suspect the modeling of the luminosity variation in the one data set with time-domain observations. Chemical models with varying luminosities combined with the full array of ALMA observations may be able to reconstruct the luminosity history, and JWST spectra will dramatically improve our knowledge of the near-source physical and chemical environment. Future models will need to consider departures from axisymmetry, given the evidence for infalling streamers.