FINDING DWARF GALAXIES FROM THEIR TIDAL IMPRINTS

The nearby Whirlpool galaxy (M51) and its companion (NGC 5195) provide a useful test case for our method. Figure 1 shows a map of the atomic hydrogen (H i) surface density distribution of M51 from THINGS (The H i Nearby Galaxy Survey; Walter et al. 2008). The companion of M51 is marked in the figure, and lies at the short arm, as can be seen from optical images (Bastian et al. 2005). At first sight, one notes the striking spiral structure in the inner regions, and the broad H i arm that extends far out and curves toward the east. Our focus here will be on what we can learn from analysis of the extended H i arm in the outskirts of M51.

Zoom In Zoom Out Reset image size

Much work has been done on modeling the dynamical interaction between M51 and its companion. Most of this work has been focused on the response of the stellar disk (Toomre 1981; Hernquist 1990b) excited by M51's companion, while recent hydrodynamical simulations (Dobbs et al. 2010) have not taken important physics into account, such as dynamical friction that are relevant for a massive perturber, like NGC 5195 (Smith et al. 1990). Moreover, to date, all such studies have attempted to only do forward modeling, i.e., to recover observed features in M51 or its companion. Our approach here is distinct from these prior studies—we show that analysis of the disturbances in the gas disk of M51 with respect to simulations allows one to address the inverse problem and quantitatively characterize M51's companion without requiring any knowledge of its optical light. We show below that this is sufficient to recover the position of the satellite and its mass. Moreover, our high-resolution global hydrodynamical simulations of M51, which model the gas, stars, and dark matter, include a sub-resolution model of star formation, energy injection from supernovae, and treat dynamical friction, are the most sophisticated simulations of M51 to date.

Figure 2 displays the projected gas surface density images from our best-fit simulation of M51 as a function of time. The best-fit time, i.e., the time at which the Fourier amplitudes of the simulation best fit those of the data, occurs at t ∼ 0.3 Gyr. As in CB09, the best-fit simulation is found by computing the complex Fourier transform,

$$\begin{equation} a_{m}(r,t)=\frac{1}{2\pi }\int _{0}^{2\pi } \Sigma (r,\phi,t)e^{-im\phi }d\phi, \; \end{equation} \tag{ 1 }$$

of the data and of the simulations as a function of time, where Σ(r, ϕ, t) is the projected gas surface density at time t. We calculate the residuals of the m = 0–4 modes of the data and the simulations for the modulus of a_m(r, t). The residuals for a given simulation are calculated as follows: S_{1 − 4} = ∑_r|[a'_{1, D}(r) − a'_{1, S}(r, t)]² + [a'_{2, D} (r) − a_{2, S}'(r, t)]² + [a'_{3, D}(r) − a_{3, S}'(r, t)]² + [a'_{4, D}(r) − a_{4, S}'(r, t)]²|. Here, a'_{m, D} and a'_m S denote the modulus of the Fourier transform for the data (D) and the simulation (S), respectively, normalized to the axisymmetric mode. The quantity S₁ = ∑_r|[a'_{1, D}(r) − a'_{1, S}(r, t)]²| is the residual of the m = 1 mode only. The best-fit time snapshot is that which minimizes S₁ and S_{1 − 4} for a given simulation. The entire simulation set is searched accordingly. The quantity Σ(r, ϕ) for the data is determined with respect to the photometric center (Walter et al. 2008). For the simulations, Σ(r, ϕ) is determined with respect to the bulge center of mass if a bulge is present, and with respect to the baryonic center of mass if a bulge is not included.

Zoom In Zoom Out Reset image size

It is clear from inspection of Figures 1 and 2 that it is at the best-fit time that the simulation produces the optimal visual match to the data image of M51 as well. The best-fit simulation has a perturber mass ratio of 1: 3 and a pericentric approach distance of 15 kpc. The mass estimate agrees closely with observation estimates for the mass of M51's companion, and with recent hydrodynamical studies of M51 (Dobbs et al. 2010; Salo & Laurikainen 2000; Smith et al. 1990). Furthermore, at the best-fit time, the perturber in the simulation is at a distance that is 1.5 times as large as the extent of the stellar spiral arms of M51, which is within a factor of ∼1.2 of the observed ratio. This ratio of distances (the location of the perturber relative to the extent of the stellar spiral arms), which is similar to an aspect ratio, is not affected by projection effects. Projection effects introduce an ambiguity into a purely radial distance, and the true deprojected distance to M51's companion is not known. As shown in Figure 1, M51's companion lies at the short arm, and this is also the case at t = 0.3 Gyr, at which time we achieve the best fit to the Fourier amplitudes. We discuss later that while the range of azimuth as determined from the phase of the m = 1 mode varies by ∼20 deg due to the variation of initial conditions and orbits, the azimuth in all cases is in the right quadrant. It is important to note that our method will work best for galaxies that are seen close to face on, which is the case for both M51 and NGC 1512. This is of course not true in general, and is a limitation of our method at present.

The Fourier amplitudes constrain both the mass and pericenter distance of the satellite. This may seem surprising as the tidal force is proportional to M/R³ (i.e., it is the first term in the approximation for the tidal force when the pericenter distance is large relative to the galaxy size) and one might expect a degeneracy between M and R. However, as first shown by CB09, one can break this degeneracy when one is not in the impulse approximation. The impulse approximation assumes that the primary galaxy does not respond along the course of the orbit of the perturber. We showed earlier (Figure 3 in CB09) that in reality, the nonlinear response of the gas in the primary galaxy along the course of the orbit of the perturber allows one to break this degeneracy, and determine M_s and R_peri independently, from the Fourier amplitudes. We also showed that the resultant Fourier amplitudes are not significantly affected by the initial conditions of the simulated primary galaxy, such as the gas fraction (provided the gas fraction is varied to within the range of typical spirals), equation of state, and the inclination of the orbit (CB09; Chang & Chakrabarti 2011).

Zoom In Zoom Out Reset image size

Shown in Figures 3(a) and (b) are the local Fourier amplitudes of the data and best-fit simulation (at the best-fit time), respectively. The Fourier amplitudes are shown here normalized to the axisymmetric (m = 0) mode. The best-fit simulation (and best-fit time) is determined by minimizing the difference between the variances of the Fourier modes of the data and the simulations. Although there is an overall level of agreement between the amplitudes of the simulations and the data, we are not seeking here to match the detailed features due to the possibility of missing diffuse emission in the interferometric H i maps. In other words, while there is general agreement for the quantities S₁ and S_{1 − 4} (which are summed over radius), there are differences at particular radii between the simulation and data Fourier amplitudes. What we find is that a perfect fit to the Fourier amplitudes is not necessary (and indeed we do not attempt to do this) to infer the mass and pericenter distance to within a factor of two, which given the other dependencies in the problem and the lack of single-dish data is of sufficient accuracy for our purposes here.

The magnitude of the m = 1 and m = 2 modes in the outskirts (r > 15 kpc) in both the simulations and the data is ∼0.5. We focus our analysis here on the outskirts (r > 15 kpc) because these regions are beyond the outer Lindblad resonance where we can more cleanly separate the effects of the perturber from the self-excited response. We estimate the pattern speed of the two-armed spiral of the primary galaxy when it is evolving in isolation using the Tremaine & Weinberg (1984) method. The outer Lindblad resonance is then given by graphically finding the radial region where Ω_p = Ω + κ/2, where κ is the epicyclic frequency, which gives r_OLR ∼ 16 kpc. This is an approximate measure, as there are multiple spiral patterns that are time dependent. We have checked that for r > 20 kpc, the Fourier amplitudes of the primary galaxy evolving in isolation are much lower (a few percent at most). Thus, we are confident that in the outskirts, we are looking at the tidal response produced by the companion. Furthermore, the outskirts are H i dominated and are thus less subject to the effects of feedback from supernovae and star formation which complicate the interstellar medium (ISM) structure (and the modeling thereof) in the inner regions of galaxies.

We note here that we are not seeking to identify a given model as particularly unique so as to allow the inference of orbital parameters. As mentioned in related papers (CB09; Chakrabarti & Blitz 2011; Chang & Chakrabarti 2011), the degeneracy due to inclination, for a given M_s, R_peri pair, does not allow us to determine the orbital history. This, coupled with the (weak) dependence on initial conditions of the simulated primary galaxy, means that we cannot clearly determine which orbital inclination, for a given mass and pericentric approach, would most uniquely fit the data. What we aim to do here is to convey the general trends in how well different simulations fit the data via their location on the variance versus variance plot.

Figure 4 displays the simulation parameter space surveyed, and also which simulations best fit the data for M51. Figure 4 is a variance versus variance plot of Fourier mode m = 1 relative to the data (S₁), plotted versus the variance of Fourier modes m = 1–4 added in quadrature (S_{1 − 4}). On such a variance versus variance plot (which shows the difference of the Fourier modes in a given simulation relative to the data at their respective best-fit times), the best fits will lie close to the origin. A given simulation is labeled by the mass ratio of the satellite and pericentric approach distance, i.e., 3R15 is M51 interacting with a 1:3 mass ratio perturber with a pericentric approach distance of 15 kpc. Table 1 lists the symbols used in Figure 4 and the simulations they refer to. The "e" inclination refers to θ₁ = 30, where θ₁ is the angle of the primary galaxy relative to the orbital angular momentum plane of the system (we adopt the same nomenclature here for the orbits as in Cox et al. 2006 and Barnes 1988). It is immediately clear from Figure 4 that low-mass perturbers (with mass ratios less than ∼1:10) cannot fit the data. The quantity S_{1 − 4} scales approximately as the square root of the mass ratio of the satellite for a fixed pericentric approach distance for low-mass companions, a result also found by Chang & Chakrabarti (2011; for a similarly defined sum of the Fourier amplitudes). The best fit to the Fourier amplitudes of M51 (where the best fit is derived by minimizing the difference in local Fourier amplitudes of the simulations and data between r = 25 and 40 kpc for all simulations as a function of time) is accomplished with a 1:3 mass ratio satellite with a pericentric approach distance of 15 kpc.

Zoom In Zoom Out Reset image size

We find that a pericentric approach distance of 15 kpc for a 1:3 mass ratio satellite yields a better fit than a 3R7 simulation (a 1:3 satellite with pericenter distance of 7 kpc). A 3R15 (black cross, and black filled circle, where the black filled circle corresponds to an isothermal equation of state and the black cross to a multiphase ISM) yields a slightly better fit than a 4R15 (red asterisk) and a significantly better fit than a 2R15 (green cross). 3R15 is significantly favored over 3R7 (points in blue) or 3R25 (orange asterisk). A 3R7 yields a better fit than a 4R7 (red cross) or a 5R7 (yellow cross). The points in blue all correspond to the "Variants of 3R7," where we have held the mass of the companion and pericenter distance constant but varied other parameters that may affect the Fourier amplitudes including orbital inclination, equation of state, bulge mass fraction, and orbital velocity. We do this to understand whether the mass and pericenter distance are truly the dominant parameters that drive the location of a given simulation on the variance versus variance plot. While the parameter space study that we show here is not exhaustive, we try to delineate the important dependencies of the problem and those that are not as significant. We see that varying orbital inclination (blue asterisk showing an inclined orbit), or changing the equation of state (blue filled circle), or modest changes in the orbital velocity (blue triangle, 3R7V+; blue square, 3R7V−), or including a bulge in the simulation (blue diamond) does not significantly affect the Fourier amplitudes. We have shown here cases where we increased (decreased) the relative orbital velocity of the companion by 30% relative to the standard case, finding that a decrease in the relative orbital velocity (3R7V−) with respect to our standard case is slightly favored, but not significantly so. We have earlier demonstrated in CB09 that varying gas fraction from ∼0.2–0.5 does not significantly affect the Fourier amplitudes. In Chang & Chakrabarti (2011), we presented an extensive study of the dependence of the Fourier amplitudes on orbital inclination, so we only show a couple of cases here to make the point that the orbital inclination of the perturber does not greatly affect the Fourier amplitudes. Inspection of the cases where the perturber mass and pericentric distance are held constant, but initial conditions, orbital inclination, equation of state, and orbital velocity are varied, shows that these other parameters affect the fit to the data only slightly. This weak dependence on all parameters aside from the mass of the perturber and pericentric approach distance allows us to solve the inverse problem. However, the dependence on other parameters does make it difficult to determine the pericentric approach distance and satellite mass to better than a factor of two.

In summary, the best-fit mass ratio for the satellite for M51 is a 1:3 satellite with a pericentric approach distance of 15 kpc. The uncertainly in the mass determination is a factor of two, i.e., a 1:5 yields a significantly poorer fit, and a 1:2 is also a poorer fit relative to the 1:3; the pericenteric approach is similarly unconstrained at the factor of two level. For massive companions (≳ 1: 10) that approach close to the galactic disk (i.e., close to the regions where the galactic disk has an exponential profile), the problem is not scale-free, and the Fourier amplitudes depend in a more complex way on the mass ratio than the square root dependence (for fixed R_peri) found by Chang & Chakrabarti (2011). We do not attempt to obtain scaling relations for the Fourier amplitudes for massive companions in this paper. It is nonetheless clear that we can constrain the satellite mass and distance of closest approach using the numerically calculated Fourier amplitudes to a factor of two, providing a quantitative (albeit approximate) means to hunt for dark-matter-dominated galaxies even when they may be optically dim.

To demonstrate the dependence of the Fourier amplitudes on the inclusion of a bulge, we explicitly show a case (depicted in the blue diamond) where we have included a bulge (with a bulge mass fraction of 0.012 and a scale length for the bulge of 0.2 relative to the disk scale length) and calculated the Fourier amplitudes from the bulge center of mass. As is clear, this simulation is not markedly different from the other simulations in S₁ or S_{1 − 4}. In principle, an offset in the center of mass (of the component with respect to which the center of mass is calculated) can introduce or suppress power in the m = 1 mode. However, when it is expressed as the global quantity S₁, we find little change with respect to center of mass variations of the gas and bulge components. Moreover, while the local Fourier amplitude of the m = 1 mode can be affected by different choices of the center of mass, the sum of the modes remains nearly constant, as the sum is equivalent to distributing power between the various modes. Blitz (1994) gives a detailed review of various choices of the center of mass of the Milky Way and concludes that all the mass concentrations near the Galactic center appear to have non-zero radial velocities. Thus, it is difficult to establish which is the correct choice. Fortunately, this consideration, as well as the inclusion of a bulge in the model, does not seriously affect the placement of a simulation in the S₁–S_{1 − 4} plane.

Simulations can also be located on the variance versus variance plots as a function of time. We consider the best-fit simulation of M51 as an example. At early times, the location of the simulated M51 is far from the origin (thus indicating a poor fit to the data); it hovers close to the origin round t ∼ 0.3 Gyr, which is the best-fit time, and at late times is also far from the origin. Thus, the time evolution of the simulated M51 on the variance versus variance plot allows us to constrain the time of encounter, which is roughly a dynamical time prior to the best-fit time. In gaseous disks, disturbances damp out on the scale of a dynamical time. Therefore, gaseous disks only possess a short-term memory of encounters, which we can utilize to infer the time of encounter. The inference of the time of encounter is a critical feature of our calculation—inferring the time of encounter is what allows us to determine the current radial and azimuthal location of the satellite.

The (relative) radial location of the perturber does not quite give us enough information to find dark (or nearly dark) galaxies. The azimuthal location is also needed. The H i image of the primary galaxy in fact provides this information. The image is composed not only of the Fourier amplitudes of the modes, but also the phase of the modes. In CB11, we employed the relative offset in the phase of the modes between the data and the simulations at the best-fit time (the time that best fit the Fourier amplitudes) to determine the perturber azimuth. This procedure is similar to visual matching of (dominant) features between the simulations and the data. We determine the azimuth of M51's companion here from the relative offset between the phase of the m = 1 mode in the simulations and the data at the best-fit time. The phase of the modes has been computed by taking the Fourier transform of the projected gas surface density:

$$\begin{equation} \phi (r,m) = \arctan \frac{[-{\rm Imag\,FFT} \Sigma (r,\phi)]}{[-{\rm Re\,FFT} \Sigma (r,\phi)]} \;. \end{equation} \tag{ 2 }$$

The phase of the modes contains information on the shape of the spiral planform. Tightly wrapped spirals produced by the self-excited spiral structure inside of the inner and outer Lindblad resonances will have sharp gradient in the phase, while open spirals produced by tidal interactions will have a flatter profile (Shu 1984). We show below that the phase of the m = 1 mode for this known tidally interacting system is indeed flat in the outskirts, as is that of the simulation at the best-fit time. It is also worth noting that the time variation of the phase independently provides a handle on the time of encounter, just as the Fourier amplitudes do. The phase in the outskirts is not always flat as a function of radius. At early times, it resembles that of an isolated galaxy (prior to the encounter), it is flat about a dynamical time after pericentric approach, and the spirals wind up at late times, producing a gradient in the phase. Thus, if the phase is flat at the same time one achieves the best fit to the Fourier amplitudes, this is compelling evidence of a tidal encounter.

We constrain the satellite's azimuth by translating from the simulation frame to the observational frame by using the relative offset in the phase of the m = 1 mode between the simulations and data. This procedure is similar to matching of dominant features (in this case the m = 1 mode) between the simulation and observed H i image. Figures 5(a) and (b) display the phase images of the m = 1 mode for the simulation of M51 and the H i data, respectively. Both of these images show considerable structure in the outskirts (r > 15 kpc) which we show in a line plot rendering in Figure 6. The angle of the perturber in the observational frame as determined by the relative offset in the phase is given by ϕ_perturber = ϕ^sim_sat − ϕ^sim₁ + ϕ₁^data, where the ϕ'₁s are the phases of the m = 1 mode in the simulations and data. This calculation is necessary as we cannot assume that the simulations and observations are exactly aligned. The quantity [ − ϕ^sim₁ + ϕ₁^data] tells us how we have to rotate the simulation coordinate system. We have denoted the azimuth of the satellite in simulation coordinates as ϕ^sim_sat, which is the angle of the satellite in the M51 best-fitting simulations, and is equal to 203 deg (on average) in center-of-mass coordinates. The relative offset between the phase of the modes, which we calculate from Figure 6(a), is given by the median of the quantity [ − ϕ^sim₁ + ϕ₁^data], evaluated from r = 15 kpc to r = 40 kpc. This calculation yields an angle for the perturber of 81 deg, which is consistent with NGC 5195 lying close to the tip of the short arm of M51 (or roughly 90 deg from the x-axis). This is a significant result—it derives simply from TA, i.e., from the calculation of the relative offset of the m = 1 mode in the simulation and the data. The range of azimuth varies by at most ∼20 deg due to the variation of initial conditions and orbits for the range of acceptable fits discussed earlier. We depict in Figure 6(b) the absolute relative difference of the inferred azimuth (ϕ_inf) from the real azimuth of M51's satellite (ϕ_real) as a function of S₁ for several simulations. There is a general trend for this relative difference to be smaller as S₁ decreases, although there is some spread in the uncertainty at a given S₁ due to the fact that the azimuth determination depends somewhat on the other modes as well. Nonetheless, the lowest S₁ values and associated azimuth values represents our best determination, which is fairly accurate for the best-fit cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The NGC 1512/1510 interacting system provides another opportunity to test our method. Figure 7 displays the Australia Telescope Compact Array (ATCA) image from the THINGS-SOUTH program of NGC 1512 by E. De Blok et al. (2011, in preparation). NGC 1512 is a barred spiral galaxy, and its companion is a blue compact dwarf galaxy, NGC 1510, which is marked by a cross in Figure 7. We initialize our simulation of NGC 1512 with the parameters inferred by Koribalski & Lopez-Sanchez 2009), and these parameters, along with the orbits are given in Section 2.

Zoom In Zoom Out Reset image size

Figure 8 displays the projected gas surface density images of the best-fit simulation for the interaction of the NGC 1512/1510 system. The best-fit time is at t = 0.18 Gyr, when the Fourier amplitudes of the simulations most closely match those of the data. The best fit is achieved for a 1:50 perturber with a pericentric approach distance of 8 kpc. The separation between the perturber and NGC 1512 at the best-fit time in the simulation is a factor of two larger than the extent of the stellar spiral arms, which is reasonably close (to within a factor of ∼2) to the observational determination. The mass inferred for NGC 1510 is close to that estimated by Koribalski & Lopez-Sanchez 2009) from the observed H i flux (to within a factor of ∼1.5), although no dynamical mass estimates are available for NGC 1510.

Zoom In Zoom Out Reset image size

Individual H i clouds in the NGC 1512/1510 system have been found out to ∼80 kpc by KS09 that they denote tidal dwarf galaxies (TDGs). They infer that the age of the stellar populations in the TDGs could be as young as 150 Myr and possibly as old as 300 Myr. These ages are coincident with the difference in time between pericentric approach and the best-fit time. As before, we use the phase of the m = 1 mode in the simulation and the data, to infer the azimuthal location of NGC 1510 (denoted by crosses in Figure 8 in the simulation frame). We find that the azimuth of NGC 1510 is in the southeast corner, in agreement with the observed image.

Figure 9 depicts the global Fourier amplitudes of the best-fit simulation of NGC 1512 as a function of time, with the data shown as a vertical stripe (in green). NGC 1512 is quite different from M51 or the Milky Way (as it is lower in mass by a factor of 10), and its companion is smaller in mass ratio than M51's. We choose to employ the global Fourier amplitudes for this data set, i.e., we define

$$\begin{equation} C_{\rm m}= \int \nolimits _D \Sigma (r,\phi,t)e^{-im\phi }r dr d\phi \;, \end{equation} \tag{ 3 }$$

where D is the integration bounds over ϕ from 0 to 2π, and over the radial range from r_min (which should be greater than the outer Lindblad resonance to select regions that are dominated by the perturber rather than the self-excited response) to R_max. We use the global rather than the local Fourier amplitudes because the NGC 1512 data set is constructed from a mosaic of individual interferometer pointings. The mosaicking leads to varying signal-to-noise levels across the field of view, which makes conclusions based on local Fourier amplitudes less reliable. The global Fourier amplitudes are less affected than the local Fourier amplitudes by the mosaicking as the global Fourier amplitudes are integrated over both radius and angle. We find the best fit to the global Fourier amplitudes for a 1:50 mass ratio satellite (which is consistent with the observational estimate by Koribalski & Lopez-Sanchez (2009) with R_peri = 7 kpc, which occurs at t = 0.18 Gyr. Figure 10 is the variance versus variance plot computed using the global Fourier amplitudes for NGC 1512. Although the global Fourier amplitudes do not give us as much information and therefore less ability to discriminate between 1:50 and 1:100 satellites (the two points shown in green and yellow), we can clearly discriminate between ∼1:50 and 1:10 (the red asterisk) to see that a 1:10 mass ratio satellite cannot have produced the disturbances in NGC 1512's H i disk. As before, we find that modest changes in the orbital velocity of the companion do not significantly affect the Fourier amplitudes (depicted by the black triangle and black square).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 11(a) and (b) display the phase m = 1 modes of the simulation and H i data of NGC 1512, respectively. The line plot rendering of these images is shown in Figure 12(a). The azimuth of NGC 1510 derived from the relative offset as discussed previously gives a value of 267 deg (with a maximum uncertainty of ∼20 deg). This is consistent with the observed image of NGC 1512 as presented in Koribalski & Lopez-Sanchez (2009), where NGC 1510 is shown to lie in the right-hand bottom quadrant. In Figure 12(b), we depict the relative uncertainty in the determination of NGC 1510's azimuth as a function of C₁ for several simulations, and again find the azimuth to be determined quite accurately ((ϕ_inf − ϕ_real)/ϕ_real ≲ 10 %) for low C₁ values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have utilized the phase of the m = 1 mode as a separate and independent constraint to infer the azimuth of the companion once the mass and pericenter distance are constrained. The phase of the higher order modes introduces an m-fold degeneracy, that were we to obtain a perfect fit to the data, may allow for more discrimination if the phase of the modes were additionally used, along with the Fourier amplitudes, to constrain the mass and pericenter distance. In this proof of principle paper, we choose to treat them as independent constraints. We find that it is particularly compelling that the shape of the phase of the m = 1 mode independently provides a handle on the time of encounter just as the Fourier amplitudes do. In principle, the phase and Fourier amplitudes are independent constraints. The model's validity (and indeed this method's) is reinforced when these two metrics independently yield the same information.