Searching for Cool Dust. II. Infrared Imaging of The OH/IR Supergiants, NML Cyg, VX Sgr, S Per, and the Normal Red Supergiants RS Per and T Per^∗

The targets were observed with SOFIA/FORCAST during Cycles 3 & 4 (OBS IDs: 03_0082, 04_0013; PI: R. M. Humphreys). FORCAST is a dual-channel mid-IR imager covering the 5–40 μm range. Each channel uses a 256 × 256 pixel blocked-impurity-band (BiB) array and provides a distortion-corrected 32 × 32 field of view with a scale of 0768 pix⁻¹. FORCAST achieves near-diffraction-limited imaging, with a PSF FWHM of ∼37 in the longest filters. We elected to image in single-beam mode to maximize throughput. The observations were obtained in the standard two-position chop-and-nod mode with the direction of the nod matching the direction of the chop (NMC). The data were reduced by the SOFIA Science Center using the FORCAST Redux pipelines version 1.0.3 (S Per), 1.0.5 (VX Sgr), 1.0.7 (NML Cyg), and 1.1.0 (RS Per, T Per). After correcting for bad pixels and droop effects, the pipeline removed sky and telescope background emission by first subtracting chopped image pairs and then subtracting nodded image pairs. The resulting positive images were aligned and merged. The details of the FORCAST pipeline are discussed in the Guest Investigator Handbook for FORCAST Data Products, Rev. A3.⁶ The FORCAST and MIRAC observations are summarized in Table 1.

Bright point sources cause cross-talk in the horizontal direction on the FORCAST array. To mitigate this effect, chop angles were selected so that the cross-talk pattern from one chop position did not overlap with the other chop position. Additionally, the FORCAST pipeline applies a correction that reduces the effect, although some of the pattern remains for some targets. The effect is strongest for the brightest IR targets, NML Cyg and VX Sgr; it is less so for S Per and is not present in the images of RS Per or T Per. However, one effect that may appear in some of the fainter targets (especially T Per) is possible coma introduced from the NMC chopping pattern.⁷ This effect may explain the asymmetries in the surface brightness profile of T Per, shown in Figure 10 and discussed in Section 3.4.

For each of the stars, observations of the asteroid 2 Pallas were used for PSF calibration in the same filters with the four-position slide in either the mirror position (for the short wavelength channel) or the open position (for the long-wavelength channel). The color temperature of the asteroid Pallas (∼160 K) is far less than the effective temperatures of the target stars (≳3200 K). This color difference is not ideal for a PSF calibrator, as the cooler source will peak at longer wavelengths, possibly resulting in a broader profile. However, Pallas was the only source observed in each SOFIA cycle under the same conditions and at each wavelength studied in this work. We analyzed another calibrator (α Aur) at 11.1 and 31.5 μm and measured a similar FWHM at each wavelength. We present the profiles of Pallas in the figures below for consistency, acknowledging the possibility that we have overestimated the size of the PSF at the longer wavelengths.

The FORCAST pipeline coadds the merged images. We use the standard deviation of the mean of fluxes extracted from the merged images (prior to coadding) as the 1σ uncertainty of the fluxes in the coadded images of each of our targets. This uncertainty is negligible compared to the 6% uncertainty that we adopt for the flux calibration, per the GI Handbook Section 4.1 (Herter et al. 2013). The bandpasses of the selected FORCAST filters are such that only small color corrections are required. Based on the F_ν ∝ ν² spectral shapes of our targets in the relevant ranges, we have applied color corrections of 1.004, 1.071, 1.004, 1.044, 1.025, and 1.025 to fluxes extracted from the F111, F197, F253, F315, F348, and F371 images, respectively. Aperture photometry was performed using the open-source Astropy (Astropy Collaboration et al. 2013) affiliated photutils⁸ package. Apertures span between 15'' and 20'', chosen to encompass the extended emission around each object. FORCAST photometry is reported in Table 2 and included in the SEDs in Section 3. Photometric error is reported as measured uncertainty in the sky background apertures.

NML Cyg, S Per, and T Per were observed with the mid-IR adaptive optics system on the MMT using the Mid-infrared Array Camera and Bracewell Infrared Nulling Cryostat (MIRAC3/4/MIRAC-BLINC; Hoffmann et al. 1998; Hinz et al. 2000; Skemer et al. 2008). The observations of NML Cyg are described in Schuster et al. (2009), and discussed here in Section 3.5. S Per was observed on UT 2006 November 05 and T Per on UT 2009 October 02 at 8.9 and 9.8 μm. MIRAC achieved Strehl-ratios close to 0.95, providing diffraction-limited imaging and stable PSFs (e.g., Biller et al. 2005). MIRAC4 employed a Si:As array with 256 × 256 pixels, and observations were made with a standard chop-and nod sequence to remove IR background emission. Cross-talk in the array electronics introduced faint artifacts in the horizontal and vertical directions, which is not completely removed by chop-and-nod subtraction. As described in Paper I, the horizontal cross-talk is mitigated during the data reduction with a code from Skemer et al. (2008). For consistency with the FORCAST photometry, we perform aperture photometry on the MIRAC images using photutils. We report the results in Table 2, include the photometry in the SEDs below, and as input to DUSTY.

To populate the mid-IR SEDs, we include IRAS photometry (and AKARI photometry when available) from point-source catalogs in the literature for RS Per and VX Sgr (Smith et al. 2004), S Per and T Per (Abrahamyan et al. 2015), and NML Cyg (Schuster 2007). The Abrahamyan et al. (2015) catalog cross-correlates IRAS point sources with WISE, the latter of which presents some issues due to its large beam-size (up to 12'' at 22 μm; Wright et al. 2010). For stars embedded in nebulosity or crowded fields, the WISE photometry can be systematically too bright.

Additionally, optical photometry is compiled from the Extended Hipparchos Compilation catalog (XHIP; Anderson & Francis 2012), the SKY2000 Master Catalog (Myers et al. 2015), or the AAVSO Photometric All Sky Survey (APASS; Henden 2016, see SEDs in Section 3). These optical data, as well as the published photometry from IRAS and AKARI, are dereddened using the extinction law from O'Donnell (1994). The values for interstellar extinction A_V chosen for each source are listed in Table 3 and in the SED captions below.

Table 3.  DUSTY Model Parameters and Mass-loss Rates

Model	AVa	Teff	Tdustb	τV	r1	$\dot{M}$ c
		K	K		au	${\dot{M}}_{\odot }\,{\mathrm{yr}}^{-1}$
Inputs				Outputs
VX Sgr
r−2.0	2.0	3200	1000	6.7	86	4.5 × 10−5
r−1.6				3.7	76	2.2 × 10−5
S Per
r−2.0	3.1	3500	1000	1.4	45	2.6 × 10−5
r−1.6				1.2	43	2.4 × 10−5
RS Per
r−2.0	1.7	3600	1000	0.3	50	3.5 × 10−5
r−1.6				0.3	50	3.9 × 10−5
T Per
r−2.0	2.1	3700	1000	0.1	30	7.5 × 10−6
r−1.6				0.1	31	8.1 × 10−6
NML Cyg
r−2.0	4.0	3300	1000	41	133	4.8 × 10−4
r−1.8				37	128	4.2 × 10−4

Notes.

^aDUSTY output models are fit to extinction-corrected SEDs with these values of A_V. ^bDust temperature at condensation radius, r₁. ^c $\dot{M}$ is computed as an average mass-loss rate over the lifetime of the shell. The outflow velocity is assumed to be 25 km s⁻¹ unless noted in the sections for the individual stars.

Download table as:  ASCIITypeset image

We also compile spectra from ISO-SWS (de Graauw et al. 1996) for all targets except T Per. S Per and RS Per spectra are from the Japanese guaranteed observing time program REDSTAR1 (PI T. Tsuji; Aoki et al. 1998), and NML Cyg and VX Sgr were observed with the AGBSTARS program (Justtanont et al. 1996; Speck et al. 2000). The color and extinction-corrected spectra are displayed in the SEDs below and are provided as near- to mid-IR photometric input to DUSTY.

We also include in our analysis the publicly available 70 and 160 μm observations made with Herschel/PACS. VX Sgr, NML Cyg, and S Per were observed as part of the Herschel key program Mass-loss of Evolved StarS (MESS; Groenewegen et al. 2011). The Herschel Interactive Processing Environment (HIPE; Ott et al. 2010)⁹ was used to download the images, but photometry was performed using photutils for consistency with the SOFIA images. Apertures span between 45 and 70'' to encompass extended emission around each object. As the PACS pixels are large on sky, we did not have enough pixels in traditional sky annuli to model the background. Instead, we first mask the star and its nebulosity and then model the background across the field of view as a two-dimensional polynomial. For each of the PACS fields, these background models were fairly flat but had high rms variation. As summarized in Table 2, this uncertainty was as high as ∼40% for VX Sgr and S Per.

The width of the PACS bandpasses requires color corrections to be applied to the 70 and 160 μm photometry from the images. In Paper I, the necessary corrections were estimated by convolving the "blue" (70 μm) filter response functions to the sources' ISO LWS spectra. However, lacking spectra for all of the sources in this work, we instead fit the mid- to far-IR photometry from SOFIA and IRAS with a power-law of the form F_ν = ν^β to represent the targets' SEDs at the PACS wavelengths. The results are modest corrections of 1.003 and 1.04 for the two bandpasses. PACS photometry is reported in Table 2 and included in the SEDs in Section 3. Photometric error is reported as measured uncertainty in the sky background models.

To estimate the mass-loss rates, mass-loss histories, and dust density distributions, we used the DUSTY radiative-transfer code (Ivezic et al. 1997) to model the observed SEDs and azimuthal average intensity profiles at each of the MIRAC and FORCAST wavelengths, in a manner similar to that used in Paper I. DUSTY solves the one-dimensional (1D) radiative-transfer equation for a spherically symmetric dust distribution around a central source. We provide as input the chosen optical properties, chemistry, size distribution of the dust grains, and a dust temperature, which fixes the inner boundary of the surrounding dust shell (the dust condensation radius, r₁). We generate a grid of models for each star with fixed stellar effective temperatures based on the spectral type of each target, fixed shell extent (1000 × r₁), and fixed dust condensation temperature (1000 K). Our grid consists of varying optical depths of the circumstellar material (0.01 < τ_V < 50) and different dust density distribution functions, described below. For a given set of inputs, DUSTY outputs a model SED and radial profiles of the dust shell at requested wavelengths.

As noted in Paper I, the spherical symmetry assumed by DUSTY fails to model the azimuthal complexities observed in the asymmetric outflows of massive stars such as VY CMa (Smith et al. 2001; Humphreys et al. 2005, 2007; Shenoy et al. 2013) and IRC +10420 (Humphreys et al. 1997; Tiffany et al. 2010; Shenoy et al. 2015). However, DUSTY allows for a consistent analysis of the dust, SEDs, and intensity profiles of the targets in this work and those in Paper I.

At a given wavelength, an output optical depth τ_λ from the model, and thus its grain opacity κ_λ, specifies the dust mass density $\rho (r)$ throughout the shell. If we assume a constant expansion rate v_exp of the outflowing dust shell, following the arguments set forth in Paper I, we can estimate the mass loss rate as:

$$\begin{eqnarray*}&&\dot{M}(t)={g}_{d}\,4\pi {r}^{2}\,\rho (r){v}_{\exp }\end{eqnarray*}$$

where radius r is a probe on timescale t as r = v_exp t, and g_d is the gas-to-dust ratio. For consistency with Paper I, we assume g_d = 100:1 (Knapp et al. 1993); however, this can be as high as 200:1 for supergiants (Decin et al. 2006; Mauron & Josselin 2011). For the dust optical properties, we use the "cool" circumstellar silicates from Ossenkopf et al. (1992), and assume the grain radii follow a Mathis, Rumpl, Nordsieck (MRN) size distribution $n(a)\propto {a}^{-3.5}\,{da}$ (Mathis et al. 1977) with a_min = 0.005 μm and a_max = 0.25 μm.

In general, the mass density distribution of the outflows can be modeled with DUSTY as a power-law $\rho (r)\propto {r}^{-q}$. An index of q = 2 is the case of constant mass-loss rate and constant expansion velocity for the shell, while q < 2 indicates a gradual decline in the mass-loss rate over the dynamical age of the expanding shell. A steeper power-law index q > 2 represents a mass distribution with more recent high mass loss, and less significant mass loss in the past. For each of our targets, the fundamental research question is how well the stars' SEDs and radial profiles in the mid-IR can be modeled with a constant mass-loss rate scenario in DUSTY.

For each of the targets in our sample, we perform three Monte Carlo experiments. In the first, we force DUSTY to use the constant mass-loss rate distribution $\rho (r)\propto {r}^{-2}$, and by varying the optical depth of the CS material, recover the best-fitting, constant mass-loss SED in the near- to mid-infrared. For the second set of simulations, we allow the power-law index of the mass distribution to vary between 1 and 3 with a step size of 0.2, deliberately excluding the r⁻² case, while also allowing the optical depth to vary. For both set of DUSTY models, we evaluate the best fit based on a reduced χ² measurement of the extinction-corrected SED and the DUSTY output spectrum. We then compare the DUSTY-predicted intensity profiles to the observed radial profiles in the SOFIA wavelengths (and MIRAC, when available). The image and profile models output from DUSTY do not account for the optics of the telescopes, so the intensity profiles are convolved with an azimuthal average of the PSF and are displayed in the figures below.

For the third and final set of DUSTY models, we select the best-fitting r⁻² model, and re-run DUSTY with those same parameters, this time enhancing the dusty density profile by a factor of 10 at 50 condensation radii (50 × r₁). An example model is shown in Figure 1. These "enhanced," piecewise-defined models explore the possibility of an extreme mass-loss event in a star's past, similar in some respects to the models from the second experiment described above where the density distribution can be shallower than r². The latter DUSTY models imply a smoothly changing mass-loss rate over the lifetime of the star, whereas the "enhanced" models simulate a single eruptive event in the mass-loss history. Similar piecewise defined density profiles were used in Paper I to model the SED of IRC +10420. Though a factor of 10 enhancement in mass-loss rate is likely an extreme case, we apply this model to explore how well the scenario of a constant mass-loss rate with a single one-off eruptive event reproduces the observed IR SED of our target stars.

Zoom In Zoom Out Reset image size

Finally, we can estimate an average mass-loss rate for the non-constant mass-loss models ($q\ne 2$) by integrating the density distribution $\rho (r)$ and multiplying by the gas-to-dust mass ratio (100:1) to compute the total mass of the shell M. We assume an average expansion velocity to estimate the dynamical age of the shell Δt = r₂/v_exp where r₂ is the outer radius of the shell predicted by a given model. The expansion velocity, v_exp, is assumed to be 25 km s⁻¹ unless specified in the sections for individual stars below. The average mass-loss rate is then $\langle \dot{M}\rangle =M/{\rm{\Delta }}t$. The specific parameters for the DUSTY models for each target in our program, as well as the output DUSTY models and computed average mass-loss rates, are summarized in Table 3, and the best-fitting SEDs to the observed photometry are shown in the figures below.

Note that in Table 3, the first row for each target star represents the best-fitting DUSTY simulations forced to evaluate the models in the constant mass-losing, r⁻² dust profile case. The second row represents the best-fitting SED with non-constant mass-loss. The columns on the left reflect the input values, and the right-hand columns are the recovered output parameters from the best-fitting models for each target and each simulation set (constant versus non-constant mass-loss rates). We do not include the enhanced, piecewise-defined models here, as the parameters were fixed to the r⁻² model for each star. Throughout the text, we will refer to the three different models as constant (r⁻²), non-constant (${r}^{-q},q\ne 2$), and enhanced (r⁻², e) mass-loss rates.

VX Sgr has a marginally resolved, nearly symmetric extended circumstellar envelope in its HST visual images (Schuster et al. 2006). Additionally, Vlemmings et al. (2005) has identified a dipole magnetic field in its ejecta mapped by its H₂O masers, which may be a clue to its mass loss mechanism. VX Sgr is also a semi-regular variable that behaves like a fundamental mode pulsator (i.e., a Mira variable), which is rare for such a luminous star. It has been observed to vary by several magnitudes with corresponding changes in its apparent spectral type from M4 to M10.

During one of its Mira-like episodes, Humphreys & Lockwood (1972) noted a decline of ∼0.5 mag out to 10 μm over a few months. Therefore, we have chosen optical photometry to align in light-curve phase with the 2MASS and IRAC photometry compiled in Smith et al. (2004). To constrain the 2–10 μm regime of the SED, we also include the photometric average over the 3.6-year cycle of the COBE DIRBE project (Price et al. 2010).

The SED is shown in Figure 2, with observed data plotted as open symbols and extinction-corrected photometry in solid. The constant mass-loss DUSTY model is overplotted with a dashed blue line, and the best-fitting power-law model with index q = 1.6 is shown by the dashed–dotted green line. Note that both models fit the 10 μm silicate feature in the ISO spectrum, but the q = 1.6 model better simulates the cool thermal dust emission out to 100 μm. However, the observed IR flux from ∼30 to 70 μm falls in between the two models. Displayed in dotted red is the "enhanced" DUSTY model, which is the r^−2.0 constant mass-loss model with a factor of 10 enhancement in dust density at 50 × r₁. This model appears to over-estimate the thermal dust emission, implying too much dust is produced to match the observations of VX Sgr.

Zoom In Zoom Out Reset image size

The derived mass-loss rates for the models are summarized in Table 3. The outflow velocity v_exp adopted for this calculation is 24.3 km s⁻¹ from the AGB/supergiant CO-line survey by De Beck et al. (2010). The mass-loss rates from DUSTY (2–5 × 10⁻⁵ M_⊙ yr⁻¹) are somewhat lower than the measurements from CO-line profiles (6.1 × 10⁻⁵ M_⊙ yr⁻¹; De Beck et al. 2010). However, the most obvious explanation for this is due to the assumed gas-to-dust ratio. Where we have assumed 100:1 for consistency with Paper I, De Beck et al. (2010) allow the gas-to-dust ratio to vary when fitting the observed outflow velocities (using GASTRoNOoM; Decin et al. 2006). Mass-loss rates scale linearly with the gas-to-dust ratio, so if we had applied a ratio of 200:1, perhaps more appropriate for RSGs (Decin et al. 2006; Mauron & Josselin 2011), our estimated mass-loss rate would be more consistent with the derived measurement from CO-line profiles.

In Figure 3, we compare the observed radial profiles to the PSF calibrator (2 Pallas) and the DUSTY output image models. The ejecta around VX Sgr is only marginally resolved above the PSF at 19.7 μm, but the envelope is more easily distinguished from the PSF at longer wavelengths. The DUSTY model profiles are convolved with the PSF in each band, and we note that the constant mass-loss rate model aligns more closely with the observed surface brightness profiles of VX Sgr. Both the shallower r^−1.6 and enhanced models over-estimate the amount of thermal dust emission observed in the FORCAST images.

Zoom In Zoom Out Reset image size

The azimuthal profiles combined with the SED modeling suggest that the mass-loss rate of VX Sgr is fairly constant with perhaps a period of elevated mass-loss in the past, as illustrated by the infrared excess emission in the 20–160 μm photometry. The model intensity profiles, though, all predict a higher surface brightness for the extended emission than what was actually observed at SOFIA wavelengths (Figure 3). One possible explanation for not observing this emission is that VX Sgr has experienced a sudden decline in mass-loss rate in very recent times. Exploring this possibility is beyond the scope of this paper, but we plan to make high-resolution 5–12 μm observations using LMIRCam and NOMIC on the LBT (Skrutskie et al. 2010; Hoffmann et al. 2014). At FWHM spatial resolutions of 012 and 029 at 5 and 12 μm, respectively, we can explore the dust shell at ∼200 au scales and combine these observations with our SOFIA data and DUSTY modeling.

We also note that although the envelope is spherically symmetric in HST optical images (Schuster et al. 2006), the H₂O masers around VX Sgr appear to align with the equatorial plane of the star's dipole magnetic field. Vlemmings et al. (2005) suggest that this alignment could create an overdensity in the circumstellar material in this plane, as modeled by Matt et al. (2000). While we do not see evidence for asymmetry in the FORCAST images, we note again that DUSTY assumes spherical symmetry in its models. As discussed further with S Per and NML Cyg below, it is likely that DUSTY may fail to accurately model stars with known asymmetric outflows and profiles.

S Per (Sp. Type M3-4e Ia) is an OH/IR source and a member of the Per OB1 association (Humphreys 1978), with a distance of 2.3 ± 0.1 kpc as determined by VLBI H₂O maser astrometry (Asaki et al. 2010). Schuster et al. (2006) present HST images showing that the star is embedded in an elongated circumstellar envelope with a position angle of ∼20° E of N with a FWHM of ∼01 (240 au). Schuster et al. (2006) speculated that the shape could be due to bipolarity in the star's ejecta or a flattened circumstellar halo, and they note this elongated structure is also seen in OH and H₂O maser observations on the same scale and with similar orientation (Richards et al. 1999; Vlemmings et al. 2001). Fitting elliptical Gaussians to S Per's MIRAC4 images yields a mean position angle of 19° ± 2° E of N, matching the orientation seen in the HST images.

The observed SED is shown in Figure 4 along with the three DUSTY models. Both the shallower q = 1.6 and enhanced DUSTY models accurately reconstruct the near- to mid-IR flux, while the constant mass-loss rate model underestimates the thermal dust emission. We note here one possible complication in our analysis. DUSTY simulations assume spherical symmetry in CS material, which could lead to underestimating the density, and thus optical depth, of the ejecta relative to the observed compact envelope seen in the Schuster et al. (2006) WFPC2 images of S Per. Our models for this star, then, may not best represent the stellar outflows and dusty envelope.

Zoom In Zoom Out Reset image size

The observed azimuthal average radial profiles from SOFIA/FORCAST for S Per are presented in Figure 5. S Per has resolvable extended emission above the PSF; however, the q = 2 and q = 1.6 profiles, once convolved with the large PSF beam of FORCAST, are virtually indistinguishable. The enhanced DUSTY model, though, predicts too much emission close to the central star.

Zoom In Zoom Out Reset image size

In Figure 6, we illustrate the surface brightness profiles at higher spatial resolution with MIRAC. Here, the observed surface brightness profile is clearly resolved above the PSF at the shorter wavelengths; however, the DUSTY models underestimate the shape of the stellar envelope. Adding a period of enhanced mass loss, the DUSTY model in dotted red, produces too much emission at the shorter wavelengths. Unfortunately, then, the radial profile models do not provide any conclusive evidence that the mass-loss history of S Per is constant versus non-constant. Note that the deviations from a smooth profile in the 31.5 μm and the two MIRAC figures are due to the asymmetry in the outflows. From the SED, though, we glean that S Per may have had a higher mass-loss rate in the past, but we acknowledge that DUSTY is not ideal for simulating stellar ejecta of stars with known bipolar/asymmetric envelopes.

Zoom In Zoom Out Reset image size

As reported in Table 3, the two DUSTY models predict mass-loss rates between ∼2 and 3 × 10⁻⁵ M_⊙ yr⁻¹. Richards et al. (1999) summarizes results from previous literature to show a range of published mass-loss rates from as low as 7 × 10⁻⁶ M_⊙ yr⁻¹ (OH 1612 MHz; Jura & Kleinmann 1990) to as high as 2 × 10⁻⁴ M_⊙ yr⁻¹ (CO-line profiles; Knapp & Morris 1985). With such a large range of published values, each measuring mass-loss rates with a different observational technique, we can only conclude that we have derived a rate within published bounds.

Fok et al. (2012) also performed DUSTY modeling on a number of Galactic RSGs, including S Per, T Per, and RS Per. However, they used a different mode, the "dusty AGB" radiatively driven wind mode, and a higher gas-to-dust ratio of 200:1. The radiatively driven wind mode in DUSTY is provided for modeling AGB star envelopes and is not necessarily appropriate for RSGs (e.g., Heras & Hony 2005). Groenewegen (2012) analyzed the systematic difference in mass-loss rates computed using DUSTY in this mode as compared to the default where the user supplies the density distribution as a power-law function. Groenewegen (2012) found that the mass-loss rates computed with the radiatively driven wind mode differ significantly from those obtained from the default mode in which the equation of radiative transfer alone is solved when applied in the context of RSGs. Thus, our results, with fixed single-component power-law distributions, are not directly comparable to the Fok et al. (2012) models. Nonetheless, Fok et al. (2012) yield a best-fitting model with $\dot{M}=1\times {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. We note that Gehrz & Woolf (1971) derived a mass-loss rate for S Per of 2.7 × 10⁻⁵ M_⊙ yr⁻¹ using an independent analysis of the 3.6–11.4 μm SED. These results are consistent with our measurements.

Departures from circular symmetry have been reported for both RS Per and T Per based on H-band interferometric imaging with CHARA by Baron et al. (2014) at an angular scale of 1.3 mas, which the authors attribute to surface asymmetries or spots. We do not see much evidence for asymmetry in the FORCAST images, and the azimuthal profiles for RS Per do not show significant excess emission above the PSF, though the angular scales for FORCAST are much larger than Baron et al. (2014) observed with the CHARA array. RS Per has the 10 μm silicate emission feature in its SED but is not a known maser source. It is likely a normal RSG that may just be entering a more active phase with enhanced mass loss, perhaps driven by surface activity like that seen in the OH/IR supergiants and VY CMa. T Per, another member of the Perseus OB1 association, similarly shows no evidence for SiO maser emission (Jiang et al. 1999). Both stars exhibit long-period variability of ∼4200 and 2500 days for RS Per and T Per, respectively (AAVSO; Kiss et al. 2006).

The SED of RS Per is shown in Figure 7. Both models fit the 10 μm silicate feature, as well as the ISO spectrum, from 2–11 μm. However, the constant mass-loss, q = 2 profile underestimates the flux for the mid- to far-IR at wavelengths larger than 20 μm, while the shallower dust density distribution q = 1.6 profile and the enhanced r^−2.0 profile better match the longer wavelength SED. However, after convolution with the FORCAST PSF, the radial profiles of both models in Figure 8 appear too similar to the constant mass-loss case beyond 2" to distinguish either as a better fit to the extended emission. For T Per, shown in Figure 9, both the constant mass-loss and enhanced models produce insufficient thermal dust emission at the longer wavelengths.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As summarized in Table 3, the mass-loss rates derived for RS Per and T Per are 4 × 10⁻⁵ and 8 × 10⁻⁶ M_⊙ yr⁻¹, respectively. Fok et al. (2012) estimates a mass-loss rate of 3.0 × 10⁻⁶ M_⊙ yr⁻¹ for RS Per and 5 × 10⁻⁷ M_⊙ yr⁻¹ for T Per using DUSTY "AGB mode." As discussed in Section 3.3, this radiatively driven wind mode is less appropriate for RSGs, and so our results are not directly comparable. Additionally, we note that the SED fits provided in their work do not extend longward of 30 μm, so we cannot qualitatively gauge which SED modeling mode (our power-law profiles versus their radiatively driven wind models) would fit best with our new FORCAST photometry through 40 μm. Particularly in the case of T Per, we note that this long-wavelength IR photometry is crucial in constraining the models.

While the r^−1.6 DUSTY model is clearly the better fit to the observed mid-IR SED, the radial profiles of RS Per in Figure 8 reveal that the SOFIA images lack the spatial resolution necessary to favor one model over the other as a better fit to the extended envelope emission. The radial profiles of T Per in Figures 10 and 11 reveal a much clearer extended profile above the PSF at 31.5 μm and around the 10 μm silicate feature with MIRAC. However, the two power-law models predict very similar output profiles. The enhanced model accurately recovers the extended emission around the 10 μm silicate feature, but the MIRAC images in the two other wavebands do not resolve emission extended above the PSF.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note a curious ripple in the 31.5 μm FORCAST profile in Figure 10. As mentioned in Section 2.1, some optical aberrations may be introduced to the images due to chopping patterns of the secondary mirror on SOFIA. One possible source of the profile shape could be coma effects that stretch out T Per along one axis. When generating azimuthally averaged surface brightness profiles, this asymmetry would cause the profile to deviate from a smooth power-law.

HST visual images of NML Cyg revealed a peculiar bean-shaped asymmetric nebula only ≈02 across and coincident with the distribution of its H₂O masers. Schuster et al. (2006) showed that its circumstellar envelope is shaped by photodissociation from the powerful nearby association Cyg OB2 inside the Cygnus X superbubble, which is relatively void of gas and dust. This configuration allows the UV radiation from the numerous luminous hot stars in Cyg OB2 to travel the ≈80 pc to NML Cyg unimpeded. Subsequent adaptive optics mid-IR imaging at 8.8, 9.8, and 11.7 μm with MIRAC3 on the MMT (Schuster et al. 2009) spatially resolve the physical structures near the star (∼240 au) responsible for its 10 μm silicate-absorption feature and an asymmetric excess at 03–05 from the star due to thermal emission from hot dust. This excess is also oriented toward the Cyg OB2 association and is attributed to the destruction of NML Cyg's dusty wind by the hot stars in Cyg OB2.

As illustrated in the SED in Figure 12, the 10 μm silicate feature is seen in absorption, rather than emission. Discussed in detail in Schuster et al. (2006, 2009), this absorption is due to NML Cyg's thick circumstellar envelope obscuring the central star. Indeed, the DUSTY models predict a large optical depth (τ_V > 40 for both models) to fit the silicate feature in absorption as well as the mid- to far-IR photometry.

Zoom In Zoom Out Reset image size

At the FORCAST wavelengths, 19.7–37.1 μm, NML Cyg appears as a point source, with no evidence of cold dust hidden or protected from the UV radiation in Cyg OB2. Additionally, no preferential extension toward Cyg OB2 is evident at the angular resolution of FORCAST. The radial profiles in Figure 13 do not show any obvious excess emission above the PSF. While we note that the shape of the observed azimuthal profile seems similar to the r⁻² model, the models seem to greatly over-predict the surface brightness flux for the large estimated optical depth.

Zoom In Zoom Out Reset image size

As noted in Schuster et al. (2009), the MIRAC images do indeed appear asymmetric along the NW–SE axis. For these images, we generate isophotes and separate our radial profiles into two axes. The major axis (NW–SE) is shown with solid points in Figure 14, and the minor axis (NE–SW) is shown with open circles. Here, we see that the constant mass-loss, q = 2, model does fit the observed major-axis brightness profile for the 9.8 μm image, though the q = 1.8 profile is similar in shape for all three MIRAC bands. We note that DUSTY cannot take into account the external radiation field from Cyg OB2, so the model profile shapes are not necessarily conclusive as to the mass-loss history of NML Cyg.

Zoom In Zoom Out Reset image size

NML Cyg has one of the highest mass-loss rates of any RSG/hypergiant—from 6.4 × 10⁻⁵ M_⊙ yr⁻¹ (Morris & Jura 1983) to as high as 1.6 × 10⁻⁴ M_⊙ yr⁻¹ (Lucas et al. 1992). We calculate even higher mass-loss rates from the DUSTY models, though, at 4–5 × 10⁻⁴ M_⊙ yr⁻¹ with average outflow velocity 23 km s⁻¹ (Schuster et al. 2009). It is likely, however, that the complex morphology of NML Cyg cannot be well-modeled with a single-component power-law dust mass distribution, similar to the difficulties in modeling S Per with its known asymmetric profile.