We present a new algorithm for identifying superbubbles in HI column density maps of both observed and simulated galaxies that has only a single adjustable parameter. The algorithm includes an automated galaxybackground separation step to focus the analysis on the galactic disk. To test the algorithm we compare the superbubbles it finds in a simulated galactic disk with the ones it finds in 21cm observations of a similar galactic disk. The sizes and radial distribution of those superbubbles are indeed qualitatively similar. However superbubbles in the simulated galactic disk have lower central HI column densities. The HI superbubbles in the simulated disk are spatially associated with pockets of hot gas. We conclude that the algorithm is a promising method for systematically identifying and characterizing superbubbles using only HI column density maps that will enable standardized tests of stellar feedback models used in galaxy simulations.
['Astrophysics - Astrophysics of Galaxies']
A New Superbubble Finding Algorithm Description and Testing
{'No Header': ', 1-10 (2024)', 'A New Superbubble Finding Algorithm: Description and Testing': 'Brock Wallin, 1 Benjamin D. Wibking, 1 ★ and G. Mark Voit 1  \n1 Department of Physics and Astronomy, Michigan State University, 567 Wilson Rd, East Lansing 48824, Michigan, USA  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'We present a new algorithm for identifying superbubbles in H I column density maps of both observed and simulated galaxies that has only a single adjustable parameter. The algorithm includes an automated galaxy-background separation step to focus the analysis on the galactic disk. To test the algorithm, we compare the superbubbles it finds in a simulated galactic disk with the ones it finds in 21 cm observations of a similar galactic disk. The sizes and radial distribution of those superbubbles are indeed qualitatively similar. However, superbubbles in the simulated galactic disk have lower central H I column densities. The H I superbubbles in the simulated disk are spatially associated with pockets of hot gas. We conclude that the algorithm is a promising method for systematically identifying and characterizing superbubbles using only H I column density maps that will enable standardized tests of stellar feedback models used in galaxy simulations.  \nKey words: galaxies: ISM - ISM: bubbles - techniques: image processing', '1 INTRODUCTION': "Spatially correlated supernovae can excavate large cavities called superbubbles within galactic disks, recognizable as hot, low-density, ionized regions surrounded by higher density shells of swept up gas (e.g., Mac Low & McCray 1988). Once gas in a galactic disk becomes self-gravitating, it gives birth to massive stars. Strong stellar winds from those massive stars then blow bubbles filled with hot GLYPH<16> 𝑇 ≳ 10 6 K GLYPH<17> gas (Weaver et al. 1977). Shortly thereafter, the massive stars end their lives with supernova explosions (e.g., Maoz et al. 2014) that locally disrupt star formation (e.g., Tacconi et al. 2020). When multiple stars in the same cluster go supernova within a few tens of millions of years of each other, their supernova remnants combine to increase the size of the bubble by adding energy to the bubble's interior. Together, clustered supernovae and strong stellar winds can thereby produce superbubbles hundreds of parsecs in size, capable of blowing through the top and bottom of a galactic disk, producing kpc-scale holes HI holes in it, and allowing the hot gas interior to the bubble to erupt into the circumgalactic medium (CGM) (e.g., Tomisaka & Ikeuchi 1986; Brinks & Bajaja 1986; Mac Low & McCray 1988).  \nSeveral recent analyses of supernova feedback and its impact on galaxy evolution have focused attention on the importance of superbubbles and how their properties depend on correlated supernovae (e.g., Orr et al. 2022; Oku et al. 2023; Vasiliev et al. 2023; Tan & Fielding 2024; Li et al. 2024). This paper is aimed at understanding whether numerical simulations of galaxy evolution produce superbubbles with properties similar to those in real galactic disks. As an example, we look at the superbubbles created by the magnetohydrodynamic (MHD) simulation of Wibking & Krumholz (2022) and compare them to the superbubbles found in a 21-cm radio observation of the face-on disk galaxy NGC 6946 (Walter et al. 2008).  \n- ★ Corresponding author: benwibking+astro@gmail.com  \nWe know of few such comparisons between simulated and observed superbubbles in the recent literature. This work is intended to bridge that gap.  \nThe simulated galaxy is Milky Way-like and its stellar feedback processes include both supernovae and ionizing photons from young stars. Star particles created in the simulation produce supernovae at a constant rate of 3 × 10 -4 SNe M ⊙ -1 Myr -1 at an age within 0 < 𝑡 age < 30 Myr. The sizes and locations of the superbubbles they produce depend on the spatial and temporal correlations of star formation in the simulated galaxy.  \nTo identify those superbubbles, we developed a simple and unbiased algorithm that can easily be applied to both simulations and observations of galactic disks. Several algorithms have previously been developed to automatically identify and measure interstellar filaments and bubbles in observed (e.g., Daigle et al. 2003; Men'shchikov et al. 2010, 2012; Collischon et al. 2021) and simulated (e.g., Williams et al. 1994) galaxies. Boomsma et al. (2008) performed an analysis of NGC 6946 and identified H I holes and bubbles using 21 cm column-density and velocity information. Similar analyses have been done for other galaxies, such as NGC 2403 (Thilker et al. 1998), IC 2574 (Walter & Brinks 1999), and IC 10 (Wilcots & Miller 1998).  \nThe algorithm presented here is simpler. It is designed to identify likely superbubble candidates using just an H I column density map with no velocity dispersion information. Simplicity was a key goal, so that human input could be kept to a minimum, thereby reducing confirmation bias in determining the superbubbles. The algorithm also determines the sizes and spatial distributions of the superbubbles, so that simulations and observations can be directly compared.  \nThe rest of this paper is organized as follows: Section 2 describes the simulated and observed galaxies used in the development and testing of the algorithm. Section 3 presents the superbubble finding algorithm in its entirety, from an input column density map to an output list of identified superbubbles and their properties. In section  \nFigure 1. Neutral hydrogen column density map of NGC 6946 using data acquired from the THINGS database (Walter et al. 2008).  \n<!-- image -->  \nFigure 2. Neutral hydrogen column density map of a face-on simulated galaxy from Wibking & Krumholz (2022).  \n<!-- image -->  \n4, we show the results of applying the algorithm to both the simulated galaxy and NGC 6946. Section 5 discusses the performance of the algorithm and suggests steps for future work to test and improve it. Section 6 summarizes our work.", '2 DATA SET': 'For development and testing of the algorithm, we used a 21-cm observation of NGC 6946 (Walter et al. 2008) and a simulated galaxy from Wibking & Krumholz (2022). Figures 1 and 2 show H I column density maps for these galaxies. Superbubbles from supernova explosions can be seen in both images as regions within the H I disk that have lower column density than the surrounding gas.', '2.1 NGC 6946': "We selected NGC 6946 for comparison with the simulated galaxy because its physical size is similar. The star formation rate of NGC 6946 is reported to be 4 . 76 M ⊙ yr -1 , and it has an H I mass of 4 . 15 × 10 9 M ⊙ (Walter et al. 2008). The superbubbles in NGC 6946 appear as fuzzy regions not much lower in column density than the rest of the galactic disk, making them challenging to identify by eye. Overcoming that challenge is one of the main goals of the algorithm. The simulated column-density map in Figure 2 is sharper because it does not include the noise and beam effects associated with the radio observation. We would like to determine superbubble properties that are independent of those effects.  \nOur algorithm takes as input H I column density maps that have pixel values in units of H cm -2 , so we first need to convert both the simulated galaxy and NGC 6946 data files into those units. Starting with NGC 6946, the FITS tables provide velocity-integrated H I 21 cm intensities ( 𝑆 Δ 𝑣 ) in units of Jy beam -1 ms -1 . To convert those intensities into velocity-integrated surface brightness temperatures, we account for the beam size using the formula  \n𝑇 B Δ 𝑣 = ( 6 . 07 × 10 5 ) 𝑆 Δ 𝑣 FWHM maj × FWHM min (1)  \nwhich gives 𝑇 B Δ 𝑣 in units of Kkms -1 when 𝑆 Δ 𝑣 is in Jy beam -1 kms -1 and the full width at half max of the beam's major axis (FWHM maj ) and minor axis (FWHM min ) are in units of arseconds (Walter et al. 2008). The H I column density is then  \n𝑁 HI = 1 . 823 × 10 18 𝑇 𝐵 Δ 𝑣 (2)  \nFinally, the outer regions of the NGC 6946 21 cm map can have unphysically small and sometimes negative flux values. Those are set to 10 18 H cm -2 .", '2.2 Simulated Galaxy': "Superbubbles in the simulated galaxy are more obvious than in NGC 6946becausetheminimumHIcolumndensitiesneartheircentersare much smaller. The simulated galaxy is evolved for 570 Myr with gas particles of mass 859 . 3 M ⊙ and total H I gas mass of 4 . 06 × 10 9 M ⊙ , and at the time of observation has a star formation rate of ∼ 2 M ⊙ yr -1 (Wibking & Krumholz 2022).  \nOur pre-processing procedure to create a synthetic H I map from the simulation follows these steps:  \nStep 1 - Density Conversion: The total gas density in the simulation is converted to hydrogen density using a hydrogen to total gas mass fraction of 𝑋 H = 0 . 71. Any unphysical gas temperature and density values are adjusted: temperature values in very low temperature regions are set to 1 K, and density values in very low density regions are set to 10 -6 g cm -3 . We also use temperature information to mask all density particles with temperatures greater than or equal to 10 4 K, in order to remove all ionized gas from the simulation.  \nStep 2 - Dust Optical Depth: The next few steps account for the ratio of neutral atomic hydrogen to molecular hydrogen, following Krumholz & Gnedin (2011) in its entirety. That ratio depends on the disk's optical depth. To determine it, we estimate the simulated disk's scale height using the formula ℎ = 𝜌 /|∇ 𝜌 | , where 𝜌 is the disk's gas density evaluated at each cell, and |∇ 𝜌 | is the magnitude of the density gradient, with the gradient calculated between all neighboring cells across the full grid. At times the gradient is very small which makes the scale height very large. To account for this, we set the maximum scale height of the simulated galaxy to be 100 pc,  \napproximately the scale height of the Milky Way's disk. This is done because large differences in scale height will have a large impact on the final column density values. The disk's column density Σ = 𝜌ℎ is then calculated and used to find the dust optical depth 𝜏 𝑐 = Σ 𝜎 𝑑 / 𝜇 H , taking 𝜎 𝑑 = GLYPH<16> 10 -21 cm -2 GLYPH<17> 𝜎 𝑑, -21 to be the dust cross section per hydrogen nucleus and 𝜇 H = 2 . 3 × 10 -24 g to be the mean mass per hydrogen nucleus.  \nStep 3 - Scaled Radiation Field: To quantify the scaled UV radiation field within the simulated disk, we start with this equation from Krumholz & Gnedin (2011):  \n𝜒 = 71 GLYPH<18> 𝜎 𝑑, -21 R -16 . 5 GLYPH<19> 𝐺 ' 0 𝑛 H , 0 (3)  \nHere, 𝜎 𝑑, -21 represents the dust cross section per hydrogen nucleus to 1000 photons, R -16 . 5 represents the rate coefficient for H 2 formation on dust grains, normalized to a Milky Way value of 10 -16 . 5 cm 3 s -1 , 𝐺 ' 0 represents the ambient UV radiation field intensity normalized to 2 . 72 × 10 -3 ergs cm -3 -1 (Parravano et al. 2003), and 𝑛 H , 0 represents the number density of hydrogen nuclei. We assume that GLYPH<0> 𝜎 𝑑, -21 /R -16 . 5 GLYPH<1> = 1 and also that the ambient UV radiation field intensity has a constant value of 1, so that  \n𝜒 = 71 𝑛 H , 0 (4)  \nStep 4 - H I Fraction: Using the dust optical depth and scaled radiation field from steps 3 and 4 we compute the 𝑠 parameter:  \n𝑠 = ln ( 1 + 0 . 6 𝜒 + 0 . 01 𝜒 2 ) 0 . 6 𝜏 𝑐 (5)  \nThis parameter is used by Krumholz & Gnedin (2011) to determine the fraction of neutral atomic hydrogen, via the fitting formula  \n𝑓 HI ≃ GLYPH<18> 3 4 GLYPH<19> 𝑠 1 + 0 . 25 𝑠 (6)  \nWe set 𝑓 HI equal to unity if this formula gives a greater value.  \nStep 5 - H I Column Density: Using the H I fraction calculated in the previous step, we find the projected column density map of the simulated galaxy using a 33 ° inclination to match the observation of NGC 6946 (Walter et al. 2008). We sum the H I density of the simulated galaxy along the line of sight at 33 ° to the plane of the galaxy to get a column density map in units of H cm -2 .  \nStep 6 - Beam Smoothing: With the simulation data in the form of a 21-cm H I column density map, we follow these steps to create a synthetic observation. We create an empty 2D array with the same size as the column density map and fill the array with a 2D Gaussian, centered at the central pixel. The 2D Gaussian array is created with a mean of 0 and a standard deviation of 𝜎 = 6 . 04 arcsec in both the x and y directions. This standard deviation is the FWHM maj of the beam used when observing NGC 6946. To combine the H I column density map and the 2D Gaussian array, we first take the 2D fast Fourier transform of each. Next, both transformed arrays are multiplied together element-wise, producing a 2D map of the combined column density and 2D Gaussian in Fourier space.  \nStep 7 - Adding Noise: To construct the synthetic 2D noise array, we again create an empty 2D array with the same size as the H I columndensity map. The synthetic noise is sampled from a 1D Gaussian with a mean of 0 and standard deviation found as follows. We take the 2D fast Fourier transform of the NGC 6946 H I map, then calculate the standard deviation of the transformed map and use this value as the standard deviation of the 1D Gaussian. Then we randomly sample from the 1D Gaussian for each element in the empty 2D array to  \nFigure 3. Pre-processed H I column density map of the simulated galaxy, inclined by 33 ° to match the inclination of NGC 6946.  \n<!-- image -->  \ncreate a synthetic 2D noise array in Fourier space.  \nBefore adding the 2D noise array to the H I map, we ensure the condition F [ 𝐼 ( 𝜔 )] = F [ 𝐼 (-𝜔 ) ∗ ] holds. This says that for any frequency of the 2D array found in Fourier space, the corresponding Fourier space element of the array with this frequency must be the complex conjugate of the Fourier space element of the array corresponding to the negative frequency. Each frequency is looped through to ensure this condition holds, and where this condition does not hold we set the negative frequency element to be the complex conjugate of the positive frequency element.  \nWedeal with the special conditions of the four corners and four edges separately. For each of the four corners, we take their absolute value and use this as their new value. We handle the edges minus the four corners in a similar way as to the rest of the image, but we only focus on either the 𝑥 frequency or 𝑦 frequency at a time, keeping the other the same. First, the element corresponding to the positive of the 𝑥 frequency is set to be the complex conjugate of the corresponding negative 𝑥 frequency keeping the 𝑦 frequency the same. We repeat this for the 𝑦 frequency, keeping the 𝑥 frequency component the same. When the Fourier transform is inverted to get a column density map in physical space, this step ensures that the noise-added map's column density values are real.  \nTo combine the 2D noise array with the synthetic beam applied H I map, both 2D arrays are in Fourier space already and we add them element-wise. Then we take the 2D inverse fast Fourier transform of the combined array to get the synthetic observation we will use in the superbubble identifying algorithm.", '2.3 Discussion of Images': 'The procedure used to make the simulated image entails many assumptions. To mitigate their effects, we applied Gaussian smoothing to both the NGC 6946 image and the pre-processed simulated image, so that small-scale structures are softened. The Gaussian kernel has a standard deviation corresponding to twice the semi-major axis FWHMvalueofthe applied beam. We use a factor of 2 here to ensure noticeable smoothing in the image without erasing the superbubblescale features we are trying to detect. Figure 3 shows the final result  \nFigure 4. Cumulative probability distribution functions and percentile thresholds for H I column density. NGC6946 is on the left and the simulation is on the right. Only the galactic disks are considered; the exterior regions are ignored. Colored lines relate the 50 th , 65 th , and 80 th area percentiles to the column density thresholds used in this analysis. All pixels with lower column density than a threshold value are assigned to candidate superbubbles by the algorithm.  \n<!-- image -->  \nfor the simulated galaxy and is the input image fed into the superbubble finding algorithm.  \nThere is a noticeable difference in minimum column density values between the simulated image and NGC 6946. Typically, NGC 6946 has larger H I column densities than the simulated image. These differences are most evident when looking at the centers of the superbubble regions, where the simulation has very low H I column densities. In NGC 6946 the central column densities within each superbubble region are larger, making the superbubbles less evident and more difficult to pick out.', '3 METHODS': 'The algorithm we have developed automatically identifies and measures H I holes, which may be physical H I superbubbles that have broken out of the galactic disk (see §4.3). The algorithm was designed with as few adjustable parameters as possible to minimize any confirmation bias. Unlike other algorithms mentioned in the introduction, we avoid using 21-cm velocity information in order to keep the model simple and usable with only a column density map. Our algorithm consists of the following steps:  \nStep 1 (Section 3.1): Given a column density map of neutral atomic hydrogen taken from a simulation or a 21-cm emission observation, first embed the image in a square array. Then, define an exterior region around the galactic disk and remove it from the analysis, leaving just the galactic disk.  \nStep 2 (Section 3.2): Using the disk map produced in Step 1, find all regions that fall below a column density threshold defined by a fixed percentile of column density by area.  \nStep 3 (Section 3.3): Using the regions identified in Step 2, mask out regions that are unlikely to be superbubbles. The remaining regions are considered candidate superbubbles.', '3.1 Exterior Removal': "Thefirst step the algorithm performs is to find and remove the exterior region around the galactic disk, so that we are left an image of just the galactic disk. By automatically removing the exterior, we can perform the algorithm on any input image.  \nStep 1a: Thealgorithmfirst selects an H I column density threshold for defining the disk. It takes the arithmetic average of the H I column density to be the threshold cut, and then calculates the image area above that threshold. Next, it computes an approximate radius 𝑟 corresponding to this area, using the formula 𝑟 ≈ √︁ 𝐴 / 2, dividing by 2 rather than 𝜋 to obtain a more conservative estimate of the disk's size. This choice decreases the possibility that some of the galactic disk could be incorrectly masked during the exterior removal process.  \nStep 1b: Next, the algorithm smooths the image with a uniform filter that has a smoothing scale equal to the radius 𝑟 from Step 1a. For each pixel, this filter averages all pixels within the smoothing scale and sets the pixel's value to that average. For pixels near the edge, pixel values within the filter that would extend over the edge of the image are taken to have the value of the nearest edge pixel.  \nStep 1c: The algorithm then uses one-fourth of the maximum columndensity in the smoothed image from Step 1b as a threshold to cut the smoothed image into two regions. Everything below the threshold is considered exterior. Everything above the threshold is considered the galactic disk. Using one-fourth of the maximum ensures that the cut will not be affected by extremely large exterior regions and will not require fine-tuning by the user. In each case we have tested, masking the image like this leaves only two contiguous regions, the galaxy region and the exterior, which can be cleanly separated.", '3.2 Identifying Candidate Superbubbles': "After the exterior is removed, the next step is identifying candidate superbubbles using isocontours of column density.  \nStep 2a: The column density threshold for a given galaxy is chosen based on a given column density percentile 𝑃 . The user chooses the  \nFigure 5. Superbubble candidates in NGC 6946 (left column) and the simulated galaxy (right column). The column density thresholds used are the 50 th (top row), 65 th (middle row), and 80 th (bottom row) area percentiles. Yellow regions are candidate superbubbles. Purple regions are parts of the galactic disk above the column-density threshold. Aqua regions are below the threshold but excluded from consideration as superbubbles because they fail one of the exclusion criteria. Crimson regions are considered exterior to the galactic disk. Each image lists the fraction of the disk's area covered with (yellow) superbubbles.  \n<!-- image -->  \nFigure 6. Bubble persistence maps for NGC 6946 (left) and the simulated galaxy (right). Each map is based on applying the algorithm seven times, with percentile thresholds going from 50 th to 80 th in steps of 5%. Colors show the number of thresholds at which a given pixel is part of a candidate superbubble.  \n<!-- image -->  \nvalue of 𝑃 for all galaxies of interest, and the algorithm calculates the corresponding column density thresholds for each galaxy. This method was chosen because of the systematic differences in typical column density values between NGC 6946 and the simulated galaxy. (Using a fixed percentile enables a direct comparison of the identified superbubbles between the two maps, even if the corresponding column density thresholds are very different.)  \nStep 2b: Using the column density threshold calculated in Step 2a, the algorithm selects all pixels in the galaxy image whose values are below the threshold, so that 𝑃 percent of the disk's area is below the threshold. Contiguous sub-threshold regions are considered candidate superbubbles. The algorithm numbers each candidate and measures its area and centroid using the regionprops function from the SciPy library.", '3.3 Excluding Unlikely Superbubbles': 'With the image processed and potential superbubble regions identified, the final step is to exclude regions that are unlikely to be real superbubbles. Once we remove unlikely superbubbles, we are left with regions that collectively cover a smaller area of the galactic disk than the selected percentile threshold. We call that fractional area the \'superbubble coverage fraction" and use that number to quantify the total area of the disk identified as superbubbles.  \nStep 3a: We exclude regions smaller than a minimum superbubble size. The minimum size is set by the size of the beam used in observing the galaxy, since features smaller than the beam size are unphysical. The algorithm uses a physical size corresponding to twice the size of the beam as the lower limit of a superbubble\'s approximate diameter. All regions that fulfill √︁ 𝐴 / 𝜋 < 𝐷 , where 𝐷 is the semi-major axis diameter of the beam, and 𝐴 is the total area of the superbubble, are masked out. We use twice the size of the beam rather than the beam size itself because we want to look for superbubbles that would be resolvable as extended features, rather  \nthan small scale fluctuations in the column density map that may be more likely to be caused by noise.  \nStep 3b: We exclude regions that are contiguous with the exterior. All pixels directly bordering the exterior, all pixels bordering these pixels, and so on, are identified. (A border with a width of 5 pixels is used to ensure all regions close to the exterior are removed.) Any candidate superbubble contiguous with the border is then excluded.', '4.1 Algorithm Output': "Here we discuss how the algorithm's results depend on choices for the area percentile, 𝑃 . We will focus on the 50 th , 65 th , and 80 th percentiles. Figure 4 highlights these area percentiles on the cumulative column density histogram for both NGC 6946 and the simulated galaxy. A particular area percentile applied to both NGC 6946 and the simulated galaxy corresponds to different column density thresholds. The column density threshold values used as the cut to determine the superbubbles for NGC 6946 and the simulated galaxy for the 50 th area percentile are 10 20 . 36 Hcm -2 and 10 20 . 13 Hcm -2 respectively. By using an area percentile rather than a specific column density enables us to directly compare these two galaxies and their superbubble properties.  \nAfter applying the algorithm to the simulated galaxy and NGC 6946, we produced images of the identified superbubbles for the selected area percentile. Figure 5 shows the results using the 50 th , 65 th , and 80 th area percentiles for both galaxies. The 50 th area percentile identifies the most superbubbles for both the simulated galaxy and NGC 6946. There is a trend in both the simulated galaxy and NGC 6946: the concentration of superbubbles moves towards the center of the galactic disk and the number of superbubbles found decreases as 𝑃 increases. Both galaxies have higher H I column density near their centers, so as the column density threshold increases, the identified superbubbles become more concentrated near the center of  \nthe galaxy. Increasing the area percentile also causes more of the galaxy to become below the threshold, allowing candidate superbubbles to be increasingly interconnected, contiguous with the exterior, and removed from consideration, as can be seen in Figure 5.  \nFigure 5 shows all the remaining superbubble regions in yellow and shows the excluded regions in aqua. The remaining purple region is the part of the galactic disk that is above the column density threshold. After both exclusion conditions are applied, the remaining superbubble candidates are labeled and saved. The superbubble coverage fraction reported in each panel represents the ratio of (nonexcluded) superbubble area (in yellow) to remaining galactic disk area (yellow, purple). It does not include the aqua parts consisting of low column density regions either contiguous with the exterior or smaller than the minimum size requirements.  \nIn the simulated galaxy, the superbubble coverage fraction decreases as 𝑃 increases, but this trend does not hold for NGC 6946. There the superbubble coverage fraction decreases as 𝑃 goes from 65 to 80 to 50. This difference arises because regions of high column density in NGC 6946 remain more contiguous as 𝑃 increases, allowing larger superbubbles to persist up to the 80 th area percentile.", '4.2 Bubble Persistence': 'We tested the stability of the algorithm by exploring what we call \'bubble persistence.\' We define bubble persistence as the number of area percentiles at which the same superbubble is found. Figure 6 provides a visual representation of bubble persistence. We take the superbubble images (yellow regions in Figure 5) and add them together. We add a total of 7 superbubble images, going from the 50 th to 80 th percentiles in steps of 5 percentile points, and count the number of times that a particular pixel belongs to a superbubble. Pixels that qualify as \'superbubble" many times have high \'bubble persistence" and indicate where a candidate superbubble is robust. There is a trend of higher bubble persistence towards the center of the galactic disk for both galaxies, which can be attributed to the generally larger column density near the center of the galaxy.  \nThe algorithm allows us to compare the radii of the identified superbubbles which we demonstrate by creating superbubble radius histograms for the 50 th , 65 th , and 80 th area percentiles using bins of 0 . 1 kpc. In Figure 7 we show superbubble radius histograms, weighted by the total area of each galaxy, to show how much of the galactic disk is covered with superbubbles at each radius. We see a shared trend: smaller superbubbles collectively cover more of a galaxy\'s disk than larger ones. The likely origin of this trend is that much more supernova energy is required to make a large bubble. Since we observe similar trends in galaxy coverage, we are encouraged that the algorithm can find superbubbles and measure their sizes consistently when given different column density maps.  \nFigure 8 presents similar histograms showing the number of superbubbles found within different ranges of galactic radius. The number of superbubbles found in each radial bin is weighted by the area of that annulus. Each bin covers 1 kpc in radius from the galactic center outward. Note that the fractional number of superbubbles decreases toward larger galactic radii.  \nFigure 9 shows the average column density profile surrounding each superbubble in both galaxies, again using the 50 th , 65 th , and 80 th area percentiles. As discussed earlier, the typical H I column density of NGC 6946 is larger than that of the simulation, especially near the bubble centers. The 80 th area percentile for the simulated galaxy is an outlier, since only 2 superbubbles are found and averaged.  \nTo calculate the average column density profiles, we selected 100 annuli with a 20 pc radial width about each superbubble\'s centroid,  \nFigure 7. Distribution of superbubble radius for NGC 6946 (blue) and the simulated galaxy (red) at the 50 th , 65 th , and 80 th area percentiles, from top to bottom.  \n<!-- image -->  \ngoing from the superbubble centroid out to 2 kpc. As expected, column densities near the centers of the superbubbles are lower than the rest of the profile, which increases to a maximum and slowly flattens out. There can be a slight peak before the average profile becomes flat. This may reflect a buildup of gas being bulldozed outward near the edges of the superbubbles.', '4.3 Superbubble Temperature Analysis': 'In order to check whether the candidate superbubbles identified in the simulation are real superbubbles, we perform an analysis using temperature information from the simulation. Starting with the superbubbles identified using the 50 th area percentile, we look at midplane gas temperatures at the locations of the superbubble centroids.  \nFigure 10 shows a midplane temperature map with the superbubble centroids plotted as dots or crosses. Here, the galaxy and its superbubbles are treated as if they are not inclined. The majority of superbubble centroids are located in regions of hot gas ( 𝑇 > 10 5 K). Those are denoted by red circles. However, some of the superbubble centroids, shown by blue crosses on the plot, are found in regions of cold ( 𝑇 < 10 5 K) midplane gas. The prevalence of red dots in the temperature contour map indicates that most of the candidate superbubbles are associated with regions of hot gas excavated by correlated supernovae.  \nWecan think of two reasons why some of the candidate superbubble centroids are located in regions of cold gas. The first is that when we identify the superbubbles we do so using the inclined galaxy image shown in Figure 3, so we first need to unincline the superbubble centroids to plot them on the midplane temperature slice map. Per-  \nFigure 8. Distribution of superbubble distance from the galactic center for NGC 6946 (blue) and the simulated galaxy (red) at the 50 th , 65 th , and 80 th area percentiles, from top to bottom.  \n<!-- image -->  \nforming this operation may cause the superbubble centroids to move off of a region of hot gas when plotted on the midplane temperature slice. The second is that the algorithm may combine several separate superbubbles into a single region because of their close proximity to one another. The centroids of the combined superbubbles will then be outside of regions of hot gas because they will be found between the individual superbubbles.  \nFigure 11 shows histograms of gas temperature across the entire midplane and at each superbubble centroid to compare their distributions. When comparing them, we see that the temperature distribution at the superbubble centroids is skewed to higher temperatures, confirming that more of the identified superbubbles are found in regions of high temperature. Most of the superbubble centroids are found in gas with temperatures above 10 6 K, similar to the theoretical values reported in Weaver et al. (1977).', '5 DISCUSSION': "Wehavedevelopedasuperbubblefindingalgorithm with only one adjustable parameter, the column density percentile threshold 𝑃 , in order to provide a robust and simple superbubble identification method and to minimize the possibility of confirmation bias by avoiding finetuning of a large number of parameters. The results of the algorithm are qualitatively insensitive to this value as long as it's between the 50 th and 80 th area percentiles, as demonstrated in Figure 6. The consistency of the results when using different area percentiles is a mark of success of our algorithm. Having the same candidate superbubble region being chosen even when using different area percentiles  \nFigure 9. Average superbubble column density profile for NGC6946 (blue) and the simulated galaxy (red) at the 50 th , 65 th , and 80 th area percentiles.  \n<!-- image -->  \ndemonstrates the robustness of the algorithm. However, the algorithm does have a few potential drawbacks which we want to discuss.  \nWe have already mentioned some improvements that we could apply to mitigate some of the known limitations, but we wish to present a few others. We will collect these into three categories: improvements that can be made to the algorithm itself, improvements to the synthetic H I column density map from the simulated data, and future work that compares simulated galaxies with observed galaxies when using the algorithm.", '5.1 Improvements to the Algorithm': 'One simplifying assumption we have made is to estimate the area of each superbubble as a perfect circle and find an approximate radius for each superbubble. This means that the algorithm is unable to determine if there are multiple superbubbles that are close together and combined into one region, because we do not measure the shape of the regions and only perform basic size and exterior connectivity masking. We could alter this and use an eccentricity measurement already calculated by the algorithm to find the semi-major and semiminor axis of each superbubble. Boomsma et al. (2008) shows an example of a method that does this in their identification of the superbubbles of NGC 6946. This improvement would preserve some of the shape information of superbubbles when analyzing trends of superbubble properties, e.g., allowing analysis of the distribution of superbubble radii along with their shape properties.  \nWe could also impose additional filters on candidate superbubbles using shape information. We expect superbubbles to be blown  \nFigure 10. Locations of candidate superbubble centroids in the simulated galaxy within a contour map of midplane temperature. A single temperature contour shows where 𝑇 = 10 5 K, and green shading shows where 𝑇 > 10 5 K. Red circles represent superbubble centroids located in hot gas. Blue crosses represent superbubble centroids located in cold gas. This set of centroids was found using the 50 th area percentile threshold, and the galaxy inclination is face-on.  \n<!-- image -->  \nFigure 11. Distribution of midplane gas temperature in the simulated galaxy for the superbubble centroids in Figure 10 (orange, right axis) and the entire galactic disk (green, left axis).  \n<!-- image -->  \nout along paths of least resistance in the ISM, meaning they are not necessarily circular as can be seen in Figure 10, but are expected to be small in scale when compared to the entire galactic disk. By including trends in real superbubble shape, the algorithm would better account for the non-standard shape of superbubbles when identifying candidate regions. With such a restriction, the algorithm may  \nbe able to more consistently find superbubbles even with changes in inclination and noise. An example of using improved requirements is shown through the Minkowski maps used in Collischon et al. (2021) when identifying bubbles. Using shape information could help the algorithm choose superbubbles more robustly, but could also introduce more confirmation bias by allowing for additional tunable parameters.', '5.2 Improvements to the Synthetic H I Map': 'The conversion from simulated galaxy data to a column density map could be made more realistic by using a spatially variable UV radiation field that is largest at the center of the galaxy and decreases radially outward. We expect this improvement may change the population of identified superbubbles near the centers of galaxies, since the very large UV radiation field in the centers of galaxies can make it more difficult for atomic hydrogen to become molecular (e.g., Krumholz & Gnedin 2011).  \nOur processing pipeline used to produce a synthetic H I map is very simple compared to the data processing pipelines used for radio interferometry observations (e.g., CASA Team et al. 2022). The most significant difference is that our synthetic H I map does not include the effect of sparse sampling of Fourier modes in the image plane due to a small number of antenna pairs. How this would affect our results is uncertain, but it is conceivable that this may explain the differences in radial profiles of column density apparent in Figure 9. Processing the simulation data with the same pipeline used for observations will be necessary in order to draw quantitative conclusions about how well the properties of simulated and observed superbubbles agree with each other.', '5.3 Future Work': 'Applying our algorithm to a wide range of galaxies, both observed and simulated, will allow us to quantify the variation in superbubble properties across galaxy populations, and will inform our expectations for how closely simulated and observed superbubble properties should agree. We used an H I map of NGC 6946 for most of our discussion above, but it is possible to apply this algorithm to all of galaxies in the THINGS H I survey. It will also be of interest to apply our algorithm to galaxies in other H I surveys in order to empirically test the robustness of our algorithm to differences in noise and beam effects from different radio telescopes and the processing algorithms applied by different radio surveys.  \nNew James Webb Space Telescope (JWST) infrared (IR) dust emission maps can also be explored with the algorithm, as these maps approximately trace H I and are of high resolution (Williams et al. 2024). Our algorithm is especially of interest for analyzing superbubbles in IR dust maps, since we only require 2D maps as input, rather than 3D position-position-velocity cubes used by some other H I superbubble algorithms. The high spatial resolution and low noise of these observations also makes them excellent candidates for comparison to galaxy simulations.  \nComparing the catalogues of superbubbles produced by our algorithm from input simulated and observed galaxies tests the realism of the implementation of supernova feedback in simulations in a direct way with only a single adjustable parameter. Because the algorithm was designed to allow for comparison of different galaxies, more analysis of a variety of observed and simulated galaxies could help test a broad range of simulation properties that could be affected by supernova feedback processes. Future work that puts error bars on superbubble properties in both simulations and observations will be necessary in determining the significance in differences between simulated and observed superbubble properties. This work will open the door to robust tests of supernova feedback models in simulations via their direct effects on the observed galactic gas.', '6 CONCLUSION': 'We have presented a new algorithm for identifying and measuring superbubbles in H I column density maps that requires little human input. We compared the candidate superbubbles identified in both a simulated galaxy and an observed galaxy (NGC 6946) to see whether the simulated galaxy created realistic superbubbles. Lessons learned during development of the algorithm motivated use of an area percentile method to find and compare superbubbles in H I maps of galaxies with different column density distributions.  \nWe found similarities in both size and radial distribution among the candidate superbubbles in the simulated galaxy and NGC 6946, indicating that the simulation may be producing realistic superbubbles. However, even with the similarities in superbubbles, there was no similar trend in superbubble coverage fraction between the simulated galaxy and NGC 6946. We use information about bubble persistence and the average column density profile to make additional comparisons. We are optimistic that applying the algorithm to a large sample of galaxy images, both simulated and observed, will improve its usefulness as a diagnostic tool of feedback models in galaxy simulations and provide a basis for a standardized comparison of superbubbles across observational and theoretical works.  \nWe hope that our work will lead to robust and reproducible tests of stellar feedback models used in simulations of galaxy formation and evolution.', 'ACKNOWLEDGEMENTS': 'This research was supported in part by the NSF through grant AAG2106575.', 'DATA AVAILABILITY': 'The data that support the findings of this paper are openly available in GitHub at https://github.com/BrockWW/ BubbleFindingAlgorithm.git .', 'REFERENCES': 'This paper has been typeset from a T E X/L A T E X file prepared by the author.'}