* Li (2021) Li, H.-B. 2021, Galaxies, 9, 41
  * Lindegren et al. (2021) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2
  * Magalhães et al. (2005) Magalhães, A. M., Pereyra, A., Melgarejo, R., et al. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 343, Astronomical Polarimetry: Current Status and Future Directions, ed. A. Adamson, C. Aspin, C. Davis, & T. Fujiyoshi, 305
  * Maharana et al. (2020) Maharana, S., Kypriotakis, J. A., Ramaprakash, A. N., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11447, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 114475E
  * Maharana et al. (2021) Maharana, S., Kypriotakis, J. A., Ramaprakash, A. N., et al. 2021, Journal of Astronomical Telescopes, Instruments, and Systems, 7, 014004
  * Marchwinski et al. (2012) Marchwinski, R. C., Pavel, M. D., & Clemens, D. P. 2012, ApJ, 755, 130
  * Martínez-Solaeche et al. (2018) Martínez-Solaeche, G., Karakci, A., & Delabrouille, J. 2018, MNRAS, 476, 1310
  * Miville-Deschênes et al. (2016) Miville-Deschênes, M. A., Duc, P. A., Marleau, F., et al. 2016, A&A, 593, A4
  * Montier et al. (2015) Montier, L., Plaszczynski, S., Levrier, F., et al. 2015, A&A, 574, A135
  * Mouschovias et al. (2006) Mouschovias, T. C., Tassis, K., & Kunz, M. W. 2006, ApJ, 646, 1043
  * Panopoulou et al. (2019a) Panopoulou, G. V., Hensley, B. S., Skalidis, R., Blinov, D., & Tassis, K. 2019a, A&A, 624, L8
  * Panopoulou & Lenz (2020) Panopoulou, G. V. & Lenz, D. 2020, ApJ, 902, 120
  * Panopoulou et al. (2016) Panopoulou, G. V., Psaradaki, I., & Tassis, K. 2016, MNRAS, 462, 1517
  * Panopoulou et al. (2019b) Panopoulou, G. V., Tassis, K., Skalidis, R., et al. 2019b, ApJ, 872, 56
  * Patat et al. (2010) Patat, F., Maund, J. R., Benetti, S., et al. 2010, A&A, 510, A108
  * Pavel et al. (2012) Pavel, M. D., Clemens, D. P., & Pinnick, A. F. 2012, ApJ, 749, 71
  * Pelgrims (2019) Pelgrims, V. 2019, A&A, 622, A145
  * Pelgrims et al. (2021) Pelgrims, V., Clark, S. E., Hensley, B. S., et al. 2021, A&A, 647, A16
  * Pelgrims et al. (2020) Pelgrims, V., Ferrière, K., Boulanger, F., Lallement, R., & Montier, L. 2020, A&A, 636, A17
  * Pelgrims et al. (2022) Pelgrims, V., Ntormousi, E., & Tassis, K. 2022, A&A, 658, A134
  * Pereyra & Magalhães (2004) Pereyra, A. & Magalhães, A. M. 2004, ApJ, 603, 584
  * Planck Collaboration Int. L (2017) Planck Collaboration Int. L. 2017, A&A, 599, A51
  * Planck Collaboration Int. XXXV (2016) Planck Collaboration Int. XXXV. 2016, A&A, 586, A138
  * Planck Collaboration Int. XXXVIII (2016) Planck Collaboration Int. XXXVIII. 2016, A&A, 586, A141
  * Planck Collaboration XI (2020) Planck Collaboration XI. 2020, A&A, 641, A11
  * Planck Collaboration XII (2020) Planck Collaboration XII. 2020, A&A, 641, A12
  * Planck Collaboration XXX (2016) Planck Collaboration XXX. 2016, A&A, 586, A133
  * Ramaprakash et al. (2019) Ramaprakash, A. N., Rajarshi, C. V., Das, H. K., et al. 2019, MNRAS, 485, 2355
  * Rezaei Kh. et al. (2020) Rezaei Kh., S., Bailer-Jones, C. A. L., Soler, J. D., & Zari, E. 2020, A&A, 643, A151
  * Rezaei Kh. & Kainulainen (2022) Rezaei Kh., S. & Kainulainen, J. 2022, ApJ, 930, L22
  * Ritacco et al. (2022) Ritacco, A., Boulanger, F., Guillet, V., et al. 2022, arXiv e-prints, arXiv:2206.07671
  * Robin et al. (2012) Robin, A. C., Luri, X., Reylé, C., et al. 2012, A&A, 543, A100
  * Romero-Shaw et al. (2022) Romero-Shaw, I. M., Thrane, E., & Lasky, P. D. 2022, arXiv e-prints, arXiv:2202.05479
  * Santos et al. (2014) Santos, F. P., Franco, G. A. P., Roman-Lopes, A., Reis, W., & Román-Zúñiga, C. G. 2014, ApJ, 783, 1
  * Skalidis et al. (2018) Skalidis, R., Panopoulou, G. V., Tassis, K., et al. 2018, A&A, 616, A52
  * Skalidis & Pelgrims (2019) Skalidis, R. & Pelgrims, V. 2019, A&A, 631, L11
  * Skalidis et al. (2021a) Skalidis, R., Sternberg, J., Beattie, J. R., Pavlidou, V., & Tassis, K. 2021a, A&A, 656, A118
  * Skalidis & Tassis (2021) Skalidis, R. & Tassis, K. 2021, A&A, 647, A186
  * Skalidis et al. (2021b) Skalidis, R., Tassis, K., Panopoulou, G. V., et al. 2021b, arXiv e-prints, arXiv:2110.11878
  * Skilling (2004) Skilling, J. 2004, in American Institute of Physics Conference Series, Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering, ed. R. Fischer, R. Preuss, & U. V. Toussaint, 395–405
  * Soler et al. (2016) Soler, J. D., Alves, F., Boulanger, F., et al. 2016, A&A, 596, A93
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132
  * Sugitani et al. (2011) Sugitani, K., Nakamura, F., Watanabe, M., et al. 2011, ApJ, 734, 63
  * Tassis & Pavlidou (2015) Tassis, K. & Pavlidou, V. 2015, MNRAS, 451, L90
  * Tassis et al. (2018) Tassis, K., Ramaprakash, A. N., Readhead, A. C. S., et al. 2018, arXiv e-prints, arXiv:1810.05652
  * Zucker et al. (2022) Zucker, C., Goodman, A. A., Alves, J., et al. 2022, Nature, 601, 334
  * Zucker et al. (2019) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, ApJ, 879, 125

##  Appendix A Mock starlight polarization data

To test and validate the inversion method presented in this paper we developed
a pipeline to generate mock samples of starlight polarization data. Our
simulation pipeline has the main following characteristics: ($i$) it makes use
of star samples extracted from the Gaia data to guarantee realistic
distributions of star distances, number density, brightness, etc., ($ii$) it
relies on a self-consistent implementation of the thin-layer dust-cloud model
discussed in Sect. 2 and presented below which includes a prescription for the
intrinsic scatter, and finally ($iii$) it computes and propagates to the
simulated stellar polarization data realistic uncertainties as expected for
the Pasiphae survey. These three parts are described in the rest of this
Appendix. We further explore our toy-model to simulate the polarization signal
from the magnetized ISM in Appendix B.

Our simulation pipeline works as follows. First, we choose a sample of stars
in a small-aperture beam and extract parallaxes, parallax uncertainties and
apparent magnitudes. Second, we simulate the stellar polarization measurements
for any chosen setup of the magnetized and dusty ISM and estimate the
polarization uncertainties. We ‘observe’ the stars by randomly drawing the
Stokes parameters and parallaxes within their respective uncertainties.

###  A.1 Star samples from the Gaia Universe Model Snapshot

We seek to create realistic samples of stellar photometry and distance, based
on the expected sky footprint of the Pasiphae survey. In this section we
discuss how we generate and choose representative samples of stars from the
GUMS database.

Figure 15: Distribution of stellar distances (left) and apparent magnitudes
(right) after the selection criteria have been applied. The distributions for
the large ($0.5^{\circ}$ radius) and small ($0.25^{\circ}$ radius) beams are
shown with black and orange lines, respectively.
Figure 16: Stellar properties and selection criteria. Left: Mean SDSS $r$-band
magnitude in bins defined on the parallax SNR vs. parallax plane. Bins with no
stars are shown in white. Right: 2D distribution of the number of stars in the
$r$ vs. $\varpi$ plane. The cut at ${\rm{SNR}}_{\varpi}=5$ is shown with a
dashed gray line in both panels. Note that the introduction of the parallax
SNR cut imposes a biased selection in the $r-\varpi$ plane. The data shown are
for the large (30-arcmin radius) beam. We do not show the full dynamic range
of $\varpi$ and ${\rm{SNR}}_{\varpi}$ for better visualization of the key
features of the dataset.

####  A.1.1 Mock stellar catalog generation

To obtain realistic samples of stellar parameters (photometry, distances) for
the high Galactic latitude sky, we use the GUMS database (Robin et al. 2012)
associated with Gaia Early Data Release 3 (EDR3, Lindegren et al. 2021), which
provides some photometric information as well as parallaxes. We begin by
selecting sightlines toward each of the HEALPix $N_{\rm{side}}=4$ pixels at
high Galactic latitude (with $|b|>60^{\circ}$). We then query the Gaia archive
within a circular region centered on each of the aforementioned pixels, with
two different radii of 0.5∘ and 0.25∘. For ease of computation, we exclude
very faint stars with ${\rm{G}}>18.5$ within the query. This brightness cut is
later superseded by a stricter cut based on our derived $r$ band photometry,
as explained next. We thus obtain 24 samples of stars for each selected beam
size.

Next, we wish to compute realistic uncertainties for the parallax and
polarization of each star. The GUMS query returns the $V-I$ color, the
$V$-band absolute magnitude $M_{V}$, the extinction $A_{V}$, the apparent
magnitude in the $G$-band, the Gaia broadband color $G_{BP}$ and the
barycentric distance, $d$, in parsec. We determine the parallax errors as a
function of $G$-band magnitude based on EDR3 performance (Lindegren et al.
2021). We use the quoted median parallax error per bin of $G$ of their table 4
(5-parameter solution sources) and create an interpolating function. Each star
is thus assigned a ‘true’ parallax by inverting the distance and converting to
units of milli-arcseconds (mas) and an uncertainty based on its $G$ band
magnitude. The estimation of uncertainties in the polarization measurements
from Pasiphae relies on (a) SDSS $r$-band photometry and (b) the expected
performance of the survey. We derive $r$-band magnitudes in the following
section and briefly discuss the dependence of the uncertainty in the stellar
Stokes parameters on $r$-band magnitude in Appendix A.4.

The stellar properties of the samples cover a broad range of apparent
magnitudes, parallaxes and SNR in parallax (Fig. 16, left, shows a
representative sample of size 1∘). We tailor the stellar samples to the
expected limiting brightness of the Pasiphae survey, by applying a cut of
$r<16\,\rm mag$. The total number of stars per sample (post-cut) varies in the
range [225, 414] and [64, 118] for the 1∘ and 0.5∘ beams, respectively. For
ease of computation, we choose one representative sample for each beam size on
which to apply our analysis. The representative sample for each beam is that
which has a total number of stars (post-brightness-cut) close to the average
number of stars of the ensemble of sightlines. We further limit the samples to
stars with a SNR in parallax higher than 5.

As shown in Fig. 15 (left), the distribution of stellar parallaxes (distances)
peaks at $\sim 1\,{\rm{mas}}$ (1000 pc), but this is brightness-dependent (see
Fig. 16, right). As one would expect, our brightness cut limits the number of
stars at small parallaxes (large distances). We also note there is a less
intuitive selection bias that results from our choice to use stars with
$\varpi/\sigma_{\varpi}>5$. Fig. 16 shows the line of
$\varpi/\sigma_{\varpi}=5$ in two different planes (SNRϖ vs. $\varpi$ and $r$
vs. $\varpi$). As seen in the right panel, imposing the SNR cut results in a
parallax-dependent brightness selection: stars to the left and top of the
dashed line are excluded from our analysis. Our samples miss some faint and
nearby stars. This is the result of the dependence of the parallax error on
stellar brightness. In other words, the limiting magnitude of star samples
analyzed in this work is parallax-dependent.

###  A.2 Derivation of stellar $r$-band photometry given GUMS outputs

We must determine the apparent magnitude in the SDSS-$r$ band. This requires
some manipulation of different color transformations of the photometric data
provided by GUMS. We begin by computing the apparent magnitude in the $V$
band:

$m_{V}=M_{V}+5\,\log(d)-5+A_{V}$ (21)

We can connect the SDSS photometry with Gaia $G$ band photometry using the
relations from Jordi et al. (2010), as initially referenced by the GUMS paper
(Robin et al. 2012):

$\displaystyle G-G_{BP}=$ $\displaystyle-0.1703-1.0813\,(r-i)$
$\displaystyle\phantom{o}-0.1424\,(r-i)^{2}+0.0271\,(r-i)^{3},$ (22)

where the uncertainty of the coefficients is $\sigma=0.1$. We solve this
equation for the color $r-i$. We can then use the color transformations by
Jordi et al. (2006):

$\displaystyle(i-I)=(0.247\pm 0.003)\,(R-I)+(0.329\pm 0.002)$ (23)
$\displaystyle(r-i)=(1.007\pm 0.005)\,(R-I)-(0.236\pm 0.003)\;.$ (24)

Solving the second expression for $(R-I)$ and substituting for $(R-I)$ in the
first expression, we obtain:

$(i-I)=0.0578+0.245\,(r-i)+0.329\;.$ (25)

Using the colour $(r-I)$, we substitute $(i-I)$ from the previous expression
and we solve for $m_{r}$:

$\displaystyle m_{r}-m_{I}=(r-I)=(r-i)+(i-I)$ (26) $\displaystyle
m_{r}=m_{I}+(r-i)+0.0578+0.245\,(r-i)+0.329.$ (27)

We finally obtain an expression for SDSS $r$-band apparent magnitude:

$m_{r}=m_{I}+1.245\,(r-i)+0.3868,$ (28)

where $m_{I}=m_{V}-(V-I)$ can be computed from the $V-I$ color provided in
GUMS.

###  A.3 ISM polarization signal: a toy model

Both the 3D orientation of the magnetic field and the dust density
distribution are expected to vary within a cloud, at least to some extent, due
to compressible MHD turbulence. These two effects induce a spread in the
polarization signal of stars in the background of a cloud (the intrinsic
scatter). This is implemented by considering the signal part to be made of two
components: a regular component, and a stochastic component. To build
intuition on the effect of intrinsic scatter on the Stokes parameters, we
simplify the problem and only implement the scatter that arises from
fluctuations in the 3D geometry of the magnetic field. We defer a treatment of
fluctuations in the dust distribution to future work. Hence, for a given
cloud, a single value of the maximum degree of polarization ($P_{\rm{max}}$)
is used to model the polarization signal of all background stars (Eq. 1). On
the contrary, both the inclination angle and the position angle of the
magnetic field are considered to vary from star to star because of
fluctuations in magnetic field geometry.

We model the total 3D magnetic field as the sum of a regular component
($\mathbf{B}_{\rm{reg}}$) and a stochastic one ($\mathbf{B}_{\rm{sto}}$). We
consider $\mathbf{B}_{\rm{reg}}$ to be uniform within a cloud and to have a
norm of unity. We model $\mathbf{B}_{\rm{sto}}$ through 3D Gaussian
realizations of white noise. That is, the stochastic component is built from
realizations of 3D isotropic random vectors obtained from sampling an
independent normal distribution for each of its three components. Assuming a
large sample of random realizations for $\mathbf{B}_{\rm{sto}}$ we estimate
the sample $rms$ and use it to normalized all the random draws. The
distributions of the norms of $\mathbf{B}_{\rm{sto}}$ have a standard
deviation of one. To reach a statistically stable normalization of
$\mathbf{B}_{\rm{sto}}$, the $rms$ evaluation should be performed for a
sufficiently large sample of realizations. We use at least 1000 realizations.
Finally, we model the 3D total magnetic field as the sum of a regular
($\mathbf{B}_{\rm{reg}}$) and a stochastic component:

$\mathbf{B}_{\rm{tot}}\propto\mathbf{B}_{\rm{reg}}+A_{\rm{turb}}\,\mathbf{B}_{\rm{sto}}\;,$
(29)

where $A_{\rm{turb}}$ quantifies the amplitude of the stochastic component
with respect to the regular one. A different realization of
$\mathbf{B}_{\rm{sto}}$ is attributed to each star in the sample. As a result,
to each star corresponds a different inclination and position angle.
Therefore, by virtue of Eq. (1), different values of the Stokes parameters are
obtained only from fluctuations in the 3D geometry of the magnetic field. We
emphasize that the parameter $A_{\rm{turb}}$ used above should not be confused
with the turbulent-to-mean magnetic field ratio used in methods to estimate
the strength of the magnetic field (e.g., Skalidis & Tassis 2021).
$A_{\rm{turb}}$ is a metric of the statistical fluctuations of the magnetic
field geometry in 3D and is model dependent. We find in Appendix B that values
in the range 0.1 to 0.3 may be representative of clouds at intermediate and
high Galactic latitudes.

To summarize, our toy model has five free parameters per cloud: the cloud
parallax ($\varpi_{\rm{C}}=1/d_{\rm{C}})$, the maximum degree of polarization
($P_{\rm{max}}$), the inclination ($\gamma_{\mathbf{B}_{\rm{reg}}}$) and
position ($\psi_{\mathbf{B}_{\rm{reg}}}$) angles of $\mathbf{B}_{\rm{reg}}$
and, finally, the relative amplitude of fluctuations in magnetic field
orientation ($A_{\rm{turb}}$). We notice that, apart from the cloud parallax,
these parameters are not the same as the model parameters entering our data
equation (Eq. 10) and, therefore, are not the parameters being sampled in the
maximization process.

Due to projection effects, and since the inclination angle is positively
defined, the mean inclination angle of the total magnetic field may generally
be larger than the inclination angle of the regular magnetic field given as
input. In general this results in a depolarization of the cloud signal which
may be stronger in one of the two polarization channels (Stokes $q_{\rm{V}}$
or $u_{\rm{V}}$) depending on the position angle $\psi_{\mathbf{B}}$.
Therefore, in addition to producing a scatter in the polarization plane, the
stochastic component in the 3D magnetic field geometry induces a bias in the
mean Stokes parameters. We further explore our toy model in Appendix B with a
particular emphasis on the bias and the covariance that magnetic-field
fluctuations produce. Here, we note that the bias is not physical and simply
originates from the construction of our model. What matters, and what should
be modeled from real observations, are the mean values of the Stokes
parameters of a cloud and the dispersion around these means.

###  A.4 Realistic noise

Once the polarization signal from the input ISM setup is attributed to each
star in the sample, we add noise in both polarization and parallax. We take
the values of parallax uncertainties from our samples of stars derived from
the Gaia data (Appendix A.1) and the star parallaxes are simply randomized
within their uncertainty range following Gaussian distributions. To give
realistic uncertainties on the Stokes parameters of starlight polarization, we
rely on current expectations of the performance of the WALOP-N instrument to
be used for the Pasiphae survey in the northern hemisphere. The observational
uncertainties, which can be estimated from the optical modeling of the WALOP-N
instrument (Anche et al. 2022), depend primarily on the magnitude of the stars
in the observation band (SDSS-$r$ band) and on the exposure time. To take into
account instrumental systematics, we add to our observational uncertainty
budget a contribution of 0.1% to the two relative Stokes values $q_{\rm{V}}$
and $u_{\rm{V}}$ (see Maharana et al. 2020; Maharana et al. 2021). Therefore,
for each star, an estimate of the observational uncertainty is obtained and
used to randomize the Stokes parameters obtained from the ISM configuration.
The randomization of $q_{\rm{V}}$’s and $u_{\rm{V}}$’s is performed
independently. The uncertainty in the Stokes parameters as a function of star
magnitude is shown in Fig. 17 for two typical exposure times.

Stellar Magnitude in SDSS-$r$ Band [mag]

Figure 17: Uncertainties on individual measurements of the relative Stokes
parameters in SDSS-r band as a function of star’s magnitude as expected from
Pasiphae’s northern instruments (WALOP-N) for 5 minutes and 15 minutes
exposure times shown with solid blue and dashed purple lines, respectively.

Examples of mock starlight polarization data obtained for a five minutes
exposure time, for a beam of $0.5^{\circ}$ circular diameter, and a median
number density are shown in Fig. 1 for the cases of one and two clouds along
the LOS.

##  Appendix B Polarization observables and intrinsic scatter

For the purpose of this work we have developed and implemented a toy model,
described in Sect. A.3, that allows us to generate mock observations based on
the multi-layer magnetized ISM paradigm described in Sect. 2. In particular,
our implementation makes it possible to account for a source of intrinsic
scatter in the simulated data. We expect such scatter to exist because of
small-scale variations in the ISM, that generically result from turbulence.
Generally fluctuations are expected in both density distribution and magnetic
field geometry but dust grain properties, dust temperature, etc. could also
vary on small scales. In our modeling, and as described in Sect. A.3, we only
consider fluctuations in 3D geometry of the magnetic field and thus ignore
fluctuations in dust density distribution or other properties that may affect
the value of the maximum degree of polarization $P_{\rm{max}}$ (see Eq. 1). In
this appendix we explore our toy-model implementation to gain intuition on the
effects of the intrinsic scatter on the Stokes parameters for different levels
of intrinsic scatter (controlled through the value of $A_{\rm{turb}}$ – see
Eq. 29) and for several combinations of the inclination angle
($\gamma_{\mathbf{B}_{\rm{reg}}}$) and position angle
($\psi_{\mathbf{B}_{\rm{reg}}}$) of the regular component of the magnetic
field. In the adopted polarization convention (IAU), $q_{\rm{V}}$ is maximum
for $\psi_{\mathbf{B}_{\rm{reg}}}=0^{\circ}$ and zero at
$\psi_{\mathbf{B}_{\rm{reg}}}=45^{\circ}$ where $u_{\rm{V}}$ is maximum.

We begin by considering the relative Stokes parameters normalized by the value
of $P_{\rm{max}}$. Therefore, in the absence of a stochastic component, the
degree of polarization varies from 0 to 1 for the case of magnetic field
pointing toward the observer ($\gamma_{\mathbf{B}_{\rm{reg}}}=90^{\circ}$) to
the case where the magnetic field lies in the POS
($\gamma_{\mathbf{B}_{\rm{reg}}}=0^{\circ}$). From geometrical considerations,
we expect that the effect of the addition of a stochastic component on the
values of the Stokes parameters $q_{\rm{V}}$, $u_{\rm{V}}$ will depend on the
chosen value of $A_{\rm{turb}}$, on $\gamma_{\mathbf{B}_{\rm{reg}}}$, and
possibly on $\psi_{\mathbf{B}_{\rm{reg}}}$. To infer such dependence, we want
to consider different orientations of the regular magnetic field as seen from
the observer.

We consider the following setup. We place an observer above the pole of the
northern hemisphere ($b\geq 0^{\circ}$) of a HEALPix map with
$N_{\rm{side}}=2$ (Górski et al. 2005). This map contains 28 pixels to which
correspond 3D orientations with co-latitude and longitude coordinates. The
former corresponds to the inclination of the magnetic field line with respect
to the POS and latter to the position angle in the POS. By construction, these
orientations are symmetric in longitude and therefore only half of the points
do show different Stokes parameters given the definition of the polarization;
$\psi_{\mathbf{B}}$ is a two-circular quantity defined in the range
$[0^{\circ},\,180^{\circ})$. Hence, we drop half of the northern hemisphere
and thus consider 14 combinations of
$(\gamma_{\mathbf{B}_{\rm{reg}}},\,\psi_{\mathbf{B}_{\rm{reg}}})$. To each,
corresponds a pair of $(q_{\rm{V}},\,u_{\rm{V}})$ at position
$\cos^{2}\gamma_{\mathbf{B_{\rm{reg}}}}\;(\cos[2\psi_{\mathbf{B_{\rm{reg}}}}],\,\sin[2\psi_{\mathbf{B_{\rm{reg}}}}])$
in the polarization plane (see Eq. 1). Those positions are shown as red
crosses in Fig. 18. These crosses correspond to the Stokes parameters that
would be observed for all background stars if the magnetic field was made of
the regular component only. Then, we add the stochastic component in magnetic
field geometry according to Eq. 29. This creates variations about those
values. Example of variations obtained in the polarization plane for
$A_{\rm{trub}}=0.1$ is shown in Fig. 18.

Figure 18: Examples of scatter produced in the $(q_{\rm{V}},\,u_{\rm{V}})$
plane by the addition of a stochastic component to a regular magnetic field
for different inclination and position angle of the latter and for
$A_{\rm{turb}}=0.1$. Fourteen configurations are shown, as explained in the
text. The red crosses indicate the $(q_{\rm{V}},\,u_{\rm{V}})$ obtained for
$A_{\rm{turb}}=0$ (regular field only). The scatter around each cross shows a
2D histogram of 10,000 realizations of $(q_{\rm{V}},\,u_{\rm{V}})$ values
obtained with the addition of the stochastic component. We note that the
$(q_{\rm{V}},\,u_{\rm{V}})$ are normalized by $P_{\rm{max}}$.

For most of the 14 configurations, the 2D scatter produced in the polarization
plane cannot be considered as resulting from two independent normal
distributions; the correlation coefficient is not zero in general. Instead,
configurations with the regular magnetic field close to the POS are
significantly asymmetric. Interestingly, depending on the position angle, the
correlation coefficient can take large values. For example, in the case of
$\gamma=0^{\circ}$ and $\psi=22.5^{\circ}$, the correlation coefficient is
close to $-1$. The more the regular magnetic field points to the observer, the
more the scatter becomes symmetric in $q_{\rm{V}}$ and $u_{\rm{V}}$.

At this stage, we notice that fluctuations in dust density distribution, or
any other ISM parameters affecting the value of $P_{\rm{max}}$, would add
radial scatter on this plot. These would then reduce the possible asymmetry
and therefore would reduce the value of the correlation coefficient. From the
shape of the 2D scatter seen in Fig. 18 we may conclude that modeling the
effect of the intrinsic scatter through the characterization of a bivariate
normal distribution with given covariance matrix is a fair approach which
should be general enough to cover a large class of models of sources of
fluctuations in the magnetized ISM. This motivates our choice while writing
the model equation in the main text (see Sect. 2). From this point, we
therefore consider that the polarization signal in the
$(q_{\rm{V}},\,u_{\rm{V}})$ plane, mean and intrinsic scatter, can be self-
consistently described by a bivariate normal distribution centered on the mean
polarization and with a covariance matrix potentially having a non-zero off-
diagonal element.

A scatter in the $(q_{\rm{V}},\,u_{\rm{V}})$ plane corresponds to scatter in
the degree of polarization and in the EVPA. The scatter in EVPA is thought to
inform on the degree of turbulence in clouds and to be related to the
amplitude of the magnetic field (e.g., Skalidis & Tassis 2021), at least when
the magnetic field lies mostly in the POS.

Figure 19: Scatter in EVPA as a function of $A_{\rm{turb}}$. For each value of
$A_{\rm{turb}}$, $10^{5}$ realizations of $\mathbf{B}_{\rm{tot}}$ are
computed. The Stokes parameters are computed and polarization position angles
determined and compared to the polarization vector corresponding to the
regular component only. The standard deviation of the polarization angle
difference is then computed and plotted against $A_{\rm{turb}}$. We carry the
analysis for $\gamma_{\mathbf{B}_{\rm{reg}}}=0^{\circ},30^{\circ}$ and
$60^{\circ}$.

We explore the dependence of the scatter in EVPA as a function of the
parameter $A_{\rm{turb}}$. Because of projection effects, we carry this
analysis for three values of $\gamma_{\mathbf{B}_{\rm{reg}}}$
($0^{\circ},\,30^{\circ},\,60^{\circ}$). Our results are reported in Fig. 19.
The EVPA scatter increases linearly at low values of $A_{\rm{turb}}$, then the
increase slows down until the EVPA scatter reaches saturation. The larger the
angle of inclination, the greater the rate of increase and thus, breakup and
saturation occur at lower values of $A_{\rm{turb}}$. Given that there is
currently a lack of starlight polarization data for diffuse sightlines, there
is no much observational constraints on the scatter in EVPA and, hence, on the
value that $A_{\rm{turb}}$ may take. In small angular regions toward denser
clouds or denser sightlines, authors report values of the EVPA scatter to
range from few degrees up to about fifteen degrees (Soler et al. 2016; Planck
Collaboration Int. XXXV 2016; Panopoulou et al. 2016; Skalidis et al. 2021b),
also consistent with numerical simulations of the ISM in sub/trans-Alvénic
regimes (e.g., Skalidis et al. 2021a). We will assume that those values are
representative to our case. Therefore, we choose values of $A_{\rm{turb}}$ in
the range 0.1 to 0.3 to be representative of degree of intrinsic scatter we
may find at intermediate and high Galactic latitude, in general, from our
future survey.

In Fig. 20, we show distributions of the Stokes parameter $q_{\rm{V}}$
(normalized to $P_{\rm{max}}$) as a function of $A_{\rm{turb}}$ and for the
same three inclination angle values studied before. We show the mean
(continuous lines) and the interval between percentiles 16 and 84 (shaded
areas) of the distributions. We see that depending on the inclination angle
there might be a non-zero difference between the mean of the distributions and
the values corresponding to the regular field only (dashed lines). This is a
bias that purely results from projection effects. This bias is larger for
small inclination angles (i.e., magnetic field close to the POS) than for
large inclination angles and vanishes when the field lines point toward the
observer. The same figure is obtained for Stokes $u_{\rm{V}}$ but with a
position angle rotated by 45∘. This result is consistent with what Hu &
Lazarian (2022) obtained from numerical simulations of incompressible MHD
turbulence in molecular cloud.

Figure 20: Distributions of $q_{\rm{V}}/P_{\rm{max}}$ corresponding to various
inclination angle and amplitude of the intrinsic scatter. We carry the
analysis out for $\gamma_{\mathbf{B}_{\rm{reg}}}=0^{\circ},\,30^{\circ}$, and
$60^{\circ}$ and $A_{\rm{turb}}\in[0,\,1]$. For each value pair, $10^{5}$
realizations of $\mathbf{B}_{\rm{tot}}$ are computed and $q_{\rm{V}}$ is
evaluated. The continuous lines show the mean of the distribution for
$\gamma_{\mathbf{B}_{\rm{reg}}}=0^{\circ},\,30^{\circ},\,\text{and}\,60^{\circ}$
in blue, orange and green, respectively. The shaded areas indicate the range
span between percentiles 16 and 84 of the distribution. The dashed lines
indicate the $q_{\rm{V}}/P_{\rm{max}}$ values obtained when no stochastic term
is added in the magnetic field.

In the following we explore the effects of the intrinsic scatter on the
element of the covariance matrix of the $(q_{\rm{V}},\,u_{\rm{V}})$ pairs that
it generates. In particular, we are interested in exploring the dependence on
the inclination and position angle of the regular magnetic field and on the
dependence with the amplitude of the stochastic component relative to the
regular one.

In Fig. 21, we present the evolution of $C_{\rm{int},qq}$ and
$C_{\rm{int},qu}$ as a function of $\gamma_{\mathbf{B}_{\rm{reg}}}$, for the
three position angles
$\psi_{\mathbf{B}_{\rm{reg}}}=0^{\circ},\,22.5^{\circ},\,\text{and}\,45^{\circ}$
and for $A_{\rm{turb}}=0.2$. We do not show $C_{{\rm{int}},uu}$, as it is
identical to $C_{{\rm{int}},qq}$ but for a position angle rotated by
$45^{\circ}$.

Figure 21: Dependence of $C_{{\rm{int}},qq}$ (blue) and $C_{{\rm{int}},qu}$
(purple) as a function of $\gamma_{\mathbf{B}_{\rm{reg}}}$ and for three
values of the position angles: $\psi_{\mathbf{B}_{\rm{reg}}}=0^{\circ}$
(dashed lines), $\psi_{\mathbf{B}_{\rm{reg}}}=22.5^{\circ}$ (dashed-dotted
lines), and $\psi_{\mathbf{B}_{\rm{reg}}}=45^{\circ}$ (dotted lines). A value
of $A_{\rm{turb}}=0.2$ is chosen here. The elements of the covariance matrix
are given in units of $P_{\rm{max}}^{2}$. A value of 0.01 corresponds to a
polarization uncertainty ($\sigma_{q}$ or $\sigma_{q}$) of 10 per cent of
$P_{\rm{max}}$.

We see that the covariance-matrix elements are non-trivial but certainly
trigonometric functions of the inclination angle. The cross term vanishes when
$\psi_{\mathbf{B}_{\rm{reg}}}=0^{\circ}$ or any multiple of $45^{\circ}$
(i.e., when either $q_{\rm{V}}$ or $u_{\rm{V}}$ is zero). If fluctuations in
magnetic field geometry were the only sources of intrinsic scatter, a detailed
analysis of the covariance-matrix elements from observational data points
could potentially lead to constraints on inclination angle of the magnetic
field with respect to the POS. This information is generally not accessible
from dust-related polarization observables and providing a direct access to it
would constitute a breakthrough in dust studies of the magnetized ISM.
However, as we have already noticed, other sources of fluctuations are
expected in the turbulent ISM and they directly affect the values of
$P_{\rm{max}}$, adding a certainly non-negligible contribution to
$C_{{\rm{int}},qq}$ and $C_{{\rm{int}},uu}$ and reducing the relative
importance of $C_{{\rm{int}},qu}$, thus bringing some complexity to the
situation depicted in Fig. 21. While future work will be necessary to clarify
whether those sources of scatter can be disentangled, we note that other
sources of complication will come from both observational noise and from the
a-priori unknown value of $A_{\rm{turb}}$.

Finally, we address the question of the significance of the intrinsic scatter
and how it compares to observational uncertainties for a large range of
$A_{\rm{turb}}$ values and for different inclination and position angles of
the regular component of the magnetic field. For this investigation we
consider a sample of 100 stars in the background of a cloud. We fix our
observational uncertainty on the polarization as
$\sigma_{p}=\sigma_{q}=\sigma_{u}=0.2\%$ and we fix a degree of polarization
of the cloud to be $p_{\rm{C}}=0.6\%$ in absence of scatter. We thus change
the value of $P_{\rm{max}}$ to compensate for change of
$\gamma_{\mathbf{B}_{\rm{reg}}}$.

For each set of ($A_{\rm{turb}}$, $\gamma_{\mathbf{B}_{\rm{reg}}}$,
$\psi_{\mathbf{B}_{\rm{reg}}}$) values we can generate simulated Stokes
parameters (i) computed only from $\mathbf{B}_{\rm{reg}}$, (ii) computed with
$\mathbf{B}_{\rm{tot}}$ for specified $A_{\rm{turb}}$ and (iii) with the
addition of the observational noise. For each simulated dataset we estimate
the total covariance matrix of the Stokes parameters and the contributions
from the observational and from the intrinsic scatter taken separately.
Denoting $C_{{\rm{tot}},xy}=C_{{\rm{obs}},xy}+C_{{\rm{int}},xy}$ with $x,\,y$
being either $q$ or $u$, we can measure the different terms in our simulations
and thus compare the relative amplitude of the different contributions. We
perform this exercise for 10,000 simulated datasets to build distributions of
the elements of the covariance matrix. The result is shown in Fig. 22 in
which, for each matrix element taken separately and for each of their
contributions (noise and intrinsic scatter), we show how the distributions
change as a function of $A_{\rm{turb}}$ by reporting the median and the
interval between 16 and 84 percentiles as shaded areas.

Figure 22: Comparison of the contributions from the observational noise and
from the intrinsic scatter to the different elements of the covariance matrix
as a function of the amplitude of the stochastic component and for different
configurations of $\mathbf{B}_{\rm{reg}}$. We show the total elements of the
covariance matrix (solid lines) along with the contribution from the noise
(dashed lines) and from the intrinsic scatter (dotted lines). Blue, green and
orange correspond to $C_{{\rm{int}},qq}$, $C_{{\rm{int}},uu}$ and
$C_{{\rm{int}},qu}$ respectively. The lines and shaded areas show the median
and 16 and 84 percentile of the values obtained while repeating the analysis
for 10,000 simulations in which both the noise and the intrinsic scatter vary.
In this simulation we have set $p_{\rm{reg}}=0.6\%$ and
$\sigma_{p}=\sigma_{q}=\sigma_{u}=0.2\%$ corresponding roughly to a SNR of 3
for the individual star polarization measurement.

The relative importance of the contribution from the intrinsic scatter with
respect to the contribution from the noise to the total scatter in the data
points is well demonstrated by the different panels of Fig. 22. Depending on
the exact setup of the ISM and amplitude of the stochastic component in the
magnetic field, the data scatter can be dominated either by the intrinsic
scatter or by the noise. Keeping in mind that this also depends on the SNR
level which we have fixed to $\approx 3$ for individual measurements, we see
that the contribution from the intrinsic scatter exceeds the noise
contribution in at least one of the polarization channels when
$A_{\rm{turb}}\gtrsim 0.3$ and that this threshold reduces when
$\gamma_{\mathbf{B}_{\rm{reg}}}$ increases. When either or both the amplitude
of the cloud polarization gets lower or the overall uncertainty on individual
star measurements gets larger, this threshold goes to lower $A_{\rm{turb}}$
values (not shown in the figure).