Three Dusty Star-forming Galaxies at z ∼ 1.5: Mergers and Disks on the Main Sequence

The majority of star-forming galaxies form a correlation between their star formation rates (SFR) and stellar masses (M${}_{\star }$; e.g., Rodighiero et al. 2011; Whitaker et al. 2012; Sargent et al. 2014; Schreiber et al. 2015). Often referred to as the main sequence of star-forming galaxies (MS), the correlation has a tight scatter of ≈0.3 dex in M${}_{\star }$ (e.g., Brinchmann et al. 2004; Daddi et al. 2007; Elbaz et al. 2007; Noeske et al. 2007; Speagle et al. 2014; Tacchella et al. 2016). The other, populated by starburst galaxies (SB), lies at SFRs a few times higher than the MS at fixed M${}_{\star }$, and does not have a tight scatter (e.g., Rodighiero et al. 2011; Sargent et al. 2014). In addition to these two star-forming populations, there is a population of quiescent galaxies that lie below the main sequence (e.g., Noeske et al. 2007; Leslie et al. 2016; Tacchella et al. 2016). The small scatter in the MS is typically interpreted to be the result of smooth star formation driven by the net gas inflow and outflow rate and the gas consumption rate (e.g., Bouché et al. 2010; Dutton et al. 2010; Tacchella et al. 2016; Scoville et al. 2017), while galaxies with enhanced SFRs are typically observed to be undergoing major mergers or short-lived starburst events (e.g., Mihos & Hernquist 1996; Di Matteo et al. 2008; Bournaud et al. 2011; Rodighiero et al. 2011; Whitaker et al. 2012; Silverman et al. 2018). Indeed, the strongest SB galaxies in the local universe are nearly all observed to be undergoing major mergers (e.g., Joseph & Wright 1985; Armus et al. 1987; Sanders & Mirabel 1996).

Dusty star-forming galaxies are a submillimeter-identified class of galaxy that are selected for their bright dust emission and characteristically elevated SFRs. They likely peak in number density around z ∼ 2 (Casey et al. 2014) and represent a crucial phase in the evolution of z = 0 giant elliptical galaxies (e.g., Swinbank et al. 2004, 2006; Engel et al. 2010; Michałowski et al. 2010; Menéndez-Delmestre et al. 2013; Toft et al. 2014). The physical cause of their bright submillimeter emission is still a matter of debate (see review by Casey et al. 2014). Those with the highest SFRs, of order a few times 10³ M${}_{\odot }$ yr⁻¹, are virtually all driven by major mergers (e.g., Swinbank et al. 2004; Greve et al. 2005; Alaghband-Zadeh et al. 2012). However, many studies conclude they are not unanimously merger driven, especially at more modest SFRs of ∼10² M${}_{\odot }$ yr⁻¹ (e.g., Genel et al. 2008; Tacconi et al. 2008; Bothwell et al. 2010, 2013; Swinbank et al. 2010, 2011; Casey et al. 2011; Hodge et al. 2012; Drew et al. 2018; McAlpine et al. 2019). The selection of DSFGs is typically based on a flux density cutoff in the submillimeter, canonically S_ν ≳ 2–5 mJy, or SFR ≳ 100 M${}_{\star }$ yr⁻¹ for galaxies at z ≳ 1 (e.g., Smail et al. 1997; Barger et al. 1998; Hughes et al. 1998; Eales et al. 1999). This corresponds to a selection for high SFRs. Some studies find DSFGs lie on the high-mass end of the main sequence at z > 2 (Michałowski et al. 2012, 2014) and suggest that this implies they are driven by gas accretion rather than merging (Michałowski et al. 2017). Given their selection, for a given bin of SFR it is stellar mass that determines their location with respect to the main sequence.

With increasing redshift, the normalization of the MS increases while the typical scatter remains the same (e.g., Daddi et al. 2007; Elbaz et al. 2007; Rodighiero et al. 2010; Speagle et al. 2014; Whitaker et al. 2014; Scoville et al. 2017). There is tension in the literature over whether mergers outside the local universe lie above or within the MS because the SFRs of SB galaxies in the local universe are comparable with those on the MS at moderate redshifts. Some studies find mergers predominantly lie in the SB regime above the MS at z > 1 (e.g., Kartaltepe et al. 2012; Hung et al. 2013; Cibinel et al. 2019), or predominantly occupy higher SFRs, irrespective of the stellar mass (e.g., Kartaltepe et al. 2010; Ellison et al. 2013). Other studies suggest that the stellar masses typically measured for the elevated SB population at z ∼ 2 are unreliable, and that most SB galaxies actually comprise the high-M${}_{\star }$ end of the MS (e.g., Michałowski et al. 2012, 2014; Koprowski et al. 2014).

The determination of galaxy classification is vital to addressing the question of what role mergers play in galaxy evolution through cosmic time. Numerous techniques have been employed to determine galaxy classification including imaging, kinematics, and nonparametric analyses (see review by Conselice 2014). Imaging studies classify galaxies based on visual signatures of mergers, interactions, or disks in images (e.g., Hubble 1926; de Vaucouleurs 1963; Abraham et al. 1996a, 1996b; Brinchmann et al. 1998; Kartaltepe et al. 2015). While large imaging studies make this kind of analysis relatively easy to perform, they may fall short in a few key ways at $z\gtrsim 1$ when trying to distinguish mergers from disks. At these redshifts, the optical waveband begins to probe rest-frame UV emission originating from younger stars. Not only is this light more highly dust obscured than at longer rest-frame wavelengths, but the UV is not a sensitive probe of the bulk of the stellar mass, which most closely resembles the distribution of the total mass of the galaxy. Additionally, kinematically regular galaxies may appear morphologically disturbed in imaging (e.g., Bournaud et al. 2008), because galaxies at z > 1 tend to be clumpier than their low-z counterparts (e.g., Abraham et al. 1996a; Elmegreen & Elmegreen 2006). This may lead to disks being misclassified as mergers or irregular galaxies. Also, surface brightness dimming may hide features of a merger, such as tidal tails (e.g., Hibbard & Vacca 1997).

Kinematic studies overcome some of the issues associated with imaging studies. They measure the motions of gas inside galaxies via observations of the Doppler shift of emission lines. Kinematic observations of disk galaxies exhibit smooth rotational fields, while mergers show more complex velocity fields (see Glazebrook 2013 for a review). Kinematic studies are still susceptible to misclassification however, especially if observed near or shortly after coalescence or if observed at low spatial resolution (e.g., Hung et al. 2015).

Nonparametric analyses (e.g., Abraham et al. 2003; Lotz et al. 2004, 2008) rely on the grouping of galaxies in parameter space to identify galaxy morphology and assembly history. The strength of this technique is that it can be applied to any type of galaxy without prior knowledge about the form the model should take. Its weakness is that the distinction between classes may not be as clear as with imaging or kinematic studies (Conselice 2014).

In this paper we present rest-frame UV and rest-frame optical imaging, a rest-frame optical emission line kinematic analysis, and a nonparametric analysis of three DSFGs at z ∼ 1.5 as case studies of high SFR-selected galaxies. Our goal is to identify their physical drivers. We then compare their location in the SFR–M${}_{\star }$ plane with their physical drivers to investigate whether merging DSFGs lie above, within, or below the MS.

Section 2 of this paper describes our observations and data reduction, Section 3 presents our imaging and kinematic analyses, Section 4 describes Gini–M₂₀ statistics, Section 5 discusses the galaxies in the context of the MS, and Section 6 summarizes. Throughout this work we adopt a Planck Λ CDM cosmology with ${H}_{0}=67.7$ km s⁻¹ Mpc⁻¹, Ω_Λ = 0.6911 (Planck Collaboration et al. 2016) and a Chabrier initial mass function (Chabrier 2003).

The spectroscopic data presented in this manuscript were obtained with the Multi-Object Spectrometer For Infra-Red Exploration (MOSFIRE; McLean et al. 2010, 2012) on Keck I as part of a spectroscopic follow-up campaign to measure the redshifts of DSFGs identified at flux densities >12.4 mJy and >2.4 mJy at 450 μm and 850 μm, respectively, with Scuba-2 in the COSMOS field (see Casey et al. 2013, 2017). The galaxies in the present paper were selected from their parent sample based on their high signal to noise (S/N) and spatially resolved Hα and [N ii] emission. The three galaxies presented in this paper are the only ones that have sufficiently resolved emission to allow for a kinematic analysis from the original 114 targets in the Casey et al. (2017) sample. One galaxy, named 850.95, is also published in Drew et al. (2018, hereafter D18) as an observational counter example to the hypothesis that galaxies at intermediate redshifts may have declining rotation curves. The rotation curve of this galaxy is flat in the outer galaxy, much like typical disk galaxy rotation curves at z = 0. In this paper, we discuss the Hα kinematics of 850.95 and refer the reader to D18 for a discussion of its dark matter content.

The three galaxies presented in this paper, 450.25, 450.27, and 850.95, have prefixes of either 450 or 850, corresponding to the wavelength at which they were initially identified in Casey et al. (2013). Only 450.27 is detected at both 450 and 850 μm. Table 1 lists basic characteristics of each galaxy. For additional details about the parent sample and its selection, see Casey et al. (2017).

Table 1.  Parameters of Each Galaxy

Source	450.25	450.27	850.95
R.A.	10:00:28.58	09:59:42.92	09:59:59.80
Decl.	+02:19:28.3	+02:21:45.1	+02:27:07.4
zspec	1.515	1.531	1.555
M⋆ (M${}_{\odot }$)	(3.4 ± 0.5) × 1011	(3.0 ± 0.6) × 1011	(3.8 ± 3.0) × 1010
LIR (${L}_{\odot }$)	(1.5${}_{-0.5}^{+0.7}$) × 1012	(4.1 ± 0.5) × 1012	(3.0${}_{-0.9}^{+1.2}$) × 1012
SFR (M${}_{\odot }$ yr−1)	157${}_{-53}^{+80}$	382${}_{-45}^{+51}$	373${}_{-90}^{+110}$
Gini	0.53 ± 0.01	0.63 ± 0.01	0.48 ± 0.01
M20	−1.725	−1.808	−1.654
Physical Driver	Merger	Merger	Disk

Note. Positions are from Casey et al. (2017). M${}_{\star }$ is estimated using MAGPHYS (da Cunha et al. 2008) with the HIGHZ extension (da Cunha et al. 2015). The errors on M${}_{\star }$ differ slightly from Casey et al. (2017) because they are estimated using a range of SFHs from a continuous star formation history to an instantaneous burst, following the procedure of Hainline et al. (2011). The errors on L${}_{\mathrm{IR}}$ and, as a direct result, SFR are estimated following the procedure of Casey (2012).

Download table as:  ASCIITypeset image

H-band spectroscopic observations of 450.25, 450.27, and 850.95 were obtained on 2013 December 31 at W. M. Keck Observatory. The FWHM of seeing was 085. Galaxy 450.25 was observed for a total integration time of 2880 s, 450.27 for 1320 s, and 850.95 for 1920 s. The slit width was set to 07 and a 15 ABBA nod pattern was used between exposures. The spectra were reduced with the MOSPY Data Reduction Pipeline,⁵ and one-dimensional spectra were extracted using the iraf⁶ package, apall. Apertures of extraction were placed on each pixel with an aperture radius of half the average seeing. Adjusting the aperture size does not significantly change the radial velocity and velocity dispersion measurements presented in the following sections. Variance weighting was used in the spectral extraction.

Figure 1 shows S/N spectra for 450.25 and 450.27, including Hα, [N ii], and continuum emission. See Figure 1 in Drew et al. (2018) for the S/N spectrum of 850.95. The white crosses denote the seeing in the vertical dimension and the combined instrument resolution and seeing in the horizontal dimension. The negative images to the north and south of the galaxy are characteristic of the data coaddition step associated with the nod pattern. Slits were randomly oriented with respect to the galaxies because MOSFIRE slit masks prevent custom orientations for individual galaxies when observing in multiplex mode.

Zoom In Zoom Out Reset image size

We simultaneously fit Gaussians to Hα, [N ii]λ6548 and [N ii]λ6583 to the extracted apertures using the iraf package, splot, forcing the centroids to have fixed spacing and the Gaussian widths to be tied. We exclude pixels contaminated with telluric emission in the fits. Errors in fit centroids and widths are derived using 1000 Monte Carlo perturbations of the data by sky noise. The bottom panels of Figures 2 and 3 show the position–velocity and position–dispersion diagrams measured from the MOSFIRE spectra. These will be discussed further in Section 3. We find no evidence of broad line emission in any of the spectra and Casey et al. (2017) finds no evidence of x-ray emission or mid-IR spectral energy distribution (SED) slopes in these galaxies that would indicate luminous active galactic nuclei are present.

Zoom In Zoom Out Reset image size

The imaging includes H-band from the UltraVISTA survey (rest-frame ∼6400 ${\rm{\mathring{\rm A} }}$; McCracken et al. 2012) and Hubble F814W (rest-frame ∼3200 ${\rm{\mathring{\rm A} }}$; Koekemoer et al. 2007) data. While Y, J, H, and Ks imaging is available for these galaxies from the UltraVISTA survey, we present the H band because it matches the band the spectra are observed in. The H-band imaging is consistent with the longest available from the UltraVISTA survey, the Ks band. These observations are the longest-wavelength high spatial resolution imaging available. The F814W imaging is the longest-wavelength Hubble Space Telescope (HST) imaging available. The images are presented in the top rows of Figures 2 and 3 and will be discussed further in Section 3.

Zoom In Zoom Out Reset image size

We classify our galaxies based on their position–velocity and position–dispersion diagrams into two categories: mergers or disks. In the position–velocity diagram we expect a disk to have smooth, symmetric velocity field about a single spatial axis. In the position–dispersion diagram we expect a peak centered on the spatial symmetry axis of the position–velocity diagram. Evidence for a merger would be a disrupted position–velocity and position–dispersion fields and/or a discontinuity in velocity between two galaxies (e.g., Glazebrook 2013 and references therein).

3.1. 450.25

3.2. 450.27

3.3. 850.95

Next, we compute the nonparametric diagnostics Gini and M₂₀ (e.g., Abraham et al. 2003; Lotz et al. 2004, 2008) for each of our galaxies. The Gini statistic quantifies how the light is distributed throughout a galaxy. A Gini value of 1 would imply all the emission originates from a single pixel, while a Gini value of 0 would imply a uniform light distribution across multiple pixels. M₂₀ measures the second moment of the galaxy's brightest 20% of pixels relative to the total second moment. These two quantities have been shown to roughly separate mergers, ellipticals, and disk galaxies (e.g., Lotz et al. 2008).

To perform the analysis, first we isolate the pixels associated with each galaxy above a threshold S/N of 8 in H-band imaging using the Python package, Photutils (Bradley et al. 2019). Thresholds are chosen to be as low as possible while still separating unrelated field galaxies in the segmentation maps. Next, we run the Photutils source deblending routine on galaxy 450.27 to separate the two perpendicular components seen in the H-band imaging (see Figure 3). Whether or not this deblending is performed does not change the Gini–M₂₀ classification of merger for this galaxy. Finally we use the Python package, Statmorph (Rodriguez-Gomez et al. 2019) to measure Gini and M₂₀ statistics from the segmentation maps. Statmorph also simultaneously fits Sérsic profiles to the pixels associated with each galaxy in the segmentation maps, the results of which we discussed in previous sections.

Figure 4 shows our galaxies in Gini–M₂₀ space. Galaxy 450.27 lies within the merger region, while galaxies 450.25 and 850.95 lie within the disk region. Our kinematic and morphological analyses conclude 450.25 is a merger, 850.95 is a disk, and 450.27 is possibly a merger, although it is a bit ambiguous. Considering this together with the location of 450.27 in the Gini–M₂₀ diagram, we conclude that it is indeed a merger. The Gini–M₂₀ plot in Figure 7 from Lotz et al. (2008) shows that the overwhelming majority of galaxies that lie in the merger region of Gini–M₂₀ space are true mergers. This figure also shows that mergers may lie in any region of Gini–M₂₀ space. Galaxy 450.25 is confirmed via imaging to be undergoing a merger or interaction (see Figure 2), but it lies within the disk region of the Gini–M₂₀ diagram. This may be caused by the fact that it is in an early merger/interaction stage so the H-band emission is not yet concentrated in a just few pixels.

Zoom In Zoom Out Reset image size

Now we consider the merger/disk classifications in the context of the main sequence of star-forming galaxies. The SFRs and stellar masses of our sample are reported by Casey et al. (2017), who measure these quantities using the high-z extension of MAGPHYS (da Cunha et al. 2008, 2015) using multiwavelength photometry (UV through submillimeter) from the COSMOS collaboration (Capak et al. 2007; Laigle et al. 2016). In the present paper, to account for systematic uncertainties on stellar mass caused by the assumption of a star formation history, we follow the procedure of Hainline et al. (2011), which was developed to estimate uncertainties on M${}_{\star }$ for similarly selected DSFGs at z ∼ 2. To summarize, we take the errors on M${}_{\star }$ to be half the difference between the stellar masses estimated using instantaneous burst histories and those estimated using continuous star formation histories. Works by Michałowski et al. (2012, 2014) demonstrate that different assumed SFHs can strongly affect the derived M${}_{\star }$. This systematic uncertainty from choice of SFH is larger than those reported by the MAGPHYS fits.

Figure 5 shows our sample plotted on the main sequence of star-forming galaxies at z ∼ 1.5 from three works in the literature. The first (Rodighiero et al. 2011) is fitted to far- and near-IR selected galaxies between 1.5 < z < 2.5 in the COSMOS and GOODS-South fields using combined UV and IR SFRs. The second (Whitaker et al. 2014) is fitted to galaxies from the 3D-HST photometric catalogs between 1.5 < z < 2.0 using combined UV and IR SFRs. The third (Koprowski et al. 2016) is fitted to 850 μm selected galaxies at z > 1.5 from the S2CLS survey (Geach et al. 2017) using UV through submillimeter SEDs.

Zoom In Zoom Out Reset image size

Galaxy 450.27 is the only galaxy to lie within the 0.3 dex scatter of the main sequence. Galaxy 850.95 is the only disk galaxy in our sample but it lies at a SFR ∼4× greater than the MS. Galaxy 450.25 is a merger with visible tidal tails but it lies at a SFR ∼4× less than the MS. The overall distribution of the galaxies is consistent with the works of Michałowski et al. (2012, 2014, 2017) and Koprowski et al. (2016), who find that DSFGs roughly comprise the high M${}_{\star }$ end of the MS. In particular, all three galaxies fit within the scatter of the samples at 1.0 < z < 1.5 and 1.5 < z < 2.0 in Figure 10 of Michałowski et al. (2017).

Michałowski et al. (2017) suggest because DSFGs comprise the high M${}_{\star }$ end of the MS major mergers are not a dominant driver of their SFRs. Mergers are short-lived events that are expected to elevate the SFRs of galaxies above the main sequence. The present paper, as well as many other works (e.g., Di Matteo et al. 2008; Hung et al. 2013; Cibinel et al. 2019), demonstrate that merging systems may lie on the main sequence. Puglisi et al. (2019) find that up to 50% of the most massive galaxies on the MS at $z\sim 1.3$ may have star formation driven by merging. Major merger activity may not enhance SFRs at every stage of a merger (e.g., Bergvall et al. 2003; Di Matteo et al. 2008; Jogee et al. 2009; Narayanan et al. 2015; Fensch et al. 2017; Silva et al. 2018). Additionally, high resolution hydrodynamical simulations by Fensch et al. (2017), show that mergers at high redshift have lower star formation efficiencies compared with those at low redshift. We caution the reader that a detailed analysis of the merger classifications of DSFGs on the high M${}_{\star }$ end of the MS needs to be performed before it can be concluded that they are not undergoing merging.

It is interesting that the SFR of 450.25, a clear early stage merger, is ≲4× the MS value at its M${}_{\star }$. Quiescence is known to correlate with galaxy compactness (e.g., Bell et al. 2012; Lee et al. 2018), but the location of 450.25 in Gini–M₂₀ space demonstrates the galaxy is not compact. It is possible that 450.25 was previously quenched but the merger has yet to fully turn star formation on again. A study of the molecular gas content of this galaxy would be illuminative.

Mergers or galaxy interactions drive the SFRs of galaxies lying at higher SFRs than the main sequence at z = 0 (e.g., Mihos & Hernquist 1996), however at higher redshift where the SFRs of galaxies at all M${}_{\star }$ are elevated, the distinction between these two populations is less clear. Kinematic observations of DSFGs in the literature reveal a mix between merging and secularly evolving galaxies (e.g., Swinbank et al. 2006; Alaghband-Zadeh et al. 2012) but such studies are limited to small sample sizes because they are observationally expensive. The determination of where DSFGs lie in the SFR–M${}_{\star }$ plane is important to help determine their importance in galaxy evolution through the cosmos. Some studies find that DSFGs sit above the main sequence (e.g., Hainline et al. 2011), but recent work by Michałowski et al. (2012, 2014, 2017) has determined the average DSFG actually comprises the high stellar mass end of the MS, and conclude that this is evidence that major mergers do not drive their SFRs.

In this paper we combine imaging, kinematic, and nonparametric analyses to determine whether three submillimeter identified DSFGs at z ∼ 1.5 have merger or secular recent SHFs. We find two to be undergoing merging or interactions and one to be an isolated disk galaxy.

Rest-frame UV and optical imaging of galaxy 450.25 shows tidal tails, which are clear evidence it is undergoing a merger or interaction. The rest-frame optical emission is not well fit by a Sérsic profile because the best-fit parameter values are unphysical. This is due to warps caused by the interaction/merger. Position–Velocity and position–dispersion diagrams reveal disorder characteristic of a merger. The Gini coefficient and M₂₀ values, G = 0.534 ± 0.006 and M₂₀ = −1.725, place 450.25 in the disk region of the Gini–M₂₀ diagram, possibly because this is an early stage merger or interaction.

Rest-frame UV imaging of galaxy 450.27 shows two distinct star-forming knots which we conclude are likely two galaxy cores close to coalescence. Rest-frame optical imaging shows one distinct emission region rather than two knots, although it is observed through seeing with a FWHM comparable with the separation between the two UV knots. The rest-frame optical emission is not well fit by a Sérsic profile because the best-fit parameter values are unphysical. The position–velocity diagram looks Keplerian while the position–dispersion diagram looks asymmetric in the magnitude of velocity dispersion on each side of the galaxy. The Gini coefficient and M₂₀ values, $G={0.633}_{-0.005}^{+0.007}$ and M₂₀ = −1.808, place 450.27 in the merger region of the Gini–M₂₀ diagram.

Rest-frame optical imaging of galaxy 850.95, published by D18, is well fit by an exponential disk profile with Sérsic index n = 1.29 ± 0.03. The rest-frame UV emission is clumpy and is offset from dust continuum emission detected with ALMA. The position–velocity and position–dispersion diagrams show clear signatures of rotation with a velocity profile well fit by an arctangent function and a centrally peaked, symmetric dispersion curve. The Gini coefficient and M₂₀ values, G = 0.481 ± 0.005 and M₂₀ = −1.654, place 850.95 in the disk region of the Gini–M₂₀ diagram.

Despite its disk classification, 850.95 sits at a SFR ∼4× above the MS, placing it in the starburst regime. Despite merging activity, 450.27 lies within the scatter of the MS, and 450.25 lies ∼4× below the MS. It is unexpected that a merging galaxy has a lower SFR than typical galaxies on the main sequence. Further investigation as to the cause of the suppressed star formation is needed.

Our sample hints that perhaps the specific SFR (SFR/M${}_{\star }$) is not a useful discriminator of the recent merger history of high-M${}_{\star }$ galaxies at z > 1. A detailed investigation of a statistical sample of DSFGs is needed to determine the recent star formation histories of DSFGs both on and off the main sequence to determine the physics driving their SFRs.

The authors thank Justin Spilker and Jorge Zavala for useful discussions, as well as the anonymous reviewer who provided valuable comments and suggestions. P.D. acknowledges financial support by a NASA Keck PI Data Awards: 2012B-U039M, 2011B-H251M, 2012B-N114M, 2017A-N136M, NSF grants AST-1714528 and AST-1814034, and the University of Texas at Austin College of Natural Sciences. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. Some of the data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. This research made use of APLpy, an open-source plotting package for Python hosted at http://aplpy.github.com, Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013), and the python packages Matplotlib (Hunter 2007), Numpy (van der Walt et al. 2011), and Pandas (McKinney 2010).