The Serkowski relation is the cornerstone of studies of starlight polarization as a function of wavelength. Although empirical its extensive use since its inception to describe polarization induced by interstellar dust has elevated the relation to the status of an indisputable law serving as the benchmark for validating interstellar dust grain models. We revisit the effects of the 3D structure of the interstellar medium ISM on the wavelength dependence of interstellar polarization. We use analytical models to show how the wavelength dependence of both the polarization fraction and direction is affected by the presence of multiple clouds along the line of sight LOS accounting for recent developments in dust distribution modelling and utilizing an expanded archive of stellar polarization measurements. We highlight concrete examples of stars whose polarization profiles are severely affected by LOS variations of the dust grain and magnetic field properties and we provide a recipe to accurately fit multiple cloud Serkowski models to such cases. We present for the first time compelling observational evidence that the 3D structure of the magnetized ISM often results to the violation of the Serkowski relation. We show that 3D effects impact interstellar cloud parameters derived from Serkowski fits. In particular the dust size distribution in single  cloud sightlines may differ from analyses that ignore 3D effects with important implications for dust modelling in the Galaxy. Our results suggest that multiband stellar polarization measurements offer an independent probe of the LOS variations of the magnetic field constituting a valuable new tool for the 3D cartography of the ISM. We caution that unless 3D effects are explicitly accounted for a poor fit to the Serkowski relation does not by itself constitute conclusive evidence that a star is intrinsically polarized.
['Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Solar and Stellar Astrophysics']
3D ISM structure challenges the Serkowski relation
{'3D ISM structure challenges the Serkowski relation': 'N. Mandarakas 1 , 2, ⋆ , K. Tassis 1 , 2 , and R. Skalidis 3  \n- 1 Department of Physics, University of Crete, Vasilika Bouton, 70013 Heraklion, Greece\n- 2 Institute of Astrophysics, Foundation for Research and Technology - Hellas, 100 Nikolaou Plastira, Vassilika Vouton, 70013 Heraklion, Greece\n- 3 Owens Valley Radio Observatory, California Institute of Technology, MC 249-17, Pasadena, CA 91125, USA  \nReceived -; accepted -', 'ABSTRACT': 'Context. The Serkowski relation is the cornerstone of studies of starlight polarizations as a function of wavelength. Although empirical, its extensive use since its inception by Serkowski et al. (1975) to describe polarization induced by interstellar dust has elevated the relation to the status of an indisputable \'law": this is the benchmark against which models of interstellar dust grains are validated. Aims. We revisit the e ff ects of the three - dimensional structure of the interstellar medium (ISM) on the wavelength dependence of  \ninterstellar polarization.  \nMethods. We use analytical models to show how the wavelength dependence of both the polarization fraction and direction is a ff ected by the presence of multiple clouds along the line of sight (LOS). We account for recent developments in dust distribution modelling and we utilize an expanded archive of stellar polarization measurements to establish the e ff ect of multiple clouds along the LOS. We highlight concrete examples of stars whose polarization profiles are severely a ff ected by LOS variations of the dust grain and magnetic field properties, and we provide a recipe to accurately fit multiple cloud Serkowski models to such cases.  \nResults. We present, for the first time, compelling observational evidence that the three - dimensional structure of the magnetized ISM often results to the violation of the Serkowski relation. We show that three - dimensional e ff ects impact interstellar cloud parameters derived from Serkowski fits. In particular, the dust size distribution in single - cloud sightlines may di ff er from analyses that ignore 3D e ff ects, with important implications for dust modelling in the Galaxy.  \nConclusions. Our results suggest that multiwavelength stellar polarization measurements o ff er an independent probe of the LOS variations of the magnetic field orientation, and thus constitute a potentially valuable new tool for the 3D cartography of the interstellar medium. Finally, we caution that, unless 3D e ff ects are explicitly accounted for, a poor fit to the Serkowski relation does not, by itself, constitute conclusive evidence that a star is intrinsically polarized.  \nKey words. ISM, polarization, etc', '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 et al. 2016). Due to the multifaceted role of interstellar dust in astrophysics and cosmology, a detailed understanding of its nature is imperative.  \nThe 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 sub-mm 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 (farIR, and sub-mm) polarized thermal radiation. Dust polarization provides one of the most reliable and widely employed methods for tracing the properties of the ISM magnetic fields (e.g., Ska-  \nlidis & Tassis 2021; Skalidis et al. 2021, 2022; Pattle et al. 2023; Pelgrims et al. 2023, 2024).  \nDust-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 wavelength λ through the formula  \nP ( λ ) = Pmax · exp GLYPH<18> -K · ln 2 λ max GLYPH<19> , (1)  \nλ  \nwith the use of three free parameters: Pmax the maximum polarization fraction observed at λ max , and K that 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.  \nSince 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 some 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), while in emis-  \nsion, the degree of polarization is relatively flat across frequencies (Hensley & Draine 2023).  \nThe Serkowski relation applies for single-cloud lines of sight (LOSs), and to date, whenever it is applied, it is implicitly assumed that the observed (integrated) P ( λ ) signal is induced by a single dominant polarization screen (cloud) along the LOS. Reality, however, is more complicated. G aia has opened up the third dimension of the sky (depth, Gaia Collaboration et al. 2018), and since then it has become clear that in the majority of the LOSs there are more than one clouds with comparable contributions in the total dust column (e.g., Green et al. 2019; Edenhofer et al. 2024). The variations of the magnetic field along the LOS are imprinted on the individual stellar polarization measurements, allowing the tomographic mapping of the ISM magnetic field (Pelgrims et al. 2024).  \nThe 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 e ff ect 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 non-constant EVPA in scenarios involving multiple clouds with di ff erent magnetic field alignment and / or di ff erent 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 di ff ering 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 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 with simplistic models, finding no firm relationship between their models and contemporary observations.  \nWhile modern researchers occasionally reference these earlier works most of the focus remains on the e ff ect of a two-cloud scenario 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 the P -λ profile 1 or the derived dust parameters.  \nIn this article, we revisit the e ff ects of the three-dimensional structure of the ISM on the wavelength dependence of ISP, accounting for recent developments in dust distribution modelling 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 e ff ects and disentangle the polarizing components within the LOS. Ultimately, we demonstrate that the LOS variations of dust and magnetic field properties a ff ect 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 e ff ects 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 a ff ected 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 in a given wavelength of an intrinsically unpolarized source lying behind NC number of clouds will be  \nqobs = NC X i = 1 qi (2)  \nand similarily for u , with qi being the q induced by the i -th dust cloud in the LOS.  \nIn order to explore this e ff ect, we created analytical models for the cases of 1 and 2 dust clouds along the LOS, for di ff erent 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 a ff ect the cumulative (observed) polarization signal. For our analytical models we used Eq.1 to calculated Pi ( λ ) for 0 . 1 µ m ≤ λ ≤ 1 . 2 µ m , selecting typical values for Pmax 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 ( λ ) = Pi ( λ ) · cos 2 θ i and ui ( λ ) = Pi ( λ ) · sin 2 θ i . For the combined model, we calculated qobs ( λ ), uobs ( λ ) from Eq. 2 and finally ˆ P ( λ ), ˆ θ ( λ ) as  \nˆ P ( λ ) = q q 2 obs ( λ ) + u 2 obs ( λ ) (3a)  \nˆ θ ( λ ) = 1 2 · tan -1 uobs ( λ ) qobs ( λ ) ! , (3b)  \nwhere we used the 2 - argument arctan function.  \nIn Fig. 1 we present some examples of P -λ and θ -λ profiles resulting from the presence of one or two clouds in the LOS, as 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 2 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.  \nAssuming a di ff erence of a factor of two in λ max between our two hypothetical clouds yields a bimodal P ( λ ), which is significantly di ff erent than the shape of a typical Serkowski curve (Fig. 1 upper right panel). Although the di ff erence appears to be significant, it would be challenging to identify it in practice, because multi-band polarization measurements are usually limited to a narrow range of frequencies.  \nModest LOS variations in λ max can result in even more complications. In our assumed case of two clouds, a 40% di ff erence  \nFig. 1: Demonstration of the e ff ect of 3D structure in the Serkowski realtion, 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 corresponds to individual clouds. Left column : The simple case of 1 cloud (Eq. 1). Middle left column : An example with 2 Clouds in the LOS, where both the combined P -λ , and θ -λ seem to follow the Serkowski relation. Middle right column : An example with 2 Clouds in the LOS, where the combined P -λ seems to follow the Serkowski relation, but θ -λ not. Right column : An example with 2 Clouds in the LOS where both P -λ and θ -λ profiles deviate from Serkowski expectations. The parameters for each model are: Left : Pmax = 1%, λ max = 0 . 6 µ m , θ = 30 o . Middle left : Pmax 1 = 1 . 5%, λ max 1 = 0 . 55 µ m , θ 1 = 30 o , Pmax 2 = 0 . 7%, λ max 2 = 0 . 5 µ m , θ 2 = 35 o . Middle right : Pmax 1 = 0 . 7%, λ max 1 = 0 . 5 µ m , θ 1 = 0 o , Pmax 2 = 0 . 8%, λ max 2 = 0 . 7 µ m , θ 2 = 60 o . Right : Pmax 1 = 1 . 3%, λ max 1 = 0 . 4 µ m , θ 1 = 0 o , Pmax 2 = 1 . 4%, λ max 2 = 0 . 8 µ m , θ 2 = 75 o . For all clouds we have used K = 1 . 68 · λ max (Wilking et al. 1980).  \n<!-- image -->  \nin λ max between the two clouds leads to nearly identical to the Serkowski curve profiles (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. In cases where 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 di ff erences in the EVPAs of two clouds lying along the LOS, such as the one in the employed models which is close to 80 degrees, can be revealed through the variation of the EVPA as a function of λ (bottom middle right and right panels in Fig. 1). However, obtaining a su ffi cient number of data points that would allow us to probe such variations is challenging 2 , and for this reason these data are sparse.  \nThese 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 multi-band polarization is 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.', '3.1. Polarimetric data selection - Statistical sample': 'A recent agglomeration of multi-band optical polarization data (Panopoulou et al. 2023), 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 or in 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 single-cloud LOSs are better fit by the Serkowski relation than stars in multiple-cloud LOSs.  \nFor acquisition of archival optical polarization measurements for the statistical comparison, we explored Tables 4, 6 of the compiled polarization catalog of Panopoulou et al. (2023), which provide polarimetry and distances from G aia EDR3 (Bailer-Jones et al. 2021) respectively for ∼ 42 , 000 sources. We excluded stars with distance greater than 1.25 kpc, as this is the limit of our dust map (see Sect. 3.2). As distance, we take the average of the values (\'r\\_med\\_geo\',\'r\\_med\\_photogeo\') provided in the catalog. We excluded measurements in filters with no well-defined e ff ective wavelength ( λ e f f ), therefore, the excluded filters are \'0", \'20", \'23", \'51" which correspond to \'No filter or unclear", \'weightedmean", 0.735-0.804 µ m , 0.5-0.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  \nq mean = P ( qi · w i ) P w i (4)  \nq std = P w i ( xi -¯ x ) 2 GLYPH<16> N -1 N GLYPH<17> · P w i , (5)  \nwith w i = 1 /σ 2 qi , and σ qi the uncertainty of the measurement in q , and similarly for u . We excluded measurements with debiased P /σ P < 2. Additionally, we did not consider measurements with polarization uncertainty σ P < 0 . 1% because they are likely underestimated 3 , as the most reliable polarization standards to date (Blinov et al. 2023) are only accurate down to 0 . 1%. 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. (2023) 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. After all the cuts, the final sample contains 223 sources with measurements in at least three optical bands.', '3.2. Extinction': 'We collected extinction data for our target stars from the Edenhofer et al. (2024) dust map. We choose this map, as it provides good angular resolution (14\'), as well as, 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 ( AV ) based on the conversion table provided by Zhang et al. (2023).  \nFor each of our targets, we examined the "mean" (Edenhofer et al. 2024) cumulative and di ff erential AV vs distance ( D ) profiles by eye, in order to identify the number of significant steps ( NC ) in the profile, which presumably correspond to distinct dust clouds. For targets where more than 3 clouds were identified, or the dust profile was very complex, we set NC > 3. The number of target stars in each subgroup with di ff erent number of identified clouds along the LOS are #( NC = 1) = 49, #( NC = 2) = 74, #( NC = 3) = 70, #( NC > 3) = 30. The targets that correspond to each NC are distributed over di ff erent positions on the sky (Fig. 2).  \nWe 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 that used by Panopoulou & Lenz (2020), who identified the number of clouds in H i data by detecting peaks in the di ff erential 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 norm C = P A i V A max V , with A i V the extinction of a given peak in the LOS, and A max V the extinction of the highest peak in the di ff erential AV vs D profile. However, this was not useful in our study for the following reasons. 1) N norm C  \nis ambiguous. For example, N norm C = 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 di ffi cult 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 a 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 e ffi ciency, magnetic field coherence, and the inclination of the magnetic field.  \nIn our attempts to devise an objective method for 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 and likely introduces some contamination into the subsamples, it is selfconsistent and allows statistical patterns to emerge in a dataset as large as ours. Fig. 3 illustrates examples of dust profiles that we assigned to di ff erent subsamples.', '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, 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.  \nAs a first step we converted Eq. 1 to the q -u space, with the introduction of θ in the model as an additional free parameter. Therefore, our model becomes  \nˆ q = Pmax · exp GLYPH<18> -K · ln 2 λ max λ GLYPH<19> · cos(2 θ ) , (6a)  \nˆ u = Pmax · exp GLYPH<18> -K · ln 2 λ max λ GLYPH<19> · sin(2 θ ) . (6b)  \nThe quantity we wished to minimize for the fits is  \nm = Lq + Lu + Lq Lu + Lu Lq (7)  \nwith  \nLq = 1 N λ · X qi -ˆ qi σ qi ! 2 , (8a)  \nLu = 1 N λ · X ui -ˆ ui σ ui ! 2 , (8b)  \n<!-- image -->  \nAmax  \nFig. 2: Positions of targets depending on their NC . Maps are in Galactic coordinates, in mollweide projection, centered in (0,0). Parallels are drawn every 30 o . Meridians are drawn every 60 o . Longitude increases towards the right-hand side. Colors correspond to the average λ max of each area. Marked areas with 'A' through 'G' contain stars that are located behind a single as well as multiple clouds. See discussion in Sect. 4.3.Fig. 3: Dust profiles for four example targets belonging in four di ff erent subsamples. From left to right: NC = 1, NC = 2, NC = 3, NC > 3. Top row: Cumulative AV vs D . Bottom row: Di ff erential AV vs D .  \n<!-- image -->  \nwhere N λ the number of di ff erent bands where we had measurements for each target, qi , ui the measurement in a given band, ˆ qi , ˆ ui our model as introduced in Eq. 6a, 6b, and σ qi , σ ui the uncertainties of the measurements in q and u respectively. The terms Lq and Lu are the reduced χ 2 metric for q and u respectively. The terms Lq / Lu and Lu / Lq were introduced to ensure a balanced contribution from Lq and Lu . For example, if the objective was  \nto minimize just Lq + Lu , the optimal value V could be achieved with Lq = V and Lu = 0, Lu = V and Lq = 0, or any combination in between. By introducing Lq / Lu and Lu / Lq , we penalize the model when it favors solutions where Lq = 0 or Lu = 0, thus promoting a more balanced fit.  \nFor the modelling procedure we used the O ptuna package in python (Akiba et al. 2019). O ptuna is an automatic hyperpa-  \nrameter 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 e ff ective in multi-parameter optimization with few data points and many parameters, such as the case of the q -u modelling, and more e ffi cient in finding the best solution than traditional algorithms such as Markov-Chain Monte-Carlo (MCMC). We used the default options of O ptuna with n\\_trials = 5000 for our fits. Moreover, we adopted the following priors. For Pmax 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 the mean 0.6 µ m and standard deviation 0.2 µ m , and ensured it was always positive. For K we imposed a Gaussian prior distribution with mean 1 . 68 · λ max and standard deviation 0 . 13 · λ max (Wilking et al. 1980), and ensured it was always positive. For θ we assumed a Gaussian prior distribution with mean the average of the polarization values of measurements in individual bands (derived from the average q and u as in Eq. 4 and converted to θ using Eq. 3b) and 10 o standard deviation .  \nWe evaluated the quality of our fits by measuring the χ 2 per degree of freedom (DoF):  \nχ 2 ν ≡ χ 2 DoF = P GLYPH<18> qi -ˆ qi σ q i GLYPH<19> 2 + P GLYPH<18> ui -ˆ ui σ u i GLYPH<19> 2 2 · N -ν , (9)  \nλ  \nwhere ν = 4 is the number of free parameters in the model ( Pmax , λ 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 used again the q -u fitting technique, following the same logic discussed in Sect. 3.3. The multiple cloud model is derived from Equations 2 and 6, by adding up the q and u components of the individual clouds.  \nWhen 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 O ptuna , we utilized an MCMC algorithm to rapidly explore the parameter space. Given that the two clouds are interchangeable, we imposed the condition λ max 1 > λ max 2 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 O ptuna algorithm, thereby significantly narrowing the search space. Finally, we ran the fit using O ptuna , imposing Gaussian priors on the parameters, based on the mean and the spread of the parameters identified by the MCMC algorithm. We chose O ptuna 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 are evaluated with χ 2 (Eq. 9).  \nν', '4.1. Statistical sample': 'We grouped the sample of sources found in Sect. 3.1 into subsamples depending on the number of clouds ( NC ) lying in their LOS, as derived from the 3D extinction map. We obtained four main subsamples: targets in LOSs with one dominant cloud, NC = 1; with two clouds, NC = 2; three clouds NC = 3; and with more than three clouds, NC > 3. For this sample we fit the traditional Serkowski curve (Eq. 1) for all targets, regardless of NC , to evaluate the capability of the Serkowski relation to fit polarimetric data for subsamples with di ff erent NC . 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 χ 2 ν for the di ff erent subsamples (Fig. 4), and we applied the Kolmogorov-Smirnov (KS) test to identify any statistical discrepancies between the χ 2 ν distributions of the NC = 1 subset and the other three subsets.  \nWe found that the χ 2 ν distribution of the NC = 1 subset is significantly di ff erent from those of subsets with NC = 2 or NC = 3. On the contrary, the χ 2 ν distributions of the NC = 1 subset and of the NC > 3 subset are not discrepant. The p -values of the KS tests for all comparisons are quoted in Table 1, along with the median χ 2 ν for the fits of each subsample. Despite the sparse nature of the data (measurements in only ∼ 3 . 5 di ff erent filters on average per target) the sample is su ffi ciently large to allow 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 multi-cloud LOSs is mainly driven by the behaviour 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 ( NC > 3), there is no statistical di ff erence 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 di ffi cult to detect, especially with limited data points across the spectrum.  \nTable 1: p - values of the KS test between the subsample of targets with NC = 1 and other subsamples. See Sect. 4.1 for details.', '4.2. Two revealing examples': 'To demonstrate observationally the impact of the LOS variations of λ max and EVPA on the Serkowski relation, we employed archival spectropolarimetric data (Bagnulo et al. 2017). We identified two revealing cases, where the multi-wavelength polarization data show clear evidence of the LOS e ff ects.  \nThe 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 e ff ects are significant. Firstly, we  \n0.8  \n0.8  \n0.6  \n0.  \n(a)  \nFig. 4: CDFs of χ 2 for di ff erent subsamples as described in each label.  \n<!-- image -->  \n(c)  \n(b)  \n(d)  \nν  \nobserve that the P -λ deviates from the typical Serkowski curve (Fig. 5 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 behaviour 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. 5). 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, hence its polarization signal carries information from all of them. We fit the data with models with di ff erent number of clouds. Although the extinction profile exhibits three distinct steps, the data can be described equally well with both a 2- and a 3-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  \nthe middle left panel of Fig. 1). We proceed with the simplest case of a 2-cloud model.  \nAs a second example, we show results for the star HD163181, which is an intermediate-size luminous supergiant of spectral class O9.5 (Fig. 6). 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, hence indicating the existence of more than one dominant clouds 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. 6). However, when we considered two clouds along this LOS, as the 3D extinction map suggest, we were able to capture the observed variability of both P and θ with wavelength (red dashed dotted line in Fig. 6 middle panel).  \n0.6  \n0.8  \n0.6  \nNc = 2  \nNc > 1  \nNc = 3  \n8  \n8  \nTable 2: Fitted parameters for the stars HD163181 and HD93222 in Sect. 4.2.  \n<!-- image -->  \n<!-- image -->  \nFig. 5: Example of a case where both the P -λ and the θ -λ profiles deviate from Serkowski expectations. We are able to fit both profiles by introducing a 2-Cloud model. Data from Bagnulo et al. (2017) for HD93222. 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 fitted 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). Dashed-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.  \n<!-- image -->  \nFig. 6: Example of a case where 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 HD163181. Panels and lines similar to Fig. 5. Parameters of the fits are quoted in Table 2  \n<!-- image -->  \n.', '4.3. On 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 e ff ects a ff ect the distributions of λ max , and Pmax that one obtains through the classical fitting of the Serkowski relation in the P -λ space; K  \nlinearly 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  \n1 N λ · X Pi -ˆ Pi σ pi ! 2 , (10)  \nFig. 8: Normalized distributions of Pmax for subsamples with di ff erent number of clouds in the LOS. See Sect. 4.3 for details. The following bins are omitted for clarity. For NC = 2: bin at Pmax / AV ≃ 23 with a height of 0.0135. For NC = 3: bin at Pmax / AV ≃ 13 with a height of 0.0145.  \n<!-- image -->  \nFig. 7: Normalized distributions of λ max for subsamples with di ff erent number of clouds in the LOS. See Sect. 4.3 for details.  \n<!-- image -->  \nassuming a one-cloud model, and imposing the same priors as discussed in Sect. 3.3. Since the e ff ect 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. 7 for λ max , Fig. 8 for Pmax / AV ).  \nPast 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 NC = 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 NC > 1 (multiple-cloud sightlines) are unimodal and peak at around 0.55 µ m, which is close to the value that is considered as the Galactic average. By performing Monte Carlo (MC) simulations we found that the unimodality in the λ max distribution, which is observed for NC > 1, can be reproduced with multiple-cloud models with values of λ max randomly drawn from the NC = 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 e ff ects. Notably, the larger spread of the distribution of λ max for the NC = 2 case is not surprising. Clarke & Al-Roubaie (1984) had theorized that in the case of two clouds in the LOS, the acquired λ max from the Serkowski fit is di ff erent from the average λ max of the two clouds when the two clouds have di ff erent parameters. Furthermore, it is worth considering the following. 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. 2). In most cases, the average λ max di ff ers 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.  \nThe distribution of Pmax / AV is also a ff ected 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 di ff erence between the EVPAs each would impart on its own.When the EVPA di ff erence is close to 0 o the polarization adds up, while when the di ff erence is close to 90 o it cancels out.  \nThis is exactly what we observe in the distributions of Pmax / AV (Fig. 8). The Pmax / AV distribution for NC = 1 is free from LOS e ff ects, hence it should resemble the "ground truth". This distribution peaks at Pmax / AV ≃ 3 -3 . 5% / mag, with only few points with Pmax / AV > 4 . 5% / mag. On the other hand, the NC = 2 distribution peaks in a lower value of Pmax / AV ≃ 2% / mag, with its tail extending to higher values of Pmax / AV than in the NC = 1 case. The shape of the NC = 2 distribution can be understood by the LOS integration e ff ects acting constructively (generating the extended tail), or destructively (generating the lower Pmax / AV peak). The distribution for NC = 3 displays a similar shape, though it is smoother, as the addition of a third cloud facilitates producing polarization closer to the average observed in the NC = 1 case. In the NC > 3 distribution, we generally observe lower values for Pmax / AV compared to the NC > 1 distribution. This is expected, as with more dust components in the LOS, variations in dust and magnetic field properties are more probable, leading to increased depolarization. Overall, both the λ max , and the Pmax / AV are significantly a ff ected by the 3D structure of the ISM, which, if ignored, can lead to erroneous Serkowski best-fit parameters.  \nOur results suggest that fitting a single Serkowski relation is invalid in LOSs with multiple clouds. There can be cases where a single Serkowski relation yields good fits, even when 3D e ff ects 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.', '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 dis-  \nTable 3: Mean λ max values in µ m and their error on the mean for di ff erent regions in the sky and di ff erent subsamples. Regions with N / A in the error of the mean are cases where there is only one star in the specific region for the subsample.  \ntributions 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. 7). 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.  \nWe propose the following two options for constraining the populations of dust in di ff erent 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, di ff erent 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 di ff erent parts of the Galaxy.', '5.2. Implications for the identification of intrinsically polarized stars': 'It is customary to iterpret 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 e ff ectively (Sect. 4.2). This suggests that deviations from the Serkowski curve alone are insu ffi cient to indicate the presence of an intrinsic polarization component. In order to validate that a star not following the Serkowski formula is indeed intrinsically polarized, polarization data could be combined with other evidence, such as detection of significant variability, or comparison with neighboring stars.', '5.3. The possibility of 3D tomography with the Serkowski relation': 'Fitting multiple-cloud models to multiwavelength (optical to nearinfrared) polarization data paves the way to constraining the LOS variations of 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 each individual cloud using the techniques discussed in Sect. 3.4. However, the following considerations should be kept in mind.  \nFirstly, 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 a LOS 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 su ffi ciently 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 direction along the LOS.  \nAnother strategy one could follow to identify cases of clouds with comparable properties is to obtain multi-band or, ideally, spectropolarimetric measurements of multiple stars at di ff erent 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 di ff erent 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, for performing ISM tomography along many LOSs (Pelgrims et al. 2024) in the context of optopolarimetric surveys, such as the P asiphae 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. Multi-band 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 e ff ect of the 3D structure of the ISM on the well-known and widely-used Serkowski relation. Our findings are summarized as follows.  \n- 1. We revisited theoretically how the Serkowski formula could become invalid in LOSs with many dust components. We showed that the variations of both the degree of polarization and EVPA with wavelength can be severely a ff ected 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 be mainly observed in the θ -λ profile.\n- 2. As a test case, we selected two di ff erent targets that appear to be behind multiple dust clouds, based on their extinction  \nprofiles, 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 fitted 2-cloud models.  \n- 3. We proposed the q -u fitting technique as the most appropriate approach to fit the Serkowski relation. This method allows for tracing both P and θ as functions of wavelength simultaneously, enabling fitting for more than one clouds along the LOS.\n- 4. We used a sample of 223 stars behind di ff erent number 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 well the data in the cases of one or more than three clouds in the LOS. The latter is an outcome of averaging.\n- 5. Our analysis of the fitted 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.\n- 6. A poor Serkowski fit by itself does not constitute conclusive evidence of intrinsic polarization in a source.\n- 7. Fitting the Serkowski formula for stars behind multiple dust clouds leads to an incorrect estimation of the fitted parameters, and, thus, of the underlying dust physics.\n- 8. Our results hint towards the possibility of performing magnetic tomography using multiband polarization data in combination with G aia distances.  \nIn conclusion, the Serkowski relation should be cautiously applied when the ISM 3D structure is complex.  \nAcknowledgements. We thank V. Pavlidou for very helpful 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.', 'References': "Akiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. 2019, arXiv e-prints, arXiv:1907.10902  \nAndersson, B. G., Lazarian, A., & Vaillancourt, J. E. 2015, ARA&A, 53, 501 André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 Bagnulo, S., Cox, N. L. J., Cikota, A., et al. 2017, A&A, 608, A146  \nBailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147  \nBastien, P. 2015, in Polarimetry of Stars and Planetary Systems, ed. L. Kolokolova, J. Hough, & A.-C. Levasseur-Regourd, 176  \nBijas, N., Eswaraiah, C., Sandhyarani, P., Jose, J., & Gopinathan, M. 2024, MNRAS, 529, 4234  \nBlinov, D., Maharana, S., Bouzelou, F., et al. 2023, A&A, 677, A144  \nCikota, A., Hoang, T., Taubenberger, S., et al. 2018, A&A, 615, A42  \nClarke, D. 2010, Stellar Polarimetry  \nClarke, D. & Al-Roubaie, A. 1984, MNRAS, 206, 729  \nCodina-Landaberry, S. & Magalhaes, A. M. 1976, A&A, 49, 407  \nCoyne, G. V. & Gehrels, T. 1966, AJ, 71, 382  \nDraine, B. T. 2024, ApJ, 961, 103  \nEdenhofer, G., Zucker, C., Frank, P., et al. 2024, A&A, 685, A82  \nGaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1  \nGinski, C., Garufi, A., Benisty, M., et al. 2024, A&A, 685, A52  \nGreen, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93  \nGuillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16  \nHensley, B. S. & Draine, B. T. 2021, ApJ, 906, 73  \nHensley, B. S. & Draine, B. T. 2023, ApJ, 948, 55  \nHildebrand, R. H. 1983, QJRAS, 24, 267  \nMessinger, D. W., Whittet, D. C. B., & Roberge, W. G. 1997, ApJ, 487, 314 Neha, S., Soam, A., & Maheswar, G. 2024, arXiv e-prints, arXiv:2406.02193  \nMartin, P. G. 1974, ApJ, 187, 461 Martin, P. G. & Whittet, D. C. B. 1990, ApJ, 357, 113 McMillan, R. S. & Tapia, S. 1977, ApJ, 212, 714 Panopoulou, G. V. & Lenz, D. 2020, ApJ, 902, 120  \nPanopoulou, G. V., Markopoulioti, L., Bouzelou, F., et al. 2023, arXiv e-prints, arXiv:2307.05752  \nPatat, F., Maund, J. R., Benetti, S., et al. 2010, A&A, 510, A108  \nPattle, K., Fissel, L., Tahani, M., Liu, T., & Ntormousi, E. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 193 Pelgrims, V., Mandarakas, N., Skalidis, R., et al. 2024, A&A, 684, A162 Pelgrims, V., Panopoulou, G. V., Tassis, K., et al. 2023, A&A, 670, A164 Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 594, A10 Serkowski, K., Mathewson, D. S., & Ford, V. L. 1975, ApJ, 196, 261 Skalidis, R., Sternberg, J., Beattie, J. R., Pavlidou, V., & Tassis, K. 2021, A&A, 656, A118  \nSkalidis, R. & Tassis, K. 2021, A&A, 647, A186  \nSkalidis, R., Tassis, K., Panopoulou, G. V., et al. 2022, A&A, 665, A77  \nTassis, K., Ramaprakash, A. N., Readhead, A. C. S., et al. 2018, arXiv e-prints,  \narXiv:1810.05652  \nTopasna, G. A., Mateja, F. M., & Kaltcheva, N. T. 2023, PASJ, 75, 269  \nTreanor, P. J. 1963, AJ, 68, 185  \nWhittet, D. C. B. 2015, ApJ, 811, 110  \nWhittet, D. C. B. 2022, Dust in the Galactic Environment (Third Edition)  \nWhittet, D. C. B., Gerakines, P. A., Hough, J. H., & Shenoy, S. S. 2001, ApJ, 547,  \n872  \nWhittet, D. C. B., Martin, P. G., Hough, J. H., et al. 1992, ApJ, 386, 562 Wilking, B. A., Lebofsky, M. J., Martin, P. G., Rieke, G. H., & Kemp, J. C. 1980, ApJ, 235, 905  \nZhang, X., Green, G. M., & Rix, H.-W. 2023, MNRAS, 524, 1855"}