A three-dimensional map of the hot Local Bubble using diffuse interstellar bands

TL;DR

This study maps the three-dimensional structure of the Local Bubble using diffuse interstellar bands, revealing the presence and resilience of DIB carriers within this hot, X-ray emitting cavity, and suggesting the existence of small-scale shielding structures.

Contribution

It provides the first 3D map of the Local Bubble using DIBs, demonstrating the stability of certain DIB carriers in harsh conditions and indicating the presence of shielding structures.

Findings

DIB carriers are present within the Local Bubble.

The $780$ DIB carrier withstands X-ray photo-dissociation.

Low $797/780$ ratio inside the Bubble suggests larger or shielded carriers.

Abstract

The Solar System is located within a low-density cavity, known as the Local Bubble, which appears to be filled with an X-ray emitting gas at a temperature of 10$^6$ K. Such conditions are too harsh for typical interstellar atoms and molecules to survive. There exists an enigmatic tracer of interstellar gas, known as Diffuse Interstellar Bands (DIB), which often appears as absorption features in stellar spectra. The carriers of these bands remain largely unidentified. Here we report the three-dimensional structure of the Local Bubble using two different DIB tracers ($\lambda$5780 and $\lambda$5797), which reveals that DIB carriers are present within the Bubble. The map shows low ratios of $\lambda$5797/$\lambda$5780 inside the Bubble compared to the outside. This finding proves that the carrier of the $\lambda$5780 DIB can withstand X-ray photo-dissociation and sputtering by fast ions, where the carrier of the $\lambda$5797 DIB succumbs. This would mean that DIB carriers can be more stable than hitherto thought and that the carrier of the $\lambda$5780 DIB must be larger than that of the $\lambda$5797 DIB. Alternatively, small-scale denser (and cooler) structures that shield some of the DIB carriers must be prevalent within the Bubble, implying that such structures may be an intrinsic feature of supernova-driven bubbles.

Equations

$$N=1.13\times 10^{17}\times \frac{W}{f{\lambda }^2}$$

$$\rho (X)=const.exp\left\{\frac{-1}{2}\left(X-X_0^T\right)C_0^{-1}\left(X-X_0\right)\right\}$$

$$\left[X-X_0\right]=C_X{G}^T{\left\{G{C}_X{G}^T\right\}}^{-1}\left\{G\left[X-X_0\right]-f\left(X\right)\right\}$$

$$N=g\left(m\right)={\int }_0^r{\rho }_{0}exp\left(m\left(r\right)-\frac{|rsin\left(b\right)|}{h_0}\right)dr$$

$${\tilde{m}}_{k+1}={\tilde{m}}_k+C_{post}\left\{G_k^T{\left(C_{obs}^{-1}\left(N_{obs}-g\left({\tilde{m}}_k\right)\right)\right)}^T-C_{prior}^{-1}\left({\tilde{m}}_k-m_0\right)\right\}$$

$$C_{post}={\left(G_k^T{C}_{obs}^{-1}{G}_k+C_{prior}^{-1}\right)}^{-1}$$

$$C_{prior}={\sigma }_m^{2}exp\left(\frac{-r^2+r^{{}^{\prime }2}-2r{r}^{{}^{\prime }}cos\left(\theta \right)}{{\zeta }^2}\right)$$

$${\chi }^2=\sum _{i=1}^n\frac{{\left(N_i^{obs}-N_i^{model}\right)}^2}{{\sigma }_{N_i}^2}$$

$$N_v=\frac{{\sum }_{i=1}^m\frac{l_i}{D}{N}_i}{{\sum }_{i=1}^m\frac{l_i}{D}}$$

Full text

{affiliations}

School of Astronomy, Institute for Research in Fundamental Sciences, 19395–5531 Tehran, Iran

Department of Physics and Astronomy, The University of Western Ontario, N6A 3K7, Canada

Lennard-Jones Laboratories, Keele University, ST5 5BG, UK

The Open University, Associate Lecturer Services (STEM), 351, Altrincham Road, Sharston, Manchester, M22 4UN, UK

The origin of the Local Bubble (LB) is unknown, but measurements of 60Fe column densities can be explained by successive explosions of massive stars (supernovae) within the Scorpius–Centaurus stellar group[13]. The high temperature within the cavity is hostile to atoms, molecules and dust grains[2, 3, 4]. Observations aimed at detecting highly ionized gas, which could be present at million-degree temperatures, reveal no such plasma within the LB[14]. The DIB carriers offer an alternative tracer of the wall and interior of the LB. Although their nature is unknown, recent studies indicate that DIB carriers are likely carbon-based organic molecules[8]. They are universal and have been detected in different environments of Milky Way (MW) and within different galaxies[5, 6, 7] and in earlier partial results from our survey[9, 10, 11]. Our 3D map of the distribution of DIB carriers within the vicinity of the Sun opens a new window on a substantial fraction of interstellar carbon that may be locked up in the molecular carriers and that may play an important role in interstellar chemistry; this may include the C${}_{60}^{+}$ (“buckyballs”) anion which may be responsible for the DIBs at $9577$ Å and $9632$ Å wavelengths[15]. Therefore, we employed two of the strongest DIB tracers at $5780$ Å and $5797$ Å wavelengths, in direction of 359 different sightlines, to mapped the LB out to a distance of $200$ pc. While the $\lambda 5780$ DIB could possibly persist under harsh conditions of the LB[16, 17], the $\lambda 5797$ DIB cannot survive in such an environment and must be shielded within the inner regions of clouds[10].

The outlines of the LB have been mapped with absorption lines of singly ionized calcium (Ca ii) and neutral sodium (Na i) in the spectra of background stars[2], as well as by the attenuation of stellar light by interstellar dust based on $E(B-V)$ measurements[3, 18, 19], revealing that the LB cavity extends out to 80 pc in the Galactic plane (GP) and up to several hundred pc perpendiculars to the GP. Some dust maps have been produced by translating the strength of $\lambda$15273 DIB into $E(B-V)$ and merging it with other $E(B-V)$ measurements as priors; thus not strictly speaking a DIB map, but really a dust map [18]. Also, a pseudo-3D map of the $\lambda$8620 DIB within 3 kpc from the Sun have been produced recently[20]. They compared the DIB map with the dust distribution and found the $\lambda$8620 DIB and dust to have a similar distribution, but the scale height of the DIB exceeds that of the dust. This is a clear indication that DIBs and dust do not necessarily trace the same interstellar material and/or conditions[9]. Their pseudo-3D map traces large-scale structure, with the whole LB being confined to one voxel.

Because of the high temperature of the LB, it is difficult to detect ordinary atoms within the LB. However, our observations show that DIB carriers are present within the LB[9, 10, 11]. An example of different DIBs within and around the LB is shown in Fig. 1. To ensure that the observed DIBs are physical features not noise, only those absorptions with a confidence level exceeding $3\sigma$ are accepted in the mapping. In the following, we present three principal slices of the 3D distribution of the DIB pseudo volume density within $200$ pc from the Sun. The maps are in Galactic coordinates, i.e. the Sun is located at the origin and the primary direction is towards the Galactic Center (GC). The slices show the maps in imaginary planes of the Galactic, meridian and rotational planes. The GP slices the MW disk, the meridian plane is perpendicular to the GP with its x-axis pointing towards the GC, and the rotational plane is perpendicular to the GP and faced towards the GC.

In upper left panel of Fig. 2 ($\lambda 5780$ in the GP) the bulk of the material lies in front of the Scutum and Aquila dark nebulae (in the 1st quadrant in the upper-right sector) at distances of 30–120 pc from the Sun. The outer layer of this dense DIB structure coincides with a dense concentration in the dust map, however, the inner region is stretched within the LB and continued to interrupt the whole LB where there is not any dust[3] (see Supplementary Fig. 1 for an RGB map of dust, DIB and Na i). A lower-density DIB filament lies in front of the Cygnus rift molecular clouds. In this direction, X-rays with energies of order keV probably affect the DIB carriers beyond a distance of approximately $100$ pc[21]. In the 2nd quadrant (upper-left sector) DIB carriers are abundant in the direction of the Taurus dark clouds and molecular. The 3rd quadrant (lower-left sector) is characterized by a general paucity of the DIB carrier, in the direction of the $\beta$ Canis Majoris (CMa) interstellar tunnel and the GSH 238+00+09 supershell ($l=260$°). However, a low-density DIB trunk defies the odds; it may be associated with the photo-ionizing effect of $\beta$CMa [2]. Finally, the 4th quadrant (lower-right sector) features the most famous cavity connected to the LB (Loop I) in the direction of the GC ($l=345$°) at a distance of $200$ pc. A narrow tunnel connects the LB with Loop I, as revealed in previous LB maps[2, 3]; our map reveals that it is filled with a DIB filament and that Loop I itself is filled with DIB material as well (see Supplementary Fig. 2 for more slices).

The vertically extended structures are seen in the meridian plane (middle left panel of Fig. 2), one towards the GC in front of the Ophiuchus and Lupus complexes and R Coronae Australis (CrA) and another in the opposite direction in front of the Taurus star formation complex in between the Hyades and Pleiades star clusters. The atomic maps of the LB in this view show a striking open-ended tunnel known as the Local Chimney, tilted at an angle of 35$\degree$ from the North Galactic Pole (NGP)[22, 2], but our map shows that the tunnel is not entirely devoid of DIB carriers. Interestingly, there are some high latitude clouds (H i shells) in this direction [23]. The tilted tunnel can also be discerned in the rotational plane, but most noticeable is the DIB complex in the 3rd quadrant and material some $150$ pc below the GP in the direction of the South Galactic Pole.

In the right column of Fig. 2, we present the $\lambda 5797$ distributions. The carrier of the $\lambda 5797$ DIB is highly susceptible to the destructive effect of energetic photons[9, 10]. The GP view shows the largest concentrations lie between $30$–$120$ pc distance in the direction of Scutum and a filament extends over $200$ pc to meet Loop I, both akin to what was seen in the $\lambda 5780$. Likewise, the Local Chimney is visible in the meridian plane and a thin filament coincides with the densest $\lambda 5780$ DIB concentration in the 3rd quadrant of the rotational plane. Apart from the similarity between the two DIB carriers there are also notable differences, with the $\lambda 5797$ DIB distribution generally more fragmented and more tenuous. The $\lambda 5780$ DIB is known to trace relatively energetic environments as compared to the $\lambda 5797$ that traces more neutral, shielded regions (see Supplementary Fig. 3 where the strength of the $\lambda 5780$ saturates at high Na i densities). Indeed, the ratio of the two DIB carrier distributions shows that the interior of the LB is depleted in the $\lambda 5797$ DIB carrier, more so than the $\lambda 5780$ DIB carrier (Fig. 3) though more neutral, shielded cloudlets seem to persist within the LB.

The occurrence of DIBs within the LB, Local Chimney and the pathway tunnels, emphasizes that the $\lambda 5780$ and $\lambda 5797$ DIBs cannot all trace dust perfectly, as some studies showed that different DIBs trace different parts of clouds[9]. Conversely, it implies that at least some DIBs will yield different information compared to dust maps.

The question arises how DIBs could survive in such a high temperature? To address this question, we should emphasize that, in general, two different processes could destroy atoms and molecules in the hot regions of the ISM. Photo-dissociation by X-ray photons from stellar winds and supernovae, and sputtering by ions and electrons in the hot plasma and behind fast non-radiative shocks. Both processes can lead to the loss of atoms from large, carbonaceous molecules (like PAH and DIBs) and causing molecular destruction. On the other hand, to better understand the LB environment, we could compare the LB with the Orion nebula and the M82 galaxy. The center of Orion nebula (Orion-S), has a moderate temperature (T = $10^{4}$ K) measured from [N ii], and a high electron density (n${}_{e}=10^{4}$ cm*-3*) estimated from the [S ii] forbidden doublet-line ratio[24], while the starburst M82 galaxy shows a high temperature (T = 5.8 $\times$ 10${}^{{}^{6}}$ K) and low density (nH = 0.013 cm*−3*) in its bipolar outflow[25]. Therefore, the LB (with nH = 0.01 cm*−3* and T=106 K) is more similar to M82 but the small colder regions within the LB (such as three sub-regions of Ca ii) could be more akin the Orion Nebula. In the galactic outflows from M82, the electron collisions destroyed PAHs with 50–200 C-atom, although the bigger molecules can survive longer, eventually, all PAHs are destroyed within several thousand years[12]. In the denser and colder regions of the Orion Nebula, the PAH erosion is not caused by electron sputtering. In these regions the ion (He) collisions damage the PAHs, yielding a timescale of 107 years for the PAH destruction[12].

Based on this discussion, in a high-temperature and low-density gas, small PAHs would be destroyed. Therefore, the survival of small molecules in such conditions requires a protective environment and/or an efficient reformation mechanism. For the small shielded regions of the LB’s interior, the DIBs can survive in the surface layers of cold clouds within the hot LB for several million years (similar to Orion)[12]. However, for non-shielded regions of high temperature, a possible reason for observing DIBs could be described with a formation mechanism. As some very low reddening is still seen within the LB (with E(B-V) $<$ 0.1 mag), a small fraction of interstellar grains [26] has still remained within the LB. When these grains are eroded by electron collisions, any DIB carriers stuck on them are released into the ISM. Therefore, we could observe weak DIBs within the LB after $\sim$14 Myr (estimated age of LB)[13] showing that before the creation of the LB, this part of the ISM was characterized by dense, cool gas. Therefore during the lifetime of the LB, DIB carriers have been released into the ISM at a steady rate, making it possible to see very weak DIBs within the LB to this day.

The evolution of multiple supernova shells is expected to instantiate dynamical instabilities within the LB. These have been linked to the formation of local cool clouds[27] while tunnels and holes in the wall may have resulted from the Vishniac instability[28]. Cold, dense structures may also form out of hot gas through thermal instabilities[29]. It thus appears plausible from our maps that DIBs such as the $\lambda 5780$ may be sensitive, and potentially unique tracers of such structures at the boundary between cool and hot ISM phases.

Our maps show the LB in a truly new light – the DIBs likely highlight the cloud surfaces and shallow clouds that are in direct contact with energetic electrons and/or photons, while the densest interstellar structures are better mapped with atomic species or dust (see Supplementary Fig. 1). Thus, the same supernovae that only a few million years ago may have given rise to the global increase in the 60Fe isotope found in deep-sea crusts, and which occurred at the same time as when Earth’s temperature started to decrease[30], may have descended upon Earth, or may do so in the future as Earth traverses through the walls of the LB.

0.1 Observations

To map DIB absorption in and around the LB, we conducted a high signal-to-noise (S/N) survey of 637 nearby early-type stars in both hemispheres (see Supplementary Fig. 4). The Southern hemisphere survey was carried out with the 3.5m New Technology Telescope (NTT) at La Silla, Chile[9], and the Northern survey has been conducted with 2.5m Isaac Newton Telescope (INT) at La Palma, Spain[10]. The Southern observations were performed during 10 nights from March 2011 to August 2012. The Faint Object Spectrograph and Camera (EFOSC2) was used at a spectral resolving power of $R=\lambda/\Delta\lambda=5500$, and covering a wavelength range from 5672 Å to 6772 Å. The Northern observations were carried out over the 35 nights from October 2011 June 2013. In these observations we used the Intermediate Dispersion Spectrograph (IDS) to cover the $\lambda\lambda 5780$, $5797$ and $5850$ DIBs in the $5750$–$6040$ Å region at a spectral resolving power of $R=\lambda/\Delta\lambda=2000$.

More than 60 flat frames (captured by quartz lamp), about 15 arc lamp exposures (CuAr+CuNe at the INT and He+Ar at the NTT), and a large number of bias frames were taken for each observing night. To obtain a typical $S/N>1000$, we observed 9–25 science frames for each target. For data reduction, first, we produced master frames of bias and flat field by combining their individual spectra. The master bias then has been subtracted from master flat and science frames, and finally, the science frames were divided by the master flat. Later the weighted optimal extracted spectra were calibrated using arc frames. Possible bright cosmic ray impacts on detector have been removed using the sigma-clipping method (see Supplementary Fig. 5 for a comparison between measurements for a select number of targets that had been observed both at the INT and NTT).

0.2 Measurements

Quantifying the equivalent widths (EW) for DIB features is extremely challenging, especially when they are weak, shallow, and/or blended with other spectral lines. Since the DIB profile is probably composed of an overlapping set of an unknown number of transitions, the overall shape of their profiles is unexplained (and differs between DIBs). Single cloud line-of-sights confirm that most of them do not have a perfect Voigt profile. For instance, high-resolution spectra of the $\lambda$5797 profile reveal substructure. Because we are using low-resolution spectrographs, the shape of $\lambda$5780 and $\lambda$5797 are matched by a Gaussian profile[31]. First, a normalized, rectified spectrum is constructed by fitting a low order Legendre polynomial, by choosing $\pm$1–2 Å around the central DIB and masking the DIB feature itself (note that in figure 1 we only plot a small range of the continuum in order to highlight the fit to the DIB profile). Then, the line widths of $\lambda$5780 and $\lambda$5797 DIBs were obtained in terms of the standard deviation ($\sigma$) of Gaussian fits, and the FWHM calculated as FWHM = 2(2ln 2)${}^{1/2}\sigma$ = 2.355$\sigma$. The EW is then calculated from the integral of the Gaussian fit.

The major source of uncertainty in DIB’s EW is their blending with other stellar and interstellar lines. For instance, the $\lambda$5780 DIB could be blend with $\lambda$5778 DIB feature[32]. For each DIB, the statistical-uncertainty computed by summing the standard deviation of the Gaussian fit’s residuals[31, 17]. However, the main source of EW uncertainty in DIB studies is the systematic-error of continuum position. This error was estimated based on fitting three different continuum lines in $\pm$2 Å range around the peak (linear fit, quadratic fit, and simultaneous fit to the DIB and a linear continuum[33]). Later, the interval between the highest and lowest values of EW was taken as the systematic-uncertainty.

0.3 Data selection

We used high-quality data from the South and North DIB surveys with absorptions stronger than $3\sigma$ (where $\sigma$ is the standard deviation of the noise). Cool stars, with spectral types later than A3, were omitted to avoid contamination from stellar spectral lines. The latter was further mitigated by only using those sightlines with equivalent spectral width exceeding $6$ mÅ and DIB full width at half maximum (FWHM) exceeding $0.5$ Å. On the other hand, we consider sightlines with no discernible feature at the position of the $\lambda 5780$ and $\lambda 5797$ DIBs as zero column density sightlines. The absence of DIB features in these sightlines shows that the DIB density is negligible in some areas and moreover helps to reveal the fragmental distribution of DIBs in the LB. Also, since the spectrum of binary systems is contaminated by the stellar spectra, to avoid any non-real estimation, we reject those in our final catalog. Based on these selection criteria, among 637 observed targets only 359 sightlines used for mapping the LB. The measured DIB equivalent widths can be described by a density which decreases exponentially with distance from the GP, with a scale height of $100$ pc[2, 3]. This can be compared to the scale height of Na i, $h_{0}=170$ pc, and that of Ca ii, $h_{0}=450$ pc[2]. The equivalent width can be converted to column density by the following equation:

[TABLE]

Here, the wavelength is in Å, $W$ is the equivalent width in mÅ, and $f$ is the oscillator strength. The structure of the carrier must be known to be able to compute the oscillator strength but the DIB carriers are unknown. The $\lambda 5797$ DIB carriers tend to be present in the dense core of neutral Na i clouds, protected from the UV background radiation[10], while the $\lambda 5780$ DIB carriers are more abundant in the skin of the Na i clouds. On the other hand, both the $\lambda 5780$ and $\lambda 5797$ DIBs exhibit a strong correlation with Na i, understood as a general gas column dependency[9, 10]. Therefore, in order to be able to assign a value to the oscillator strength of the DIBs, we suppose that the DIB column densities are in the same range as the Na i column densities. By taking $f=0.3$ the Na i and $\lambda 5780$ DIB column densities lie in the same range (as shown in Supplementary Fig. 3). Thus we derived pseudo column densities for the carriers of the $\lambda 5780$ and $\lambda 5797$ DIBs.

All distances to the target stars are measured from parallaxes from the second Gaia data release (GDR2)[34, 35]. However, reliable distances for the majority of stars in the GDR2 cannot be obtained only by inverting the parallax. A correct inference procedure must instead be used to account for the nonlinearity of the transformation and the asymmetry of the resulting probability distribution considering the low signal-to-noise ratios of many measured parallaxes in the GDR2. Therefore we used the catalog of Bailer-Jones et al. (2018)[36] who provide purely geometric distance estimates for GDR2 sources.

0.4 Inverse method

In general, data values can be calculated directly from a given model (the ‘forward’ problem). On the other hand, to reconstruct a model from a set of measurements one must use an inverse method to estimate values at the positions where no direct measurements exist. It is possible to find the probability density of posterior estimates by assuming a Gaussian distribution[37] for both the data ($N$) and a priori model parameters ($m$). By defining the model parameters and data as a set of $X=[N,m]^{T}$ ($T$ is the transpose matrix) vectors, the a priori distribution and the model form a cloud in the parameter space. This cloud is centered on the observed data and the mean of a priori parameters. The shape of the cloud is described by the covariance matrix, $C_{X}$, which on its diagonal contains the covariance matrix of the data, $C_{obs}$, and the covariance of a priori values, $C_{prior}$[38]. The probability distribution takes the following form:

[TABLE]

One approach to estimating the optimal density distribution is to find the maximum likelihood point of $\rho(X)$ on the surface of $f(X)=0$. To estimate this likelihood point one should maximize $\rho(X)$ with the $f(X)=0$ constraint, while minimizing the exponential argument of $\rho(X)$, instead of determining the whole probability distribution[37]. This equation could be solved with Lagrange multiplier equations:

[TABLE]

Where $G$ is the matrix of partial derivatives $G=\partial g\left(m\right)/{\partial m}$. In the above equation, the variable $X$ appears on both sides and $f$ is a function of $X$; this makes it difficult to solve explicitly[37]. But the equation could be generalized in an iterative process, starting with some initial trial solution, $m_{0}$. If this a priori guess was close enough to the maximum likelihood point, then the successive approximations will converge to the true solution, $m_{est}$, within a few iterations.

In the case of constructing a 3D volume density distribution from a set of 2D column densities, $N=[N_{1},N_{2},...,N_{n}]^{T}$, with an inverse method along different sightlines, the most probable likelihood point could be found with an iterative Newtonian method[37, 38]. The volume density, $\rho$, is related to column density, $N$, by a simple integral equation, and the gas volume density in the ISM is known to approximately vary as a function of distance and Galactic latitude[2, 3] as $\rho_{0}(r)\exp(-|rsin(b)|/h_{0})$. Hence, to parameterize $\rho$ in each point of space, $N$ can be expressed as follows:

[TABLE]

Therefore, the observed column densities are related to volume densities by a model $N=g(m)$. To define the $g(m)$ function, the optimal $m$ parameter must be found with an inverse method. Using a Newtonian method one could estimate the optimal $m$ parameter with the following equations:

[TABLE]

[TABLE]

Here, $m_{0}$ is the a priori value of the parameter, $G^{T}$ is the transpose operator of $G$, the subscript $k$ refers to the iteration order, and $\widetilde{m}_{k}$ is the estimate of the m parameter in the $k$-th iteration. The value $C_{prior}$ represents the covariance between point $r$ and $r^{\prime}$, which in its simplest analytical form is thus[37]:

[TABLE]

Here, $\theta$ is the angle between the i-th and j-th sightline, $\zeta$ is the correlation length, and $\sigma$ is the a priori uncertainty for each point. At each iteration, a $\chi^{2}$ minimization criterion is used to control the algorithm convergence:

[TABLE]

The error on the observed column density ($\sigma_{N}$) is computed using similar methods as in previous works[2, 3].

The smoothing length implies that spatial detail below this characteristic scale will be smoothed out. The smoothing length used in this study is $\zeta=30$ pc which means that, if at a given point there is a deviation from the a priori model of given sign and magnitude, we want the deviation in the a priori model to be smoothed based on the column densities in the neighborhood of 30 pc boxes. The value for $\zeta$ was determined on the basis of the average distance between targets.

0.5 Priors

The inverse method is an ill-posed problem with a non-unique solution for each iteration as the estimated solution will depend on the initial guess. To determine the dependency of the distribution to the a priori guesses, we checked three different prior estimates. Firstly, we used a clumpy a priori guess for the initial value. We chose a clump as dense as the observed column density and change its position along the sightline at a fraction of $0.25$, $0.5$, and $0.75$ of the star’s distance. By using these prior values the iteration could not converge to an acceptable $\chi^{2}$ (the best-achieved $\chi^{2}=12.4$). In the second attempt, we divided each sightline into several equal segments and for each of those adopted the same initial value of the $m_{0}$ parameter based on the measured column density. With this initial guess, the inversion converged to $\chi^{2}=2.67$. However, some small dense clouds accumulated near the Sun despite a general accumulation of zero column density sightlines at these short distances. In the final attempt, we used an a priori Gaussian distribution centered at the middle of the sightline with the density distributed around $m_{0}$ (estimated based on the measured column density along that sightline) and a Gaussian width of $\sigma=2$, resulting in $\chi^{2}=0.76$. Considering our method is an iterative procedure it produces models that approach the data successively closer with each iteration. On the other hand, the ideal value for $\chi^{2}$ is unity, which would indicate that the measured densities and model are in agreement within the error variance. However, in cases where the errors on the data have been overestimated the $\chi^{2}$ would be less than unity. Actually, to prevent the algorithm from getting stuck in local minima instead of the global minimum, we deliberately increased the noise on the measured column densities. Therefore, our $\chi^{2}$ is a little below unity. If the $\chi^{2}$ value did not change in a subsequent iteration the algorithm was considered to have converged. In previous studies, the reduced $\chi^{2}$ for Na i was reported to be $1.5$, and $1.8$ for H i[2]. We varied the Gaussian width in the range of $\sigma\pm 2$ and found exactly the same results.

0.6 Approach

There are two approaches to solve equation 2. The first method is to divide the entire 3D space into equal voxels and estimate the model parameters in every single voxel. Another way is to estimate the model parameters along each sightline without dividing the space into separate voxels. In the first method, the un-smoothed distribution of stars within space, and the low number of observed sightlines in the case of small voxels cause some voxels not to be crossed by any star sightline. Therefore the number of constraints would be less than the number of voxels. In that case, the iteration will not converge to a correct answer. This problem can be circumvented by increasing the size of the voxels so as to put at least one star within every single voxel, to maximize the number of constraints. This is done at the expense of map resolution, which can result in the oversight of fine structures. The second method is the non-blocking method[37], which, instead of dividing the entire volume into voxels, only divides each observed sightline into segments and then estimates the volume density at every single point by inversion. In our case, we divided each sightline into $30$ segments while for the furthest point in a $540$ pc box we could estimate the volume density at each $30$ pc. Although most target stars are scattered at distances below $150$ pc, allowing a volume density estimate at each $5$ pc, by dividing the 3D space into $30^{3}$ pc voxels the final volume density in each voxel is the weighted mean value of all densities passing through the voxel. The weighted mean density in each voxel is computed based on the following equation:

[TABLE]

Here, $N_{i}$ and $l_{i}$ are the volume density and length of i-th segment, respectively, and $D$ is the space diagonal of the voxel. In cases where a zero column density sightline passes through a voxel crossed by some other non-zero column densities, the volume density within the voxel decreases. In the returned DIB map there are some strong absorptions seen through regions with low volume density; in simplistic map-making methods, these sightlines would cause the “fingers of God” radial striping.

Considering the correlation function represents the relation of each pair of points, the convergence of the algorithm sensitively depends on the acceptable fluctuation level of the model around the best estimate ($\sigma_{m}$). This value depends on the distribution and quantity of the data and is determined based on trial and error. To obtain the optimal value for $\sigma_{m}$ we tested different values and found that by using $\sigma_{m}=0.08$ the inversion converges to the minimum $\chi^{2}=0.76$ in only 5 iterations. If we choose a larger value for $\sigma_{m}$ the algorithm may diverge in the first iterations, while by picking very small values for $\sigma_{m}$ the process may need more than $100$ iterations before converging which in the case of non-linear inversions is prohibitive. After converging to an acceptable $\chi^{2}$ value, if the estimated column density was less than $0.5N_{obs}$ or greater than $2.5N_{obs}$ we considered it an outlier[2]. Additionally, we chose a very low a priori density, $\rho_{0}=10^{-10}$ cm*-3*, which is one dex below the a priori value for corresponding Na i atoms.

To check the reliability of the inversion, we simulated the result obtained by applying the method on an artificial distribution of clouds which is described in the Supplementary file.

0.7 Test the inverse method

To check the reliability of the inversion, we simulated the result obtained by applying the method on an artificial distribution of clouds. In this simulation we produced three rectangular clouds with different densities and separated from one another in order to check the inversion’s ability to determine the clouds’ edges. In the first test, we placed multiple stars in each voxel. In this case the number of constraints exceeds the number of voxels and the returned map is similar to original clouds in term of distribution and densities. We then estimated the cloud densities by using randomly distributed stars, fewer in number than voxels. The lack of constraints as a result of empty voxels renders the inversion method unable to accurately determine the shape of the clouds – especially at the cloud edges (See Supplementary Fig. 6). As another test, we randomly varied the observed column densities within their error range (based on a random Gaussian distribution) and applied the inversion method on the new dataset. The returned results were, in essence, identical to those from the previous test, and the map did not show significant differences. The outer layers of clouds showed very small changes which were, however, smoothed out once the final voxel median smoothing had been applied.

{addendum}

The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.