TY - JOUR
AU1 - Rodriguez-Gomez,, Vicente
AU2 - Snyder, Gregory, F
AU3 - Lotz, Jennifer, M
AU4 - Nelson,, Dylan
AU5 - Pillepich,, Annalisa
AU6 - Springel,, Volker
AU7 - Genel,, Shy
AU8 - Weinberger,, Rainer
AU9 - Tacchella,, Sandro
AU1 - Pakmor,, Rüdiger
AU1 - Torrey,, Paul
AU1 - Marinacci,, Federico
AU1 - Vogelsberger,, Mark
AU1 - Hernquist,, Lars
AU1 - Thilker, David, A
AB - Abstract We have generated synthetic images of ∼27 000 galaxies from the IllustrisTNG and the original Illustris hydrodynamic cosmological simulations, designed to match Pan-STARRS observations of log10(M*/M⊙) ≈ 9.8–11.3 galaxies at |$z$| ≈ 0.05. Most of our synthetic images were created with the skirt radiative transfer code, including the effects of dust attenuation and scattering, and performing the radiative transfer directly on the Voronoi mesh used by the simulations themselves. We have analysed both our synthetic and real Pan-STARRS images with the newly developed statmorph code, which calculates non-parametric morphological diagnostics – including the Gini–M20 and concentration–asymmetry–smoothness statistics – and performs 2D Sérsic fits. Overall, we find that the optical morphologies of IllustrisTNG galaxies are in good agreement with observations, and represent a substantial improvement compared to the original Illustris simulation. In particular, the locus of the Gini–M20 diagram is consistent with that inferred from observations, while the median trends with stellar mass of all the morphological, size and shape parameters considered in this work lie within the ∼1σ scatter of the observational trends. However, the IllustrisTNG model has some difficulty with more stringent tests, such as producing a strong morphology–colour relation. This results in a somewhat higher fraction of red discs and blue spheroids compared to observations. Similarly, the morphology–size relation is problematic: while observations show that discs tend to be larger than spheroids at a fixed stellar mass, such a trend is not present in IllustrisTNG. methods: numerical, techniques: image processing, galaxies: formation, galaxies: statistics, galaxies: structure 1 INTRODUCTION The past few years have seen substantial progress in our understanding of galaxy formation and evolution from a statistical perspective, both on the observational side with surveys such as the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (Grogin et al. 2011) and the Sloan Digital Sky Survey (SDSS, York et al. 2000), and on the theoretical side with large-scale hydrodynamic cosmological simulations such as Illustris (Vogelsberger et al. 2014a,b; Genel et al. 2014), EAGLE (Schaye et al. 2015; Crain et al. 2015), and Horizon-AGN (Dubois et al. 2014). However, galaxy formation is complex and fully understanding it will require increased synergy between observations and theory. One way of achieving such increased interaction between theory and observations is by means of so-called ‘forward modelling’ of simulation data. In particular, the generation and subsequent analysis of synthetic images from hydrodynamic simulations is quickly becoming a powerful tool to connect galaxy formation theory with observations (e.g. Lotz et al. 2008a; Jonsson, Groves & Cox 2010; Snyder et al. 2015a,b; Bignone et al. 2017; Bottrell et al. 2017a; Trayford et al. 2017). Nevertheless, the generation of realistic synthetic images is not without challenges of its own. While some aspects of synthetic images are relatively well understood and standardized by now, such as the use of stellar population synthesis models (e.g. Bruzual & Charlot 2003), other aspects have been a consistent challenge for theorists, such as the detailed modelling of dust absorption and scattering. Some simple models that estimate the overall, spatially unresolved effects of dust on starlight have been found to work well in a variety of situations (e.g. Charlot & Fall 2000). However, quantifying the spatially resolved effects of dust on galaxy structure is particularly challenging because this requires modelling a very large ensemble of photons as it travels through the interstellar medium (ISM), considering wavelength-dependent absorption and scattering events along the way. This complex problem is usually known as radiative transfer in astrophysics, and is closely related to the concept of ray tracing in other fields, such as computer science. Thankfully, there are many publicly available codes designed for radiative transfer calculations in astrophysics, such as sunrise (Jonsson 2006), hyperion (Robitaille 2011), and skirt (Baes et al. 2011; Camps & Baes 2015). In this paper, we create realistic synthetic images from the state-of-the-art IllustrisTNG simulation (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018) using the radiative transfer code skirt (Baes et al. 2011; Camps & Baes 2015), designing them to match observations of low-redshift galaxies from the Pan-STARRS |$3{\pi}$| Steradian Survey (Chambers et al. 2016). Then, we quantify various structural parameters of both the simulated and observed galaxies using the same morphology code, which allows us to make a fair comparison between theory and observations. This constitutes the first systematic study of image-based galaxy morphology from a hydrodynamic cosmological simulation while including the effects of (spatially resolved) dust attenuation and scattering. This paper is organized as follows. In Section 2, we introduce the simulations used for this work and define the observational and simulated galaxy samples. In Section 3, we explain how synthetic images are created using the radiative transfer code skirt. In Section 4, we describe the newly developed statmorph code, which we use to measure various image-based morphological diagnostics. Our main results are presented in Section 5. Finally, we discuss these results and present our conclusions in Section 6. 2 METHODOLOGY 2.1 The IllustrisTNG simulation suite The IllustrisTNG Project (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018b; Springel et al. 2018) is a suite of magneto-hydrodynamic cosmological simulations run with the moving-mesh code arepo (Springel 2010; Pakmor, Bauer & Springel 2011; Pakmor et al. 2016), featuring an updated version of the Illustris galaxy formation model (Vogelsberger et al. 2013; Torrey et al. 2014) described in Weinberger et al. (2017) and Pillepich et al. (2018a). In this study, we use the highest resolution version of ‘TNG100’ (hereafter referred to as the IllustrisTNG simulation), which follows the evolution of 2 × 18203 resolution elements within a periodic cube measuring 75h−1 ≈ 110.7 Mpc per side, which translates to an average mass of the baryonic resolution elements of |$1.39 \times 10^6 \, {\rm M}_{\odot }$|⁠. The gravitational softening length of the DM and stellar particles, which essentially sets the spatial scale of the simulation, is 0.5h−1 ≈ 0.74 kpc. This value is ideal for a comparison to the Pan-STARRS observations described in Section 2.2. The IllustrisTNG model was roughly tuned to match the following observables: (i) the global star formation rate (SFR) density at |$z$| = 0–8, (ii) the galaxy mass function at |$z$| = 0, (iii) the stellar-to-halo mass relation at |$z$| = 0, (iv) the black hole-to-stellar mass relation at |$z$| = 0, (v) the halo gas fraction at |$z$| = 0, and (vi) galaxy sizes at |$z$| = 0. Note, however, that the model was not tuned to match galaxy morphology, much less using image-based statistics such as the ones presented in this paper. We also use data from the original Illustris simulation (Genel et al. 2014; Vogelsberger et al. 2014a,b; Nelson et al. 2015; Sijacki et al. 2015), which follows the evolution of 2 × 18203 resolution elements within a periodic cube measuring 75h−1 ≈ 106.5 Mpc per side, with a baryonic mass resolution of |$1.26 \times 10^6 \, {\rm M}_{\odot }$| and a gravitational softening length of 0.5h−1 ≈ 0.71 kpc. The main updates to the original Illustris galaxy formation model (Vogelsberger et al. 2013; Torrey et al. 2014) implemented in IllustrisTNG, besides the inclusion of ideal magneto-hydrodynamics, consist of a new active galactic nucleus (AGN) feedback model that operates at low accretion rates (Weinberger et al. 2017) and various modifications to the galactic winds, stellar evolution, and chemical enrichment schemes (Pillepich et al. 2018a). In both the IllustrisTNG and original Illustris simulations, galaxies were identified with the subfind halo-finding algorithm (Springel et al. 2001b; Dolag et al. 2009). Unless noted otherwise, the cosmological parameters used throughout this work are Ωm = 0.3089, Ωb = 0.0486, ΩΛ = 0.6911, σ8 = 0.8159, ns = 0.9667, and h = 0.6774, consistent with recent Planck measurements (Planck Collaboration XIII 2016).1 2.2 The observational galaxy sample We consider galaxies observed with the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS), specifically the first completed component, Pan-STARRS1 (PS1), which consists of a 1.8 m telescope equipped with a 1.4 gigapixel camera, producing images with a pixel scale of 0.258 arcsec pixel−1. The PS1 instrument has been used in a number of distinct sky surveys, in particular the |$3{\pi}$| Steradian Survey, which mapped the Northern hemisphere using five photometric filters (PS1 g, r, i, |$z$|⁠, andy) at an angular resolution of ∼1 arcsec and with a 5σ point source limiting sensitivity of up to ∼23 mag (Chambers et al. 2016). The final, stacked image products have a pixel scale of 0.25 arcsec pixel−1, slightly smaller than the native pixel scale, as a result of the warp pixel processing stage (Waters et al. 2016). A catalogue of galaxies from the |$3{\pi}$| Survey can be found in the PS1 Optical Galaxy Survey (POGS), an ongoing project that aims to create a value-added, multiwavelength atlas of galaxies in the nearby Universe (Vinsen & Thilker 2013; Thilker et al., in preparation). For each galaxy included in the POGS catalogue, we use its i-band stacked image (along with its associated segmentation map, stellar mask, weight map, and point spread function, PSF) for the morphological analysis presented in Section 5. We obtain the stellar masses and redshifts of POGS galaxies using the NASA Sloan Atlas (NSA),2 keeping only galaxies that are successfully matched: those with angular separations between POGS and NSA coordinates smaller than 2.5 arcsec. From an initial sample of 62 061 POGS galaxies, our coordinate matching criterion results in a sample of 56 614 objects with NSA counterparts, which is presented in Fig. 1. This figure shows that, by construction, POGS galaxies have stellar masses |$M_{\ast } \gt 10^{9.5} h^{-2} \, {\rm M}_{\odot }$| and redshifts |$z$| < 0.05. Figure 1. View largeDownload slide The distribution of stellar masses and redshifts of 56 614 galaxies from the POGS sample. The red rectangle indicates the redshift and stellar mass range considered in our study. The colour scale indicates the number of galaxies per 2D bin. Figure 1. View largeDownload slide The distribution of stellar masses and redshifts of 56 614 galaxies from the POGS sample. The red rectangle indicates the redshift and stellar mass range considered in our study. The colour scale indicates the number of galaxies per 2D bin. In order to make a straightforward comparison to a single snapshot from the IllustrisTNG and original Illustris simulations, we only consider galaxies in the redshift interval |$z$| = 0.045–0.05. Within this redshift range, we find that the completeness of the POGS sample (relative to the NSA) is high, but falls below 90 per cent at |$M_{\ast } \gtrsim 10^{11} h^{-2} \, {\rm M}_{\odot }$|⁠. Therefore, we only consider POGS galaxies with stellar masses in the range M* = 109.5–|$10^{11} h^{-2} \, {\rm M}_{\odot }$|⁠, or M* ≈ 109.8–|$10^{11.3} \, {\rm M}_{\odot }$|⁠. We also note that the spatial resolution of PS1 observations for this redshift range is ∼1 kpc, which is ideal for a comparison with the simulations described in Section 2.1. Our final POGS sample, which consists of 10 829 galaxies, is indicated with a red rectangle in Fig. 1. 2.3 The simulated galaxy sample On the simulation side, we consider a single simulation snapshot at |$z$| = 0.0485 (snapshot 95 in IllustrisTNG, or snapshot 131 in Illustris original), consistent with the observational redshift range mentioned in Section 2.2. Within this snapshot, we consider all simulated galaxies with |$M_{\ast } \gt 10^{9.5} \, {\rm M}_{\odot }$|⁠, a stellar mass range that encompasses that of the POGS sample. Such galaxies are typically resolved with more than ∼2250 (∼2500) stellar particles in the IllustrisTNG (Illustris original) simulation. The resulting mock galaxy catalogue, which includes both central and satellite galaxies, consists of 12 470 (14 721) galaxies in IllustrisTNG (Illustris original). 3 GENERATING SYNTHETIC IMAGES In order to generate reasonably realistic synthetic images from our simulations, we not only predict the light distribution produced by stellar populations and star-forming regions, but also model the effects of dust attenuation and scattering as the starlight traverses the ISM. This requires radiative transfer calculations, which can be very computationally expensive for large images (especially above a thousand pixels on each dimension), and not very worthwhile for galaxies with negligible amounts of star-forming gas, which we use as a proxy for dust (see Section 3.2.2). Therefore, in order to reduce the overall computational cost of our synthetic images, we have developed two synthetic image generation pipelines, with and without dust radiative transfer, which we refer to as our ‘skirt pipeline’ and ‘galaxev pipeline’, respectively. As their names suggest, our pipelines use the skirt radiative transfer code (Baes et al. 2011; Camps & Baes 2015) and the galaxev stellar population synthesis code (Bruzual & Charlot 2003). We process galaxies with negligible fractions of star-forming gas using our inexpensive galaxev pipeline, while the rest are treated with full radiative transfer using our skirt pipeline. The threshold is defined at a star-forming gas fraction fgas,sf = 1/100 (with respect to the total baryonic mass), as illustrated in Fig. 2. This is approximately the value at which we find that dust effects begin to have a noticeable impact on galaxy morphology. Figure 2. View largeDownload slide The fraction of star-forming gas as a function of stellar mass for all galaxies in snapshot 95 from IllustrisTNG, which corresponds to |$z$| = 0.0485. The 2D histogram shows the overall galaxy distribution, while the black line shows the median trend. Synthetic images of galaxies above the magenta line are created with our skirt pipeline (Section 3.2), while galaxies below the line are processed with our galaxev pipeline (Section 3.3). Galaxies with zero star-forming gas are shown at fgas,sf ≈ 10−6. Figure 2. View largeDownload slide The fraction of star-forming gas as a function of stellar mass for all galaxies in snapshot 95 from IllustrisTNG, which corresponds to |$z$| = 0.0485. The 2D histogram shows the overall galaxy distribution, while the black line shows the median trend. Synthetic images of galaxies above the magenta line are created with our skirt pipeline (Section 3.2), while galaxies below the line are processed with our galaxev pipeline (Section 3.3). Galaxies with zero star-forming gas are shown at fgas,sf ≈ 10−6. Throughout this section, we express some distances in comoving coordinates (which is how distances are handled internally by the arepo code), which we indicate with a ‘c’ prefix. In particular, we refer to comoving kpc as ‘ckpc’. 3.1 General details Several aspects of the synthetic images are the same for both our skirt and galaxev pipelines. Each galaxy is observed from a viewing angle perpendicular to the xy-plane of the simulation volume, which essentially corresponds to a random orientation for each galaxy. The field of view of each image is equal to 15 times the (3D) stellar half-mass radius of the corresponding galaxy in IllustrisTNG, and 10 times in Illustris original. These values were chosen empirically, so that the final image products are large enough for the morphological measurements described in Section 4, after adding an adequate amount of background noise (Section 3.4.2). The number of pixels is chosen so that the resulting pixel scale matches that of the observations (0.25 arcsec pixel−1, which corresponds to |$0.174 \, h^{-1}$| ckpc pixel−1 at |$z$| = 0.0485). For each stellar particle, we define an adaptive smoothing scale equal to the (3D) distance to the Nth nearest stellar particle. This approach was described in Torrey et al. (2015), although we use N = 32 instead of N = 16. We find that our morphological measurements are not particularly sensitive to the number of neighbours used, at least for choices ranging from N = 4 to N = 64. Unless noted otherwise, all smoothing calculations in this section are carried out through a convolution with a ‘standard’ smoothed particle hydrodynamics (SPH) spline kernel (Monaghan & Lattanzio 1985; Hernquist & Katz 1989; Monaghan 1992) redefined over the interval [0, hsml] (Springel, Yoshida & White 2001a), where hsml is some appropriately chosen smoothing length: \begin{eqnarray*} W(r, h_{\rm sml}) = \frac{2^{\nu } \sigma }{h_{\rm sml}^{\nu }} \left\lbrace \begin{array}{lcc}1 - 6 q^2 + 6 q^3 & \quad {\rm if} & \quad 0 \le q \le \displaystyle \frac{1}{2}, \\ 2(1-q)^3 & \quad {\rm if} & \quad \displaystyle \frac{1}{2} \lt q \le 1, \\ 0 & \quad {\rm if} & \quad q \gt 1, \end{array}\right. \end{eqnarray*} (1) where q ≡ r/hsml, r is the radial distance, ν is the number of dimensions, and σ is a normalization constant with the values \begin{equation*} \frac{2}{3}, \frac{10}{7{\pi} }, \frac{1}{{\pi} } \end{equation*} in one, two, and three dimensions, respectively. 3.2 skirt pipeline In this section, we describe our skirt synthetic image generation pipeline, based on the radiative transfer code skirt (Baes et al. 2011; Camps & Baes 2015), which we use to produce synthetic images including the effects of dust attenuation and scattering. 3.2.1 Stellar sources Each stellar particle in the simulation represents a coeval stellar population. The spectral energy distributions (SEDs) of stellar particles older than 10 Myr are modelled with the galaxev population synthesis code (Bruzual & Charlot 2003), as implemented internally in skirt. The user must simply provide the initial mass (i.e. neglecting mass loss due to stellar evolution), metallicity, and stellar age of each stellar particle. More details about galaxev can be found in Section 3.3. On the other hand, stellar populations younger than 10 Myr are essentially treated as starbursting regions, and their SEDs are modelled with the mappings-iii photoionization code (Groves et al. 2008), which includes emission from H ii regions and their surrounding photodissociation regions as well as absorption by gas and dust in the ‘birth clouds’ associated with such star-forming regions. The mappings-iii models require the following five parameters for each star-forming region: (i) the star formation SFR, assumed to be constant over the last 10 Myr, (ii) the metallicity, (iii) the compactness parameter |$\mathcal {C}$|⁠, (iv) the ISM pressure P0, and (v) the cloud covering fraction, fPDR. For every young stellar particle, we assume that the SFR is given by its initial mass divided by 10 Myr, and use its nominal metallicity (inherited from its parent gas cell). Since changes in the compactness parameter |$\mathcal {C}$| or in the ISM pressure P0 only have a noticeable effect on the far-infrared regime of the SEDs, which we do not consider in this work, we simply use ‘typical’ fixed values of |$\log _{10} \mathcal {C} = 5$| and |$\log _{10} [(P_{0}/k_{\rm B}) / {\rm cm^{-3} \, K}] = 5$| (Groves et al. 2008). Following Jonsson et al. (2010), we adopt a value of fPDR = 0.2 for the cloud covering fraction. We note that, due to the coarse sampling of star formation in cosmological simulations, the spatial distribution of young stellar populations becomes sensitive to stochastic noise (see Trayford et al. 2015; Nelson et al. 2018, for a discussion). This motivated Camps et al. (2016) and Trayford et al. (2017), using the eagle cosmological simulation (Crain et al. 2015; Schaye et al. 2015), to implement a post-processing resampling technique to model the expected frequency of young stellar populations at a higher resolution. For simplicity, we do not resample the young stellar populations, and defer a detailed quantification of the impact of such resampling techniques to future work. Finally, as mentioned in Section 3.1, every stellar particle is associated with a smoothing scale equal to the distance to its 32nd nearest neighbour. This scale is used internally by skirt to smooth the stellar mass density distribution, using the 3D version of the SPH kernel introduced in equation (1). 3.2.2 Dust modelling One remarkable feature of skirt is that it provides the option of performing the radiative transfer calculations directly on a 3D Voronoi mesh (Camps, Baes & Saftly 2013), which is particularly well suited for hydrodynamic simulations run with arepo (Springel 2010). This means that we can reconstruct the gas (and metal) density distributions in a self-consistent fashion, i.e. exactly as implemented in the hydrodynamic solver (except for the cell gradients) in order to dynamically evolve the system. Therefore, for each simulated galaxy, we consider a cubical region around it with the same dimensions as the field of view (Section 3.1). The positions of the gas cells – which are actually mesh-generating points – are used by skirt to uniquely reconstruct the Voronoi mesh inside this volume.3 This, along with the recorded density value for each gas cell, completely determines a gas density distribution for each galaxy. In order to convert the gas density distribution into a dust density distribution, we assume that the diffuse dust content of each galaxy is traced by the star-forming gas.4 The latter consists of gas cells above a given density threshold (⁠|$n_{\rm H} \sim 0.1 \, {\rm cm}^{-3}$|⁠), which are stochastically converted into stellar particles (for details, see Springel & Hernquist 2003; Vogelsberger et al. 2013; Pillepich et al. 2018a). The dust densities of gas cells that are not star forming, as well as gas cells that are not part of the galaxy in question, are set to zero. For the remaining gas cells, we assume a constant dust-to-metal mass ratio of 0.3 (Camps et al. 2016). This completely determines the overall dust density distribution. Finally, skirt offers several options to model the dust composition. Following Camps et al. (2016) and Trayford et al. (2017), we choose the multicomponent dust mix of Zubko, Dwek & Arendt (2004) including graphite grains, silicate grains, and polycyclic aromatic hydrocarbons. We note that the diffuse dust density distribution described in this section is separate from the dust content in birth clouds, which is already included in the mappings-iii models for star-forming regions (see Section 3.2.1) before any radiative transfer calculations are done. Since birth clouds have a finite lifetime (∼10 Myr), starlight from older stellar populations is only expected to be absorbed and scattered by diffuse dust in the ambient ISM, which we model using skirt. 3.2.3 Running skirt skirt is a Monte Carlo radiative transfer code, which means that the radiation field is represented as a discrete, but very large, number of ‘photon packages’ (also referred to as ‘rays’). Typical values for the number of rays (per wavelength) range from 106 to 108. However, since the sizes of our synthetic images are proportional to the sizes of the galaxies (see Section 3.1), choosing any fixed number of rays for the entire galaxy sample would not be optimal. Instead, we define a constant ‘ray density’ of 50 rays pixel−1, which we find to be more than sufficient to carry out the morphological measurements described in Section 4. Similar values for the ray density have been used implicitly in other recent studies (e.g. Snyder et al. 2015b; Trayford et al. 2017). Our wavelength grid consists of 50 logarithmically spaced wavelengths between 0.35 and 0.95 |$\mu$|m, an interval that encloses the rest-frame and observer-frame PS1 g, r, i, |$z$| bands. Although this spectral resolution is not fine enough to adequately resolve some of the narrow emission lines included in the mappings-iii models (Section 3.2.1), we find that it yields converged results for the purposes of this paper, such that doubling the resolution of the wavelength grid does not change any of our measurements appreciably. Note that since we are only interested in optical and near-infrared frequencies, thermal dust emission is not included in our analysis. 3.2.4 Broad-band integration The main data product generated by skirt is a 3D data cube for each galaxy, consisting of a full rest-frame SED for each pixel. We then assume that the source is located at |$z$| = 0.0485 and generate the data cube that would be measured by a local observer, taking cosmological effects into account (e.g. surface brightness dimming). Each SED is then multiplied by each of the PS1 g, r, i, |$z$| filter curves (Tonry et al. 2012) and integrated over the full wavelength range. This procedure results in idealized synthetic images for the PS1 g, r, i, |$z$| filters as would be hypothetically observed by an instrument with a point-like PSF (i.e. a Dirac delta function) and an infinite signal-to-noise (S/N) ratio. Fig. 3 shows examples of such idealized synthetic images for randomly selected IllustrisTNG galaxies, showing the light distribution that would be obtained without and with dust modelling (first two columns), and separating the dust-modelled light distribution from the second column into non-scattered and scattered rays (third and fourth columns, respectively). Figure 3. View largeDownload slide Idealized synthetic images (using the PS1 g, r, i filters) of randomly selected |$M_{\ast } \approx 10^{10.5} \, {\rm M}_{\odot }$| galaxies from the IllustrisTNG simulation at |$z$| = 0.0485, showing the different light components. From left to right, the different columns show (a) the light distribution that would be observed without any dust effects, (b) the light distribution obtained with the full dust modelling, (c) the contribution of rays that have not been scattered, and (d) the contribution from rays that have been scattered at least once. The field of view in each panel corresponds to 10 stellar half-mass radii. Note that the inclusion of dust tends to reduce the light concentration near the galactic centre. Figure 3. View largeDownload slide Idealized synthetic images (using the PS1 g, r, i filters) of randomly selected |$M_{\ast } \approx 10^{10.5} \, {\rm M}_{\odot }$| galaxies from the IllustrisTNG simulation at |$z$| = 0.0485, showing the different light components. From left to right, the different columns show (a) the light distribution that would be observed without any dust effects, (b) the light distribution obtained with the full dust modelling, (c) the contribution of rays that have not been scattered, and (d) the contribution from rays that have been scattered at least once. The field of view in each panel corresponds to 10 stellar half-mass radii. Note that the inclusion of dust tends to reduce the light concentration near the galactic centre. 3.3 galaxev pipeline Since radiative transfer calculations are computationally intensive, and not particularly worthwhile for galaxies with little dust content, we have also developed a synthetic image generation pipeline based solely on the stellar population synthesis code galaxev (Bruzual & Charlot 2003). In the absence of dust, this pipeline produces results that are effectively indistinguishable to those of our skirt pipeline, and yet is at least two orders of magnitude faster due to the absence of radiative transfer calculations. 3.3.1 Stellar sources In order to model the spectra of stellar particles in the simulation, we consider the ‘default’ simple stellar population (SSP) models included in galaxev, which were computed using the ‘Padova 1994’ evolutionary tracks and a Chabrier (2003) initial mass function. These models provide the rest-frame luminosity per unit wavelength of an SSP, Lλ(λ, t, Z), as a function of wavelength λ, age t, and metallicity Z (the luminosity is normalized to a mass of |$1 \, {\rm M}_{\odot }$|⁠). In practice, this quantity is provided as a ‘data cube’ sampled at 221 unequally spaced SSP ages between 0 and 20 Gyr, seven metallicities between 10−4 and 0.1,5 and 1221 wavelengths between |$91 \, {\rm \mathring{\rm A} }$| and |$160 \, \mu$|m.6 Since the spectra of the stellar sources in the skirt pipeline were also calculated using galaxev (Section 3.2.1), the only difference between the two pipelines at this stage lies in the modelling of the young stellar populations (<10 Myr), which is done using mappings-iii models (instead of galaxev) in the skirt pipeline. However, such young stellar populations are (by construction) practically absent in those gas-depleted galaxies where the galaxev pipeline replaces the skirt pipeline (see Fig. 2). Optionally, a simple Charlot & Fall (2000) model can be implemented at this stage in order to estimate the effects of an unresolved dust distribution (e.g. Torrey et al. 2015; Trayford et al. 2015; Nelson et al. 2018). According to this prescription, the luminosity per unit wavelength of each SSP is attenuated as |$L_{\lambda }^{\rm CF00} = L_{\lambda } \exp (-\tau _{\lambda })$|⁠, where \begin{eqnarray*} \tau _{\lambda } = \left\lbrace \begin{array}{lcc}\tau _{1} (\lambda / 5500 \, {\rm \mathring{\rm A} })^{\delta } & \quad {\rm if} & \quad t \le t_{\rm BC}, \\ \tau _{2} (\lambda / 5500 \, {\rm \mathring{\rm A} })^{\delta } & \quad {\rm if} & \quad t \gt t_{\rm BC}. \end{array}\right. \end{eqnarray*} (2) For simplicity, we adopt the recommended values for all parameters: τ1 = 1, τ2 = 0.3, δ = −0.7, and |$t_{\rm BC} = 10 \, {\rm Myr}$|⁠. We emphasize, however, that the inclusion of a Charlot & Fall (2000) model in our galaxev pipeline is optional, and in fact we only use it in Section 5 for comparison purposes. In either case, once the rest-frame luminosity per unit wavelength of a given SSP is known, we proceed to calculate the observer-frame flux per unit wavelength, taking cosmological effects into account, assuming that each SSP is located at |$z$| = 0.0485. Finally, given the observer-frame flux per unit wavelength of any SSP, we use appropriate broad-band filter response functions (Tonry et al. 2012) to obtain the integrated, apparent g, r, i, |$z$| magnitudes that would be measured by the observer. For every broad-band filter, the magnitudes are stored in a ‘table’ (sampled at the same age and metallicity values as the original galaxev data cube) that can be readily interpolated and normalized in order to obtain the apparent g, r, i, |$z$| magnitudes of any stellar particle, given its initial mass, age, and metallicity.7 3.3.2 Adaptive smoothing For each galaxy, an image is created by adding the flux contributions of all stellar particles, assuming that each of them has a spatial distribution given by the 2D version of the spline kernel given by equation (1). The adaptive smoothing length, pixel scale, and field of view are the same as those used in our skirt pipeline, and have already been described in Section 3.1. This yields idealized synthetic images that, in the absence of dust, are equivalent to those produced by our skirt pipeline (Section 3.2), but are very inexpensive to make. 3.4 Post-processing As described in the previous sections, both the skirt and galaxev pipelines produce idealized synthetic images for the PS1 g, r, i, |$z$| broad-band filters. Since the filter curves from Tonry et al. (2012) are given as capture cross-sections with units of |${\rm m}^2 \, e^{-} \, {\rm photon}^{-1}$|⁠, the units of the resulting wavelength-integrated images are |$e^{-} \, {\rm s}^{-1} \, {\rm pixel}^{-1}$| (electron counts per second per pixel). However, before performing any morphological measurements on these images, we must carry out the post-processing stages described next. 3.4.1 Convolving with a PSF Astronomical images are affected by telescope optics and, in the case of ground-based instruments, by atmospheric noise. These factors can be modelled in our synthetic images by convolving with an appropriate PSF, which for simplicity we choose to be an azimuthally symmetric Gaussian function. We match the full width at half-maximum (FWHM) of the PSF to the median seeing of the PS1 instrument, which is 1.11 arcsec (1.31 arcsec) in the i band (g band). 3.4.2 Modelling noise: σsky and G Including sky background noise is important not just for the sake of realism, but also because most of the morphological measurements described in Section 4, as well as the segmentation image described in Section 3.4.3, are built upon the notion that the ‘galaxy’ pixels are surrounded by ‘sky’ pixels with random values drawn from some probability distribution. We model sky background noise as a Gaussian random variable that is added to each pixel flux value, with a standard deviation σsky that is constant across the entire image. While σsky could be redefined for each individual image based on a ‘target’ S/N, this would be somewhat inconsistent with data from most galaxy surveys, which by design do not change the exposure time of individual observations arbitrarily. Therefore, we assume the same value of σsky for all of our simulated galaxies. We find that typical values for the least noisy quarter of POGS galaxies are σsky = 1/12 and |$1/15 \,\, e^{-} \, {\rm s}^{-1} \, {\rm pixel}^{-1}$| in the i band and g band, respectively. However, we find that directly applying this level of background noise to our synthetic images results in a somewhat large fraction (up to ∼30 per cent) of ‘bad’ measurements (see Section 4), especially in the g band, at low stellar masses, and for galaxies from the original Illustris simulation, which tend to have lower surface brightness (see also Bottrell et al. 2017a, who encountered similar issues while analysing the Illustris simulation). Since it would be instructive to include these galaxies in the analysis, at the cost of sacrificing some realism in the amount of background noise applied to the synthetic images, we choose somewhat lower values for σsky. For IllustrisTNG, we adopt values of σsky = 1/15 and |$1/25 \,\, e^{-} \, {\rm s}^{-1} \, {\rm pixel}^{-1}$| in the i band and g band, respectively, and half of those values for the original Illustris simulation. Besides an adequate amount of background noise, some of the morphological measurements presented in Section 4 (namely, the S/N ratio per pixel and the Sérsic fit parameters) also require a ‘weight map’ that captures the uncertainty in each pixel value. This includes contributions from both the random sky background noise and from Poisson statistics in the number of electron counts measured by the charge-coupled device camera, as well as other sources of error (e.g. the ‘read noise’, which for simplicity we will consider to be included in the background noise). While such a weight map could be constructed for each synthetic image and used as input for the morphology code, a simpler approach is to define a ‘gain’ factor that is used internally by the code to generate such a weight map, according to equation (33) from the galfit user’s manual.8 With a slight abuse of terminology, we define the gain G as a scalar factor that converts any image units (in this case |$e^{-} \, {\rm s}^{-1} \, {\rm pixel}^{-1}$|⁠) into |$e^{-} \, {\rm pixel}^{-1}$|⁠.9 In the case of our mock POGS images, the gain G is simply the total exposure time, for which we adopt typical values of 1200 and 750 s in the i band and g band, respectively. 3.4.3 Creating a segmentation image A segmentation image (also known as segmentation map) is a 2D, integer-valued array with the same shape as the original image, where pixels belonging to the same source or ‘segment’ are labelled with the same integer value. A value of zero is reserved for the background. Segmentation maps are often created with the popular software SExtractor (Bertin & Arnouts 1996), which includes deblending of extended, overlapping sources, as well as a neural network algorithm to separate stars from galaxies. However, since our synthetic images contain a single galaxy each, simpler tools will suffice, as described next. For each synthetic image, we create a segmentation image using the photutils photometry package.10 We estimate the sky background by measuring the median of all pixel values, clipped iteratively at 3σ until convergence, and then assuming that all pixel values 1.5σ above the sky median belong to some source. We only keep the largest source since by construction there should be a single galaxy in each image. Finally, we smooth the main source segment by convolving the segmentation image with a uniform ‘boxcar’ filter measuring 10 pixels in each dimension. 4 CALCULATING OPTICAL MORPHOLOGIES To facilitate the analysis of our synthetic and survey images, we have developed statmorph, a python package for calculating non-parametric morphologies of galaxy images, which we make publicly available.11 This tool is largely based on the IDL code described in Lotz, Primack & Madau (2004) and Lotz et al. (2006, 2008a,b) for calculating the Gini–M20 (Lotz et al. 2004) and concentration–asymmetry–smoothness (CAS, Conselice 2003) statistics, with some additional features such as calculating multimode–intensity–deviation statistics (MID, Freeman et al. 2013) and fitting 2D Sérsic profiles. The code can handle images with a single source each, which is the mode used in this work, as well as large ‘mosaic’ images with hundreds or thousands of sources. We have run statmorph on all of our Pan-STARRS galaxy images (Section 2.2), as well as on the synthetic images described in Section 3. The results are presented in Section 5. Fig. 4 shows some of the measurements carried out by statmorph for the i-band component of the simulated galaxy shown in the top panel of Fig. 3, including full dust modelling. However, before describing in detail the measurements shown in Fig. 4 (Sections 4.1–4.8), we give a brief overview of the code next. Figure 4. View largeDownload slide Some of the morphological measurements performed by statmorph, shown for the galaxy in the top row of Fig. 3 (using the full dust model). Top left: a cut-out from the original galaxy image, with pixel intensities scaled as f(x) = log (ax + 1)/log (a + 1), where a = 104. The blue point represents the point that minimizes the asymmetry (which is used as the galactic centre for most purposes), while the blue lines indicate the orientation and extent of the elliptical aperture that contributes half of the total flux. Top row, second from left: the fitted Sérsic model, with background noise added for realism. The Sérsic index, n, is indicated in the lower left-hand corner. The red point and lines show basic shape and size measurements derived from the Sérsic fit. Top row, third from left: the residual image obtained from subtracting the Sérsic model from the original image. Top right: the residual image obtained from subtracting the rotated image (with respect to the asymmetry centre) from the original image. Bottom left: the original image, with a black contour indicating the ‘galaxy’ pixels identified by the original segmentation map. The blue square shows the ‘skybox’ used to measure the properties of the background. Bottom row, second from left: the original image, with a black contour showing the pixels included in the calculation of the Gini–M20 quantities. The text labels show the values of the Gini–M20 and CAS statistics. Bottom row, third from left: the watershed segmentation image described in Section 4.6.2. Each coloured region indicates the pixels associated with each local maximum of the brightness distribution (by following the pixels’ maximum gradient paths). The blue and red points indicate the brightest and second brightest maxima of the image. The text label shows the values of the MID statistics. Bottom right: the detection mask used to calculate the shape asymmetry, AS. The cyan line shows the ‘minimal’ circle that encloses the detection mask, which is used to define rmax (Section 4.7.2). Figure 4. View largeDownload slide Some of the morphological measurements performed by statmorph, shown for the galaxy in the top row of Fig. 3 (using the full dust model). Top left: a cut-out from the original galaxy image, with pixel intensities scaled as f(x) = log (ax + 1)/log (a + 1), where a = 104. The blue point represents the point that minimizes the asymmetry (which is used as the galactic centre for most purposes), while the blue lines indicate the orientation and extent of the elliptical aperture that contributes half of the total flux. Top row, second from left: the fitted Sérsic model, with background noise added for realism. The Sérsic index, n, is indicated in the lower left-hand corner. The red point and lines show basic shape and size measurements derived from the Sérsic fit. Top row, third from left: the residual image obtained from subtracting the Sérsic model from the original image. Top right: the residual image obtained from subtracting the rotated image (with respect to the asymmetry centre) from the original image. Bottom left: the original image, with a black contour indicating the ‘galaxy’ pixels identified by the original segmentation map. The blue square shows the ‘skybox’ used to measure the properties of the background. Bottom row, second from left: the original image, with a black contour showing the pixels included in the calculation of the Gini–M20 quantities. The text labels show the values of the Gini–M20 and CAS statistics. Bottom row, third from left: the watershed segmentation image described in Section 4.6.2. Each coloured region indicates the pixels associated with each local maximum of the brightness distribution (by following the pixels’ maximum gradient paths). The blue and red points indicate the brightest and second brightest maxima of the image. The text label shows the values of the MID statistics. Bottom right: the detection mask used to calculate the shape asymmetry, AS. The cyan line shows the ‘minimal’ circle that encloses the detection mask, which is used to define rmax (Section 4.7.2). The statmorph code requires the following data: |$\tt image$|⁠: An image containing the object(s) of interest, given as a 2D array. The code assumes that the input image is already background subtracted. |$\tt segmap$|⁠: The corresponding segmentation map, given as a 2D array (of the same size as the image) with different integer values used to label different sources. A value of zero is reserved for the background. |$\tt weightmap$| or |$\tt gain$|⁠: The weight map is a 2D array (of the same size as the image) representing the 1σ variation of each pixel value, including the contribution from the sky background.12 If the weight map is not provided by the user, a ‘gain’ parameter must be provided instead, i.e. a multiplicative factor that converts the image units into electron counts per pixel. This is used internally by the code to estimate the weight map assuming Poisson statistics. The following two input arguments are optional, but their use is strongly recommended, if such data are available: |$\tt mask$|⁠: A 2D array (of the same size as the image) with boolean values indicating the pixels that should be excluded from the calculation (e.g. pixels contaminated by foreground stars). |$\tt psf$|⁠: A 2D array (usually smaller than the image) representing the PSF of the observations. This is only used for the Sérsic profile fitting. Given the above information, statmorph will perform the morphological measurements described in Sections 4.1–4.8. Before doing so, however, we remove any ‘bad pixels’ (as caused by cosmic rays or instrumental noise) from the image (Lotz et al. 2004), which we identify as pixels with flux values that deviate by more than 10σ (this value can be modified by the user) from their neighbouring pixels. In general, all of the ‘fixed’ parameters implicit in the morphological measurements defined in Sections 4.1–4.8 (e.g. the radius of the circular aperture used to measure the CAS statistics) can be modified by the user. For simplicity, however, we recommend leaving most parameters to their default values, which correspond to the ones most widely used in the literature. After the calculations are done, statmorph will notify the user if there was a problem by means of the following two ‘bad measurement’ flags, where values of ‘0’ and ‘1’ indicate good and bad measurements, respectively: |$\tt flag$|⁠: This indicates if there was a problem with the basic measurements, which can happen for a variety of reasons (e.g. when there is some artefact, foreground star, or secondary source in the image that was not properly masked, or when the image is not background subtracted). Although the code will attempt to perform all of the morphological measurements regardless, we recommend only trusting the morphologies of galaxies that have |$\tt flag == 0$|⁠. |$\tt flag\_sersic$|⁠: This indicates if there was a problem during the Sérsic profile fitting, which can happen for irregular or merging galaxies. Generally, it is not necessary to enforce this flag, unless one is actually interested in Sérsic profile fits. 4.1 Shape measurements For every source labelled in the segmentation map, some basic shape measurements are performed. First, the centroid or barycentre (xc, yc) is calculated in the usual way: \begin{eqnarray*} x_{\rm c} = \frac{\sum _{i,j}I_{ij} x_j}{\sum _{i,j} I_{ij}}, y_{\rm c} = \frac{\sum _{i,j}I_{ij} y_i}{\sum _{i,j} I_{ij}}, \end{eqnarray*} (3) where Iij is the pixel value at the ith row and jth column of the image I, (xj, yi) are the coordinates at the centre of the pixel, and the sum is carried out over all the pixels belonging to the source, as specified by the input segmentation map. In practice, the centroid does not always correspond to the ‘visual’ centre of a galaxy, especially for merging systems or objects with extended tidal features. Therefore, the centroid is mostly used as an initial guess for more robust measures of the galactic centre, such as the point that minimizes the asymmetry index A, which is described in Section 4.5.2. For either choice of the galactic centre, all central moments up to second order are calculated: \begin{eqnarray*} \mu _{pq} = \sum _{i,j} I_{ij} (x_j - x_{\rm c})^p (y_i - y_{\rm c})^q, \end{eqnarray*} (4) where p, q ∈ {0, 1, 2}. In particular, the second-order moments are used to construct the covariance matrix of the light distribution I: \begin{eqnarray*} {\rm cov}({\bf I}) = \frac{1}{\mu _{00}} \left(\begin{array}{cc}\mu _{20} & \quad \mu _{11} \\ \mu _{11} & \quad \mu _{02} \end{array}\right) \equiv \left(\begin{array}{cc}\overline{x^2} & \quad \overline{xy} \\ \overline{xy} & \quad \overline{y^2} \end{array}\right). \end{eqnarray*} (5) This completely determines a 2D Gaussian distribution with second moments equal to those of the original source. The orientation of the source, i.e. the angle between the x-axis and the major axis of the corresponding Gaussian distribution, is given by \begin{eqnarray*} {\rm orientation} = \frac{1}{2} \arctan \left(\frac{2 \, \overline{xy}}{\overline{x^2} - \overline{y^2}}\right). \end{eqnarray*} (6) Similarly, if λ1 and λ2 are the eigenvalues of the covariance matrix (5), where λ1 ≥ λ2, then the elongation and ellipticity of the source can be calculated as \begin{eqnarray*} {\rm elongation} = \frac{a}{b} \end{eqnarray*} (7) and \begin{eqnarray*} {\rm ellipticity} = 1 - \frac{b}{a}, \end{eqnarray*} (8) where |$a \equiv \sqrt{\lambda _1}$| and |$b \equiv \sqrt{\lambda _2}$|⁠. Equivalently, closed-form solutions for a and b can be found in equations (24–25) from the SExtractor user’s manual.13 Finally, we note that an alternative shape measurement included in statmorph is the ellipticity derived from Sérsic profile fitting, which is described in Section 4.8. 4.2 Size measurements The main size measurements carried out by statmorph are the half-light radius rhalf and the Petrosian radius rpetro (Petrosian 1976), although other size measurements are necessary for some of the morphological statistics described next, and are therefore included as well. This includes r20 and r80 (Section 4.5.1), rmax (Section 4.7.2), and the Sérsic half-light radius (Section 4.8). Most of our size measurements are carried out for both circular and elliptical apertures (using the orientation and ellipticity from Section 4.1). For the sake of readability, we will often use the term ‘radius’ both for the radius of a circle and for the semimajor axis of an ellipse, unless noted otherwise. 4.2.1 Half-light radius The half-light radius, rhalf, contains half of the light emitted by the galaxy. We calculate it both for circular and elliptical apertures, setting the centre at the point that minimizes the asymmetry index A (Section 4.5.2). There is some ambiguity in defining the outer radius, which is assumed to contain the galaxy’s total flux. In an attempt to capture as much of the galaxy’s light as possible, we define the outer radius at rmax (Section 4.7.2). However, we note that this definition of rhalf is sensitive to the limiting surface brightness and S/N of the observations. As an alternative, we describe r50 in Section 4.5.1, which represents a measure of the half-light radius consistent with r20 and r80 from the CAS statistics. Finally, in Section 4.8 we present the Sérsic half-light radius, which contains half of the light of a fitted 2D Sérsic profile. 4.2.2 Petrosian radius The Petrosian radius (Petrosian 1976), rpetro, is the radius at which the mean surface brightness is equal to some fraction η of the mean surface brightness within rpetro. The Petrosian radius is largely insensitive to variations in the limiting surface brightness and S/N of the observations, which makes it ideal for comparing galaxy properties across different redshifts. Typically, the Petrosian ‘ratio’ is set to η = 0.2. We calculate the Petrosian radius both for circular and elliptical apertures, assuming that the centre coincides with the point that minimizes the asymmetry index A (Section 4.5.2). 4.3 Noise and flux measurements 4.3.1 Sky background properties It is generally useful to quantify some properties of the sky background, which is defined as the set of pixels that do not belong to any source (i.e. pixels with a value of zero in the segmentation map) and that are also not masked. For most purposes, it is sufficient to simply search for a square region in the image that contains only background pixels (however, see Shi et al. 2009, for a more sophisticated approach). Our code finds such a ‘skybox’ automatically, searching first for regions measuring 32 × 32 pixels and iteratively halving the dimensions of the skybox candidates until a blank region of the image is found. Once an appropriate skybox has been determined, the mean, median, and standard deviation of its pixel flux values are calculated. The skybox is also used in Section 4.5 to measure the asymmetry and smoothness of the background. 4.3.2 S/N per pixel As mentioned at the beginning of Section 4, the weight map quantifies the amount of variation in each pixel value. Therefore, we calculate the average S/N per pixel, 〈S/N〉, as the average ratio between the image and the weight map over some aperture. For consistency with Lotz et al. (2004), such an aperture is given by the Gini segmentation map defined in Section 4.4.1, which roughly corresponds to an elliptical aperture at the Petrosian radius. 4.3.3 Petrosian flux In order to estimate galaxy magnitudes and colours, it is useful to output the sum of the pixel flux values over some aperture. Here, we follow the SDSS convention and measure the Petrosian flux over an aperture with a radius equal to twice the Petrosian radius (Stoughton et al. 2002). We measure the Petrosian flux both for circular and elliptical apertures. 4.4 Gini–M20 statistics The Gini–M20 classification system (Lotz et al. 2004) has been used extensively to quantify galaxy morphologies, including not just ‘normal’ galaxies but also irregular and merging galaxies. In particular, Lotz et al. (2008b) defined a simple partitioning of the Gini–M20 diagram that can be used to roughly separate early-type, late-type, and merging galaxies at 0.2 < |$z$| < 0.4. A generalization of this partitioning scheme leads to the definition of the Gini–M20 ‘bulge statistic’ (Section 4.4.3) and the Gini–M20 ‘merger statistic’ (Section 4.4.4). 4.4.1 Gini coefficient The Gini coefficient (G) is a statistic traditionally used in economics to quantify the wealth inequality in a population, which has recently found an application in astronomy as well (Abraham, van den Bergh & Nair 2003; Lotz et al. 2004). For a set of n pixel flux values Xi, where i = 1, 2, ..., n, the Gini coefficient can be computed as (Glasser 1962) \begin{eqnarray*} G = \frac{1}{\bar{X} n (n-1)} \sum _{i=1}^{n} (2i - n - 1) X_i, \end{eqnarray*} (9) where |$\bar{X}$| represents the mean over the pixel values. A value of G = 1 is obtained when all of the flux is concentrated in a single pixel, while a homogeneous brightness distribution yields G = 0. In practice, the Gini coefficient is sensitive to the choice of pixels that are included in the calculation. Originally, Abraham et al. (2003) considered pixels that lie above a constant surface brightness threshold, which makes a comparison between low- and high-redshift galaxies difficult. Therefore, Lotz et al. (2004) proposed a method to construct a ‘Gini’ segmentation map that depends only on the Petrosian radius, which is therefore insensitive to the (1 + |$z$|⁠)4 surface brightness dimming of distant galaxies. This segmentation map is created in the following way. First, the galaxy image is convolved with a Gaussian kernel with σ = rpetro/5, where rpetro is the elliptical Petrosian ‘radius’ (i.e. the semimajor axis). Then, the mean surface brightness at rpetro is used to define a flux threshold, so that pixels with flux values above this threshold are assigned to the galaxy. Note that while the shape of the Gini segmentation map is quite smooth by construction, G is always calculated for the original, unsmoothed image.14 Lotz et al. (2004) also noted that the traditional definition of the Gini coefficient, as given by equation (9), becomes unreliable for very noisy galaxy images, especially at 〈S/N〉 ≲ 3, due to the inclusion of increasingly negative pixel flux values. They found that a first-order correction consists in calculating G for the absolute values of the pixel flux values, effectively redefining the Gini coefficient as \begin{eqnarray*} G = \frac{1}{\overline{|X|} n (n-1)} \sum _{i=1}^{n} (2i - n - 1) |X_i|, \end{eqnarray*} (10) where |$\overline{|X|}$| is the mean of the absolute values |Xi|. This formula can recover the ‘true’ Gini coefficient to within 10 per cent for galaxies with 〈S/N〉 = 2–3 (Lotz et al. 2004). 4.4.2 M20 The M20 statistic (Lotz et al. 2004) measures the second moment of a galaxy’s brightest regions, containing 20 per cent of the total flux, relative to the total second-order central moment, μtot. The latter is calculated as \begin{eqnarray*} \mu _{\rm tot} = \sum _{i=1}^n \mu _i \equiv \sum _{i=1}^n I_i \left[(x_i - x_{\rm c})^2 + (y_i - y_{\rm c})^2\right], \end{eqnarray*} (11) where Ii are the pixel flux values, (xc, yc) is the galaxy’s centre, and the sum is carried out over all the pixels identified by the Gini segmentation map (Section 4.4.1). The centre (xc, yc) is the point that minimizes μtot, which corresponds to the centroid of the pixels labelled in the segmentation map. The second moment of the galaxy’s brightest regions is obtained by sorting the pixels by flux and summing μi over the brightest pixels until the sum of the brightest pixels equals 20 per cent of the galaxy’s total flux. The result is then normalized by μtot, so that M20 is ultimately defined as (Lotz et al. 2004) \begin{eqnarray*} M_{20} \equiv \log _{10}\left(\frac{\sum _i \mu _i}{\mu _{\rm tot}}\right), \,\,\,\, {\rm while} \,\,\,\, \sum _i I_i \,\,\lt\,\, 0.2 I_{\rm tot}, \end{eqnarray*} (12) where Itot is the total flux of the pixels identified by the segmentation map. 4.4.3 The bulge statistic Lotz et al. (2008b) classified 0.2 < |$z$| < 0.4 merging galaxies as objects on the G–M20 diagram such that G > −0.14M20 + 0.33, and defined early-type galaxies (including E/S0/Sa Hubble types) as non-merging galaxies that also satisfy G > 0.14M20 + 0.80. Note that the corresponding lines on the G–M20 diagram have slopes of −0.14 and 0.14, respectively, and that they intersect at \begin{eqnarray*} G_0 = 0.565, M_{20,0} = -1.679. \end{eqnarray*} (13) The Gini–M20 ‘bulge statistic’ (Snyder et al. 2015b) is defined as the position along a line (with values increasing towards bulge-dominated systems) with origin15 at (G0, M20,0) that is perpendicular to the line separating early-type and late-type galaxies, scaled by a factor of 5: \begin{eqnarray*} F(G, M_{20}) = -0.693 M_{20} + 4.95 G - 3.96. \end{eqnarray*} (14) This quantity is strongly correlated with other measures of galactic bulge strength, such as the concentration index (Section 4.5.1) and the Sérsic index (Section 4.8). 4.4.4 The merger statistic Similarly, the Gini–M20 ‘merger statistic’ (Snyder et al. 2015a) is the position along a line (with values increasing towards merging systems) with origin at (G0,M20,0) that is perpendicular to the line that divides merging and non-merging galaxies: \begin{eqnarray*} S(G, M_{20}) = 0.139 M_{20} + 0.990 G - 0.327. \end{eqnarray*} (15) We note that Lotz et al. (2004) originally defined slightly different division lines, tuned to observations of |$z$| = 0 galaxies (instead of the more distant 0.2 < |$z$| < 0.4 galaxies considered in Lotz et al. 2008b), which would perhaps be a better choice for this work (we consider Pan-STARRS galaxies at |$z$| ∼ 0.05). Instead, we adopt the classification scheme of Lotz et al. (2008b) in order to provide a consistent, ‘general-purpose’ definition of the G–M20 bulge and merger statistics that can be applied over a wider range of redshifts (Snyder et al. 2015a,b). 4.5 CAS statistics Another set of widely used non-parametric morphological indicators consists of the CAS statistics. Although galaxy concentration has been measured as early as Morgan (1958, 1959), while the asymmetry and smoothness indices have predecessors in the works of Rix & Zaritsky (1995) and Takamiya (1999), respectively, these quantities have evolved over the years, until reaching their current form as described in Conselice (2003). 4.5.1 Concentration The concentration index (C) is usually defined as (Bershady, Jangren & Conselice 2000; Conselice 2003) \begin{eqnarray*} C = 5\log _{10}\left(\frac{r_{80}}{r_{20}}\right), \end{eqnarray*} (16) where r20 and r80 are the radii of circular apertures containing 20 and 80 per cent, respectively, of the galaxy’s light.16 There is some ambiguity in defining the aperture that contains the ‘total’ flux of a galaxy. Bershady et al. (2000) defined the total flux as the flux contained within an aperture with radius equal to 2rpetro. However, more recent studies (Conselice 2003; Lotz et al. 2004) measure the total flux within 1.5rpetro, which is also the definition that we adopt here. The centre of the aperture corresponds to the point that minimizes the asymmetry index A, described next. 4.5.2 Asymmetry The asymmetry index (A) is obtained by subtracting the galaxy image rotated by 180° from the original image (Schade et al. 1995; Abraham et al. 1996; Conselice, Bershady & Jangren 2000): \begin{eqnarray*} A = \frac{\sum _{i,j}|I_{ij} - I_{ij}^{180}|}{\sum _{i,j}|I_{ij}|} - A_{\rm bgr}, \end{eqnarray*} (17) where Iij and |$I_{ij}^{180}$| are the pixel flux values of the original and rotated images, respectively, and Abgr is the average asymmetry of the background. The sum is carried out over all pixels within 1.5rpetro of the galaxy’s centre, which is determined by minimizing A. 4.5.3 Smoothness The ‘clumpiness’ or smoothness index (S) is obtained by subtracting the galaxy image smoothed with a boxcar filter of width σ from the original image (Conselice 2003): \begin{eqnarray*} S = \frac{\sum _{i,j}I_{ij} - I_{ij}^{S}}{\sum _{i,j}I_{ij}} - S_{\rm bgr}, \end{eqnarray*} (18) where Iij and |$I_{ij}^{S}$| are the pixel flux values of the original and smoothed images, respectively, and Sbgr is the average smoothness of the background. The sum is carried out over all pixels at distances between σ and 1.5rpetro from the point that minimizes the asymmetry index A. The central region is excluded because most galaxies are highly concentrated. In addition, pixels such that |$I_{ij} - I_{ij}^{S}\,\, \lt\,\, 0$| are excluded from the numerator in equation (18), so that only high-frequency features that are brighter than the local average contribute to S. Following Lotz et al. (2004), we set σ = 0.25rpetro (instead of the original choice of σ = 0.3rpetro). Note that larger values of S actually correspond to galaxies that are less smooth (i.e. more ‘clumpy’). 4.6 MID statistics The MID statistics (Freeman et al. 2013; Peth et al. 2016) were introduced as an alternative to the Gini–M20 and CAS statistics that is plausibly more sensitive to recent mergers. However, these new statistics have not been tested as extensively as the Gini–M20 and CAS statistics, especially using hydrodynamic simulations (e.g. Lotz et al. 2008a, 2010a,b; Bignone et al. 2017). The MID statistics are calculated over a segmentation map that can be considered to be a generalization of the ‘Gini’ segmentation map defined in Section 4.4.1, removing the assumption of ellipticity. It is constructed by finding a surface brightness threshold such that its value equals a fraction η of the mean surface brightness of the pixels above the threshold, where typically η = 0.2. The main source is then identified as the set of connected pixels (including the brightest pixel) with flux values above the threshold. The shape of the segmentation map is then regularized slightly with a 3 × 3 boxcar filter. 4.6.1 Multimode The multimode statistic (M) measures the ratio between the areas of the two most ‘prominent’ clumps within a galaxy. Its calculation mostly consists in finding such substructures, which is done as follows. First, all pixels within the MID segmentation map are sorted by brightness. Then, for a given quantile q (between 0 and 1), the set of all pixels with flux values above the qth quantile will generally consist of n groups of contiguous pixels, which are sorted by area (largest first): Aq, 1, Aq, 2, ..., Aq, n. Finally, M is defined as the quantile q that maximizes the area ratio17 between the two largest groups (Peth et al. 2016): \begin{eqnarray*} M = \max _{q}\left(\frac{A_{q,2}}{A_{q,1}}\right). \end{eqnarray*} (19) 4.6.2 Intensity The intensity statistic (I) measures the ratio between the two brightest subregions of a galaxy. To calculate it, the galaxy image is first slightly smoothed using a Gaussian kernel with σ = 1 pixel. Then, the image is partitioned into pixel groups according to the watershed algorithm: each distinct subregion consists of all the pixels such that their maximum gradient paths lead to the same local maximum. We perform this operation using a function from the scikit-image image-processing package.18 A watershed segmentation image is illustrated in Fig. 4. Once the pixel groups are defined, their summed intensities are sorted into descending order: I1, I2, etc. The intensity statistic is then defined as (Freeman et al. 2013) \begin{eqnarray*} I = \frac{I_2}{I_1}. \end{eqnarray*} (20) 4.6.3 Deviation The deviation statistic (D) measures the distance between the image centroid, (xc, yc), calculated for the pixels identified by the MID segmentation map, and the brightest peak found during the computation of the I statistic, |$(x_{I_1}, y_{I_1})$|⁠. This distance is normalized by |$\sqrt{n_{\rm seg}/{\pi} }$|⁠, where nseg is the number of pixels in the segmentation map, which represents an approximate galaxy ‘radius’ (Freeman et al. 2013): \begin{eqnarray*} D = \sqrt{\frac{{\pi} }{n_{\rm seg}}} \sqrt{(x_{\rm c} - x_{I_1})^2 + (y_{\rm c} - y_{I_1})^2}. \end{eqnarray*} (21) 4.7 Other non-parametric morphologies 4.7.1 Outer asymmetry The outer asymmetry (AO) is similar to the standard asymmetry, A (Section 4.5.2), but its calculation excludes pixels within the inner elliptical aperture that contributes 50 per cent of the galaxy’s light (Wen, Zheng & An 2014). More precisely, the flux integration is carried out over an elliptical annulus with inner and outer semimajor axes ahalf and amax, which are the elliptical generalizations of rhalf (Section 4.2.1) and rmax (Section 4.7.2), respectively. The orientation and ellipticity of the annulus are calculated as described in Section 4.1. Wen et al. (2014) originally rotated each image around the centroid of the outer elliptical annulus. However, for similar reasons to those explained at the end of Section 4.7.2, we instead rotate each galaxy image around the point that minimizes the standard asymmetry, A. 4.7.2 Shape asymmetry The shape asymmetry (AS) is calculated using the same mathematical expression as the standard asymmetry (equation 17). However, the measurement is done on the binary detection mask (i.e. the segmentation map) rather than the galaxy image. This increases the sensitivity to low surface brightness features in the galaxy’s outskirts (Pawlik et al. 2016). The calculation of AS mostly consists in carefully separating the ‘galaxy’ pixels from the background. The binary detection mask or ‘shape asymmetry segmentation map’ is created as follows. First, the sky background level is estimated over a circular annulus with inner and outer radii equal to two and four times the Petrosian semimajor axis (Section 4.2.2).19 Within this aperture, the sky background level is estimated by iteratively clipping the histogram of the flux pixel values at 3σ until convergence of the ‘mode’ defined in Bertin & Arnouts (1996): mode = 2.5 × median − 1.5 × mean. Then, a brightness threshold is defined at 1σ above the mode and the galaxy image is smoothed slightly with a 3 × 3 boxcar (mean) filter. The binary detection mask is defined as the contiguous group of pixels above the threshold that includes the brightest pixel in the galaxy image. An example is shown in the bottom right-hand panel of Fig. 4. Once the detection mask has been created, we define rmax as the distance between the point that minimizes the standard asymmetry (A) and the most distant pixel that belongs to the detection mask.20 This defines a ‘minimal’ circular aperture that encloses the detection mask.21 As shown by Pawlik et al. (2016; their fig. 2), even a circular aperture at 2rpetro does not capture the low surface brightness features of some morphologically disturbed galaxies, which makes rmax a more useful size measurement in this context. Finally, AS is calculated by applying equation (17) to the binary detection mask (note that the background correction term disappears since it is zero by definition). However, instead of finding the point that minimizes AS, the rotation is done around the point that minimizes the standard (i.e. flux-weighted) asymmetry, A. This ensures that the galaxy is rotated around its core, rather than around some arbitrary point that could be far from the ‘visual’ centre of the galaxy, especially for strongly disturbed systems. 4.8 Sérsic profile fitting Although statmorph was originally created to calculate non-parametric morphologies, we have decided to include Sérsic profile fits as well, mostly due to their importance in astronomy. The Sérsic profile (Sérsic 1968) is a useful parametrization of the brightness profile of a galaxy: \begin{eqnarray*} I(r) = I_{e} \exp \left\lbrace -b_n \left[ \left(\frac{r}{r_e}\right)^{1/n} - 1\right] \right\rbrace , \end{eqnarray*} (22) where n is known as the Sérsic index, bn is a coefficient chosen so that a circular aperture with radius re contains half of the galaxy’s flux, and Ie is the surface brightness at r = re. Note that the exponential profile (which approximates the brightness profile of spiral galaxies) and the de Vaucouleurs profile (typical of elliptical galaxies) correspond to the special cases n = 1 and n = 4, respectively. Although equation (22) defines a 1D profile (i.e. an azimuthally symmetric brightness profile), it can also be used to represent a 2D profile – with elliptical isophotes – through an appropriate transformation of the coordinate system. In total, a 2D Sérsic profile has seven free parameters: Ie, re, n, the centre of the profile (which contributes two free parameters), and the ellipticity and orientation of the elliptical isophotes. There are highly sophisticated algorithms available for fitting Sérsic profiles to astronomical images, such as gim2d (Simard 1998), galfit (Peng et al. 2002), imfit (Erwin 2015), and profit (Robotham et al. 2017), among others. These codes typically provide the option of fitting other mathematical models besides the Sérsic profile, and support fitting multiple components simultaneously. However, since we are only interested in fitting a single-component Sérsic profile to each galaxy image, simpler tools are sufficient for our purposes, as summarized next. We fit a 2D Sérsic profile to each galaxy image using the astropy22 modelling package, which in turn uses fitting functions from other open source projects, in particular scipy.23 The fitting is done with the Levenberg–Marquardt algorithm, assigning a weight to each individual pixel given by the weight map described at the beginning of Section 4. During each step of the fitting, the modelled Sérsic profile is convolved with the PSF of the observations, if provided. The initial ‘guessed’ values for the fitted parameters are based on the non-parametric measurements presented in Sections 4.1 through 4.7, and are generally not far from the final, optimized values. 5 RESULTS Throughout this section, we present the results of performing the morphology measurements described in Section 4 on the synthetic images described in Section 3, as well as on our sample of real Pan-STARRS images (Section 2.2). We only show those measurements that satisfy the following properties: They should be classified as ‘good’ measurements, i.e. they should satisfy flag == 0 (Section 4). For parameters derived from Sérsic fitting we also require flag_sersic == 0. The mean S/N per pixel, 〈S/N〉, should be higher than 2.5 (Lotz et al. 2006). The radius of the circular aperture containing 20 per cent of the galaxy’s flux, r20, should be larger than half of the FWHM of the PSF. We find that ∼90 per cent of our galaxy sample satisfies these requirements in any stellar mass range, except for real Pan-STARRS galaxies at |$M_{\ast } \gtrsim 10^{11} \, {\rm M}_{\odot }$|⁠, where the percentage drops to ∼80 per cent. All of the morphological measurements presented throughout this section were obtained in the PS1 i band. We highlight that the mean FWHM of PS1 i-band observations is 1.11 arcsec, which corresponds to ∼1 kpc at |$z$| ∼ 0.05. On the other hand, the spatial resolution of the Illustris and IllustrisTNG simulations, essentially set by the gravitational softening length, is slightly below ∼1 kpc, which is consistent with the observational sample. This also happens to be approximately the spatial resolution limit at which structural measurements such as the Gini coefficient (G), M20, the concentration index (C), the asymmetry index (A), and the clumpiness index (S) are still considered to be reliable (Lotz et al. 2006). 5.1 Overview Fig. 5 shows the Gini–M20 diagram for our observational Pan-STARRS sample (left-hand panel), as well as for galaxies from the IllustrisTNG (middle panel) and original Illustris (right-hand panel) simulations. All panels show objects with log10(M*/M⊙) ≈ 9.8–11.3 at |$z$| ≈ 0.05. The bin colours indicate the median concentration index (C) of the galaxies in each bin, while the contours contain 68 and 95 per cent of the galaxy population. The orange and blue lines, introduced in Lotz et al. (2008b), are commonly used to classify 0.2 < |$z$| < 0.4 galaxies as early-type, late-type, or merger candidates, as indicated in the left-hand panel. These lines are used to define the G–M20 ‘bulge’ and ‘merger’ statistics (Snyder et al. 2015a,b). Figure 5. View largeDownload slide The Gini–M20 diagram for our galaxy sample, consisting of objects with stellar masses log10(M*/M⊙) ≈ 9.8–11.3 at |$z$| ≈ 0.05. The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and original Illustris simulations, respectively. The Gini–M20 diagram is coloured according to the median value (in each 2D bin) of the concentration index, C. The black and dark grey contours contain 68 and 95 per cent of the galaxy distribution. The orange and blue lines were defined in Lotz et al. (2008b) to roughly separate merger candidates (above the orange line) from early-type (above the blue line) and late-type (below the blue line) galaxies at 0.2 < |$z$| < 0.4. These lines are also used to define the G–M20 ‘bulge’ and ‘merger’ statistics (Sections 4.4.3 and 4.4.4). Overall, the locus of the G-M20 diagram in IllustrisTNG and the associated C values are consistent with observations. The apparent lack of merging systems in IllustrisTNG is simply a consequence of the synthetic image generation procedure, which considers a single subfind object at a time, effectively biasing against early-stage mergers. Figure 5. View largeDownload slide The Gini–M20 diagram for our galaxy sample, consisting of objects with stellar masses log10(M*/M⊙) ≈ 9.8–11.3 at |$z$| ≈ 0.05. The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and original Illustris simulations, respectively. The Gini–M20 diagram is coloured according to the median value (in each 2D bin) of the concentration index, C. The black and dark grey contours contain 68 and 95 per cent of the galaxy distribution. The orange and blue lines were defined in Lotz et al. (2008b) to roughly separate merger candidates (above the orange line) from early-type (above the blue line) and late-type (below the blue line) galaxies at 0.2 < |$z$| < 0.4. These lines are also used to define the G–M20 ‘bulge’ and ‘merger’ statistics (Sections 4.4.3 and 4.4.4). Overall, the locus of the G-M20 diagram in IllustrisTNG and the associated C values are consistent with observations. The apparent lack of merging systems in IllustrisTNG is simply a consequence of the synthetic image generation procedure, which considers a single subfind object at a time, effectively biasing against early-stage mergers. Evidently, the locus of the Gini–M20 diagram in IllustrisTNG is in very good agreement with the one derived from observations. In particular, the median concentration of IllustrisTNG galaxies within the inner contour (containing 68 per cent of the sample) ranges from C ∼ 2.5 to C ∼ 4, just like in the observations. On the other hand, the bulk of the galaxy population in the original Illustris population has significantly lower median concentrations (C ∼ 2–3), while the G–M20 locus lies at lower Gini values (i.e. the pixel flux distribution becomes more homogeneous) and higher M20 values (i.e. the brightest 20 per cent of the galaxy’s light becomes more spatially extended). This range of G and M20 values is roughly consistent with the measurements of Bignone et al. (2017) for g-band images of Illustris galaxies at |$z$| = 0. Note, however, that Fig. 5 is not directly comparable to fig. 2 from Snyder et al. (2015b), which shows galaxies at a higher mass range. The improved morphological realism of IllustrisTNG galaxies evident from Fig. 5 is primarily due to two reasons. First, the reparametrization of the galactic winds (Pillepich et al. 2018a) results in galaxies with roughly correct sizes at the low-mass end, while also producing thinner discs. Secondly, the new AGN feedback model (Weinberger et al. 2017) is able to quench galaxies more efficiently at the high-mass end, allowing dissipationless processes – in particular dry mergers – to have an impact on the morphologies of galaxies that are already quenched (Rodriguez-Gomez et al. 2017). Nevertheless, Fig. 5 also reveals the existence of IllustrisTNG galaxies with somewhat low M20 values (and correspondingly high Gini coefficients and concentrations) compared to observations. This discrepancy primarily affects massive galaxies. While a detailed investigation of this issue is beyond the scope of this paper, we offer two possible explanations. The first one is related to the numerical implementation of the galaxy formation model. When the kinetic mode AGN feedback is active, momentum is imparted in a weighted fashion to a fixed number of neighbouring gas cells (256 in the fiducial TNG100 run). It is possible that this ‘feedback injection region’ is too large, such that even if the overall galaxy becomes quenched, there is not enough quenching at the smallest galactocentric distances, leaving some residual star formation at the centres of massive galaxies. The second possibility involves the natural processes considered to be important in the formation of massive elliptical galaxies, such as dry mergers (e.g. Khochfar & Burkert 2003; Naab, Khochfar & Burkert 2006; Rodriguez-Gomez et al. 2017). The importance of these mechanisms might be exaggerated in the IllustrisTNG model due to the sharp transition between the two AGN feedback modes. This picture is supported by Fig. 2, which shows a sudden drop in the gas fractions of galaxies with |$M_{\ast } \gtrsim 10^{10.5} \, {\rm M}_{\odot }$|⁠. Another notable difference between the G–M20 diagrams of the Pan-STARRS and IllustrisTNG galaxies is that the latter seems to lack merging systems. This is not an indication of missing mergers in IllustrisTNG but is simply a consequence of how the synthetic images were created. By construction, each image contains a single galaxy identified by the subfind substructure finder, effectively biasing our simulated sample against early-stage mergers. This was done for the sake of simplicity and corresponds to the choice made in Snyder et al. (2015b). Figs 6 and 7 show median trends – as a function of stellar mass – for a wide selection of morphological, shape and size measurements (Section 4), as indicated by the y-axis label in each panel. The black lines and grey shaded regions indicate the median and the 16th to 84th percentile range (at a fixed stellar mass) for the Pan-STARRS observations. The blue and red lines correspond to the IllustrisTNG and original Illustris simulations, respectively, with different line styles indicating different models: the solid lines represent our fiducial skirt implementation, including full dust modelling and emission lines from star-forming regions, while the dashed and dotted lines correspond to simpler models (without radiative transfer) including effects from an unresolved dust distribution (i.e. Bruzual & Charlot 2003 combined with Charlot & Fall 2000) and without dust (i.e. a ‘vanilla’ Bruzual & Charlot 2003 implementation), respectively. Figure 6. View largeDownload slide Median trends of various morphological parameters, as indicated in the y-axis labels, as a function of stellar mass. The top, middle, and bottom rows show the Gini–M20, CAS, and MID statistics, respectively. In each panel, the black line shows the median trend for the observational Pan-STARRS sample, while the grey shaded region indicates the 16th to 84th percentile range of the same data. The blue and red lines show results from the IllustrisTNG and original Illustris simulations, respectively, with different line styles indicating model variations: the solid lines show results obtained with the skirt radiative transfer code including full dust modelling, while the dashed and dotted lines represent simpler models with an unresolved dust distribution and without dust, respectively. This figure shows that the optical morphologies of IllustrisTNG galaxies are in better overall agreement with observations than those of the original Illustris simulation. Figure 6. View largeDownload slide Median trends of various morphological parameters, as indicated in the y-axis labels, as a function of stellar mass. The top, middle, and bottom rows show the Gini–M20, CAS, and MID statistics, respectively. In each panel, the black line shows the median trend for the observational Pan-STARRS sample, while the grey shaded region indicates the 16th to 84th percentile range of the same data. The blue and red lines show results from the IllustrisTNG and original Illustris simulations, respectively, with different line styles indicating model variations: the solid lines show results obtained with the skirt radiative transfer code including full dust modelling, while the dashed and dotted lines represent simpler models with an unresolved dust distribution and without dust, respectively. This figure shows that the optical morphologies of IllustrisTNG galaxies are in better overall agreement with observations than those of the original Illustris simulation. Figure 7. View largeDownload slide Same as Fig. 6, but this time showing Sérsic parameters (top row), size measurements (left-hand column), shape measurements (top centre and middle panels), the Gini–M20 merger statistic (centre right-hand panel), and alternative definitions of the asymmetry index (bottom centre and bottom right-hand panels). Figure 7. View largeDownload slide Same as Fig. 6, but this time showing Sérsic parameters (top row), size measurements (left-hand column), shape measurements (top centre and middle panels), the Gini–M20 merger statistic (centre right-hand panel), and alternative definitions of the asymmetry index (bottom centre and bottom right-hand panels). Figs 6 and 7 clearly show that the IllustrisTNG simulation is in better overall agreement with observations than the original Illustris simulation, with median trends in all of the morphological measurements lying within the ∼1σ scatter of the Pan-STARRS data. It is worth noting, however, that the observational data (solid black lines) always display more monotonic trends with stellar mass, which could be caused by uncertainties (usually around ∼0.2 dex) in the stellar mass estimates. In IllustrisTNG (blue lines), the transition seen in many panels around |$M_{\ast } \approx 5 \times 10^{10} \, {\rm M}_{\odot }$| is possibly related to the onset of the kinetic mode AGN feedback, which approximately happens for galaxies above that mass (Weinberger et al. 2017). Naturally, we would find smoother trends in the simulations if we randomly perturb the stellar masses of the simulated galaxies by some amount reflecting the uncertainty in the observations. In agreement with Snyder et al. (2015b), we find that the original Illustris simulation produces qualitatively correct trends for parameters such as the Gini coefficient (G), M20, and the concentration index (C), although galaxies in general tend to be too large and disc like compared to observations. This is also consistent with the results of Bottrell et al. (2017b) for the Illustris simulation. It is also worth noting that the F(G,M20) bulge statistic, concentration index (C), and Sérsic index (n) display very similar trends with stellar mass. In fact, these quantities will be considered to be essentially interchangeable in the following sections. Finally, Figs 6 and 7 also show that the inclusion of dust effects has a non-negligible effect on galactic structure, making galaxies at the low-mass end less concentrated near their centres and more spatially extended. 5.2 Morphology and galaxy colour After demonstrating that the optical morphologies of galaxies in the IllustrisTNG simulation are essentially correct to first order for the stellar masses considered, we proceed to explore the connection between morphology and other galaxy properties, such as colour. Fig. 8 shows the F(G, M20) bulge statistic as a function of stellar mass, where the colour of each bin indicates the SDSS g − i colour index. The thick and thin black lines show the median and the 16th to 84th percentile range of F(G,M20) at a fixed stellar mass. For the observational sample, the SDSS colours were obtained by matching to the NSA catalogue as described in Section 2.2. For the simulation sample, we generated mock SDSS images following the same methodology of Section 3 (this time using SDSS broad-band filter curves and a pixel scale of 0.396 arcsec pixel−1) and then calculated the ratio between the Petrosian fluxes (integrated over a circular aperture with radius equal to twice the Petrosian radius, following the SDSS convention) in the g and i bands. Figure 8. View largeDownload slide The Gini–M20 bulge statistic, F(G, M20), as a function of stellar mass. The colour of each 2D bin corresponds to the median g − i colour index of the galaxies that fall into that bin. The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In each panel, the thick and thin lines show the median and the 16th to 84th percentile range at a fixed stellar mass. As evidenced by the colour gradients, galaxy morphology and colour are strongly correlated both in observations and in the original Illustris simulation. In the case of IllustrisTNG, however, the colour gradient becomes nearly horizontal, which means that galaxy colour is more strongly correlated with stellar mass than with morphology. The same applies to other morphological indicators, such as concentration (C) or Sérsic index (n). Figure 8. View largeDownload slide The Gini–M20 bulge statistic, F(G, M20), as a function of stellar mass. The colour of each 2D bin corresponds to the median g − i colour index of the galaxies that fall into that bin. The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In each panel, the thick and thin lines show the median and the 16th to 84th percentile range at a fixed stellar mass. As evidenced by the colour gradients, galaxy morphology and colour are strongly correlated both in observations and in the original Illustris simulation. In the case of IllustrisTNG, however, the colour gradient becomes nearly horizontal, which means that galaxy colour is more strongly correlated with stellar mass than with morphology. The same applies to other morphological indicators, such as concentration (C) or Sérsic index (n). As already mentioned in the previous section, the overall F(G,M20) values of simulated galaxies show better agreement with observations in IllustrisTNG than in Illustris original. Fig. 8 shows that this is true not only for the median trend, but also for the 1σ scatter in F(G, M20) at a fixed stellar mass, which is about twice as large in the original Illustris simulation than in the observations, while the discrepancy between IllustrisTNG and Pan-STARRS is ∼30 per cent. Fig. 8 also shows that bulge strength is very strongly correlated with galaxy colour in the observations (left-hand panel), as evidenced by the nearly vertical gradient of the g − i colour index in the M*−F(G, M20) plane. A strong correlation between morphology and colour at a fixed stellar mass is also apparent for the original Illustris simulation (right-hand panel), despite the fact that it overproduces galaxies with late-type morphologies. Surprisingly, the aforementioned colour gradient is approximately horizontal in the IllustrisTNG simulation (middle panel), which suggests that galaxy morphology and colour are relatively independent from each other in the newer model. Fig. 9 further explores the connection between the F(G,M20) bulge statistic and the g − i colour index by plotting these quantities directly against each other, as shown by the grey scale and contours in the lower panels. This figure includes the full galaxy sample, with stellar masses log10(M*/M⊙) ≈ 9.8–11.3. The histograms in the upper panels show the relative distribution of galaxy colours for early-type and late-type galaxies, defined as those with F(G, M20) ≥ 0 and F(G,M20) < 0, respectively. While such distributions are clearly different from each other in the observations and in the original Illustris model, they are quite similar in IllustrisTNG. We reach the same qualitative trends if we individually examine narrower stellar mass ranges (specifically, log10(M*/M⊙) = 9.8–10.3, 10.3–10.8, and 10.8–11.3), as well as if we replace F(G, M20) with the concentration index (C) or Sérsic index (n). These trends are also qualitatively consistent with the fact that the specific star formation rate (sSFR) and effective stellar mass surface density (Σe) are less strongly correlated in IllustrisTNG than in observations (Habouzit et al. 2018, their figs 8 and 12). Figure 9. View largeDownload slide Distribution of the F(G,M20) bulge statistic as a function of the g − i colour index for the full galaxy sample (log10(M*/M⊙) ≈ 9.8–11.3). The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In the bottom panels, the grey scale and contours show the overall galaxy distribution, while the histograms in the upper panels show the marginal distributions of early-type and late-type galaxies, defined as those with F(G,M20) ≥ 0 and F(G, M20) < 0, respectively. While IllustrisTNG features a clear colour bimodality, it also produces a relatively large proportion of red (g − i ≥ 1) discs and blue (g − i < 1) spheroids relative to observations. We obtain qualitatively similar results if we replace F(G,M20) with the concentration index (C) or Sérsic index (n). Figure 9. View largeDownload slide Distribution of the F(G,M20) bulge statistic as a function of the g − i colour index for the full galaxy sample (log10(M*/M⊙) ≈ 9.8–11.3). The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In the bottom panels, the grey scale and contours show the overall galaxy distribution, while the histograms in the upper panels show the marginal distributions of early-type and late-type galaxies, defined as those with F(G,M20) ≥ 0 and F(G, M20) < 0, respectively. While IllustrisTNG features a clear colour bimodality, it also produces a relatively large proportion of red (g − i ≥ 1) discs and blue (g − i < 1) spheroids relative to observations. We obtain qualitatively similar results if we replace F(G,M20) with the concentration index (C) or Sérsic index (n). The trends described above result in a higher fraction of red discs and blue spheroids relative to observations. Defining a red/blue cut at g − i = 1, the percentages of red discs and blue spheroids relative to the total number of discs (F(G,M20) < 0) and spheroids (F(G,M20) ≥ 0) are 35 and 21 per cent for the Pan-STARRS data, compared to 39 and 48 per cent in IllustrisTNG, respectively. If we instead define a threshold between early types and late types at F(G,M20) = −0.1, which is perhaps better suited for our low-redshift sample, the aforementioned percentages of red discs and blue spheroids would change to 28 and 24 per cent in the observations, and 37 and 49 per cent in IllustrisTNG. The formation of red discs within the IllustrisTNG framework will be studied in detail in Tacchella et al. (in preparation). 5.3 Morphology and galaxy size Fig. 10 once again shows the F(G, M20) bulge statistic as a function of stellar mass. In this case, however, the 2D bin colours indicate the median galaxy size, parametrized by the semimajor axis of an ellipse containing half of the total luminosity (Section 4.2.1), with redder colours indicating larger galaxies. The different panels show Pan-STARRS observations of galaxies at |$z$| ∼ 0.05 and the corresponding synthetic images from the IllustrisTNG and Illustris simulations. The thick and thin black lines indicate the median and 16th to 84th percentile range (at a fixed stellar mass) of F(G, M20), which are identical to those previously shown in Fig. 8. Figure 10. View largeDownload slide The Gini–M20 bulge statistic, F(G,M20), as a function of stellar mass. The colour of each 2D bin indicates the median size (specifically, the semimajor axis of an ellipse containing half of the total flux) of the galaxies that fall into that bin. The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In each panel, the thick and thin black lines show the median and the 16th to 84th percentile range at a fixed stellar mass. Both in the observations and in the original Illustris simulation, discs tend to be larger than spheroids at any given stellar mass, while the trend is less clear in IllustrisTNG. The same is true if we replace F(G, M20) with the concentration index (C) or Sérsic index (n), or if we use Sérsic half-light radii. Figure 10. View largeDownload slide The Gini–M20 bulge statistic, F(G,M20), as a function of stellar mass. The colour of each 2D bin indicates the median size (specifically, the semimajor axis of an ellipse containing half of the total flux) of the galaxies that fall into that bin. The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In each panel, the thick and thin black lines show the median and the 16th to 84th percentile range at a fixed stellar mass. Both in the observations and in the original Illustris simulation, discs tend to be larger than spheroids at any given stellar mass, while the trend is less clear in IllustrisTNG. The same is true if we replace F(G, M20) with the concentration index (C) or Sérsic index (n), or if we use Sérsic half-light radii. The left-hand panel of Fig. 10 shows a strong correlation between the F(G, M20) bulge statistic and galaxy size, with disc-like systems being larger, on average, than spheroids in any stellar mass range. This is a well-known observational trend, which is also reproduced by the original Illustris simulation (right-hand panel), despite the fact that Illustris galaxies exhibit systematically large sizes compared to observations. On the other hand, such a well-defined trend is absent in the IllustrisTNG simulation (middle panel). In Fig. 11, we further explore the relationship between morphology and galaxy size. The bottom panels show the overall distribution of galaxies in the morphology-size plane, while the upper panels show the size distribution of early-type (F(G,M20) ≥ 0) and late-type (F(G,M20) < 0) galaxies separately. We consider the full galaxy sample at |$z$| ∼ 0.05 with stellar masses log10(M*/M⊙) ≈ 9.8–11.3. This figure confirms that late-type galaxies in the observational sample (left-hand panel) and in the original Illustris simulation (right-hand panel) tend to be larger than early-type galaxies, on average, although there is significant overlap between the two populations. Interestingly, in IllustrisTNG the trend is reversed: early-type galaxies tend to be larger than late-type ones. Figure 11. View largeDownload slide Distribution of the F(G,M20) bulge statistic as a function of galaxy size (namely, the semimajor axis of an ellipse containing half of the total flux) for the full galaxy sample (log10(M*/M⊙) ≈ 9.8–11.3). The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In the bottom panels, the grey-scale and contours show the overall galaxy distribution, while the histograms in the upper panels show the marginal distributions of early-type and late-type galaxies, defined as those with F(G,M20) ≥ 0 and F(G,M20) < 0, respectively. This figure shows that the correlation between galaxy morphology and size in IllustrisTNG is in tension with the observations. We obtain similar results if we replace F(G,M20) with the concentration index (C), or if we plot the Sérsic index (n) versus the Sérsic half-light radius. Figure 11. View largeDownload slide Distribution of the F(G,M20) bulge statistic as a function of galaxy size (namely, the semimajor axis of an ellipse containing half of the total flux) for the full galaxy sample (log10(M*/M⊙) ≈ 9.8–11.3). The left-hand panel shows galaxies observed with Pan-STARRS, while the middle and right-hand panels show galaxies from the IllustrisTNG and Illustris simulations, respectively. In the bottom panels, the grey-scale and contours show the overall galaxy distribution, while the histograms in the upper panels show the marginal distributions of early-type and late-type galaxies, defined as those with F(G,M20) ≥ 0 and F(G,M20) < 0, respectively. This figure shows that the correlation between galaxy morphology and size in IllustrisTNG is in tension with the observations. We obtain similar results if we replace F(G,M20) with the concentration index (C), or if we plot the Sérsic index (n) versus the Sérsic half-light radius. We find that these qualitative trends continue to hold if we individually examine narrower stellar mass ranges (log10(M*/M⊙) = 9.8–10.3, 10.3–10.8, and 10.8–11.3), as well as if we replace F(G, M20) with other structural measurements such as the concentration index (C) or Sérsic index (n), or if we use size measurements derived from Sérsic profile fitting (Section 4.8) instead of the non-parametric size measurements currently shown. Recently, Genel et al. (2018) explored galaxy sizes in the IllustrisTNG simulation and found that the sizes of main-sequence galaxies are systematically larger than those of quenched galaxies, in agreement with observations, except at |$M_{\ast } \lesssim 10^{9.5} \, {\rm M}_{\odot }$|⁠. We have verified that we can reproduce those results in IllustrisTNG if we replace the F(G,M20) bulge statistic in Fig. 10 with the sSFR predicted by the simulation. Similarly, replacing F(G, M20) with the colour index g − i in Fig. 11 would demonstrate that blue galaxies tend to be larger than red ones (within the stellar mass range considered in this paper). In other words, while the morphology–size relation in IllustrisTNG is somewhat problematic, the sSFR-size and colour-size relations are consistent with observations. 6 DISCUSSION AND CONCLUSIONS We have made a robust, quantitative comparison between the optical morphologies of galaxies from the state-of-the-art TNG100 simulation (from the IllustrisTNG suite) and those of galaxies observed with Pan-STARRS. By repeating our analysis on the original Illustris simulation, we have assessed the degree of improvement that the IllustrisTNG galaxy formation model has achieved with respect to its predecessor in the context of image-based galaxy morphology. We generated ∼27 000 synthetic images of |$M_{\ast } \gt 10^{9.5} \, {\rm M}_{\odot }$| galaxies from the IllustrisTNG and Illustris simulations at |$z$| = 0.0485, designed to match Pan-STARRS observations of galaxies at |$z$| = 0.045–0.05. Most synthetic images were created with the radiative transfer code skirt, including the effects of dust attenuation and scattering. The radiative transfer calculations were performed on a reconstruction of the same Voronoi mesh that was used by the arepo code to dynamically evolve the system, which means that the gas density distribution – used as a proxy for the dust content of the simulated galaxies – was determined in a self-consistent fashion. Some examples of idealized synthetic images generated with skirt are presented in Fig. 3. We analysed both the synthetic images and the real Pan-STARRS images with the newly developed statmorph code, which calculates non-parametric morphological diagnostics such as the Gini–M20 (Lotz et al. 2004) and CAS (Conselice 2003) statistics, and also performs 2D Sérsic fits. Fig. 4 illustrates some of the various morphological measurements carried out for each galaxy. We find that the locus of the Gini–M20 diagram in IllustrisTNG is consistent with observations (Fig. 5), while galaxies from the original Illustris simulation occupy a region of the diagram with lower Gini and higher M20 values, which is also characterized by lower concentrations. This shows that bulge-dominated galaxies are underproduced by the original Illustris model, in agreement with Bottrell et al. (2017b). Overall, the optical morphologies of IllustrisTNG galaxies are in good agreement with Pan-STARRS observations. The median trends (at a fixed stellar mass) for the 18 morphological, shape and size parameters considered in this study are always within the ∼1σ scatter of the observations in any stellar mass range, which represents a substantial improvement with respect to the original Illustris simulation (Figs 6 and 7). We reiterate that the IllustrisTNG model was not tuned to match image-based morphologies. We also find that the inclusion of attenuation and scattering by a spatially resolved dust distribution has a non-negligible effect on the resulting galaxy morphologies, especially at the low-mass end, making galactic nuclei less concentrated and therefore increasing the galaxies’ spatial extent. However, further inspection shows that the morphology–colour (Figs 8 and 9) and morphology–size (Figs 10 and 11) relations in IllustrisTNG are qualitatively inconsistent with observations. Galactic bulge strength, quantified by the G–M20 bulge statistic, is very weakly correlated with galaxy colour in IllustrisTNG, which results in a somewhat larger fraction of red discs and blue spheroids relative to observations. Furthermore, late-type galaxies in IllustrisTNG tend to be smaller than early types, which is the opposite of what is seen in observations. (Nevertheless, we do find sSFR–size and colour–size relations that are qualitatively consistent with observations, in agreement with Genel et al. 2018.) We find similar trends with other structural measurements, such as the concentration index (C) or Sérsic index (n). The lack of a strong connection between galaxy morphology and colour in IllustrisTNG is interesting. Using the Illustris simulation, Rodriguez-Gomez et al. (2017) showed that mergers play a dominant role in shaping galaxy morphology, although only for objects with high stellar masses, which tend to be quenched (this result is reproduced in IllustrisTNG as well). This raises the question of whether mergers also play an important role in quenching galaxies, by triggering powerful starbursts and AGN activity through nuclear gas inflows (e.g. Springel, Di Matteo & Hernquist 2005). Recently, using the ‘TNG300’ simulation from the IllustrisTNG suite (carried out on a larger box than TNG100, measuring ∼300 Mpc per side, but at a lower mass resolution), Weinberger et al. (2018) found that galaxy mergers do not play an important role in quenching galaxies within the IllustrisTNG framework. Previously, however, by running zoom-in simulations of Illustris mergers, Sparre & Springel (2016) showed that the mass resolution of the Illustris simulation (which is ∼8 times higher than that of TNG300) is not enough to produce starbursts, which were observed in the merger zoom-ins (at a 40 times better mass resolution). In addition, a higher mass resolution could also benefit the current AGN feedback model. By construction, quenching is only achieved in the IllustrisTNG model when the black hole (BH) mass is above |${\sim }10^8 \, {\rm M}_{\odot }$| (Weinberger et al. 2017). At a higher resolution, a gas-rich merger would also result in more accretion to small galactocentric radii and hence in more BH growth, making the BH more likely to cross the mass threshold needed in the model for it to cause quenching. Therefore, until we have enough computing power to produce large-scale cosmological simulations at a resolution high enough to capture such starbursts and nuclear gas inflows self-consistently, the implementation of a tunable ‘merger boost factor’ into the galaxy formation model – i.e. the inclusion of some mechanism that leads to increased quenching via the merger channel – might be worth exploring. This would also have implications for the morphology–size relation since gas-rich mergers tend to make galaxies smaller (by means of a very dense star-forming core), besides contributing to bulge formation and quenching. This work highlights the importance of making fair comparisons to observations through forward modelling of simulation data, and using the same tools to analyse both simulations and observations. An interpretative framework based on hydrodynamic cosmological simulations, such as the one presented here, will be extremely valuable in order to understand galaxy formation in the era of next-generation instruments such as the Large Synoptic Survey Telescope and the James Webb Space Telescope. ACKNOWLEDGEMENTS We thank Chris Hayward, Kate Rowlands, Dries Van De Putte, Liza Sazonova, and Mike Fall for useful comments and discussions, as well as Maarten Baes and Peter Camps for making the skirt code public. VRG, JL and GS acknowledge support from the National Science Foundation (NSF) under Grant No. AST-1517559. The POGS catalogue was created with support from NSF grant AST-1412596. This work used the Extreme Science and Engineering Discovery Environment (XSEDE; Towns et al. 2014), which is supported by NSF grant ACI-1548562. The XSEDE allocation TG-AST160043 utilized the Comet and Data Oasis resources provided by the San Diego Supercomputer Center. The IllustrisTNG flagship simulations were run on the HazelHen Cray XC40 supercomputer at the High Performance Computing Center Stuttgart (HLRS) as part of project GCS-ILLU of the Gauss Centre for Supercomputing (GCS). Ancillary and test runs of the project were also run on the compute cluster operated by HITS, on the Stampede supercomputer at TACC/XSEDE (allocation AST140063), at the Hydra and Draco supercomputers at the Max Planck Computing and Data Facility, and on the MIT/Harvard computing facilities supported by FAS and MIT MKI. The original Illustris simulations were run on the Harvard Odyssey and CfA/ITC clusters, the Ranger and Stampede supercomputers at TACC/XSEDE, the Kraken supercomputer at ORNL/XSEDE, the CURIE supercomputer at CEA/France as part of PRACE project RA0844, and the SuperMUC computer at the Leibniz Computing Centre, Germany, as part of project pr85je. The Flatiron Institute is supported by the Simons Foundation. Footnotes 1 However, when showing results from the original Illustris simulation, we use 9-yr Wilkinson Microwave Anisotropy Probe observations (Hinshaw et al. 2013): Ωm = 0.2726, ΩΛ = 0.7274, Ωb = 0.0456, σ8 = 0.809, ns = 0.963, and h = 0.704. 2 http://www.nsatlas.org/data 3 Behind the scenes, this is done using the voro++ open source library for computing Voronoi tessellations (Rycroft 2009). 4 The simulations used in this work do not track dust explicitly, although impressive progress has been made in this direction (McKinnon, Torrey & Vogelsberger 2016; McKinnon et al. 2017, 2018). 5 We use a relatively recent release of galaxev that extends the upper limit of the metallicity range from 0.05 to 0.1. 6 We use the lower resolution BaSeL 3.1 spectral library instead of the higher resolution STELIB measurements, mostly for consistency with the skirt implementation of galaxev. We find that the mean variation in the PS1 g, r, i, |$z$| integrated fluxes as a result of using different spectral libraries is below a thousandth of a magnitude. 7 We clip any metallicity values that fall outside the range covered by the galaxev models. 8 https://users.obs.carnegiescience.edu/peng/work/galfit/README.pdf 9 The ‘true’ gain is a factor that converts photon counts into electron counts, which is of order one for PS1 observations. 10 https://photutils.readthedocs.io 11 https://statmorph.readthedocs.io 12 We caution that, in general, different authors adopt different definitions of the weight map, which is sometimes expressed as the variance or the inverse variance (instead of the standard deviation) of the pixel values. Furthermore, it may or may not include the contribution from the background noise, and in some cases it may even refer to an exposure map (given in seconds). 13 https://www.astromatic.net/software/sextractor 14 We also point out that the Gini segmentation map is generally expected to be continuous: when this is not the case, the ‘bad measurement’ flag is activated, although G will be calculated anyway for the main individual segment. 15 The original definition by Snyder et al. (2015b) placed the origin at a slightly different position: G0 = 0.533, M20, 0 = −1.75. However, since we are only interested in relative values and because F(G, M20) is a linear combination of G and M20, the location of the origin is irrelevant for our purposes. 16 For convenience, we also calculate r50 in an analogous fashion, which can be considered as an alternative to the definition of rhalf presented in Section 4.2.1. 17 In the original definition by Freeman et al. (2013), the area ‘ratio’ included an additional factor of Aq,2, which introduces a dependence on the size of the galaxy. This was changed by Peth et al. (2016), which is the definition followed in this work. 18 https://scikit-image.org 19 Instead, Pawlik et al. (2016) used an annulus with inner and outer radii equal to 20 and 40 times the FWHM of the galaxy profile. 20 We differ slightly from Pawlik et al. (2016), who measured rmax with respect to the brightest pixel instead of the asymmetry centre. 21 For completeness, we also calculate the elliptical generalization of rmax, i.e. the semimajor axis of the minimal elliptical aperture (with shape defined as in Section 4.1) that contains the binary mask. 22 http://www.astropy.org 23 https://www.scipy.org REFERENCES Abraham R. G. , Tanvir N. R. , Santiago B. X. , Ellis R. S. , Glazebrook K. , van den Bergh S , 1996 , MNRAS , 279 , L47 https://doi.org/10.1093/mnras/279.3.L47 Crossref Search ADS Abraham R. G. , van den Bergh S. , Nair P. , 2003 , ApJ , 588 , 218 https://doi.org/10.1086/373919 Crossref Search ADS Baes M. , Verstappen J. , De Looze I. , Fritz J. , Saftly W. , Vidal Pérez E. , Stalevski M. , Valcke S. , 2011 , ApJS , 196 , 22 https://doi.org/10.1088/0067-0049/196/2/22 Crossref Search ADS Bershady M. A. , Jangren A. , Conselice C. J. , 2000 , AJ , 119 , 2645 https://doi.org/10.1086/301386 Crossref Search ADS Bertin E. , Arnouts S. , 1996 , A&AS , 117 , 393 https://doi.org/10.1051/aas:1996164 Crossref Search ADS Bignone L. A. , Tissera P. B. , Sillero E. , Pedrosa S. E. , Pellizza L. J. , Lambas D. G. , 2017 , MNRAS , 465 , 1106 https://doi.org/10.1093/mnras/stw2788 Crossref Search ADS Bottrell C. , Torrey P. , Simard L. , Ellison S. L. , 2017a , MNRAS , 467 , 1033 https://doi.org/10.1093/mnras/stx017 Crossref Search ADS Bottrell C. , Torrey P. , Simard L. , Ellison S. L. , 2017b , MNRAS , 467 , 2879 https://doi.org/10.1093/mnras/stx276 Crossref Search ADS Bruzual G. , Charlot S. , 2003 , MNRAS , 344 , 1000 https://doi.org/10.1046/j.1365-8711.2003.06897.x Crossref Search ADS Camps P. , Baes M. , 2015 , Astron. Comput. , 9 , 20 https://doi.org/10.1016/j.ascom.2014.10.004 Crossref Search ADS Camps P. , Baes M. , Saftly W. , 2013 , A&A , 560 , A35 https://doi.org/10.1051/0004-6361/201322281 Crossref Search ADS Camps P. , Trayford J. W. , Baes M. , Theuns T. , Schaller M. , Schaye J. , 2016 , MNRAS , 462 , 1057 https://doi.org/10.1093/mnras/stw1735 Crossref Search ADS Chabrier G. , 2003 , PASP , 115 , 763 https://doi.org/10.1086/376392 Crossref Search ADS Chambers K. C. et al. , 2016 , preprint(arXiv:1612.05560) Charlot S. , Fall S. M. , 2000 , ApJ , 539 , 718 https://doi.org/10.1086/309250 Crossref Search ADS Conselice C. J. , 2003 , ApJS , 147 , 1 https://doi.org/10.1086/375001 Crossref Search ADS Conselice C. J. , Bershady M. A. , Jangren A. , 2000 , ApJ , 529 , 886 https://doi.org/10.1086/308300 Crossref Search ADS Crain R. A. et al. , 2015 , MNRAS , 450 , 1937 https://doi.org/10.1093/mnras/stv725 Crossref Search ADS Dolag K. , Borgani S. , Murante G. , Springel V. , 2009 , MNRAS , 399 , 497 https://doi.org/10.1111/j.1365-2966.2009.15034.x Crossref Search ADS Dubois Y. et al. , 2014 , MNRAS , 444 , 1453 https://doi.org/10.1093/mnras/stu1227 Crossref Search ADS Erwin P. , 2015 , ApJ , 799 , 226 https://doi.org/10.1088/0004-637X/799/2/226 Crossref Search ADS Freeman P. E. , Izbicki R. , Lee A. B. , Newman J. A. , Conselice C. J. , Koekemoer A. M. , Lotz J. M. , Mozena M. , 2013 , MNRAS , 434 , 282 https://doi.org/10.1093/mnras/stt1016 Crossref Search ADS Genel S. et al. , 2014 , MNRAS , 445 , 175 https://doi.org/10.1093/mnras/stu1654 Crossref Search ADS Genel S. et al. , 2018 , MNRAS , 474 , 3976 https://doi.org/10.1093/mnras/stx3078 Crossref Search ADS Glasser G. J. , 1962 , J. Am. Stat. Assoc. , 57 , 648 https://doi.org/10.1080/01621459.1962.10500553 Crossref Search ADS Grogin N. A. et al. , 2011 , ApJS , 197 , 35 https://doi.org/10.1088/0067-0049/197/2/35 Crossref Search ADS Groves B. , Dopita M. A. , Sutherland R. S. , Kewley L. J. , Fischera J. , Leitherer C. , Brandl B. , van Breugel W. , 2008 , ApJS , 176 , 438 https://doi.org/10.1086/528711 Crossref Search ADS Habouzit M. et al. , 2018 , preprint(arXiv:1809.05588) Hernquist L. , Katz N. , 1989 , ApJS , 70 , 419 https://doi.org/10.1086/191344 Crossref Search ADS Hinshaw G. et al. , 2013 , ApJS , 208 , 19 https://doi.org/10.1088/0067-0049/208/2/19 Crossref Search ADS Jonsson P. , 2006 , MNRAS , 372 , 2 https://doi.org/10.1111/j.1365-2966.2006.10884.x Crossref Search ADS Jonsson P. , Groves B. A. , Cox T. J. , 2010 , MNRAS , 403 , 17 https://doi.org/10.1111/j.1365-2966.2009.16087.x Crossref Search ADS Khochfar S. , Burkert A. , 2003 , ApJ , 597 , L117 https://doi.org/10.1086/379845 Crossref Search ADS Lotz J. M. , Primack J. , Madau P. , 2004 , AJ , 128 , 163 https://doi.org/10.1086/421849 Crossref Search ADS Lotz J. M. , Madau P. , Giavalisco M. , Primack J. , Ferguson H. C. , 2006 , ApJ , 636 , 592 https://doi.org/10.1086/497950 Crossref Search ADS Lotz J. M. , Jonsson P. , Cox T. J. , Primack J. R. , 2008a , MNRAS , 391 , 1137 https://doi.org/10.1111/j.1365-2966.2008.14004.x Crossref Search ADS Lotz J. M. et al. , 2008b , ApJ , 672 , 177 https://doi.org/10.1086/523659 Crossref Search ADS Lotz J. M. , Jonsson P. , Cox T. J. , Primack J. R. , 2010a , MNRAS , 404 , 575 https://doi.org/10.1111/j.1365-2966.2010.16268.x Crossref Search ADS Lotz J. M. , Jonsson P. , Cox T. J. , Primack J. R. , 2010b , MNRAS , 404 , 590 https://doi.org/10.1111/j.1365-2966.2010.16269.x Crossref Search ADS Marinacci F. et al. , 2018 , MNRAS , 480 , 5113 https://doi.org/10.1093/mnras/sty2206 McKinnon R. , Torrey P. , Vogelsberger M. , 2016 , MNRAS , 457 , 3775 https://doi.org/10.1093/mnras/stw253 Crossref Search ADS McKinnon R. , Torrey P. , Vogelsberger M. , Hayward C. C. , Marinacci F. , 2017 , MNRAS , 468 , 1505 https://doi.org/10.1093/mnras/stx467 Crossref Search ADS McKinnon R. , Vogelsberger M. , Torrey P. , Marinacci F. , Kannan R. , 2018 , MNRAS , 478 , 2851 https://doi.org/10.1093/mnras/sty1248 Crossref Search ADS Monaghan J. J. , 1992 , ARA&A , 30 , 543 https://doi.org/10.1146/annurev.aa.30.090192.002551 Crossref Search ADS Monaghan J. J. , Lattanzio J. C. , 1985 , A&A , 149 , 135 https://doi.org/10.1017/CBO9781107415324.004 Morgan W. W. , 1958 , PASP , 70 , 364 https://doi.org/10.1086/127243 Crossref Search ADS Morgan W. W. , 1959 , PASP , 71 , 394 https://doi.org/10.1086/127415 Crossref Search ADS Naab T. , Khochfar S. , Burkert A. , 2006 , ApJ , 636 , L81 https://doi.org/10.1086/500205 Crossref Search ADS Naiman J. P. et al. , 2018 , MNRAS , 477 , 1206 https://doi.org/10.1093/mnras/sty618 Crossref Search ADS Nelson D. et al. , 2015 , Astron. Comput. , 13 , 12 https://doi.org/10.1016/j.ascom.2015.09.003 Crossref Search ADS Nelson D. et al. , 2018 , MNRAS , 475 , 624 https://doi.org/10.1093/mnras/stx3040 Crossref Search ADS Pakmor R. , Bauer A. , Springel V. , 2011 , MNRAS , 418 , 1392 https://doi.org/10.1111/j.1365-2966.2011.19591.x Crossref Search ADS Pakmor R. , Springel V. , Bauer A. , Mocz P. , Munoz D. J. , Ohlmann S. T. , Schaal K. , Zhu C. , 2016 , MNRAS , 455 , 1134 https://doi.org/10.1093/mnras/stv2380 Crossref Search ADS Pawlik M. M. , Wild V. , Walcher C. J. , Johansson P. H. , Villforth C. , Rowlands K. , Mendez-Abreu J. , Hewlett T. , 2016 , MNRAS , 456 , 3032 https://doi.org/10.1093/mnras/stv2878 Crossref Search ADS Peng C. Y. , Ho L. C. , Impey C. D. , Rix H.-W. , 2002 , AJ , 124 , 266 https://doi.org/10.1086/340952 Crossref Search ADS Peth M. A. et al. , 2016 , MNRAS , 458 , 963 https://doi.org/10.1093/mnras/stw252 Crossref Search ADS Petrosian V. , 1976 , ApJ , 209 , L1 https://doi.org/10.1086/182253 Crossref Search ADS Pillepich A. et al. , 2018a , MNRAS , 473 , 4077 https://doi.org/10.1093/mnras/stx2656 Crossref Search ADS Pillepich A. et al. , 2018b , MNRAS , 475 , 648 https://doi.org/10.1093/mnras/stx3112 Crossref Search ADS Planck Collaboration XIII , 2016 , A&A , 594 , A13 https://doi.org/10.1051/0004-6361/201525830 Crossref Search ADS Rix H.-W. , Zaritsky D. , 1995 , ApJ , 447 , 82 https://doi.org/10.1086/175858 Crossref Search ADS Robitaille T. P. , 2011 , A&A , 536 , A79 https://doi.org/10.1051/0004-6361/201117150 Crossref Search ADS Robotham A. S. G. , Taranu D. S. , Tobar R. , Moffett A. , Driver S. P. , 2017 , MNRAS , 466 , 1513 https://doi.org/10.1093/mnras/stw3039 Crossref Search ADS Rodriguez-Gomez V. et al. , 2017 , MNRAS , 467 , 3083 https://doi.org/10.1093/mnras/stx305 Crossref Search ADS Rycroft C. H. , 2009 , Chaos: Interdiscip. J Nonlinear Sci. , 19 , 041111 https://doi.org/10.1063/1.3215722 Crossref Search ADS Schade D. , Lilly S. J. , Crampton D. , Hammer F. , Le Fevre O. , Tresse L. , 1995 , ApJ , 451 , L1 https://doi.org/10.1086/309677 Crossref Search ADS Schaye J. et al. , 2015 , MNRAS , 446 , 521 https://doi.org/10.1093/mnras/stu2058 Crossref Search ADS Sérsic J. L. , 1968 , Atlas de Galaxias Australes . Cordoba , Argentina Shi Y. , Rieke G. , Lotz J. , Perez-Gonzalez P. G. , 2009 , ApJ , 697 , 1764 https://doi.org/10.1088/0004-637X/697/2/1764 Crossref Search ADS Sijacki D. , Vogelsberger M. , Genel S. , Springel V. , Torrey P. , Snyder G. F. , Nelson D. , Hernquist L. , 2015 , MNRAS , 452 , 575 https://doi.org/10.1093/mnras/stv1340 Crossref Search ADS Simard L. , 1998 , in Albrecht R. , Hook R. , Bushouse H. , eds, ASP Conf. Ser. Vol. 145, Astronomical Data Analysis Software and Systems VII . Astron. Soc. Pac , San Francisco , p. 108 Snyder G. F. , Lotz J. , Moody C. , Peth M. , Freeman P. , Ceverino D. , Primack J. , Dekel A. , 2015a , MNRAS , 451 , 4290 https://doi.org/10.1093/mnras/stv1231 Crossref Search ADS Snyder G. F. et al. , 2015b , MNRAS , 454 , 1886 https://doi.org/10.1093/mnras/stv2078 Crossref Search ADS Sparre M. , Springel V. , 2016 , MNRAS , 462 , 2418 https://doi.org/10.1093/mnras/stw1793 Crossref Search ADS Springel V. et al. , 2018 , MNRAS , 475 , 676 https://doi.org/10.1093/mnras/stx3304 Crossref Search ADS Springel V. , 2010 , MNRAS , 401 , 791 https://doi.org/10.1111/j.1365-2966.2009.15715.x Crossref Search ADS Springel V. , Hernquist L. , 2003 , MNRAS , 339 , 289 https://doi.org/10.1046/j.1365-8711.2003.06206.x Crossref Search ADS Springel V. , Yoshida N. , White S. D. , 2001a , New Astron. , 6 , 79 https://doi.org/10.1016/S1384-1076(01)00042-2 Crossref Search ADS Springel V. , White S. D. M. , Tormen G. , Kauffmann G. , 2001b , MNRAS , 328 , 726 https://doi.org/10.1046/j.1365-8711.2001.04912.x Crossref Search ADS Springel V. , Di Matteo T. , Hernquist L. , 2005 , ApJ , 620 , L79 https://doi.org/10.1086/428772 Crossref Search ADS Stoughton C. et al. , 2002 , AJ , 123 , 485 https://doi.org/10.1086/324741 Crossref Search ADS Takamiya M. , 1999 , ApJS , 122 , 109 https://doi.org/10.1086/313216 Crossref Search ADS Tonry J. L. et al. , 2012 , ApJ , 750 , 99 https://doi.org/10.1088/0004-637X/750/2/99 Crossref Search ADS Torrey P. , Vogelsberger M. , Genel S. , Sijacki D. , Springel V. , Hernquist L. , 2014 , MNRAS , 438 , 1985 https://doi.org/10.1093/mnras/stt2295 Crossref Search ADS Torrey P. et al. , 2015 , MNRAS , 447 , 2753 https://doi.org/10.1093/mnras/stu2592 Crossref Search ADS Towns J. et al. , 2014 , Comput. Sci. Eng. , 16 , 62 https://doi.org/10.1109/MCSE.2014.80 Crossref Search ADS Trayford J. W. et al. , 2015 , MNRAS , 452 , 2879 https://doi.org/10.1093/mnras/stv1461 Crossref Search ADS Trayford J. W. et al. , 2017 , MNRAS , 470 , 771 https://doi.org/10.1093/mnras/stx1051 Crossref Search ADS Vinsen K. , Thilker D. , 2013 , Astron. Comput. , 3 , 1 https://doi.org/10.1016/j.ascom.2013.10.001 Crossref Search ADS Vogelsberger M. et al. , 2014a , MNRAS , 444 , 1518 https://doi.org/10.1093/mnras/stu1536 Crossref Search ADS Vogelsberger M. et al. , 2014b , Nature , 509 , 177 https://doi.org/10.1038/nature13316 Crossref Search ADS Vogelsberger M. , Genel S. , Sijacki D. , Torrey P. , Springel V. , Hernquist L. , 2013 , MNRAS , 436 , 3031 https://doi.org/10.1093/mnras/stt1789 Crossref Search ADS Waters C. Z. et al. , 2016 , preprint(arXiv:1612.05245) Weinberger R. et al. , 2017 , MNRAS , 465 , 3291 https://doi.org/10.1093/mnras/stw2944 Crossref Search ADS Weinberger R. et al. , 2018 , MNRAS , 479 , 4056 https://doi.org/10.1093/mnras/sty1733 Crossref Search ADS Wen Z. Z. , Zheng X. Z. , An F. X. , 2014 , ApJ , 787 , 130 https://doi.org/10.1088/0004-637X/787/2/130 Crossref Search ADS York D. G. et al. , 2000 , AJ , 120 , 1579 https://doi.org/10.1086/301513 Crossref Search ADS Zubko V. , Dwek E. , Arendt R. G. , 2004 , ApJS , 152 , 211 https://doi.org/10.1086/382351 Crossref Search ADS © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - The optical morphologies of galaxies in the IllustrisTNG simulation: a comparison to Pan-STARRS observations
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/sty3345
DA - 2019-03-01
UR - https://www.deepdyve.com/lp/oxford-university-press/the-optical-morphologies-of-galaxies-in-the-illustristng-simulation-a-7uaPhU1nN5
SP - 4140
VL - 483
IS - 3
DP - DeepDyve
ER -