© The Authors 2025

1 Introduction

Dust is ubiquitous in the interstellar medium (ISM) of our Galaxy and plays a significant role in many astrophysical processes. The thermal (gray-body) emission of interstellar dust spans many frequencies (100-1000 GHz), making dust a powerful tracer of the morphology of ISM structures as well as a diagnostic of their physical properties, such as temperature (Hildebrand 1983; André et al. 2010). At the same time, dust emission poses a challenge for cosmological experiments because it veils the cosmic microwave background signal over a wide range of frequencies (Planck Collaboration X 2016). Due to the multifaceted role of interstellar dust in astrophysics and cosmology, a detailed understanding of its nature is imperative.

The interaction of aspherical dust grains with the magnetic field permeating the ISM leads to their alignment, resulting in the polarization of light, which can be observed from ultraviolet to submillimeter wavelengths (Andersson et al. 2015). Short-wavelength (UV, optical, near-IR) polarization is induced by dichroic extinction of the starlight (which usually starts out unpolarized) by aligned dust grains. The same grains emit long-wavelength (far-IR and submillimeter) polarized thermal radiation (Whittet 2022; Hensley & Draine 2023). Dust polarization provides one of the most reliable and widely employed methods of tracing the properties of the ISM magnetic fields (e.g., Skalidis & Tassis 2021; Skalidis et al. 2021, 2022; Pattle et al. 2023; Pelgrims et al. 2023, 2024).

Dust-induced stellar polarization follows a complex empirical relation with wavelength, known as the “Serkowski relation” (Serkowski et al. 1975). It expresses the polarization fraction, P, as a function of the wavelength, λ, through the formula $$P(\lambda )=P_{max}\cdot \exp (-K\cdot {\ln }^2\frac{{\lambda }_{max}}{\lambda }),$$(1)

with the use of three free parameters: Pmax, the maximum polarization fraction observed at λ_max, and K, which quantifies the spread of the profile. The polarization angle (also referred to as electric vector position angle, EVPA, denoted here as θ) is implicitly assumed to be constant with wavelength.

Since its discovery by Serkowski et al. (1975), this relation has been extensively studied and confirmed by a number of authors (Wilking et al. 1980; Whittet et al. 1992, 2001; Bagnulo et al. 2017; Cikota et al. 2018). It is so well established that it is often dubbed the “Serkowski law.” It is even a standard practice in the literature to use it to extrapolate the polarization fraction in one particular wavelength from measurements made at another wavelength (Bijas et al. 2024; Neha et al. 2024; Ginski et al. 2024). Most importantly, the Serkowski relation with its complicated shape provides stringent constraints for dust grain models (Guillet et al. 2018; Hensley & Draine 2021; Draine 2024).

The Serkowski relation, by construction, describes the polarization induced by a single cloud. For this reason, in the bibliography, whenever the formula is applied, it is implicitly assumed (even if not stated directly) that the observed (integrated) P(λ) signal is induced by a single dominant polarization screen (cloud) along the LOS. Reality, however, is more complicated. Gaia has revealed the third dimension of the sky (depth, Gaia Collaboration 2018), and since then it has become clear that in the majority of the LOSs there is more than one cloud with a comparable contribution in the total dust column (e.g., Green et al. 2019; Edenhofer et al. 2024). The variations in the magnetic field along the LOS are imprinted on the individual stellar polarization measurements, allowing for the tomographic mapping of the ISM magnetic field (Pelgrims et al. 2024).

The impact of multiple polarizing clouds along the LOS on stellar polarization has been discussed in the literature, though not extensively. Most of the focus has been on the effect on the EVPA, with less attention given to the polarization fraction spectral profile. During the early stages of interstellar polarization (ISP) studies, when the notion of a wavelength-independent EVPA emerged from observational evidence, Treanor (1963) briefly mentioned the possibility of a nonconstant EVPA in scenarios involving multiple clouds with different magnetic field alignment and/or different grain sizes. Coyne & Gehrels (1966) attributed the observed variability in EVPA with wavelength among their targets to the potential presence of multiple clouds along the LOS. Martin (1974) was the first to construct mathematical models to describe such scenarios, involving multiple clouds with differing properties along the LOS. Some authors focused on the impact of multiple dust components in the LOS on the parameters derived from Eq. (1). For example, Codina-Landaberry & Magalhaes (1976) explored the behavior of the K parameter in a medium with a changing grain alignment, while Clarke & Al-Roubaie (1984) investigated how the presence of two clouds along the LOS would influence the parameters K and λ_max, finding no firm relationship between their models and contemporary observations.

While modern researchers occasionally discuss scenarios of multiple clouds, most of the focus remains on the EVPA (e.g., McMillan & Tapia 1977; Messinger et al. 1997; Clarke 2010; Whittet 2015; Bagnulo et al. 2017; Whittet 2022). Little attention is given to the broader consequences of such scenarios, in either the P - λ profile¹ or the derived dust parameters.

In this article, we revisit the effects of the 3D structure of the ISM on the wavelength dependence of ISP, accounting for recent developments in dust distribution modeling and utilizing an expanded archive of polarimetric measurements. We use analytical models to show how both the P - λ and θ - λ profiles are expected to behave in scenarios with multiple clouds along the LOS. For the first time, we present compelling observational evidence to demonstrate these effects and disentangle the polarizing components within the LOS. Ultimately, we demonstrate that the LOS variations in dust and magnetic field properties affect the Serkowski relation in two ways: (1) they challenge the universality of the relation, because Eq. (1) can be a poor fit to the data when 3D effects are significant; and (2) even when Eq. (1) fits the data well, the best-fit parameters may not be representative of dust physics. Conversely, deviations of stellar polarization spectral profiles could potentially be utilized as a diagnostic of the 3D complexity of the magnetized ISM.

2 Theoretical expectations

In the case of multiple dust clouds along the LOS, initially unpolarized star light is affected sequentially by all of them, until it reaches the observer. In the limit of low polarization, the normalized Stokes parameters, q = Q/I, u = U/I, of the transversing light, are additive (e.g. Patat et al. 2010). Therefore, the observed q, u at a given wavelength for an intrinsically unpolarized source lying behind N_C number of clouds will be $$q_{obs}=\sum _{i=1}^{N_{\mathrm{C}}}{q}_i,$$(2)

and similarly for u, with q_i being the q induced by the i-th dust cloud in the LOS.

In order to explore this effect, we created analytical models for the cases of one and two dust clouds along the LOS, for different parameters of the clouds. Given that the Serkowski relation is an intrinsic property of individual ISM clouds, variations in grain sizes would be manifested by variations in λ_max, which is proportional to the average size of aligned grains (Draine 2024). Variations in the plane-of-sky (POS) morphology of the magnetic field are traceable through the polarization angle, θ. Variations in λ_max and in θ between two or more clouds along the LOS significantly affect the cumulative (observed) polarization signal. For our analytical models, we used Eq. (1) to calculate P_i (λ) for 0.1 μm ≤ λ ≤ 1.2 μm, selecting typical values for P_max and λ_max and fixing K = 1.68·λ_max (Wilking et al. 1980). We also assigned a θ_i value (constant with wavelength) representing the average POS magnetic field direction within the cloud. Thus, for each cloud, i, we have qi (λ) = P_i (λ)·cos 2θ_i and u_i(λ) = P_i(λ)·sin2θi. For the combined model, we calculated q_obs(λ), u_obs(λ) from Eq. (2) and finally $\hat{P}(\lambda ),\hat{\theta }(\lambda )$ as $$\hat{P}(\lambda )=\sqrt{q_{obs}^2(\lambda )+u_{obs}^2(\lambda )}$$(3a) $$\hat{\theta }(\lambda )=\frac{1}{2}\cdot {\tan }^{-1}\left(\frac{u_{obs}(\lambda )}{q_{obs}(\lambda )}\right),$$(3b)

where we used the two-argument arctan function.

In Fig. 1, we present some examples of P - λ and θ-λ profiles resulting from the presence of one or two clouds in the LOS, explored by our analytical models. The simple case of a single dust cloud (Serkowski relation - Eq. (1)) is shown in the left panel of Fig. 1 for reference. Three cases with two clouds in the LOS are shown in the middle left, middle right, and right panels of Fig. 1. Specifically, the middle left panel demonstrates an example where both the P - λ and θ - λ curves resemble the Serkowski relation, although the profiles are constructed by two components. The middle right panel demonstrates an example where the P - λ curve resembles the Serkowski relation, yet the θ - λ profile deviates from constancy. The right panel demonstrates an example where both the P-λ and the θ-λ profiles deviate from the Serkowski curves.

Assuming a difference of a factor of two in λ_max between our two hypothetical clouds yields a bimodal P(λ), which is significantly different than the shape of a typical Serkowski curve (Fig. 1, upper right panel). Although the difference appears to be significant, it would be challenging to identify itin practice without the use of spectropolarimetry. Current polarimetric research on the ISM predominantly employs wide-filter polarimetry, which has limited resolving power. Moreover, studies with multiband polarization measurements (such as the ones used in this work - see Sect. 3.1) are usually limited to a narrow range of frequencies.

Modest LOS variations in λ_max can result in even more complications. In our assumed case of two clouds, a 40% difference in λ_max between the two clouds leads to a profile that is nearly identical to the single cloud Serkowski curve profile (middle right panel in Fig. 1), although the cumulative signal by construction consists of two underlying profiles. In this case, erroneously assuming a single cloud along the LOS would lead to inaccurate constraints for the parameters of the Serkowski relation, and hence for the grain properties. The imprint of a second cloud in these cases can be noticed by the nonconstant behavior of θ with wavelength.

In cases in which the λ_max and θ parameters are similar between clouds (middle left panel in Fig. 1), the resulting profiles of both polarization and EVPA are indistinguishable from the single cloud case within observational uncertainties, given the current accuracy of optical polarimeters. Significant differences in the EVPAs of two clouds lying along the LOS, such as the ones employed in the middle right and rightmost models of Fig. 1, can be revealed through the variation in the EVPA as a function of λ. However, obtaining a sufficient number of data points that would allow us to probe such variations is challenging², and for this reason these data are sparse.

These analytic calculations pose a challenge in the interpretation of P(λ) observations: while the Serkowski relation would still hold for individual clouds, the majority of the results extracted from multiband polarization are likely contaminated by the 3D structure of the ISM, leading to inaccurate dust-grain constraints. The upside of this problem is that it opens up new possibilities for studying the magnetized ISM in 3D.

Fig. 1 Demonstration of the effect of 3D structure in the Serkowski relation, using analytical models for P - λ and θ-λ curves. In each case, the solid black curve corresponds to the combined model and the dashed red and dotted blue lines correspond to individual clouds. Left column: simple case of one cloud (Eq. (1)). Middle left column: example with two clouds in the LOS, in which both the combined P-λ, and θ-λ seem to follow the Serkowski relation. Middle right column: example with two clouds in the LOS, in which the combined P-λ seems to follow the Serkowski relation, but θ - λ not. Right column: Example with two clouds in the LOS in which both P-λ and θ - λ profiles deviate from Serkowski expectations. The parameters for each model are: Left: Pmax = 1%, λmax = 0.6μιη. θ = 30°. Middle left: Pmax1 = 1.5%, λmax1 = 0.55 μm, θ1 = 30°, Pmax2 = 0.7%, λmax2 = 0.5 μm, θ2 = 35°. Middle right: Pmax1 = 0.7%, λmax1 = 0.5 μm, θ1 = 0°, Pmax2 = 0.8%, λmax2 = 0.7 μm, θ2 = 60°. Right: Pmax1 = 1.3%, λmax1 = 0.4 μm, θ1 = 0°, Pmax2 = 1.4%, λmax2 = 0.8 μm, θ2 = 75°. For all clouds, we have used K = 1.68·λmax (Wilking et al. 1980).

3 Methods

3.1 Polarimetric data selection - Statistical sample

A recent agglomeration of multiband optical polarization data (Panopoulou et al. 2025), which is the largest to date, allows us, for the first time, to study in a statistical fashion the impact of the 3D ISM structure of the Galaxy on the Serkowski relation, and the consequences for magnetic field and dust physics studies. In order to fit multiple-cloud models to polarization data, it is optimal to use polarization measurements across a broad optical wavelength range. However, for most stars in the literature, measurements exist only in a single band or a few bands. Even so, it is possible to evaluate the validity of the Serkowski relation in a statistical manner, by studying an ensemble of sources. We designed such an experiment to test whether stars in singlecloud LOSs are better fit by the Serkowski relation than stars in multiple-cloud LOSs.

To acquire archival optical polarization measurements for the statistical comparison, we explored Tables 4 and 6 of the compiled polarization catalog of Panopoulou et al. (2025), which provide polarimetry and distances from Gaia EDR3 (Bailer-Jones et al. 2021), respectively, for ∼42000 sources. We excluded stars with distances greater than 1.25 kpc, as this is the limit of our dust map (see Sect. 3.2). As a distance, we took the average of the values (‘r_med_geo’, ‘r_med_photogeo’) provided in the catalog. We excluded measurements in filters with no well-defined effective wavelength (λ_eff); therefore, the excluded filters are “0,” “20,” “23,” and “51,” which correspond to “No filter or unclear,” “weighted-mean,” 0.735-0.804 μm, and 0.50.85 μm, respectively, in the catalog. For targets with multiple measurements in the same filter, we used the weighted mean and standard deviation of the q measurements, according to $$q_{\text{mean}}=\frac{\sum (q_i\cdot w_i)}{\sum w_i}$$(4) $$q_{\text{std}}=\frac{\sum w_i(x_i-\overline{x}{)}^2}{\left(\frac{N-1}{N}\right)\cdot \sum w_i},$$(5)

with $w_i=1/{\sigma }_{qi}^2$, and σ_qi the uncertainty of the measurement in , and similarly for u. We excluded measurements with debiased P/σ_P < 2. Additionally, we did not consider measurements with a polarization uncertainty of σ_P < 0.1% because they are likely underestimated. After all the cuts, we selected targets that have surviving measurements in at least three bands. We note that we did not impose any cut for intrinsically polarized stars (although the Panopoulou et al. (2025) catalog provides this information), as these targets may have been flagged as intrinsically polarized simply because they do not follow the Serkowski relation. However, the final sample does not contain any star that is flagged as intrinsically polarized. The final sample contains 223 sources with measurements in at least three optical bands.

The reasoning for excluding data with σ_P < 0.1% is presented below. Initially, we considered not imposing any restrictions on the level of σ_P. However, we found that error bars below ∼0.1% are likely underestimated on average. If the error bars accurately represented the uncertainties in the measurements, the ${\chi }_{\nu }^2$ would remain independent of the error bar level when analyzing targets within the same subsample (i.e., those behind the same number of clouds). To test this hypothesis, we examined the targets in the subgroup with N_C = 1. We chose this particular subsample, as we know not only that the ${\chi }_{\nu }^2$ should be flat with respect to the error bars, but also that it should be ∼ 1 for targets behind a single cloud. We fit the model described in Sect. 4.1 and, for each target, plotted the average error bar value across different filters against the corresponding ${\chi }_{\nu }^2$. We demonstrate in the left panel of Fig. 2 that targets with smaller error bars exhibit higher ${\chi }_{\nu }^2$, revealing a clear, nonconstant correlation between the two quantities. χ2 seems to approach the expected value of ∼1 for ∼σ_P ≥ 0.1%. Based on the above findings, we gradually excluded points with lower reported uncertainties and repeated the fit described in Sect. 4.1, until the expected behaviors were satisfied. We found that the dependence practically disappears and ${\chi }_{\nu }^2$ approaches the expected value of ∼ 1 for targets with average σ_P > 0.1% (Fig. 2, right). Therefore, we discarded all measurements with σ_P < 0.1%.

Fig. 2 Left: average reported uncertainties in polarization, $\overline{{\sigma }_P}$, versus ${\chi }_{\nu }^2$ for targets that survive the cuts described in Sect. 3.1, but without any restriction for the error bars. The black line is a linear fit for all the data in log-log space. The dashed blue line is similar, but excluding the points in the top left of the plot. Right: similar to left, but with the restriction for the error bars of σP ≥ 0.1%.

3.2 Extinction

We collected extinction data for our target stars from the Edenhofer et al. (2024) dust map. We chose this map, as it provides a good angular resolution (14′) as well as a parsec-scale distance resolution. The map is available in two versions, one that extends out to 1.25 kpc from the Sun and one that extends up to 2 kpc, but uses less data for the dust profile reconstruction. We used the 1.25 kpc version, as it should be more accurate, and we did not want to introduce any additional sources of noise. Extinction values in the Edenhofer et al. (2024) map are expressed in the units used by Zhang et al. (2023), the work upon which the dust map was based. We converted the extinction values to represent the extinction in the V band (A_V) based on the conversion table provided by Zhang et al. (2023).

For each of our targets, we examined the “mean” (Edenhofer et al. 2024) cumulative and differential A_V versus distance (D) profiles by eye, in order to identify the number of significant steps in the profile, which presumably correspond to distinct dust clouds (N_C). For targets for which more than three clouds were identified, or the dust profile was very complex, we set N_C > 3. The number of target stars in each subgroup with a different number of identified clouds along the LOS are #(N_C = 1) = 49, #(N_C = 2) = 74, #(N_C = 3) = 70, #(N_C > 3) = 30. The targets that correspond to each N_C are distributed over different positions on the sky (Fig. 3).

We considered the possibility of using a more objective metric to characterize the number of dust clouds along each LOS, as each step in the dust profiles is typically not equally prominent, and in many cases numerous small steps contribute to the cumulative extinction. Therefore, we explored a method similar to the one used by Panopoulou & Lenz (2020), who identified the number of clouds in H I data by detecting peaks in the differential H I column density profile and then normalizing each peak by the height of the tallest peak among them. In our case, the normalized number of clouds in each LOS would be $N_{\mathrm{C}}^{norm}=\sum \frac{A_V^i}{A_V^{max}}$ with $A_V^i$ the extinction of a given peak in the LOS, and $N_{\mathrm{C}}^{norm}$ the extinction of the highest peak in the differential A_V versus D profile. However, this was not useful in our study for the following reasons. 1) $N_{\mathrm{C}}^{norm}$ is ambiguous. For example, $N_{\mathrm{C}}^{norm}=2$ could signify two equally prominent peaks, or one prominent peak plus two peaks half its height, or even one prominent peak plus three peaks a third of its height, and so on. This lack of clarity reduces its utility. 2) The profiles in the Edenhofer et al. (2024) map are sometimes irregular, with many minor peaks that may not represent real features but instead result from sampling noise. Additionally, peaks can appear very close together, making it difficult to distinguish between separate clouds and variations within the same cloud. Although we considered smoothing the profiles to address these issues, we found that the number and relative height of the resulting peaks were highly sensitive to the parameters chosen for smoothing, introducing further ambiguity. 3) Even when one peak appears higher than another, this does not necessarily correlate with increased polarization. The dust-induced polarization of aligned grains depends largely on factors beyond peak height, such as dust density, grain alignment efficiency, magnetic field coherence, and the inclination of the magnetic field with the LOS.

In our attempts to devise an objective method of peak counting, we found that the resulting numbers were often ambiguous or lacked physical meaning. Consequently, we opted for a visual approach. Although this method has its uncertainties, it is self-consistent and allows statistical patterns to emerge in a dataset as large as ours. Figure 4 illustrates examples of dust profiles that we assigned to different subsamples.

Fig. 3 Positions of targets depending on their NC. Maps are in Galactic coordinates, in mollweide projection, centered on (0,0). Parallels are drawn every 30°. Meridians are drawn every 60°. Longitude increases toward the right-hand side. Colors correspond to the average λmax of each area. Areas marked with “A” through “G” contain stars that are located behind a single cloud and multiple clouds. See discussion in Sect. 4.3.

Fig. 4 Dust profiles for four example targets belonging in four different subsamples. From left to right: NC = 1, NC = 2, NC = 3, NC > 3. Top row: cumulative AV vs. D. Bottom row: differential AV vs. D.

3.3 Fitting the Serkowski relation

In order to evaluate the performance of the Serkowski relation in our statistical sample, we used an extremely powerful, if unconventional in Serkowski studies, technique: we fit our data in the Stokes parameters space q - u, exploiting its power to simultaneously trace changes in the profiles of both the degree of polarization and EVPA. Traditionally, the EVPA is not taken into account while fitting the Serkowski formula, and thus part of the information imprinted in the polarization signal is ignored. Although the variability of the EVPA with wavelength has been hypothesized to be connected to the existence of multiple clouds (e.g. McMillan & Tapia 1977; Whittet 2015), q - u fitting has never been employed to trace the various components (clouds) contributing to the P - λ and θ - λ profiles.

As a first step, we converted Eq. (1) to the q - u space. This step allows us to include θ in the model, which is ignored in the Serkowski formula. Therefore, our model becomes $$\hat{q}=P_{max}\cdot \exp (-K\cdot {\ln }^2\frac{{\lambda }_{max}}{\lambda })\cdot \cos (2\theta ),$$(6a) $$\hat{u}=P_{max}\cdot \exp (-K\cdot {\ln }^2\frac{{\lambda }_{max}}{\lambda })\cdot \sin (2\theta ).$$(6b)

The quantity we wished to minimize for the fits is $$m=L_q+L_u+\frac{L_q}{L_u}+\frac{L_u}{L_q}$$(7)

with $$L_q=\frac{1}{N_{\lambda }}\cdot \sum {\left(\frac{q_i-{\hat{q}}_i}{{\sigma }_{q_i}}\right)}^2,$$(8a) $$L_u=\frac{1}{N_{\lambda }}\cdot \sum {\left(\frac{u_i-{\hat{u}}_i}{{\sigma }_{u_i}}\right)}^2,$$(8b)

where N_λ is the number of different bands for which we had measurements for each target. q_i, u_i is the measurement in a given band, and $\widehat{q_i}$, $\widehat{u_i}$ our model introduced in Eqs. (6a), (6b). σ_qi, σ_ui are the uncertainties of the measurements in q and u, respectively. The terms L_q and L_u are the reduced χ² metric for q and u, respectively. The terms L_q/L_u and L_u/L_q were introduced to ensure a balanced contribution from L_q and L_u . For example, if the objective was to minimize just L_q + L_u, the optimal value V could be achieved with L_q = V and L_u = 0, L_u = V and L_q = 0, or any combination in between. By introducing L_q/L_u and L_u/L_q, we penalize the model when it favors solutions where L_q = 0 or L_u = 0, and thus promote a more balanced fit.

For the modeling procedure, we used the OPTUNA package in Python (Akiba et al. 2019). OPTUNA is an automatic hyperparameter optimization software framework, originally developed for machine learning hyperparameter-tuning purposes. However, it can also address general optimization problems through the use of a Bayesian optimization framework. We found it highly effective in multiparameter optimization with few data points and many parameters, such as the case of the q - u modeling, and more efficient in finding the best solution than traditional algorithms such as Markov chain Monte Carlo (MCMC). We used the default options of OPTUNA with n_trials = 5000 for our fits. Moreover, we adopted the following priors. For P_max, we initially assumed a rather loose Gaussian prior distribution with a mean of 3% and standard deviation of 3%, and enforced it to be positive. For λ_max, we assumed a Gaussian prior distribution with a mean of 0.6 μm and standard deviation of 0.2 μm, and ensured it was always positive. For K, we imposed a Gaussian prior distribution with a mean of 1.68·λ_max and standard deviation of 0.13·λ_max (Wilking et al. 1980), and ensured it was always positive. For θ, we assumed a Gaussian prior distribution with the mean the value derived from the average q and u as in Eq. (4) and converted to θ using Eq. (3b), and a 10^° standard deviation.

We evaluated the quality of our fits by measuring the χ² per degree of freedom (DoF): $${\chi }_{\nu }^2\equiv \frac{{\chi }^2}{DoF}=\frac{\sum {\left(\frac{q_i-{\hat{q}}_i}{{\sigma }_{q_i}}\right)}^2+\sum {\left(\frac{u_i-{\hat{u}}_i}{{\sigma }_{u_i}}\right)}^2}{2\cdot N_{\lambda }-\nu },$$(9)

where ν = 4 is the number of free parameters in the model (P_max, λ_max, K, θ). We used 2·N_λ as the number of points in Eq. (9) as we accounted for the measurements in both q and u to perform the simultaneous fits. We note that we cannot use m directly (Eq. (7)) to evaluate the fits, as m is not a normalized metric, and thus cannot be used for comparing fits to each other.

3.4 Fitting multiple clouds

In order to fit multiple clouds (as in the examples of Sect. 4.2), we again used the q - u fitting technique, following the same logic discussed in Sect. 3.3. The multiple cloud model is derived from Eqs. (2) and (6), by adding up the q and u components of the individual clouds.

When trying to determine the correct parameters for a model with two or more clouds, the number of parameters is increased (four parameters per cloud), resulting in a substantial number of possible combinations and a vast parameter space to explore. In such situations, identifying the optimal combination of parameters that minimize the cost function can be challenging. To address this, we employed the following procedure. Initially, instead of directly using OPTUNA, we utilized an MCMC algorithm to rapidly explore the parameter space. Given that the two clouds are interchangeable, we imposed the condition λ_max1 > λ_max2 to focus the algorithm on a single solution. Subsequently, we visually examined the corner plots of the parameters generated by the MCMC algorithm to establish reasonable constraints for the OPTUNA algorithm, thereby significantly narrowing the search space. Finally, we ran the fit using OPTUNA, imposing Gaussian priors on the parameters, based on the mean and the spread of the parameters identified by the MCMC algorithm. We chose OPTUNA to derive the final set of parameters for the model, as it directly provides the optimal set of parameters by minimizing Eq. (7). The fits were evaluated with ${\chi }_{\nu }^2$ (Eq. (9)).

4 Results

4.1 Statistical sample

We grouped the sample of sources found in Sect. 3.1 into subsamples depending on the number of clouds (N_C) lying in their LOS, as derived from the 3D extinction map. We obtained four main subsamples: targets in LOSs with one dominant cloud, N_C = 1; two clouds, N_C = 2; three clouds N_C = 3; and more than three clouds, N_C > 3. For this sample, we fit the traditional Serkowski curve (Eq. (1)) for all targets, regardless of N_C, to evaluate the capability of the Serkowski relation to fit polarimetric data for subsamples with different N_C . We evaluated the performance of the Serkowski relation in all of the targets as described in Sect. 3.3. We constructed the cumulative density functions (CDFs) of ${\chi }_{\nu }^2$ for the different subsamples (Fig. 5), and we applied the Kolmogorov-Smirnov (KS) test to identify any statistical discrepancies between the ${\chi }_{\nu }^2$ distributions of the N_C = 1 subset and the other three subsets.

We found that the ${\chi }_{\nu }^2$ distribution of the N_C = 1 subset is significantly different from those of subsets with N_C = 2 or N_C = 3. In contrast, the ${\chi }_{\nu }^2$ distributions of the N_C = 1 subset and of the N_C > 3 subset are not discrepant. The p-values of the KS tests for all comparisons are quoted in Table 1, along with the median ${\chi }_{\nu }^2$ for the fits of each subsample. Despite the sparse nature of the data (measurements in only ∼3.5 different filters on average per target), the sample is sufficiently large to allow for the statistical distinction of LOSs with two or three clouds from those with only a single cloud. The deterioration of the Serkowski fits in the multicloud LOSs is mainly driven by the behavior of the EVPA with wavelength, and is made clear thanks to our q - u space fits. In contrast, if we perform the fits only on P - λ, it is not possible to discriminate between LOSs with and without multiple clouds with current data. When the structure in the LOS is too complicated (N_C > 3), there is no statistical difference in the polarization profiles compared to the single cloud case. This is not surprising, as with an increasing number of clouds the variations induced by individual screens tend to average out. Even deviations in the EVPA, which should be more prominent, would become smeared and difficult to detect, especially with limited data points across the spectrum.

Fig. 5 CDFs of ${\chi }_{\nu }^2$ for different subsamples, as is described in each label.

4.2 Two revealing examples

To demonstrate the impact of the LOS variations in λ_max and EVPA on the Serkowski relation, we employed archival spec-tropolarimetric data (Bagnulo et al. 2017). We identified two revealing cases, where the multiwavelength polarization data show clear evidence of the LOS effects. The examples discussed here are not part of our statistical sample, as their distances are higher than the limitations imposed in Sect. 3.1.

The LOS toward the main sequence star HD93222 with spectral type O7 represents a prominent example of what the P - λ, and θ-λ profiles look like when 3D effects are significant. First, the polarization as a function of wavelength deviates from the typical Serkowski curve (Fig. 6 middle panel). Fitting the P - λ data with the Serkowski relation requires unphysical values for both λ_max, and K (Table 2), which still yield a merely adequate fit. In addition, θ varies smoothly with wavelength, λ, which is unexpected when the polarization is dominated by a single cloud with a well-defined mean magnetic field. Thus, both P and EVPA show an atypical behavior with λ for this star. Using 3D dust extinction maps (Edenhofer et al. 2024), we verified the existence of multiple clouds along the LOS toward the HD93222 star. The clouds are identified as abrupt increments in extinction as a function of distance (left panel in Fig. 6). The average extinction profile shows three prominent clouds located at around 0.25, 1.35, and 1.8 kpc. The target star, whose estimated distance is 2.44 kpc (Bailer-Jones et al. 2021), lies beyond these clouds, and hence its polarization signal carries information from all of them. We fit the data with models with different number of clouds. Although the extinction profile exhibits three distinct steps, the data can be described equally well with both a two- and a three-cloud model. This may indicate that two out of three clouds have similar properties, and their constructive profile resembles that of one cloud (such as the case of the middle left panel of Fig. 1). We proceeded with the simplest case of a two-cloud model.

As a second example, we show results for the star HD163181, which is an intermediate-size luminous supergiant of spectral class O9.5 (Fig. 7). In this case, the observed degree of polarization follows the Serkowski formula, but the EVPA varies with wavelength. The 3D extinction profile has two significant steps, indicating the existence of more than one dominant cloud along this LOS. Assuming the standard Serkowski curve with a single dominant cloud can yield good fits to the data in the P - λ space (Bagnulo et al. 2017), but it fails to explain the observed variability of the EVPA; we verified this by fitting a one cloud model in the q - u plane (dashed blue line in Fig. 7). However, when we considered two clouds along this LOS, as the 3D extinction map suggests, we were able to capture the observed variability of both P and θ with wavelength (dash-dotted red line in Fig. 7, middle panel).

4.3 The parameters of the Serkowski formula

Here, we demonstrate how the derived parameters of the Serkowski relation can be incorrect, when erroneously assuming a single cloud in the LOS. We explored how 3D effects affect the distributions of λ_max, and P_max that one obtains through the classical fitting of the Serkowski relation in the P-λ space; K linearly correlates with λ_max (Wilking et al. 1980), and thus can be omitted from this analysis. We fit the data only in P - λ space (Eq. (1)) by minimizing $$\frac{1}{N_{\lambda }}\cdot \sum {\left(\frac{P_i-\widehat{P_i}}{{\sigma }_{p_i}}\right)}^2,$$(10)

assuming a one-cloud model, and imposing the same priors as discussed in Sect. 3.3. For these fits, we used the debiased P values. Since the effect of the 3D dust structure is more prominent in θ - λ space, fits in P-λ were comparably good for all four subsamples. We then explored the normalized distributions of the derived parameters for all subsamples (Fig. 8 for λ_max, Fig. 9 for P_max/AV).

Past studies of the Serkowski parameters show that the distribution of λ_max is singly peaked, with an average value close to 0.5 μm (Martin & Whittet 1990). In contrast, we find that for N_C = 1 (single-cloud sightlines) the intrinsic distribution of λ_max is bimodal, with peaks at approximately 0.28 and 0.62 μm. On the other hand, all distributions with N_C > 1 (multiple-cloud sightlines) are unimodal and peak at around 0.55 μm, which is close to the value that is considered to be the Galactic average. By performing MCMC simulations, we found that the unimodality in the λ_max distribution, which is observed for N_C > 1, can be reproduced with multiple-cloud models with values of λ_max randomly drawn from the N_C = 1 distribution. The intrinsic Galactic distribution of λ_max may thus be bimodal, and the currently accepted value for the Galactic average of λ_max (∼0.5 μm) is likely contaminated by LOS effects. Notably, the larger spread of the distribution of λ_max for the NC = 2 case is not surprising. Clarke & Al-Roubaie (1984) theorized that in the case of two clouds in the LOS, the acquired λ_max from the Serkowski fit is different from the average λ_max of the two clouds when the two clouds have different parameters. Furthermore, it is worth considering that most areas containing stars behind a single cloud also contain stars behind multiple clouds. We calculated the average λ_max for targets in distinct areas on the sky (Fig. 3). In most cases, the average λ_max differs significantly between subsets (see Table 3) located in the same area. This suggests that using the λ_max of a star from an area with NC > 1 as a characteristic value for the LOS dust properties could be misleading.

The distribution of P_max/AV is also affected by the number of clouds along the LOS. We can immediately see this by assuming, without loss of generality, that two clouds lie along a LOS. The polarization induced by the two clouds can either (partially or totally) add up, or (partially or totally) cancel out, depending on the difference between the EVPAs each would impart on its own. When the EVPA difference is close to 0^o , the polarization adds up, while when the difference is close to 90° it cancels out. On the other hand, the extinction, AV, induced by each cloud is always additive. Therefore, with more clouds, we would expect to have in general lower P_max /AV than in the N_C = 1 cases. However, since P_max is a derived parameter rather than a direct measurement, the fit can be poor in cases with N_C > 1, and P_max may not represent any physical truth. In such scenarios, P_max, and thus P_max/AV, could take on arbitrary values. This is exactly what we observe in the distributions of P_max/AV (Fig. 9). The P_max/AV distribution for N_C = 1 is free from LOS effects, and hence it should resemble the “ground truth.” This distribution peaks at $P_{max}/A_V\simeq 3-3.5\mathrm{\%}/\mathrm{mag}$, with only a few points with P_max/AV > 4.5%/ mag. On the other hand, the N_C = 2 distribution peaks at a lower value of $P_{max}/A_V\simeq 2\mathrm{\%}/\mathrm{mag}$, its tail extending to higher values of P_max /AV than when N_C = 1. The shape of the N_C = 2 distribution exactly follows the expectations discussed above. Similarly, for the distributions of more components, P_max /AV consistently has lower values than the case of N_C = 1.

The polarization-to-extinction ratio has garnered significant interest within the ISM community (e.g. Hiltner 1956; Serkowski et al. 1975; Panopoulou et al. 2019; Planck Collaboration XII 2020; Angarita et al. 2023). While here we present our findings for P_max/AV , discussions typically focus on P/AV in the context of measured polarization in a specific band, rather than on a fit parameter such as P_max . Consequently, our results for the distributions with N_C > 1 cannot be fully integrated into this broader discussion if the fit parameter, P_max , does not accurately reflect the true value for a given field. On the contrary, the findings for the N_C = 1 case should be free from bias. Within this subsample, we identify three targets with extreme P_max /AV values: two with $P_{max}/A_V\simeq 10-11\mathrm{\%}$/mag and one with $P_{max}/A_V\simeq 6\mathrm{\%}/\mathrm{m}\mathrm{a}\mathrm{g}$. These values correspond to $P_{max}/E(B-V)\simeq 31-34\mathrm{\%}/\mathrm{m}\mathrm{a}\mathrm{g}$ and $P_{max}/E(B-V)\simeq 18.6\mathrm{\%}/\mathrm{m}\mathrm{a}\mathrm{g}$, respectively, assuming that RV = 3.1 with RV = AV/E(B – V) (Fitzpatrick 2004). These values exceed the highest P_max/E(B – V) reported by Angarita et al. (2023) to date. However, the observed maximum P/E(B -V) has consistently been revised upwards over time: from P/E(B – V) ≤ 9%/mag (Hiltner 1956; Serkowski et al. 1975), to P/E(B – V) ≤ 13%/mag (Panopoulou et al. 2019; Planck Collaboration XII 2020), and most recently to P/E(B -V) < 16%/mag (Angarita et al. 2023). While our values may appear extreme, they may well reflect true physical conditions, and further research is warranted.

Overall, both the λ_max and the P_max/A_V are significantly affected by the 3D structure of the ISM, which, if ignored, can lead to erroneous Serkowski best-fit parameters. Our results suggest that fitting a single Serkowski relation is invalid in LOSs with multiple clouds. There can be cases in which a single Serkowski relation yields good fits, even when 3D effects are important (Sect. 4.2). However, in these cases the best-fit parameters of a single Serkowski curve are misleading, and not representative of the dust grain physics.

All the targets used in the statistical sample are presented in Table 4 with their properties. The complete table can be found only in electronic form (Sect. 6).

Fig. 6 Example of a case in which both the P - λ and the θ - λ profiles deviate from Serkowski expectations. We were able to fit both profiles by introducing a 2-Cloud model. Data from Bagnulo et al. (2017) for HD 93222. Left: extinction vs. distance profile taken from the 2kpc version of the map by Edenhofer et al. (2024). The vertical red line corresponds to the distance of the source. Middle: data in the P - λ space, together with three fit models. Solid black line: Fit using Eq. (1) in P - λ space only. Dashed blue line: Best-fit 1-Cloud model in the q - u space (Eq. (6)). Dash-dotted red line: Best-fit 2-Cloud model in the q - u space (Eq. (6)). Dotted orange lines correspond to the two individual components of the 2-Cloud model. Right: data in θ - λ space together with the three models and the individual components of the 2-Cloud model. Parameters of the fits are quoted in Table 2.

Fig. 7 Example of a case in which the P – λ profile follows the Serkowski relation, but the EVPA does not. By employing a model with two components in the LOS, it is possible to fit both P and θ. Data from Bagnulo et al. (2017) for HD 163181. Panels and lines similar to Fig. 6. Parameters of the fits are quoted in Table 2.

Fig. 8 Normalized distributions of λmax for subsamples with different number of clouds in the LOS. See Sect. 4.3 for details.

Fig. 9 Normalized distributions of Pmax for subsamples with different number of clouds in the LOS. See Sect. 4.3 for details. The following bins are omitted for clarity. For NC = 2: bin at $P_{max}/A_V\simeq 23$ with a height of 0.0135. For NC = 3: bin at $P_{max}/A_V\simeq 13$ with a height of 0.0145.

5 Discussion

5.1 Implications for dust modeling

The findings presented in this study have significant implications for dust modeling. Dust models and grain alignment theories must satisfy the constraints obtained from Serkowski fits (Martin & Whittet 1990; Andersson et al. 2015; Draine 2024). However, our results indicate that the Serkowski parameters derived from fitting the observed polarization data in complex interstellar environments may not accurately represent the underlying dust properties, leading to erroneous conclusions about grain size distributions and composition. The results from Sect. 4.3 suggest that the intrinsic distribution of λ_max in the clouds of our Galaxy may actually be bimodal (Fig. 8). This finding contradicts past studies (Martin & Whittet 1990), which likely stem from (incorrectly) fitting the Serkowski formula in regions with multiple clouds. Although, our employed sample is the largest to date, the statistics are still limited, and more data are required to further explore the intrinsic distribution of λ_max in our Galaxy.

We propose the following two options for constraining the populations of dust in different ISM environments: 1) fit the Serkowski formula only to targets behind a single cloud, or 2) decompose the cumulative polarization signal into its individual components for targets behind multiple clouds to derive the characteristics of each. Moreover, different regions of the Galaxy (e.g., the disk vs. the halo, or varying distances from the Galactic center) may exhibit distinct dust physics that we are currently missing due to the misapplication of the Serkowski relation. By studying the polarization of targets behind multiple clouds and decomposing the signal to derive individual parameters, it may become possible to uncover dust clouds with diverse and potentially novel properties in different parts of the Galaxy.

5.2 Implications for the identification of intrinsically polarized stars

It is customary to interpret deviations from the Serkowski relation in the P(λ) behavior of point sources as indicative of intrinsic polarization mechanisms, such as the presence of circumstellar disks (e.g. Topasna et al. 2023). However, the wavelength dependence of polarization in such cases is known to lack a typical pattern (Bastien 2015). In contrast, we have demonstrated that fitting simple two-cloud models can account for deviations from the Serkowski formula quite effectively (Sect. 4.2). This suggests that deviations from the Serkowski curve alone are insufficient to indicate the presence of an intrinsic polarization component. In order to confirm that a star not adhering to the Serkowski formula is intrinsically polarized, other evidence, such as the detection of variability in the polarization or a difference in polarization between the star and its neighbors, is required.

5.3 The possibility of 3D tomography with the Serkowski relation

Fitting multiple-cloud models to multiwavelength (optical to near-infrared) polarization data paves the way for constraining the LOS variations in the plane-of-the-sky magnetic field morphology (tomography). If the number of clouds in the LOS is known, for instance from 3D extinction maps (e.g., Edenhofer et al. 2024; Green et al. 2019), it is possible to extract the parameters of individual clouds using the techniques discussed in Sect. 3.4. However, the following considerations should be kept in mind. It is possible that two or more clouds in the LOS have similar parameters. In such cases, the polarization data may be adequately fit with fewer clouds than what the 3D extinction maps suggest. For example, if two clouds along aLOS have comparable λ_max and θ, their polarization contributions will combine, resulting in a profile resembling that of a single cloud but with an increased degree of polarization. This might explain our target star HD93222 (Sect. 4.2), where the extinction profile displays at least three distinct steps, yet the data are sufficiently well fit by a 2-Cloud model. Another possibility is that the third cloud has a low abundance of aligned grains, because of, for example, a magnetic field with a direction along the LOS.

Another strategy one could follow to identify cases of clouds with comparable properties is to obtain multiband or, ideally, spectropolarimetric measurements of multiple stars at different distances along the same LOS. The polarization signal would vary with distance, as each star would only be influenced by the foreground dust column. Thus, stars lying behind the same clouds would show similar polarization trends, while distinct changes would only occur for stars behind different clouds. Following such methods could complement the BISP-1 algorithm (Pelgrims et al. 2023), which employs single-band optical polarization data with parallaxes from Gaia. This algorithm is particularly relevant for performing ISM tomography along many LOSs (Pelgrims et al. 2024), especially in the context of optopolarimetric surveys, such as the Pasiphae survey (Tassis et al. 2018). While the method described by Pelgrims et al. (2023) can determine the number of clouds along a LOS and their distances, it does not fully capture the dust physics within each cloud. Multiband polarimetric measurements of a few stars at varying distances along the LOS can provide insights into the dust physics within each cloud. Thus, multiwavelength polarization data promises to provide an independent method for performing ISM magnetic field tomography, which could usher in the era of high-precision 3D magnetized ISM cartography.

6 Summary and conclusions

In this paper, we showcase the effect of the 3D structure of the ISM on the well-known and widely used Serkowski relation. Our findings are summarized as follows:

1. We have revisited theoretically how the Serkowski formula could become invalid in LOSs with many dust components. We have shown that the variations in both the degree of polarization and EVPA with wavelength can be severely affected in some stars. In other stars, the polarization profile may misleadingly appear to follow the Serkowski relation, despite the contribution from multiple dust clouds. In such cases, the imprint of the 3D structure can mainly be observed in the θ – λ profile;

2. As a test case, we selected two different targets that appear to be behind multiple dust clouds, based on their extinction profiles, and for which archival spectropolarimetric data exist. We demonstrated that the single-cloud model, implicitly assumed for Serkowski fits, fails in these cases, while we successfully fit two-cloud models;

3. We propose the q – u fitting technique as the most appropriate approach to fit the Serkowski relation. This method allows one to trace both P and θ as functions of wavelength simultaneously, enabling fitting for more than one cloud along the LOS;

4. 4. We used a sample of 223 stars behind different numbers of clouds with polarimetric measurements in three or more bands, to fit and evaluate statistically the performance of the Serkowski formula. We found that the Serkowski formula is not a good model for LOSs with two or more clouds, while it fits the data well in the cases of one or more than three clouds in the LOS. The latter is an outcome of averaging;

5. Our analysis of the fit parameters of the 223 stars suggests that the intrinsic distribution of the λ_max parameter in our Galaxy may be bimodal, something that was veiled until now, due to treating sightlines with many clouds as a single dust component;

6. A poor Serkowski fit does not by itself constitute conclusive evidence of intrinsic polarization in a source;

7. Fitting the Serkowski formula for stars behind multiple dust clouds leads to an incorrect estimation of the fit parameters, and thus of the underlying dust physics;

8. Our results hint at the possibility of performing magnetic tomography using multiband polarization data in combination with Gaia distances.

In conclusion, the Serkowski relation should be cautiously applied when the ISM 3D structure is complex. Finally, our findings underscore the critical role of spectropolarimetry in driving future advancements in ISM research. While wide-field polarimetry offers valuable insights, its limitations become evident when considering the physics of dust. Therefore, we strongly advocate for the widespread adoption of spectropo-larimetry within the ISM community for future research.

Data availability

Table 4 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/698/A168

Acknowledgements

The authors thank V. Pavlidou and G. Panopoulou for their constructive feedback. The authors also thank B. Hensley and A. Readhead for very fruitful discussions. NM and KT were supported by the European Research Council (ERC) under grant agreements No. 7712821. This work was supported by NSF grant AST-2109127. N.M. was funded by the European Union ERC-2022-STG – BOOTES – 101076343. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. Support for this work was provided by NASA through the NASA Hubble Fellowship grant # HST-HF2-51566.001 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

References

All Tables

All Figures

Fig. 1 Demonstration of the effect of 3D structure in the Serkowski relation, using analytical models for P - λ and θ-λ curves. In each case, the solid black curve corresponds to the combined model and the dashed red and dotted blue lines correspond to individual clouds. Left column: simple case of one cloud (Eq. (1)). Middle left column: example with two clouds in the LOS, in which both the combined P-λ, and θ-λ seem to follow the Serkowski relation. Middle right column: example with two clouds in the LOS, in which the combined P-λ seems to follow the Serkowski relation, but θ - λ not. Right column: Example with two clouds in the LOS in which both P-λ and θ - λ profiles deviate from Serkowski expectations. The parameters for each model are: Left: Pmax = 1%, λmax = 0.6μιη. θ = 30°. Middle left: Pmax1 = 1.5%, λmax1 = 0.55 μm, θ1 = 30°, Pmax2 = 0.7%, λmax2 = 0.5 μm, θ2 = 35°. Middle right: Pmax1 = 0.7%, λmax1 = 0.5 μm, θ1 = 0°, Pmax2 = 0.8%, λmax2 = 0.7 μm, θ2 = 60°. Right: Pmax1 = 1.3%, λmax1 = 0.4 μm, θ1 = 0°, Pmax2 = 1.4%, λmax2 = 0.8 μm, θ2 = 75°. For all clouds, we have used K = 1.68·λmax (Wilking et al. 1980).

In the text

Fig. 2 Left: average reported uncertainties in polarization, $\overline{{\sigma }_P}$, versus ${\chi }_{\nu }^2$ for targets that survive the cuts described in Sect. 3.1, but without any restriction for the error bars. The black line is a linear fit for all the data in log-log space. The dashed blue line is similar, but excluding the points in the top left of the plot. Right: similar to left, but with the restriction for the error bars of σP ≥ 0.1%.

In the text

Fig. 3 Positions of targets depending on their NC. Maps are in Galactic coordinates, in mollweide projection, centered on (0,0). Parallels are drawn every 30°. Meridians are drawn every 60°. Longitude increases toward the right-hand side. Colors correspond to the average λmax of each area. Areas marked with “A” through “G” contain stars that are located behind a single cloud and multiple clouds. See discussion in Sect. 4.3.

In the text

	Fig. 4 Dust profiles for four example targets belonging in four different subsamples. From left to right: NC = 1, NC = 2, NC = 3, NC > 3. Top row: cumulative AV vs. D. Bottom row: differential AV vs. D.
In the text

	Fig. 5 CDFs of ${\chi }_{\nu }^2$ for different subsamples, as is described in each label.
In the text

Fig. 6 Example of a case in which both the P - λ and the θ - λ profiles deviate from Serkowski expectations. We were able to fit both profiles by introducing a 2-Cloud model. Data from Bagnulo et al. (2017) for HD 93222. Left: extinction vs. distance profile taken from the 2kpc version of the map by Edenhofer et al. (2024). The vertical red line corresponds to the distance of the source. Middle: data in the P - λ space, together with three fit models. Solid black line: Fit using Eq. (1) in P - λ space only. Dashed blue line: Best-fit 1-Cloud model in the q - u space (Eq. (6)). Dash-dotted red line: Best-fit 2-Cloud model in the q - u space (Eq. (6)). Dotted orange lines correspond to the two individual components of the 2-Cloud model. Right: data in θ - λ space together with the three models and the individual components of the 2-Cloud model. Parameters of the fits are quoted in Table 2.

In the text

	Fig. 7 Example of a case in which the P – λ profile follows the Serkowski relation, but the EVPA does not. By employing a model with two components in the LOS, it is possible to fit both P and θ. Data from Bagnulo et al. (2017) for HD 163181. Panels and lines similar to Fig. 6. Parameters of the fits are quoted in Table 2.
In the text

	Fig. 8 Normalized distributions of λmax for subsamples with different number of clouds in the LOS. See Sect. 4.3 for details.
In the text

	Fig. 9 Normalized distributions of Pmax for subsamples with different number of clouds in the LOS. See Sect. 4.3 for details. The following bins are omitted for clarity. For NC = 2: bin at $P_{max}/A_V\simeq 23$ with a height of 0.0135. For NC = 3: bin at $P_{max}/A_V\simeq 13$ with a height of 0.0145.
In the text