HI vs. H$\alpha$ - Comparing the Kinematic Tracers in Modeling the Initial Conditions of the Mice

1.8 MB·

by  S. Alireza Mortazavi

MNRAS000,1–12(2016) Preprint17January2017 CompiledusingMNRASLATEXstylefilev3.0
HI vs. Hα - Comparing the Kinematic Tracers in Modeling
the Initial Conditions of the Mice
S. Alireza Mortazavi (cid:63), 1 Jennifer M. Lotz, 2 Joshua E. Barnes, 3,4 George C. Privon, 5
Gregory F. Snyder 2
1Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles St., Baltimore, MD 21218, USA
2Space Telescope Science Institute, 3700 San Martin Dr., Baltimore, MD 21218, USA
3Institute of Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA
7 4Yukawa Institute for Theoretical Physics, Kyoto University, Kitashirakawa Oiwakecho, Sakyo-ku, Kyoto 606-8502, Japan
1 5Instituto de Astrof´ısica, Pontificia Universidad Cato´lica de Chile, Vicun˜a Mackenna 4860, 7820436 Macul, Santiago, Chile
0
2

n AcceptedXXX.ReceivedYYY;inoriginalformZZZ
a
J

3 ABSTRACT
1 Weexplorethe effectofusingdifferentkinematictracers (HIandHα)onreconstruct-
ing the encounter parameters of the Mice major galaxy merger (NGC 4676A/B). We
]
observed the Mice using the SparsePak Integral Field Unit (IFU) on the WYIN tele-
A
scope, and compared the Hα velocity map with VLA HI observations. The relatively
G
high spectral resolution of our data (R ≈ 5000) allows us to resolve more than one
. kinematic component in Hα emission line of some of the fibers. We separate the Hα
h
emission of the star forming regions from that of shock-heated regions based on the
p
- [N II]/Hα line ratio and the velocity dispersion of the components. We show that the
o kinematicsofstar-formingregionsagreewiththatofthecoldgas(HI),particularlyin
r thetidaltailsofthesystem.WeputconstraintstheencounterparametersoftheMice
t
s by matching the morphology and kinematics of these tidal tails with an equal-mass
a
galaxy merger simulation. We use an automated modeling method based on the Iden-
[
tikitsoftware package.In ourmethod,we quantify the goodness of the match and the

1 uncertaintyoftheresultingencounterparameters.Mostoftheinitialconditionsrecon-
v structed using Hα and HI kinematics are consistent within the estimated uncertainty,
9 and qualitatively agree with the results of previous works based on visual matching
1 techniques. For example, we find 210±50 Myrs, and 180±50 Myrs for the time since
8 40 40
pericenter, when modeling Hα and HI kinematics, respectively. This confirms that in
3
somecases,HαkinematicscanbeusedinsteadofHIkinematicsforreconstructingthe
0
initialconditionsofgalaxymergers,andourautomatedmodelingmethodisapplicable
.
1 to some merging systems.
0
7 Key words: galaxies:individual:NGC4676 – galaxies:interactions –galaxies:kine-
1 matics and dynamics
:
v
i
X
r
a 1 INTRODUCTION modeling puts new constraints on our understanding of
galaxy evolution and cosmology. For example, isolated hy-
Mergersarekeyprocessesingalaxyformationandevolution.
drodynamicaldisc-discgalaxymergershaveshownthatthe
They are one of the major contributors to the mass assem-
initial conditions such as pericentric separation and initial
bly of galaxies, they induce starbursts in galaxies, and they
orientationofdiscsaffectsthetimingandamountofmerger
are likely to be responsible for the transformation of disk-
inducedstarformation(Coxetal.2008,Snyderetal.2011).
dominated,rotation-supportedgalaxiestobulge-dominated,
Dynamical modeling constrains these parameters indepen-
dispersionsupportedones(Toomre&Toomre1972,Barnes
dentofthemeasuredstarformationhistoryinthesystem,so
& Hernquist 1996, Mihos & Hernquist 1996).
one can use them as independent tools for testing mod-
Reconstructing the encounter parameters of galaxy
els of merger-induced star-formation. In addition, we have
mergers (including their initial conditions, orbital param-
learned that the initial orientation of interacting discs cor-
eters, and observer-dependent parameters) via dynamical
relates with whether the remnant will be a fast or a slow
rotator.(Boisetal.2011;Naabetal.2014).Besides,wecan
(cid:63) [email protected] put constraints on cosmological dark matter simulations by
(cid:13)c 2016TheAuthors
2 Mortazavi et al.
measuring the orbital parameters (e.g. eccentricity) for a throughchaoticprocesses,anditishardtoreproducethem
statistical sample of galaxy mergers. Distribution of orbital notonlywithcollisionlesstestparticlesofIdentikit(Barnes
parametersofgalaxymergerscanbecomparedwiththatof &Hibbard2009;Mortazavietal.2016),butalsoinhydrody-
dark matter halo mergers in cosmological simulations (e.g. namicalsimulationsincludinggaseouscomponents.Itisalso
see Benson 2005, Khochfar & Burkert 2006). expensive to measure the kinematics of cold gas in tidally
Interactinggalaxiesoftenexperiencestrongtidalforces interacting galaxies with a spatial resolution that resolves
when they pass by each other, producing pronounced tidal the velocity gradient across the tails. The third option is to
features in discs. These features are usually strongest after measure the velocity of star forming regions using nebular
the first pericenter, though their shape and strength is a lineemission(e.g.Hα).Interactioninducestarformationin
complex function of the initial parameters of the orbit and gas-rich disc galaxies, and one often finds H II regions in
the structure of the galaxies (Dubinski et al. 1995, Springel the tidal tails and bridges (Jog & Solomon 1992; de Grijs
& White 1999, Barnes 2016). The sensitivity of the shape et al. 2003). Measuring Hα emission is a lot less expensive
and velocity of these features to the initial conditions make than measuring the velocity of stars or the cold gas (HI).
themstrongtoolsformodelingthedynamicsofgalaxymerg- Nonetheless,gasdissipationalsoaffecttheseregions.Inad-
ers and their initial conditions (Barnes & Hibbard 2009). dition, before using Hα as velocity tracer, we must make
Dynamical modeling of galaxy mergers is possible by sure that the ionized gas resides with the bulk of baryonic
finding the most similar simulation to the morphology and matter and is not displaced in position or velocity by non-
kinematics the data. Here,“most similar” is a vague term. gravitationalphenomenasuchashighvelocityshocksdriven
Mostpreviousattemptstomodelgalaxymergershaveused by supernovae (SNe) or active galactic nuclei (AGN).
qualitative, subjective matching criteria obtained by visual Here, we explore the effect of using different kinematic
inspection of the model and data (e.g. Toomre & Toomre tracers on the reconstructed encounter parameters. IFU
1972; Hibbard & Mihos 1995; Barnes & Hibbard 2009). In galaxy surveys like SAMI (Croom et al. 2012), CALIFA
Mortazavietal.(2016)wedevelopedanautomatedmethod (Sa´nchezetal.2012),andMaNGA(Bundyetal.2015)mea-
based on Identikit 2 (Barnes 2011). In this method we use suretheresolvedHαkinematicsofarelativelylargenumber
collisionless massless particles to reproduce tidal features. of galaxies, including major mergers. Consistency between
Our method is not only less subjective than the visual modeling Hα kinematics and modeling the more extended
matching techniques, but also provides well-defined error- HI emission shows that our method can be applied on the
bars for the initial parameters. data from these large surveys. As a result, we will have dy-
For modeling a galaxy merger, we need to know the namical models of statistically significant samples, and we
lineofsightvelocityoftidalfeatures.Therehavebeensome canestimatethedistributionoforbitalparametersofmajor
attempts to model merger systems without implementing galaxy mergers in the nearby Universe.
velocity information, only trying to match the morphology In this work, we focus on modeling one of the most fa-
of the model with data (Shamir et al. 2013, Holincheck mousgalaxymergersystems,NGC4676(a.k.a.Arp242,the
et al. 2016). Velocity measurement requires resolved spec- Mice).NGC4676isanearlystagegalaxymergeratredshift
troscopywhichtendstobemoreexpensiveandlessavailable z=0.02205. It has a very distinctive morphology consisting
than observing morphology. Optical morphology of many oftwostrongtidaltailsinthenorthandsouthofthesystem
galaxy mergers are/will be available through all sky imag- resemblingtwoplayingMice.Thestraighttailinthenorth-
ing surveys such as Sloan Digital Sky Survey (SDSS, York erngalaxyindicatesthatthisgalaxyisalmostedge-on.The
et al. 2000) and Large Synoptic Survey Telescope (LSST, southern disc seems to be slightly tilted, though still close
Abelletal.2009).However,thereismoredegeneracyinthe to an edge-on view. In this work, we model the Mice using
matched solutions when one does not utilize velocity infor- the kinematics of two different components: cold gas (HI),
mation(Hibbardetal.1994,Barnes2011),andmorphology- and star-forming regions (Hα). In §2 we describe reduction
only matching may result in a best-fit model that is incon- and analysis of the SparsePak IFU data for obtaining the
sistent with the observed velocity gradient across the tidal Hα kinematics of the Mice. In §3 we briefly mention the
tails and bridges (e.g. Borne & Richstone 1991,Hibbard & characteristicsoftheHIdata.In§4wedescribethemethod
Mihos 1995). we use for modeling the initial conditions of the Mice using
Different velocity tracers can be used to measure the bothHαandHIline-of-sightvelocitymaps,andwepresent
kinematics of tidal features. Velocities of stars are usu- the modeling results. In §5 we discuss these results, com-
ally ideal to match with test-particle and collisionless self- paring them to some previous reconstructions of encounter
consistent simulations. We may assume that stars have had parameters of the Mice in the literature and demonstrating
enough time to redistribute as collisionless particles, if the some of the their implications.
stellar population is formed long before the encounter be-
gins.Nonetheless,measuringthevelocityofstarsinthefaint
tidal tails and bridges requires a high signal to noise ratio
in the continuum of spectra which is expensive to obtain. 2 SPARSEPAK IFU DATA
Another option is to measure the velocity of cold neutral
2.1 Hα-[N II]Observations
hydrogengas(21cmHIemission).Neutralhydrogenisusu-
ally more extended than stars in disc galaxies, and produce We observed the Mice using the SparsePak Integral Field
stronger tidal features when discs interact in prograde or- Unit (IFU) on the WIYN telescope at Kitt Peak National
bits.Strongtidalfeaturesareusefulforconstrainingthedy- Observatory(KPNO)(Bershadyetal.2004)inMarch2008.
namicalmodel.However,oneshouldkeepinmindthatcold OurgoalwastomeasurethekinematicsofHαemissionline.
gas is dissipative. Some dissipative structures are produced Wedonotrequireauniformcoverageofvelocityinformation
MNRAS000,1–12(2016)
HI vs Hα 3
2.2 Data Reduction
Wefit1componentand2componenttripleGaussiancurves
30°41'26.425"
on top of a line with free slope as the background (contin-
uum) over the wavelength range 6665-6755 where Hα and
[N II] emission lines appear for the Mice system. The triple
Gaussian curves have free amplitudes, with centers sepa-
26.420" rated by the wavelength difference between [N II][λ6548],
Hα[λ6563], and [N II][λ6584]. The ratio of [N II][λ6584]/[N
II][λ6548] was fixed to the theoretical value of 2.95 (Acker
et al. 1989). The signal to noise of the continuum is not
high enough to fit the stellar model properly. However, in
26.415"
our analysis we take the underlying Hα absorption into ac-
count. We use an F-test to decide which fit to the emission
linesisfavoredbyourdata.Twocomponentfitsareusually
preferred in the central regions where the signal to noise
26.410" is higher and parts of the system with different kinematics
arelikelytooverlap.Figure4bshowsfibersinwhichdouble
component fits were favored. In each of these fibers, nar-
rowandbroadcomponentsaredeterminedbytheirvelocity
26.405" dispersion. In some of the very central fibers even a 2 com-
ponentfitisnotenoughtoproperlymodeltheshapeofthe
emission line. Figure 2 shows an example of a fiber where 2
componentmodelispreferred.Thisfiberismarkedwiththe
redcrossinFigure1.InFigure2bwecomparetheCALIFA
26.400"
spectrumoftheexactsamelocationontheMice.Themulti-
191°33'59.45" 59.46" plecomponentsareclearlywashedoutduetolowerspectral
resolution.
Figure 1. The HI surface density contours are plotted over the
SDSSr-bandimage.SkypositionsofthefibersinfourSparsePak
pointings are shown separately, with green circles, red triangles, 2.3 Shocked vs. Star Forming Regions
blue stars, and yellow squares. The red cross indicates the fiber
for which the spectrum is shown in Figure 2. The centers of the Hα emission originates from ionized gas regardless of how
twogalaxiesandthenortherntailmatchinr-bandimageandHI thegaswasionized.InnormalHIIregions,atomichydrogen
map.Thesoutherntail,however,seemtobegeneratingtwoself- is photo-ionized by the UV emission from O/B stars. Pho-
gravitation clouds of cold gas (indicated by the brown arrows). toionizationdoesnotsignificantlyaffecttheoverallkinemat-
These self-gravitating features can not be reproduced with test- icsofthegasrelativetooldstarsandtheneutralgasinthe
particlesinIdentikit.Weusethemorphologyofther-bandimage
vicinity. We expect the velocity obtained from Hα emission
alongwiththekinematicsoftheHIgas.
ofthesephoto-ionizedsourcestoberelativelysimilartothe
velocityofthebulkofthebaryons,whichismostlygoverned
by gravitational forces. In dynamical modeling of a galaxy
merger using collisionless test particles only gravitational
effects are to be considered. So, Hα kinematics of normal
HII regions should be the ideal for our modeling method.
over the system. Often, just the velocity of a few fibers in On the other hand, high speed stellar winds, SN remnants,
the tail regions is enough to break the degeneracy in the and feedback from Active Galactic Nuclei (AGN) produce
mergerparameterspace.SparsePakisespeciallysuitablefor shocks and heat up the interstellar gas. Shock-heated gas
this purpose as it has a relatively large field of view (∼1(cid:48)), alsoemitHα,butunlikephoto-ionizedgas,oftenitacquires
attheexpenseofmissingareasbetweenthesparselyplaced adifferentvelocitythroughanon-gravitationalprocess.The
fibers. We observed the Mice with 4 SparsePak pointings. kinematicsobtainedfromthisHαemissioncansignificantly
Figure1showsthelayoutofthefiberpositionsontheMice. disagree from that of the bulk of the baryons. Therefore, it
ForSparsePakobservationsweusedthebenchspectro- is important to distinguish emission from the photo-ionized
graphandthe860lines/mmgratingblazedat30.9◦ inorder and shocked regions.
2, obtaining a dispersion of 0.69 ˚A/pixel (FWHM) in the One way to separate the star-forming and shocked re-
wavelength range of 6050-7000˚A. We obtain a velocity res- gions is to utilize the Baldwin, Phillips & Terlevich (BPT)
olution of ∼ 31 km s−1. Our spectral coverage is less than diagram which uses [O III]λ5007 / Hβ, [N II]λ6583/Hα, [S
current and ongoing galaxy surveys such as CALIFA, and II]λλ6716, 6731/Hα, and [O I]λ6300/Hα flux ratios (Bald-
MaNGA, But our spectral resolution is higher. In the red win et al. 1981, Kewley et al. 2006). Our SparsePak obser-
band, the dispersion of CALIFA, and MaNGA surveys are vations,however,werelimitedtowavelengthrangeof6050-
2.0 ˚A/pixel and 0.83˚A/pixel, respectively. Higher spectral 7000˚AanddidnotincludeHβand[OIII]emissionlines.We
resolution enables us to resolve multiple emission line com- need a different method to distinguish shocks from normal
ponents,usuallyappearinginthecentralregionsofgalaxies star-forming regions.
where multiple gaseous components overlap. Theshockedregionstendtohavebothhigher[NII]/Hα
MNRAS000,1–12(2016)
4 Mortazavi et al.
0.2
1400 2
Compnent 1: Compnent 2:
1200 HvFNWeIαIlλ=oH6c3M5i8t=85y53=5.=8.65±37198±11696508.5000....66.238±82 ±77.. 278161k1.3.m2/s HvFNWeIαIlλ=oH6c1M5i6t=83y63=2.=2.67±15688±11471050..004..3..1123±90± 5511..43 8222k..29m/s 0.0
1000
ux(count) 860000 [NII]/Hα 00..42 3
Fl
400
0.6
1
200
0 0.8
6660 6670 6680 6690 6700 6710 6720 6730 6740 6750 0 50 100 150 200 250
λ( ) velocity dispersion σ (km/s)
(a)
Figure3.Plotof[NII]/Hαvs.velocitydispersionforallcompo-
nentswithS/N>3.Thedownwardarrowsindicatesthatthevalue
of[NII]/Hαisanupperlimit.WedidnottaketheunderlyingHα
1.4 absorptionintoaccount.(Seetext.)Datapointsarevisuallyclas-
2)csec− 1.2 HvFNWeIαIlλ=oH6c6M5i.t1=8y±3=1=1006...11404 75.43±.±40000..±11..110000..44 km/s sehimifigeihtdeteridn[tNofrIo3Im]g/rHonuαoprams.naGdlrvsoetualoprc-f1ioty(rlmodwiinsepgrelrresefigto)inoantrshe.atGnakrgoernuopuaps2c1o(.tmTopph)oenyheaanvrtees
ar1.0 likelytobeemittedfromshock-heatedgas.Group3(lowerright)
1−
have large velocity dispersion and lower [N II]/Hα than group
2−0.8 2 suggesting projection of multiple unresolved components. The
m
c vertical and horizontal dashed blue lines show respectively the
1s−0.6 limits of σ =90 km/s and [N II]/Hα=−0.2 for shocked regions
rg intheliterature.(Seetext)
e0.4
17−
0
(10.2
Fλ three groups. Group 1 has low [N II]/Hα and low velocity
0.0 dispersion, group 2 has high [N II]/Hα and high velocity
6660 6670 6680 6690 6700 6710 6720 6730 6740 6750
λ( ) dispersion,andgroup3haslow[NII]/Hαandhighvelocity
dispersion. We chose the dividing lines shown in Figure 3
(b)
toseparatethesethreegroups.InthisFigure,wealsoshow
Figure 2. (a) An example of the SparsePak spectrum in one of with dashed blue lines, σ = 90 km/s (limit of shocks in
thecentralfibers(indicatedwiththeredcrossinFigure1).The
Rich et al. 2011) and [N II]/Hα= −0.2 (limit of compos-
blackpointswitherrorbarsshowthespectrumdatapoints.The
itegalaxiesinstandardBPTdiagrams;Kewleyetal.2006).
solid blue curves show the fitted triple Gaussians (three compo-
In our spectral analysis we did not fit the stellar model to
nents)andthethefittedbackgroundline.Theredsolidlineshows
thespectra,andweusedalinewithfreeslopetomodelthe
the sum of fitted functions. One can visually confirm that more
thatoneoreventwocomponentsarenecessarytofullydescribein continuumaroundHα-[NII]triplet.Wedidnottaketheun-
theemissionlines.Inthiswork,wetestwhethertwocomponents derlyingHαabsorptionfeatureofthestellarcontinuuminto
are resolved in each fiber. (b) CALIFA spectrum from the exact account. This means that The Hα equivalent width (EW)
sameskyposition.Thepointsandsolidlinesarethesameas(a). is underestimated in our analysis, and [N II]/Hα is overes-
ThelowerspectralresolutionofCALIFAinstrumentwashesout timated. As a result, the data points in Figure 3 are upper
themultiplecomponents. limits.
Group 1 components are taken to be the components
originating from normal star forming regions. We use these
and usually exhibit higher velocity dispersion. Rich et al. components to make the Hα velocity map for dynamical
(2011) showed that the histogram of velocity dispersion of modeling. (See Figure 4a). Most of the components in this
emission line components of their high spectral resolution group are both below and to the left of the horizontal and
IFU data of some galaxy mergers reveals a bi-modality. vertical dashed blue lines, respectively. The rest, are within
Theyconcludethatcomponentswithσ<90km/sareemit- the error-bar distance to these lines. The smooth rotation
ted from HII and turbulent star forming regions, and com- seenacrosscirclesinFigure4a,andtheiragreementwiththe
ponents with σ > 90 km/s are emitted from low velocity velocityofcoldgasconfirmsthechoiceofthesecomponents
shocks.CombiningthisconclusionwiththetraditionalBPT as the velocity of star forming regions.
diagrams, we use a plot of [N II]/Hα vs. the velocity dis- Group2componentsareemittingfromtheshockedre-
persion to separate the shocked regions from normal star gions. The spatial position of these components are shown
forming ones. inFigure4c.Theylieonthetwosides,aboveandbelowthe
Figure3showstheplotof[NII]/Hαvs.velocitydisper- discsnearthecoresofbothgalaxies,whichisconsistentwith
sionforallkinematiccomponentswithS/N>3.Fromvisual the bi-cones of shocked material in Wild et al. (2014). This
inspection, we find that the data points are clustered into suggeststhattheprocessproducingtheshock-ionizedgasis
MNRAS000,1–12(2016)
HI vs Hα 5
happening in the core of the galaxies, which can be merger 4 DYNAMICAL MODELING
inducedcentralstar-burstorAGN.Wildetal.(2014)argue
FordynamicalmodelingoftheMice,weusetheautomated
thatbasedonthespatialdistributionof[OIII]/Hαratioin
pipeline that we developed based on Identikit (Barnes &
CALIFA data and the X-ray observations, NGC 4676A do
Hibbard 2009; Barnes 2011; Mortazavi et al. 2016). Iden-
nothostanAGNandsourceofionizationisfastoutflowing
tikitisasoftwarepackageformodelingtheinitialcondition
shocks, probably from a central starburst. In NGC 4676B
of major galaxy mergers. It uses a precompiled library of
they do not rule out the possibility of an AGN being the
N body simulations of galaxy mergers with different initial
source of bi-conal structure. The agreement of the spatial
conditions. The discs are modeled by test particles. This
distributionofregionswithhigh[NII]/Hαbetweenourdata
facilitates the run of multiple discs at the same time, and
and CALIFA analysis, for which stellar model is taken into
improvesthepaceofsearchintheparameterspace.Interac-
account, suggests that the underlying Hα absorption does
tivevisualmatchingwithIdentikitwasusedbyPrivonetal.
not significantly change our measurements of [N II]/Hα.
(2013) to model the Mice along with 3 other systems.
Figure4dshowsthespatialdistributionofgroup3com-
ponents. They are in the central regions and, their high ve- Theisolatedgalaxiesconsistofcollisionlessmassiveand
locitydispersioncanbeexplainedbyprojectionofmorethan massless test particles. massive particles are distributed in
2 components that are not resolved. spherically symmetric fashion representing the mass of the
Usingthesegroups,wecanestimatethefractionoftotal dark matter halo, the disk, and the bulge. In these models,
Hα flux emitted from shocked gas. To do so, we divide the the luminous mass fraction, fL =(Md+Mb)/(Mh+Md+
sumoftheHαfluxfromgroup2componentstothetotalHα Mb)is0.2andtheconcentrationparameterofthehalomass
flux. Group 3 components can be both unresolved group 1 profile, ch is 4. The discs are represented by test particles
and2components.Weuse[NII]/Hα=−0.2limittoseparate that are initially in circular orbits. The scale length of the
them for the purpose of this estimate. We obtain a fraction disc , ad, is 1/3 of the scale radius of the dark matter halo,
of26±2%ofHαfluxbeingemittedfromshock-heatedgas. ah. Barnes (2016) explored the effect of different structure
The error of this fraction is obtained from the variance of parameters in the isolated galaxies on the morphology of
shock fraction in 100 bootstrapped samples of data points tidalfeatures.Inthisworkwekeepworkingonasinglemass
in [N II]/Hαvs. σ space (i.e. Figure 3). model for the sake of simplicity.
Inthisworkweuseanautomatedmatchingfeaturethat
is introduced in Identikit II (Barnes 2011). In Identikit II,
thesimilaritybetweenmodelanddataisquantifiedwithan
informal measure called score. The user places small boxes
3 JVLA HI DATA in the same position as the tidal tails and bridges of the
system, extended in two spatial directions and one velocity
Thekinematicsofcoldgasareavailablefromradiointerfer-
direction (LOSV). Identikit calculates the scores based on
ometricobservationsofHI21cmemissionlineperformedby
thenumberofdisctestparticlesresidinginthesephasespace
Hibbard&vanGorkom(1996).Theobservationstookplace
boxes. The models that better reproduce the tidal features
in May of 1991 and May of 1992 using the D and C con-
inthedatawillplacemoreparticlesintheseboxes,andthey
figurations of then Very Large Array (VLA), respectively.
willgainabetterscore.IdentikitIIhasbeenappliedwitha
The HI data has a spatial resolution of 20”(≈ 9 kpc) and
limitedrangeofparametersontheMicegalaxiesbyBarnes
velocityresolutionof43.1km/s.ThetotalHImassisinthe
&Privon(2013).Theyinspectedthemodelswithhighscore
systemis7.5×109M (Hibbard&vanGorkom1996;Wild
(cid:12) visually and confirmed that they show good, or fair visual
et al. 2014).
match with the HI data.
In the Mice, there are a couple of self-
InMortazavietal.(2016)wedescribedthemethodthat
gravitating/dissipative features in the southern tail
we have developed for automated searching of parameter
with no counterparts in the optical r-band image. They
spaceandreconstructingtheinitialparameterswithrobust
are indicated by the brown arrows in Figure 1. After a
uncertainties.Inthismethodwefirstputtheboxesontidal
few unsuccessful trials to reproduce these features with
features of the galaxy merger using a semi-automated tech-
Identikit, we decided to use the morphology of the optical
nique. The boxes are placed randomly in regions that are
image(SDSS,r-band)alongwiththevelocityoftheHIgas,
far enough from the center of galaxies. Boxes near the cen-
same as Privon et al. (2013).
ter of the galaxies, where test particles have higher number
density, are always populated, and they do not put strong
constraint on the initial conditions of the merger. So, we
put the boxes over the ends of the tidal tails to enforce the
3.1 Kinematics of Hα vs. H I
similarityintheoverallshapeandthevelocitygradient.We
WeusethekinematicsoftheHαemissionfromnormalstar- ignore self-gravitating features like stellar clusters or blobs
formingregionsformodelingthedynamicsoftheMice.Fig- of gas, because test particles in Identikit do not reproduce
ure 4a shows the Hα velocity map the mice galaxies over- them. In cases where no velocity information is available
plotted on the H I velocity map. In this color plots the (e.g. the absence of Hα in the southern tail of the Mice),
striking agreement between the two velocity maps is evi- we put boxes that are only extended in spatial directions,
dent. This agreement is especially important in the tidal anddonotputconstraintonvelocity.Thesizeoftheboxes
tailswhichplaythemostimportantroleindynamicalmod- are selected to match the spatial resolution of the obser-
eling, as we try to find the model that best matches the vation. In case of the HI data The boxes are 8 kpc wide,
morphology and kinematics of the tidal features. almostmatchingthespatialresolutionofJVLAobservation
MNRAS000,1–12(2016)
6 Mortazavi et al.
400 400
320 320
240 240
160 160
s) s)
80 m/ 80 m/
K K
0 y ( 0 y (
cit cit
o o
80 el 80 el
v v
160 160
240 240
320 320
400 400
(a) (b)
0.2 250
225
0.1
200
0.0 s)
175m/
k
0.1 150n (
o
0.2NII]/Hα 125spersi
[ di
0.3 100city
o
75 el
0.4 v
50
0.5
25
0.6 0
(c) (d)
Figure 4.(a)Velocityofcomponentsingroup1ofFigure3,whichhavelowvelocitydispersionandlow[NII]/Hα,isover-plottedon
theHIVelocitymap.Thereisanexcellentagreementintheoverallvelocitygradientandthevelocityoftidaltails.Theblueboxshows
thezoom-inareaoftheotherpanels.(b)ThefibersinwhichtwocomponentfitispreferredbytheF-testoveronecomponentfit.The
velocityofthetwocomponentsareshownbythecolormap.Theupperhalfofthecirclesshowthevelocityofthenarrowcomponents,
and their lower half show the velocity of the broad components. (c) The fibers with components in group 2 of Figure 3. The spatial
distribution of this fibers is consistent with the bi-cone shocked structure seen in Wild et al. (2014), and suggests a central process is
responsibleforshocks.(d)Thefiberswithcomponentsingroup3ofFigure3.Thesecomponentshavelargevelocitydispersionandlow
[N II]/Hα. The color of points show their velocity dispersion. They are in the central regions where multiple gaseous component are
projectedinthesamelineofsight.NotethatthebroadcomponentinthespectrumofFigure2belongstothisgroup.
of Hibbard & van Gorkom (1996). When modeling the Hα sincefirstpass,velocityscalingandvelocityoffset.Forfuture
kinematics,thesizeoftheboxeswheresettobeclosetothe references we call these five parameters ”external parame-
sizeofSparsePakfibers(4kpc).Inplaceswheretheerrorin ters”. Calculating the score for each of these models, Iden-
Hα velocity measurement was larger that the extent of the tikit constrains viewing angle, initial orientation of discs,
box in velocity direction, we increased the size of the box length scaling and position offset. We call these parameters
proportionally.Figure6showsasampleofselectedboxesin “internalparameters”.Bycalculatingthescoreforalllibrary
morphology and velocity directions. members,weobtaina5dimensionalscoremapforexternal
parameters. Each point of this score map corresponds to a
Table1showstherangeofencounterparametersweex-
set of calculated internal parameters. The model with the
ploredinthiswork.WeproducedalibraryofIdentikitmod-
maximumscoreisthebest-fitmodelandthecorresponding
elswithavarietyofeccentricity,pericentricseparation,time
MNRAS000,1–12(2016)
HI vs Hα 7
ParameterClass Parameter RangeTested
orbitalparameters eccentricity [0.80-1.10]
pericentricdistance [0.03125-1.0000]×Rvir
massratio 1
observerdependent timesincepericenter fromfirstpasstosecondpass
parameters viewingangle viewedfrom320evenlydistributed
directionsonasphere
position setbylockingthecenters
lengthscalingL setbyviewingangleandlockingthecenters
velocityoffset [-0.1,0.1]
velocityscalingV [-0.500-+0.500]∗
initialorientation 1280evenlydistributedorientations
ofdiscs
Table 1. Range of encounter parameters explored. The total number of parameter space points for which we calculated score is ∼
3.2×1010.∗Thevelocity(length)scaling,V (L),relatesthedimensionlessvelocity(length)inIdentikittothephysicalvelocity(length).
external and internal parameters are the best-fit encounter suggests that the reconstructed eccentricity is less than or
parameters. equal to this value. The best-fit pericentric separation is
The uncertainty of the best-fit parameters are calcu- 17±12 kpc for Hα modeling, and is 24±8 kpc for HI mod-
6 13
lated using a Bootstrap statistical method. Score, by itself, eling.Themergerstagedefinedasthetimesincepericenter
doesnotprovidestatisticalprobabilityasothermeasuresof divided by the time from first passage to second passage,
goodness of the fit like χ2. So, in order to estimate the un- is 0.19±0.20 and 0.14±0.70 for Hα and HI modelings, re-
0.01 0.01
certaintyofthescoreswere-dotherandomboxpositioning spectively.Thebest-fittimesincepericenteris200±40 Myrs
30
procedure described above several times (10 times in this for Hα and 170±90 Myrs for HI. The best viewing angles
50
work), and calculate the score each time. This process ran- in spherical coordinates (θ, φ) relative to the orbital plane
domly moves the position of boxes on the tidal tails and is (76◦±3◦,0◦±61◦) and (67◦±4◦,52◦±9◦ ) for Hα and HI
9◦ 41◦ 9◦ 66◦
bridges of the system. The variation in the scores for each models, respectively. For each Identikit model, we viewed
model measures the error of the score, which eventually, the system from 320 equally separated viewing angles, so
translates into the uncertainty of the best-fit parameters. the resolution of our search in viewing angle was about
The models with scores that are within 1, 2, or 3 standard 11◦. The initial orientation (i,ω) of NGC 4676A is found
deviation from the best score, respectively, determine the to be (21◦±11◦,29◦±33◦) and (56◦±31◦,50◦±55◦) in mod-
5◦ 38◦ 24◦ 96◦
extent of 1σ, 2σ, or 3σ separation in the parameter space eling Hα and HI kinematics, respectively. For initial ori-
from the best-fit model. entation of NGC 4676B we find (38◦±21◦,187◦±20◦) and
10◦ 24◦
In Mortazavi et al. (2016) we applied our modeling (69◦±23◦,162◦±25◦) for Hα and HI modeling.
27◦ 61◦
methodononedge-onviewsofhydrodynamicalsimulations
AllofthereconstructedparametersformodelingHIand
of prograde disc-disc interactions (similar to the Mice). We
Hαkinematicsareconsistentwithin1σuncertainty.Thiscan
evaluated the bias in the reconstructed interaction parame-
beseenbycomparingFigure5(a)and(b)inwhichthedark
ters comparing them to the known input parameters of the
areas with high scores mostly overlap. Figure 5 also shows
hydrodynamical simulations. We showed that when testing
the model by Privon et al. (2013) which was obtained with
young stars in the simulation (representing the Hα data),
visual matching of HI kinematics of the Mice to Identikit
theaveragebiasisnotsignificantineccentricity,R ,and
peri modelsandconfirmingthemodelwithafullyself-consistent
merger stage (< 1σ). When testing the cold gas of the
N-bodysimulation.Privonetal.(2013)onlyexploredmod-
same simulation, the reconstructed physical R was al-
peri elswithparabolicorbits(eccentricity=1).Excepteccentric-
most twice the correct value, but the reconstruction of ec-
ity, our parameters for modeling Hα kinematics are consis-
centricity and time was only biased by less that 1σ.
tent to within 1σ of the Privon et al. (2013) parameters.
TheparametersofourHImodelareabitfurtheraway,but
stillwithin2σ ofPrivonetal.(2013).Figure6showsoneof
4.1 Summary of Results
the best-fit models to the Hαkinematics. This figure can be
In Figure 5 one can see 3 slices of the 5 dimensional score compared to Figure 5 of Privon et al. (2013).
mapatthebest-fitparameterpoint,acrosseccentricity,peri-
centricseparation,andmergerstage.Intheseplotsthebest-
fitmodelcanbecomparedwiththebest-fitmodelbyPrivon
5 DISCUSSION
et al. (2013).
The best-fit parameters for Hα and HI kinematics are The peculiar morphology of the Mice has made it a curi-
presentedinTable2.Thebest-fiteccentricityis≤0.80±0.05 ous case for dynamical modeling. There has been several
−−−
and ≤ 0.80±0.15 in modeling Hα and HI kinematics, re- attemptstomodeltheinitialconditionsofthissystem.The
−−−
spectively. The reason for the ≤ sign is that the range of reported parameters in some of these attempts are shown
eccentricitiesexploredinthisworkwasfrom0.8to1.1.Ob- in Table 3. In all but one of these dynamical models the
tainingthebesteccentricityontheedgeofoursearchrange mass ratio is presumed to be 1. This assumption is made
MNRAS000,1–12(2016)
8 Mortazavi et al.
fractional physical e mergerstage physicaltime (θ,φ) NGC4676A NGC4676B
Rperi Rperi ∆T/Tfirstto (Myr) degrees◦ (i1,ω1) (i2,ω2)
(Rvir) (kpc) secondpass∗ degrees◦ degrees◦
Hα 0.38±0.16 17±12 ≤0.80±0.05 0.19±0.21 200±40 (76±3,0±61)∗ (21±11,29±33) (38±21,187±20)
0.09 6 −−− 0.01 30 9 41 5 38 10 24
HI 0.50±0.03 24±8 ≤0.80±0.15 0.14±0.70 170±90 (67±4,52±9 ) (56±31,50±55) (69±23,162±25)
0.44 13 −−− 0.01 50 9 66 24 96 27 61
Table 2. Reconstructed parameters from modeling Hα and HI kinematics. ∗ The resolution of automated search in viewing angle was
∼11◦ (viewdepthparameter=3inIdentikit;seeBarnes2011).Theerrorsreportedherearefromthewidthofdistributionofviewing
anglesinmodelswithhighscore.Notethatnegativeazimuthalanglesshouldbeinterpretedasa360◦ rotation.Theseparametersmatch
thedefinitionsof(Privonetal.2013)andcanbecomparedwiththeirvalues.
1.2 1.2 1.15
1.0 1.0 1.10
1.05
0.8 0.8
y
()RRperivir 00..46 ()RRperivir 00..46 ccentricit001...990050
E
0.2 0.2
0.85
0.0 0.0 0.80
0.2 0.2 0.75
0.750.800.850.900.951.001.051.101.15 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Eccentricity Tafterpericenter/∆Tfirsttosecondpass Tafterpericenter/∆Tfirsttosecondpass
5.0 σ 2.5 σ best fit
(a)
1.2 1.2 1.15
1.0 1.0 1.10
1.05
0.8 0.8
y
()RRperivir 00..46 ()RRperivir 00..46 ccentricit001...990050
E
0.2 0.2
0.85
0.0 0.0 0.80
0.2 0.2 0.75
0.750.800.850.900.951.001.051.101.15 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Eccentricity Tafterpericenter/∆Tfirsttosecondpass Tafterpericenter/∆Tfirsttosecondpass
5.0 σ 2.5 σ best fit
(b)
Figure 5. (a) The average score map for modeling Hα kinematics of the Mice galaxy merger. Score map is 5-dimensional parameter
space.Hereweareshowing3slicesofthescoremap,throughthebest-fitparameter,acrosseccentricity,pericentricseparation,andtime
sincepericenter.Thedarkareasindicatebetterscore.Thegraycolorisscaledwiththestandarddeviationofscoresinthebest-fitmodel.
shownwiththecyanbox.Theredboxshowsthebest-fitparametersofPrivonetal.(2013)whichusedvisualmatchingtoHIkinematics.
Inthemiddlepanelitishiddenbehindthecyanbox.Themagentaboxshowsthebest-fitparametersoftheHImodel.(b)Theaverage
scoremapformodelingHIkinematicsoftheMice.Coloredsquaresarethesame.
because of the similar stellar mass ratio (1.5:1, Wild et al. modelingarebothconsistentwiththemedianofthesemod-
2014) and the strong tidal tails in both galaxies (Toomre els. In both of our models the Mice is found to be at a rel-
& Toomre 1972; Barnes 2016). We also adopt equal mass atively early stage. The galaxies are currently at about 1/5
merger models in this work. (1/8) of time from the first passage to the second passage
fortheHα(HI)model.Takingtheuncertaintyintoaccount,
Time is the best-constrained encounter parameter for bothvaluesareconsistentwithPrivonetal.(2013)(≈1/5).
the Mice galaxies so far. Among the reported values in Ta- OurHαandHImodelspredictthatthesecondpassagewill
ble 3, the lowest value for time since pericenter is 120 Myr occurin530±480 Myrsand530±370 Myrs,respectively,also
220 220
for Toomre & Toomre (1972). The largest is 430 ±190 Myr consistent with model of Privon et al. (2013) for which the
for Holincheck et al. (2016), and the median is 175 Myr for second passage will happen in ≈ 430 Myrs.
Privon et al. (2013). Our reconstructions from HI and Hα
MNRAS000,1–12(2016)
HI vs Hα 9
Y 10 kpc Y 50 km/s
10 kpc 10 kpc
X V
(a) (b)
V 10 kpc Z 10 kpc
50 km/s 10 kpc
X X
(c) (d)
Figure 6. (a) The reconstructed Identikit model (green and red test particles) are over-plotted on the SDSS r-band image of the
Micegalaxies.Thegreen,redandblueboxesindicatethephasespaceregionsassociatedwithNGC4676A,NGC4676B,andbothdisks,
respectively.TheyareusedbyIdentikittocalculatescoreasameasureofthegoodnessofthematch.Theyellowcirclesshowtheposition
of SparsePak fibers used to measure the Hα kinematics. (b) and (c) show the the line of sight kinematics of the reconstructed model
versus the y and x spatial directions, respectively. The position of yellow rectangles show the Hα velocity measured at each fiber, and
theirextentindicatethedispersionintheHαemissionline.Inthesepanels,theboxesshowthevelocityofphasespaceregionsdetermined
by the measured Hα velocities on tidal tails. Some boxes in panel (a) are placed where there is no Hα velocity information (no fiber),
therefore they do not have counterparts in panels (b) and (c), and are only used to constrain the morphology of the model. (d) shows
themodelviewedfromthedirectionoftheyaxisinpanel(a).ThisfigurecanbecomparedtoFigure5ofPrivonetal.(2013).
Bothofourmodelspreferellipticalorbits,whichiscon- Physical pericentric distance in our Hα and HI models
sistentwithToomre&Toomre(1972)andHolinchecketal. are 17±12 kpc and 24±8 kpc which are consistent with all
6 13
(2016).Nonetheless,inMortazavietal.(2016)itwasshown previous reported values except that of Barnes (2004) (8.9
that our modeling method underestimates the eccentricity kpc).R ishardtoconstrain,asitseffectonmorphology
peri
by ∼0.1 in most tests on hydrodynamical merger simula- of tides is correlated with mass profile of isolated galaxies,
tions with parabolic orbits. The large error-bars in the re- i.e. disc scale-length and halo scale-radius, as well as the
constructed eccentricity in Mortazavi et al. (2016) indicate
that our modeling method is not sensitive enough to ec-
centricity, and other constraints may be required to better
constrain the initial energy of orbits in galactic encounters.
MNRAS000,1–12(2016)
10 Mortazavi et al.
dynamicalmodel µ e Rperi physicaltimesince kinematics
(kpc) pericenter(Myrs)
Toomre&Toomre(1972)1 1.0 0.6 ... 120 Hαonlycenters
Mihosetal.(1993)2 1.0 0.6 23 180 Hαlongslit
Barnes(2004)3 1.0 1.0 8.9 170 HImap∗
Privonetal.(2013)2 1.0 1.0 14.8 175 HImap∗
Holinchecketal.(2016)1 0.59±0.52 0.70±0.27 15.9±5.0 430±190 novelocityused
Table 3. Previous dynamical models of the Mice. ∗ From Hibbard & van Gorkom (1996). 1 Test particles with central potential. 2
Collisionlessselfconsistentsimulations.3 Hydrosimulationsincludinggasandstarformation.
luminous mass fraction1 (Barnes 2016). In Mortazavi et al. the Hubble time. Otherwise, one may ask why the galax-
(2016) we showed that our method overestimates R in ies did not merge during the previous pericenter. Assuming
peri
mostofthetestsonhydrodynamicalsimulations,suggesting that the mass of the dark matter halo of the larger partner
thatourreconstructionofR oftheMiceisanupperlimit. is 1011−1012M the period of the orbit will be 10-3 Gyrs.
peri (cid:12)
Most recently, Holincheck et al. (2016) used judge- This is less than the age of the Universe and suggests that
ment of citizen scientists (Galaxy Zoo, Lintott et al. 2010) the orbit has become bound due to recent processes that
to match the morphology of a relatively large sample of happenedinthelast3-10Gyrs,e.g.recentmassgrowthdue
merger systems. In addition to parameters explored in this to minor mergers.
work, they varied mass ratio parameter, and they managed Merger induced star formation rate (SFR) is affected
to present best-fit models for a relatively large sample of bytheencounterparametersinhydrodynamicalsimulations.
62 tidally interacting galaxy pairs. For each galaxy pair Cox et al. (2008) show that mergers with larger pericentric
∼ 4.6×106 parameter space points were rejected initially distanceinducestarburstsatlatertimes.Moreover,therel-
for not displaying tidal features, and ∼ 50,000 parameter ative inclination of the discs with respect to orbital plane
spacepointswereviewedandjudgedvisuallybycitizensci- correlateswiththetotalamountofmergerinducedstarfor-
entists.However,asseeninTable3,theuncertaintyoftheir mations(alsoseeSnyderetal.2011).Ontheotherhand,the
measurementsarerelativelylargecomparedtotheonespre- starformationinthesesimulationsdependsonpresumptions
sented in this work. Lack of kinematic information in their about sub-grid physics, the physical processes that happen
analysis and fewer parameter space points they explored inscalessmallerthantheresolutionofthesimulation.Even
contributetolargererrorbars.Inthiswork,wecalculatethe the most recent simulations of galaxy formation do not yet
scorefor3.2×1010 pointsinparameterspaceelaboratedin resolve particles at scales smaller than stars. To implement
Table 1. star formation, most current simulations assume a density
Thestrongtidaltailsandtheedge-onviewoftheMice threshold above which stars are formed. The exact value of
makesitastraight-forwardcasefortestingagalaxymerger thisthresholdisdifferentindifferentsimulationsandvaries
model. However, one may ask how common the encounter with the resolution of the hydrodynamical simulation (e.g.
parameters of the Mice are among galaxy mergers. While seeBrooks&Teyssier2015),andtheyareusuallycalibrated
the answer to this question requires modeling a larger sam- to reproduce the Kennicutt-Schmidt relationship between
pleofmergers,wecancompareourvaluesforreconstructing gasmasssurfacedensityandstarformationrate(Kennicutt
theMice,withthedistributionoforbitalparametersofdark 1998). Barnes (2004) proposed a shock-induced star forma-
matterhalomergersincosmologicalsimulations.Theeccen- tion recipe in order to reproduce extended merger induced
tricity of the Mice found in this work (e ≤ 0.80), is lower star formation observed in tidal tails of interacting galaxies
than the most frequent eccentricity (e = 1.0) among dark (e.g. see de Grijs et al. 2003). In addition to sub-grid star
matterhalomergers(e.g.seeFigure11ofBenson2005and formation recipe, the model for feedback from SNe, stel-
Figure 6 of Khochfar & Burkert 2006). Parabolic (e = 1) lar winds, and AGN adds to the complications of star for-
orbits have zero total mechanical energy, and are expected mation model in hydrodynamical simulations (e.g. see Cox
to be seen if the protogalaxies are initially at rest and far et al. 2006; Oppenheimer et al. 2010). Finding independent
fromeachother.Theuncommonlyloweccentricitywefound constraints on star formation history in galaxy mergers im-
maybeduetoabiasinourmodelingmethod.Ontheother proves our understanding of the sub-grid physics in hydro-
hand, it may indicate that orbits of interacting galaxies de- dynamical simulations.
cayearlierthanthedarkmatterhalostheylivein.Theperi-
Our estimate of the merger stage and other encounter
centric separation found for the Mice (0.38±0.16 R and
0.09 vir parameters of the Mice is independent of star formation
0.50±0.03 R for Hα and HI, respectively) are both un-
0.44 vir rate. It only depends on the dynamics of tidal features. As
certain. They both fall into the most common values in the
a result, by looking at a“Mice-like”simulation of a galaxy
distributionofR sindarkmatterhalomergersoftheMil-
peri mergerwiththesameencounterparametersandatthesame
lenniumcosmologicalsimulation(seeFigure4ofKhochfar&
timestepwefoundinthiswork,wecanputanindependent
Burkert2006).Asanitycheckforthisansweristocalculate
constraint on the sub-grid physics affecting merger induced
the period for this elliptical orbit to see if it is longer than
starformation.Todoso,wecancomparethetotalenhance-
mentofstarformationinthehydrodynamicalsimulationto
themeasuredstarformationintheMice.Unfortunately,the
1 LuminousmassfractionfL= MdM+dM+bM+bMh.InIdentikitmodels measurements of the SFR in the Mice is uncertain and de-
usedinthisworkfL=0.2. pends on the observed star formation indicator. Wild et al.
MNRAS000,1–12(2016)