© ESO, 2017

1. Introduction

The light of stars often shows a degree of polarization correlated to interstellar extinction, up to a degree of a few percent per A_V magnitude. This phenomenon has long been recognized as the effect of non-spherical cosmic dust grains that are aligned with interstellar magnetic field lines (Hall 1949; Hiltner 1949). Dust extinguishes starlight, and the component of the electric field parallel to the longer axis of a grain is more extinguished than the orthogonal field component. Furthermore, interstellar dust grains align their shorter axes with the interstellar magnetic field, so that they are not generally randomly oriented (e.g., Andersson et al. 2015, and references therein). The overall result is that the dusty magnetized interstellar medium (ISM) polarizes the starlight that was not originally polarized. The polarization fraction p and the polarization angle ψ of starlight therefore provide information on both interstellar dust and the Galactic magnetic field, or at least the component of the field that is parallel to the plane of the sky.

It is mainly the large grains that are aligned (e.g. Kim & Martin 1995) and, in typical ISM conditions, their thermal emission mainly falls in the far-infrared (FIR) and submillimeter (submm) range. This emission is also polarized, since emission is more efficient for the electric field component parallel to the longer axis, and it is an important complement to observations of polarized extinction in the optical and near-infrared (NIR; e.g., Hildebrand 1988; Planck Collaboration Int. XIX 2015). It should be noted that since radiation polarized parallel to the longer axis of a grain is least intense in extinction and most intense in emission, we expect ψ in the submm to be orthogonal to ψ in the optical.

The main factors that determine the polarization fraction p are the optical properties of the dust, the alignment efficiency, and the orientation of the magnetic field lines. Polarization can be expressed as (Lee & Draine 1985)

$\mathit{p}\mathrm{=}{\mathit{p}}_{\mathrm{0}}\hspace{0.17em}\mathit{R}\hspace{0.17em}{\cos }^{\mathrm{2}}\mathit{\gamma ,}$where p₀ is the highest possible polarization given the dust properties, the parameter R, which is between 0 and 1, accounts for the effects of imperfect alignment¹, and γ is the angle between the magnetic field lines and the plane of the sky.

The wavelength dependence of polarization in extinction usually follows the so-called Serkowski curve (Serkowski et al. 1975),

$\mathit{p}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}{\mathit{p}}_{\max }\mathrm{·}\exp \mathrm{(}\mathrm{-}\mathit{K}\mathrm{·}\ln \mathrm{(}\mathit{\lambda }\mathit{/}{\mathit{\lambda }}_{\max }{\mathrm{)}}^{\mathrm{2}}\mathrm{)}\mathit{,}$where the polarization has a maximum p_max at a wavelength λ_max, usually falling in the visible; the value of the parameter K, tied to the inverse of the FWHM, is usually around unity. Since the polarization efficiency of a grain peaks at ~2π times its size (Kim & Martin 1995), λ_max traces the typical size of aligned grains: variations in λ_max between lines of sight may indicate a change in grain size distribution, in the dependence of alignment on size, or both (e.g. Andersson & Potter 2007, hereafter AP07). The problem is further complicated because λ_max also shows some dependence on the magnetic field angle γ (Voshchinnikov et al. 2016).

Since polarization depends on many factors at once, interpreting it is a degenerate and difficult problem. This is especially true of the dense and complex environments that are molecular clouds, where magnetic field orientation, grain alignment, and dust properties are expected to change on small scales. Despite this, studies on dust polarization are often focused on constraining the grain alignment efficiency (e.g. AP07) or the structure of the magnetic field (e.g. Planck Collaboration Int. XX 2015) without accounting for other variables. One way of confronting the problem is to construct a cloud model that includes dust evolution, magnetic field structure and grain alignment, but such models are complex and very demanding computationally. The judicious combination of observational data can also provide interesting insights on dust physics, while requiring far lighter calculations.

One example of such labor-saving combinations is the complementary use of extinction data and polarized dust thermal emission. This last has the dimensions of an intensity and it is usually observed in the FIR and submm. All-sky surveys in the submm like Planck (Planck Collaboration Int. XIX 2015) are opening up new possibilities for this type of multiwavelength analysis. The idea that extinction and emission combined can be more informative than either of them alone is explored in Planck Collaboration Int. XXI (2015), which examines the P₃₅₃/p_V ratio² in the diffuse ISM, P₃₅₃ being the polarized intensity in emission in the Planck 353 GHz channel, and p_V being the starlight polarization degree in the V band. This ratio is measured in MJy sr^-1, and since P₃₅₃ and p_V have (at first approximation) the same dependence on alignment and γ, it should provide strong constraints on the properties of aligned dust grains. Among the new results made possible by Planck is the determination of average P₃₅₃/p_V in the diffuse ISM: 5.4 MJy sr^-1, which is about 2.5 times higher than predicted by previous dust models (Planck Collaboration Int. XXI 2015).

The present paper aims to extend the study of P₃₅₃/p_V to denser environments, namely translucent lines of sight in molecular clouds, using an updated method and a dust model optimized for the high P₃₅₃/p_V found by Planck (see Guillet et al. 2016). The paper is organized as follows: Sect. 2 introduces the observational data we used, mostly from translucent lines of sight, and the selection that had to be made so that the comparison of extinction and emission would be meaningful. Section 3 presents the dust model we used, which was especially constructed to reproduce the P₃₅₃/p_V ratio and classic dust observables. The dust model was created for the diffuse ISM, while we study translucent clouds where dust evolution may be taking place: the modification that is needed for the model to fit the data suggest what the nature of dust evolution might be be in these areas. The model results are compared to observations in Sect. 4, and the meaning of this comparison is discussed in Sect. 5. Finally, Sect. 6 contains the conclusions and perspectives.

2. Data

Our work combines measures of starlight polarization in the near-ultraviolet (NUV) to NIR range that we recovered from the literature, and Planck measurements of total and polarized dust emission at 353 GHz for the same lines of sight. We use a total of 132 objects, which are reduced to 70 after selection (Sect. 2.3).

2.1. Extinction in the NUV to NIR

Most of our data points are from AP07, who include data on several clouds from the literature (Chamaeleon, Musca, Ophiuchus, R Coronae Australis (RCrA), and Taurus, Covino et al. 1997; Arnal et al. 1993; Vrba et al. 1993; Whittet et al. 1992, 2001, see Table 1) in a polarimetric study of the Coalsack nebula.The AP07 study is particularly useful for this purpose because the authors did not employ the Serkowski fit results from the literature, but they used the photometric and polarimetric data therein to conduct their own fits, thus minimizing the systematic effects from different fitting procedures. To increase the statistics, we also included data from Martin et al. (1992), who provide more lines of sight in Ophiuchus, and Anderson et al. (1996), who provide lines of sight that are for the most part³ not associated with the aforementioned clouds.

From the literature we obtained the Serkowski polarization parameters p_max and λ_max for all stars, as well as the polarization angle ψ_ext. The values of K are not calculated consistently in the literature: AP07 often use K = 1.15 and only fit K as a free parameter if it constitutes a statistically significant improvement, while Anderson et al. (1996) impose that K be a linear function of λ_max. For this reason, we chose not to include K in our work. The AP07 data retrieved from the Vizier online database⁴ do not include the uncertainties on ψ_ext for Ophiuchus, therefore we complemented the data with the original article (Vrba et al. 1993): we used the average on ψ in the various bands as the value of ψ_ext and their standard deviation as the uncertainty. We excluded stars with standard deviations greater than 7°, which we interpreted as stars where angles in different bands are not compatible. We also obtained the value of the V-band polarization p_V for most of the AP07 stars by adopting the data directly from the references (Covino et al. 1997; Arnal et al. 1993; Vrba et al. 1993; Whittet et al. 1992, 2001). The data from Anderson et al. (1996) are in the form of (polarized) spectra rather than multiband photometry; for their lines of sight we took the polarization at λ = 545 nm as the value for p_V. We did not use p_V for the stars in Martin et al. (1992), who fit data from multiple sources and thus provide non-unique values for each band. Finally, we obtained the extinction parameters for most of the stars: A_V and E(B−V) for the AP07 stars, E(B−V) for the Martin et al. (1992) and Anderson et al. (1996) stars.

2.2. Emission: Planck and IRAS submm maps

Our submm data consist of the Planck 353 GHz (850 μm) maps for the I, Q, and U Stokes parameters, from which the polarized intensity P₃₅₃ and angle ψ₃₅₃ were obtained. We did not use any other frequencies because their signal-to-noise ratio (S/N) was lower. For selection purposes we also used the all-sky submm dust opacity maps created by Marc-Antoine Miville-Deschênes using Planck and IRAS data (Planck Collaboration XI 2014).

We used the second Planck public data release⁵, which consists of HEALPix all-sky maps of ten quantities: the Stokes parameters I, Q, and U, the number of hits, the variances of the Stokes parameters II, QQ, and UU, and the covariances IQ, IU, and QU. Maps are in NESTED ordering and Galactic coordinates; they have a pixelization N_side = 2048 for a total of 12·2048² = 50 331 648 pixels with 1′.7 side lengths, so that the beam of the instrument (FWHM ~ 5′) is well sampled. The maps are in units of K_CMB and are converted into MJy/sr with a conversion factor of 287.45 at 353 GHz (Planck Collaboration IX 2014). To obtain the value of I, Q, and U at the position of each star and increase the S/N, we employed the same technique as Planck Collaboration Int. XXI (2015): we averaged the values for the Stokes parameters on a Gaussian point-spread function (PSF) centered on the star coordinates and with a FWHM of 5′, bringing the effective resolution to ~7′. In the case of Q and U, since we are working on a flat map recovered from a spherical map, we had to account for the fact that the direction of north changes from pixel to pixel. We did this by rotating the doublet (Q, U) until it was on the equator in the local reference frame.

With the value of the submillimeter Q and U for all the lines of sight, we calculated the polarized intensity in emission, ${\mathit{P}}_{\mathrm{353}}\mathrm{=}\sqrt{{\mathit{Q}}^{\mathrm{2}}\mathrm{+}{\mathit{U}}^{\mathrm{2}}}$, and the polarization angle ${\mathit{\psi }}_{\mathrm{353}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\arctan \mathrm{\left(}\mathit{U,Q}\mathrm{\right)}$ using the HEALPix angle convention (where the relative signs of Q and U are inverted with respect with the IAU convention). Because it is a quadratic function of measures with finite noise, P₃₅₃ has a positive bias, and it was debiased with the conventional formula (Wardle & Kronberg 1974): ${\mathit{P}}_{\mathrm{deb}}\mathrm{=}\sqrt{{\mathit{P}}_{\mathrm{bias}}^{\mathrm{2}}\mathrm{-}{\mathit{\sigma }}_{\mathit{P}}^{\mathrm{2}}}$. We did not apply the more recent debiasing method (e.g. Plaszczynski et al. 2014; Montier et al. 2015) because after the smoothing, the environments we studied have a high S/N, and therefore a low bias. There was no need to apply corrections for the cosmic microwave and infrared backgrounds, which are negligible at this wavelength and for our dataset.

For each star we also calculated the polarization angle dispersion function S (Planck Collaboration Int. XIX 2015), which is a tracer of disorder in polarization and therefore in the magnetic field orientation. The (I,Q,U) triplet was again smoothed to increase S/N, using a Gaussian PSF with 5’ FWHM and bringing the maps to a ~7′ resolution. The maps thus obtained are oversampled (four pixels per beam), so we also degrade the pixelization of the Q and U maps to N_side = 1024 to approach the Nyquist criterion. The dispersion function S was then computed for the pixel containing the star, with a lag δ = 5′.

2.3. Selection

Fig. 1 Illustration of line-of-sight and beam effects when comparing extinction and emission, from Planck Collaboration Int. XXI (2015). The measured extinction is E(B−V), the extinction obtained from conversion of submillimeter emission is E(B−V)S.

Since this paper compares different phenomena (dust extinction and thermal emission) that are observed at different wavelengths (NUV to NIR vs. submm), we need to ensure that the comparison is meaningful (see e.g. Planck Collaboration Int. XXI 2015). This means that first, we must be sure that we are observing the same type of grain at the two wavelengths and that second, the two wavelengths must probe the same volumes of ISM.

The first condition is met, in first approximation, as a consequence of dust physics: only large grains contribute to polarization, and this is true both in extinction and in emission. Large grains are also the main contributor to submm emission as well as visual and IR extinction (small grains are important in the UV). It should be kept in mind, however, that polarization and overall extinction do not necessarily trace the same grains: a population of large grains that are spherical or unaligned would contribute to extinction and emission, but not to polarization.

The second condition is not trivial. Extinction measured on background stars, as is the case of the data described in Sect. 2.1, only probes the matter in front of the star itself; emission has no such limitation, especially in the submm, where the ISM is optically thin (see Fig. 1). In presence of a background to the star, the total intensity I measures systematically more dust than what is observed in the optical. The effect of background on polarization is more complex, and we detail it in Sect. 2.4. To compare extinction and emission, we therefore need to discard those lines of sight that have significant dust emission from behind the star. As shown in Planck Collaboration Int. XXI (2015), this selection can be based on three criteria:

• Galactic latitude: all stars close to the Galactic planeare very likely to have significant background; wetherefore only keep stars stars with Galactic latitude| b | ≥ 2°. This forces us to exclude some well-studied clouds, such as the Coalsack nebula (AP07), which is located on the Galactic plane.

• Polarization angle: dust polarization in extinction should be orthogonal to that in emission. In lines of sight where the angles ψ_ext and ψ₃₅₃ are not orthogonal, extinction and emission do not come from the same dust; we exclude these lines of sight with a tolerance of 3σ or 10°, whichever is smaller. For stars in Anderson et al. (1996), whose ψ_ext are given without uncertainties, we assume σ(ψ_ext) = 0. Since Planck angle uncertainties are usually larger than V-band angle uncertainties, this should not make a large difference.

• Column density: the dust submm optical depth τ₃₅₃ can be converted into an expected A_V or E(B−V) and compared to the actual extinction measured; lines of sight where the τ₃₅₃-derived extinction shows an excess have significant background. For the conversion we used the empirical factor E(B−V) /τ₃₅₃ = 1.49 × 10⁴ that was obtained by Planck Collaboration XI (2014) for the diffuse ISM. In molecular clouds, however, τ₃₅₃/E(B−V) is known to increase by a factor ~2–3 compared to the diffuse ISM (e.g. Stepnik et al. 2003; Ysard et al. 2013; Planck Collaboration XI 2014). We therefore decided to relax this condition and we kept all lines of sight where the τ₃₅₃-derived E(B−V) is lower than three times the measured value.

While none of these selection procedures intrinsically exclude all of the contaminated sightlines, combining them and using them together with the selection made by AP07 gives more reliable results.

We improved the data quality by another selection. We only kept those stars for which we had an S/N greater than 3 in P₃₅₃ and greater than 5 in λ_max. For a few of the stars in AP07 the quality of the Serkowski fit was low and the Serkowski parameters were not an adequate representation of the polarization curve; we recovered the observational data from Covino et al. (1997), Arnal et al. (1993), Vrba et al. (1993), Whittet et al. (1992, 2001) and excluded those stars that, at a visual inspection, do not follow Serkowski. Finally, we excluded those stars that according to Anderson et al. (1996) are likely to have intrinsic polarization. The combined selection left us with the values of λ_max, p_max, E(B−V), I₃₅₃, and P₃₅₃ for 70 lines of sight, 56 of which also have information on A_V and p_V.

2.4. Line-of-sight and beam depolarization

The magnetic field in the ISM has a non-negligible disordered (or “meandering”) component that introduces a confounding variable called “depolarization”. When there is confusion between polarized sources with different orientation angles, it is possible for the orthogonal components of the polarizations to cancel each other out, so that the overall polarization that is observed may be lower than that of each source taken separately. Depolarization may occur if the interstellar magnetic field changes orientation along a line of sight (line-of-sight depolarization), or if an instrument has a finite observational beam and the magnetic field changes orientation on scales smaller than said beam (beam depolarization). In most polarization studies, the two effects are subsumed under the name of “beam depolarization” or simply “depolarization”. However, since the two types of depolarization have different effects on the extinction and the emission, we treat them separately here.

The line-of-sight depolarization, at first approximation, has the same effect on extinction and emission if extinction and emission probe the same ISM. Complications arise if there is significant emission from the background to the star: if the magnetic field orientation is very different in the foreground and in the background, depolarization in emission may be very different from that in extinction. This would give unreliable measurements of the ratio P₃₅₃/p_V, for example. The selection described in Sect. 2.3, if effective, should ensure that the line-of-sight depolarization affects extinction and emission in the same way. We remark that having a uniform magnetic field orientation on the line of sight is not equivalent to having no background emission: since P₃₅₃ is additive, in this case we would observe an excess of polarization in emission as compared to extinction, and overestimate P₃₅₃/p_V.

The beam depolarization affects observations that have finite beam size. This is usually the case at FIR and submm wavelengths: the Planck beams measure 5′ or more. Extinction observations on stars, on the other hand, are point-like and suffer no beam effects, so that beam depolarization affects only polarized emission. Unlike line-of-sight depolarization, beam effects are unaffected by our selection criteria. However, the amount of beam depolarization can be estimated from observational data, such as the function S (Sect. 2.2), which measures field disorder.

2.5. Final observables

Fig. 2 Observational data to reproduce. Different colors indicate different clouds, different symbols indicate different references (circles: AP07; triangles: Anderson et al. 1996, squares: Martin et al. 1992). Top: the “classical” pV/τ vs. λmax plot. This panel does not contain the stars without AV and pV measurements. Center: comparison of λmax in extinction with the polarization P353/I353 in emission. Bottom: the ratio of polarization P353/pmax as a function of λmax.

Our observational data that are to be compared to a dust model are plotted in Fig. 2. The top panel shows the normalized polarization in extinction, p_V/τ_V, as a function of λ_max. The two quantities have a clear negative correlation; we also see that the values for polarization are very widely scattered, and their upper limit marks an envelope, as is typical of polarized observations, whose shape may be partly determined by line-of-sight depolarization (e.g. Planck Collaboration Int. XX 2015). A very similar behavior can be seen in the submm polarization fraction P₃₅₃/I₃₅₃ as a function of λ_max, shown in the central panel. The bottom panel of Fig. 2 shows a different observable: P₃₅₃/p_max. This quantity, like the ratio P₃₅₃/p_V used in Planck Collaboration Int. XXI (2015), is meant to trace the optical properties of grains by normalizing out the effects of alignment and magnetic field orientation that affect both emission and extinction. We decided to use p_max in the construction of this ratio, as opposed to p_V, to avoid introducing spurious correlations: many of our stars have high values of λ_max (>0.6 μm), and p_V and λ_max are going to be negatively correlated in that range. The bottom panel of Fig. 2 shows that even if the scatter in P₃₅₃/p_max is quite large, it is still small compared to the scatter observed in p_V/τ_V and P₃₅₃/I₃₅₃; the dependence of P₃₅₃/p_max on λ_max is also much less pronounced. This is consistent with our expectations that this quantity be nearly independent of alignment and magnetic field orientation.

3. Model: DUSTEM with polarization

The dust model we use should ideally have the three following characteristics: it should predict polarization in both extinction and emission, it should be compatible with the latest results from the Planck mission (especially the P₃₅₃/p_max ratio, which is underpredicted by pre-Planck models), and it should allow to modify the dust properties to simulate dust evolution. Unfortunately, while models that allow detailed dust evolution exist, to our knowledge, they either do not predict polarization (e.g. Jones et al. 2013) or are calibrated on extinction alone and cannot be expected to reproduce the correct P₃₅₃/p_max ratio (e.g. Hirashita & Voshchinnikov 2014). We decided to use a model optimized for fitting the latest Planck data instead, at the cost of a simplified treatment of dust evolution where only grain size is accounted for.

We adopted the dust model recently developed by Guillet et al. (2016), called “Model A”, which is a modified version of the model by Compiègne et al. (2011). The computation is performed with the DustEM Fortran numerical tool⁶ and its IDL wrapper⁷. The model populations and size distributions are chosen to minimize the number of free parameters; the parameters themselves are calculated by fitting the observables typical of the low-latitude diffuse ISM (| b | < 30°): the extinction curve, the polarization in extinction up to 4 μm, and the spectral energy distribution (SED) in emission, polarized and non, updated with Planck results. The model therefore reproduces the average observations for the diffuse ISM, including ratios such as I₃₅₃/A_V and P₃₅₃/p_V.

The model includes three grain types (see Table 2 and Fig. 3): a population of neutral polycyclic aromatic hydrocarbons (PAHs) with a lognormal size distribution, plus two populations of large grains (amorphous carbon and silicates) that are distributed as power laws: dn/ da = a^α. We are mainly interested in observables where the large grain contribution is dominant, therefore the model, unlike the model of Compiègne et al. (2011), has no separate population for very small carbonaceous grains: the very small grains are included in the amorphous carbon population, which is why the power law for carbons is weighted toward small sizes. Large grains are prolate spheroids with an axial ratio of 3 (oblate grains of the same axial ratio cannot reproduce the high P₃₅₃/I₃₅₃ observed by Planck: Guillet et al. 2016). The neutral PAHs and the amorphous carbon grains have the same compositions as their counterparts in Compiègne et al. (2011); the silicate grains have the same composition as in Weingartner & Draine (2001), with added porosity: 20% of their volume consists of vacuum inclusions. The porosity of the silicate grains is essential in increasing their P₃₅₃/p_V ratio to the value observed by Planck.

Fig. 3 Size distribution for the model we used (model A from Guillet et al. 2016).

Fig. 4 Top: grain alignment function for different values of aalig. Bottom: corresponding polarization curves in extinction, normalized to NH = 1021 cm-2. Dots indicate the value of aalig (top) and that of λmax from free-K fits (bottom).

In the model silicate grains are aligned according to the phenomenological alignment function provided by DustEM: $\mathit{f}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{f}}_{\max }\hspace{0.17em}\left(\mathrm{1}\mathrm{+}\tanh \left(\frac{\ln \mathrm{(}\mathit{a}\mathit{/}{\mathit{a}}_{\mathrm{alig}}\mathrm{)}}{{\mathit{p}}_{\mathrm{stiff}}}\right)\right)\mathrm{,}$(1)where f(a) is the fraction of grains aligned as a function of the equivalent radius⁸a, f_max is the maximum alignment efficiency, a_alig is the size threshold for grain alignment and the parameter p_stiff denotes the width of the transition. This alignment function is designed to increase monotonically with size, since small grains are generally unaligned. This parametric function is not designed to test any particular alignment process, and while its shape resembles the typical result of the radiative torque model (Lazarian & Hoang 2007), it is also compatible with magnetic alignment for grains with superparamagnetic inclusions (Mathis 1986; Voshchinnikov et al. 2016). The variation of the alignment function for varying a_alig and its effects on the polarization curve in extinction, are shown in Fig. 4.

In addition to causing polarization, grain alignment affects dust extinction and emission as well. The resulting correction is very small and it is generally ignored in dust models; nonetheless, DustEM provides the option of including alignment effects in extinction and emission. We here chose to ignore these effects to allow a more direct comparison to the results of Guillet et al. (2016), where minor alignment effects were ignored as well.

3.1. Model results: fitting

DustEM provides the full extinction curve and dust emission for the model, including polarization; most of the parameters introduced in Sect. 2.5 have to be obtained from a fit to the DustEM output.

We took the extinction curve interpolated at 550 nm to be A_V of the model, and the thermal emission at 353 GHz, plus a color correction to account for the spectral response of the corresponding Planck band, as the value I₃₅₃ of the SED. The same interpolation, operated on the polarized extinction and emission, gives us the model prediction for p_V and P₃₅₃. The p_max, λ_max, and K of the model were calculated by interpolating the model polarized extinction at the effective central wavelengths of the UBVRIJH photometric bands (Table 3) and fitting a Serkowski function to the synthetic observations thus obtained, keeping K as a free parameter.

We also fit the model SED with a modified blackbody: I_λ = B_λ(T)·τ₀·(λ/λ₀)^{− β}. The fit was performed on the model emission interpolated and color-corrected at the wavelengths of the Planck HFI bands (350, 550, 850, 1380, and 2100 μm) and the IRAS 100 μm band. Since emission at those wavelengths is dominated by large grains, integrating the modified blackbody over wavelength provides the radiance ℛ, or emitted power, of the large grain populations (see Planck Collaboration XI 2014).

3.2. Model variations

The model so far described was developed to fit the average observables in the diffuse ISM. Inside molecular clouds, however, evolution alters the properties of dust considerably: the main alteration is grain growth through accretion and coagulation (Boulanger et al. 1990; Stepnik et al. 2003; Köhler et al. 2012). Translucent lines of sight, such as those studied in this paper, typically probe regions at the onset of such an evolution: for instance, Stepnik et al. (2003) and Ysard et al. (2013) found that coagulation takes place where A_V is greater than 2 or 3. This raises the question of whether our sample can be explained by a dust model designed for the diffuse ISM. To answer this question, we studied how the model output is affected by grain alignment efficiency, magnetic field orientation, and grain size. By comparing these results to the observations we can determine whether the data can be explained by the variation of alignment and field orientation alone, using the same dust as in the diffuse ISM, or if dust growth is necessary, and to what extent. The modifications to our model are purely phenomenological and are not based on simulations of grain growth, dust alignment or magnetic field structure; however, they are useful for estimating the variations that physical models will have to reproduce.

3.2.1. Variations in alignment function

Loss of grain alignment inside molecular clouds is sometimes invoked to explain the decrease in polarization efficiency observed at high A_V (e.g. Andersson et al. 2015, and references therein). This weakening of polarization is also in qualitative agreement with some alignment theories, e.g., the radiative torque model, which predicts that only the largest grains are aligned inside molecular clouds.

The alignment efficiency in DustEM is a function of the three parameters a_alig, p_stiff, and f_max, as shown in Eq. (1). We simulated different alignment efficiencies by running the model with different values of a_alig, keeping the same p_stiff for simplicity. We did not change f_max (equal to 1 in our model) because it has the exact same effects on p_V and P₃₅₃, and it could have no effect on P₃₅₃/p_V; the parameter a_alig, on the other hand, affects p_V and P₃₅₃ similarly but not identically, since the polarization cross-section has a different size-dependent behavior in the visible (where scattering is dominant) and in the submm.

Figure 4 shows the change in f(a) as a function of a_alig and the effect that this has on the polarization curve; higher a_alig corresponds to lower p_max because fewer grains are aligned, and to higher λ_max because the average aligned grain is larger.

3.2.2. Variations in magnetic field orientation

An ordered magnetic field forming an angle γ with the plane of the sky introduces a factor cos²γ in the polarized intensity, as shown in Sect. 1. Dust polarization models would greatly benefit from measures of the angle γ; unfortunately, dust only traces the magnetic field component parallel to the plane of the sky, which means that this information is not usually available. The angle γ is therefore another variable parameter in our model: we ran the model for γ = 0°, 30°, 45°, and 60°.

The Galactic magnetic field is the sum of an ordered component and a disordered, or meandering, component: this latter component causes the phenomena known as line-of-sight and beam depolarization, as explained in Sect. 2.4. Our model does not include a disordered magnetic field component and therefore it cannot predict depolarization; however, the polarization angle dispersion S (Sect. 2.2) can be used as a measure of field disorder (see e.g. Planck Collaboration Int. XX 2015, which compare the observed S with MHD simulations). Using this we were able to assess some of the effects of the disordered magnetic field, as shown in Sect. 4.1.

3.2.3. Variations in grain size distribution

Gas accretion on grain surfaces (e.g. Jones et al. 2013) and formation of aggregates are known to increase grain sizes inside molecular clouds. This growth is supported by theoretical studies (Köhler et al. 2012; Hirashita & Voshchinnikov 2014) and it is consistent with observed phenomena such as the flattening of extinction curves in dense environments (e.g. Fitzpatrick 1999; Weingartner & Draine 2001) and the coreshine observed in the NIR (Pagani et al. 2010).

As we described, no dust models currently treat dust evolution realistically while reproducing the emission-to-extinction polarization ratio revealed by Planck. We therefore opted for a model that is compatible with Planck data (model A from Guillet et al. 2016) at the cost of a simplified treatment of dust evolution. Specifically, we focused on a single aspect affected by dust evolution: the size distribution of grains. In our model, the size distribution for large grains is a power law defined by three parameters: the minimum and maximum grain sizes a_min and a_max, and the power-law index α.

Most of the mass in silicates is in the large grains (see Fig. 3), meaning that the size distribution is most sensitive to a_max; furthermore, none of the observables we use are in the UV where the contribution of small grains is important. Therefore, we decided to vary the a_max of silicates between 350 nm and 1 μm (the standard value is ~500 nm) while keeping α and a_min fixed. Although we chose to fix α mainly as a matter of convenience, we note that this is consistent with the model of grain growth by Hirashita & Voshchinnikov (2014), where the slope of the size distribution does not change much during evolution, and the largest variation is in the upper size cutoff.

The case of carbon grains is different, as their distribution is weighted toward small sizes: the mass available for large grains is now also dependent on a_min and on the amount of PAHs, so that a realistic model becomes a necessity. For this reason we only varied the size distribution of silicates while leaving that of carbon grains fixed. While this choice does not give realistic results for the variation of extinction and emission with size distribution, it still allows us to predict the dust polarization, which in our model depends on silicates alone.

3.2.4. Multiparameter study: Monte Carlo simulation

The phenomena described in Sects. 3.2.1 to 3.2.3 are all expected to occur in molecular clouds, so that variations of a single model parameter at a time are not realistic, even if studying their effect can be instructive. We decided to use a Monte Carlo simulation to explore the effects of simultaneous variations of many parameters. As explained in the previous section, we can only vary the size distribution of silicates, which gives realistic results for polarization but not for unpolarized observables; therefore our Monte Carlo results can only be compared to polarization observables. The model was run one thousand times, and the values of a_alig and a_max for silicates were uniformly distributed within the ranges $\begin{array}{ccc}& & \mathrm{50}\hspace{0.17em}\mathrm{nm}\mathit{<}{\mathit{a}}_{\mathrm{alig}}\mathit{<}\mathrm{300}\hspace{0.17em}\mathrm{nm}\\ & & \mathrm{350}\hspace{0.17em}\mathrm{nm}\mathit{<}{\mathit{a}}_{\max }\mathit{<}\mathrm{1000}\hspace{0.17em}\mathrm{nm}\mathit{.}\end{array}$The operation was repeated for four values of γ (0°, 30°, 45°, and 60°), bringing the simulations to a total of 4000.

We found that not all combinations of a_alig, a_max and γ in the ranges chosen are realistic: the synthetic observables obtained for some such combinations have values that are never observed. We set out to find a realistic range of parameters by imposing that model results have a range as close as possible to that of actual observations (e.g., that 0.35 <λ_max< 0.8 μm and 0.5 <K< 1.5 as per Voshchinnikov & Hirashita 2014). Restricting a_alig and a_max to the following ranges eliminates most of the unrealistic values for λ_max and K: $\begin{array}{ccc}& & \mathrm{75}\hspace{0.17em}\mathrm{nm}\mathit{<}{\mathit{a}}_{\mathrm{alig}}\mathit{<}\mathrm{150}\hspace{0.17em}\mathrm{nm}\\ & & \mathrm{350}\hspace{0.17em}\mathrm{nm}\mathit{<}{\mathit{a}}_{\max }\mathit{<}\mathrm{800}\hspace{0.17em}\mathrm{nm}\mathit{,}\end{array}$while the same four values for γ are kept. This selection left us with 844 Monte Carlo iterations.

4. Results

4.1. Alignment efficiency and magnetic field orientation

Fig. 5 Observational data compared with our variable-alignment model (75 ≤ aalig ≤ 150 nm). The color and symbol scheme are the same as in Fig. 2. The black square marks the position of the standard model. The top panel does not contain the stars without AV and pV measurements.

The effects of dust alignment efficiency and magnetic field orientation are shown in Fig. 5, which compares the model results with observational data. Dots represent the observed values of p_V/τ_V (top), P₃₅₃/I₃₅₃ (middle), and P₃₅₃/p_max (bottom) as a function of λ_max. Curves represent the model; within each curve, a_alig varies between 75 and 150 nm, and the four curves correspond to the four values of γ, 0°,30°,45°, and 60°.

The combination of variable alignment and magnetic field orientation can reproduce most of the observations in the case of p_V/τ_V and P₃₅₃/I₃₅₃, both the general trends and the dispersion. The curve for γ = 0° coincides roughly with the highly polarized, low-λ_max lines of sight in our sample. For higher values of γ the polarization decreases, but the dispersion in λ_max caused by the variation of a_alig increases, pushing the maximum λ_max to higher values: as a result, the model predicts that weakly polarized lines of sight that can have either small or large λ_max, which is indeed the trend found in the observational data. The relation between γ and λ_max described by Voshchinnikov et al. (2016) is evident in Fig. 5: although λ_max is mainly affected by the alignment size threshold a_alig, the model curves with a higher γ are clearly shifted to higher values of λ_max. The figure, however, reveals something more: the strength of the γ-λ_max relation itself increases with a_alig. The leftmost tips of the model curves, corresponding to a_alig = 75 nm, all have very similar λ_max. In contrast, the rightmost tips, corresponding to a_alig = 150 nm, show wide differences in λ_max, comparable to the differences that are due to a_alig itself. It should be noted that our model assumes a uniform magnetic field, and therefore it does not include line-of-sight or beam depolarization. If these effects were important, it might mean that the ordered component of the field is closer to the plane of the sky than our model predicts.

Figure 5 also shows that alignment and magnetic field orientation have very little effect on the model results for P₃₅₃/p_max, as was indeed expected: the different curves are close to each other. In fact, the dispersion observed in the observed P₃₅₃/p_max is much larger than predicted by the model, which suggests that variations in the dust optical properties occur in translucent lines of sight. Again, one possible confounding factor in this interpretation is the depolarization caused by meandering of the magnetic field: we will now attempt to assess its effects.

Fig. 6 Correlation between the observables we studied and the magnetic field meandering as measured by S. Top: polarization fraction P353/I353 as a function of λmax. Bottom: polarization ratio P353/pmax as a function of λmax. The color and symbol scheme are the same as in Fig. 2.

While our model assumes a uniform magnetic field, the disorder in the field lines can be estimated from the angle dispersion function S (Sect. 2.2). The top panel of Fig. 6 shows an anticorrelation between S and the polarization fraction in emission, a well-known effect that is commonly attributed to line-of-sight depolarization caused by meandering of the field (e.g. Planck Collaboration Int. XX 2015). Beam depolarization is also a possible cause, but we can see that beam effects appear negligible in our sample: the bottom panel of Fig. 6 shows that there is no clear influence of field meandering on the polarization ratio P₃₅₃/p_max. Beam depolarization, if present, should affect emission but not extinction, introducing an anticorrelation between S and P₃₅₃/p_max. The fact that we see no such correlation suggests that beam depolarization is negligible in our sample: this supports the idea that we are probing dust in relatively homogeneous regions, which is an important assumption in extinction/emission comparison. The effect of line-of-sight depolarization on the polarization fraction is unfortunately impossible to assess without more advanced modeling, but even if present, it should have little effect on P₃₅₃/p_max following our selection (Sects. 2.3 and 2.4).

4.2. Grain size distribution

Fig. 7 Observational polarization ratio P353/pmax compared with a variable size distribution model. The solid line shows the evolution of the model for 350 ≤ amax ≤ 103 nm (left to right). The dashed line shows the effect of variable power-law index: −3.82 ≤ α ≤ −2.82 (left to right). The black square marks the position of the standard model. The color and symbol scheme are the same as in Fig. 2.

As explained in Sect. 3.2.3, we ignored variations of size distribution in carbonaceous grains, which affect I₃₅₃ and A_V, but not P₃₅₃ and p_max. Because of this, our results on the effects of grain growth on observations are unrealistic for P₃₅₃/I₃₅₃ and p_V/τ_V, but realistic for P₃₅₃/p_max, and we show only this last observable.

Figure 7 compares the observations and the model results for 350 ≤ a_max ≤ 10³ nm (solid black line). The effect of a variable a_max on λ_max is small compared to the effect of grain alignment; on the other hand, a_max has a strong effect on P₃₅₃/p_max and is a plausible contributor to the large dispersion observed in this quantity. For comparison, the picture also shows the model results for fixed a_max and α varying of ±0.5 around the standard value, that is, −3.82 ≤ α ≤ −2.82 (dashed line). The variations in α have a modest effect on λ_max comparable to that of a_max; however, unlike a_max, the parameter α has nearly no influence on the P₃₅₃/p_max ratio.

4.3. Multiparameter Monte Carlo

Fig. 8 Observational polarization ratio P353/pmax compared with Monte Carlo iterations of the model (grey). The black square marks the position of the standard model. The color and symbol scheme are the same as in Fig. 2.

We mentioned above that, realistically, all the model parameters we use will simultaneously vary inside a molecular cloud, and we chose to represent this with a Monte Carlo simulation (Sect. 3.2.4). As in the previous section, our modeling does not account for the size distribution of carbon grains, and we only show polarization results.

Figure 8 compares the Monte Carlo results to the observational data. The Monte Carlo simulation was run for γ = 0°, 30°, 45°, and 60°; the other model parameters are uniformly distributed in the regions 75 nm ≤ a_alig ≤ 150 nm and 350 ≤ a_max ≤ 800 nm. The combined variation of dust alignment, field orientation, and grain size allows the model to reproduce most of the observations. However, some lines of sight, spanning the full range of observed λ_max, have a lower P₃₅₃/p_max than the model can reproduce. This may be the result of variations in the dust polarization properties that are due to factors other than size (e.g., grain shape, porosity, or chemical composition), which are not considered in our model.

5. Discussion

Fig. 9 P353/pmax ratio as a function of the modified blackbody temperature (top panel) and as a function of the normalized radiance ℛ /AV (bottom panel). The bottom panel does not include stars without AV and pV measurements. The grey curves show the model for varying intensity of the radiation field (0.1 ≤ G0 ≤ 3).The grey square marks the position of the standard model. The color and symbol scheme for the observarions are the same as in Fig. 2.

We have shown that the grain alignment and the orientation of the magnetic field both have important effects on our observables: a combination of variable a_alig and variable γ reproduces the general trends and the dispersions for p_V/τ_V and P₃₅₃/I₃₅₃ as a function of λ_max; in fact, we find the familiar envelope in the distribution of points (Fig. 5). There are some stars, however, mainly in the Musca cloud, with a higher polarization than the model can reproduce. This is because the model was made to reproduce a p_V/τ_V of ~3%, the usually accepted value for the diffuse ISM (e.g. Serkowski et al. 1975; Andersson et al. 2015), while the outlier stars in Fig. 5 have p_V/τ_V ~ 4%. We should point out, however, that the measures in question have relatively large error bars and most stars are within ~1σ above the envelope.

Figure 8 also shows that a combination of variable alignment, magnetic field orientation, and grain growth reproduces the general trend and most of the data scatter in the P₃₅₃/p_max vs. λ_max relation. Again, some data are beyond the model range; in this case it is the stars with low P₃₅₃/p_max. These lines of sight may belong to regions with different dust properties: while our model reproduces the high P₃₅₃/p_max ratio from Planck Collaboration Int. XXI (2015), most models give lower values. Another possibility is that the low P₃₅₃/p_max is an effect of beam depolarization, which would lower the observed value of P₃₅₃ without affecting p_max; however, we showed in Sect. 4.1 that beam depolarization does not seem important in our sample. A low P₃₅₃/p_max could also be due to low dust heating, because P₃₅₃ is proportional to dust emission: our model results are calculated for the typical interstellar radiation field intensity⁹ in the diffuse ISM, G₀ = 1, but we would expect the radiation field to be less intense in clouds. However, in the long-wavelength range of the Planck function, which is the case of Planck observations of dust, emission tends to depend linearly on dust temperature, meaning that the effect of heating on dust SED may be weak.

A proper estimation of heating effects would be no trivial task, as there is no univocal relation between the observed A_V and the extinction actually experienced by dust (see AP07 for a discussion), but a preliminary analysis is shown in Fig. 9. The gray curves in the figure represent the model results for 0.1 ≤ G₀ ≤ 3. The top panel shows no visible correlation between P₃₅₃/p_max and dust temperature, suggesting that G₀ only has modest effects (it should be kept in mind, though, that temperature is also influenced by dust evolution: e.g., Stepnik et al. 2003; Planck Collaboration Int. XVII 2014). The same conclusion seems supported by the bottom panel: this shows that P₃₅₃/p_max is not correlated with the normalized radiance ℛ /A_V. Since the radiance ℛ is the bolometric emission of large dust grains (Sect. 3.1), ℛ /A_V is a measure of the power emitted (and therefore absorbed) by large grains, averaged on the line of sight¹⁰.

6. Conclusions and perspectives

Dust polarization depends on the efficiency of grain alignment, the orientation and meandering of the magnetic field, and the optical properties of the dust itself. Many studies have attempted to explain polarization in molecular clouds in terms of one of these factors (see, e.g., AP07 for a study focused on alignment and Planck Collaboration Int. XX (2015) for one focused on magnetic field meandering); nonetheless, it is likely that all factors are at play at once.

In this paper we used dust model A from Guillet et al. (2016) and varied its alignment efficiency, grain size, and magnetic field orientation to reproduce the diverse conditions that one may find in molecular clouds. The results were compared to extinction and emission data from both bibliographic sources and the Planck survey, with particular attention to the polarization observables p_V/τ_V, P₃₅₃/I₃₅₃, and P₃₅₃/p_max as a function of λ_max.

We find that none of the model parameters employed can explain the full set of observations on its own. Monte Carlo simulations show that most of the data can be reproduced by letting a_alig vary between 75 and 150 nm, a_max vary between 350 nm and 800 nm, and γ vary between 0° and 60°. This means that any studies of polarization in molecular clouds need to take into account all these aspects to explain the full range of the data. In particular, the ratio P₃₅₃/p_max is very useful in reducing the contributions of alignment and magnetic field orientation, and highlighting variations in dust properties, and especially size distribution. Within the context of our model, the variations observed in P₃₅₃/p_max can be partly explained by varying the maximum grain size a_max, while the power-law index α has little effect.

Nonetheless, some of the observations fall beyond the range of the model results, most notably some lines of sight with very low values of P₃₅₃/p_max that are found over the full range of observed λ_max. This probably indicates variations in dust properties other than size distribution, such as shape, structure, or chemical composition. Dust evolution that is not related to size may also influence our estimate of magnetic field orientation: for instance, lines of sight with a low P₃₅₃/I₃₅₃, which our model attributes to magnetic field lines nearly orthogonal to the plane of the sky (large γ), may be instead explained by dust with a low polarization cross-section. The lines of sight with low P₃₅₃/p_max can in principle be explained without a need for dust evolution: a dim radiation field due to extinction, or the beam depolarization due to a disordered magnetic field, would lower the value of polarization in emission without affecting extinction. However the analysis we conducted in the previous sections suggests that neither dust heating nor beam depolarization have much effect on our sample.

It is possible that the width of the magnetic field orientation range found by our model (0° to 60°) is overestimated. In addition to the aforementioned influence of dust properties, field meandering, which is absent in our model, could introduce line-of-sight depolarization, the effects of which are degenerate with increasing angle γ of (the ordered component of) the magnetic field. Including field meandering in an ISM model would make for a very interesting follow-up study, but a polarized radiation transfer code is needed for this. It would also be useful to independently determine γ in future research, for example, using MHD simulations or where available, measures of the line-of-sight magnetic field such as the Zeeman polarization.

Not all potentially interesting cases were considered in this paper: this is a first application of this technique to the study of polarization in molecular clouds. The full implications of this technique will become clearer with more detailed modeling and more observational data. A continuation of this work would benefit from extending the dataset to NIR extinction and polarization. Observations in the NIR can probe denser lines of sight than those in the visible; furthermore, increasing the wavelength of observation means that we are closer to the Rayleigh limit, so that the observables are better trackers of the overall mass of aligned grains and less dependent on the details of alignment and size distribution. This makes NIR observations an interesting complement to observations in the optical. Finally, different

types of observations could improve our constraints on the model variables: maps of molecular lines and elemental depletion could be useful in constraining grain growth processes such as accretion and coagulation.

Acknowledgments

We would like to thank S. N. Shore for his stimulating discussion and insightful remarks on a variety of subjects, from ISM physics to data analysis to writing, during the editing of this article. We are also indebted to N. V. Voshchinnikov for his helpful comments.

References

All Tables

All Figures

	Fig. 1 Illustration of line-of-sight and beam effects when comparing extinction and emission, from Planck Collaboration Int. XXI (2015). The measured extinction is E(B−V), the extinction obtained from conversion of submillimeter emission is E(B−V)S.
In the text

Fig. 2 Observational data to reproduce. Different colors indicate different clouds, different symbols indicate different references (circles: AP07; triangles: Anderson et al. 1996, squares: Martin et al. 1992). Top: the “classical” pV/τ vs. λmax plot. This panel does not contain the stars without AV and pV measurements. Center: comparison of λmax in extinction with the polarization P353/I353 in emission. Bottom: the ratio of polarization P353/pmax as a function of λmax.

In the text

	Fig. 3 Size distribution for the model we used (model A from Guillet et al. 2016).
In the text

	Fig. 4 Top: grain alignment function for different values of aalig. Bottom: corresponding polarization curves in extinction, normalized to NH = 1021 cm-2. Dots indicate the value of aalig (top) and that of λmax from free-K fits (bottom).
In the text

	Fig. 5 Observational data compared with our variable-alignment model (75 ≤ aalig ≤ 150 nm). The color and symbol scheme are the same as in Fig. 2. The black square marks the position of the standard model. The top panel does not contain the stars without AV and pV measurements.
In the text

	Fig. 6 Correlation between the observables we studied and the magnetic field meandering as measured by S. Top: polarization fraction P353/I353 as a function of λmax. Bottom: polarization ratio P353/pmax as a function of λmax. The color and symbol scheme are the same as in Fig. 2.
In the text

	Fig. 7 Observational polarization ratio P353/pmax compared with a variable size distribution model. The solid line shows the evolution of the model for 350 ≤ amax ≤ 103 nm (left to right). The dashed line shows the effect of variable power-law index: −3.82 ≤ α ≤ −2.82 (left to right). The black square marks the position of the standard model. The color and symbol scheme are the same as in Fig. 2.
In the text

	Fig. 8 Observational polarization ratio P353/pmax compared with Monte Carlo iterations of the model (grey). The black square marks the position of the standard model. The color and symbol scheme are the same as in Fig. 2.
In the text

Fig. 9 P353/pmax ratio as a function of the modified blackbody temperature (top panel) and as a function of the normalized radiance ℛ /AV (bottom panel). The bottom panel does not include stars without AV and pV measurements. The grey curves show the model for varying intensity of the radiation field (0.1 ≤ G0 ≤ 3).The grey square marks the position of the standard model. The color and symbol scheme for the observarions are the same as in Fig. 2.

In the text