# BUBBLY: a method for detecting and characterizing interstellar bubbles using Fabry-Perot spectroscopy

A. Camps-Farina, J. Zaragoza-Cardiel, J. E. Beckman, J. Font, B.

Garcia-Lorenzo, S. Erroz-Ferrer, P. Amram

► To cite this version:

A. Camps-Farina, J. Zaragoza-Cardiel, J. E. Beckman, J. Font, B. Garcia-Lorenzo, et al., BUBBLY: a method for detecting and characterizing interstellar bubbles using Fabry-Perot spectroscopy. Monthly Notices of the Royal Astronomical Society, 2015, 447 (4), pp.3840-3848. 10.1093/mnras/stu2623. hal-01434366

HAL Id: hal-01434366

https://hal.science/hal-01434366

Submitted on 15 Sep 2021

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# BUBBLY*: a method for detecting and characterizing interstellar bubbles using Fabry-Perot spectroscopy

A. Camps-Fariña,<fnref n="1" />,<fnref n="2" />† J. Zaragoza-Cardiel,<fnref n="1" />,<fnref n="2" /> J. E. Beckman,<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> J. Font,<fnref n="1" />,<fnref n="2" />,<fnref n="4" />

B. García-Lorenzo,<fnref n="1" />,<fnref n="2" /> S. Erroz-Ferrer<fnref n="1" />,<fnref n="2" /> and P. Amram<fnref n="5" />

Accepted 2014 December 9. Received 2014 December 9; in original form 2014 November 5

## ABSTRACT

We present a new method for the detection and characterization of expansion in galaxy discs based on Hα Fabry-Perot spectroscopy, taking advantage of the high spatial and velocity resolution of our instrument (GHαFoS). The method analyses multipeaked emission line profiles to find expansion along the line of sight on a point-by-point basis. At this stage we have centred our attention on the large-scale structures of expanding gas associated with HII regions which show a characteristic pattern of expansion velocities, of the order of 100 km s⁻<fnref n="1" />, as a result of both bubble shape and projection effects. We show an example of the expansion map obtained with our method from a superbubble detected in the Antennae galaxies. We use the information obtained from the method to measure the relevant physical parameters of the superbubbles, including their ages which can be used to date young star clusters.

**Key words:** ISM: bubbles – HII regions – ISM: kinematics and dynamics – galaxies: starburst.

## 1 INTRODUCTION

The study of multicomponent emission lines is a topic of major interest as they are produced by many systems in astrophysics. It has long been known that the spectrum of an AGN has two components to its emission, characterized by a narrow and a broad peak each associated with two parts of the system each having very different kinematics (Antonucci 1993). The presence of an underlying broad supersonic component, which we now know corresponds to unresolved kinematic components, is also characteristic of bright HII regions (Arsenault & Roy 1986; Chu & Kennicutt 1994; Yang et al. 1996). In this case the observed emission peak may in fact be a combination of the main peak and two unresolved peaks associated with an expanding shell. Multiple peaks also appear in CO observations when independent clouds overlap in the line of sight (Blitz & Stark 1986; Stutzki & Guesten 1990), and outflows of gas have been observed in ionized and molecular gas (Bally & Lada 1983; Lada 1985; Bland & Tully 1988; Heckman, Armus & Miley 1990). It is apparent that the study of these multiple components of an

emission line, related to different kinematic phenomena, is relevant to a wide range of astrophysical problems.

Massive stars in OB associations and star-forming regions inject large quantities of mechanical energy into the surrounding interstellar medium (ISM) in the form of intense winds and supernova explosions that sweep up the material around them. This creates an expanding shell of dense material leaving a cavity of hot low-density gas inside. The individual bubbles created by each star combine quickly into one big superbubble associated with the region. A definition by Chu (2008) by which a bubble is blown by the stellar winds while a superbubble also has input from supernova explosions, can be used to make a formal separation between the two. Note that Meaburn (1980) called large shells around bursts of star formation ‘Giant Shells’. These are the ‘Superbubbles’ of Chu (2008) and not to be confused with the ‘Super-Giant Shells’ of Meaburn (1980). As the purpose of this article is to explain our detection technique we will not go more deeply into the subject here.

In this paper, we have used observations of the Antennae galaxies, where the linear resolution is such that objects smaller than a few tens of parsecs are not readily detectable, so that all our detections in this article are of superbubbles or giant shells, although the method is equally applicable to bubbles or expanding shells, given data on galaxies sufficiently close. Both phenomena have been studied theoretically and observationally for a long time (Dyson & de Vries

<footnote n="1">¹Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="2">²Department of Astrophysics, University of La Laguna, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="3">³CSIC, E-2806 Madrid, Spain</footnote>

<footnote n="4">⁴National Observatory of Athens, Lofos Nymfon, Thiesio, Athens, 11810, Greece</footnote>

<footnote n="5">⁵Laboratoire d'Astrophysique de Marseille, Aix Marseille Université, CNRS, F-13388 Marseille, France</footnote>

<footnote>*The code used in this paper is hosted at GitHub, it can be found here: https://github.com/artemic/Bubbly.git
†E-mail: artemic@iac.es</footnote>
1972; Castor, McCray & Weaver 1975; Weaver et al. 1977; Heiles 1979; McCray & Snow 1979; McCray & Kafatos 1987), and are a key aspect of feedback in galaxies. There have been detailed studies analysing their kinematics, mainly in relatively nearby giant H II regions such as 30 Doradus (Meaburn 1984; Chu & Kennicutt 1994) or NGC 604 (Yang et al. 1996), and other regions in the Magellanic clouds and several nearby galaxies (Meaburn 1980; Chu & Kennicutt 1986; Martin 1996; Nazé et al. 2001) as well as inside the Milky Way (Meaburn 1970; Chu et al. 1982).

As a superbubble expands and sweeps up the gas it disrupts the dense clouds where the stars are forming, so that star formation in the cluster may terminate with the formation of the superbubble. The effects on the surrounding ISM, however, can trigger further star formation in neighbouring regions (Gerola & Seiden 1978; McCray & Kafatos 1987; Palous, Tenorio-Tagle & Franco 1994), which implies that the superbubbles have a dual effect of quenching further star formation in the cluster that formed them and, at the same time, spreading star formation through the disc as each new star-forming region triggered by the superbubble will create more superbubbles in turn. This means mechanical feedback from H II regions plays a major role in galaxy evolution as it heavily regulates star formation in the disc; simulations of these effects provide support for this scenario as shown in Dalla Vecchia & Schaye (2008) and Ceverino & Klypin (2009).

One important consequence of the cut-off in time of star formation in the cluster as the ISM is swept into a superbubble is that the age of the superbubble coincides closely with the age of the star cluster. This means that if we can measure the dynamical parameters of the superbubble and compare them with an accurate model for its evolution the age of the cluster can be found in a straightforward way, as an alternative to (or complementing) spatially resolved observations of the cluster and fitting models for the stellar radiation such as STARBURST99 (Leitherer et al. 1999; Sharma et al. 2011), or using the Hertzsprung–Russell diagram (which requires resolving individual stars). It also allows studies of cluster ages to be carried out for objects at greater distances where the spatial resolution is limited.

Here we present a method for finding and mapping expanding components of the ionized gas throughout a galaxy or, for very nearby galaxies, within a sub-field. We use a Fabry-Perot spectrometer to map the Hα line profile throughout the galaxy and then search for the double peaked signature of expansion symmetrically spaced in velocity with respect to the main peak. It has been found in previous studies (Relaño & Beckman 2005; Relaño et al. 2007) that Hα profiles of bright H II regions have wings on both sides of the peak emitted by expanding approaching and receding parts of a shell or bubble surrounding the region.

With our instrument we achieve high-velocity resolution (~8 km s⁻¹) which allows us to separate the different kinematic components of the profile, and thus to detect the expansion signature on each individual pixel of a two-dimensional field.

In this paper, we describe the method, giving examples of the results and of derived parameters to showcase its potential.

In Section 2, we describe the data to which we have applied the method, giving the relevant technical specifications of the telescope and the instrument and a brief description of the reduction process. The method itself and its mathematical foundation is explained in Section 3, while in Section 4 we explain the results that can be obtained from it and the possible applications. Finally, in Section 5 we discuss the advantages and shortcomings of the method and explore possible improvements and applications.

## 2 DATA DESCRIPTION

Our code is designed to work on data cubes of galaxies with two spatial dimensions and a third spectral or velocity direction. It is optimized for a spectral line profile consisting of a single blended spectral line, to consider each resolved peak a kinematic feature, but there is no inherent problem in using it with broader spectra apart from the requirement to identify the spectral line component corresponding to each contributing peak.

The usefulness of the analysis is heavily dependent on the resolution that can be achieved with the observations, both spatial and spectral, in relation to the typical velocities of the target phenomena. In the case of bubbles or superbubbles, the velocities are generally in the range 30–150 km s⁻¹ (Chu & Kennicutt 1994; Martin 1996; Yang et al. 1996) for ionized gas and 5–50 km s⁻¹ for atomic gas (Kamphuis, Sancisi & van der Hulst 1991; Chakraborti & Ray 2011; McClure-Griffiths et al. 2002), requiring very good velocity resolution in both cases. Because of this it is clear that high-resolution instruments which take spectra in two dimension, such as Fabry-Perot detectors or aperture synthesis radio telescopes such as VLA and ALMA are ideal in supplying inputs for the kind of multipeaked profile analysis presented here.

Other methods capable of performing analysis of kinematic components, such as echelle spectrometers can take advantage of the algorithms described here, but those reach their true usefulness with two spatial dimensions, as the number of spatial elements to be analysed grows geometrically. Instruments with one spatial dimension, while also useful for these studies, do not allow a full characterization of the bubble, only along a line usually cutting through the centre of the region. Stepped longslit or long multislit spectra, without cross-dispersion, and fed by fibre bundles produce similar 3D line profiles to the Fabry-Perot described here, with the advantage of a much wider spectral range. However, the instruments available have fields which, at most, cover an area an order of magnitude smaller than a Fabry-Perot. For this reason when profiling a single emission line the latter is normally used.

An excellent example of the use of echelles in this context can be found in Meaburn’s article on 30 Doradus (Meaburn 1984).

The data used here to develop and apply the method were taken with the integral field spectrometer Galaxy Hα Fabry-Perot System (GHαFaS; see Fathi et al. 2008; Hernandez et al. 2008) mounted at the Nasmyth focus of the 4.2 m William Herschel Telescope (WHT) at the Observatorio del Roque de los Muchachos (ORM, La Palma). It produces a 3.4 × 3.4 arcmin data cube with seeing-limited angular resolution (spaxel size ~0.2 arcsec) and a spectrum at each point on the object, with 48 channels covering a spectral range which depends on the central wavelength of the emission but is usually ~400 km s⁻¹, yielding a velocity resolution of ~8 km s⁻¹ for Hα.

GHαFaS does not use an optical derotator so a correction is applied by software, the process is described in detail in Blasco-Herrera et al. (2010). Velocity calibration, phase correction, sky removal and adaptive binning are performed using the procedure described in Daigle et al. (2006) to produce the final cube which is used in our study.

An important feature of GHαFaS system is the IPCS (Image Photon Counting System) detector, which includes a CCD and a photon intensifier. The main advantage this provides, which under many circumstances more than counter balances the lower quantum efficiency, is that it has zero readout noise, an important property when studying extended features of galaxies with low surface brightness which is especially useful when considering very faint secondary peaks.
## 3 MULTIPEAKED EMISSION ANALYSIS

We aim to detect the presence of multiple components in the emission line profile for each spatial position in the data cube, so we need a way to automatically estimate and fit the components for it to be practical. Otherwise, we would have to oversee and interact with a very large number of points (for GHαFAS, with a 1024 × 1024 field it is over a million spaxels), which is impractical. To this end we have developed algorithms included in BUBBLY capable of performing this work automatically which yield reliable and precise results.

### 3.1 The second derivative method

In order to obtain the number of Gaussians for each profile and good estimates of their parameters, which is necessary to avoid local minima when fitting, we take advantage of the properties of the second derivative of Gaussians, following Goshtasby & O'Neill (1994).

A single Gaussian, or a normal distribution, $f(x)$, is defined by

$$f(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(x-\mu)^2}{2\sigma^2}}, \qquad (1)$$

where $\mu$ and $\sigma$ are the mean and the standard deviation of the distribution, respectively. The second derivative, $f''(x)$, is

$$f''(x) = \frac{1}{\sqrt{2\pi}\sigma^5} e^{-\frac{(x-\mu)^2}{2\sigma^2}} (\mu^2 - \sigma^2 + 2\mu x). \qquad (2)$$

The second derivative will take the value zero at $x_1 = \mu - \sigma$, and $x_2 = \mu + \sigma$, so, we can estimate the mean and the standard deviation as

$$\mu = \frac{x_1 + x_2}{2} \qquad (3)$$

$$\sigma = \frac{x_2 - x_1}{2}. \qquad (4)$$

This can be generalized to a function comprising $N$ Gaussians:

$$f(x) = \sum_{i=1}^{i=N} \frac{1}{\sqrt{2\pi}\sigma_i} e^{-\frac{(x-\mu_i)^2}{2\sigma_i^2}}, \qquad (5)$$

where $\mu_i$ and $\sigma_i$ are the mean and standard deviation of the ith Gaussian. The second derivative of equation (5) is

$$f''(x) = \sum_{i=1}^{i=N} \frac{1}{\sqrt{2\pi}\sigma_i^5} e^{-\frac{(x-\mu_i)^2}{2\sigma_i^2}} (\mu_i^2 - \sigma_i^2 + x^2 - 2\mu_i x). \qquad (6)$$

In order to obtain an initial value of the means and the standard deviations of the Gaussians, and to estimate the ranges of their possible values, we assume that, for any components, i, j, of the sum of Gaussians of comparable standard deviations, $\sigma_i \sim \sigma_j$, their exponential terms in equation (6) satisfy

$$\frac{1}{\sqrt{2\pi}\sigma_j} e^{-\frac{(x-\mu_j)^2}{2\sigma_j^2}} \gg \frac{1}{\sqrt{2\pi}\sigma_i} e^{-\frac{(x-\mu_i)^2}{2\sigma_i^2}} \qquad (7)$$

for $i \neq j$ and for values of $x$ close to $\mu_j$ and far enough from $\mu_i$. This implies, but only as an initial approximation from which to start the process, that the Gaussians do not overlap, i.e. that they do not contribute significantly to each other's profile. Then,

$$\mu_i \approx \frac{x_{i,1} + x_{i,2}}{2} \qquad (8)$$

$$\sigma_i \approx \frac{x_{i,2} - x_{i,1}}{2}, \qquad (9)$$

where $x_{i,1}$ and $x_{i,2}$ are the pair of zeros associated with the ith Gaussian component. This approximation clearly works better as the Gaussians are better resolved, but we want to be able to reach the limit of our spectral resolution, so we use them to estimate initial values and possible ranges used in the subsequent fit.

The limits for the mean are the values for the zeros themselves, $x_{i,1}$ and $x_{i,2}$, preventing each component from overlapping completely with the rest:

$$\mu_i^{\text{final}} \in [x_{i,1}, x_{i,2}]. \qquad (10)$$

We can also set limits to the standard deviations, but we implement them differently depending on whether the approximate value is greater than one spectral (or velocity) resolution element, $\delta$. If it is greater we let it vary between half and double the estimated value, while if it is lower we use the resolution instead:

$$\sigma_i^{\text{final}} \in \begin{cases} \left[ \frac{\sigma_i}{2}, 2\sigma_i \right] & \text{if } \sigma_i > \delta \\ \left[ \frac{\delta}{2}, 2\delta \right] & \text{if } \sigma_i < \delta. \end{cases} \qquad (11)$$

For the intensity we take the initial value as the intensity of the profile at the estimated mean and limit it to between zero intensity and 1.5 times the estimated value.

At this point we have a complete set of detected components of the profile with their respective estimated intensity, offset velocity and standard deviation, which we can feed to the fitting algorithms to obtain the final results.

### 3.2 Fitting and results

We use the properties described above to perform the detection and fitting of the Gaussian components for each spatial point. BUBBLY sweeps through the points which have sufficiently high level of emission, defined by the noise, and performs the following steps: detecting zeros of the second derivative of the profile, choosing the valid pairs of zeros, estimating the input parameters for the fit and finally performing the fit and saving the results.

A real line profile will have noise which, given its fluctuations, will produce a large number of zeros of the second derivative of the profile. There is also the possibility of an out of range peak extending partially into the range, producing an unpaired zero that would throw off the determination of all the components. To account for this, we select the valid pair of zeros by imposing a condition that the in first zero the second derivative goes from positive to negative and vice-versa for the second which is characteristic of a true peak (see Fig. 1). We also weed out all pairs between which the profile does not reach an intensity above a threshold signal-to-noise ratio (SNR).

Once we have the final choice of zero pairs we estimate the initial values and ranges as explained in the previous section and feed them to the algorithm for multiple Gaussian fitting. Our algorithm has been constructed using modified versions of functions specifically developed to fit Gaussians created by Adam Ginsburg for his GAUSSFITTER PYTHON program. This program is a wrapper with multiple options (such as upper and lower limits to the parameters) to the MPFIT module for Levenberg-Marquardt least-squares fitting.<fnref n="1" />

As a result, we obtain multicomponent line emission surface brightness, velocity and velocity dispersion maps. Since the process is done individually on each spatial point the number of detected components will vary over the field.

<footnote n="1">¹ Both GAUSSFITTER and MPFIT PYTHON packages can be found at http://code.google.com/p/agpy/</footnote>
<figure><img src="image_3.png" /></figure>

**Figure 1.** Graph illustrating the effectiveness of the second derivative method to separate Gaussian components. The points correspond to the zeros of the second derivative. Note that the centre of any Gaussian lies between two zeros of the second derivative, a positive to negative zero and a negative to positive zero. This holds true regardless of the direction in which you sweep the profile to search for these zeros and is used to find and determine the initial parameters for each Gaussian component.

These maps can be highly useful for many applications, in this paper we will focus on superbubbles, but there have been a number of different types of studies requiring the analysis of multiple emission peaks. A particular topic of interest is the detection and characterization of Seyfert galaxies based on the presence of multiple components in the emission line profiles. In Arribas, Mediavilla & Garcia-Lorenzo (1996), García-Lorenzo (2013) the double-peaked emission line profiles of NGC 1068 are sought and used to analyse the active nucleus of the galaxy, while Arribas, Colina & Clements (2001) perform a similar study on Arp 220. Liu et al. (2010), Ge et al. (2012) and Pilyugin et al. (2012) explore the possibility of using multiple peak analysis on a sample of Sloan Digital Sky Survey (SDSS) galaxies, where automated methods similar to ours would be of great value. Another application of the presence of multiple peaks in data cubes is found in Williams, de Geus & Blitz (1994), where a routine to separate the blended emission of molecular gas clumps is presented. Yet another example is the detection of hydrogen gas transfer between a pair of interacting galaxies in Font et al. (2011).

As an example of the application of our method in particular, in Zaragoza-Cardiel et al. (2013) a local inspection of the multicomponent velocity maps around the nuclear starburst in NGC 3396 was made. The authors detected, using the secondary components of the velocity, an outflow of ionized gas.

# 4 APPLICATION TO EXPANDING STRUCTURES

We can easily apply the results from the multipeak fit to search for expanding structures. The kinematics along the line of sight for an expanding shell around an H II region present three peaks, one in the middle produced by the bulk of the region, one for the gas approaching us because of the expansion and a corresponding peak receding from us emitted by the far side of the shell. Because of this we define the expansion signature that we will seek in the profiles as a bright main emission peak, coming from the bulk of an H II

<figure><img src="image_4.png" /></figure>

**Figure 2.** Example of expansion detection in an H II region of the Antennae (see Fig 3). The individual components (green, red and blue) are plotted along with the profile (black) and the fit result (crosses). The brightest, main peak corresponds to the emission from the bulk of the H II region, and the two smaller, symmetrically placed peaks on each side correspond to the approaching and receding parts of the superbubble along the line of sight in that particular spatial point. Further peaks can be observed on the profile, but they are not taken into consideration here as they are not symmetric in velocity about the central peak, which they would be in the case of expansion. They represent a complexity of kinematic structure which is clearly present, but which we are, for the present, not pursuing in detail.

region, and two secondary peaks placed symmetrically, one at each side of the main peak. An example of an expansion signature from a superbubble is shown in Fig. 2.

The inhomogeneity of the ISM, the effect of the density profile of the gas decreasing perpendicularly to the disc and instrumental error will generally cause these two velocities to differ slightly from perfect symmetry, so we allow a maximum deviation from symmetry with respect to the central peak of twice the spectral resolution of our instrument, $2\delta v \sim 16 \text{ km s}^{-1}$. Tests were performed comparing the results obtained using different allowed deviations. Initially, we used a deviance of $\delta v$, the instrumental error for these measurements, but we found that the $2\delta v$ maps contained the same features as the $\delta v$ maps, i.e. no additional macroscopic features, and the structures detected for each feature appeared more complete. Further tests using $3\delta v$ did not yield more coherent structures but did introduce false signatures, so the $2\delta v$ window was used throughout this work.

Once it has the multiple component maps BUBBLY determines the main peak of the line profile for each spatial point, using flux as the criterion and then checks the other detected peaks to find a pair that satisfies the conditions imposed. It then creates a map containing the mean separation between these peaks (the expansion velocity at that point) and another map with their mean intensity. For the example shown here we used instead the relative intensity with respect to the total intensity in order to use fluxes obtained from narrow-band images ($I_b = 0.5*(I_1/I + I_2/I)$). BUBBLY offers the option of imposing additional conditions on the secondary peaks to accept them as an expansion signature, namely having similar intensities and velocity dispersions for the symmetric pairs.

Another result of the overall mapping process is a velocity map of the main peak, which gives the overall kinematics of the galaxy. This map is an alternative to the usual method for
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison between the Hα narrow-band image (left) of the Antennae galaxies and the corresponding expansion map (right). The circles represent the estimated bubble radii and the white line signals the position for the Hα profile shown in Fig. 2. The scale is in velocity with respect to the emission from the bulk of the HII region.</figcaption></figure>

obtaining a kinematic map in Hα which uses the moments method. The latter calculates the moments of the line profile to find the central ‘characteristic’ velocity, while our algorithm takes the velocity of the brightest peak of each profile. In general, our algorithm gives a more reliable result, because if there are secondary emission peaks these will weight the moments method and give a distorted value for the velocity of the spaxel concerned. However, in kinematically complex environments, where there may be a number of peaks with comparable intensities, one can obtain artificial jumps in the velocity map using our algorithm, because it selects the brightest of these peaks, regardless of whether this represents the overall underlying velocity, while the moments method finds an average, and smooths the differences. In this section, we have included an example of the application of BUBBLY to a data cube and the results that originate from it to demonstrate the detection capabilities of the program. The selected object is the Antennae pair of galaxies, which exhibits very interesting results on superbubbles originating from powerful current star formation. A paper analysing the superbubbles in depth and comparing the detections to previous work on them is in preparation.

## 4.1 Finding the expanding structures

After processing the data cube with BUBBLY the main products we obtain are the expansion velocity map, the intensity map and velocity, the intensity and dispersion maps for the main peak. The first two are the relevant ones for the study of superbubbles, as they are what we will use to detect and characterize them. In Fig. 3, we show the expansion map of the Antennae galaxies and its corresponding Hα narrow-band image as an example of the products of our method. It can be seen in the image that the centres of the superbubbles detected coincide well with the centres of the brightest HII regions. A paper analysing the properties of the superbubbles in this object will be submitted in the near future.

We detect the superbubbles by inspecting the expansion map to find regions with coherent patterns in the expansion velocity that

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Detail of a superbubble detected in the brightest HII region in the Antennae galaxies. The scale is velocity with respect to the emission from the bulk of the HII region and the white circle represents the estimated superbubble radius. As expected from an expanding structure, we can see the highest velocity at the centre, falling to lower values towards the edge. The spatial coherence and effect of projection on the observed velocity can be well seen in this example.</figcaption></figure>

correspond to the expected radial distribution of observed expansion velocity characteristic of superbubbles. These arise from the fact that the velocities we measure are the components directed along the line of sight due to the Doppler effect, so a spherically expanding shell will show projection effects on the expansion map. Therefore, observed expansion velocity will be highest at the centre of the bubble where the velocity vector aligns with the line of sight, decreasing radially until it is zero at the edge, so we look for radially decreasing circular patterns in the map.

A particularly clear case of a detected superbubble from the Antennae is shown in Fig. 4 with the estimated bubble size and position shown by a white circle. This example is the clearest one
on the object, because usually we do not detect expansion in the very centre of the shell (which is typically optically thin) but only an annulus concentric with the H II region.

The absence of observable expansion components in the central zone of a superbubble is due to a combination of phenomena. One is that superbubbles of this size enter a blowout phase as they grow beyond the vertical density scaleheight of the disc, so that the gas of the central part thins out as it expands freely into the halo; another is the extinction caused by dust which may obscure the receding component of the bubble signature, as that emission has to travel through the entire region to get to us and the centre of the far face of the superbubble is the part with the largest optical path to the observer.

Both of these factors are combined with the aforementioned fact that the observed column density for a spherical shell is lowest at the centre because of the projection effect.

Naturally, the assumption of sphericity is inexact at best. There have been many theoretical and simulation studies of the evolution of such superbubbles (Mac Low & McCray 1988; Zaninetti 2004), all of them concluding that they soon acquire an elongated shape, with the long axis perpendicular to the plane of the galaxy, as the ambient density decreases in that direction, morphing into a peanut shape for bigger and older superbubbles. The projection effects described still apply, but it can be expected, for example, that the observed projected expansion velocity will be affected. As the superbubble has a peanut shape, its observed edge as seen from above will not have a zero value of the vertical velocity, so we still see the gradient, but the expansion velocity will decrease more slowly and we will stop detecting the bubble before it approaches zero. The expansion velocity profiles we have found do agree better with the elongated case.

### 4.1.1 Multiple expanding features

In some cases around the brightest H II regions, profiles exhibiting multiple pairs of symmetrically separated peaks can be observed (Fig. 5). These additional features follow the same general shape

<figure><img src="image_7.png" /><figcaption>Figure 5. Profile observed in the Antennae galaxies where two pairs of symmetrically spaced secondary peaks are detected, implying the possible presence of two superbubbles in one region. An alternative possibility is that it could be a projection effect with two H II regions overlapping on the line of sight which may indicate more complex kinematics. As in Fig. 2, we show separated peaks only for detected components that are symmetrical to the main peak.</figcaption></figure>

as the main detected superbubbles, forming rings around the H II regions, but they are not as common and the analysis and determination of their properties is more difficult.

However, they represent a very interesting possibility as their presence may imply multiple bursts of star formation in the same regions or perhaps a system of more local bubbles inside the large region and its associated superbubble. Another possibility to explain the presence of multiple expansion signatures in a spaxel would be that we are detecting kinematics which are considerably more complex than that of an assumed superbubble.

Deeper spectral imaging should help to characterize these and other complex features.

## 4.2 Bubble properties

The amount of information about a bubble or shell that we can extract is critical for our method to be truly useful. We have four primary quantities from which we may derive the physical parameters of a particular bubble: the radius, the H II region luminosity, the shell-to-region luminosity ratio, and the expansion velocity of the shell. We can calculate the bubble luminosity, density, mass, kinetic energy and age. The basic equation we use is the luminosity of a shell (Relaño & Beckman 2005)

$$L_{H\alpha}(\text{shell}) = 4\pi R^2 \Delta R n_e^2 \alpha(H, T) h\nu. \quad (12)$$

Among the parameters appearing in this equation are the effective recombination coefficient for hydrogen $\alpha_{H\alpha}^{\text{eff}}(H, T)$, which can be taken as $4 \times 10^{-13}$ for H II regions (Osterbrock 1989), $h\nu_{H\alpha}$ is the energy of an H$\alpha$ photon and the shell width $\Delta R$ which we take as ~2pc. We do not expect the width to vary much from case to case and when calculating the mass and subsequent quantities its contribution varies as $\Delta R^{-1/2}$, so it does not affect the final result very strongly. Observational uncertainties give a much larger possible variation so we can safely consider it as constant.

Since we already know the luminosity of the shell from the product of the region luminosity and the shell-to-region luminosity ratio we are especially interested in the density to calculate the total mass of the shell and its energy. We can also make a zero-order estimate of the age of the bubble by dividing the radius by the expansion velocity, which is equivalent to assuming the region has expanded at constant velocity. Thus

$$n_e = \sqrt{\frac{L_{H\alpha}(\text{shell})}{4\pi R^2 \Delta R \alpha_{H\alpha}^{\text{eff}}(H, T) v_{H\alpha}}} \quad (13)$$

$$M = 4\pi R^2 \Delta R n_e m_p \quad (14)$$

$$E_k = \frac{1}{2} M * v_{\text{exp}}^2 \quad (15)$$

$$\text{Age} = R/v_{\text{exp}}. \quad (16)$$

In Table 1, there is a list of these measured parameters for some of the superbubbles we have detected in the Antennae. The errors for the parameters are determined by propagation of the errors in radius, expansion velocity and H$\alpha$ luminosity. The first two parameters are estimated by us directly by inspection of the expansion map, so we took a conservative error of 10 per cent for both radius and expansion velocity. For the luminosity the most important source of error is, by far, the uncertainty in the radius of the region, though it does
<table><caption>**Table 1.** Properties of some detected bubbles in the Antennae galaxies.</caption><thead><tr><th>Radius<br>(pc)<br>(1)</th><th><i>v</i><sub>exp</sub><br>(km s<sup>-1</sup>)<br>(2)</th><th><i>L</i><sub>Hα</sub>(shell)<br>10<sup>38</sup> erg s<sup>-1</sup><br>(3)</th><th><i>n</i><sub>e</sub><br>(cm<sup>-3</sup>)<br>(4)</th><th>Mass<br>(10<sup>35</sup> Kg)<br>(5)</th><th>Kinetic energy<br>(10<sup>51</sup> erg)<br>(6)</th><th>Age<br>(Myr)<br>(7)</th></tr></thead><tbody><tr><td>334 ± 33</td><td>115 ± 12</td><td>3.04 ± 0.25</td><td>6.4 ± 0.7</td><td>7.36 ± 1.5</td><td>48.7 ± 13.6</td><td>2.8 ± 0.4</td></tr><tr><td>304 ± 30</td><td>82 ± 8</td><td>0.44 ± 0.06</td><td>2.7 ± 0.3</td><td>2.55 ± 0.5</td><td>8.6 ± 2.6</td><td>3.6 ± 0.5</td></tr><tr><td>319 ± 32</td><td>73 ± 7</td><td>1.83 ± 0.21</td><td>5.2 ± 0.6</td><td>5.45 ± 1.2</td><td>14.5 ± 4.5</td><td>4.3 ± 0.6</td></tr><tr><td>415 ± 42</td><td>82 ± 8</td><td>0.77 ± 0.10</td><td>2.6 ± 0.3</td><td>4.59 ± 0.9</td><td>15.4 ± 4.4</td><td>4.9 ± 0.7</td></tr><tr><td>269 ± 27</td><td>70 ± 7</td><td>0.16 ± 0.02</td><td>1.8 ± 0.2</td><td>1.36 ± 0.3</td><td>3.3 ± 0.9</td><td>3.8 ± 0.5</td></tr><tr><td>383 ± 38</td><td>57 ± 6</td><td>0.80 ± 0.10</td><td>2.9 ± 0.3</td><td>4.33 ± 0.9</td><td>7.0 ± 2.0</td><td>6.6 ± 0.9</td></tr></tbody></table>

*Notes.* All parameters refer solely to the bubble itself, not the underlying H II region. The kinetic energy (7) is given in units of $10^{51}$ erg, which is approximately the energy one supernova injects into the medium. In this manner, the energy can be taken to be in units of supernovae necessary to drive the expansion in that particular region. The age (8) is calculated as $R(2)/v_{exp}(3)$.

not propagate in a mathematically straightforward way as it affects the pixels over which we sum the detected flux and depends on the physical characteristics of each region. What we do is measure three sets of fluxes for each region, two of which correspond to the high and low extreme values of the radius as allowed by the radius error and take their mean deviation from the central value as the error for the flux.

### 4.2.1 Age dating

A precise determination of the age of a superbubble is very important as it is equivalent to the age of the associated H II region and can therefore be used to give a short-term map of the star-forming history of the galaxy. In Table 1, we list ages for each superbubble, but these are merely intended as rough approximations given the crude way in which they are calculated, as simply $R/v_{exp}$. The superbubble evolves as a complex interaction of the winds, supernovae and ionization pressure with the non-isotropic ISM, and the expansion velocity will vary throughout its lifetime.

A better approach to the issue is therefore necessary and entails models of the interaction. We have experimented with analytic models for the expansion of superbubbles with moderate success; comparison with observations of line equivalent widths in clusters to produce their stellar ages gives compatible ages in most cases, especially for the younger regions. The models assume homogeneous distribution of gas in all directions, while in reality and on the scales we are studying, it becomes necessary to take into account the 3D shape of the distribution. We are working with detailed physical simulations to provide a better determination of the ages and other parameters.

## 5 DISCUSSION

We have presented our method and associated software designed to detect and fit multiple components to integral field spectroscopy of emission lines, with an application to superbubbles using our Fabry-Perot instrument.

In Section 3, we have explained in detail the method for the line profile decomposition into a set of Gaussians and the way it is implemented in *BUBBLY*. The next section applied the results from the decomposition to produce expansion maps that can be used to search for and characterize superbubbles. In this regard with a relatively simple approach we already obtain reliable detection and characterization of superbubbles.

The presence of multiple components in an emission line is a powerful way to study the fainter and more difficult to detect, kine-

matic phenomena in a galaxy, as well as numerous other applications such as separating the emission of different regions or clouds or even galaxies that overlap or AGN and violent galactic outflows (see Section 3.2). Our algorithm has the advantage of being both simple and easy to use and very flexible with regards to the input data, which makes it a potent tool to study galaxy kinematics.

In this paper, we have centred on the application of multicomponent analysis to the search and characterization of superbubbles, where the method proves especially useful allowing us to detect several superbubbles clearly associated with giant H II regions. The fact that it does not need purpose-programmed observations to work means that it offers added value to integral field spectroscopy observations of galaxies.

Current alternative methods for the detection of superbubbles are to use morphological features observed through narrow H α filters such as cavities and arches around OB associations, shock-fuelled forbidden emission lines such as S II and O III and morphological and kinematic features in H I. X-ray observations of extended emission interior to the observed shells is a complementary diagnostic.

All morphology-based methods require very high spatial resolution and shells which are bright with respect to their surroundings in order to detect them, which restricts the methods to nearby galaxies and/or observations with very high angular resolution. This is not the case for our method, where the shell emission is separated from the H II region emission kinematically, although of course with better angular resolution we will be able to detect more and smaller bubbles.

Except for H I cube observations the other methods do not provide kinematics of the shell which prevents a detailed characterization of the bubble and its dynamic properties, limiting the utility of the data. It would be possible to complement the observations with high-resolution spectra to detect the expanding components but the total observing time required would be very high, considering that each region would require an individual spectrum.

The most direct competition for our method when it is applied to superbubbles is, therefore, observing H I using the classical 21 cm method. The long wavelength of this H I emission makes it very difficult to achieve the spatial resolution necessary to distinguish these features properly, (even though the velocity resolution is high), making it very difficult to extend these studies beyond the closest galaxies in the Local Group, such as M101 where several shells have been detected in this manner (Kamphuis et al. 1991; Chakraborti & Ray 2011).

In any case the shells detected with H I normally correspond to the late stages of the superbubble, after the massive stars which
ionized the shell have disappeared and the superbubble has decelerated greatly. Because of this the typical expansion velocities lie in the range 5–50 km s⁻¹ (Kamphuis et al. 1991; McClure-Griffiths et al. 2002; Chakraborti & Ray 2011) and in many cases it is difficult or impossible to find the OB association corresponding to the superbubble. This is not the case for our method which traces the superbubbles from the ionized gas, so we find young superbubbles clearly associated with the cluster. In this manner, it can be said that H₁ observations and ours are complementary, each detecting superbubbles in different stages of their evolution.

Of course, the method for multipeak detection can be applied to H₁ data cubes without problems, but to detect superbubbles we would have the problem of lacking a main peak associated with the emission of the region and therefore a reference velocity for the secondary peaks. This could be solved using the velocity field of the galaxy for comparison, which is in fact a normal product of the analysis of these data.

It would certainly be possible to adapt our code to be run on H₁ data cubes to apply the method and yield the expansion maps.

We are limited in the range of application by the faint nature of the secondary components we try to detect and the spatial resolution necessary to see the structure of the superbubble, which is required in order to decouple the shell signature from complex non-expanding kinematic structure along the line of sight.

The first problem prevents expansion detection at some points where one of the components is not observed, though careful examination of the profiles of the regions clearly shows that in most of these points one of the pair of secondary peaks that make up the expansion signature appears, or even both, with one of them having a lower intensity and falling below the minimum SNR required by the program. These appear at the same velocities as in the points with double detection, so it is evident that they correspond to the same expanding structure.

The second problem (spatial resolution) impedes a more complete and statistical study of all detected expansions, as we need to see the characteristic projection structure of the velocity. Otherwise we could distinguish real signatures and their properties at much smaller scales, yielding many more bubbles/shells per galaxy and allowing a much more detailed study of the feedback effects within the galaxy.

The values for the parameters we obtain (see Table 1 for some examples) appear to be reasonable and consistent among the regions. The results should be used with some care, however, as they rely on several simplifications such as assuming a perfect spherical shell and a homogeneous medium, which we know to be only moderate approximations. We expect these deviations to be important but not critical, so the properties we have obtained should be an acceptable approximation for many purposes.

By far the biggest source of uncertainty when deriving these parameters is the radius of the shell as it influences strongly the determination of the flux and it also appears as an important parameter for the calculations of the physical parameters. This is because the determination is performed by inspection and suffers from a lack of spatial coherence in many cases, so we chose a conservative margin of error to be safe. The value of the radius is also subject to some uncertainty because the true shape of the superbubble is non-spherical. We expect superbubbles to be elongated in the direction perpendicular to the disc, and therefore, we are usually measuring the lesser radius of an ellipsoid. This would mean we are underestimating the kinetic energy by a factor depending on how elongated the structure is, but another effect acts in the opposing sense. In an elongated superbubble the expansion along the major axis (perpendicular to

the disc) is faster than along the minor axis (parallel to the disc), this means that the expansion velocity we use is an upper limit and thus we are overestimating the kinetic energy. The combination of the two effects should tend to cancel somewhat, but in general younger and smaller superbubbles will be more spherical and thus provide more reliable energy and age determinations.

We are working to use hydrodynamic simulations to constrain these effects and obtain more reliable parameter determinations.

## 5.1 Further developments and applications

Throughout the text the current shortcomings of the method and possible developments to address them and further improve the results are commented. Here is a summary of the main possible developments and applications.

(i) Write additional modules to take advantage of the multiple component emission detected in order to search for wind streams, inflows and outflows within a galaxy.

(ii) Improve the detection capabilities of the code by implementing pattern recognition algorithms in order to correlate the signatures between themselves and detect lower emission signatures, allowing us to better characterize bubbles/shells and to expand the study to smaller scales.

(iii) Obtain a precise determination of the age of the expanding structure and, subsequently, the underlying cluster by means of hydrodynamic simulations.

(iv) Explore the possibility of the presence of multiple expansions within a region and the implications with regards to the conditions determining star formation.

(v) Apply results obtained from our method to study feedback and evolution in galaxies.

The feasibility of this program is based on the effectiveness of BUBBLY as a superbubble detector.

## ACKNOWLEDGEMENTS

This research has been supported by the Spanish Ministry of Economy and Competitiveness (MINECO) under the grant AYA2007-67 625-CO2-01, and by the Instituto de Astrofísica de Canarias under project P/308 603. JEB acknowledges financial support to the DAGAL network from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/ under REA grant agreement number PITN-GA-2011-289 313.

The WHT is operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.

We thank the referee, John Meaburn, for very helpful comments which have enabled us to make a significant improvement in the paper.

## REFERENCES

* Antonucci R., 1993, ARA&A, 31, 473

* Arribas S., Mediavilla E., Garcia-Lorenzo B., 1996, ApJ, 463, 509

* Arribas S., Colina L., Clements D., 2001, ApJ, 560, 160

* Arsenault R., Roy J.-R., 1986, AJ, 92, 567

* Bally J., Lada C. J., 1983, ApJ, 265, 824

* Bland J., Tully B., 1988, Nature, 334, 43

* Blasco-Herrera J. et al., 2010, MNRAS, 407, 2519

* Blitz L., Stark A. A., 1986, ApJ, 300, L89
Castor J., McCray R., Weaver R., 1975, ApJ, 200, L107

Ceverino D., Klypin A., 2009, ApJ, 695, 292

Chakraborti S., Ray A., 2011, ApJ, 728, 24

Chu Y.-H., 2008, in Bresolin F., Crowther P. A., Puls J., eds, Proc. IAU Symp. 250, Massive Stars as Cosmic Engines. Cambridge Univ. Press, Cambridge, p. 341

Chu Y.-H., Kennicutt R. C., Jr, 1986, ApJ, 311, 85

Chu Y.-H., Kennicutt R. C., Jr, 1994, ApJ, 425, 720

Chu Y.-H., Troland T. H., Gull T. R., Treffers R. R., Kwitter K. B., 1982, ApJ, 254, 562

Daigle O., Carignan C., Hernandez O., Chemin L., Amram P., 2006, MNRAS, 368, 1016

Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431

Dyson J. E., de Vries J., 1972, A&A, 20, 223

Fathi K. et al., 2008, ApJ, 675, L17

Font J. et al., 2011, ApJ, 740, LL1

García-Lorenzo B., 2013, MNRAS, 429, 2903

Ge J.-Q., Hu C., Wang J.-M., Bai J.-M., Zhang S., 2012, ApJS, 201, 31

Gerola H., Seiden P. E., 1978, ApJ, 223, 129

Goshtasby A., O'Neill W., 1994, Graph. Models Image Process., 56, 281

Heckman T. M., Armus L., Miley G. K., 1990, ApJS, 74, 833

Heiles C., 1979, ApJ, 229, 533

Hernandez O. et al., 2008, PASP, 120, 665

Kamphuis J., Sancisi R., van der Hulst T., 1991, A&A, 244, L29

Lada C. J., 1985, ARA&A, 23, 267

Leitherer C. et al., 1999, ApJS, 123, 3

Liu X., Shen Y., Strauss M. A., Greene J. E., 2010, ApJ, 708, 427

Mac Low M.-M., McCray R., 1988, ApJ, 324, 776

McClure-Griffiths N. M., Dickey J. M., Gaensler B. M., Green A. J., 2002, ApJ, 578, 176

McCray R., Snow T. P., Jr, 1979, ARA&A, 17, 213

McCray R., Kafatos M., 1987, ApJ, 317, 190

Martin C. L., 1996, ApJ, 465, 680

Meaburn J., 1970, Nature, 228, 1069

Meaburn J., 1980, MNRAS, 192, 365

Meaburn J., 1984, MNRAS, 211, 521

Nazé Y., Chu Y.-H., Points S. D., Danforth C. W., Rosado M., Chen C.-H. R., 2001, AJ, 122, 921

Osterbrock D. E., 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei. University Science Books, Mill Valley, CA

Palous J., Tenorio-Tagle G., Franco J., 1994, MNRAS, 270, 75

Pilyugin L. S., Zinchenko I. A., Cedrés B., Cepa J., Bongiovanni A., Mattsson L., Vílchez J. M., 2012, MNRAS, 419, 490

Relaño M., Beckman J. E., 2005, A&A, 430, 911

Relaño M., Beckman J. E., Daigle O., Carignan C., 2007, A&A, 467, 1117

Sharma S., Corbelli E., Giovanardi C., Hunt L. K., Palla F., 2011, A&A, 534, A96

Stutzki J., Guesten R., 1990, ApJ, 356, 513

Weaver R., McCray R., Castor J., Shapiro P., Moore R., 1977, ApJ, 218, 377

Williams J. P., de Geus E. J., Blitz L., 1994, ApJ, 428, 693

Yang H., Chu Y.-H., Skillman E. D., Terlevich R., 1996, AJ, 112, 146

Zaninetti L., 2004, PASJ, 56, 1067

Zaragoza-Cardiel J. et al., 2013, MNRAS, 432, 998

This paper has been typeset from a *TEX/LATEX* file prepared by the author.