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

The far-infrared (FIR) opacity of dust in dark clouds within the Ophiuchus molecular cloud is investigated through multiwavelength infrared observations from UKIDSS, Spitzer, and Herschel. Employing the infrared color excess technique with both near-infrared and mid-infrared photometric data, a high-resolution extinction map in the K band (A K ) is constructed for three dark clouds: L1689, L1709, and L1712. The derived extinction map has a resolution of 1′ and reaches a depth of A K ∼ 3 mag. The FIR optical depths τ 250 at a reference wavelength of 250 μm are obtained by fitting the Herschel PACS and SPIRE continuum data at 100, 160, 250, 350, and 500 μm using a modified blackbody model. The average dust opacity per unit gas mass at 250 μm, r κ 250, is determined through a pixel-by-pixel correlation of τ 250 with A K , yielding a value of approximately 0.09 cm2 g−1, which is about 2–3 times higher than the typical value in the diffuse interstellar medium. Additionally, an independent analysis across 16 subregions within the Ophiuchus cloud indicates spatial variations in dust opacity, with values ranging from 0.07 to 0.12 cm2 g−1. Although the observed trend of increasing dust opacity with higher extinction implies grain growth, our findings indicate that rapid grain growth has clearly not yet occurred in the dark clouds studied in this work.


INTRODUCTION
Dark clouds within the Milky Way represent a subset of molecular clouds that remain visually opaque due to their high gas and dust density (Lynds & Sandage 1963).These clouds, particularly their cores, are often the birthplaces of new stars and accompanying planetary systems.Precise knowledge of the initial conditions, including density and temperature, within these dark clouds can considerably enhance our understanding of star formation processes.However, molecular hydrogen, a major component in star formation, remains undetectable under conditions of low temperature and high density.Additionally, molecules such as CO are often depleted within the cores of dark clouds, and the associated molecular emission (i.e. 12 CO) is usually optically thick.As a result, the near-infrared (NIR) extinction of dust and its far-infrared (FIR) emission serve as important indicators of physical conditions within the dark cloud cores prior to its collapse into a star (Lada et al. 2007).Moreover, the study of dust emission properties is of significant value in obtaining more reliable mass estimates for molecular clouds and further refining the initial mass function (IMF) of star formation.Consequently, the investigation of dust properties within dark clouds is fundamental to understanding the processes governing star formation.
A fundamental property of dust is their ability to absorb short-wavelength photons and re-emit energy at longer wavelengths (Li 2005).This absorption and emission behavior is commonly characterized by the 'opacity' parameter κ ν , expressed in units of cm 2 g −1 .This parameter plays a crucial role in estimating dust mass by fitting the FIR and millimeter continuum emission (Fanciullo et al. 2020;Li et al. 2022).However, achieving accurate measurements has been constantly challenging.Historical measurements of the opacity have shown discrepancies spanning more than two orders of magnitude (Clark et al. 2019;Bianchi et al. 2022), and spatial variations have been observed across different environments (Roy et al. 2013;Martin et al. 2012).Moreover, the widely adopted models for dust evolution studies, such as those proposed by Ossenkopf & Henning (1994) and Ormel et al. (2011), present absorption coefficients that differ by over threefold.The assumptions underpinning these models may not be universally applicable across complex interstellar environments, ranging from cold, starless cloud cores to heated region in star clusters.Recent studies have suggested that the opacity provided by dust evolution models does not adequately account for the FIR emission of cloud cores (Webb et al. 2017).As a result, it becomes essential to constrain the opacity through observational data.
The opacity can be investigated by considering the ratio of extinction in the optical/NIR bands to FIR emission (Bianchi et al. 2003;Kramer et al. 2003;Shirley et al. 2011;Suutarinen et al. 2013).On the one hand, the Herschel Space Observatory (Pilbratt et al. 2010) provides photometric observations of dust emission spanning wavelengths between 70-500 µm.This enhances our ability to probe the physical properties of dust at relatively high angular resolution and across great depths in molecular clouds, although the results are somewhat mitigated by degeneracies in dust optical depth and temperature.On the other hand, extinction is independent of temperature and serves as an effective proxy for tracing the column density of dust.However, the construction of an extinction map relies heavily on the statistical analysis of the number density of background stars.As a result, in dense dark cloud cores, the measurement of extinction maps requires high sensitivity and photometric depth (Webb et al. 2017).Previous extinction maps derived from the 2MASS survey typically have a spatial resolution of around 3 ′ (Lombardi et al. 2008;Juvela & Montillaud 2016), which is insufficient for resolving the internal structure of dark cloud cores.Although many studies have reported an increase in opacity with increasing extinction, opacity in dense cloud cores tends to be underestimated due to the limited depth of A V < 10 mag (Martin et al. 2012;Roy et al. 2013;Juvela et al. 2015;Lewis et al. 2022).To better study the dust properties of clouds and cores, it is essential to construct an extinction map with a spatial resolution on par with Herschel observations (i.e.< 1 ′ ).The deeper and more sensitive NIR surveys from Visible and Infrared Survey Telescope for Astronomy (VISTA) and the United Kingdom Infrared Telescope (UKIRT) can help to derive such high-resolution extinction maps (Chu & Hodapp 2021;Zhang & Kainulainen 2022).
In this work, we aim to investigate dust opacity by correlating NIR extinction with FIR emission within the Ophiuchus cloud, focusing specifically on the three dark clouds L1689, L1709, and L1712.These clouds are located at a distance of ∼139 pc (Ortiz-León et al. 2018) with relatively clean background, making it a suitable target for studying dust evolution in molecular cloud.Previous works have shown evidence of dust growth in Ophiuchus cloud based on the infrared extinction law (Chapman et al. 2009;Li & Chen 2023).Here, we focus on the FIR dust opacity by combining NIR extinction and FIR emission.
This paper is organized as follows.Section 2 describes the data briefly.Section 3 presents the data analysis approach to derive extinction map and FIR optical depths.Section 4 provides the determination of dust opacity and discussions.Section 5 is the summary.

UKIDSS Data
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°×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 1 .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.

Spitzer Data
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 2 .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.

Herschel Data
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 µm, and an approximate resolution of 12 ′′ at 160 µm.In contrast, the SPIRE data provide spatial resolutions of 18.1 ′′ , 25.2 ′′ , and 36.9 ′′ at 250, 350, and 500 µ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 3 .

NIR Extinction mapping
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 λ1 − m λ2 ) is the color excess in the two bands of λ 1 and λ 2 , which is the difference between the observed color (m λ1 − m λ2 ) obs and the intrinsic color (m λ1 − m λ2 ) 0 .The term C el is the coefficient to convert color excess E(m λ1 − m λ2 ) to extinction A K , which depends on the assumed extinction law.We adopt the results of C 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 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.

Intrinsic Colors
2 https://irsa.ipac.caltech.edu/data/SPITZER/C2D/ 3http://archives.esac.esa.int/hsa/whsa/Normally, intrinsic colors are approximated by averaging the observed color of background stars from a nearby lowextinction region (also known as a control field).Nonetheless, selecting a low-extinction region from Spitzer observation is complicated because it covers the main structure of the Ophiuchus cloud.To address this, we choose to leverage the TRILEGAL model in simulating the stellar distribution in the Milky Way4 (Girardi et al. 2005).In the TRILEGAL simulations, we designated a field of view with 2 square degrees directed towards the Ophiuchus cloud, set the extinction to zero, and keep all other parameters (i.e. the initial mass function, binary fraction, and Galaxy component) at their default values.We then selected the desired photometric system and set the H band limiting magnitude to be 20 mag.Consequently, a stellar catalog was generated containing details such as stellar mass, temperature, and magnitudes in the corresponding bands.We then computed the average and 1σ dispersion of model colors as the estimates of the intrinsic colors and their associated uncertainties, respectively.For comparative purposes, we selected a lowextinction region from the UKIDSS observations, and estimated the average colors (J − H) 0 ≈ 0.495 ± 0.293 mag and (H − K) 0 ≈ 0.163 ± 0.266 mag.These values closely match the estimates from TRILEGAL simulations within 1σ uncertianty, as revealed in Table 1.This indicates that the TRILEGAL simulation can yield a reliable approximation of intrinsic colors.
Figure 1 presents the color-color diagram (CCD) of observed J − H versus H − K for three dark clouds: L1689, L1709, and L1712.This figure also exhibits the distribution of stars as derived from the TRILEGAL simulations, as well as stars from the 'low-extinction' field of the UKIDSS dataset.A noticeable shift of the stars from the dark clouds along the reddening vector in the CCD provides a strong indication of significant extinction present within these clouds.

Map Construction
Upon calculating the color excess and extinction A K for each individual star, we employed A K to generate an extinction map.The final extinction map is produced by utilizing a Gaussian filter with a Full Width at Half Maximum (FWHM) of 1 arcmin and a pixel size of 0.5 arcmin.The extinction value at each pixel is determined by the Gaussian weighted mean of extinction values for sources located within a radius three times the standard deviation (3σ) of the Gaussian beam.Thereafter, we calculated the errors associated with the mean extinction, represented as σ A K , using the conventional formula for a weighted average.
Figure 2 shows the resultant extinction map for the dark clouds L1689, L1709, and L1712 with a color scale showing A K up to 1.6 map. Figure 3 illustrates the relationship between the extinction error (σ A K ) and A K .For A K < 0.5 mag, σ A K is relatively constant at approximately 0.02 mag, but it begins to increase for A K > 0.5 mag, reaching up to roughly 0.2 mag.Although the depth of extinction map reaches A K ≈ 3 mag, it should be noted that regions with extinction A K > 2 mag are likely not resolved in Figure 2. Previous measurements of the Ophiuchus extinction map derived from 2MASS data typically demonstrated a resolution of approximately 3 arcmin (e.g., Lombardi & Alves 2001;Juvela & Montillaud 2016).The extinction map presented in this study improve the resolution by a factor of 3, thereby enabling the detection of core structures within dark clouds and filaments.

Zero-point Correction
Due to the unknown thermal background of the instruments on Herschel, the flux measurements obtained from Herschel represent relative fluxes rather than absolute fluxes.In this work, we used the Level 2.5 data products for the Herschel PACS at 100 µm and 160 µm, as well as Level 3 data products for the Herschel SPIRE at 250 µm, 350 µm  (2) The intercept a ν from the linear fit offers an estimate for the absolute photometric calibration of Herschel observations, which represent the zero point of Herschel flux.Ideally, the slope b ν of the linear fit should be close to 1.We estimate the errors of a ν using a bootstrap method and present the results for the observations of L1689, L1709, and L1712 in Table 2.There are substantial zero-point offsets in the PACS 100 µm and 160 µm images, ranging from several tens to hundreds of MJy/sr.This method has also been applied to the SPIRE data that have been previously corrected for zero-point offsets.The derived zero-point values for the SPIRE data were found to be near zero (less than 6 MJy/sr), demonstrating the effectiveness of this method in zero-point calibration.

Maps of Optical Depth and Dust Temperature
The zero-point corrected Herschel images were convolved to a resolution of 1 arcmin and regridded to a pixel size of 0.5 arcmin, consistent with the A K extinction map obtained in Section 3.1.Given that dust thermal emission at the FIR and submillimeter wavelengths is typically optically thin, the surface brightness I ν (in units of MJy/sr) of dust emission can be described by: where B ν (T d ) is the Planck function for dust temperature T d and frequency ν. τ ν is the dust optical depth, which can be expressed as: where N H is the total hydrogen column density N H = 2N (H 2 ) + N (H), µ is the mean molecular weight (with µ = 1.4,assuming 10% He), and r is the dust-to-hydrogen mass ratio M d /M H , which is presumed to be 1/100 and believed to vary minimally from diffuse to dense environment.κ ν is the dust opacity, frequently refered to as the mass absorption coefficient or sometimes the emissivity.Another commonly characterized quantity is σ ν ≡ τ ν /N H = µm H rκ ν , also called the opacity (the emission cross-section per H nucleon).In subsequent discussions, dust opacity is expressed as rκ ν , the opacity per unit gas mass.This can be converted to opacity per unit dust mass by simply dividing by the dust-to-hydrogen mass ratio r.

Li et al.
Adopting λ 0 = 250 µm (ν 0 = 1200 GHz) as the reference wavelength, Equation 3 can be rewritten as a modified blackbody: I ν = B ν (T d )τ 250 (ν/1200 GHz) β .Here we fix β at 2 while treating T d and τ 250 as free parameters.A χ 2 minimization technique is then applied to fit the SED for each pixel in the calibrated Herschel maps.The errors in SED fitting parameters are estimated using the Monte Carlo (MC) method.All SEDs are weighted by σ −2 , where σ is the sum of the squares of (i) the uncertainty of the zero correction obtained from Planck model, (ii) the mean Root Mean Square (RMS) noise at a resolution of 1 arcmin (RMS 1 ′ ) estimated from the Herschel noise maps.Calibration errors and color correction errors are accounted for in the MC analysis.Table 2 lists the relevant parameters used in SED fitting, including flux calibrations, color corrections, zero point (a ν ) and RMS 1 ′ .A 7% uncertainty in flux calibration for all five bands is adopted according to the PACS and SPIRE manuals (PACS Observer's Manual, Version 2.4; SPIRE Observers'Manual, Version 2.4).The color corrections and uncertainties are refered to Pezzuto et al. (2012) and Sadavoy et al. (2013).Figure 4 illustrates a representative SED and best-fit model.Figure 5 shows the optical depth (τ 250 ) and dust temperature (T d ) maps for the dark clouds L1689, L1709 and L1712.Generally, regions with lower temperatures correlate with regions of higher optical depths.The derived typical values of the 1σ relative errors of T d and τ 250 are 3.5 % and 15 %, respectively.

Determination of Opacity
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 τ 250 outlined in Section 3, and with the employment of an empirically calibrated ratio that connects A K to N H , the dust opacity rκ 250 can be determined as follows: which is dependent on the empirical ratio of N H /A K .In the typical diffuse ISM, the ratio of N H to A V is estimated to be approximately 1.9 × 10 21 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 H /A K is referred to be about 1.7 × 10 22 cm −2 mag −1 .For dense clouds, Vuong et al. (2003) reports N H /A J ≈ (5.6 − 7.2) × 10 21 cm −2 mag −1 toward ρ Ophiuchi cloud, while Martin et al. (2012) gives N H /E(J − K) = (11.5 ± 0.5) × 10 21 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 H /A K range of approximately (1.4 − 1.8) × 10 22 cm −2 mag −1 .Given this context, N H /A K = 1.7 × 10 22 cm −2 mag −1 is adopted in this work.Figure 6 illustrates a clear correlation between A K and τ 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 τ 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κ 250 , a linear regression is conducted on the observed correlation between τ 250 and A K .For this purpose, the Python package lts linefit6 (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 τ 250 /A K = 0.0037, corresponding to a mean opacity value of rκ 250 = 0.093 cm 2 g −1 throughout the entire region.Considering a power-law variation of τ 250 and A K found in the literature, we fitted a power-law function to all data points, yielding the relationship τ 250 = 0.004A K 1.20 + 0.003.This result demonstrates that τ 250 increases with A K (or N H ), to a power of approximately 1.2, implying that rκ 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 ± 0.01 in Orion A, while Okamoto et al. (2017) discovered an index of 1.32±0.04for Perseus.More recently, Hayashi et al. (2019) reported a value of 1.21±0.04for 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 τ 250 and A K is relatively weaker in comparison to a linear relationship.Note-(a) The positions of sub-regions are outlined in Figure 5.
It is notable that these studies as mentioned above primarily used 2MASS extinction maps that have a comparatively low resolution of approximately 3 ′ − 5 ′ , 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 τ 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 τ 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 ∼ 15 ′ , corresponding to a physical scale of ∼0.6 pc.In order to derive the opacity rκ 250 via Equation 5, a linear fit is performed to the τ 250 versus A K relation in each sub-region by lts linefit.This fitting approach provides robust τ 250 -A K correlations for each sub-region.The resultant FIR opacities rκ 250 are listed in Table 3, spanning a range from 0.07 to 0.12 cm 2 g −1 .Table 3 also presents the derived σ 250 , along with the mean dust temperature ⟨T d ⟩ and the average extinction ⟨A K ⟩ for each sub-region.

Comparison with Previous Studies
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 µm, denoted as rκ 250 .In Section 4.1, we derived an average opacity rκ 250 ≈ 0.093 cm 2 g −1 encompassing the entire region under investigation.When analyzing the sub-regions individually, we observed rκ 250 values ranging from 0.07 cm 2 g −1 to 0.12 cm 2 g −1 .The diffuse ISM at high latitudes typically exhibits an opacity around τ 250 /N H ≈ 1 × 10 −25 cm 2 H −1 (Boulanger et al. 1996).Observations from the Planck mission suggest τ 250 /N H values within the range of (0.6 − 1.6) × 10 −25 cm 2 H −1 (Planck Collaboration et al. 2011Collaboration et al. , 2014b,c),c), which correspond to rκ 250 values between 0.025 and 0.068 cm 2 g −1 .These findings align with the predicted rκ 250 = 0.047 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κ 250 ≈ 0.07 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 σ 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κ 250 = 0.08±0.01cm 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κ 250 ≈ 0.07 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 ∼20% in A K and consequently in rκ 250 .In addition, a lower spectral index β assumed in SED fitting results in lower τ 250 , which subsequently lead to decreasing rκ 250 by approximately 25% for β=1.8 compared to β=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 (∼0.2 pc) biases their opacity estimate towards nearby translucent ISM.As a result, the rκ 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.

Opacity as a Tracer of Dust Evolution
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(Köhler et al. , 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κ 250 ∼ 0.06 cm 2 g −1 , which increases to 0.13 cm 2 g −1 when dust grains evolve over a coagulation time of 10 5 years at a gas density of 10 5 cm −3 .Under similar conditions, the updated models by Ormel et al. (2011) predict opacities reaching 0.14-0.17cm 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 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 τ 250 /N 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 τ 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 τ 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 starforming 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κ 250 exceeding 0.1 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.

SUMMARY
In this work, we use multi-wavelength infrared observations from UKIDSS (1-3 µm), Spitzer (3-8 µm), and Herschel (100-500 µm) to expore dust opacity in three dark clouds (L1689, L1709, and L1712) within the Ophiuchus cloud complex.Our primary conclusions are as follows: 1.By integrating the NICE and NICER extinction measurement methods and leveraging NIR and MIR photometry, we generate extinction maps for the three dark clouds at a resolution of 1 arcmin (equivalent to 0.04 pc) and depth up to A K ∼ 3 mag.This resolution and depth exceed previous of prior 2MASS-based extinction maps by more than threefold, thereby enabling the resolution of distinct dense cores and filaments within the clouds.
2. Comparing the A K extinction with Herschel FIR emission allows us to determine the overall dust opacity of the Ophiuchus cloud to be rκ 250 ≈ 0.09 cm 2 g −1 on average, which is approximately 2-3 times higher than that of the diffuse ISM.This is consistent within errors with previous measurements based on 2MASS extinction map (e.g.Martin et al. 2012;Roy et al. 2013).It can be expected that observations by the James Webb Space Telescope (JWST), combined with submillimeter observations will extend this method and determine the dust opacity of interstellar clouds with unprecedented precision.
3. It is found that the dust opacity rκ 250 increases with density (A K or N H ) as a power law with an index of ∼0.2, suggesting grain growth.An analysis of sub-regions indicates that rκ 250 also exhibits spatial variations, generally increasing with higher environmental density.Although there are indications of grain growth, the evidence does not suggest rapid grain growth has taken place in the dark clouds examined in this study.
We would like to thank the anonymous referee for the very helpful comments and suggestions.We thank Dr.      Figure 8.The opacity rκ250 as a functin of average extinction ⟨AK ⟩ for the 16 sub-regions listed in Table 3.The dashed line represents a linear fit to the data, offering direct evidence of grain growth.

ν
and the observed Herschel fluxes convolved to the 5 arcmin resolution, S Herschel ν : S Herschel ν = a ν + b ν S Planck ν .

Figure 1 .Figure 2 .
Figure 1.Color-color density distribution (J − H vs. H − K) for stars within the regions L1689, L1709, and L1712.Stars from the control (zero-extinction) region and the TRILEGAL simulation are shown by grey and red markers, respectively.The blue star symbol represents the adopted intrinsic color (J − H)0 and (H − K)0.The black solid line presents the reddening law as measured by Li & Chen (2023).For comparison, the reddening laws predicted by WD01 RV = 3.1 and RV = 5.5 models (Weingartner & Draine 2001) are presented by the grey dashed and dot-dashed lines, respectively.

Figure 3 .
Figure 3. σA K vs. AK distribution for the extinction map of Figure 2. The individual mean extinctions per pixel are shown in black dots.The blue dots with errorbars denote the median values and standard deviations, respectively.

Figure 4 .Figure 5 .
Figure 4.A SED example for a single spatial pixel within the Ophiuchus cloud using Herschel observations.The black solid line indicates the best-fit model with β = 2 as described in Section 3.2.The derived optical depth at 250 µm is τ250 = 0.0018 and dust temperature is T d = 17.5 K.The grey shaded area illustrates the 90% confidence interval in MC analysis.

Figure 6 .Figure 7 .
Figure 6.Correlation between FIR optical depth τ250 and NIR extinction AK in three dark clouds L1689, L1709 and L1712.The green dashed line presents a linear fit to the data point using the lts linefit routine.The grey points shows the outliers as determined by lts linefit, which are excluded from the linear fit.The red dot-dashed indicates the power-law correlation obtained by fitting all the data points.The best-fit parameters are displayed in the plot.

Table 1 .
Parameters Employed in the Estimation of Extinction Using the NICE and NICER Methods.Intrinsic colors from TRILEGAL simulations.Errors are the standard deviation of the colors.(b) The coefficients for converting color excess E(m λ 1 − m λ 2 ) to extinction AK by adopting the extinction law from Li & Chen (2023).

Table 2 .
Lombardi et al. 2014;Singh & Martin 2022)ons, zero points (aν )and RMS noise at 1 ′ resolution used for SED fitting.The SPIRE Level 3 data were already zero-point corrected, however the PACS level 2.5 data have not undergone zero-point corrections.Therefore, before fitting the Spectral Energy Distribution (SED), it is imperative to correct the zero points of Herschel PACS data.Following the methodology commonly employed in prior studies (e.g.Lombardi et al. 2014;Singh & Martin 2022), we use the Planck data products derived from the Planck 2013 all-sky model of thermal dust emission (Planck Collaboration et al. 2014a) 5 .These data provide three maps of dust temperature, (T d ), dust emissivity index (β) and optical depth at 353 GHz (τ 353 ), all of which were produced from an all-sky SED fit by a modified black-body model.The resolution of the T d and τ 353 maps is 5 arcmin, while the resolution of β map is 30 arcmin.We then compute the fluxes at Herschel bands by convolving the Planck model SED with the Herschel filters at a resolution of 5 arcmin, denoted as S Planck ν. Subsequently, we perform a linear fit to the S Planck

Table 3 .
Basic Properties of the 16 Sub-regions and Measured Dust Opacities.
Yuping Tang for helpful discussions.This work is supported by the National Key R&D program of China (2022YFA1603102), the National Natural Science Foundation of China (12133002).H.Z. is funded by the China Postdoctoral Science Foundation (No. 2022M723373) and the Jiangsu Funding Program for Excellent Postdoctoral Talent.X.C. thanks to Guangdong Province Universities and Colleges Pearl River Scholar Funded Scheme (2019).