ABSTRACT

We combine parallax distances to nearby O stars with parsec-scale resolution three-dimensional dust maps of the local region of the Milky Way (within 1.25 kpc of the Sun) to simulate the transfer of Lyman continuum photons through the interstellar medium (ISM). Assuming a fixed gas-to-dust ratio, we determine the density of ionized gas, electron temperature, and Hα emissivity throughout the local Milky Way. There is good morphological agreement between the predicted and observed Hα all-sky map of the Wisconsin Hα Mapper. We find that our simulation underproduces the observed Hα emission while overestimating the sizes of H ii regions, and we discuss ways in which agreement between simulations and observations may be improved. Of the total ionizing luminosity of 5.84×1050 photonss1, 15 per cent is absorbed by dust, 64 per cent ionizes ‘classical’ H ii regions, 11 per cent ionizes the diffuse warm ionized medium, and 10 per cent escapes the simulation volume. We find that 18 per cent of the high-altitude (|b|>30) Hα arises from dust scattered rather than direct emission. These initial results provide an impressive validation of the three-dimensional dust maps and O-star parallaxes, opening a new frontier for studying the ionized ISM’s structure and energetics in three dimensions.

1 INTRODUCTION

Understanding the creation and transport of Lyman continuum (LyC) photons from massive stars is a key element in determining the nature of galaxies as star formation engines. This ionizing radiation is mostly reprocessed into hydrogen and helium recombination and thermal dust emission, while some fraction of it can escape galaxies to produce the intergalactic radiation field. Measurement of the effects of this radiation allows us to quantify the star formation rate of galaxies and then compare to the structural and evolutionary factors that may regulate star formation (see, for example, Kennicutt & Evans 2012).

Observations of the local Milky Way provide an opportunity to study this process in detail because we can resolve individual ionizing sources and the structure of the surrounding gas and dust. In the 1950s, a primary focus of these investigations lay in determining the spiral structure of the Milky Way by using wide-area Hα sky surveys (Sharpless & Osterbrock 1952; Gum 1955; Sharpless 1959; Rodgers, Campbell & Whiteoak 1960) to identify nearby H ii regions and measuring the distances to the exciting stars (Morgan, Whitford & Code 1953).

The original idealization of these H ii regions as ‘Stromgren spheres’ (Strömgren 1939) in a uniform medium has been modified over the subsequent decades as our understanding of the complexity of the interstellar medium (ISM) increased. Massive stars are observed to form as part of clusters in dense molecular gas, but the combination of stellar winds, ionizing radiation, and supernovae leads to a complex density structure for the ISM, potentially filling the volume with networks of hot (T106 K) bubbles and superbubbles, interstellar chimneys, and ‘galactic fountain’ flows (Cox & Smith 1974; Shapiro & Field 1976; McKee & Ostriker 1977; Bregman 1980; Norman & Ikeuchi 1989).

This framework had to be modified when pulsars in globular clusters were used to demonstrate a significant mass of warm ionized gas extending a few kiloparsecs above the Galactic plane (Reynolds 1989).1 This layer of warm ionized medium (WIM), sometimes called the ‘Reynolds layer’, or the diffuse ionized gas (DIG), exists in the volume beyond the classical H ii regions, and analogous layers were detected in other edge-on spiral galaxies (Dettmar 1990; Rand, Kulkarni & Hester 1990). In a series of papers, Reynolds outlined several mysteries associated with this gas: what powers it, what ionizes it, what supports it at heights greater than a thermal scale height for 104 K gas, what fraction of ionizing photons emitted escape the galaxy, and so on (see Haffner et al. 2009). These mysteries, combined with the realization that WIM emission was a significant source of foreground emission for cosmological studies (Reynolds 1992), led to the construction of the Wisconsin Hα Mapper (WHAM), which obtained the first velocity-resolved all-sky map of optical line emission from the DIG of the Milky Way (Haffner et al. 2003).

Although numerous investigations have addressed the issues above, the unknown values of the filling factors and density structure of different phases have limited our ability to make conclusive statements about the energization and physics of the ISM. Previous modelling efforts of LyC transport (Miller & Cox 1993; Dove & Shull 1994; Bland-Hawthorn & Maloney 1999; Wood & Reynolds 1999; Dove, Shull & Ferrara 2000) explored the effects of different density structures, while Zurita et al. (2002) used H i maps of NGC 157 as an input to explore the role of H ii region LyC leakage in ionizing the WIM. However, recent development of maps of the three-dimensional (3D) distribution of dust in the local ISM (Schlafly et al. 2014; Green et al. 2015, 2019; Lallement et al. 2019, 2022; Leike, Glatzle & Enßlin 2020; Vergely, Lallement & Cox 2022; Edenhofer et al. 2024) informed by parallax (Gaia Collaboration 2016, 2023) and extinction (Zhang, Green & Rix 2023) measurements for millions of stars make it possible to study the transport of Lyman continuum (and other) photons, turning the local ISM into a laboratory for the physical processes regulating the structure and dynamics of the Galactic ISM.

A realistic model of Lyman continuum transport would provide a 3D grid of free-electron density and hence local contributions towards dispersion measure along any sightline, crucial in interpreting dispersions from extragalactic fast radio bursts. The free-electron structure is also required for interpretation of surveys of Faraday rotation measure in constraining the 3D structure of the local magnetic field. Our simulations also estimate interstellar contributions of photoionized species with commonly used absorption lines (Mg, Ca, Na, etc.). They will help constrain Lyman continuum escape fractions with relevance to the study of cosmic reionization, give realistic fractions of ionizing photons lost due to dust absorption, and also help to identify the extent to which ionizing flux is able to escape individual H ii regions. An estimate of direct versus scattered emission on every line of sight is also possible, as well as the ability to study the morphology and nature of individual regions of interest. It is also possible to determine contributions from individual ionizing sources, and constrain zones of influence of individual stars. Generating synthetic images also opens opportunities to fit for ISM and stellar parameters such as dust-to-gas ratio, ionizing luminosity, and distances to every O star in the local volume.

In Section 2, we describe our simulations. In Section 3, we compare the resulting simulated Hα sky maps to observations and derive some key, albeit preliminary, properties of the WIM. In Section 4, we discuss numerous avenues for future investigations.

2 SIMULATION SET-UP

2.1 Monte Carlo photoionization and Hα dust scattering

We use the Monte Carlo photoionization code cmacionize (Vandenbroucke & Wood 2018; Vandenbroucke & Camps 2020), which is based on the earlier version of Wood, Mathis & Ercolano (2004). cmacionize simulates the transport of photons in a 3D density grid and calculates the equilibrium ionization and temperature structure for a given set of ionizing sources. We have updated the radiative and dielectronic recombination rates and include collisional ionization rates from the CHIANTI data base (Del Zanna et al. 2021). For each of the 10 iterations to convergence in the photoionization simulation, we utilize 109 Monte Carlo photon packets.

Having computed the 3D temperature and ionization structure, we calculate the Hα emissivity using the tables of Osterbrock & Ferland (2006). The resulting Hα emissivity grid is the input for a Monte Carlo simulation of dust scattering and absorption using a modified version of the code described in Wood & Reynolds (1999), with dust properties at Hα from Hensley & Draine (2023). The opacity of the dust grains at Hα per hydrogen nucleon is set as 2.32×1022 cm2H1, the dust albedo was 0.72, and the Henyey–Greenstein factor for forward scattering was 0.44. The dust-to-gas ratio for both the photoionization simulation and dust scattering simulation was set as 0.006 05 (Hensley & Draine 2023). For the dust scattering simulation, we use 4×109 Monte Carlo photon packets.

2.2 Density structure

The 3D density structure was derived from the differential extinction maps within 1.25 kpc of the Sun (Edenhofer et al. 2024). This work releases 12 realizations of their 3D dust map, alongside the mean and standard deviation of the samples. Here, we utilize the mean map. This work thus does not take into account the uncertainties in the dust map.

The maps are converted to total hydrogen density using the same technique as O’Neill et al. (2024) (based on the method in Zucker et al. 2021). We first use the extinction curves of Zhang et al. (2023) to translate the Edenhofer et al. (2024) map (in total extinction per parsec – E/pc) into Gaia G-band extinction per parsec (AG/pc) via multiplication by a factor of 2.04. Then, using the extinction to column density factor of AG/NH=4×1022 from Draine (2009) gives NH/pc of 5.1×1021 cm2pc1×E, equivalently 1652 cm3×E. The resulting density was interpolated on to a regular 10243 Cartesian grid using the interpolation script shipped with the data from Edenhofer et al. (2024). This script interpolates the log-density structure and returns the exponential of the result.

The simulation box is a cube of side 2.5 kpc with the Sun at the centre.

2.3 Ionizing sources

Our photoionization simulation has 87 known O stars with good spectral classifications within 1.25 kpc of the Sun. We used the Galactic O-Star Catalog (Maíz Apellániz et al. 2013)2 and additional 10 O stars from a compilation of H ii regions and associated ionizing sources by K. Jardine.3 The Wolf–Rayet (WR) star contained within γ2 Velorum in the Gum Nebula was also included. Barring nearby sources that have been spectrally misclassified or deeply embedded sources, we anticipate that we have a reasonably complete sample of WR and O stars within 1.25 kpc.

Distances were obtained from Gaia parallaxes (Gaia Collaboration 2021) and earlier catalogues (Fabricius et al. 2002; van Leeuwen 2007; Zacharias et al. 2013). Of sources within 1.25 kpc, the mean statistical distance uncertainty was 9 per cent with only nine sources with statistical uncertainties exceeding 20 per cent. We discuss the effects of distance uncertainties on our results in Section 4.

The ionizing luminosities and effective temperatures of each source were determined by interpolating tables 4–6 of Martins, Schaerer & Hillier (2005).

The effective temperatures were used to select an appropriate spectrum for each ionizing source from the WMBasic library (Pauldrach, Hoffmann & Lennon 2001). The ionizing luminosity and spectral shape of the WR star were taken from Crowther (2007). The effects of any missing O stars, B stars, and other sources are discussed in Section 4. Uncertainties in ionizing luminosities would change both the size of ionized regions and the total Hα intensity generated. Variations in spectral type will have less effect on the simulated Hα sky, but might alter the structure of other ionized species, which are not analysed in this preliminary work.

3 RESULTS

3.1 The Hα sky: the inside and outside view

Fig. 1 shows the simulated all-sky map of Hα from ionized gas within 1.25 kpc, alongside the total Hα emission measured by WHAM (Haffner et al. 2003, 2010). These survey data are accessible via whampy (Krishnarao 2019). We have applied a velocity cut to the WHAM data at ±12 kms1 to eliminate emission from outside the simulated volume. This is a helpful first-order approximation, but a full knowledge of the thermal and kinetic components in the regions near the boundary of the region is required to completely remove spatially non-local emission. Also shown is the intensity of Hα arising from dust scattering. The Hα excess/deficit is also displayed using a quality factor, Q, defined as

(1)

for simulated intensity ISIM and WHAM observed intensity IOBS. This is such that Q=0 represents a perfect match to the observational data, positive values represent a multiplicative factor of overestimation of Hα, and negative values represent a multiplicative factor of underestimation. Several regions are annotated in the WHAM image.

Mollweide projections (centred on $l=180{}^{\circ }$) of the observed H$\alpha$ sky from WHAM (upper left), the simulated H$\alpha$ sky (upper right), the intensity of simulated H$\alpha$, which originates from dust scattering rather than direct H$\alpha$ emission (lower left), and the quality factor of simulated H$\alpha$ surface brightness as described in Section 3.1 (lower right). Individual regions are circled in the H$\alpha$ maps by their numbers in the catalogues by Gum (G), Sharpless (S), and Jardine (J) with the exception of the following named regions: the Orion-Eridanus superbubble (OE), the Gum nebula (Gum), the California nebula (CN), the North American Nebula (NAN), the $\zeta$ Ophiuchi H ii region ($\zeta$ Oph) and the $\lambda$ Orionis H ii region ($\lambda$ Ori). Blue rectangles on the upper panels highlight regions of high H$\alpha$ in the WHAM map originating from outside our 1.25 kpc radius simulated volume, and thus do not appear in our simulated sky.
Figure 1.

Mollweide projections (centred on l=180) of the observed Hα sky from WHAM (upper left), the simulated Hα sky (upper right), the intensity of simulated Hα, which originates from dust scattering rather than direct Hα emission (lower left), and the quality factor of simulated Hα surface brightness as described in Section 3.1 (lower right). Individual regions are circled in the Hα maps by their numbers in the catalogues by Gum (G), Sharpless (S), and Jardine (J) with the exception of the following named regions: the Orion-Eridanus superbubble (OE), the Gum nebula (Gum), the California nebula (CN), the North American Nebula (NAN), the ζ Ophiuchi H ii region (ζ Oph) and the λ Orionis H ii region (λ Ori). Blue rectangles on the upper panels highlight regions of high Hα in the WHAM map originating from outside our 1.25 kpc radius simulated volume, and thus do not appear in our simulated sky.

There is a good morphological agreement between the observations and simulations. Of the 87 O stars, 65 can be associated with 45 Hα-selected H ii regions from WHAM, 38 of which are reproduced in our simulations, including the ζ Oph (Sh 2–27) and λ Ori (Sh 2–264) H ii regions, the California (Sh 2–220) and North American Nebula (NGC 7000/Sh 2–117), the Orion-Eridanus superbubble, and the Gum nebula. A zoom-in image of the Orion-Eridanus superbubble region (Fig. 2) shows good agreement, reproducing several arcs and shells, although Barnard’s Loop is less prominent in the simulations than the observations. Of the seven nearby H ii regions that are not reproduced by the simulation, most are in the inner Galaxy, |l|<30; the dust maps and O-star distances may merit re-examination in this direction.

Zoom-in of the Orion-Eridanus superbubble in H$\alpha$ from the WHAM sky (left) and our simulated sky (right). The solid blue contours on both images show lines of 10 R brightness in the WHAM sky, and purple dashed contours are 100 R on the WHAM image. The H$\alpha$ emission in the simulated view of Orion-Eridanus is slightly more extended than observed.
Figure 2.

Zoom-in of the Orion-Eridanus superbubble in Hα from the WHAM sky (left) and our simulated sky (right). The solid blue contours on both images show lines of 10 R brightness in the WHAM sky, and purple dashed contours are 100 R on the WHAM image. The Hα emission in the simulated view of Orion-Eridanus is slightly more extended than observed.

Most of the features in the WHAM map that do not appear in the simulations are for understandable reasons. Several high-latitude H ii regions are ionized by B stars, while other regions lie beyond the 1.25 kpc distance range of the 3D dust map, notably the l=1030 Sagittarius ‘spur’ (Kuhn et al. 2021) and the l=280300 Carina Arm ‘tangency’ region. Because the emission lines from these regions are bright and broad, they contaminate low-velocity local emission. Gaussian fitting is required to remove this distant emission.

The Hα excess/deficit shows that over most of the sky the simulation underpredicts the observed Hα emission, typically by a factor of 2. There are three principal reasons why this might be the case. First, we only consider ionization by local O stars. Including emission from B stars, subdwarf OB stars, hot white dwarfs, and other sources of LyC radiation will increase the predicted Hα emissivity. Secondly, especially at lower latitudes, our simulations are missing Hα emission from beyond 1.25 kpc. Similarly, any photons that would enter the box from outside the simulation volume are not included, so the Hα in the box would be expected to be in deficit, in particular towards the edges and top/bottom of the simulation box.

Some regions of very low Hα surface brightness are present in the simulated sky. In a future study, these low emission measure sightlines will be compared with the free-electron column density from pulsar dispersion measures.

Given the good morphological agreement of the 3D simulation with the Hα sky view, one can create reliable views of the Hα emission from different vantage points. Fig. 3 provides the face-on Hα view of our location in the Milky Way as seen by an external observer. Some H ii regions are compact and in dense and dusty regions; others are larger and more ring-like, including two prominent examples of broken rings near the Sun: the Orion-Eridanus superbubble and the Gum Nebula that appear to be venting into Galactic supershell GSH 238+00+09. The face-on view also shows regions of diffuse ionization, i.e. the warm ionized medium. These images are reminiscent of views of H ii regions in nearby galaxies, e.g. Hodge & Kennicutt (1983), Greenawalt et al. (1998), or Groves et al. (2023).4

Simulated top down view of the galactic plane around the sun in H$\alpha$ (left) and simulated neutral hydrogen column density (right). The circled regions are the same as in Fig. 1, with the addition of the blue ellipse outlining supershell GSH 238+00+09 (labelled GSH238). The light blue circles on the right panel show our O-star locations, with the sizes denoting the ionizing luminosities in $Q_{46}$ ($\rm {s}^{-1}/10^{46}$). The yellow crosses denote the location of the Sun in both panels. The galactic centre is to the right. In this H$\alpha$ projection, it appears that Orion-Eridanus and the Gum nebula are venting into GSH 238+00+09. For a 3D interactive figure, visit https://wim.meikleobney.net.
Figure 3.

Simulated top down view of the galactic plane around the sun in Hα (left) and simulated neutral hydrogen column density (right). The circled regions are the same as in Fig. 1, with the addition of the blue ellipse outlining supershell GSH 238+00+09 (labelled GSH238). The light blue circles on the right panel show our O-star locations, with the sizes denoting the ionizing luminosities in Q46 (s1/1046). The yellow crosses denote the location of the Sun in both panels. The galactic centre is to the right. In this Hα projection, it appears that Orion-Eridanus and the Gum nebula are venting into GSH 238+00+09. For a 3D interactive figure, visit https://wim.meikleobney.net.

3.2 The Hα sky: individual regions

Although the simulated all-sky map shows good agreement with WHAM, an examination of individual H ii regions points the way to future improvements. We highlight five individual regions that have been extensively studied, and refer readers to their respective entries in the two-volume Star Formation Handbook (Reipurth 2008a, b) for more information.

ζ Ophiuchi H ii region (Sh 2–27): This is the nearest O-star ionized H ii region (at a distance of 135 pc), with slightly elliptical R=5 bright Hα emission centred at l=5.8,b=+23.7. ζ Oph is a ‘runaway’ O star associated with Sco OB2 (Preibisch & Mamajek 2008). From the 3D grids of emission measure, we find that it is ionization-bounded (all the ionizing photons are absorbed) and no photons escape to ionize the low-density DIG in agreement with prior modelling by Wood et al. (2005).

The simulated H ii region has a lower intensity Hα halo extending 3 (7 pc) beyond the observed radius, but with similar ellipticity. This larger extent may be due to (1) too high an ionizing luminosity from the interpolated tables, (2) under-resolving higher density substructure in the nebula (as noted in Zucker et al. 2021), or (3) a drop in the mean dust-to-gas ratio around ζ Oph, resulting in an underestimate of the gas density. Potentially in support of the final option, we note that the 3D dust maps of Edenhofer et al. (2024) show a cavity of dust density centred on the expected location of ζ Oph; however, this could also be a physical cavity due to the dynamical effect of ionizing photons from ζ Oph. Uncertainties in the distance to ζ Oph could also contribute.

Orion-Eridanus superbubble: This large structure at l=190, b=40 (Bally 2008) is morphologically very similar to observations. The shape of the Eridanus loop at the negative latitude end of the bubble closely resembles the WHAM map, suggesting that the 3D structure of the superbubble in the dust map is reliable. A previous 3D analysis focused on Barnard’s Loop (Foley et al. 2023), but without radiative transfer. Our simulations do not convincingly reproduce this structure, possibly as a result of uncertainties in the stellar distances.

We find that photons from O stars in the superbubble are in general trapped within the low-density void (nH103 cm3), and do not ionize through the dense walls (nH10 cm3). However, many holes and channels in the walls of the superbubble are identifiable and allow a significant amount of radiation to leak and produce nearby DIG. The face-on view in Fig. 3 shows that the bubble appears to have broken open on the left side as seen from the Sun (l=220), potentially ‘venting’ into the Galactic supershell GSH 238+00+09 as predicted by Heiles (1998).

λ Orionis Hii region (Sh 2–264): Projected within the Orion-Eridanus complex (l=165, b=12), this H ii region (Mathieu 2008) is nearly circular in the WHAM data. As with ζ Oph, it is present in the simulation, but slightly larger than observed. In 3D, this H ii region appears as a hemispherical outward protrusion in the back wall of the Orion-Eridanus superbubble. The highest density in the dust map along the line of sight to λ Ori is at a distance of 420 pc, within 1σ of the distance to λ Ori: 386±60 pc. Slightly increasing the stellar distance would locate this star in a denser environment and produce a smaller, brighter H ii region. One could imagine using the combination of the high spatial resolution of the dust maps combined with Hα images to refine the positions of the ionizing sources within this volume.

Gum Nebula: This nebula (Pettersson 2008), centred on (l=260,b=10), consists of a lower higher density half (the IRAS Vela Shell) and an upper lower density ionized cap with a larger expansion velocity. It also contains a significant number of ‘bright-rimmed clouds’, dense clouds with Hα rims ionized by ζ Pup and γ2 Velorum. This nebula is recognizable in our simulation, but the upper filament has a lower Hα surface brightness than in the WHAM map. The simulated Hα in the region mostly outlines the boundaries of the dense shell and clouds, recently investigated by Gao et al. , in preparation. Like the Ori-Eri superbubble, in the face-on view in Fig. 3 this nebula appears to be fractured on the side that abuts the supershell GSH 238+00+09.

North American Nebula (Sh 2–117): The O star in our simulation with the highest ionizing luminosity is the O3.5 iii 2MASS J20555125+4352246 (referred to as the Bajamar star in Maíz Apellániz et al. 2016), with an ionizing luminosity of 5.1×1049 photonss1, representing more than 8 per cent of the total in our simulation. This source is responsible for ionizing the North American Nebula (NGC 7000 or Sh 2–117), and produces the highest Hα brightness anywhere on the simulated sky. This is consistent with WHAM observations where it is brightest region within 1.25 kpc.

3.3 Constraints on LyC transport and the Galactic WIM

Our Monte Carlo codes allow us to distinguish four different photon termination scenarios: absorbed by dust, absorbed by gas of n>0.1 cm3, absorbed by gas of n<0.1 cm3, or escape from the simulation.

We find that 15 per cent of photons are absorbed by dust, 64 per cent ionize the dense gas, 11 per cent ionize the low-density gas, and 10 per cent escape the simulation box (in all directions). Not all of the escaping photons will reach the intergalactic medium, as we are missing the absorption from the remainder of the galactic disc. Our estimate that 11 per cent of the ionizing photon budget from O stars goes to the ionization of the DIG shows good agreement with the observationally motivated estimates of around 15 per cent (Reynolds 1984; Haffner et al. 2009).

To ensure that these results are well converged with resolution, a lower resolution run was carried out using a Cartesian grid of 5123 grid cells. The estimates for Lyman continuum escape fraction, dust losses, dense/diffuse gas contributions, and high-latitude dust scattering are found to agree within 2 per cent of the total photon budget going from 10243 to 5123. The morphological match is also retained at lower resolution, although some grid artefacts become apparent in structures, which are very local to the Sun.

By volume, 41 per cent of our simulation grid is ionized (neutral fraction <0.01), 51 per cent neutral (neutral fraction >0.99), and 8 per cent partially ionized. The volume filling factor of DIG (density less than 0.1 cm3) increases with height above the galactic mid-plane. The DIG filling factor is 0.3 in the mid-plane, rising to 0.75 at z=500 pc and 0.5 at z=+500 pc. Additional ionizing sources and time-dependent recombination are required to accurately model the vertical extent of the DIG (McCallum et al. 2024).

4 CONCLUSIONS AND FUTURE WORK

We have presented our initial simulations of the 3D Hα sky, made possible by the dust extinction maps of Edenhofer et al. (2024), the parallax distances to nearby O stars (Gaia Collaboration 2016), and the radiative transfer methods developed in Wood & Reynolds (1999), Wood et al. (2004), Vandenbroucke & Wood (2018), and Vandenbroucke & Camps (2020). The simulated sky bears a striking resemblance to the observed Hα sky (Finkbeiner 2003; Haffner et al. 2003), particularly the substructure seen in the Orion-Eridanus superbubble and the Gum Nebula, although the overall Hα emission is less than expected. Three quarters of the observed O-star ionized H ii regions within 1.25 kpc appear in our simulations, with varying levels of morphological agreement. Given the promising agreement with the data, we provided initial estimates for the ultimate fate of Lyman continuum photons emitted in the local Milky Way and the first estimates of the observed Hα produced by dust scattering within this new realistic density structure.

Our simulations provide an impressive validation of the dust extinction maps, and add two important new elements to our evolving 3D understanding of the local ISM: the spatial variation in the local hydrogen and helium ionization fractions and the extreme ultraviolet radiation field. This opens numerous avenues to test, refine, and extend these initial simulations, as discussed below.

Further tests of the simulations: Our preliminary assessment of the agreement between observed and simulated classical H ii regions can be greatly expanded. We find that modelled H ii regions are in general larger than observed. Obtaining better agreement with sizes, brightnesses, and morphology of observed regions will involve experimenting with the uncertainties in stellar positions, Lyman continuum production, the sharpness of density interfaces, density filling factors, the gas-to-dust ratio (and variations thereof), and parameters determining the dust scattering and absorption. Our experiments so far give us reason to believe that we can improve constraints on all of these factors by comparing simulations with data. With availability of the 12 underlying realizations of the dust map from Edenhofer et al. (2024), this also leaves the possibility of doing a Monte Carlo uncertainty estimation in the Hα sky by running a number of independent simulations, which sample the various dust map realizations and O-star positions.

The agreement between simulations and observations may be further tested by comparing integrated columns of ionized gas with dispersion measures of pulsars with parallax measurements (e.g. Ocker, Cordes & Chatterjee 2020). In addition, the detection of pulsar scintillation arcs (Stinebring et al. 2022; Ocker et al. 2024), which have been used as a probe of foreground ionized ‘screens’, may be compared with expected locations of ionized screens created by bubble boundaries in the simulations.

Although we have compared our simulation to the Hα sky maps, the same simulations may be used to predict the expected far-ultraviolet continuum from two-photon hydrogen recombination emission, radio free–free emission, radio recombination lines, and low-frequency thermal absorption of radio synchrotron emission from the Galactic disc and halo. This may provide insight, for example, regarding a long-standing puzzle about apparent inconsistency in inferred WIM properties based on Hα versus free–free emission (e.g. Dong & Draine 2011).

Further refinements of these simulations: Given the underproduction of Hα emission over much of the sky, it would be valuable to include the effect of ionizing sources in dust maps that extend to larger distances across the Galactic disc. Although these available maps have lower resolution than our adopted map, with Gaia Data Release 4 on the horizon, the quality and extent of the 3D dust maps are expected to increase. More distant ionizing sources will also contribute to the ionization of the high-altitude gas in the solar neighbourhood. Adding information about the vertical extent of the ISM above a kiloparsec would also allow for the measurement of the LyC flux as a function of height, improving our understanding of LyC escape from the Galaxy.

A second avenue of refinement would be the addition of other sources contributing to the ionization of the WIM, namely B stars, hot low-mass evolved stars (Flores-Fajardo et al. 2011), including subdwarf OB stars, hot white dwarfs, post-asymptotic giant branch (AGB) stars, and AGB manqué stars. Gaia-based catalogues for several of these classes of sources are now available, e.g. Pantaleoni González et al. (2021), Culpan et al. (2022), and Gentile Fusillo et al. (2021), although the Lyman continuum production of these different classes of sources is in general not well constrained.

While the simulations in this letter focused on the Hα sky, the equilibrium ionization state of metals is also tracked in the 3D volume. Synthetic emission line maps can be produced for collisionally excited optical emission lines (e.g. [N ii] λ6584 Å and [O iii] λ5007 Å). This capability will be of particular use for interpreting data from the Local Volume Mapper’s high-resolution optical spectroscopic survey of the Galaxy (Drory et al. 2024) and may also be compared to observations of the DIG of nearby face-on spiral galaxies.

Further extensions and ramifications of these simulations: As the work outlined above improves our understanding of the 3D distribution of the electron density and transport of ionizing radiation, it will become possible to use this improved understanding to examine other phases of the Galaxy. In particular, a reliable map of the 3D electron density may be used to better understand the magnetic field structure of the Galaxy as traced by Faraday rotation of pulsars and extragalactic sources (e.g. Dickey et al. 2022) and Faraday synthesis (Brentjens & de Bruyn 2005). In addition, one may explore the effect of other sources of diffuse ionization on spectral diagnostics, in particular the X-ray radiation from hot bubbles in the ISM and ionization by low-energy cosmic rays.

Perhaps the most remarkable result is encapsulated in Fig. 3, which shows the face-on view of the local part of the Milky Way disc. The 3D dust maps that have enabled this research have a resolution of roughly 1 pc. When viewed from a distance of 2 Mpc, this resolution would be equivalent to 0.1 arcsec. This means that the resolution of our face-on map of the local Galaxy is higher than could be achieved by space-based telescopes for any galaxy except for those within our Local Group. This extraordinarily high-resolution view of the Milky Way, which can be probed in 3D, will be a valuable resource for understanding the unresolved details in the star-forming discs of other galaxies.

ACKNOWLEDGEMENTS

LM acknowledges financial support from a UK-STFC PhD studentship. DK acknowledges support from HST-AR-17060 provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of 389 Universities for Research in Astronomy, Inc., under NASA 390 contract NAS5-26555. WHAM and LMH are supported by U.S. National Science Foundation award AST 2009276. We thank Michelangelo Pantaleoni for the helpful discussions relating to our sample of O stars.

KW, RAB, and LMH acknowledge the support and encouragement of Ron Reynolds in their early careers, which helped to prepare them for the discoveries presented here. Although Ron sadly passed away in April 2024, we fondly remember him as we continue to explore the mysteries of the warm ionized medium of the Galaxy.

DATA AVAILABILITY

3D cubes of electron density and Hα emissivity, a 2D sky map of Hα (at 15 arcmin resolution), and a table of O-star properties and positions are available for download at 10.5281/zenodo.15041318.

Footnotes

1

The existence of this gas had been posited to explain thermal absorption of low-frequency radio emission (Hoyle & Ellis 1963), but the height of this gas was observationally unconstrained.

4

For a particularly striking comparison, compare Fig. 3 to this amateur astronomer panorama of the Large Magellanic Cloud: https://www.cielaustral.com/galerie/photo95.htm.

REFERENCES

Bally
 
J.
,
2008
, in
Reipurth
 
B.
, ed.,
Handbook of Star Forming Regions: Volume I, The Northern Sky
, Vol.
4
.
Astron. Soc. Pac
,
San Francisco
, p.
459
 

Bland-Hawthorn
 
J.
,
Maloney
 
P. R.
,
1999
,
ApJ
,
510
,
L33
 

Bregman
 
J. N.
,
1980
,
ApJ
,
236
,
577
 

Brentjens
 
M. A.
,
de Bruyn
 
A. G.
,
2005
,
A&A
,
441
,
1217
 

Cox
 
D. P.
,
Smith
 
B. W.
,
1974
,
ApJ
,
189
,
L105
 

Crowther
 
P. A.
,
2007
,
ARA&A
,
45
,
177
 

Culpan
 
R.
,
Geier
 
S.
,
Reindl
 
N.
,
Pelisoli
 
I.
,
Gentile Fusillo
 
N.
,
Vorontseva
 
A.
,
2022
,
A&A
,
662
,
A40
 

Del Zanna
 
G.
,
Dere
 
K. P.
,
Young
 
P. R.
,
Landi
 
E.
,
2021
,
ApJ
,
909
,
38
 

Dettmar
 
R. J.
,
1990
,
A&A
,
232
,
L15

Dickey
 
J. M.
 et al. ,
2022
,
ApJ
,
940
,
75
 

Dong
 
R.
,
Draine
 
B. T.
,
2011
,
ApJ
,
727
,
35
 

Dove
 
J. B.
,
Shull
 
J. M.
,
1994
,
ApJ
,
430
,
222
 

Dove
 
J. B.
,
Shull
 
J. M.
,
Ferrara
 
A.
,
2000
,
ApJ
,
531
,
846
 

Draine
 
B. T.
,
2009
, in
Henning
 
T.
,
Grün
 
E.
,
Steinacker
 
J.
, eds,
ASP Conf. Ser. Vol. 414, Cosmic Dust – Near and Far
.
Astron. Soc. Pac
,
San Francisco
, p.
453
 

Drory
 
N.
 et al. ,
2024
,
AJ
,
168
,
198
 

Edenhofer
 
G.
,
Zucker
 
C.
,
Frank
 
P.
,
Saydjari
 
A. K.
,
Speagle
 
J. S.
,
Finkbeiner
 
D.
,
Enßlin
 
T. A.
,
2024
,
A&A
,
685
,
A82
 

Fabricius
 
C.
,
Høg
 
E.
,
Makarov
 
V. V.
,
Mason
 
B. D.
,
Wycoff
 
G. L.
,
Urban
 
S. E.
,
2002
,
A&A
,
384
,
180
 

Finkbeiner
 
D. P.
,
2003
,
ApJS
,
146
,
407
 

Flores-Fajardo
 
N.
,
Morisset
 
C.
,
Stasińska
 
G.
,
Binette
 
L.
,
2011
,
MNRAS
,
415
,
2182
 

Foley
 
M. M.
 et al. ,
2023
,
ApJ
,
947
,
66
 

Gaia Collaboration
 
2021
,
A&A
,
649
,
A1

Gaia Collaboration
,
2016
,
A&A
,
595
,
A1
 

Gaia Collaboration
,
2023
,
A&A
,
674
,
A1
 

Gentile Fusillo
 
N. P.
 et al. ,
2021
,
MNRAS
,
508
,
3877
 

Green
 
G. M.
 et al. ,
2015
,
ApJ
,
810
,
25
 

Green
 
G. M.
,
Schlafly
 
E.
,
Zucker
 
C.
,
Speagle
 
J. S.
,
Finkbeiner
 
D.
,
2019
,
ApJ
,
887
,
93
 

Greenawalt
 
B.
,
Walterbos
 
R. A. M.
,
Thilker
 
D.
,
Hoopes
 
C. G.
,
1998
,
ApJ
,
506
,
135
 

Groves
 
B.
 et al. ,
2023
,
MNRAS
,
520
,
4902
 

Gum
 
C. S.
,
1955
,
Mem. RAS
,
67
,
155

Haffner
 
L. M.
 et al. ,
2009
,
Rev. Mod. Phys.
,
81
,
969
 

Haffner
 
L. M.
,
Reynolds
 
R. J.
,
Madsen
 
G. J.
,
Hill
 
A. S.
,
Barger
 
K. A.
,
Jaehnig
 
K. P.
,
Mierkiewicz
 
E. J.
,
Percival
 
J. W.
,
2010
, in
American Astronomical Society Meeting Abstracts #215
. p.
415.28

Haffner
 
L. M.
,
Reynolds
 
R. J.
,
Tufte
 
S. L.
,
Madsen
 
G. J.
,
Jaehnig
 
K. P.
,
Percival
 
J. W.
,
2003
,
ApJS
,
149
,
405
 

Heiles
 
C.
,
1998
,
ApJ
,
498
,
689
 

Hensley
 
B. S.
,
Draine
 
B. T.
,
2023
,
ApJ
,
948
,
55
 

Hodge
 
P. W.
,
Kennicutt
 
R. C. J.
,
1983
,
AJ
,
88
,
296
 

Hoyle
 
F.
,
Ellis
 
G. R. A.
,
1963
,
Aust. J. Phys.
,
16
,
1
 

Kennicutt
 
R. C.
,
Evans
 
N. J.
,
2012
,
ARA&A
,
50
,
531
 

Krishnarao
 
D.
,
2019
,
J. Open Source Softw.
,
4
,
1940
 

Kuhn
 
M. A.
 et al. ,
2021
,
A&A
,
651
,
L10
 

Lallement
 
R.
,
Babusiaux
 
C.
,
Vergely
 
J. L.
,
Katz
 
D.
,
Arenou
 
F.
,
Valette
 
B.
,
Hottier
 
C.
,
Capitanio
 
L.
,
2019
,
A&A
,
625
,
A135
 

Lallement
 
R.
,
Vergely
 
J. L.
,
Babusiaux
 
C.
,
Cox
 
N. L. J.
,
2022
,
A&A
,
661
,
A147
 

Leike
 
R. H.
,
Glatzle
 
M.
,
Enßlin
 
T. A.
,
2020
,
A&A
,
639
,
A138
 

Maíz Apellániz
 
J.
 et al. ,
2013
,
preprint
()

Maíz Apellániz
 
J.
 et al. ,
2016
,
ApJS
,
224
,
4
 

Martins
 
F.
,
Schaerer
 
D.
,
Hillier
 
D. J.
,
2005
,
A&A
,
436
,
1049
 

Mathieu
 
R. D.
,
2008
, in
Reipurth
 
B.
, ed.,
Handbook of Star Forming Regions, Volume I: The Northern Sky
, Vol.
4
.
Astron. Soc. Pac
,
San Francisco
, p.
757

McCallum
 
L.
,
Wood
 
K.
,
Benjamin
 
R.
,
Peñaloza
 
C.
,
Krishnarao
 
D.
,
Smith
 
R.
,
Vandenbroucke
 
B.
,
2024
,
MNRAS
,
530
,
2548
 

McKee
 
C. F.
,
Ostriker
 
J. P.
,
1977
,
ApJ
,
218
,
148
 

Miller
 
W. W.
 III,
Cox
 
D. P.
,
1993
,
ApJ
,
417
,
579
 

Morgan
 
W. W.
,
Whitford
 
A. E.
,
Code
 
A. D.
,
1953
,
ApJ
,
118
,
318
 

Norman
 
C. A.
,
Ikeuchi
 
S.
,
1989
,
ApJ
,
345
,
372
 

O’Neill
 
T. J.
,
Zucker
 
C.
,
Goodman
 
A. A.
,
Edenhofer
 
G.
,
2024
,
ApJ
,
973
,
136
 

Ocker
 
S. K.
 et al. ,
2024
,
MNRAS
,
527
,
7568
 

Ocker
 
S. K.
,
Cordes
 
J. M.
,
Chatterjee
 
S.
,
2020
,
ApJ
,
897
,
124
 

Osterbrock
 
D. E.
,
Ferland
 
G. J.
,
2006
,
Astrophysics of Gaseous Nebulae and Active Galactic Nuclei
,
University Science Books
,
Sausalito, California

Pantaleoni González
 
M.
,
Maíz Apellániz
 
J.
,
Barbá
 
R. H.
,
Reed
 
B. C.
,
2021
,
MNRAS
,
504
,
2968
 

Pauldrach
 
A. W. A.
,
Hoffmann
 
T. L.
,
Lennon
 
M.
,
2001
,
A&A
,
375
,
161
 

Pettersson
 
B.
,
2008
, in
Reipurth
 
B.
, ed.,
Handbook of Star Forming Regions, Volume II: The Southern Sky
, Vol.
5
.
Astron. Soc. Pac
,
San Francisco
, p.
43

Preibisch
 
T.
,
Mamajek
 
E.
,
2008
, in
Reipurth
 
B.
, ed.,
Handbook of Star Forming Regions, Volume II: The Southern Sky
, Vol.
5
.
Astron. Soc. Pac
,
San Francisco
, p.
235
 

Rand
 
R. J.
,
Kulkarni
 
S. R.
,
Hester
 
J. J.
,
1990
,
ApJ
,
352
,
L1
 

Reipurth
 
B.
, ed.,
2008a
,
Handbook of Star Forming Regions, Volume II: The Southern Sky
, Vol.
4
.
Astron. Soc. Pac
,
San Francisco

Reipurth
 
B.
, ed.,
2008b
,
Handbook of Star Forming Regions, Volume II: The Southern Sky
, Vol.
5
.
Astron. Soc. Pac
,
San Francisco

Reynolds
 
R. J.
,
1984
,
ApJ
,
282
,
191
 

Reynolds
 
R. J.
,
1989
,
ApJ
,
339
,
L29
 

Reynolds
 
R. J.
,
1992
,
ApJ
,
392
,
L35
 

Rodgers
 
A. W.
,
Campbell
 
C. T.
,
Whiteoak
 
J. B.
,
1960
,
MNRAS
,
121
,
103
 

Schlafly
 
E. F.
 et al. ,
2014
,
ApJ
,
789
,
15
 

Shapiro
 
P. R.
,
Field
 
G. B.
,
1976
,
ApJ
,
205
,
762
 

Sharpless
 
S.
,
1959
,
ApJS
,
4
,
257
 

Sharpless
 
S.
,
Osterbrock
 
D.
,
1952
,
ApJ
,
115
,
89
 

Stinebring
 
D. R.
 et al. ,
2022
,
ApJ
,
941
,
34
 

Strömgren
 
B.
,
1939
,
ApJ
,
89
,
526
 

van Leeuwen
 
F.
,
2007
,
A&A
,
474
,
653
 

Vandenbroucke
 
B.
,
Camps
 
P.
,
2020
,
A&A
,
641
,
A66
 

Vandenbroucke
 
B.
,
Wood
 
K.
,
2018
,
Astron. Comput.
,
23
,
40
 

Vergely
 
J. L.
,
Lallement
 
R.
,
Cox
 
N. L. J.
,
2022
,
A&A
,
664
,
A174
 

Wood
 
K.
,
Haffner
 
L. M.
,
Reynolds
 
R. J.
,
Mathis
 
J. S.
,
Madsen
 
G.
,
2005
,
ApJ
,
633
,
295
 

Wood
 
K.
,
Mathis
 
J. S.
,
Ercolano
 
B.
,
2004
,
MNRAS
,
348
,
1337
 

Wood
 
K.
,
Reynolds
 
R. J.
,
1999
,
ApJ
,
525
,
799
 

Zacharias
 
N.
,
Finch
 
C. T.
,
Girard
 
T. M.
,
Henden
 
A.
,
Bartlett
 
J. L.
,
Monet
 
D. G.
,
Zacharias
 
M. I.
,
2013
,
AJ
,
145
,
44
 

Zhang
 
X.
,
Green
 
G. M.
,
Rix
 
H.-W.
,
2023
,
MNRAS
,
524
,
1855
 

Zucker
 
C.
 et al. ,
2021
,
ApJ
,
919
,
35
 

Zurita
 
A.
,
Beckman
 
J. E.
,
Rozas
 
M.
,
Ryder
 
S.
,
2002
,
A&A
,
386
,
801
 

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.