Dynamics and Shocks from H Emission of Nearby Galaxy Mergers
We examine the dynamical properties of interacting galaxies and the properties of shocked gas produced as a result of the interaction. We observed 22 galaxy mergers using the SparsePak IFU at Kitt Peak National Observatory (KPNO). The goal of the observations was to obtain the H velocity maps over the entire luminous parts of the galaxies including the faint tidal tails, and to find extended shocks and outflows. Our sample consists of major and minor galaxy mergers with mass ratios . We fit multiple kinematic components to the H and [N II] emission lines, develop an MCMC code to robustly estimate the error of fit parameters, and use the F-test to determine the best number of kinematic components for each fiber. We use [N II]/H and velocity dispersion of components to separate star-forming (HII) regions from shocks. We use the kinematics of the H emission from HII regions and an automated modeling method to put the first ever constraints on the encounter parameters of one of the observed systems. Besides, we estimate the fraction of shocked H emission, , and examine the spatial distribution of shocks. We find that close galaxy pairs have, on average, a higher shock fraction than wide pairs, and coalesced mergers have the highest average . In addition, galaxy pairs with more equal mass ratio tend to have a higher . Combining the dynamical models from the literature and this work, we inspect trends between and dynamical encounter parameters. Our findings are generally consistent with shocks being produced either by direct collision of the ISM or by the chain of events provoked by the tidal impulse during the first passage.
keywords:galaxies: interactions – galaxies: kinematics and dynamics – galaxies: peculiar
Galaxies have evolved in size, shape, kinematic properties, and star formation rate (SFR) over the past few billion years (Bell et al., 2004; Faber et al., 2007; Brown et al., 2007), and galaxy mergers may have played a significant role in this evolution (Toomre & Toomre, 1972; Narayanan et al., 2008; Lotz et al., 2011; Rodriguez-Gomez et al., 2016). In most galaxy mergers, the gravitational tidal force between interacting galaxies and the orbital energy that is converted into internal energy of the remnant can redistribute stars, gas, and dark matter particles in the phase space, leading to the transformation of the shape and kinematic properties of galaxies over time (Toomre, 1977). Moreover, interaction affects the star formation rate (SFR) both in the core and in the outskirts of gas rich galaxies (Barnes & Hernquist, 1996; Sanders & Mirabel, 1996; Mihos & Hernquist, 1996; de Grijs et al., 2003; Cox et al., 2008). In case of certain encounter geometries, gas clouds of two or more interacting galaxies directly collide, and due to ram pressure and tidal stripping and shock-heating, the fraction, temperature, and density of the gas in the remnant are affected (Cox et al., 2004a; Soto & Martin, 2012; Vollmer et al., 2012). Studying properties of gas in interacting galaxies allows us to evaluate the comparative significance of each process in the overall evolution of galaxies.
A case study can be the evolution of the sources of nebular emission as a merger proceeds, particularly the evolution of the fraction of shock-heated gas. Merger induced shocks in nearby tidally interacting galaxies are not only observed, but also predicted by hydrodynamical simulations of binary galaxy mergers (e.g. see Monreal-Ibero et al. 2010; Rich et al. 2011; Belfiore et al. 2016 for observations and Cox et al. 2004a for simulations). As we get closer to coalescence, we expect significant shock heating to take place as a result of the mixing and collision of random gas orbits as well as the return of gas removed previously through tidal striping (Cox et al., 2004b). Though, shocks have also been observed in early stage mergers (Soto & Martin, 2012). In the earlier stages, shocks may be produced through different modes. In one mode, each of the interacting discs experiences a strong tidal force during the first passage. The maximum of this tidal force on each galaxy depends on their separation at pericenter, as well as the mass of the companion(s). Tidal force results in a gas flow toward the center of each disc (Mihos & Hernquist, 1994; Barnes & Hernquist, 1996). The inflowing gas may collide with the gas at lower radii at a high enough velocity to induce shocks. The inflowing gas may also reach the core and trigger (or enhance) starburst or Active Galactic Nucleus (AGN) generating strong outflows (Ellison et al., 2011; Scudder et al., 2012). These outflows produce shocks as they blow into the interstellar medium (ISM). In this mode, more shocks are expected if the tidal force is stronger, so we anticipate to find higher shock emission when either the companion is more massive or the pericentric separation is smaller. In a different mode, collision of galaxies with certain geometries can result in a fast direct encounter of gaseous clouds of the ISM in the two discs, creating widespread shocks (Vollmer et al., 2012; Medling et al., 2015). For both of these modes, The chain of events would begin after the first passage, but the timeline would depend on the mode. In the first mode, a merger-driven central starburst or AGN would produce shocks a few hundreds of Myrs after the first passage (Ellison et al., 2011; Scudder et al., 2012). In the second mode, in which shock are due to the collision of gas clouds in head-on encounters (e.g. in UGC 12914/5), shock production would be immediate. Observing a sample of interacting galaxies and inspecting the relationship between encounter parameters and shock fraction would put new observational constraints on these speculations.
For example, Rich et al. (2015) used galaxy-wide H emission in a sample of 17 (U)LIRGS to show the relation of shock fraction and merger stage. They binned galaxies (mergers) into four groups based on their separation: isolated galaxies, wide pairs with separation between 10-100 kpc, close pairs with separation10 kpc, and coalesced mergers. Removing unambiguous AGN hosts from their sample, they found that going along the sequence from isolated galaxies to coalesced mergers, an increasing fraction of H is emitted from shock-heated sources, which must be a result of increasing turbulence as merger proceeds. However, the projected separation of interacting galaxies does not necessarily represent the merger stage. An interacting pair of galaxies just after the first passage are as close to each other as a pair near coalescence. Also, two galaxies at large separation appear to be a close pair from certain viewing angles. In order to properly constrain the merger stage, and other possible encounter parameters affecting shocks, we need a dynamical model of the merger system.
Dynamical modeling constrains merger stage and other encounter parameters (e.g. pericentric distance, eccentricity, etc.) utilizing simulations that best reproduce tidal tails and bridges. Simulations that are used for this purpose are usually collisionless N-body simulations, because incorporating gas physics extensively increases the computational demand (Barnes & Hibbard, 2009; Holincheck et al., 2016). We often need to measure and model the line of sight velocity of tidal features, in addition to their apparent shape, to ensure the uniqueness of the model (Barnes, 2011). As a kinematic tracer, H emission is relatively easy to measure for star-forming galaxies even in low surface brightness regions. Galaxy mergers often induce star formation both in the center and in the tidal tails (Jog & Solomon, 1992; Hattori et al., 2004; Whitmore & Schweizer, 1995; de Grijs et al., 2003; Scudder et al., 2012). If we can separate H emission originating at star-forming regions from shocked gas, H observations may be used not only to measure the velocity of baryons in tidal features but also to detect shocks and their distribution throughout a merger.
Mortazavi et al. (2018) proposed a method for separating shocks from star-forming regions, which only uses the kinematics and ratio of [N II] and H emission lines, for their data did not include [OIII] and H lines. By separating star-forming regions from shocks, they found an excellent match between the velocity of star-forming H and cold HI gas in the Mice galaxy merger. They used the automated dynamical modeling method developed in Mortazavi et al. (2016) to model the Mice using both H and HI kinematics, finding consistent results. Also, they showed that the properties of shocked H emission obtained from their method is consistent with the data from CALIFA galaxy survey, which includes [OIII] and H lines (Wild et al., 2014; Sánchez-Menguiano et al., 2016).
Galaxy-wide nebular emission in a sample of tidally interacting galaxies can be used to improve our understanding of the physical processes during the merger. In order to find the encounter parameters of galaxy mergers, the easiest kinematic tracer is the H emission line. Additionally, the H and [N II] emission lines contain valuable information about the source of ionization, especially shocks. We can measure what fraction of H flux is emitted from shock-ionized gas, and where the shocks are spatially located. We may use the reconstructed encounter parameters and the measured shocks to investigate how shocks evolve with different merger parameters. Statistical approach would help us to isolate the effect of each parameter on shocks production mechanisms.
In this work, we present H observations of a modest sample of 22 galaxy mergers using the SparsePak Integral Field Unit (IFU; Bershady et al. 2004) on the WYIN telescope at Kitt Peak National Observatory (KPNO). Observations of 21 systems are new to this work, while the data for the Mice galaxy merger was presented before in Mortazavi et al. (2018). We present the analysis of emission lines and shock detection for all of the observed systems and a new dynamical model for one of them. We examine correlation between encounter parameters and shocks. In §2 we describe the observational setup, the instrument, and the target selection. In §3 we present the analysis of emission lines including fitting multiple kinematic components and an MCMC code to measure the error of the fit parameters. In §4 we discuss how to separate emission originating at shocked gas from that arising out of star-forming regions, and we inspect the relationship between H shock fraction and encounter parameters. In §5 we present our attempts to model the dynamics of equal mass mergers in our sample, and in §6 we inspect the correlations between shock fraction and some of the encounter parameters. In §7 we discuss our results and implications about how merger encounter parameters affect merger-induced shocks. We present some notes on a couple of observed systems in the Appendix.
We observed 22 galaxy merger systems using the SparsePak IFU on the WIYN telescope at KPNO. The observations took place from March 2008 to May 2013 in five observing runs consisting of a total of 14 nights, with the dome being closed three full nights and parts of some other nights due to weather conditions. Our primary goal was to measure the kinematics of H and [N II] emission lines, so we kept observing in non-photometric conditions, and occasionally thin clouds were covering our targets during the observation. Resultingly, we did not attempt to flux calibrate our data.
2.1 The Instrument
The WIYN SparsePak IFU consists of a grid of 82 sparsely packed fibers each 4.687” in diameter, covering a Field of View (FoV) of 72”71.3”. Seven sky fibers are placed at distance from the science grid on the north and east sides (Bershady et al., 2004).
The bench spectrograph and 860 lines/mm grating blazed at 30.9 in order 2 was used, obtaining a dispersion of 0.69 Å/pixel (FWHM) in the wavelength range of 6050-7000Å. (R4500, and velocity resolution 30 km/s near the H line) Our spectral coverage is less than current and ongoing galaxy surveys such as CALIFA (Sánchez et al., 2012), MaNGA (Bundy et al., 2015) and SAMI (Croom et al., 2012). Our spectral resolution, however, is higher than CALIFA and MaNGA. In red band, CALIFA, and MaNGA surveys have spectral resolution, R, of 850 and 2000, respectively. Higher spectral resolution enables us to resolve multiple emission line components, usually appearing in the central regions of galaxies where multiple gaseous components overlap.
SparsePak is especially suitable for the purpose of finding velocity maps for
dynamical modeling of galaxy mergers, as it has a sparsely packed grid that covers a relatively large FoV. Modeling the dynamics of interacting galaxies does not require a uniform velocity coverage on the system. Rather, velocity information that reveals the direction of rotation near the center and the large scale velocity gradient across the tidal tails and bridges is enough for dynamical modeling. Previous dynamical models of interacting galaxy pairs have used HI velocity maps with lower spatial resolution than is available with SparsePak (Hibbard &
Mihos, 1995; Privon et al., 2013). In this work, in order to observe each part of the targets, we mostly pointed SparsePak only once with no dithering.
SparsePak, however, is not ideal for inspecting shocks in galaxies. For studying shocks, full coverage is more desirable as we may miss regions with extensive shocks between the sparsely placed fibers. In addition, by pointing the denser center of SparsePak toward the center of galaxies, in most systems, we introduce a bias in the galaxy-wide shock fraction toward the value near the center.
2.2 Target Selection
Primarily, we selected binary systems of interacting galaxies with strong tidal features. In addition, our sample includes three coalesced systems with strong tidal features, suggesting that they have also experienced a recent major merger. Tidal features are one of the best indicators of a recent merger event of rotation-supported disc galaxies. These features would not have survived if the interacting galaxies were dispersion-supported. Besides, these features can be reproduced by test particle simulations (Toomre & Toomre, 1972; Dubinski et al., 1999), and their shape and velocity is sensitive to the encounter parameters, so they are utilized for reconstructing encounter parameters (Barnes & Hibbard, 2009; Barnes, 2011; Mortazavi et al., 2016; Holincheck et al., 2016). However, requiring strong tidal features introduces a selection bias to our galaxy merger sample. A prograde merger, for which the initial angular momentum of the discs are aligned with the orbital angular momentum, makes stronger tidal effects compared to a retrograde (polar) merger ,where discs are initially anti-aligned with (perpendicular to) the orbital angular momentum. As a result, the prograde mergers are more likely to appear in our sample.
All but one of the observed systems were in the footprint of Sloan Digital Sky Survey (SDSS York et al. 2000)
We used three different catalogs to find the merger systems that match our requirements. Our targets are chosen from the catalog of isolated pairs of galaxies in the northern hemisphere (Karachentsev et al., 1985), the Arp atlas of peculiar galaxies (Arp, 1966), and the Galaxy Zoo morphological classification catalog (Lintott et al., 2010). In the Galaxy Zoo catalog, we visually inspected galaxies with merger probability 0.5 and selected proper candidates.
Targets are selected to have a field of view (FoV) similar to that of SparsePak science grid (72”71.3”). Most of the observed systems and their tidal tails fit in two SparsePak pointings, one for each galaxy. Some smaller systems fit in one pointing, and some larger ones require three or four pointings. All of the systems analyzed in this work are shown in Figure 1 along with the SparsePak grid pointings. Table 1 shows the sky position of the targets, observation date, redshift, number of SparsePak pointings, and exposure time on each pointing.
|UGC 12914||0.4171||23.4898||Oct 2012||0.0146||3||65, 30, 125|
|Arp 256||4.7104||-10.3693||Oct 2012||0.0272||2||40, 65|
|VV 433||9.8322||13.1064||Oct 2012||0.0353||2||60, 105|
|UGC 480||11.6472||36.3286||Oct 2012||0.0374||2||97, 90|
|UGC 1063||22.2881||11.1360||Oct 2012||0.0193||2||65, 65|
|Arp 273||35.3778||39.3660||Oct 2012||0.0251||3||37, 50, 32|
|NGC 1207||47.0034||38.3769||Oct 2012||0.0160||2||50, 25|
|NGC 2623||129.6001||25.7545||Mar 2008||0.0185||3||185, 86, 90|
|Arp 283||139.3624||41.9970||Oct 2012||0.0060||3||55, 65, 95|
|Arp 181||157.1193||79.8182||May 2013||0.0326||1||30, 82|
|NGC 3509||166.0981||4.8286||Mar 2008||0.0257||1||110|
|Arp 87||175.1850||22.4379||May 2013||0.0237||2||25, 58|
|NGC 3921||177.7786||55.0788||Mar 2008||0.0197||3||120, 120, 120|
|UGC 07593||187.0612||44.4532||Apr 2012||0.0230||1||95|
|NGC 4676||191.5443||30.7271||Mar 2008||0.0220||4||150, 120, 120, 120|
|Arp 238||198.8870||62.1269||May 2013||0.0308||2||25, 35|
|NGC 5257/8||204.9805||0.8354||Apr 2012||0.0227||2||85, 90|
|NGC 5278/9||205.4237||55.6722||Apr 2012||0.0252||2||65, 65|
|Arp 84||209.6492||37.4391||May 2013||0.0116||4||105, 30, 35, 30|
|UGC 11695||318.0418||-1.4857||Oct 2012||0.0323||2||65, 125|
|UGC 12589||351.2615||0.0096||Oct 2012||0.0338||2||45, 125|
|Arp 284||354.0750||2.1557||Oct 2012||0.0093||4||39, 35, 55, 121|
2.3 Observation Setup
In order to determine the exposure time required at each pointing, we estimated the continuum flux density near the H ( Å) line using the SDSS r-band images. The surface brightness in the faint tidal features of our systems ranged from 8.4 to 3.8 erg s. We planned to have exposure time required to achieve S/N5 in the continuum in fibers placed on the tails. However, we lost part of our observing time due to bad weather conditions, so in order to observe more systems, we reduced the exposure time when we could see strong H lines in outskirt fibers in the first couple of exposures. Table 1 shows the number of pointings and the exposure time applied on each pointing.
For accurate wavelength calibration, CuAr and ThAr calibration lamps were used either before or after the science exposures on most pointings. We made three or more science exposures at each pointing to correct for the night sky lines and the cosmic rays. We took dome and twilight flats and zero exposures at the beginning/end of the night. Sky spectra was obtained simultaneously with either the seven SparsePak sky fibers or science fibers placed on blank sky. If the outskirts of the galaxies were far so that they covered the sky fibers, we oriented the IFU to put some of the science fibers on blank sky.
3 Emission Line Analysis
We did not make any attempt to fit the stellar model to the continuum because the tidal tails are usually too faint. The goal of theses observations were to obtain the kinematics of H and [N II] emission lines throughout the systems, including the faint tidal tails. Thus, we only fit emission line models to the H-[N II][, ] triplet. The absence of stellar model affects the H line flux measurements because we do not take the underlying H absorption into account. Our measurement of H flux is underestimated, so our [N II]/H is overestimated. We will discuss the effect of this in §3.3.
The relatively high spectral resolution of our data allows us to explore fitting more than one emission line component. This step is essential for separating the shocks using the velocity dispersion of emission lines. Sometimes, two gaseous components with different line of sight velocities lay on the same fiber, and their emission line profiles blend. One-component fit may measure a velocity that is offset from the velocity of both components. It may also measure a velocity dispersion that is broader than the velocity dispersion of individual blending components. In the following sections we describe the method we used for fitting one and two velocity components and estimating the errors of the fit parameters. We also address the issue of the underlying H absorption.
3.1 One-Component Fit
We first fit a triple Gaussian (three Gaussian functions with the same free velocity and velocity dispersion, , and free normalization factor, H flux and [N II] fluxes, with centers separated by the wavelength difference between [N II]6549, H6563, and [N II]6585) and a straight line with a free slope representing the background. We use the theoretical value of 2.95 for [N II]/[N II] ratio (Acker et al., 1989). Along with the slope and the offset of the background we fit a total of 6 parameters. The fit range is from 6518Å to 6608Å in the target rest frame. For velocity dispersion measurement, we subtract the intrinsic FWHM of our SparsePak observations (1.5 Å) from the fit FWHM of the lines in quadrature.
We use a Markov Chain Monte-Carlo (MCMC) method to estimate the uncertainty of the fit parameters. The initial MCMC step is determined with the least square minimization. We use the emcee python package, varying the the 6-dimensional parameter with the following prior constraints.
The velocity here is with respect to the reported redshift for each system. The log of the likelihood function is
where and are the spectrum data points, and is the error of the flux density measurement. After we reach the equilibrium distribution, we use the 16 and 84 percentile values of a 2000-step MCMC sample for the low and high error limits.
3.2 Two-Component Fit
Visual inspection of our one-component fits makes it clear that in some of the fibers the emission lines have a more complicated profile than a simple Gaussian. Sometimes the three emission lines in the H-[N II] triplet show a consistent deviation from a single Gaussian profile, indicating that the emitting gas projected into the fiber has more than one kinematic components. Figure (a)a shows an example of a fiber for which one-component triple Gaussian cannot well-describe the shape of the emission line triplet.
In order to model two-component emission lines, we fit two triple-Gaussian functions and a background line to the H-[N II] triplet. This is, basically, an extra triple-Gaussian added to the one-component model. Again, we perform MCMC to estimate the error of the fit parameters. We only attempt to fit the second component when we find H flux S/N3 in the one-component fit. The same prior constraints of Equation 1 are used for both components. However, we add the following constraints to make sure that the two components have either different redshift or different FWHM. Without these constraints, two components with the same kinematic properties can share the flux arbitrary between each other, reproducing a fit that is similar to a single component fit and providing no new information about the line profile. In order to prevent MCMC walkers from wasting time in these answers, we exclude them in the prior:
The likelihood function is the same as Equation 2. Fitting two components involves 10 parameters, so more MCMC time steps are required to achieve an equilibrium distribution. Similar to the one-component fit, after reaching equilibrium, we used the 16 and 84 percentile values of a 2000-step MCMC sample for the low and high error limits.
We use an F-test (Lomax, 2007) to determine whether a two-component model is preferred over the one-component one. In general, increasing the number of parameters improves the statistics because with more parameters the model can reproduce more subtle features in data, even if they are statistically insignificant. The F-test utilizes both the increase in the number of free parameters and the improvement in the to determine which model is preferable. As an example, for the fit in Figure 2 the F-test prefers the two-component fit, because the second component significantly improves the . We perform the F-test only when the H flux of the second component has a S/N larger than 3, otherwise the one-component fit is preferred without performing the F-test.
3.3 Underlying H absorption
By taking a straight line as the background, we are ignoring the underlaying stellar H absorption. Therefore, our measurement of H flux is higher than the actual value, and our [N II]/H is a higher limit. In Mortazavi et al. (2018) we discussed this effect for the SparsePak spectra of the Mice galaxies (NGC 4676). For the Mice galaxies CALIFA data with a stellar continuum model is available (Sánchez-Menguiano et al., 2016). Using CALIFA line ratios, we showed that the underlying H absorption only decrease the log([N II]/H) by less than 0.2 dex (See Figure 5a of Mortazavi et al. 2018). We showed that even with this overestimation we find shocks in the same regions of the Mice system as they were found in the CALIFA data (Wild et al., 2014).
As it will be discussed in §4.1, we use log([N II]/H) for shock detection, so we want to minimize the effect of the underlying H absorption on our measurement of log([N II]/H). In order to do so, for the purpose of our shock analysis, we remove fibers with H EW Å
Removing the fibers with low H EW takes out a significant number of detected components with S/N (); however, most of the H flux in these galaxies come from the luminous fibers with high H EW, so the flux in the removed fibers are not significant, and they do not much affect the measurements of shocked H fraction. In order to demonstrate this, in Figure 3 we show the both unweighted and H flux weighted histograms of H EW for all of the 22 observed systems. In this Figure, the cut of EW=7Å is shown with the vertical red dashed line. Even though a significant fraction of unweighted histogram (blue) is to the left of the EW cut, the weighted histogram (green) clearly displays that most of the flux belongs to the fibers with EW7Å. Note that our data is not flux calibrated, so for weighting this histogram we took the total H flux in electron counts and weighted them by the exposure time. We realize that this is not a proper flux calibration method because of CCD variations and non-photometric observing conditions. We only use method for qualitative estimates similar to Figure 3.
4 Source of Ionization
H is emitted from ionized hydrogen gas, and ionization can be the result of various processes. The most common source of ionization is star formation. In star-forming regions gas is photo-ionized by the UV light from young O and B stars with energies 13.6 eV. The photons from these stars are typically not energetic enough to ionize metals (e.g. N and O). Moreover, If the star formation is not happening in a turbulent region of the ISM, the velocity dispersion of this photo-ionized gas is only a few tens of kilometers per second.
The H emission from star-forming regions provide a proper kinematic tracers for dynamical modeling of disc galaxy mergers. The process of photo-ionization does not affect the kinematics of the gas significantly. If the gas was moving with the bulk of the old stars before photo-ionization, we can safely assume that it would continue to do so after photo-ionization. Consequently, the star-forming H velocity would provide the velocity of the bulk of the baryons. In rare situations, when the gas was already displaced from the bulk of old stars before forming new stars, the star-forming H velocity would not match the velocity of stars. We find evidence for this phenomenon in UGC 12914 (See Figure 15) .
Other forms of photo-ionization such as post-Asymptotic Giant Branch (post-AGB) stars and AGNs also ionize metals along with the neutral hydrogen. The diffuse radiation of post-AGB stars can be responsible for low-ionization emission lines like [N II], [S II], and [O I] (Belfiore et al., 2016). The hard UV and/or X-ray radiation from AGN with energies of up to a few keVs may also produce high-ionization emission lines, such as [O III].
Ionization can also be the result of conversion of kinetic energy into heat rather than photo-ionization. The kinetic energy originating from the colliding gas flows produce shocks. Shocks are disturbances that move faster that the speed of sound in the ambient medium. Shock fronts dissipate the kinetic energy and convert it into heat. The shock-heated gas cools via radiation which is itself a powerful source of ionizing photons. Photons from the post-shocked gas may travel upstream and photo-ionize the pre-shock material known as the precursor. The emission from the shock and the precursor results in high emission line ratios, particularly among the low-ionization species, such as [N II]. Brightness of the radiative shocks and its precursor scales with the rate of dissipation of kinetic energy; resultingly, one expects to see higher velocity dispersion in the emission lines from brighter shocks (Dopita & Sutherland, 1995).
Unlike photo-ionized regions around young stars, shocked gas is not a good kinematic tracer for dynamical modeling of the gravitational features like tidal tails. Processes like high speed stellar wind, supernova (SN) feedback, and AGN-driven outflows, produce shock-heated ionized gas with velocities that can significantly differ from the stars and rest of baryons in their vicinity. The H velocity of shocks is not suitable for modeling the dynamics of tidal tails, especially with collisionless N-body simulations, in which gravity is the only player in the formation of tidal tails. Hence, we need to separate shocks from the photo-ionized gas near young stars. In this section, we describe our method for doing so. We investigate the how the fraction of shock-ionized gas evolve with merger stage. We also explore the correlation of shock fraction with merger mass ratio and the mass of the companion.
4.1 Separating Shocked and Star-Forming Regions
The diagrams first proposed by Baldwin et al. 1981, known as Baldwin, Phillips & Terlevich (BPT) diagram, are used to separate nuclear star formation from AGNs/Seyferts and low-ionization nuclear emission line regions (LINERs, Kauffmann et al. 2003; Kewley et al. 2006). In some galaxies, the emission line ratios are consistent with a combination of ionization by star formation and AGN, so they are called the composite galaxies. Integral Field Spectroscopy (IFS) has made it possible to investigate emission line ratios not only in the center but also throughout the galaxies.
IFS observations have shown that high low-ionization line ratios such as [N II]/H may be found in extra-nuclear regions. Using the data from MaNGA IFU survey (Bundy et al., 2015), Belfiore et al. (2016) argued that most of the extra-nuclear LINER-like emission is due to ionization by post-AGB stars in old stellar populations. Post-AGB ionized gas typically has low H EW (3Å). In this work, however, we disregard all components with emission lines 7Å in order to alleviate the issue of the underlying H absorption (see §3.3), so it is less likely to find post-AGB ionized gas in our systems. Besides, Belfiore et al. (2016) showed that high [N II]/H is also indicative of shocks, and in order to distinguish them from other sources of ionization we need to look at the velocity dispersion of emission lines.
It has been shown that given high enough velocity resolution, the velocity dispersion of emission lines is also indicative of the source of ionization, particularly of shocks. Monreal-Ibero et al. (2006) showed that in a shock-heated gas, the velocity dispersion of emission lines is correlated with low-ionization line ratios particularly with [O I]/H and [N II]/H. Monreal-Ibero et al. (2010) and Rich et al. (2011) confirmed these results, arguing that velocity dispersion can be used as an independent indicator of shocks. Rich et al. (2011) showed that the flux weighted histogram of velocity dispersion of emission lines in their high spectral resolution IFU data of LIRGs reveal a bi-modality. It displays a peak at low velocity dispersion of a few kilometers per second corresponding to photo-ionization by stars, and a bump at velocity dispersions higher that 100 km/s. Rich et al. (2014) showed that composite emission in the extra-nuclear regions of galaxies may originate from the combination of shocks and star formation rather than AGN and star formation. Based on these observations, Rich et al. (2015) proposed a limit of km/s for velocity dispersion of emission from star-forming regions. They suggested that components with km/s are emitted from low velocity shocks.
As described in §2 our observations were carried out in the wavelength range 6050-7000Å and did not include H and [O III] emission lines, so we can not use the BPT diagnostic diagrams to separate the emission lines. On the other hand, our SparsePak observations have a relatively high velocity resolution of 30 km s, making it possible to use the width of emission lines to separate star formation and shocks.
In the top plot of Figure 4 we show the H flux weighted histogram of velocity dispersion for all 956 high-S/N (3) components in high EW fibers (7Å) in all of the observed systems. Most of the H flux from our sample is emitted at velocity dispersions below 90 km/s, which has two peaks that may correspond to normal and turbulent star formation. There is also a small bump in the H flux around the velocity dispersion of km/s.
Combining the available information of [N II]/H and width of emission line, we use plots of log([N II]/H) vs. velocity dispersion to separate the shocked regions from the star-forming ones. Based on the curve suggested by Kewley et al. (2006) for separating star-forming galaxies from composites, spectra with log([N II]/H)-0.2 are more likely to be star-forming, and spectra with log([N II]/H)-0.2 are more likely to be composite or AGN. Rich et al. (2014) argued that extra-nuclear composite-like spectra may be due to shocks. On the other hand, for velocity dispersion, Rich et al. (2015) proposed a limit of 90 km/s to separate star formation and shocks. Like the limit of line-ratio, this limit is not rigid. On the low end, components are more likely to be star-forming, and on the high end the opposite is true. In Mortazavi et al. (2018), we used the plot of log([N II]/H) vs. velocity dispersion for separating emission lines observed in the Mice galaxies. We divided the emission line components visually into three groups. The lines separating these groups are plotted in the lower panel of Figures 4, as well as Figures 9 and 11. Group 1 has low [N II]/H and low velocity dispersion. Both criteria indicate that Group 1 components are likely to be emitted from the star-forming regions. In Mortazavi et al. (2018), we used these components to produce the H velocity map that was in excellent agreement with HI velocity map of the Mice galaxies. Also, using the [O III] and H lines obtained from the CALIFA data of the Mice (Sánchez et al., 2016) we showed that these components are indeed in the star-forming regions of the BPT diagram (Kewley et al., 2006). Group 2 components have both higher [N II]/H and velocity dispersion, suggesting that shocks are their likely source of ionization. Again, In Mortazavi et al. (2018), we used CALIFA data to show that the components of this group are more likely to be in the composite and AGN regions of [N II] BPT diagram. Group 3 components have high velocity dispersion, but relatively lower [N II]/H. In Mortazavi et al. (2018), we took them as overlapping unresolved components of Groups 1 and 2. In this work, however, we consider Group 3 components as shocks along with Group 2. In most of the maps of [N II]/H in Figure 6, we find Group 3 components in extra-nuclear regions, where it is less likely that they originate from multiple unresolved velocity components. Moreover, by visually inspecting the emission lines, we confirm that most of them are indeed single high velocity dispersion components originated at shocks.
We take Group 1 as star-forming regions and Groups 2 and 3 as shocked gas. Selection of star-forming regions is fairly conservative, but for shocks one should keep in mind that high low-ionization line ratio or broad emission line can be the result of other ionizing processes such as hard AGN UV/X-ray emission and post-AGB stars (Belfiore et al., 2016). In §4.3 we will discuss how shocks are established as the source of ionization using the [N II]/H maps of Group 2 and 3 components.
Table 2 shows some of the quantities measured for the observed systems. The light ratios are mostly obtained from the K-band magnitudes of 2MASS survey (Skrutskie et al., 2006). The projected separation between the cores of the two galaxies is measured using the redshift distance and the angular separation. In this table, we also present the minimum H EW among fibers with H S/N3. We use an EW low cut of 7 Å in all systems as discussed in §3.3. In addition, we present H shock fraction ,. For each system, this value is obtained by dividing the sum of the H flux of Groups 2 and 3 by the total H flux of all fibers. This measurement of is an approximate estimate as the dividing lines between the groups are not rigid. Moreover, our spatial coverage is not complete, and we may have missed regions with effectively high or low shock fraction laying in between the sparsely positioned fibers. Furthermore, SparsePak has a denser grid in the center which is usually placed on the cores of the galaxies, so our measurement of H shock fraction is biased toward its value near the cores of galaxies, though most of the H flux is usually coming from the core anyway. We estimate the uncertainty of () by bootstrap sampling of data points in the log([N II]/H)-velocity dispersion space, based on the measured uncertainty of log([N II]/H) and velocity dispersion. For each sample, we add up the H flux of Groups 2 and 3 and divide by the total H flux. The standard deviation of the H shock fraction in 100 samplings is reported as .
|system name||light ratio||separation (kpc)||minimum H EW(Å)|
4.2 Velocity Maps
Figure 5 shows maps of the velocity of star-forming components in all of the observed systems. Velocity of H emission is not much affected by the underlying absorption, so we plot these velocity maps for all components with with S/N without the H EW cut. In 13 out of 19 separate pairs we find smooth rotation in both galaxies confirming that these components move with the bulk of the baryon in their vicinity and can be used to trace velocity of material governed by gravity. Four systems only show rotation in one galaxy. The secondary companion in UGC 480, Arp 273, and Arp 181 do not display many star-forming regions. In Arp 283 the primary galaxy shows no H emission in the outskirts, and has a disturbed velocity field in the center. NGC 5278/9 shows a disturbed rotation, and UGC 12914 displays a messy velocity field in a bridge of star-forming regions between the two galaxies, consistent with HI and CO observations. We will discuss some interesting features observed in two of the observed systems in the Appendix.
In some of the velocity fields one can see a void in the center of the galaxy. This can be due to high velocity dispersion, or [N II]/H of components near the center that has classified them in Groups 2 or 3. Also low H EW in the center of some of the galaxies indicate an old or post-starburst stellar population. Among the three coalesced systems, we only see smooth rotation in NGC 3509. This system has a tiny H shock fraction. In contrast, the other coalesced systems, NGC 2623 and NGC 3921, display very high shock fraction, few star-forming components, and no clear rotation.
4.3 Indications for Galaxy-Wide Shocks from [N II]/H Maps
IFS observations of systems with unambiguous AGN usually display a gradient of emission line ratios from center to edge. Hard ionizing radiation from AGN affects the gas in its vicinity more than the gas that is kpcs away in the disc or above it. In the presence of an AGN in the center one expects to find higher [N II]/H near the center than in the outskirts (Davies et al., 2014; Leslie et al., 2014). In case of shocks, on the other hand, depending on the shock producing mechanism the regions of high velocity dispersion and [N II]/H are not necessarily near the center, but can be anywhere in the galaxy. Resultingly, if we see no gradient of [N II]/H toward the center, we can argue that the source of ionization is likely to be shocks. Nevertheless, even if we find a gradient of [N II]/H toward the center, it is not inevitably due to hard ionizing radiation from an AGN. If the source of shocks are the superwinds or outflows from a central processes such as starbursts or AGNs, they can also produce a gradients in the observed [N II]/H. The mere existence of a gradient in [N II]/H is not an indication of an AGN.
Based on the arguments above, we look for indications of shocks in [N II]/H maps. None of these systems are unambiguous AGNs according to the literature.
4.4 Shocked Gas and Merger Sequence
Figure (a)a shows the plot of H shock fraction versus projected separation between galaxies in each pair. Assuming that the separation is an indication of merger stage, close pairs are at a later stage, closer to coalescence. One may recognize the general trend of increasing shock fraction with decreasing projected separation. Similar to Rich et al. (2015) we bin our systems into close pairs, wide pairs, and coalesced mergers. The close (wide) pairs are non-coalesced systems with projected separations less (more) than 25 kpc. We select the 25 kpc limit to have significant statistics in both bins. Figure (b)b shows the trend in these bins. In the bins of close and wide pairs, we also separate systems based on their stellar mass ratio (K-bank light ratio). We can see that not only are the close pairs more shocked than the wide pairs, but also more equal mass mergers reveal a higher fraction of shocks altogether. Similar to Rich et al. (2015) we find the highest shock fraction in the coalesced systems.
This trend is consistent with the flux-weighted histograms of velocity dispersions presented in Figure (a)a. The coalesced systems show the highest bump at high velocity dispersion, and the close pairs emit more high velocity dispersion H flux compared to the wide pairs.
Another way to see this phenomenon is to look at the plot of log([N II]/H) vs. velocity dispersion, separating components from the close and wide pairs as in Figure 9. We see a significant increase in the fraction of components in the regions of Groups 2 and 3 in the close pairs compared to the wide pairs, a signature of enhanced shock ionization, which confirms the results of Figure (b)b. However, this plot provides extra information; we see that compared to the wide pairs in the close pairs not only is the percentage of components in Groups 2 and 3 higher, but also the components in Group 1 generally reveal a higher log([N II]/H) and velocity dispersion, indicating that the environment of the ionized gas in the close pairs is more turbulent even in star-forming regions. Different light ratio bins are also separated in Figure9 with red, orange, and yellow points referring to equal-mass, non-equal-mass major, and minor mergers, respectively. One can see that particularly in the close pairs more equal mass components (red points) occupy a slightly higher region compared to the other two light ratio bins. This will be discussed in §4.5
We should keep in mind that a small projected separation does not necessarily imply a late stage merger. There are other factors than can affect the projected separation. One is the geometry of observation, for we may find a pair of widely separated galaxies that appear to be close to each other because of the geometry of our line of sight. In addition, two interacting galaxies just after the first passage look as close to each other as a pair that is about to coalesce. Moreover, distance between two interacting galaxies at the pericenter of their orbit, , affects how much separated they appear to an observer. For better constraints on the actual merger stage we need to model the dynamics of the merger. We will discuss this in §7.
4.5 Shocked Gas and Merger Mass Ratio
One of the benefits of having a sample of early stage mergers before coalescence is that we can explore how properties of the companion affect the other galaxy. We can find out about the properties of isolated galaxies from their cores which, supposedly, have not lost much of its stellar mass, morphology, and kinematics after the first passage. One of the properties that is relatively easy to measure for non-coalesced mergers is mass ratio.
The mass ratio of a galaxy merger can be estimated from the stellar mass ratio, assuming that the halo mass grows with stellar mass. This assumption is not accurate all the time (Behroozi et al., 2013), but for the sake of simplicity we hold on to it throughout this paper. We use the Ks-band magnitude from 2MASS survey to estimate stellar mass ratio between the galaxies of each pair, so in this paper, we use the terms light ratio and mass ratio interchangeably. Figure (a)a displays the plot of shock fraction vs. light ratio, , in the non-coalesced pairs of our sample. Similar to projected separation, we separate galaxies into 3 bins of light ratio shown by vertical dashed lines: representing the equal mass mergers, representing major galaxy merger of non-equal mass, and representing minor galaxy mergers. The dividing values are chosen to have a significant number of galaxies in each bin. Figure (b)b shows the average shock fraction in each bin. Average shock fraction in both galaxies as well as in the primary and the secondary galaxies are shown separately. It is slightly but significantly enhanced in more equal mass mergers; though, we should be cautious about this interpretation as in such low statistics either removing one or two systems or slightly changing the limits of the bins could affect the results. Furthermore, in Figure (b)b the error-bars are too large to provide a robust Conclusion about the difference of shock fraction in the primary (more massive) vs. secondary (less massive) companions. Though it appears that in minor mergers the shock fraction in the primary galaxies are significantly higher than the secondary ones.
Figure 11 shows the components in light ratio bins on plots of [N II]/H vs. velocity dispersion. Similar to Figure 9, as we go from right to left not only does the fraction of points in Groups 2 and 3 increase but also the points in Group 1 reveal a slightly higher log([N II]/H) and velocity dispersion. This suggests that the environment of more equal mass mergers is, on average, more turbulent even in star forming regions. Comparing the panels of Figure 9 we find that this effect is stronger in the close pairs compared to the wide pairs, for in the panel of close pairs, more equal mass mergers (red points) reveal a significantly higher log([N II]/H) and velocity dispersion. In Figure 11 black and cyan points display the components of the primary and secondary companions, respectively. In the minor merger (right panel) the black points appear to spread significantly higher than the cyan points even in the Group 1 region, revealing a slightly higher log([N II]/H). Consistent with Figure (b)b. This suggests that the environment of the primary companions are more violent that the secondary ones in minor mergers. We will discuss these results in §7.
5 Encounter Parameters of Equal Mass Mergers
The morphology of the observed systems provides some insights for their encounter parameters. All of these mergers are the result of interaction between rotation-supported disc-dominated galaxies, for all of them display strong elongated tidal features. These features would not have survived if the interacting galaxies were dispersion-supported. We take all pairs in our sample as being observed after the first passage, because tidal features seen in these systems must have been produced by the tidal impulsive force that the galaxies experienced during a close passage. Tidal forces produce deformations that grow with time. Mass striping by tidal features and dynamical friction remove angular momentum from the orbit, causing it to decay rapidly. The galaxies coalesce in a much shorter time after the second passage, compared to the time between the first and second passages. As a result, we can confidently argue that most of the non-coalesced pairs in our sample are between the first and the second passages.
We use the method described in Mortazavi et al. (2016) for measuring the encounter parameter of galaxy mergers. This method uses identikit (Barnes & Hibbard, 2009; Barnes, 2011), which is a software package for modeling the dynamics of interacting pairs of disc galaxies. identikit models of isolated galaxies contain a spherically symmetric distribution of massive collisionless particles that reproduces the potential of dark matter halo, and stellar bulge and disc. In addition, in order to represent the visible morphology and kinematics of the discs, the models include massless test particles initially in circular orbits. As the galaxies interact massless test particles move through the potential made by the massive particles. Test particles have been shown to be successful in producing the large scale tidal features in the disc-disc galaxy merger (Dubinski et al., 1999). Since they are massless, multiple discs in a single galaxy will not affect each other, so we can model different disc orientations in a single simulation run. This facilitates a relatively rapid exploration of encounter parameter space to find the best model that matches the data. Mortazavi et al. (2016) only tested equal mass identikit models, so in this work we only apply this method on the systems in the first light ratio bin in §4.5 with Ks-band light ratio 1.85.
We often need to measure and model the line of sight velocity of these tidal features, in addition to their apparent shape, to ensure the uniqueness of the model (Barnes, 2011). Collisionless stars are the ideal components to match with collisionless N-body simulations such as identikit. Nevertheless, velocity of stars are hard to measure in the faint tidal tails, as they require high continuum S/N. Another option is cold neutral gas, which is usually more extended in the discs of isolated galaxies and produce stronger tidal features in interactions. However, detecting cold gas through HI 21 cm emission line at high enough resolution requires long exposure times on large interferometers like JVLA, which is also expensive. Nebular emission, on the other hand, is relatively easy to measure for star-forming galaxies even in low surface brightness regions, and galaxy mergers usually enhance star formation in both the center and the tidal features (Jog & Solomon, 1992; Hattori et al., 2004; Whitmore & Schweizer, 1995; de Grijs et al., 2003). Mortazavi et al. (2018) obtained consistent results with velocities derived from H and HI emission lines for the dynamical modeling of the Mice merger system. However, it is important to separate star-forming regions from shocks before using the H derived velocities.
Velocity of H emission originating in star-forming regions is more suitable for dynamical modeling compared to the velocity of shocks. The baryons in the tidal features of interacting galaxies are mainly affected by the gravitational forces within and between the two galaxies. For dynamical modeling using collisionless N-body simulations, we need to match our model with the velocity of the bulk of the baryons, which are usually the collisionless stars. H, however, traces the ionized gas. Assuming that the original neutral gas had been moving along with the stars before ionization, the ionization process determines weather or not the ionized gas is still at the same velocity as the stars. In normal star-forming regions, the velocity of the ionized gas is not much affected by the ionization process, for UV emission from young and hot stars photo-ionize the surrounding ISM, making the HII regions. So, most of the time, Hemission from star-forming regions should be a proper kinematic tracer for dynamical modeling with N-body simulations. On the other hand, shock-ionization usually changes the velocity of neutral gas significantly. Starbursts and SNe produce winds that heat up the ISM with shocks and produce H emission at velocities different from the bulk of the baryons in their vicinity (Sharp & Bland-Hawthorn, 2010). Moreover, colliding gaseous clouds of merging galaxies induce shocks in the ISM (Soto & Martin, 2012). Galaxy mergers can also trigger or enhance AGNs (Ellison et al., 2011), which produce outflows and spread shocks (Wild et al., 2014). Shock excitation results in a motion of the ionized gas that is not exclusively gained by gravitational forces, and we can not model the velocity of shocked gas with simple N-body simulations. Consequently, we only use the velocity of star-forming regions presented in Figure 5 for dynamical modeling.
identikit evaluates the match between model and data by calculating a parameter called “score”. In order to calculate the score we place phase space boxes on the tidal features of the interacting galaxies. identikit calculates the score for each model based on the number of test particles that populate these boxes. In this work, we use the same algorithm as in Mortazavi et al. (2018) for placing the boxes, and their size are determined by the diameter of the SparsePak fibers ().
Mortazavi et al. (2016) tested their modeling method against a set of hydrodynamical simulations of isolated galaxy mergers with known encounter parameters (Cox et al., 2008). The test results were classified into three levels of convergence based on the distribution of reconstructed parameters. In some tests the reconstructed encounter parameters had either little or controllable systematics and were reliable (good convergence). However, some other tests resulted in either large or unpredictable systematics, and unreliable reconstructed model parameters were obtained (fair and poor convergences). The quality of convergence was determined by the distribution of the reconstructed initial orientations of the discs. When this distribution is a narrow, it is concluded that a subset of identikit models with close initial orientations have resulted in the best scores. On the other hand, a flat distributions of reconstructed disc orientations suggests that models with different disc orientations have all resulted in good scores, and consequently, the reconstructed encounter parameters are not unique. Mortazavi et al. (2016) found good convergence in most of the tests on discs that were almost in a prograde orbit. In contrast, they found poor and fair convergence when any of the two interacting discs were in a retrograde or polar orbit.
|system name||light||separation||previous models|
|NGC 5257/8||1.17||33.8||Privon et al. (2013)|
|UGC 12914||1.61||14.8||Vollmer et al. (2012)|
Table 3 shows the systems we selected for dynamical modeling. Here we exclude NGC 4676, the Mice, which we modeled and discussed extensively in Mortazavi et al. (2018). The systems in Table 3 are the most equal mass systems based on their 2MASS K-band magnitude (Skrutskie et al., 2006), which is used as a crude estimate of stellar mass. We select the galaxies with a light ratio less that 1.85, and make an attempt to model them with equal mass identikit models. For two of these systems we found prior dynamical models in the literature, which are mentioned in Table 3.
We found a first ever model for UGC 07593, though the dynamical modeling method we used did not well converge for the other systems. Here we present the encounter parameters we found for UGC 07593 and discuss briefly the possible reasons for the failure of our method for some of the other systems.
Figure 12 shows slices of the map of scores across three encounter parameters, eccentricity, pericentric distance, and time since pericenter. The slices are taken at the best-fit parameters shown by the cyan box. This system is in a relatively late stage at of the time between the first and the second passages. The reconstructed pericentric distance and eccentricity are and , respectively. Taking the physical length and velocity scaling into account we derive the encounter parameters with physical units (See Mortazavi et al. 2016). The reconstructed time since the first passage is Myr. The physical pericentric distance is kpc. Both galaxies are almost half-way between a prograde and a polar orbit with the inclinations being at about from the orbital plane.
UGC 12914: The velocity map of this system, shown in Figure 5, point in the direction of why our method cannot reconstruct the observed features. Fibers with strong H line emission are scattered in the region between the two galaxies, and the velocity of Group 1 components do not show an smooth rotation. The gas in the bridge contains a significant amount of shocks. In addition, powerful HI, CO, radio continuum, and X-ray emission is detected in the bridge, suggesting that the gas originally in the two discs have been stripped off as a result of a virtually head on encounter (Condon et al., 1993; Braine et al., 2003; Gao et al., 2003; Vollmer et al., 2012; Appleton et al., 2015). It does not accompany the rest of the baryons, particularly the collisionless stars, anymore and it can not be reconstructed by identikit test particles.
NGC 5257/8: NGC 5257 has a prominent tail in HI that is absent in the optical image. This tail played an important role in the modeling by Privon et al. (2013). Without the information in the extended tail we do not have the same constraints to pin down the encounter parameters. Furthermore, the model presented by Privon et al. (2013) indicates that NGC 5258 is in a nearly polar orbit. According to Mortazavi et al. (2016), the modeling method we used does not result in a good convergence when any of the galaxies is in a polar or retrograde orbit, so assuming that the model of Privon et al. (2013) is correct we should have expected a similar failure.
Arp 87: This system consists of an edge-on and a nearly face-on galaxy, so even without a reconstructed model of the orbit, we can be certain that one of the galaxies is closer to a polar or retrograde orbit rather than a prograde one, and the limitation of our method in modeling mergers with polar and retrograde orbits applies. Moreover, even though we find a smooth rotation in the H velocity near the center of both discs, the tidal tails do not display any detectable H emission. Lack of velocity information in the tails introduces more degeneracy to the reconstructed dynamical models.
6 Shock Fraction vs. Reconstructed Encounter Parameters
|system name||source||time (Myr)||t (Myr)||(kpc)|
|NGC 5257/8||Privon et al. (2013)||230.0||1200.0||21.0||0.15|
|NGC 4676||Mortazavi et al. (2018)||190.0||775.0||18.0||0.26|
|UGC 12914||Vollmer et al. (2012)||26.0||-||1.2||0.37|
|Arp 284||Struck & Smith (2003)||170.0||-||6.5||0.39|
|UGC 07593||this work||27.0||12.0||2.5||0.40|
|NGC 2623||Privon et al. (2013)||220.0||-80.0||1.0||0.90|
In this work, we measured the H kinematics of a sample of 22 galaxy mergers. Nineteen of these mergers are in the stage between the first passage and the coalescence, and the other three are already coalesced. We examine the source of ionization in these mergers using the kinematics and flux ratio of H and [N II] emission lines. We find that the close pairs with projected separation km/s have a higher H shock fraction compared to the wide pairs with projected separation km/s. The average of the shock fraction in the coalesced systems is even higher; shocks are responsible for an average of about 50% of H emission in these systems. These results are consistent with results of Rich et al. (2015) suggesting that if the sequence of wide pairs, close pairs, and coalesced mergers is a time sequence, then the shock fraction is growing as the merger proceeds. However, dynamical modeling, which incorporates the complex morphology and kinematics of tidal features, is required to constrain encounter parameters in most galaxy mergers. Determining merger stage only using the separation between the pairs is not alway correct. The geometry of the encounter and observer can result in a close projected separation for a physically wide pair. A pair right after the first passage appears as close as one near coalescence. The impact parameter of the encounter also affects the appearance of projected separation during the encounter. In order to obtain a robust correlation between shock fraction and merger stage we need reconstructed dynamical models of galaxy mergers.
In addition to the model of UGC 07593 from this work, we can find the dynamical models of five other systems in our sample. The Mice galaxies (NGC 4676) are extensively modeled in the literature (Toomre & Toomre, 1972; Barnes, 2004; Privon et al., 2013; Mortazavi et al., 2018). NGC 5257/8 and NGC 2623 are among the four systems modeled in Privon et al. (2013). We also find a dynamical model of UGC 12914 and Arp 284 in Vollmer et al. (2012) and Struck & Smith (2003), respectively. The reconstructed time since the first passage, time left to coalescence, and pericentric separation of these models are shown in Table 4 along with the measured shocked H fraction from this work.
Table 4 is sorted by shock fraction. One can see that time since pericenter varies with no obvious trend, but time till coalescence, t, displays a trend; systems with long times until coalescence have smaller shock fraction compared to the systems near coalescence or after coalescence. Vollmer et al. (2012) and Struck & Smith (2003) do not provide t in their model. Pericentric separation is also almost sorted in Table 4. Systems with wide pericentric distance have less fraction of shocked H than systems with small pericentric distance. This trend and the trend of time till coalescence can be seen in Figure 13.
Production of shocks is the result of complex processes and can not be reduced to simple correlations like the ones in Figure 13. In particular, during the coalescence, violent relaxation, random orbital crossings, and the return of tidally striped material may agitate the ISM so much that the shocked gas in the remnant loses memory of the orbit of merging galaxies. However, before coalescence, between the first and second passages, the processes that produce shocks can generally be categorized into two modes (Cox et al., 2004b; Soto & Martin, 2012).
The first mode is initiated by a tidal impulse. The tidal impulse at the time of the first pericenter derives the rotating gas in the discs to flow in toward the center of the interacting galaxies (Mihos & Hernquist, 1994; Barnes & Hernquist, 1996). The inflowing gas may collide with the gas at lower orbits, and if the collision is fast enough, shocks can be produced. The inflowing gas may also reach the core of the galaxies and trigger or enhance starbursts or AGN, which produce fast outflowing material that can collide with the ISM on their way out and produce shocks. In this mode shocks are produced a few dynamical times after the tidal impulse kicks in. This mode is a side effect of the gravitational tidal impulse during a close passage, so factors that affect the strength of the tidal force should affect the shock production through this mode.
In this mode of shock production, pericentric distance and initial disc orientation should strongly affect the amount shocks. We know that the tidal response is proportional to the inverse cube of the distance to the perturber, so a factor of two in pericentric separation changes the tidal impulse by a factor of eight. In prograde encounters tidally induced inflows are even more prominent, because the tidal force pushes or pulls the rotating material in the same direction during the pericentric passage, and the resulting effect of tides is maximized. Therefore, in prograde mergers with small pericentric separation, we expect to find stronger shocks produced through the first mode, and the amount of shock should increase with time since pericenter.
In a second mode, sometimes, the geometry of galactic encounter allows the ISM of the two gas rich discs to directly collide with each other. For this to take place, the two discs should have small enough distance at the pericenter, so that the discs and not just the dark matter halos cross during the first passage. The shocks are produced immediately, and decay as they cool down through radiation. Unlike the other mode of shock production, this one is less likely to happen in prograde discs, because the overlapping rotating material moves in the same direction in prograde discs. Therefore, In this mode of shock production, one expects to find more shocks in retrograde and polar discs interacting with small pericentric separation, and these shocks should become less prominent with time since pericenter.
7.1 Shock Production and Encounter Parameters
Our sample in this work was selected based on the visibility of strong tidal features (see §2.2), so it should be biased toward prograde mergers, and the first mode of shock production should be dominant. We predict higher shock fraction in systems that are closer to coalescence, because not only is there more time for the gas inflows to enhance central starbursts/AGNs and produce shocks through outflows, but also violent relaxation, random orbital crossings, and the return of tidally striped material during coalescence may generate shocks. This trend can be seen in the left panel of Figure 13. Notice that in this panel we do not have the data for UGC 12914 and Arp 284.
From the spacial distribution of the H emission, it seems that in UGC 12914 shocks are produced through the second mode. The presence of the continuous gaseous bridge observed on different bands (Condon et al., 1993; Braine et al., 2003; Gao et al., 2003; Vollmer et al., 2012; Appleton et al., 2015) implies that gas has been striped from the two galaxies because of a recent head-on encounter. According to the model presented by Vollmer et al. (2012), the relative velocity of galaxies at the time of the encounter was 1000 km/s. The gas clouds in the ISM of the two discs have collided at this speed and dragged themselves away from the stars toward the bridge. If Vollmer et al. (2012) had provided the time until final coalescence, we would expect it to be a long time, for the galaxies are still receding from each other. This would be inconsistent with the trend in the left panel of Figure 13 between shock fraction and time left to coalescence, and it is an additional evidence that the shocks in this galaxy are produced through a fundamentally different procedure.
In both proposed modes of shock production, pericentric separation affects the shock fraction in the same way. The smaller the pericentric separation is, the stronger the induced tidal forces are, and the more shocks via the first mode are expected. Similarly, a small enough pericentric separation results in direct collision of gas in the ISM of the two discs, producing shock. Both of these are consistent with the trend seen in the right panel of Figure 13. We should keep in mind that this description only applies to the time between the first and the second passages. During coalescence, shock production is more complicated and may not be characterized by these simplified modes. In fact, strong shocks may be found in the remnant of a recent coalesced merger, even if it had a large separation at the first passage.
7.2 The Effect of Mass Ratio
Assuming that the first mode of shock production is dominant in our sample, the shocks should have been created mostly as a result of the mass flows generated by tidal forces during the first passage; consequently, we should find stronger shocks in a galaxy for which the companion is more massive, for tidal force depends on the mass of the perturber. Disregarding other factors, this leads to two predictions about the effect of mass ratio on shock fraction in binary galaxy mergers: First, when galaxies have comparable mass, tidal force on both galaxies should be relatively strong compared to the the self-gravity of individual galaxies. In comparison, when one of the companions is much less massive, it may not induce a strong enough tidal force to overcome the self-gravity of its companion and affect the orbit of rotating disc particles, so only one of the companions (the secondary) would experience strong enough tides to create shocks. Therefore, more equal mass mergers are expected to have a higher overall shock fraction. Second, with a similar argument for non-equal mass mergers, we should have seen a higher shock fraction in the smaller companion as it experiences a larger tidal impulse from a more massive perturber.
The results of §4.5 confirms the first prediction, but for the second prediction, the evidence points to the contrary. It is important to note that the predictions above would have been valid only if all other factors affecting the amount of shocks were the same for all galaxies, and all shocks were created from the first mode. However, as discussed earlier, at least for UGC 12914 most of the shocks are not likely to be the consequence of a tidal impulse. Besides, there are other important factors such as gas content and its gravitational stability that may affect the fraction of shocks produced by tidal forces. Without an independent measurement of the gas content of each galaxy we cannot evaluate the effect of these other factors. From a statistical point of view, the scaling relationship between the size, baryonic mass, and gas fraction of nearby spiral galaxies suggest that those with a total baryonic mass have the most gas mass reservoir of all star-forming disc galaxies (Wu, 2018). On the other hand, gravitational instability of the gas disc may be a critical player in shock production when the tidal impulse kicks in. According to the Toomre’s stability criterion, disc stability against collapse is inversely proportional to gas surface density (Safronov, 1960; Toomre & A., 1964), and in the nearby Universe, gas surface density is generally maximum for discs with total baryonic mass (Wu, 2018). The total baryonic masses of our galaxies are likely to be on either side of these maximum values, so without an accurate constraint on the mass of each galaxy, we cannot evaluate if the effects of gas content explain the low shock fraction of the secondary companions in our minor mergers. Measuring the baryonic mass of the observed galaxies is out of the scope of this work.
It is important to emphasize that in this work, all of the relationships between shock fraction and encounter parameters are based on either a small sample of 22 observed systems or an even smaller sample of 6 systems with available dynamical models. In order to independently investigate the effect of each parameter and possibly find more complicated relations, we need a larger statistical sample of galaxy mergers observed with optical IFU instruments. Such observations could be used to understand the possible role of merger-induced shocks in the overall quenching of star formation in the Universe. Both modes of shock production discussed here, remove or heat up the cold gas reservoir in the discs which is the fuel for star formation. Hence, in a simplified picture, we expect galaxies to quench star formation after a period of spreading shocks throughout their gas content. Whether shocks are made by strong galactic winds or by the direct collision of gas in the discs, they drag the cold gas away into the halo (Cox et al., 2006; Oppenheimer et al., 2010). This gas may not be able to return to its cold state due to the gravitational shock heating (Birnboim & Dekel, 2011). Ongoing and future IFU galaxy surveys (e.g. SAMI: Croom et al. 2012, MaNGA: Bundy et al. 2015, etc.) provide a promising area for further investigation in near future.
In this work we observed 22 galaxy mergers with strong tidal features, using the SparsePak IFU on the WIYN telescope at KPNO. We reduce the data and analyze the emission lines with one- and two-component Gaussians. We use an MCMC code to estimate the uncertainty of fit parameters, and select the best number of components using the F-test. Relatively high spectral resolution of our data allows us to use velocity dispersion of emission lines along with [N II]/H line ratios to separate H originating at star-forming regions from that arising from shocked gas. We use emission line maps to confirm that most of high velocity dispersion and high [N II]/H components are galaxy wide shocks, likely to be induced as a result of interaction. We found that the fraction of H emission from shocked gas is correlated with the separation of galaxies in pairs. Close pairs have higher shock fraction than wide pairs. The three coalesced systems show the highest average shocked H fraction.
We use the modeling method developed in Mortazavi et al. (2016) along with the velocity maps of star-forming regions to model the dynamics of equal mass pairs in our sample. We find the first ever constraints on the encounter parameters of UGC 07593, but obtain poor convergence in the other four attempts. We find dynamical model of five other systems in our sample from the literature, and explore the correlations between some encounter parameters and shock fraction. In these systems time left to coalescence and pericentric separation appear to be correlated with the fraction of shocked H. We suggest two modes of shock production that are responsible for most of the shocks in the early stages of a merger, after the first passage and before the coalescence.
The authors would like to thank Colin Norman, Joshua Barnes, David Law, Gregory Snyder, Susan Kassin, Ron Allen, and Luciana Bianchi for valuable discussions throughout this work. This project was supported in part by the Space Telescope Science Institute (STScI) Director’s Discretionary Research Fund (DDRF). This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575 (see Towns et al. 2014). We used the DSS scan of Arp 273 in this work. The DSS was produced at the STScI under U.S. Government grant NAG W-2166. The images of these surveys are based on photographic data obtained using the Oschin Schmidt Telescope on Palomar Mountain and the UK Schmidt Telescope. The plates were processed into the present compressed digital form with the permission of these institutions.
Appendix A Notes on Individual Systems
In the Appendix we discuss some of the interesting features we found in two of the systems observed in this work.
Arp 284: In the secondary galaxy (NGC 7715, east) we do not detect H emission in the western side of the edge-on disc which is toward the primary galaxy. This can be seen is shown in Figure (a)a. There are HII regions in the tidal bridge between the galaxies. Smith & Wallin (1992) found an HI bridge in the same position as the optical bridge. The primary galaxy (NGC 7714, west) displays smooth rotation in the star-forming components, which is scattered throughout the tilted disc.
This system is studied extensively in the literature. Delgado et al. (1998) used HST/GHRS ultraviolet spectroscopy and ground based radio, optical, and X-ray observations to perform a spectral synthesis modeling on central 300 pc of the primary galaxy. They suggest that the center of this galaxy contains a very young starburst (4.5 Myrs) along with an older stellar population with ages around tens of Myrs or older. Dynamical modeling by Struck & Smith (2003) indicates that a significant amount of gas has been transferring from the secondary companion to the the primary galaxy since the first passage, which occurred Myr ago. This gas transfer is fueling the starburst in the core of NGC 7714. It seems that tidal striping of gas and lack of inflows is why the secondary galaxy has quenched star formation. Smith et al. (2005) observed X-ray emission in the core of the the primary galaxy attributing it to mechanical energy injected into to the ISM by either supernovae or High Mass X-ray Binaries (HMXBs), both as results of the central starburst. In addition, Infrared spectroscopy with Spitzer space telescope shows no evidence of an obscured AGN in the center of the primary, confirming that it is a young unobscured starburst galaxy (Brandl et al., 2004).
In the primary galaxy, we find a vertical cone-like structure made mostly by Group 3 and some Group 2 components, which is shown in Figure (b)b. At first glance, we speculate that this structure and its low overall [N II]/H ratio indicates that they are the outflows from the central starburst proposed by Smith et al. (2005). However, by plotting the difference in the velocity of broad and narrow components in fibers where two-component fit is preferred, we come up with a new conclusion. Figure (c)c shows the velocity of narrow components subtracted from the velocity of broad components, so blue (yellow) circles indicate that the broad component is blueshifted (redshifted) with respect to the narrow component. Assuming that the narrow component corresponds to the star-forming region in the disc and that we only see the broad component that is physically in front of the narrow component, we can argue that blueshifted broad components show outflows and redshifted ones show inflows. In Figure (c)c we see blueshifted broad components on the wide end of the cone confirming our speculation of them showing outflows. However, on the narrow end of the cone closer to the center, we see redshifted broad components, which can be explained by the inflowing gas. This could be the gas that is striped from the secondary galaxy leaving it quiescent on its west side. This observation is consistent with the dynamical model of Struck & Smith (2003) which suggested that the secondary galaxy is depleting gas into the center of the primary one, causing a starburst and generating the outflows toward the west. This observation demonstrates that by properly separating kinematic components in high spectral resolution nebular emission data, not only outflows but also inflows may be detected.
UGC 12914/5: This system is also known as the “Taffy galaxies”, because HI and 4.86 GHz observations by Condon et al. (1993) showed a joined HI and radio continuum bridge between the two galaxies. Condon et al. (1993) suggest that cosmic rays, magnetic fields, and HI gas have been striped from the two galaxies as the result of a nearly head-on collision, about 20 Myrs ago. In the bridge, they found double peaks in the HI profile separated by 100-300 km/s. Braine et al. (2003) finds a significant amount of molecular gas emission (CO) in the bridge, suggesting that 18-35% of the total gas mass in the system sits in the bridge. The CO emission does not reveal the double peak and only the high-velocity HI component has a CO counterpart.
Vollmer et al. (2012) modeled the dynamics of this system using a model that includes collisionless halo+stars particles and collisional sticky gas particles. They distinguish the molecular gas from neutral hydrogen using a prescription for total gas density. They found the best-fit model at about 26 Myr after the first pass. Their model reproduced the morphology and kinematics of stars and gas in this system including the prominent gaseous bridge. They also reproduced the double-component HI profile and the single component CO profile in the bridge region. Changing the cloud-cloud collision parameter in their simulation, they argued that the double component profile is produced mainly by the distortion caused by the collisional nature of ISM, and not by the tidal distortions.
Appleton et al. (2015) used the Chandra observatory to show that the bridge also emits a significant amount of X-ray. They showed that 19% of the X-ray luminosity of the the system comes from the bridge. They also used Herschel Far-IR data to estimate the SFR in the bridge, concluding that SFR is too low to account for X-ray emission via outflows. Moreover, they showed that the peak of the diffuse X-ray emission does not match with the peak of the radio continuum, ruling out a direct connection between the X-ray and synchrotron emission caused by cosmic rays in the bridge. They suggest that the main source of X-ray emission in the bridge is shock heating due to collision of the ISM in the two galaxies as was suggested earlier by Struck (1997). The shocked gas can be heated up to K during the collision and cool down to K in 35 Myr and to K within less than 100 Myr. Based on their estimated time of Myr after the pericenter, the temperature of the observed soft X-ray is consistent with shocks with speed range of 430-570 km/s.
Figure (c)c shows the velocity difference between the broad and narrow component in fibers where two component fit is preferred by the F-test. It shows that most of the fibers in the bridge region prefer a two-component fit. In most of these fibers, the narrow component meets the criteria for a normal star-forming region, and the broad component is considered as shocked. In most of these two-component fits the broad component is redshifted with respect to the narrow component ( km/s), suggesting that outflows caused by star formation are not responsible for producing the shocks in this region. This velocity difference is consistent with the velocity difference between peaks in the HI profile in the bridge. Based on the velocity difference we can argue that the narrow star-forming components shown in Figure (a)a correspond to the HI components with corresponding CO lines, and the broad shocked components shown in Figure (b)b correspond to the HI components without corresponding CO lines. The velocity difference is a bit lower but consistent with shock speeds required to produce the X-ray emission in the bridge region. Therefore, our data confirms that the shocked components in the bridge are produced as a result of the fast direct collision of gas in the ISM at the time of the first passage.
In Figure (a)a we also find a void of star forming components in the center of UGC 12914. Appleton et al. (2015) found a slightly extended ultra-luminous X-ray source in this location, hinting to the presence of an obscured AGN. Though, the velocity dispersion of the shocked components are not as high as what is expected for AGN hosts (Rich et al., 2014).
- pubyear: 2017
- pagerange: Dynamics and Shocks from H Emission of Nearby Galaxy Mergers–A
- In some of the coalesced systems we dithered in order to find a better coverage and better handle on the spatial distribution of outflows.
- Arp 273, also known as the Rose Galaxies, is not in SDSS footprint. However, the Hubble space telescope (HST) image taken in 2010 resolves individual HII regions in both galaxies.
- We keep these fibers for velocity measurements as the underlying absorption does not significantly affect the velocity of the emission lines.
- Some galaxies like NGC 2623 are AGN candidates, but for all of them the controversy stays.
- Acker A., Köppen J., Samland M., Stenholm B., 1989, The Messenger, 58, 44
- Appleton P. N., et al., 2015, The Astrophysical Journal, 812, 118
- Arp H., 1966, Astrophysical Journal Supplement, 14, 1
- Baldwin J. A., Phillips M. M., Terlevich R., 1981, Publications of the Astronomical Society of the Pacific, 93, 5
- Barnes J. E., 2004, Monthly Notices of the Royal Astronomical Society, 350, 798
- Barnes J. E., 2011, Monthly Notices of the Royal Astronomical Society, 413, 2860
- Barnes J. E., Hernquist L., 1996, Astrophysical Journal v.471, 471, 115
- Barnes J. E., Hibbard J. E., 2009, The Astronomical Journal, 137, 3071
- Behroozi P. S., Wechsler R. H., Conroy C., 2013, The Astrophysical Journal, 770, 57
- Belfiore F., et al., 2016, Monthly Notices of the Royal Astronomical Society, 461, 3111
- Bell E. F., et al., 2004, The Astrophysical Journal, 608, 752
- Bershady M. A., Andersen D. R., Harker J., Ramsey L. W., Verheijen M. A. W., 2004, Publications of the Astronomical Society of the Pacific, 116, 565
- Birnboim Y., Dekel A., 2011, Monthly Notices of the Royal Astronomical Society, 415, 2566
- Braine J., Davoust E., Zhu M., Lisenfeld U., Motch C., Seaquist E. R., 2003, Astronomy and Astrophysics, 408, L13
- Brandl B. R., et al., 2004, The Astrophysical Journal Supplement Series, 154, 188
- Brown M. J. I., Dey A., Jannuzi B. T., Brand K., Benson A. J., Brodwin M., Croton D. J., Eisenhardt P. R., 2007, The Astrophysical Journal, 654, 858
- Bundy K., et al., 2015, The Astrophysical Journal, 798, 7
- Condon J. J., Helou G., Sanders D. B., Soifer B. T., 1993, Astronomical Journal (ISSN 0004-6256), 105, 1730
- Cox T. J., Primack J., Jonsson P., Somerville R. S., 2004a, The Astrophysical Journal, 607, L87
- Cox T. J., Primack J., Jonsson P., Somerville R. S., 2004b, The Astrophysical Journal, 607, L87
- Cox T. J., Jonsson P., Primack J. R., Somerville R. S., 2006, Monthly Notices of the Royal Astronomical Society, 373, 1013
- Cox T. J., Jonsson P., Somerville R. S., Primack J. R., Dekel A., 2008, Monthly Notices of the Royal Astronomical Society, 384, 386
- Croom S. M., et al., 2012, Monthly Notices of the Royal Astronomical Society, 421, 872
- Davies R. L., Kewley L. J., Ho I. T., Dopita M. A., 2014, Monthly Notices of the Royal Astronomical Society, 444, 3961
- Delgado R. M. G., Garc\’\ia-Vargas M. L., Goldader J., Leitherer C., Pasquali A., 1998, arXiv.org, pp 707–719
- Dopita M. A., Sutherland R. S., 1995, Astrophysical Journal v.455, 455, 468
- Dubinski J., Mihos J. C., Hernquist L., 1999, The Astrophysical Journal, 526, 607
- Ellison S. L., Patton D. R., Mendel J. T., Scudder J. M., 2011, Monthly Notices of the Royal Astronomical Society, 418, 2043
- Faber S. M., et al., 2007, The Astrophysical Journal, 665, 265
- Gao Y., Zhu M., Seaquist E. R., 2003, The Astronomical Journal, 126, 2171
- Hattori T., et al., 2004, The Astronomical Journal, 127, 736
- Hibbard J. E., Mihos J. C., 1995, Astronomical Journal v.110, 110, 140
- Holincheck A. J., et al., 2016, Monthly Notices of the Royal Astronomical Society, 459, 720
- Jog C. J., Solomon P. M., 1992, The Astrophysical Journal, 387, 152
- Karachentsev I., Lebedev V., Shcherbanovski A., 1985, Bulletin d’Information du Centre de Donnees Stellaires, 29, 87
- Kauffmann G., et al., 2003, Monthly Notices of the Royal Astronomical Society, 346, 1055
- Kewley L. J., Groves B., Kauffmann G., Heckman T., 2006, Monthly Notices of the Royal Astronomical Society, 372, 961
- Leslie S. K., Rich J. A., Kewley L. J., Dopita M. A., 2014, Monthly Notices of the Royal Astronomical Society, 444, 1842
- Lintott C., et al., 2010, Monthly Notices of the Royal Astronomical Society, 410, 166
- Lomax R. G., 2007, Statistical Concepts: A Second Course (2007), http://scholar.google.com/scholar?q=related:XxhhnoJj_UEJ:scholar.google.com/&hl=en&num=20&as_sdt=0,5
- Lotz J. M., Jonsson P., Cox T. J., Croton D., Primack J. R., Somerville R. S., Stewart K., 2011, The Astrophysical Journal, 742, 103
- Medling A. M., et al., 2015, arXiv.org, pp 2301–2311
- Mihos J. C., Hernquist L., 1994, The Astrophysical Journal, 425, L13
- Mihos J. C., Hernquist L., 1996, Astrophysical Journal v.464, 464, 641
- Monreal-Ibero A., Arribas S., Colina L., 2006, The Astrophysical Journal, 637, 138
- Monreal-Ibero A., V\’\ilchez J. M., Walsh J. R., Muñoz-Tuñón C., 2010, Astronomy and Astrophysics, 517, A27
- Mortazavi S. A., Lotz J. M., Barnes J. E., Snyder G. F., 2016, Monthly Notices of the Royal Astronomical Society, 455, 3058
- Mortazavi S. A., Lotz J. M., Barnes J. E., Privon G. C., Snyder G. F., 2018, Monthly Notices of the Royal Astronomical Society, 474, 3423
- Narayanan D., et al., 2008, The Astrophysical Journal Supplement Series, 176, 331
- Oppenheimer B. D., Davé R., Kereš D., Fardal M., Katz N., Kollmeier J. A., Weinberg D. H., 2010, Monthly Notices of the Royal Astronomical Society, 406, 2325
- Privon G. C., Barnes J. E., Evans A. S., Hibbard J. E., Yun M. S., Mazzarella J. M., Armus L., Surace J., 2013, The Astrophysical Journal, 771, 120
- Rich J. A., Kewley L. J., Dopita M. A., 2011, The Astrophysical Journal, 734, 87
- Rich J. A., Kewley L. J., Dopita M. A., 2014, The Astrophysical Journal, 781, L12
- Rich J. A., Kewley L. J., Dopita M. A., 2015, The Astrophysical Journal Supplement Series, 221, 28
- Rodriguez-Gomez V., et al., 2016, Monthly Notices of the Royal Astronomical Society, 458, 2371
- Safronov V. S., 1960, Annales d’Astrophysique, 23, 979
- Sánchez-Menguiano L., et al., 2016, arXiv.org, p. arXiv:1610.00440
- Sánchez S. F., et al., 2012, Astronomy and Astrophysics, 538, A8
- Sánchez S. F., et al., 2016, Astronomy & Astrophysics, 594, A36
- Sanders D. B., Mirabel I. F., 1996, Annual Review of Astronomy and Astrophysics, 34, 749
- Scudder J. M., Ellison S. L., Torrey P., Patton D. R., Mendel J. T., 2012, Monthly Notices of the Royal Astronomical Society, Volume 426, Issue 1, pp. 549-565., 426, 549
- Sharp R. G., Bland-Hawthorn J., 2010, The Astrophysical Journal, 711, 818
- Skrutskie M. F., et al., 2006, The Astronomical Journal, 131, 1163
- Smith B. J., Wallin J. F., 1992, Astrophysical Journal, 393, 544
- Smith B. J., Struck C., Nowak M. A., 2005, The Astronomical Journal, 129, 1350
- Soto K. T., Martin C. L., 2012, The Astrophysical Journal Supplement, 203, 3
- Struck C., 1997, The Astrophysical Journal Supplement Series, 113, 269
- Struck C., Smith B. J., 2003, The Astrophysical Journal, 589, 157
- Toomre A., 1977, Evolution of Galaxies and Stellar Populations, pp 401–
- Toomre A., A. 1964, The Astrophysical Journal, 139, 1217
- Toomre A., Toomre J., 1972, Astrophysical Journal, 178, 623
- Towns J., et al., 2014, Computing in Science & Engineering, 16, 62
- Vollmer B., Braine J., Soida M., 2012, Astronomy and Astrophysics, 547, A39
- Whitmore B. C., Schweizer F., 1995, The Astronomical Journal, 109, 960
- Wild V., et al., 2014, Astronomy and Astrophysics, 567, A132
- Wu P.-F., 2018, Monthly Notices of the Royal Astronomical Society, 473, 5468
- York D. G., et al., 2000, The Astronomical Journal, 120, 1579
- de Grijs R., Lee J. T., Clemencia Mora Herrera M., Fritze-v Alvensleben U., Anders P., 2003, New Astronomy, 8, 155