HALOGAS: H i OBSERVATIONS AND MODELING OF THE NEARBY EDGE-ON SPIRAL GALAXY NGC 4244

The simplest model involves a single disk with a slight warp perpendicular to the line of sight (deviating by 4°–5° from the inner disk (Figure 8) and beginning at a radius of 12 kpc). For such a model, all other quantities aside from column densities and rotational velocities are held constant with radius. The inclination is found to be 88° for the best fit to the vertical profile. For the final model, the scale height is about 600 pc. This simple model (referred to as "1C" in figures) cannot be immediately eliminated based on the fit to the vertical profile and zeroth-moment map (Figure 9). However, one may note that in the bv plots in Figure 6, the one-component model lacks sufficient flux on the systemic relative to the terminal side in each panel. As can be seen in the channel maps in Figure 7, shown most prominently in the plots closest to the systemic velocity and noted with an arrow, the edge of the diagram away from the center of the galaxy comes to a definitive point at a positive minor axis offset on the approaching half, and a negative offset on the receding, which is not seen in the data. No other combination of inclination and scale height improves this.

Adding a projection of the warp along the line of sight (referred to as "W" in figures), by gradually decreasing the inclination in the outermost rings (shown in Figure 8), matches the observed "T" shape from the data in the central bv diagrams in Figure 6 significantly better than the 1C model. This is the result of bringing flux from outer rings with lower projected velocity into higher latitude regions. If the inclination decreases too quickly, the "T" shape becomes excessive.

In contrast to the one-component model, the channel maps in Figure 7 show that the flux is pushed upward and at high z is pushed away from the center in the major axis direction, in panels with velocities close to the systemic velocity in the warp model, which matches the data more closely. However, if too much of a warp is present in the model, then analogous to the "T" shape verging on becoming too "Y"-shaped in the bv diagrams, there will be an indentation on this edge. This is not seen in the data.

The best-fit line of sight (LOS) warp component deviates by 8° from the inner disk and as in the position-angle warp, starts at a radius of 12 kpc (Figure 8).

Upon initial inspection of the zeroth-moment maps, it may be noted that the edges in the outer radii taper significantly in minor axis extent. Because of this, a substantial flare can be ruled out immediately as it would cause the zeroth-moment map to become too boxy and fail to fit this aspect of the data. However, this is not to say for certain that a modest flare cannot be present. For this reason, a representative flaring model (referred to as "F" in figures) is created by starting with the one-component model and increasing the scale height outward radially. This increase begins at the optical radius (10.4 kpc; Table 1) and increases linearly, approximately doubling by a radius of 15 kpc. Comparison of zeroth-moment maps from the data (Figure 9) and model indicates that a modest flare yields a slight improvement compared to the one-component model. Unlike the warp model, the modest flare does not eliminate the rounding on the systemic side of the bv plots (Figure 6) or the thinning on the outer edges of the channel maps seen in the one-component model. A more extreme and less physical flare (abruptly doubling the scale height at the optical radius and then holding it constant at larger radii) helps to reduce these problems in the bv plots and channel maps, but once again, causes the zeroth-moment map to become boxy and an obvious mismatch to the data (not shown).

The pure flare model clearly does not match the data in the bv plots. However, adding the same modest flare to the warp model (referred to as "WF" in figures) does not appear to have a substantial effect; in fact, it is almost identical to the warp model as shown in the zeroth-moment map in Figure 9 as well as the bv plots in Figure 6. This indicates that such a flare could be present in NGC 4244, but that the component of the warp along the line of sight is dominant and effectively hides most of its contribution. Since we are unable to determine whether or not the flare exists due to the minimal contribution it would have compared to the warp, and for the sake of clarity, it is omitted from further models.

H i halos have been observed in several galaxies, with NGC 891 among the most notable (Oosterloo et al. 2007). One of the primary goals of HALOGAS is to determine the frequency of occurrence and origins of such halos. We therefore examine the possibility of such a halo in NGC 4244.

The addition of a halo to the model is made via combining a second, thicker exponential component to the first, keeping the inclination of both at 88° (referred to as "2C" in figures). Upon inspection of the vertical profile and zeroth-moment map, there is no immediate motivation for a second component, as reasonable fits can be made with less complicated models. For both halves of NGC 4244, multiple ratios of scale height and flux in each component were tested in attempts to match the vertical profile, but none was successful. In fact, when considering only the vertical profile, a two-component model cannot be made to match the data, as during the modeling process, the components tend to converge to a single scale height, resulting in a one-component model.

However, since the existence of a halo is of primary concern, we have made substantial attempts to optimize this model. Since the vertical profile and zeroth-moment map provide no useful constraints for the relative scale heights of the two components or for the percentage of flux contained within the halo, we use illustrative fixed values, with the scale height of the extended component twice that of the thin component, and with 25% of the flux residing in the extended component. These values are chosen because they represent a reasonable scenario compared to other galaxies with halos (Oosterloo et al. 2007).

Through optimization of this two-component model, it is found that a reasonable fit to the data is impossible without including a warp component along the line of sight. A halo will only add to the complexity of the model and if a warped model without a halo is sufficient (or in this case, slightly better), then the simpler model should be used. Thus, we reject the two-component model with or without such a warp. We do however retain a two-component model with a warp component along the line of sight ("2CW") in figures throughout this paper for illustrative purposes only.

At this point, prior to the addition of a lag, we place our confidence in the single-component line of sight warp model. However, as may be seen in Figures 6 and 7, the optimized models with no lag do not exactly match the data. Adding a lag at this stage will best illustrate its contribution to the models. We continue to include the one-component model, the flaring model with a warp component along the line of sight, as well as the two-component model with a warp component along the line of sight for illustrative purposes as well as to demonstrate the reliability of the result concerning the lag.

The columns in Figure 6 denoted with "L"s in the titles show the models as before, but with optimal lags for each model. Upon inspection of these plots, it becomes evident that the addition of a lag is an almost universal improvement. By adding a lag, flux in the models is drawn to the systemic side, thus duplicating the aforementioned "T" and "V" shapes in the data. Furthermore, the edge of the systemic side becomes flatter in almost all of the panels as in the data, whereas in the no lag models, only the models including a line-of-sight warp component appear relatively flat. Additionally, the channel maps are improved. They narrow in the z-direction at velocities further from systemic, more closely matching the data.

The warped model with the addition of a lag of −9⁺³ _{− 2} km s⁻¹ kpc⁻¹ and −9  ±  2 km s⁻¹ kpc⁻¹ in the approaching and receding halves, respectively, is a reasonable fit. Error bars are estimated via visual inspection of the range of lags which could potentially match the data and under the (inexact) assumption that the models have no additional uncertainties (e.g., uncertainty in column density or rotational velocity), as making such estimates would be unwieldy.

For the two-component warp model, a lagging halo is simulated by adding a lag in only the thick component. The addition of the lag substantially improves the model for the receding half with a lag as steep as −9 km s⁻¹ kpc⁻¹, but fails to adequately duplicate the "V" shape in Figure 6 in the approaching half. Even with a lag as steep as −28 km s⁻¹ kpc⁻¹ this shape cannot be duplicated. Given this difficulty, a lag of −9 km s⁻¹ kpc⁻¹ is used in the final model to optimize the fit for the overall shape of the data.

The improvement due to the lag is present in all models. This attests to the robustness of the lag in that, regardless of the optimal morphological model, there are features which cannot be re-created without it. Once again, the one-component warp model with the addition of a lag appears to be the best match to the data.

Upon establishing the existence of a lag, it must be characterized by not only the magnitude of the gradient in the z-direction but also its radial variation. Using the one-component warped model, bv diagrams farther from the center of the galaxy than those shown in Figure 6 are examined in Figure 10 using a 30'' × 30'' cube, with the range of slice locations corresponding to the dashed lines in Figure 2. The data, at these large radii, are best fit with a shallowing lag, decreasing in magnitude to −5 ± 2 km s⁻¹ kpc⁻¹ and −4 ± 2 km s⁻¹ kpc⁻¹ near a radius of 8' (10 kpc) in the approaching and receding halves, respectively. This is most easily seen on the terminal side of each panel: the data contours become more rounded, whereas the model with a constant lag retains too much of a "V" shape, with data contours extending farther on the terminal side than those of the model. However, on the negative offset side of the receding half, the shallowing of the lag may be even more substantial. This is best seen in panels corresponding to −82 and −9'. Uncertainties are still quantified based on the positive offset side, and −4 km s⁻¹ kpc⁻¹ is only given as an upper limit for the side with negative offset. There is no strong indication of any radial variation at radii smaller than approximately 6' (7.5 kpc). It should be noted that the signal-to-noise ratio (S/N) of the data is considerably diminished at and beyond a radius of 10' and determining the lag this far from the center approaches the limits of our modeling capabilities. The error estimates are once again by eye.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Although the lag in the H i layer of NGC 4244 does not fit the trend with star formation rate seen by Heald et al. (2007) for DIG layers, we can ask whether such a trend is apparent within NGC 4244: the lag is comparable on the two sides of the galaxy, but is this also true of the star formation rate? If the Heald et al. (2007) trend continues, then a large contrast in star formation rate would be unexpected. We consider the star formation in each half separately by using Hα (Hoopes et al. 1999) and 24 μm MIPS data (NASA/IPAC Infrared Science Archive). Using these, we apply the method described in Kennicutt et al. (2009) to correct for extinction. We estimate the uncorrected Hα flux in each half by measuring the percentage of the total Hα flux in each and then using the total Hα luminosity for NGC 4244 found in Kennicutt et al. (2008) of 10⁴⁰ erg s⁻¹. Accounting for the slight difference in distances used, these estimates are 5.4 × 10³⁹ and 4.4 × 10³⁹ erg s⁻¹ in the approaching and receding halves, respectively. (This ordering will remain consistent throughout this section.) Assuming a monochromatic flux (L = λL_λ) for the 24, 70, and 160 μm data, the total infrared (TIR) fluxes are found to be 1.1 × 10⁴² and 1.0 × 10⁴² erg s⁻¹. By assuming a Salpeter initial mass function and then applying a conversion between line fluxes and star formation rates supplied by Kennicutt et al. (2007), the star formation rate for each half may be determined:

$$\begin{equation} {\rm SFR}_{{M}_{\odot }\, {\rm yr}^{-1}}=9\times 10^{-42}[L_{\rm H\alpha }+aL_{\rm TIR}](\hbox{erg} \hbox{ s}^{-1}), \end{equation} \tag{ 1 }$$

where a  =  0.0024. These are found to be 0.064 and 0.054 M_☉ yr⁻¹. The star formation rate in the approaching half is approximately 15% higher than that in the receding half. However, the uncertainties in both the lag magnitude and the local star formation rates do not allow for an unambiguous conclusion regarding the trends seen by Heald et al. (2007), but there is no indication of deviation from it.

Now we consider the radially outward decrease in the magnitude of the lag, which has yet to be fully explained. This decrease begins just within the optical radius (near 7' or 9 kpc) and continues until just past the start of the warp. Thus, it may be concluded that this is an effect independent of the warp. By comparing the span of this radial variation in lag with the distribution of star formation as traced in Hα and 24 μm emission (Figure 1), the onset of the decreasing lag magnitude coincides with a sharp decline in star formation, again going against the tentative trend for extraplanar DIG rotation found by Heald et al. (2007). However, it is premature to assume that star formation is the only process affecting the magnitude of the lag and other factors may be the cause. Simple geometric arguments predict a shallowing of the lag at large radii. The most basic interpretation of the geometry, without any consideration of gas dynamical effects, considers the relationship between the rotational velocity and the gravitational potential: gas rotating at high z is farther from the gravitational center than the gas at the same radius in the midplane, resulting in a slower rotational velocity. This effect is most noticeable at smaller radii, where a certain height above the disk is a larger fraction of the radial distance from the center, resulting in a larger relative change in the gravitational potential, and thus a steeper lag.

Additionally, a radial pressure gradient could affect the lag and its radial variation. According to Benjamin (2002), a radial pressure gradient directed inward may cause a steepening of the lag. In contrast, an outward pressure gradient may result in a shallowing of the lag or even an increase in rotation velocity. Such gradients could also vary with radius. Given the uncertainty in the pressure gradients in question, as well as the potentially substantial impacts even small pressure gradients may have, it is difficult to say to what degree these will alter the magnitude of any given lag.

Observationally, a radially varying lag is also seen in NGC 891 (Oosterloo et al. 2007), where the gradient in the inner regions is −43 km s⁻¹ kpc⁻¹ and decreases to −14 km s⁻¹ kpc⁻¹ in outer radii. This shallowing of 2.5 km s⁻¹ kpc⁻² is comparable to what the models indicate for NGC 4244, although the uncertainties (as judged by eye) make this a difficult comparison. Furthermore, the shallowing of the lag in NGC 891 was seen closer to the center in both Hα (Kamphuis et al. 2007) and H i. Thus, the origins of the shallowing of the lags in NGC 4244 and NGC 891 may be different.

Also displaying radial variation in the lag is the DIG in one quadrant of NGC 4302 (Heald et al. 2007), but the trend found in that galaxy was opposite from what we detect in the H i in NGC 4244. NGC 4302 displayed a substantial steepening of the lag, going from approximately −23 km s⁻¹ kpc⁻¹ to almost −60 km s⁻¹ kpc⁻¹ at radii outward of 4.25 kpc.

Ultimately, a larger sample of galaxies where this phenomenon is observed and modeled will greatly aid in determining its cause and nature.

First, we consider the shell-like feature in the approaching half. This feature is above and slightly offset radially from a region of star formation (Figure 11), which extends between 35 and 46 along the major axis in the Hα and 24 μm data. The feature's proximity to a region of star formation as well as the lack of nearby external features indicate an internal origin. To estimate the energy required to produce it, the number density must first be found. For the calculations below, we assume this feature is indeed a shell created by multiple supernovae in the disk. Following the approach taken by Rand & van der Hulst (1993), we measure the peak flux in the limbs of the presumed shell to acquire a column density. To obtain the number density, we then consider the path length along the line of sight (l) to be

$$\begin{equation} 2\sqrt{d^{2}+2dr}, \end{equation} \tag{ 2 }$$

where d is the shell wall thickness and r is the inner radius of the shell. Since the shell walls are not well resolved, we set d equal to the beam resolution parallel to the disk of about 17''. From this, a value of 1280 pc is obtained for l, which yields a number density of 0.15 cm⁻³.

Now using Figure 12, an estimate for the expansion velocity (V_exp) is 25 km s⁻¹. Together with the estimate for the shell's radius in parsecs (R_shell = 750 pc) and assuming a constant expansion, its age is approximately 2.9 × 10⁷ years.

We now calculate the energy required to produce such a shell from Chevalier (1974) where n₀ is the H i number density in cm⁻³, R_shell is the radius of the shell in parsecs, and V_exp is in km s⁻¹ as described above:

$$\begin{equation} E_{\rm E}=5.3\times 10^{43}n_{\rm 0}^{1.12}R_{\rm shell}^{3.12}V_{\rm exp}^{1.4} \hbox{erg}. \end{equation} \tag{ 3 }$$

Given these estimates, we find that a single burst of energy of roughly 5.4 × 10⁵³ erg or the equivalent of 540 supernovae would be required to produce this shell (assuming an energy of 10⁵¹ ergs per supernova). If the shell instead formed from a continuous supply of energy as described in McCray & Kafatos (1987),

$$\begin{equation} E_{\rm C}=1.16\times 10^{41}n_{\rm 0}R_{\rm shell}^{5}t_{\rm 7}^{-3} \hbox{ erg}, \end{equation} \tag{ 4 }$$

where t₇ is the age in units of 10⁷ years, this would only require 1.7 × 10⁵³ ergs or 170 supernovae.

Now considering the Hα luminosity of this region we again use the total Hα luminosity for NGC 4244 found in Hoopes et al. (1999). The fraction of the total Hα emission originating from this region is approximately 4%. Thus, its luminosity, uncorrected for extinction, is found to be 4.0 × 10³⁸ erg s⁻¹.

To correct for extinction, we examine the 24 μm data. Summing the flux over the same region as for the Hα data and once again assuming a monochromatic flux, the 24 μm luminosity is approximately 4.8 × 10³⁹ erg s⁻¹.

Using the method described in Kennicutt et al. (2007), we determine the star formation rate, this time with a = 0.038 and replacing L_TIR with the 24 μm luminosity to account for localized regions.

The star formation rate for this region is found to be 0.0046 M_☉ yr⁻¹, which would account for 3.8% of the ongoing star formation in NGC 4244. However, Kennicutt et al. 2007 caution that this conversion is intended for entire galaxies rather than localized regions, which limits its physical meaning.

Using this star formation rate (assuming it is constant) as well as the estimated age for the shell, a calculation akin to that presented in McKee & Williams (1997) yields an estimated 650–700 supernovae to have gone off in the region, rendering supernovae a very plausible source.

Attending now to the curved feature in the receding half shown in Figure 13, we examine its potential connection to star formation. It is likely that there is some correlation between this feature and the underlying star forming region seen in the both Hα and 24 μm data. The best resolved emission is seen in the Hα data, which shows a curious extension of emission from possibly blown-out gas (Figure 13, bottom panel, highlighted in boxes) which to some extent fills in the gap defined by the H i contours, indicating that the H i gas may have been forced outward in that direction. This structure does not appear in the 24 μm data (not shown). When examining individual channel maps, it can be noted most prominently at 294 km s⁻¹ and 298 km s⁻¹, a gap in the H i forms, once again suggesting a channeling of the gas away from the midplane and into higher z in this area (Figure 14). This hole in the H i emission may be analogous to those seen in other galaxies such as NGC 4559 (Barbieri et al. 2005). However, it is possible that the features are coincidental.