2.1. Infrared Data

The source IRAS 08589−4714 was observed with the Herschel space telescope³ (^{Pilbratt et al. 2010}) in the bands of the Photodetector Array Camera and Spectrometer (PACS, 70 and 160 μm; ^{Poglitsch et al. 2010}) and of the Spectral and Photometric Imaging Receiver (SPIRE, 250, 350, and 500 μm; ^{Griffin et al. 2010}), with nominal FWHM of about, 5'', 12'', 18'',25'', and 36'' in the five bands, respectively. The data were gathered as part of the Herschel Infrared Galactic Plane Survey (Hi-GAL, ^{Molinari et al. 2010}), which mapped the Galactic plane in mosaics of ≈ 2.3° x 2.3°, taken in the five bands simultaneously. IRAS 08589−4714 is located in the Hi-GAL field centered at [l, b] = [268.3°, −01.25°]. These mosaics were obtained in the parallel mode, at a scan velocity of 60''/s on 2012 November 13. The raw images were reduced to Level 1, using the HIPE software with standard photometric pipelines. The final maps in Level 2 were processed with the software Scanamorphos, version 24 (^{Roussel 2013}). Finally, we applied appropriate flux correction factors and color corrections to the PACS and SPIRE maps, respectively, according to the corresponding manuals⁴.

To extract the fluxes in each of the Herschel bands we applied the aperture photometry method, using the HIPE software. We used aperture sizes of 20'' and 25'' in the 70 and 160 μm bands, respectively, and of ≈ 36'' in all three SPIRE bands. The background was estimated over a ring with internal and external radii of ≈ 50'' and 90'', respectively, centered on the peak emission position. The photometric uncertainty assigned to each flux is given by the standard deviation of the measured values after aperture correction. These values are taken at different positions around the central source, without background subtraction, using the same aperture as for the flux measurement of the source (^{Balog et al. 2013}). For IRS 1, 4, 5 and 6 fluxes relative errors are between 10 - 25 %. The remaining sources (IRS 2, 3, and 7) are contaminated by the emission from the environment in which they are immersed and their fluxes for λ > 160 μm have relative errors of ≈ 50 %, preventing us from modeling these sources.

This region was also observed in the bands of the Wide-field Infrared Survey Explorer (WISE) satellite⁵ (^{Wright et al. 2010}), at W1(3.4 μm), W2(4.6 μm), W3(12 μm), and W4(22 μm), with FWHM of about 6.1'', 6.4'', 6.5'' and 12.0'', respectively. The WISE images and fluxes were obtained from the IRSA (NASA/IPAC Infrared Science Archive⁶) database.

2.2. Infrared Emission and Identification of YSOs

The bottom panel of Figure 1 is a composite image of the WISE images at 4.6 μm (blue) and 12 μm (green) and the Herschel image at 70 μm (red) of the region IRAS 08589−4714. A vast region dominated by the 70 μm emission is visible on the west. This emission is likely produced by warm dust, while on the east the emission is dominated by colder dust, detected at λ > 100 μm. A transition zone from a hotter to a cooler region is also apparent in this panel (see the highlighted area over the upper panel of Figure 1). A well-defined curved edge that crosses the region from north to south, seen in the W3(12 μm) WISE band (in green) in the bottom panel, is likely produced by the emission of warm dust and polycyclic aromatic hydrocarbons (PAHs; ^{Tielens 2008}). For more details on this curved and elongated structure see § 4.

In the bands of the Herschel telescope we detected 7 IR sources (IRS, labeled from 1 to 7 in the bottom panel of Figure 1), four of which (IRS 1, 2, 3, and, 7) were also detected in the bands of the WISE telescope. From the fluxes in these bands, we determined the color indexes of the sources. These four WISE detections are potential Class I or II objects, according to the criteria by ^{Koenig et al. (2012) 7}. The left panel of Figure 2 shows the position of these sources in the WISE color-color diagram. The sources without WISE counterparts (IRS 4, 5, and 6) are probably much younger objects (^{Chambers et al. 2009}). IRS 1 and 4 are projected onto two dust clumps (blue contours in Figure 1) identified by ^{Beltrán et al. (2006)} at 1.2 mm. Table 1 lists the sources detected in the WISE and Herschel images, with their coordinates and fluxes, and their correlation with the dust clumps detected at 1.2 mm.

Fig. 2 Left Panel: WISE W1−W2 vs W2−W3 color-color diagram. The dotted lines limit the region where Class I and Class II objects lie according to the ^{Koenig et al. (2012)}’s criteria. 

Green open squares on the left panel of Figure 2 correspond to faint WISE sources that lie in the region of the WISE color-color diagram contaminated by PAH emission, according to ^{Koenig et al. (2012) 8}. PAH emission has prominent lines in the WISE band W1 and W3, producing color excesses (^{Wright et al. 2010}). These sources are listed in Table 2. Interestingly, these objects are spatially located around the arc-shaped structures near IRS 1, 5 and 6. Figure 3 shows the bottom panel of Figure 1 (4.6 μm in blue, 12 μm in green, and 70 μm in red), indicating with red triangles the positions of these faint WISE sources. None of these sources has been detected in the Herschel bands, with the exception of those coincident with IRS 5 and 6. All sources in Table 2 are potential or candidate YSOs that require deeper observations to confirm their nature and evolutionary status.

The right panel of Figure 2 shows the 2MASS near-IR color-color J − H vs H − K diagram of the IRAS 08589−4714 region. IRS 2, 3 and 7 have near IR color excesses. Another three sources with near IR color excess lie in the direction towards IRS 1. All these 2MASS sources, located to the right of the reddening band, are listed in Table 3. The 2MASS sources coincident with IRS 5 and 6 do not have color excess in these bands (they are not included in Table 3). These sources are marked in Figure 3 with red squares. Color excesses in the near and middle IR could be associated with warmer dust in the inner part of the envelope, which could be due to a circumstellar disk (^{Sicilia-Aguilar et al. 2006}). This would confirm that these sources are young.

Fig. 3 Composite image of the WISE images at 4.6 μm (blue) and 12 μm (green), in combination with the Herschel image at 70 μm (red) of the IRAS 08589−4714 region. Red triangles show the positions of faint WISE sources listed in Table 2. Red squares correspond to sources with near-IR excesses listed in Table 3. The color figure can be viewed online. 

3.1. The Spectral Energy Distributions

In this Section, we analyze the evolutionary status of some of the sources identified in the Herschel images (see § 2.2) by means of their SEDs. To construct the SEDs we use the fluxes in the Herschel bands at 70, 160, 250, 350, and 500 μm. Since protostars lie embedded in dense and opaque material we can derive the main characteristics of the envelopes surrounding the central objects. Fluxes for λ < 20μm are excluded from the modeling of the whole set of sources since the code used (DUSTY, see § 3.3) underestimates the emission at these wavelengths. This is commonly observed in spherically symmetric models that do not include cavities or inhomogeneities in the envelope to allow radiation at these wavelengths to escape.

In the case of IRS 1, we model the SED and the 70 μm radial intensity profile (see § 3.2) simultaneously, while for IRS 4, 5, and 6 we model the corresponding SEDs only, since they are too faint to obtain reliable profiles at 70 μm. In the SED corresponding to IRS 1 we include the 1.2 mm flux from ^{Beltrán et al. (2006)} and the flux at 22 μm from WISE, while in the case of the SED for IRS 4 we only include the flux at 1.2 mm. This source does not have a mid-IR counterpart. The SED of IRS 1 is shown in the left panel of Figure 4, and those of IRS 4, 5, and 6, in Figure 5.

Fig. 4 SED (left panel) and intensity profile at 70 μm (right panel) of IRS 1. The filled diamonds correspond to the fluxes from the Herschel and WISE telescopes. 

Fig. 5 SEDs of IRS 4 (left top panel), 5 (right top panel), and 6 (left bottom panel). The filled diamonds belong to the five Herschel bands. The empty diamond in the SED for IRS 4 is an upper limit to the 70 μm flux. The error bars are included. The relative errors in the fluxes are between 10 and 25%. 

3.2. The 70 μm Intensity Profile for IRS 1

The average azimuthal intensity profile is commonly employed to compare modeled and observed images. This is used to complement the analysis of SEDs, as it allows to reduce the number of free parameters, and thus minimizes the degeneration of the solutions (^{Crimier et al. 2010}). We present the modeling of the intensity profile at 70 μm corresponding to IRS 1.

Since the PSF (point spread function) at 70 μm is rather elongated (due to the scan speed, which for 60''/s gives a PSF of 5.83'' x 12.12'', ^{Poglitsch et al. 2010}), to determine the radial intensity profile several one-dimension radial profiles are obtained in different directions from the center of the source. In the case of IRS 1 (see Figure 4), we use only the hemisphere in the opposite direction to the apex of the curved structure seen in Figure 1 to avoid contamination from the surrounding diffuse and extended emissions. The uncertainty in the profile is given by the noise in the image and the non-circularity of the beam of the source. The latter effect is taken into account by the standard deviation of the average azimuthal flux. The profile is normalized by its maximum value.

3.3. The Fitting Method

The SEDs of the four sources, as well as the 70 μm intensity profile of IRS 1, were modeled with the DUSTY⁹ code of ^{Ivezic & Elitzur (1997)}. The code solves the radiative transfer in one dimension for a dusty shell that surrounds a central source, which emits as a black body and whose radiation is absorbed, scattered and re-emitted by the dust in the envelope. This envelope has a density profile that follows a power law (n ∝ r^−p), and it is parameterized by the relative size given by Y = r_ext /_rin, where r_ext and r_in are the external and internal radii of the shell. The temperature of the central source (T_star) and the temperature at the internal radius of the envelope (T_in) are fixed, with values of 15000 and 300 K, respectively. The results are rather insensitive to the values of T_star.

We tested the DUSTY SEDs for T_star between 5000 and 50000K, and no distinguishable differences in the results were found. On the other hand, the inner radius was relatively unconstrained by current data resolution (5'' corresponds to ≈ 10 000 AU at 2 kpc), as is T_in. As a matter of fact, T_in determines the radius at which the code starts the calculations. Several SEDs models of massive stars, as well as mid-infrared spectra rule out inner radius temperatures consistent with dust sublimation temperatures (≈ 1500 K). Radiation pressure, stellar winds and/or shocks are usually invoked to produce a large cavity depleted of dust and to limit the amount of short wavelength emission (^{Churchwell et al. 1990}; ^{Faison et al. 1998}; ^{Hatchell et al. 2000}). For T_in we fixed a value in a manner similar to previous work, adopting an inner radius larger than would be expected from dust sublimation (^{Hatchell et al. 2000}; ^{Jørgensen et al. 2002}; ^{Crimier et al. 2009}, ²⁰¹⁰; ^{Vehoff et al. 2010}; ^{Hirsch et al. 2012}). Finally, for very young sources it is usual to assume that the envelope is composed of particles surrounded by a thin layer of ice. For this reason, we adopted the opacity corresponding to a density of 10⁶ cm⁻³ from ^{Ossenkopf & Henning (1994)}.

The DUSTY model needs to be scaled to the distance and the bolometric luminosity (L_bol) of the central source so as to compare the models and observed data. The bolometric luminosity is, however, an output parameter, since it is calculated by integrating the modeled SED. Consequently, the luminosity is re-evaluated as input parameter iteratively until the difference between the models and observations is minimized. We adopted a distance of 2.0 kpc (see § 1), assuming that all the sources belong to the IRAS 08589−4714 region.

On the other hand, the intensity profiles modeled by DUSTY show an intense and very narrow central component (width < 1''). This component is the stellar radiation attenuated by a dusty envelope, whose width is proportional to the stellar radius (^{Ivezić & Elitzur 1996}). Given that the pixel size at 70 μm is 3′′.2, we re-binned the model profile using a slightly larger bin of 4′′. Finally, this profile is convolved with the instrumental profile at 70 μm, taken as the average intensity profile of PACS calibrated sources, such as α Tau, the Red Rectangle, IK Tau and the Vesta asteroid (^{Lutz 2012}). The first two sources are used to model the core and the last two the extended wings of the profile (^{Aniano et al. 2011}).

To systematically compare the models with the observations we built a grid of 127 100 DUSTY models. The grid was constructed for 31 values of the power index (p) in the density profile, ranging from 0.0 to 3.0 in steps of 0.1, 82 values of the envelope relative size (Y) in the range 100-910 with steps of 10, and 50 values of the optical depth (τ₁₀₀), from 0.1 to 5.0, in steps of 0.1. The space of the parameters explored is similar to that used by ^{Crimier et al. (2010)}, restricted to a more limited range of each parameter but including typical values for young stellar objects. Table 4 summarizes the ranges of the free parameters. The best fit is obtained by means of the goodness-of-fit criterion. We employed the weighted least squares method to fit the SED. The weights are inversely proportional to the square of the errors assigned to each observed data point.

In the case of IRS 1, to find a model that fits both the SED and the intensity profile simultaneously we followed the procedure used by ^{Crimier et al. (2010)}. We first fixed the value of the opacity (τ₁₀₀) and, by fitting the intensity profile, we obtained an initial value for the power law index (p) and for the envelope relative size (Y). Then we recalculated the value of τ₁₀₀ by fitting the SED, using the values of p and Y derived in the previous step. We adopted the new value of τ₁₀₀ in the next iteration, and the process was repeated. Convergence was achieved in a few steps. This procedure considers the fact that the intensity profile is strongly dependent on the size of the envelope and the density distribution, whereas the SED depends mainly on the column density, that, in turn, depends on the opacity (^{Crimier et al. 2010}).

3.4. Derived Parameters

The solid red line in Figures 4 and 5 shows the best DUSTY models for IRS 1, 4, 5, and 6. With green continuous line we show all DUSTY models comprised within the error bars of the fluxes (20 in the case of IRS 1 and 200 for IRS 4, 5 and 6). Table 5 summarizes the parameters that best reproduce the data: the luminosity L, the opacity τ₁₀₀, the power law index p, and the relative size of the envelope Y. We include in this table the physical parameters of the envelope: the inner r_in and external, r_ext radii, the temperature of the dust at r_ext, T_dust, and the mass, M_env and density in the inner part of the envelope n(r_in).

Since the DUSTY code does not provide error determinations for the modeled parameters, Table 6 lists average values and standard deviations corresponding to all DUSTY models (calculated for the same luminosity) that fall within the error bars of the fluxes (green lines in Figs. 4 and 5). Best model parameters listed in Table 5 are not always within the standard deviations given in Table 6. However, they are not significantly different from the average values. In the case of IRS 1, p and Y are derived by modeling simultaneously the SED and the 70 μm intensity profile. For this reason, the differences between the values listed in Table 5 and 6 are larger.

For the SED of IRS 1, the fluxes at 3.4, 4.6, and 12 μm are not well fitted by the model, as we have mentioned before. These fluxes are likely contaminated by shocks and PAH emission features detected in these WISE bands (^{Wright et al. 2010}). Moreover, ^{Saldaño et al. (2016)} detected an outflow associated with IRS 1, which likely produces a cavity in the envelope. Dust in the walls of the cavity is heated by UV and optical radiation coming from the stellar embryo/s, re-irradiating in the IR (see, for example, ^{Bruderer et al. 2009}; ^{Zhang et al. 2013}). Thus mid-IR fluxes are considered to be upper limits. In addition, the observed intensity profile is well fitted by the modeled profile out to 20'', but not beyond.

The density profile for the envelope of IRS 1 derived with the DUSTY code is constant, and consistent with the supposition introduced by ^{Ghosh et al. (2000)} in their model of the SED of this source. In addition, the other parameters determined by these authors (M = 55 M_⊙, R_max = 0.2 pc, L = 2.4 ×10³ L_⊙) are similar or of the same order as the parameters obtained through DUSTY. However, we derived an opacity for the envelope four times larger than the value obtained by ^{Ghosh et al. (2000)}. The difference is probably due to the fact that these authors used silicate and graphite dust grains without a layer of ice on the surface, producing less attenuation to the radiation. The inclusion of a layer of ice gives a more realistic model (^{Ossenkopf & Henning 1994}).

3.5. Virial Mass

Assuming that: (1) the parent cloud is similar to that one used in the DUSTY code; (2) the source is in hydrostatic equilibrium; and (3) the volume density of the gas decays as a power law (ρ ∝ r^−p), we estimate the envelope mass using the following expression for the virial theorem:

where G is the gravitational constant and σ is the one-dimensional velocity dispersion averaged over the entire system. The k factor is approximately equal to 1 for Y ≫1 and p<2.5. Using the parameters calculated by the DUSTY code, such as p and r _ext , and 8 lin2σ , where ∆v is the full-width at half-maximum (FWHM) of the C¹⁸O line detected towards IRS 1 (^{Saldaño et al. 2016}), we estimate σ = 0.65 ± 0.05 km s⁻¹ and a virial mass for IRS 1 listed in Table 7. For the other three starless cores (IRS 4, 5 and 6) we adopt an average velocity dispersion ≈ 0.55 ± 0.25 km s⁻¹, averaging the values determined by ^{Sánchez-Mongeet al. (2013)} and ^{Battersby et al. (2014)} for multiple starless cores observed in NH₃ molecular transitions. Table 7 lists the calculated values. Sources 1, 5 and 6 have virial masses roughly in agreement with their DUSTY mass envelope, confirming that they are gravitationally bound. On the other hand, IRS 4 has a virial mass smaller than the envelope mass, indicating that this source may be collapsing.