Spatial Variations of Dust Opacity and Grain Growth in Dark Clouds: L1689, L1709 and L1712

The primary objective of this study is an investigation of dark clouds L1689, L1709 and L1712 within the Ophiuchus cloud. We use photometry in the $J$, $H$, and $K$ bands, obtained from the Galactic Clusters Survey (GCS) of the UKIRT Infrared Deep Sky Survey (UKIDSS) Release 11 (Lucas et al., 2008). The observed region encompasses an area of approximately 2.6°$\times$2.3°. Detailed descriptions of the photometric system and its calibration are available in Hewett et al. (2006) and Hodgkin et al. (2009), respectively.

The UKIDSS observations provide superior resolution and sensitivity in comparison to 2MASS survey, achieving limiting magnitudes that are about 3.5 mag deeper across each band. The UKIDSS data can be retrieved from the WFCAM Science Archive¹¹1http://surveys.roe.ac.uk/wsa/index.html. To ensure high data quality and keep only stellar sources, the photometric errors for all $JHK$ bands are required to be less than 0.2 mag and values of ‘PSTAR’ flag larger than 0.9. This results in limiting magnitudes of 20.9, 19.4, and 19.0 mag for the $J$, $H$, and $K$ bands, respectively.

We collected observations on Ophiuchus molecular cloud from the Spitzer Space Telescope, which is a part of the Spitzer Legacy Project “From Molecular Cores to Planet-Forming Disks” (c2d) (C2D Team, 2020). We use data across four mid-infrared (MIR) IRAC bands, namely [3.6], [4.5], [5.8], and [8.0]. Detailed information about data processing can be found in Evans et al. (2003) and Padgett et al. (2008). The Spitzer c2d observations of Ophiuchus cloud cover an area of about 8 square degrees. Access to the Spitzer c2d data is facilitated through the archive on the Spitzer Science Center’s website²²2https://irsa.ipac.caltech.edu/data/SPITZER/C2D/. To obtain a high-quality data sample, the signal-to-noise ratio (SNR) of photometry is required to be larger than 5, and the ‘QType’ flag to be exclusively ‘star’. Consequently, the photometric limiting magnitudes for the four IRAC bands were 17.9, 17.1, 16.3, and 15.2 mag, respectively.

The FIR continuum emission data are obtained from the Herschel Space Observatory (Pilbratt et al., 2010) by its two primary instruments: the Photodetector Array Camera and Spectrometer (PACS) (Poglitsch et al., 2010) and the Spectral and Photometric Imaging Receiver (SPIRE) (Griffin et al., 2010). The PACS data offer resolutions of 6^′′ at both 70 and 100 $\mu$m, and an approximate resolution of 12^′′ at 160 $\mu$m. In contrast, the SPIRE data provide spatial resolutions of 18.1^′′, 25.2^′′, and 36.9^′′ at 250, 350, and 500$\,{\rm\mu m}$, respectively. It is crucial to acknowledge that these resolutions can vary depending on the specific observation mode and the conditions under which the observations were made. The datasets are publicly accessible via the Herschel Science Archive³³3http://archives.esac.esa.int/hsa/whsa/.

The Ophiuchus cloud is populated with various translucent clouds, dense cores, and young stellar objects (YSOs).To probe the central dense regions of these clouds, we incorporated both UKIDSS NIR and Spitzer MIR photometric data. Following the methodology outlined by Román-Zúñiga et al. (2009), we employed a combined approach utilizing the Near-Infrared Color Excess (NICE) method proposed by Lada et al. (1994) and the NICER method revised by Lombardi & Alves (2001) to estimate extinction. Both techniques necessitate the prior knowledge of intrinsic colors. The NICE method requires only one observed color (i.e., the difference between photometric magnitude across two wavelength bands), while the NICER method needs two observed colors. Generally, we use the following formula to calculate the extinction in the $K$ band (i.e., $A_{K}$):

where $E(m_{\lambda_{1}}-m_{\lambda_{2}})$ is the color excess in the two bands of $\lambda_{1}$ and $\lambda_{2}$, which is the difference between the observed color $(m_{\lambda_{1}}-m_{\lambda_{2}})_{\rm obs}$ and the intrinsic color $(m_{\lambda_{1}}-m_{\lambda_{2}})_{0}$. The term $C_{\rm el}$ is the coefficient to convert color excess $E(m_{\lambda_{1}}-m_{\lambda_{2}})$ to extinction $A_{K}$, which depends on the assumed extinction law. We adopt the results of $C_{\rm el}$ from Li & Chen (2023) who studied the extinction law within the Ophiuchus cloud using the same data (i.e., UKIDSS and Spitzer data) as that of this work. The corresponding coefficients $C_{\rm el}$ are listed in Table 1.

The two photometry catalogs from UKIDSS and Spitzer are merged within a tolerance radius of 1 arcsec using the TOPCAT software. The UKIDSS-Spitzer common area catalog is analyzed separately and merged partially: only Spitzer sources are used when the UKIDSS sources are either unavailable or unreliable. The deatiled procedures are as follows: (1) if all three $JHK$ photometry of a source are available, the extinction is estimated using the NICER method; (2) if a source lacks $J$ photometry but has $HK$ photometry, the extinction is evaluated through $E(H-K)$ using the NICE method; (3) if a source lacks both $J$ and $H$ photometry, the extinction is estimated by $E(K-[3.6])$. This process is continued until only $E([4.5]-[5.8])$ is used.

Interstellar dust is homogeneously mixed with gas across a variety of environments, implying that the dust-to-gas ratio does not exhibit significant variation with respect to density or the sightline direction. The ratio of dust extinction or reddening to gas column density remains relatively constant across a wide range of conditions (Zhu et al., 2017). According to Equation 4, by correlating the two independent measures of $A_{K}$ and $\tau_{250}$ outlined in Section 3, and with the employment of an empirically calibrated ratio that connects $A_{K}$ to $N_{\rm H}$, the dust opacity $r\kappa_{250}$ can be determined as follows:

which is dependent on the empirical ratio of $N_{\rm H}/A_{K}$. In the typical diffuse ISM, the ratio of $N_{\rm H}$ to $A_{V}$ is estimated to be approximately $1.9\times 10^{21}\,\rm cm^{-2}\,mag^{-1}$ (Bohlin et al., 1978). By applying the extinction law with $R_{V}=3.1$ from the models of Weingartner & Draine (2001) (hereafter WD01), where $A_{K}/A_{V}=0.11$, the corresponding $N_{\rm H}/A_{K}$ is referred to be about $1.7\times 10^{22}\,\rm cm^{-2}\,mag^{-1}$. For dense clouds, Vuong et al. (2003) reports $N_{\rm H}/A_{J}\approx(5.6-7.2)\times 10^{21}\,\rm cm^{-2}\,mag^{-1}$ toward $\rho$ Ophiuchi cloud, while Martin et al. (2012) gives $N_{\rm H}/E(J-K)=(11.5\pm 0.5)\times 10^{21}\,\rm cm^{-2}\,mag^{-1}$ in the Vela cloud. Using the ratio $A_{J}/A_{K}=2.5$ (Indebetouw et al., 2005), these ratios suggest an $N_{\rm H}/A_{K}$ range of approximately $(1.4-1.8)\times 10^{22}\,\rm cm^{-2}\,mag^{-1}$. Given this context, $N_{\rm H}/A_{K}=1.7\times 10^{22}\,\rm cm^{-2}\,mag^{-1}$ is adopted in this work.

Figure 6 illustrates a clear correlation between $A_{K}$ and $\tau_{250}$ across the entire region, after excluding data points with SNRs below 3. Notably, there is a larger scatter for $A_{K}>0.6$ mag, yet the overall correlation between $A_{K}$ and $\tau_{250}$ remains apparent. This increased dispersion at higher extinctions may suggest greater variations in dust opacity along the line of sight, particularly toward the densest core regions. To estimate the overall dust opacity $r\kappa_{250}$, a linear regression is conducted on the observed correlation between $\tau_{250}$ and $A_{K}$. For this purpose, the Python package lts_linefit⁶⁶6https://pypi.org/project/ltsfit/(Cappellari et al., 2013) was used, which takes into account errors in both variables, potential outliers and intrinsic scatter. As shown in Figure 6, the grey points represent the outliers as determined by lts_linefit. The regression yielded a slope of $\tau_{250}/A_{K}=0.0037$, corresponding to a mean opacity value of $r\kappa_{250}=0.093\,{\rm cm^{2}\,g^{-1}}$​ throughout the entire region. Considering a power-law variation of $\tau_{250}$ and $A_{K}$ found in the literature, we fitted a power-law function to all data points, yielding the relationship $\tau_{250}=0.004{A_{K}}^{1.20}+0.003$. This result demonstrates that $\tau_{250}$ increases with $A_{K}$ (or $N_{\rm H}$), to a power of approximately 1.2, implying that $r\kappa_{250}$ increases with $A_{K}$ to the power of 0.2 (see Equation 5). This dependency indicates potential grain growth as the cloud density increases. In comparative literature, Roy et al. (2013) reported a power-law index of $1.28\pm 0.01$ in Orion A, while Okamoto et al. (2017) discovered an index of 1.32$\pm$0.04 for Perseus. More recently, Hayashi et al. (2019) reported a value of 1.21$\pm$0.04 for Chamaeleon. Although the power-law index of 1.20 derived in this work is close to these findings in other molecular clouds, the power-law relationship between $\tau_{250}$ and $A_{K}$ is relatively weaker in comparison to a linear relationship.

It is notable that these studies as mentioned above primarily used 2MASS extinction maps that have a comparatively low resolution of approximately $3^{\prime}-5^{\prime}$, rendering the resolving of smaller structures challenging. Despite our extinction mapping achieving a threefold improvement in resolution, the overall results exhibit no significant alterations. A potential cause could be that the large dispersion among data points with high extinction leads to inaccuracies in the fitting results. One might hypothesize that the resolution of finer details would uncover more pronounced density gradients, potentially leading to larger deviations in the $\tau_{250}$ and $A_{K}$ relationship from a simple linear approximation. However, the continuous presence of the shallow power-law, even at higher resolutions, suggests a smooth transition from diffuse to dense conditions over scales larger than those captured by our current observations. Another consideration involves beam dilution in both FIR emission and extinction maps which could average out inherent small-scale variations. Further high-resolution mapping of approximately 10^′′ would reveal localized opacity variations. Nevertheless, the consistency between our $\tau_{250}$-$A_{K}$ results and those from previous lower-resolution studies suggests that the inferred dust evolutionary trends remain reliable across an extensive scale range.

Considering the complexity of various environments within the Ophiuchus molecular cloud, the entire region is divided into 16 distinct sub-regions to analyze the spatial variations. The positions of these sub-regions are illustrated in Figure 5. The typical sizes of these sub-regions are $\sim 15^{\prime}$, corresponding to a physical scale of $\sim$0.6 pc. In order to derive the opacity $r\kappa_{250}$ via Equation 5, a linear fit is performed to the $\tau_{250}$ versus $A_{K}$ relation in each sub-region by lts_linefit. This fitting approach provides robust $\tau_{250}$-$A_{K}$ correlations for each sub-region. The resultant FIR opacities $r\kappa_{250}$ are listed in Table 3, spanning a range from 0.07 to 0.12$\,{\rm cm^{2}\,g^{-1}}$. Table 3 also presents the derived $\sigma_{250}$, along with the mean dust temperature $\langle T_{\rm d}\rangle$ and the average extinction $\langle A_{K}\rangle$ for each sub-region.

Dust opacity can be expressed in a variety of forms across diverse wavelengths. In the present work, we represent the opacity in terms of the opacity at 250 $\mu\rm m$, denoted as $r\kappa_{250}$. In Section 4.1, we derived an average opacity $r\kappa_{250}\approx 0.093\,{\rm cm^{2}\,g^{-1}}$ encompassing the entire region under investigation. When analyzing the sub-regions individually, we observed $r\kappa_{250}$ values ranging from $0.07\,{\rm cm^{2}\,g^{-1}}$ to $0.12\,{\rm cm^{2}\,g^{-1}}$. The diffuse ISM at high latitudes typically exhibits an opacity around $\tau_{250}/N_{\rm H}\approx\rm 1\times 10^{-25}\,cm^{2}\,H^{-1}$ (Boulanger et al., 1996). Observations from the Planck mission suggest $\tau_{250}/N_{\rm H}$ values within the range of $\rm(0.6-1.6)\times 10^{-25}\,cm^{2}\,H^{-1}$ (Planck Collaboration et al., 2011, 2014b, 2014c), which correspond to $r\kappa_{250}$ values between 0.025 and 0.068$\,{\rm cm^{2}\,g^{-1}}$. These findings align with the predicted $r\kappa_{250}=0.047\,{\rm cm^{2}\,g^{-1}}$ from the dust model of Draine & Li (2007). The mean value derived in this work is approximately twice the value for the diffuse ISM. In regions 2 and 4 with lowest extinction of $A_{K}<0.2$ mag, we determined $r\kappa_{250}\approx 0.07\,{\rm cm^{2}\,g^{-1}}$, larger than that in the diffuse ISM. The NIR extinction map in this work is insensitive to the diffuse region with $A_{K}<$ 0.1 mag. As indicated in Figure 3, the sensitivity of the derived extinction map, denoted as $\sigma_{A_{K}}$, ranges from 0.02 to 0.1 mag.

Instead of extinction mapping, Suutarinen et al. (2013) employed the use of smoothed color excess $E(H-K)$ derived from NTT/SOFI observations targeting the prestellar core CrA C. They found an opacity of $r\kappa_{250}=0.08\pm 0.01\,{\rm cm^{2}\,g^{-1}}$. This value does not demonstrate a significant enhancement when compared to the diffuse ISM, to some extent, it is less than most of the results obtained in this work.

Lewis et al. (2022) conducted a comparison between the 2MASS extinction map with Planck emission at a resolution of 5 arcmin for the Ophiuchus cloud. Their derived opacity $r\kappa_{250}\approx 0.07\,{\rm cm^{2}\,g^{-1}}$ is lower than our average value and even the minimum value among the sub-regions. It is worth noting that the adopted extinction law could influence the resultant opacity. Flatter curves tend to yield higher extinctions from color excess measurements. For instance, varying between the $R_{V}$=3.1 and $R_{V}$=5.5 extinction laws from WD01 introduces uncertainties up to $\sim$20% in $A_{K}$ and consequently in $r\kappa_{250}$. In addition, a lower spectral index $\beta$ assumed in SED fitting results in lower $\tau_{250}$, which subsequently lead to decreasing $r\kappa_{250}$ by approximately 25% for $\beta$=1.8 compared to $\beta$=2.2 (Suutarinen et al., 2013). However, the primary factor contributing to the lower opacity in Lewis et al. (2022) is the limited spatial resolution. Their resolution of 5 arcmin ($\sim$0.2 pc) biases their opacity estimate towards nearby translucent ISM. As a result, the $r\kappa_{250}$ value that we derive is more reflective of values in dense environments where dust grain evolution is expected to be significant. Our results emphasize the necessity of high resolution in both extinction and emission maps to accurately characterize opacity variations related to local density enhancements.

Dust grains at various stages of evolution display distinct FIR opacity characteristics. Specifically, an increase in FIR opacity often indicates grain growth, as larger dust grains contribute more significantly to the opacity at longer wavelengths (Köhler et al., 2012, 2015; Ossenkopf & Henning, 1994; Ormel et al., 2011). The grain growth is also often observed in denser regions of the ISM and in protoplanetary disks, where dust grains can accumulate and coagulate (Testi et al., 2014). According to the models proposed by Ossenkopf & Henning (1994), ISM-like dust without ice mantles and coagulation has an opacity $r\kappa_{250}\sim 0.06$ $\,{\rm cm^{2}\,g^{-1}}$, which increases to 0.13 $\,{\rm cm^{2}\,g^{-1}}$ when dust grains evolve over a coagulation time of $10^{5}$ years at a gas density of $10^{5}$ $\rm cm^{-3}$. Under similar conditions, the updated models by Ormel et al. (2011) predict opacities reaching 0.14-0.17 $\,{\rm cm^{2}\,g^{-1}}$ for fully ice-coated aggregates (ic-sil, ic-gra) or spatially mixed aggregates (is-sil+gra). Moreover, with greater coagulation time and density, both models can yield even higher opacities (up to 0.2 $\,{\rm cm^{2}\,g^{-1}}$ or more). Therefore, we deduce that the opacities observed in the Ophiuchus cloud can be consistent with these models under specific conditions.

Numerous studies have previously identified the trend of increasing FIR dust opacity with higher environment density by correlating the FIR optical depth with extinction. For instance, Kramer et al. (2003) observed increasing opacity with decreasing temperature in IC5146. Ysard et al. (2013) found a factor of two higher opacity in the center compared to the outskirts of the L1506 filament in Taurus. Stepnik et al. (2003) utilizing PRONAOS observations of Taurus measured a 3.5 times increase in $\tau_{250}/N_{\rm H}$. Schnee et al. (2008) also reported higher emissivities for colder grains. Dust coagulation can account for the opacity enhancements and temperature drops in dense regions. Radiative transfer models indicate this opacity rise cannot be explained by changes in the radiation field alone, but instead reflects variations in intrinsic dust properties (Juvela et al., 2015; Ysard et al., 2013; Suutarinen et al., 2013). While many observations indicate grain growth in dense regions, Suutarinen et al. (2013) did not find clear opacity enhancement in a prestellar core in CrA. This suggests that opacity alone does not definitively trace evolution. Therefore, complimentary constraints from dust scattering properties are also needed to quantify grain growth (Juvela et al., 2020).

Figure 7 presents the correlations between $\tau_{250}$ and $A_{K}$ across all 16 sub-regions as defined in Figure 5. The majority of these sub-regions exhibit a tight linear relationship of $\tau_{250}$ and $A_{K}$. However, in some regions with higher dynamic range of $A_{K}$ (e.g. regions 11-16), the lts_linefit linear fitting technique has detected a significant number of outliers. These outliers primarily deviate above the best-fit linear relationship at higher $A_{K}$ values, which may suggest a potential increase in dust opacity that could be indicative of grain growth. It is also plausible that these outliers are attributed to the underestimation of extinction toward the dense cores, potentially due to a limited number of background stars. Moreover, these dense cores are too small to be resolved with the 1 arcmin beam in our extinction mapping, leading to larger dispersion when $A_{K}>1.5$ mag. The localized opacity increase toward dense cores merits further investigation through high angular resolution observations that can resolve the small-scale structure. To distinguish potential opacity variations from mapping artifacts, it is crucial to match the resolution between the FIR emission and NIR extinction measurements.

As shown in Figure 8, the overall variation in opacity across the region is within a factor of two. Nevertheless, there is a discernible general trend of increasing opacity with dust column density, particularly towards cloud cores where grain growth appears enhanced. Despite the complexities introduced by including both lower density and dense star-forming cores, the opacities remain relatively uniform, only exhibiting a factor of two variation across the entire region. The lack of strong opacity enhancements even toward dense cores implies rapid grain growth has not yet dramatically altered the dust properties in most areas. The handful of sub-regions with $r\kappa_{250}$ exceeding 0.1$\,{\rm cm^{2}\,g^{-1}}$ (see Figure 7 and 8) likely reflect grain coagulation in the cold dense cores. But overall, the narrow distribution suggests the dust populations in the Ophiuchus dark clouds are likely at similar evolutionary stages, with no clear signs of rapid grain growth.