The Solar System through dust maps. Encounter statistics with interstellar gas

Introduction

During its orbit around the center of the Milky Way, the Solar System crosses the interstellar medium (ISM), encountering environments characterized by varying gas and dust densities, elemental composition, temperatures, ionization states, and velocities (see e.g., [51, 17, 34]). The ISM traversed by the Sun has received significant scientific attention for several reasons. For example, the density of the surrounding ISM and its relative velocity with respect to the Sun affect the size and structure of the heliosphere (see e.g., [73, 50, 47]), which in turn plays a key role in modulating the Galactic cosmic-ray flux (see e.g., [72, 65]). The encountered ISM also gives rise to an influx of interstellar dust into the Solar System, as first detected by the Ulysses spacecraft (see e.g., [26, 25, 40]) and later confirmed by Galileo and Cassini (see e.g., [5, 2]). Evidence that such material can enter the Earth’s atmosphere was presented by Taylor et al. [68]. Encounters with different ISM environments along the Sun’s Galactic orbit have also been linked to anomalies in geological records (see e.g., [35, 37, 71]), possible gravitational influence on the Oort cloud (see e.g., [31, 32, 69]), and, more tentatively, Earth’s climate (see e.g., [62, 30, 55, 54]).

At present, the Solar System resides within the Local Bubble, a large ($R \sim100$–$200\ \mathrm{pc}$), hot ($T \sim10^{6}\ \mathrm{K}$), and low-density ($n \sim0.01\ \mathrm{cm}^{-3}$) cavity of the ISM (see e.g., Cox & Reynolds 1987; Galeazzi et al. [20]; Breitschwerdt et al. [10]; Zucker et al. [76]; O’Neill et al. [53]). More specifically, the Sun is currently crossing a complex of diffuse ($n_{\mathrm{H}} \sim0.1 - 0.2\ \mathrm{cm}^{-3}$), warm ($T \sim5{,}000 - 10{,}000\ \mathrm{K}$), and partially ionized clouds that are located inside the Local Bubble. These clouds, known as the Cluster of Local Interstellar Clouds (CLIC; e.g., Lallement & Bertin [42]; Redfield & Linsky [59]; Vannier et al. [70]), are thought to have formed from a supernova occurring within the pre-evacuated Local Bubble (Zucker et al. [77]), are associated to an outflow originating from the Scorpius-Centaurus OB association (Piecka et al. [57]), and have also been linked to the recent deposition of $^{60}\mathrm{Fe}$ detected in Antarctic snow and ice-cores (Koll et al. [37, 38]).

In its recent ($\lesssim20\ \mathrm{Myr}$) Galactic orbit, the Sun has not always resided within the Local Bubble. Approximately 14 Myr ago, while still outside this cavity, whose expansion started around the same epoch (see e.g., Zucker et al. [76]), the Sun was traversing the Radcliffe Wave, a coherent $\sim2.7\ \mathrm{kpc}$-long massive gaseous structure (Alves et al. [3]), in the Orion region (Maconi et al. [45]). The Solar System may subsequently have crossed the Orion–Eridanus superbubble (see e.g., Schulreich et al. [61]) before entering the Local Bubble about 5-6 Myr ago (see e.g., Zucker et al. [76]). Since then, it has traveled through the Local Bubble interior, possibly crossing the Lynx cloud (Opher et al. [54]), and entered the CLIC between $\sim240{,}000 - 160{,}000\ \mathrm{yr}$ ago (Vannier et al. [70]), reaching its current location.

Reconstructing the Solar trajectory further back in time becomes increasingly uncertain due to potential gravitational scattering and the difficulty to reconstruct past gas structures. An alternative approach is to statistically characterize the environments encountered by the Sun, for example by estimating the frequency of crossings at a given density. Different studies have addressed this question using either analytical estimates applied to observational constraints (see e.g., Talbot & Newman [67]) or numerical simulations (see e.g., Kokaia & Davies [36]).

In this paper, we aim to characterize the ISM crossed by the Sun using the three-dimensional (3D) dust maps of the Solar neighborhood by Leike et al. [44] and Edenhofer et al. [16]. These maps, together with other reconstructions available in the literature, were made possible by the European Space Agency’s Gaia mission (Gaia Collaboration et al. [18]) together with the advent of deep wide-field photometric surveys (e.g., Pan-STARRS1, Chambers et al. 2019), and the advancement of statistical techniques (see e.g., Edenhofer et al. [14]). Collectively, these reconstructions represent the first detailed 3D view of the Solar neighborhood, revealing previously unknown Galactic structures (see e.g., Lallement et al. [41]; Alves et al. [3]; Kormann et al. [39]).

Our approach consists of rigidly translating the maps along a circular orbit and sampling their density fields along the trajectory of a point particle moving on simulated Solar orbits. From the resulting time series, we estimate the frequency and duration of passages through different ISM density regimes, the fraction of time spent above different density thresholds, and the associated gas and dust exposure. The main advantage of this approach is that the sampled density field is based on observationally informed 3D dust reconstructions of the local ISM. Its main limitation is that the maps provide a present day snapshot of a constantly evolving ISM, cover a local portion of the Galaxy, and are not very sensitive to gas densities above $n_{\mathrm{H}} \sim100 - 200\ \mathrm{cm}^{-3}$.

This paper is organized as follows. In Sect. 2, we describe the dust maps used for this study; in Sect. 3, we present the integration method and maps sampling procedure, define the encounter statistics and influx estimates; and in Sect. 4, we present and discuss our results. Finally, we summarize our results in Sect. 5.

Data

Dust maps

We use the 3D interstellar dust extinction maps of the Solar neighborhood by Leike et al. [44], hereafter L20) and Edenhofer et al. [16], hereafter E24) to characterize the ISM traversed by the Solar System along its Galactic orbit. The two maps are not combined into a single density field, but are analyzed independently throughout this work. We chose these maps because they currently offer the best distance resolution (see also the discussion in Zucker et al. [78]), which is essential for resolving high-density regions. Additionally, we used the two maps as they differ in spatial extent, resolution, and spatial discretizations, allowing us to test how these quantities affect the results.

Leike et al. (2020)

The L20 dust map is based on the stellar parameter catalog by Anders et al. [4], which combines observational data from Gaia DR2, AllWISE, Pan-STARRS, and 2MASS. The map is based on a Cartesian grid centered on the Sun, with a voxel resolution of 1 pc and a reported functional resolution of 2 pc. The grid extends 740 pc in the $xy$-plane and 540 pc in the $z$-direction.

The L20 map is provided as a set of 12 posterior samples of the extinction density. We used all the 12 samples for our analysis. The map is expressed in terms of differential extinction, given as the Gaia G-band optical depth, $\tau_G$, per parsec. To convert the differential extinction to gas volume density, we follow the approach of Zucker et al. [75] and provide further details in Appendix A.

Edenhofer et al. (2024)

The E24 dust map is based on the catalog of stellar distance and extinction estimates by Zhang et al. [74], hereafter ZGR23), which is primarily based on Gaia BP/RP spectra from Gaia DR3 (Gaia Collaboration et al. [19]). E24 modeled the 3D dust extinction within 1.25 kpc of the Sun in the HEALPix (Górski et al. [23]) spherical coordinate system, with 516 logarithmically spaced distance bins from 69 pc to 1250 pc, using the NIFTy statistical library (Edenhofer et al. [15, 14]). The map has an angular resolution of 13.7′ ($N_{\mathrm{side}} = 256$) and a distance resolution ranging from 0.4 pc at 69 pc to 1.25 kpc. Within 180 pc from the Sun the distance resolution is 1 pc or better, and remains below 2 pc out to 360 pc. However, these values describe only the binning of the map, rather than its effective spatial resolution. The effective resolution remains on parsec scales, but depends on distance and on the number of stars constraining the reconstruction in each region. The E24 map excludes the innermost 69 pc around the Sun, where extremely low extinction and previous reconstructions (Leike et al. [44]) indicate a negligible dust content.

E24 reconstructed 12 samples of the inferred 3D dust extinction distribution, along with the posterior mean. In our analysis, we used the 12 samples. The reconstruction is provided in terms of the dimensionless extinction defined in ZGR23. As before, we refer to Appendix A for the conversion of the map to gas volume density.

The authors also provide the E24 map interpolated onto a Cartesian grid, together with the interpolation script. However, we perform our analysis on the original HEALPix grid. This choice avoids additional smoothing of high-density regions introduced by the interpolation. Moreover, the Cartesian interpolation does not fully remove the effect of the original distance-dependent bin sizes, since the interpolated grid still retains an imprint of the underlying HEALPix discretization.

For both maps, the conversion procedure described in Appendix A yields the total volume density of hydrogen nuclei, independent of ISM phase ($n = n_{\mathrm{H\,I}} + 2n_{\mathrm{H}_{2}}$). The L20 and E24 maps are each subject to different sources of uncertainty, which propagate into our crossing statistics. We refer to the respective papers for a detailed discussion. Additional limitations relevant to our analysis include the potential underestimation of extinction toward dense clouds, due to the lack of stellar data at high densities to constrain the 3D reconstruction [see also 75] and, for the E24 map, the increasing volume bin size with increasing distance. Additionally, converting dust extinction into gas volume density introduces an additional uncertainty, as the dust-to-gas conversion factor is not universal and may vary with the local ISM conditions (see e.g., O’Neill et al. [52]). These limitations will impact our results, as later discussed.

In Figure 1, the projections of the two maps onto the $xy$-plane are shown for one representative posterior sample each. The L20 map captures the local ISM within 370 pc of the Sun, and local structures such as the Local Bubble are resolved in details. However, due to its limited extent, it cannot unveil the large-scale organization of the ISM, including kpc-long structures such as the Radcliffe Wave (Alves et al. [3]) and other superclouds (Kormann et al. [39]), which instead become recognizable in the E24 map. In the same figure, we also show the probability density functions of the gas density of the maps. They are broadly consistent over most of the range relevant for the encounter statistics analyzed here. Differences appear mainly at very low densities, $n_{\mathrm{H}} \lesssim10^{-2}\,\mathrm{cm}^{-3}$, where the reconstructions are more uncertain.

Methods

Integration of the orbits

To derive the crossing statistics and gas flux, we integrated the orbit of the Solar System together with that of the map. At each time step, the map is rigidly translated along a circular orbit at the Sun’s current Galactocentric radius, with its orientation fixed toward the Galactic center, in order to preserve the orientation of the Galactic structures. This is not intended to model the full dynamical evolution of the ISM, but rather to account for the mean rotation of the local Galactic disk. In Appendix B, we also tested the case of a random map orientation. The setup described above is illustrated schematically in Figure 2. Throughout this work, when referring to the Sun or the Solar System in the context of the orbit integration, we mean a test particle initialized with the given initial conditions and integrated in a Milky Way-like gravitational potential.

For the orbits integration, we used the Galactic dynamics package galpy (Bovy [8]). We adopted the MWPotential2014 model for the Milky Way potential which includes a bulge, disk, and dark-matter halo component (see Bovy [8], for details). The velocity of the LSR was set to $v_{\mathrm{LSR}} = 220\ \mathrm{km}\,\mathrm{s}^{-1}$ (see e.g., Bovy et al. [9]). The map follows a circular orbit at $R = R_{\odot}$ in the Galactic midplane, with a velocity relative to the Galactic center of $(U,V,W)=(0,220,0)\ \mathrm{km}\,\mathrm{s}^{-1}$.

In the reference integration, we positioned the Sun at the Galactocentric radius of $R_{\odot}=8.122\ \mathrm{kpc}$ (GRAVITY Collaboration et al. [24]) and at a height above the disk $z_{\odot}=20.8\ \mathrm{pc}$ (Bennett & Bovy [7]). The solar velocity relative to the LSR was set to $(U_{\odot},V_{\odot},W_{\odot})=(11.1,12.24,7.25)\ \mathrm{km}\,\mathrm{s}^{-1}$ (Schönrich et al. [60]).

To study how the ISM structure and the vertical motion of the Sun affect the crossing statistics, we also performed a set of controlled experiments. In these tests, we modified the initial $z$ and $v_z$ values of the Sun to sample different maximum heights, or modified the ISM density field within a given height range by randomly shuffling the density values among the map cells. This preserves the density distribution, but removes the spatial correlations. These experiments are described in Sect. 4.2. Additional tests using different initial conditions are presented in Appendix B.

In all cases, the integrations were performed such that the map is sampled by the Sun for at least 2 Gyr.¹ We adopted a timestep of 0.001 Myr, ensuring that the corresponding spatial sampling is always smaller than the minimum maps resolution of 0.4 pc. This results in a typical spatial step of $\sim0.2\ \mathrm{pc}$ in the Galactocentric reference frame and of $\sim0.05\ \mathrm{pc}$ relative to the maps.

To account for the statistical uncertainties in the Sun’s position and velocity, we repeated the orbit integration 100 times, each with different initial conditions drawn via Monte Carlo sampling of the associated uncertainties, reported in the referenced literature. Statistical uncertainties in the dust maps were addressed by applying this procedure to each of the 12 samples. In all cases, we considered the center of the maps to be positioned at the Galactocentric height of 20.8 pc, assuming the Sun’s height by Bennett & Bovy [7] as reference.

Sampling of the dust maps

We tracked the positions of the Sun and the map throughout the integration. At each time step for which the Sun is located inside the map boundaries, we recorded its position and velocity relative to the map together with the gas density. This procedure produces a time series of sampled densities, an example of which is shown in Figure 2. We then considered density thresholds ($n_{\mathrm{th}}$) between 1 and $500\ \mathrm{cm}^{-3}$ and derived three main encounter statistics, described in §3.2.1, and gas flux and mass accretion estimates, as outlined in §3.2.2.

Encounter statistics

For each density threshold, we identified contiguous time intervals during which the encountered density satisfies $n \ge n_{\mathrm{th}}$. A crossing is defined as an event beginning when the encountered density rises above the threshold and ending when it drops below the threshold again. We estimated the following encounter statistics: crossing frequency, $f_{\mathrm{cross}}(\ge n_{\mathrm{th}})$, crossing duration, $\tau_{\mathrm{cross}}(\ge n_{\mathrm{th}})$, and residence time fraction, $\eta_{\mathrm{res}}(\ge n_{\mathrm{th}})$. Similar encounter statistics have been used in previous studies (see e.g., Talbot & Newman [67]; Kokaia & Davies [36]; Pfalzner et al. [56]), although their precise definitions are not always stated explicitly. We define the quantities used in this work as follows.

The crossing frequency is computed to quantify how often the Solar System encounters regions above a given density. For each threshold, we counted the number of crossings, $N$, above a density threshold. The crossing frequency is then defined as

where $N(\ge n_{\mathrm{th}})$ is the number of crossings above the threshold and $T_{\mathrm{inside}}$ is the total time for which the Solar orbit is inside the map volume. For each dust map sample, we combined the 100 orbit realizations by summing both the number of crossings and the corresponding sampled time. This gives one long-time estimate of the crossing frequency for each dust map sample. We then reported the median over the 12 dust map samples, together with the 16th and 84th percentiles.

The crossing duration quantifies the typical time spent in a single contiguous encounter above a given density threshold. We measured the duration of each individual crossing, $\Delta t_{\mathrm{cross},i}(\ge n_{\mathrm{th}})$, across the dust map samples and orbit realizations. The characteristic crossing duration, $\tau_{\mathrm{cross}}(\ge n_{\mathrm{th}})$, is defined as the median of this distribution, and the 16th and 84th percentiles are used to describe the spread of individual crossing durations.

The residence time fraction quantifies the fraction of time that the Sun spends in regions whose density exceeds a given value. It is defined as

where $T_{\mathrm{inside}}(\ge n_{\mathrm{th}})$ is the total time spent above the threshold. Following the same procedure as for the crossing frequency, we report the median and the 16th to 84th percentile range over the 12 dust map samples.

These quantities, $f_{\mathrm{cross}}$, $\tau_{\mathrm{cross}}$, and $\eta_{\mathrm{res}}$, are not independent. To first order, they are related by

They therefore provide complementary descriptions of the encounter history and are sensitive to different properties of the density field and of the orbit, as discussed in Sect. 4.

We also computed the volume filling factor of the gas,

where $V(\ge n_{\mathrm{th}})$ is the volume occupied by gas above the threshold and $V_{\mathrm{tot}}$ is the total spatial volume considered. When studying the vertical structure of the maps, we computed the same quantity in bins of $\lvert z\rvert$. The values of $f_{\mathrm{Vol}}$ shown later in the figures correspond to the median over the 12 dust map samples.

Gas flux and interstellar dust accretion estimates

Our setup allows us to estimate, to first order, the time-averaged gas flux encountered by the heliosphere, $\Phi_{\mathrm{gas}}^{\mathrm{tot}}$ [cm$^{-2}$ · yr$^{-1}$], as

where $n(t)$ and $v_{\mathrm{rel}}(t)$ are the gas density and relative velocity between the Sun and the map at the integration time $t$, respectively. Since the dust maps are static, $v_{\mathrm{rel}}(t)$ only accounts for the relative motion between the Solar orbit and the translated map, and does not include intrinsic gas motions within the ISM. We also computed the relative flux contribution of different gas density $\Phi_{\mathrm{gas}}(\ge n_{\mathrm{th}})$ to the total flux.

As an order of magnitude calculation, we convert the gas flux into a simple estimate of the interstellar dust mass accretion rate onto Earth, $\dot{M}_{\mathrm{dust}}$ [kg · yr$^{-1}$]. We approximate this rate as

where $\rho_{\mathrm{dust}}(t)$ is the dust mass density, $\sigma_{\mathrm{eff}}(t)$ the effective cross section of the Earth, and $\epsilon(t)$ a dimensionless efficiency factor. This factor accounts for processes that regulate the delivery of interstellar dust to Earth, including gravitational focusing, heliosphere filtering, and grain-size- and charge-dependent transport of interstellar dust. The modeling of these processes is beyond the scope of this work. For a more in depth discussion, we refer the reader to Slavin et al. [63], Krüger et al. [40], Godenkov & Izmodenov (2024), and references therein.

Equation (6) can be simplified by additionally assuming that the cross section and filtering factor do not depend on time, and that the ISM is well mixed. From this, the dust mass density can be derived as $\rho_{\mathrm{dust}}(t) = D\,\mu_{\mathrm{H}}\,m_{\mathrm{H}}\,n(t)$, with $D$ being the dust-to-gas mass ratio, $\mu_{\mathrm{H}}$ the mass contribution from helium and metals per hydrogen nucleus, and $m_{\mathrm{H}}$ the hydrogen mass. Equation (6) finally reduces to

which is the first order approximation we used in this work. For the parameters, we assumed $D = 0.01$, $\mu_{\mathrm{H}} = 1.4$, $m_{\mathrm{H}} = 1.67 \cdot10^{-27}$ kg, and $\sigma_{\mathrm{eff}} = \sigma_{\oplus}$, with $\sigma_{\oplus}$ being the geometrical cross section of the Earth.

Results and Discussion

Encounter statistics

The Solar System encounters a variety of Galactic environments along its orbit, spanning a broad range of gas densities. In Figure 2, we show an example of the density crossed by the Sun along its trajectory obtained by using the set-up described in Sect. 3.1. The trajectory is characterized by intermittent passages through higher-density regions (typically of order $\sim10 - 10^{2}$ cm$^{-3}$), separated by longer intervals spent in lower-density environments.

We quantified this encounter history using the statistics defined in Sect. 3.2.1. The results are shown in Figure 3 for the reference integration of both the L20 and E24 maps. In Table 1, we report the corresponding central values at representative density thresholds.

$n$ [cm$^{-3}$]	$f_{\mathrm{cross}}$ [Gyr$^{-1}$] L20	$f_{\mathrm{cross}}$ [Gyr$^{-1}$] E24	$\tau_{\mathrm{cross}}$ [kyr] L20	$\tau_{\mathrm{cross}}$ [kyr] E24	$\eta_{\mathrm{res}}$ L20	$\eta_{\mathrm{res}}$ E24
$\ge1$	1198	627	64	142	0.18	0.47
$\ge3$	540	539	59	138	0.06	0.26
$\ge10$	156	299	54	138	0.01	0.1
$\ge30$	36	112	51	135	0.002	0.03
$\ge100$	5	15	46	124	0.0002	0.002
$\ge200$	1	2	44	113	0.0001	0.0003
$\ge300$	0.51	0.6	44	110	2e-05	0.0001

We find that the crossing frequency, $f_{\mathrm{cross}}$, decreases smoothly with increasing density threshold. For densities higher than $n \ge1$ cm$^{-3}$, the inferred crossing frequency is greater than 600 Gyr$^{-1}$ for both maps, corresponding to one crossing every $\sim0.8 - 1.6$ Myr. For $n \ge10$ cm$^{-3}$, the inferred crossing frequency is $\sim155$ Gyr$^{-1}$ using L20 and $\sim300$ Gyr$^{-1}$ using E24, corresponding to one such crossing every $\sim6$ Myr for L20 and $\sim3$ Myr for E24, respectively. Crossings of very dense regions are rare: at $n \sim300$ cm$^{-3}$, the inferred frequency decreases to $\sim0.5$ Gyr$^{-1}$. Nevertheless, both maps predict at least one crossing per Gyr for thresholds around 200 cm$^{-3}$.

A direct comparison with previous estimates of $f_{\mathrm{cross}}$ in the literature is not straightforward, because earlier studies generally count encounters with assumed cloud populations rather than defining contiguous crossings above a fixed density threshold, as done here. In an early estimate based on a simplified population of interstellar clouds, Talbot & Newman [67] found that the Sun should encounter around 31 Gyr$^{-1}$ clouds with $n_{\mathrm{H}} \gtrsim100$ cm$^{-3}$ and about 4 Gyr$^{-1}$ clouds with $n_{\mathrm{H}} \gtrsim10^{3}$ cm$^{-3}$. More recent estimates have focused mainly on encounters with giant molecular clouds (GMCs). Kokaia & Davies [36] modeled the Galactic distribution of GMCs using a spiral arm pattern based on Hou & Han [29] and found an average encounter rate of 1.6 $\pm$ 1.3 Gyr$^{-1}$ for the Sun. Similarly, Moreira et al. [49] found that Sun-like particles encounter a GMC approximately every 700 Myr. These values are comparable to our crossing frequencies only if the relevant GMC densities correspond to thresholds of order $n_{\mathrm{H}} \sim100$ cm$^{-3}$. For this reason, we interpret the agreement with GMC-based studies as qualitative rather than one-to-one.

The estimated median duration of individual crossings, $\tau_{\mathrm{cross}}$, is of the order of several tens to several hundreds of kyr over a broad range of density thresholds. For the L20 map, the median duration of crossings above 1 cm$^{-3}$ is $\sim64^{+140}_{-40}$ kyr; this decreases to $\sim54^{+80}_{-34}$ kyr for $n \geq10$ cm$^{-3}$. At higher density, $\tau_{\mathrm{cross}}$ becomes nearly constant, with a median values around 45 kyr. For the E24 map, the median duration of crossings above 1 cm$^{-3}$ is $\sim142^{+530}_{-100}$ kyr. At higher densities, the median duration remains approximately constant and reaches $\sim124^{+140}_{-90}$ kyr at 100 cm$^{-3}$, with little variation for higher thresholds.

This flattening should not be interpreted as a physical saturation of cloud sizes, but rather as a consequence of the finite spatial resolution of the maps. In principle, crossings above higher density thresholds are expected to become shorter, because denser structures occupy progressively smaller volumes. In the reconstructed L20 and E24 density fields, however, the highest density structures are not fully resolved. For the L20 map, which is defined on a uniform Cartesian grid, it is significant that the inferred crossing duration approaches the characteristic time required for the Sun to cross a 1 pc resolution element at the average relative velocity between the Sun and the map, as shown in Figure 3. The E24 map is additionally affected by its native spatial discretization. E24 uses a set of concentric spherical shells, logarithmically spaced in radius and with a HEALPix pixelization (see Sect. 2). The volume of the spatial bins therefore increases with distance from the map center, which affects the inferred crossing duration. Figure 4 illustrates this effect by

For the residence time fraction, $\eta_{\rm res}$, the experiment with the L20 map shows that the Sun spends about 18% of its time in regions with densities above 1 cm$^{-3}$, and about 1% of its time above 10 cm$^{-3}$. The residence fraction decreases further at higher thresholds, consistent with the fact that dense gas occupies a progressively smaller volume fraction. For E24, the corresponding residence fractions are larger: about 47% of the time is spent in regions above 1 cm$^{-3}$, and about 10% of its time above 10 cm$^{-3}$. These values are also affected by the bin volume limitation described above. These results imply that in the reconstructed density field the Sun spends most of its time in low density environments. In particular, using L20 as reference because of its uniform spatial grid, the Sun remains below 10 cm$^{-3}$ for approximately 99% of the sampled time, and below 1 cm$^{-3}$ for approximately 80% of the time. Although this is a first order estimate, it is consistent with the presence of several hundredparsec-sized voids in the solar neighborhood, as observed from the top down view of the reconstructed maps (see Figure 1), and as similarly seen in nearby galaxies (see e.g., Lee et al. [43]; Barnes et al. [6]).

A related estimate for $\eta_{\rm res}$ was obtained by Pfalzner et al. [56], who studied the fraction of time spent by interstellar objects inside GMCs using the simulation set-up of Kokaia & Davies [36]. They found that typical interstellar objects in the Solar neighborhood, with velocities of order $30\ \mathrm{km\ s^{-1}}$, spend only about $0.1$ to $0.2\%$ of their time inside GMCs, while slower ISOs, with velocities of order $5\ \mathrm{km\ s^{-1}}$, spend about $1$ to $2\%$ of their time inside GMCs. Although these values are not directly comparable to our residence fractions, because they refer to ISO trajectories through an idealized GMC population rather than to the Solar orbit through reconstructed density fields, they provide a useful reference point for interpreting residence fractions in dense environments.

In general, the derived crossing statistics are subject to several limitations. First, high density structures are not fully captured by the maps (see Sect. 2). However, because such regions occupy a small volume fraction, their contribution is expected to be limited, although they may dominate the statistics of rare high density encounters. Second, the results are affected by the native spatial discretization of the maps, in particular for E24, as discussed above. Third, the effective resolution and smoothing of the maps can dilute unresolved dense structures below their peak densities, leading to an underestimate of the crossing statistics at high density thresholds. Smoothing can also broaden density structures, affecting the inferred crossing durations and residence fractions at lower thresholds. Finally, the intrinsic motions of the ISM are not captured by these maps. Gas motions could modify the true duration of an encounter by changing the relative velocity between the Sun and the gas, either shortening or lengthening individual crossings depending on the direction of the motion. We therefore interpret these quantities as estimates derived from the reconstructed density fields, rather than as direct measurements of the true physical ISM values.

Despite these limitations, the crossing statistics provide useful insight into the encounter history of the Solar System. It is notable that the two maps yield broadly similar trends despite their different spatial extents. The L20 map reconstructs the nearby $\sim400\ \mathrm{pc}$ region around the Sun, whereas E24 extends to $\sim1.25\ \mathrm{kpc}$ and includes larger Galactic structures (see e.g., Kormann et al. [39]). This similarity suggests that the central values of the encounter statistics are primarily sensitive to the density distribution sampled by the orbit, rather than to the detailed large scale organization of the ISM. This interpretation is also supported by a test in which we repeated the calculations using a random orientation of the map, which gives similar results (Appendix B). We investigate this point further in Sect. 4.2, where we study how the ISM structure and the Solar orbit affect these quantities. In Appendix B, we also show that the encounter statistics are robust to changes in the Solar initial conditions and that restricting E24 to the L20 volume improves the agreement between the two maps.

Vertical ISM stratification and the solar orbit

In this section, we investigate how the encounter statistics are affected by the vertical stratification of the ISM and by the Solar orbit, focusing in particular on the Sun’s vertical motion. In details, we aim to identify the density regimes and ranges of vertical heights sampled by the Sun in which the statistics are mainly controlled by the gas density distribution itself, and those in which the coupling between the ISM structure and the Solar motion becomes important. For this part of the analysis, we use the L20 map, taking advantage of its uniform Cartesian grid.

The different encounter statistics are sensitive to different properties of the gas density and the Sun’s orbit. The crossing frequency and crossing duration depend not only on the gas distribution, which sets the probability of encountering gas above a given threshold, but also on the spatial organization of the gas. A density field in which dense gas is concentrated into a small number of coherent structures produces fewer and longer crossings than one in which the same dense gas fraction is distributed among many smaller structures. Both quantities depend on the relative velocity between the Sun and the map, since larger relative velocities shorten individual crossings and can increase the number of structures crossed per unit time. By contrast, the residence fraction is less sensitive to the detailed spatial organization of the gas, as it measures the total fraction of time spent above a given density threshold, rather than the number or duration of individual crossings. All the quantities are sensitive to the heights sampled by the Sun and the gas stratification of the disk.

To investigate this dependence, we first examine the vertical structure of the ISM in the L20 map. Figure 5 shows this analysis through a set of diagnostic plots. The side projection of the map (upper left panel) shows that the gas is organized into structures which appear to be uniformly distributed over a vertical range, roughly within $|z| \sim100 - 130\ \mathrm{pc}$, rather than being concentrated around $z = 0\ \mathrm{pc}$. This behavior is also visible in the normalized volume filling factor, $f_{\rm vol}$, shown in the lower left panel for representative low, intermediate, and high density thresholds ($n_{\rm H} \geq5$, $50$, $500\ \mathrm{cm^{-3}}$). For the lowest threshold, the filling factor is largest near the Galactic midplane and decreases gradually with increasing $|z|$, broadly consistent with the expected vertical stratification of Galactic dust (see e.g., Misiriotis et al. [48]; Guo et al. [27]). At higher density thresholds, however, dense structures appear to populate more uniformly the vertical heights before declining rapidly at heights with $|z| \geq100 - 130\ \mathrm{pc}$. The same trend is reflected in the vertical distribution of the density bins shown in the right panel. This behavior may partially be explained by the limited volume sampled by the L20 map, rather than reflect the global vertical structure of the Galactic ISM, but it might also suggest that different density regimes may have different effective vertical distributions. Low density gas may be primarily governed by gravity, while denser structures may be shaped by feedback processes. It is worth noting that a similar behavior is also found when using the larger E24 map, although in this case the interpretation is affected by the bin size distance dependence of the native grid. We briefly discuss this in Appendix § 9. However, a detailed investigation of the scale heights dependency on density thresholds is beyond the scope of this work.

To test the role of the Solar vertical motion, we modified the initial $z$ and $v_z$ of the test particle, while keeping the remaining phase space coordinates fixed to the reference values. The resulting encounter statistics are shown in Figure 6 for representative density thresholds. This experiment isolates how the statistics change when the orbit samples different portions of the vertically stratified ISM.

Among the three quantities considered, $\tau_{\rm cross}$ shows the clearest dependence on the maximum height sampled by the Sun: it decreases with increasing vertical amplitude for all density thresholds. This trend can be explained by the fact that orbits reaching larger $|z|$ have larger vertical velocities when traversing the disk. As discussed in Sect. 4.1, $\tau_{\rm cross}$ approaches the resolution limit of the map, corresponding to the time required to cross a spatial resolution element. Because of this limited resolution, we cannot investigate in detail whether the sampling angle of the orbit might affect the inferred crossing durations.

The crossing frequency and the residence time show a more complex behavior, which depends on the adopted density threshold. Since these statistics are approximately related by relation 3, their trends are not completely independent. At low density thresholds, $f_{\mathrm{cross}}$ remains approximately constant as long as the Sun samples heights within $|z|\lesssim100$ pc, and then appears to decline at larger vertical amplitudes. Instead, $\eta_{\mathrm{res}}$ decreases more gradually. At intermediate density thresholds, corresponding to regimes where the vertical stratification of the ISM begins to deviate from a simple monotonic decline, both $f_{\mathrm{cross}}$ and $\eta_{\mathrm{res}}$ remain relatively stable while the Sun samples the region within $|z|\lesssim100$–$130$ pc. A clearer decrease appears only when the orbit reaches heights beyond this range. At the highest density thresholds, the statistics become more sensitive to the particular locations of dense structures within the finite volume of the L20 map. In this regime, neither the crossing frequency nor the residence fraction shows a simple monotonic decrease with maximum vertical height.

The results suggest that, within the local volume sampled by the L20 map, the probability of encountering dense ISM structures does not necessarily decrease smoothly with distance from the Galactic midplane, as found in other studies (see e.g., Kokaia & Davies [36]). However, in Kokaia & Davies [36], the vertical decrease follows directly from the assumed distribution of GMCs, and the sampled range of heights extends well beyond the one considered here. When we consider orbits reaching larger vertical amplitudes ($|z|\geq130$ pc), we also find that the encounter probability decreases. Our findings also suggest that the Solar vertical oscillation through the Galactic midplane does not necessarily translate into a periodic perturbation of the Solar System by dense ISM structures, as sometimes invoked in relation to periodicities in the terrestrial impact cratering record (see e.g., Rampino [58]). The result that dense structures do not appear to follow a simple monotonic decline with height may be also relevant for studies of stellar cluster dissolution. For example, Moreira et al. [49] successfully modeled the open cluster scale height using vertically decreasing encounter probabilities with molecular clouds. It would therefore be interesting to explore whether alternative assumptions, for example a more uniform probability within a finite vertical range (e.g., a sigmoid-like vertical dependence), affect the inferred evolution of open cluster populations.

The residence fraction provides the clearest way to isolate the coupling between the vertical structure of the ISM and the Solar orbit. Unlike the crossing frequency and crossing duration, it does not depend directly on how the gas is spatially organized. Instead, it measures only the total fraction of time spent above a density threshold. If this is correct, the residence fraction can be approximated as

where $\rho_{\odot}(z)$ is the fraction of time spent at a given $z$, $f_{\mathrm{vol}}(\geq n_{\mathrm{th}}|z)$ is the volume filling factor of gas above the threshold at $z$, and $z_{\min}$ and $z_{\max}$ are the minimum and maximum heights reached by the Solar orbit. If $f_{\mathrm{vol}}$ is approximately independent of $z$ over the vertical range sampled by the Sun, this expression reduces to

In this limiting case, the residence fraction is determined by the gas distribution within the sampled region, rather than by the detailed vertical coupling between the orbit and the ISM.

We test this interpretation using the fiducial integration setup and two additional experiments. In the first experiment, we randomly shuffle the L20 density field within $|z|\leq130$ pc. This preserves the gas distribution within the vertical range sampled by the Solar orbit, while removing the spatial correlations of the gas. In the second experiment, we use initial conditions for which the Sun reaches larger heights, where the vertical stratification of the disk changes more significantly. The results are presented in Figure 7.

We find that Eq. (8) successfully predicts the residence fraction in all three cases. In the fiducial and shuffled experiments, the residence fraction is also well approximated by the volume filling factor computed for the sampled vertical range, thus satisfying the approximation from Eq. (9). This indicates that, for the present Solar orbit, the $\eta_{\mathrm{res}}$ is mainly controlled by the gas distribution. In the shuffled experiment, however, $\tau_{\mathrm{cross}}$ becomes nearly constant, reflecting the loss of coherent density structures. In the experiment with a larger vertical amplitude, the volume filling factor computed over the sampled region is no longer sufficient to reproduce $\eta_{\mathrm{res}}$. In this case, the coupling between the Solar vertical motion and the vertical stratification of the ISM becomes important.

These tests show that the residence fraction is primarily controlled by the combination of the vertical volume filling factor and the fraction of time spent by the Sun at each height. By contrast, the crossing duration and crossing frequency depend more strongly on the spatial organization of the gas, because they are sensitive to how structures are organized.

Interstellar gas flux and dust exposure

We next estimate the amount of interstellar material encountered by the Solar System along its orbit. This quantity depends not only on the gas density sampled along the trajectory, but also on the relative velocity between the Sun and the surrounding ISM. We therefore compute the time averaged gas flux, $\Phi_{\mathrm{gas}}^{\mathrm{tot}}$, using Eq. 5, and quantify both the contribution from individual density ranges and the corresponding accreted interstellar dust mass. The estimates here presented are derived from the L20 map.

We obtain a time averaged gas flux along the sampled Solar orbit of $\Phi_{\mathrm{gas}}^{\mathrm{tot}} = (6.0 \pm0.9) \cdot10^{13}\ \mathrm{cm}^{-2}\cdot\mathrm{yr}^{-1}$. However, the contribution to this flux is not uniform across density ranges and in Figure 8 we report the differential and cumulative contribution to $\Phi_{\mathrm{gas}}^{\mathrm{tot}}$ by density. Most of the contribution to the flux is coming from relatively low density environments: 50% of the total flux comes from gas with $n_{\mathrm{H}} \lesssim3 - 4\ \mathrm{cm}^{-3}$, while 90% is produced by gas with $n_{\mathrm{H}} \lesssim48\ \mathrm{cm}^{-3}$. Gas in the lowest density bin, $0 < n_{\mathrm{H}} \leq1\ \mathrm{cm}^{-3}$, already accounts for about 20% of the total flux. This contribution should be interpreted with caution as such low density environments are among the most difficult regimes to reconstruct accurately with dust extinction maps. As a simple test, we added a random isotropic velocity component to the gas. Each Cartesian component was drawn from a Gaussian distribution with zero mean and dispersion $\sigma_{\mathrm{ISM}} = 8.9\ \mathrm{km\,s}^{-1}$ ([46]). This increases the total flux slightly, to $\Phi_{\mathrm{gas}}^{\mathrm{tot}} \simeq7.3 \times10^{13}\ \mathrm{cm}^{-2}\,\mathrm{yr}^{-1}$, while leaving the relative contribution of different density ranges essentially unchanged. Since a fully consistent treatment of gas motions would also modify the cloud crossing times, we use the static map calculation as our reference case and interpret this test only as an estimate of the possible effect on adding a random isotropic velocity.

The undisturbed ISM flow outside the heliosphere carries both gas and dust, but it is strongly modified as it interacts with the Solar wind. In particular, the transport of interstellar dust into the Solar System depends also on grain size, charge, and heliospheric conditions (see e.g., [64, 40]). Nevertheless, a fraction of interstellar dust grains can enter the Solar System and the Earth’s atmosphere (see e.g., [26, 68, 1]). We therefore convert the time averaged gas flux into an estimate of the accreted interstellar dust mass under the simplifying assumptions described in Sect. 3.2.2 and using Eq. (7). Using the reported parameters values, we obtain

This value represents the orbit averaged estimate over the sampled Solar trajectory. However, the Solar System is currently moving through a low density ISM environment, known as the CLIC (see § 1 and, e.g., [59]). For comparison with present day conditions, we therefore also compute the dust mass accretion rate using only the flux measured for $n_{\mathrm{H}} < 1\ \mathrm{cm}^{-3}$. In this case, we find

which corresponds to $\dot{M}_{\mathrm{dust}} \approx\epsilon430\ [\mathrm{kg}\cdot\mathrm{yr}^{-1}]$. Matching a present day interstellar dust accretion rate on Earth of order $100\ \mathrm{kg}\,\mathrm{yr}^{-1}$ (see e.g., (Strub et al. 2021)) would therefore require an effective delivery efficiency of $\epsilon\approx0.23$. As discussed in Sect. 3.2.2, $\epsilon$ absorbs unresolved focusing, heliospheric filtering, and grain size and charge dependent transport effects and should therefore be viewed as a very first order estimate.

Finally, we examined how the incident flux and accreted dust mass vary with distance from the Galactic midplane. The incident flux reflects the combined effect of the density distribution and the relative velocity along the orbit. As shown in the left panel of Fig. D.1, the mean instantaneous flux is slightly enhanced near the midplane, where the vertical velocity is larger. The fractional contribution to the total accreted dust mass depends additionally on the amount of time spent by the Sun at each height. As a result, the heights associated with the largest instantaneous flux are not necessarily those that dominate the total accreted mass over the full orbit. This is illustrated in the right panel of Fig. D.1, where the largest contribution comes from heights close to the vertical turning points of the Solar orbit.

Caveats

The encounter statistics presented in this work are subject to several limitations:

• the maps represent local, present day, static conditions ISM;

• the formation, evolution, and dispersal of clouds are neglected;

• the maps have finite spatial resolution and do not resolve compact high density regions (see Sect. 2);

• the spatial grid on which each map is defined can influence the inferred encounter statistics, especially the crossing durations (see Sect. 4.1);

• intrinsic gas motions and transient events, such as SNe, are not included.

Despite these limitations, the analysis provides a useful first order, data driven estimate of the density environments encountered by the Sun and of the corresponding gas and dust exposure along its Galactic orbit. The results should therefore be interpreted as observationally informed estimates based on the reconstructed local ISM, rather than as an accurate reconstruction of the Solar System’s actual encounter history.

Conclusions

In this work, we used 3D dust maps of the Solar neighborhood to characterize the interstellar environments encountered by the Solar System along its Galactic orbit. We integrated the Solar orbit through the reconstructed density fields and measured the frequency, duration, and residence time fraction of passages above different gas density thresholds. We also estimated the corresponding interstellar gas flux and dust exposure. Our main results can be summarized as follows:

• The Solar System samples a broad range of interstellar densities along its orbit, but spends most of its time in low density environments. The crossing frequency decreases rapidly with increasing density threshold. For $n_{\mathrm{H}} \geq1\ \mathrm{cm}^{-3}$, both maps give $f_{\mathrm{cross}} > 600\ \mathrm{Gyr}^{-1}$. For $n_{\mathrm{H}} \geq10\ \mathrm{cm}^{-3}$, the frequency decreases to $\sim155\ \mathrm{Gyr}^{-1}$ in L20 and $\sim300\ \mathrm{Gyr}^{-1}$ in E24. Crossings of very dense gas are rare, with $f_{\mathrm{cross}} \sim0.5\ \mathrm{Gyr}^{-1}$ at $n_{\mathrm{H}} \geq300\ \mathrm{cm}^{-3}$, although both maps predict one crossing per Gyr for thresholds around $200\ \mathrm{cm}^{-3}$. The median duration of individual crossings is typically several tens of kyr, for example $\sim64\ \mathrm{kyr}$ in L20 and $\sim142\ \mathrm{kyr}$ in E24 for $n_{\mathrm{H}} \geq1\ \mathrm{cm}^{-3}$. At higher density thresholds, these durations become resolution limited and should be interpreted as likely upper limits on the true crossing times of unresolved compact structures. For the residence fraction, we found that the Sun spends $\sim18\%$ of the sampled time above $1\ \mathrm{cm}^{-3}$ in L20 and $\sim47\%$ in E24, while the corresponding fractions above $10\ \mathrm{cm}^{-3}$ are about $1\%$ and $10\%$. The numerical differences between L20 and E24 reflect their different spatial discretizations, with E24 being more affected as it is defined on a spherical grid.

• The vertical gas stratification in the maps depends on the density threshold considered. Low density gas is more likely

to be found near the Galactic midplane, with its volume filling factor decreasing gradually at larger $|z|$. Higher density structures, instead, appear to populate more uniformly the vertical heights before declining rapidly at heights with $|z| \ge100 - 130$ pc. Consistently with this, the crossing frequency of dense structure remains relatively stable for Solar orbits confined within this vertical range. This suggests that Galactic midplane passages do not necessarily lead to periodic perturbations of the Solar System.

– The residence time fraction reflects the coupling between the vertical structure of the ISM and the vertical motion of the Sun. We find that it is well reproduced by weighting the volume filling factor of gas above a given density threshold at each height by the fraction of time that the Sun spends at that height. By contrast, the crossing duration and crossing frequency also depend on the spatial organization and coherence of the gas structures along the orbit.

– The gas and dust exposure of the Solar System is dominated by relatively low density environments. We find a time averaged gas flux of $\Phi_{\mathrm{gas}}^{\mathrm{tot}} = (6.0 \pm0.9) \times10^{13}\ \mathrm{cm}^{-2}\ \mathrm{yr}^{-1}$. About half of this flux is contributed by gas with $n_{\mathrm{H}} \lesssim3 - 4\ \mathrm{cm}^{-3}$, and about 90% comes from gas with $n_{\mathrm{H}} \lesssim48\ \mathrm{cm}^{-3}$. We also converted the gas flux into a first order estimate of the dust mass accretion rate. For comparison with present day low density conditions, and under very basic assumptions, we obtain $\dot{M}_{\mathrm{dust}} \approx\epsilon\ 430\ [\mathrm{kg}\cdot\mathrm{yr}^{-1}]$. This corresponds to an effective delivery efficiency of $\epsilon\approx0.23$, to match the present day $100\ \mathrm{kg}\cdot\mathrm{yr}^{-1}$.

These results provide a data driven estimate of the density environments sampled by the Solar System in the local Galactic disk. They should be interpreted as a first order approximation, as they are limited by the finite extent and resolution of the maps, the assumption of a static density field, and the absence of transient events such as nearby supernovae. Nevertheless, to our knowledge, this work provides the first estimate of these quantities based on 3D observational data.

Future work based on upcoming dust-map reconstructions and numerical simulations will allow these statistics to be better constrained. Such studies will also make it possible to characterize additional environmental quantities encountered by the Sun, including temperature and magnetic field strength. These findings could be complementary to heliospheric and geological studies aimed at understanding how changes in the ISM affect the heliosphere and may leave signatures in terrestrial records.

Acknowledgements. EM, JA, and SH were co-funded by the European Union (ERC, ISM-FLOW, 101055318). This work has made use of data from the European Space Agency (ESA) mission Gaia https://www.cosmos.esa.int/gaia, processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

References

$n$ [cm$^{-3}$]	$f_{\mathrm{cross}}$ [Gyr$^{-1}$] L20	$f_{\mathrm{cross}}$ [Gyr$^{-1}$] E24	$\tau_{\mathrm{cross}}$ [kyr] L20	$\tau_{\mathrm{cross}}$ [kyr] E24	$\eta_{\mathrm{res}}$ L20	$\eta_{\mathrm{res}}$ E24
$\ge1$	814	609	61	139	0.1	0.35
$\ge3$	306	422	56	137	0.03	0.17
$\ge10$	73	189	52	137	0.01	0.05
$\ge30$	14	48	48	132	0.001	0.01
$\ge100$	1	3	44	114	0.0001	0.0004
$\ge200$	0.33	0.3	45	107	1e-05	4e-05
$\ge300$	0.12	0.06	45	97	4e-06	1e-05

From differential extinction to gas volume density

The maps by Leike et al. [44] (hereafter L20) and Edenhofer et al. [16] (hereafter E24) are both provided in terms of the differential extinction: L20 as the Gaia G-band optical depth per parsec ($s_{\mathrm{G}} = \frac{\tau_{\mathrm{G}}}{1\ \mathrm{pc}}$) and E24 as the dimensionless extinction ($A'_{\mathrm{ZGR23}} = \frac{\mathrm{d}A_{\mathrm{ZGR23}}}{1\ \mathrm{pc}}$) defined in Zhang et al. [74] (hereafter ZGR23). To convert the differential extinction data to gas volume density, we follow the procedure of Zucker et al. [75], O’Neill et al. [53], and O’Neill et al. [52], summarized as follows.

The gas volume of a given bin is

where # is the number of particles, $\mathrm{d}v$ the bin volume, $\mathrm{d}a$ is the area of the bin, and $N_{\mathrm{H}}$ is the hydrogen column density. Assuming a constant ratio between the hydrogen column density and extinction in a given band (see e.g., Draine [13]), it is possible to write

with $\alpha_X$ a coefficient that can be derived from the curve in the work by Draine [13]. Using the Gaia G-band extinction, with a mean wavelength of $\lambda= 673\ \mathrm{nm}$ (Jordi et al. [33]), the $\alpha$ coefficient is equal to $4 \times10^{-22}$.

For the L20 map, given the relationship between extinction $A_{\mathrm{G}}$ and optical depth ($A_{\mathrm{G}} = 1.086 \times\tau_{\mathrm{G}}$), Eq. A.1 becomes

For the E24 map, we know that the dimensionless differential extinction $A'_{\mathrm{ZGR23}}$ can be converted to a given extinction band as,

where $m_X$ is a coefficient from the extinction curve by ZGR23. For the Gaia G-band extinction, $m_X = 2.036$. From this, Eq. A.1 can be written as

Recently, O’Neill et al. [52] directly compared dust extinction with H I emission and found that the conversion factor previously adopted by O’Neill et al. [53], $1650\ \mathrm{cm}^{-3}$, is underestimated by a factor of $\sim1.8$. They recommend a conversion factor of $3002\ \mathrm{cm}^{-3}$, corresponding to the median of the distribution of possible values. This reflects the fact that this conversion is not described by a unique best value, but instead has a calibration uncertainty. For a more detailed discussion, we refer to O’Neill et al. [52] and their Fig. 16. We therefore assume that the same correction applies to the L20 map and multiply the conversion factor used by Zucker et al. [75] by a factor of 1.82.

The use of the updated conversion factors affects the absolute gas density and therefore leads to changes in the crossing statistics. However, the main conclusions presented in the paper are not significantly affected. For comparison, Table A.1 reports the crossing statistics obtained using the pre-correction conversion factors, providing an equivalent version of Table 1. The quantity obtained using this procedure represents the total volume density of hydrogen nuclei, independent of phase ($n = n_{\mathrm{H I}} + 2n_{\mathrm{H_2}}$).

Validation tests

In this Section, we present the validation tests for the encounter statistics presented in the main part of the paper. We first tested how the map orientation and the different initial conditions for the Sun’s orbit influence the results.

For the map orientation test, we reorient the map randomly each time the Sun re-enters the map volume during the integration. This test is designed to assess how sensitive the encounter statistics are to the large scale organization of the ISM. The E24 map is particularly useful for this purpose, because it covers a sufficiently large volume of the Solar neighborhood to include linear structures (see e.g., Kormann et al. [39]). We find that the encounter statistics are not significantly affected by the map orientation, as shown in Fig. B.1. This suggests that, for the statistical quantities considered here, the local ISM organization into clouds, the gas distribution, and the Solar orbit are the main drivers of the results, rather than the absolute orientation of large scale structures within the map. This conclusion is valid for the local Galactic volume probed by the present maps, but may not hold on larger scales, where the orientation of coherent and larger structures such as spiral arms might become important.

For the Solar initial condition tests, we estimated the encounter statistics using two different setups. In the first case, we kept the Solar position as in the reference case described in Sect. 3.1, but assumed a Solar velocity relative to the LSR of $(U_\odot, V_\odot, W_\odot) = (10.1, 4.0, 6.7)\,\mathrm{km}\,\mathrm{s}^{-1}$ (Hogg et al. [28]). As shown in Fig. B.1 (under Sun initial condition 1), the encounter statistics remain mostly unchanged. The main exception is the crossing duration in the L20 map, shown in the left middle panel, for which the median value is larger than in the reference case. This is likely due to the smaller Solar velocity adopted from Hogg et al. [28] compared with the value from Schönrich et al. [60] used in the reference integration, which increases the time required to cross individual density structures.

In the second case, we changed the Solar position, adopting a Galactocentric radius of $R_\odot= 8.33\,\mathrm{kpc}$ (Gillessen et al. [21]) and a vertical position of $z_\odot= 27\,\mathrm{pc}$ (Chen et al. [11]), while keeping the Solar velocity relative to the LSR fixed to the Schönrich et al. [60] value. Also in this case, shown in Fig. B.1 (under Sun initial condition 2), the resulting crossing statistics are consistent with those obtained in the reference case.

These tests show that the encounter statistics are robust against the Solar initial condition variations considered here. This is expected because these changes do not substantially modify the Solar orbit: in all cases, the Sun samples a similar range of vertical heights, where the ISM vertical stratification does not vary strongly, and the Sun velocity changes only moderately. This is also in agreement with the study by Kokaia & Davies [36], where they show that the $xy$ kinematic components of the stars do not affect the clouds hit rate.

We also tested the case in which the E24 volume is restricted to the volume covered by the L20 map, in order to verify whether the encounter statistics derived from the two maps become more consistent. As a first approximation, we repeated the analysis using only the E24 bins whose geometric centers lie within the L20 map boundaries. In this test, the E24 bins are still affected by the variation of bin volume with distance from the map center, although this effect is reduced. The radial size of the selected E24 bins ranges from 0.4 to 3.3 pc. As shown in Fig. B.2, the encounter statistics derived from the two maps become more similar when the same volume is considered, although differences remain. The agreement improves in particular for densities below $\sim100\,\mathrm{cm}^{-3}$, where the density regimes are equally well resolved in L20 and E24. This test therefore provides a useful consistency check, showing that the differences between the two maps are not produced by errors in the analysis pipeline, but mainly reflect their spatial discretizations.

Vertical stratification of the E24 map

In this section, we analyze the vertical stratification of the E24 map, following the same approach used for the L20 map in Sect. 4.2. Since the detailed vertical distribution of dust is not the main focus of this work, we perform a limited analysis.

Figure C.1 shows the side projection of the E24 map, together with the vertical distribution of the density bins and the volume filling factor for different density thresholds. We show these quantities for different radial selections of the map: the full E24 volume, the region within 625 pc, the region beyond 625 pc, and the region beyond 1000 pc. These distance cuts are adopted because the E24 map covers a large volume of the Solar neighborhood and is therefore increasingly affected by distance dependent resolution effects.

We find that, in all four cases, the low density gas follows the expected decrease with increasing $|z|$. For higher density gas, the behavior is more complex. As found for the L20 map, the E24 map restricted to within 625 pc shows a distribution that cannot be described by a simple monotonic decline. Dense regions appear to be distributed relatively uniformly within $|z| \lesssim100 - 130$ pc, before declining rapidly at larger heights. This behavior is not recovered in the outer regions of the E24 map, where the dense gas distribution becomes more strongly centered around $z = 0$. In these outer regions, however, the highest density structures are also largely absent, as shown by the vertical distributions of the density bins. This is likely a consequence of the loss of spatial resolution with distance, which smooths compact dense structures and reduces their peak reconstructed densities.

A more detailed study of how the dust scale height depends on density threshold would therefore be valuable. However, such analysis is beyond the scope of the present work.

Height dependence of the incident flux and dust accretion

In Fig. D.1, we show how the incident gas flux and the inferred accreted dust mass depend on the Galactic vertical height sampled by the Sun. The histogram shows the normalized distribution of time spent by the Sun at each height. The two panels show the same height distribution, but are color coded by different quantities.

In the left panel, the colors indicate the mean instantaneous gas flux in each vertical bin, normalized by the orbit averaged gas flux. This quantity measures the typical flux experienced by the Sun when it is located at a given height, independently of how much time the Sun spends there. It therefore reflects the combined effect of the gas density sampled at that height and the relative velocity between the Sun and the map.

In the right panel, the colors indicate the fractional contribution of each vertical bin to the total accreted dust mass. For fixed dust to gas ratio, effective cross section, and filtering efficiency, this quantity is equivalent to the fraction of the total integrated gas flux accumulated in each height bin. It therefore depends not only on the instantaneous gas flux, but also on the residence time of the Sun at that height.

This comparison illustrates why the regions of largest instantaneous exposure do not necessarily dominate the total accumulated dust mass. In particular, the instantaneous flux is enhanced closer to the Galactic midplane, whereas the total accreted mass receives a larger contribution from heights where the Sun spends more time, close to the vertical turning points of its orbit.