IMPACT OF BARYONIC PHYSICS ON INTRINSIC ALIGNMENTS

The intrinsic shapes and orientations of galaxies are correlated with each other and the large-scale density field. This intrinsic alignment of galaxies is an important astrophysical systematic in weak lensing measurements (Croft & Metzler 2000; Heavens et al. 2000; Catelan et al. 2001; Jing 2002; Hirata & Seljak 2004) of upcoming surveys such as the Large Synoptic Survey Telescope⁶ (LSST; LSST Science Collaboration et al. 2009) and Euclid⁷ (Laureijs et al. 2011). Ignoring intrinsic alignments in weak lensing analysis can significantly bias the constraints on cosmological parameters such as the dark energy equation of state parameter (Krause et al. 2016). Therefore, intrinsic alignments have been studied with analytical models and also cosmological simulations including N-body and hydrodynamic simulations which can help in mitigating this contaminant signal. Analytically, intrinsic alignments have been modeled with a linear alignment model (Catelan et al. 2001; Hirata & Seljak 2004) and modifications of this model which include the nonlinear evolution of the density field (Bridle & King 2007; Blazek et al. 2015). However, it is difficult to analytically describe the alignments of a galaxy's stellar component by accurately considering the physics of galaxy formation. There are also limitations to the use of N-body simulations as one has to populate halos with galaxies by assigning a random orientation (Heymans et al. 2006) or employ semi-analytic methods (Joachimi et al. 2013). Recently, intrinsic alignments of galaxies in large-volume hydrodynamic simulations have been extensively studied with simulations of galaxy formation such as MassiveBlack-II (Khandai et al. 2015), Horizon-AGN (Dubois et al. 2014), EAGLE (Schaye et al. 2015) and Illustris (Genel et al. 2014; Vogelsberger et al. 2014b, 2014a).

Cosmological hydrodynamic simulations of galaxy formation are an important tool to study intrinsic alignments as it is possible to directly measure the shape and orientation of the stellar component of galaxies in the simulations. In a precursor of this paper, Tenneti et al. (2015b) studied the galaxy shapes and two-point statistics in the MassiveBlack-II cosmological hydrodynamic simulation. This study was extended to compare the galaxy alignments based on their morphological type in the MassiveBlack-II and Illustris simulations (Tenneti et al. 2015a). Chisari et al. (2015) used the Horizon-AGN simulation, an adaptive mesh refinement based hydrodynamic simulation of galaxy formation to study intrinsic alignments of spirals and elliptical galaxies. The redshift and luminosity evolution of alignments in the same simulation was studied in Chisari et al. (2016). Recently, Hilbert et al. (2016) studied the mass and redshift dependence of intrinsic alignments in the Illustris simulation and their dependence on stellar mass, luminosity, redshift, and photometric type. Qualitatively, the properties of galaxy shapes and alignments have a similar trend with mass across different simulations. However, differences have been noted in the amplitude of galaxy alignments and morphological fraction of disk galaxies in MassiveBlack-II and Illustris (Tenneti et al. 2015a), as well as qualitative differences in the comparison of alignments of spirals with the over-density and the redshift dependence of intrinsic alignments in Horizon-AGN (Chisari et al. 2015, 2016). Given the differences in the models of subgrid physics adopted in these simulations and also the numerical implementations of hydrodynamics, it is important to understand the details of the subgrid physics responsible for changes in the galaxy alignments and to explore the robustness of the simulation results.

In a previous study, Velliscig et al. (2015) studied intrinsic alignments using the EAGLE suite of simulations with variations in the strength of feedback. Here, we undertake a parameter space study of the subgrid model adopted in the MassiveBlack-II simulation using a suite of small-volume simulations with box size of 25 h⁻¹ Mpc on a side. We vary the free parameters in the feedback models of the simulation and test the robustness of the galaxy shapes, orientations, and two-point statistics of shape correlations to variations in these parameters. Since high-resolution hydrodynamic simulations of large volume are computationally expensive, we also test the usefulness of using small-volume simulations to capture the sensitivity of intrinsic alignment statistics to variations in the feedback parameters.

This paper is organized as follows. In Section 2, we describe the simulations used in this study along with a brief overview of the feedback models adopted in the MassiveBlack-II simulation. Section 3 provides the details of the methods adopted to calculate shapes and intrinsic alignment statistics studied in this paper. In Section 4 we compare the results from the suite of small-volume simulations with the fiducial MBII model and different amplitudes of the DC mode with those of the original 100 h⁻¹ Mpc box-size MBII simulation. The intrinsic alignment statistics in the small-volume runs with different feedback parameters are compared with those from the fiducial model in Section 5. We provide a summary of our conclusions in Section 6.

In this paper, we use the MassiveBlack-II (MBII) simulation (Khandai et al. 2015), a high-resolution cosmological hydrodynamic simulation performed in a box of volume (100 h⁻¹ Mpc)³, which includes galaxy formation physics as our base model. We complement MassiveBlack-II with smaller-volume simulations of size 25 h⁻¹ Mpc, in which we vary the key parameters of the star formation and stellar and active galactic nucleus (AGN) feedback model. We denote the smaller-volume simulations as MBII-25. The simulations are performed with the TreePM-smoothed particle hydrodynamics code, P-Gadget, a modified version of GADGET2 (Springel et al. 2005). The same version of the code was used earlier to perform the large volume MBII simulation (Khandai et al. 2015). The simulations include the wide range of physical effects thought to be crucial for properly modeling galaxy formation, such as multiphase interstellar medium (ISM), star formation, supernova and stellar wind feedback, as well as black hole accretion and feedback. Radiative cooling and heating are included as in Katz et al. (1996), along with photoheating due to an imposed ionizing UV background.

Initial conditions are generated at z = 159 and simulations are evolved to z = 0 with an equal initial number of gas and dark matter particles. The cosmological parameters are chosen with the WMAP7 cosmology (Komatsu et al. 2011): h = 0.701, Ω_m = 0.275, Ω_b = 0.046, Ω_Λ = 0.725, σ₈ = 0.816, spectral index, η_s = 0.968 The mass of each dark matter particle is 1.1 × 10⁷ h⁻¹ M_⊙. The smaller-volume simulations are performed with the same mass and spatial resolution as the original simulation. Accordingly, the initial numbers of gas and dark matter particles are equal to 2 × 1792³ and 2 × 448³ in the 100 h⁻¹ Mpc and 25 h⁻¹ Mpc box-size simulations, respectively. We note that all the small-volume simulations are started with the same initial conditions at z = 159. The details of the star formation and feedback models of the simulation and the changes adopted in the small-volume runs are described below.

2.1. Star Formation and Stellar and AGN Feedback

The star formation and feedback model adopted in the simulation is based on an earlier multiphase ISM model of Springel & Hernquist (2003). Specifically, if the local gas density ρ is greater than a critical density threshold ρ_th, a multiphase ISM consisting of cold clouds in pressure equilibrium with a hot ambient gas is assumed. The effective pressure P_eff is defined as ${P}_{\mathrm{eff}}=(\gamma -1)({\rho }_{h}{\mu }_{h}+{\rho }_{c}{\mu }_{c})$ (Springel & Hernquist 2003), where ρ_c, ρ_h are the local densities of the cold and hot phases respectively, $\rho ={\rho }_{c}+{\rho }_{h}$, and μ_h and μ_c are specific energies of the hot and cold components. The threshold density ρ_th is determined self-consistently by requiring that the effective pressure be a continuous function of density.

Star formation is modeled by spawning individual stellar particles stochastically from the cold clouds. The rate of star formation is given by

$$\begin{eqnarray}&&\displaystyle \frac{d{\rho }_{* }}{{dt}}=\displaystyle \frac{{\rho }_{c}}{{t}_{* }}-\beta \displaystyle \frac{{\rho }_{c}}{{t}_{* }}\end{eqnarray} \tag{ 1 }$$

where β = 0.1 is the mass fraction of short-lived stars and t_* is the star formation timescale with density dependence given by

$$\begin{eqnarray}&&{t}_{* }(\rho )={t}_{0}^{* }{\left(\displaystyle \frac{\rho }{{\rho }_{\mathrm{th}}}\right)}^{-0.5},\end{eqnarray} \tag{ 2 }$$

where ${t}_{0}^{* }=2.1\,\mathrm{Gyr}$.

The energy released by supernovae heats the ambient gas and the heating rate is set by the energy balance condition

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}({\rho }_{h}{\mu }_{h})=\beta \displaystyle \frac{{\rho }_{c}}{{t}_{* }}({\mu }_{\mathrm{SN}}).\end{eqnarray} \tag{ 3 }$$

Here ${\mu }_{\mathrm{SN}}=\tfrac{3}{2}{{kT}}_{\mathrm{SN}}$ where T_SN is the equivalent supernova temperature which is equal to 10⁸ K in the fiducial model.

2.2. Wind Feedback

Galactic winds are implemented with the wind velocity given by

$$\begin{eqnarray}&&{v}_{w}=\sqrt{\displaystyle \frac{2\beta \chi {\mu }_{\mathrm{SN}}}{\eta (1-\beta )}},\end{eqnarray} \tag{ 4 }$$

where χ = 1.0 is the fraction of supernova energy carried by the wind and η = 2.0 is the wind loading factor. For a given time step Δt, a gas particle is added to the wind probabilistically with the probability

$$\begin{eqnarray}&&{p}_{w}=1-\exp \left[-\displaystyle \frac{\eta (1-\beta )x{\rm{\Delta }}t}{{t}_{* }}\right].\end{eqnarray} \tag{ 5 }$$

2.3. AGN Feedback

The simulations also include the physics of black hole (BH) accretion and feedback, based on the models of Springel et al. (2005) and Di Matteo et al. (2005). BHs are treated as collisionless particles introduced into halos of mass greater than 5.0 × 10¹⁰ h⁻¹ M_⊙ at regular time intervals, separated by ${\rm{\Delta }}\mathrm{log}(a)=\mathrm{log}(1.25)$. The densest particle is converted into a seed BH of mass ${M}_{\mathrm{BH},\mathrm{seed}}=5\times {10}^{5}\,{h}^{-1}\,{M}_{\odot }$ which grows in mass by black hole accretion and mergers. The BH accretion rate is given by the modified Bondi rate formula

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{BH}}=\displaystyle \frac{4\pi \alpha {G}^{2}{M}_{\mathrm{BH}}^{2}\rho }{{({c}_{{\rm{s}}}^{2}+{v}_{\mathrm{BH}}^{2})}^{3/2}},\end{eqnarray} \tag{ 6 }$$

where ρ is the local gas density, c_s is the local speed of sound, v is the velocity of BH relative to the gas. The accretion rate is limited to two times the Eddington rate, ${\dot{M}}_{\mathrm{Edd}}$. A dimensionless parameter α is set to 100; that value has been found experimentally to approximately correct for the gas density close to the BH, which is reduced in the effective sub-resolution model of the ISM.

The AGN feedback is modeled by coupling 5% (the value chosen to match the slope in the observed ${M}_{\mathrm{BH}}\mbox{--}\sigma$ relation; Springel et al. 2005) of the bolometric luminosity radiated from the BH,

$$\begin{eqnarray}&&{L}_{\mathrm{bol}}={\epsilon }_{r}{\dot{M}}_{\mathrm{BH}}{c}^{2},\end{eqnarray} \tag{ 7 }$$

with the radiation efficiency _r = 0.1. The energy is deposited isotropically to the 64 nearest gas particles within the BH particle kernel.

2.4. Parameter Space Study

In this section, we describe the method adopted to calculate shapes and the also provide details of the intrinsic alignment statistics explored in this paper.

3.1. Calculation of Shapes

The 3D shapes of the dark matter and stellar components in subhalos are determined using the the eigenvalues and eigenvectors of the reduced inertia tensor given by

$$\begin{eqnarray}&&{\widetilde{I}}_{{ij}}=\displaystyle \frac{{\sum }_{n}{m}_{n}({x}_{{ni}}{x}_{{nj}})/{r}_{n}^{2}}{{\sum }_{n}{m}_{n}},\end{eqnarray} \tag{ 8 }$$

where the summation is over particle index n, and

$$\begin{eqnarray}&&{r}_{n}^{2}=\displaystyle \frac{{x}_{n0}^{2}}{{a}^{2}}+\displaystyle \frac{{x}_{n1}^{2}}{{b}^{2}}+\displaystyle \frac{{x}_{n2}^{2}}{{c}^{2}}.\end{eqnarray} \tag{ 9 }$$

Here a, b, and c are half-lengths of the principal axes of the ellipsoid.

The eigenvectors of the inertia tensor are ${\hat{e}}_{a},{\hat{e}}_{b},{\hat{e}}_{c}$ with corresponding eigenvalues λ_a > λ_b > λ_c. The eigenvectors represent the principal axes of the ellipsoid, with the half-lengths of the principal axes (a, b, c) given by $(\sqrt{{\lambda }_{a}},\sqrt{{\lambda }_{b}},\sqrt{{\lambda }_{c}})$. The 3D axis ratios are b/a and c/a.

Similarly, in 2D, the projected shapes are calculated by projecting the positions of the particles onto the XY plane and modeling the shapes as ellipses. Here, we denote the eigenvectors as ${\hat{e}}_{a}^{\prime },{\hat{e}}_{b}^{\prime }$ with corresponding eigenvalues ${\lambda }_{a}^{\prime }\gt {\lambda }_{b}^{\prime }$. The lengths of the semimajor and semiminor axes are $a^{\prime} =\sqrt{{\lambda }_{a}^{\prime }}$ and $b^{\prime} =\sqrt{{\lambda }_{b}^{\prime }}$ with the axis ratio b'/a'.

The details of the iterative method for measuring axis ratios can be found in Tenneti et al. (2015b). In the first iteration, we start with the half-lengths of the principal axes all equal to 1 and determine the eigenvalues and eigenvectors of the ellipsoid. After each iteration, the lengths of the principal axes of ellipsoids are rescaled such that the enclosed volume is constant and particles outside the ellipsoidal volume are discarded. This process is repeated until convergence is reached such that the fractional change in axis ratios is below 1%.

In addition to the distribution of the axis ratios, b/a and c/a of the stellar components of subhalos, we are also interested in the orientation of the major axis of the stellar shape with the shape of dark matter in subhalos. So, we compute the probability distribution of the misalignment angle

$$\begin{eqnarray}&&{\theta }_{m}=\arccos (| {\hat{e}}_{{da}}\cdot {\hat{e}}_{{ga}}),\end{eqnarray} \tag{ 10 }$$

where ${\hat{e}}_{{da}}$ and ${\hat{e}}_{{ga}}$ are the major axes of the shapes defined by the dark matter and stellar matter components respectively.

3.2. Two-point Statistics

In this paper we quantify the intrinsic alignments of galaxies with the large-scale density field using the ellipticity–direction (ED) and the projected shape–density (${w}_{\delta +}$) correlation functions.

The ED correlation function cross-correlates the orientation of the major axes of the 3D shapes of dark matter or stellar component of galaxies with the large-scale density field. Consider a subhalo centered at position x with the major axis direction ${\hat{e}}_{a}$. Let the unit vector in the direction of a tracer of the matter density field at a distance r be $\hat{{\boldsymbol{r}}}({\boldsymbol{x}})$. Based on the notation in Lee et al. (2008), the ED correlation function is given by

$$\begin{eqnarray}&&{\omega }_{\delta }(r)=\langle | {\hat{e}}_{a}({\boldsymbol{x}})\cdot \hat{{\boldsymbol{r}}}({\boldsymbol{x}}){| }^{2}\rangle -\displaystyle \frac{1}{3},\end{eqnarray} \tag{ 11 }$$

which is zero for randomly oriented galaxies in a uniform distribution. In the simulations the matter density field is traced using the positions of dark matter particles.

The projected shape correlation function, ${w}_{\delta +}$ is directly related to the correlation function measured in observations. Following the notation of Mandelbaum et al. (2006), we define the the matter–intrinsic shear correlation function ${\hat{\xi }}_{\delta +}({r}_{p},{\rm{\Pi }})$ and the corresponding projected two-point statistic ${w}_{\delta +}$. In this paper, r_p is the comoving transverse separation of a pair of galaxies in the XY plane and Π is their separation along the Z direction.

The components of the projected ellipticities of a galaxy are given by

$$\begin{eqnarray}&&({e}_{+},{e}_{\times })=\displaystyle \frac{1-{(b^{\prime} /a^{\prime} )}^{2}}{1+{(b^{\prime} /a^{\prime} )}^{2}}[\cos (2\phi ),\sin (2\phi )],\end{eqnarray} \tag{ 12 }$$

where b'/a' is the axis ratio of the projected shape of the stellar component of a galaxy, and ϕ is the position angle of the major axis with respect to the reference direction (position of the dark matter particle). Here, e₊ refers to the radial component and e_× is the component rotated at 45°. The matter–intrinsic shear correlation function is given by,

$$\begin{eqnarray}&&{\hat{\xi }}_{\delta +}({r}_{p},{\rm{\Pi }})=\displaystyle \frac{{S}_{+}D}{{RR}}\end{eqnarray} \tag{ 13 }$$

where S₊ represents the "shape sample," selected on the basis of a binning in subhalo mass and the "density sample" labeled by D consists of the dark matter particles used to trace the matter density field. S₊D is given by the following sum over all galaxy–dark matter particle pairs with separations r_p and Π:

$$\begin{eqnarray}&&{S}_{+}D=\displaystyle \sum _{i\ne j| {r}_{p},{\rm{\Pi }}}\,\displaystyle \frac{{e}_{+}(j| i)}{2{ \mathcal R }},\end{eqnarray} \tag{ 14 }$$

where ${e}_{+}(j| i)$ is the + component of the ellipticity of a galaxy (j) from the shape sample relative to the direction of a dark matter particle (i) selected from the density sample. Here, ${ \mathcal R }=(1-{e}_{\mathrm{rms}}^{2})$ is the shear responsivity that converts from distortion to shear (Bernstein & Jarvis 2002), with e_rms being the rms ellipticity per component of the shape sample. The RR term in Equation (13) refers to the expected number of randomly distributed pairs in a particular (r_p, Π) bin around galaxies in the shape sample.

The projected shape correlation function ${w}_{\delta +}({r}_{p})$ is given by

$$\begin{eqnarray}&&{w}_{\delta +}({r}_{p})={\int }_{-{{\rm{\Pi }}}_{\max }}^{+{{\rm{\Pi }}}_{\max }}{\hat{\xi }}_{\delta +}({r}_{p},{\rm{\Pi }})\,d{\rm{\Pi }}.\end{eqnarray} \tag{ 15 }$$

We calculate the matter–intrinsic shear correlation function over the whole length of the box, L_box with Π_max = L_box/2, where the length of the box is 100 $\,{h}^{-1}\,\mathrm{Mpc}$ or 25 $\,{h}^{-1}\,\mathrm{Mpc}$. The projected correlation functions are obtained via direct summation.

To study the effects of modifying baryonic feedback parameters, we use small-volume simulations, as larger-simulation volumes would not be feasible at present. Smaller-volume simulations, however, will be a subject to larger cosmic variance, and so may be biased relative to the larger box.

In order to estimate the error we are going to incur by using smaller boxes, we use the DC mode formalism (Sirko 2005) that allows one to approximately quantify the effect of the missing large-scale power. Ideally, one would need to run a whole ensemble of the simulations with randomly chosen DC modes. However, due to limited computational resources, we only perform three independent realizations of the 25 h⁻¹ Mpc box with the amplitude of DC mode set to zero and to ±Δ₀, where Δ₀ is the rms density fluctuations in a cubic 25 h⁻¹ Mpc box at z = 0. For any of our statistical measures we then can use the spread between the three realizations as an, admittedly crude, estimate of the uncertainty due to the limited simulation volume.

For the WMAP7 cosmological parameters and the box size of 25 h⁻¹ Mpc box size Δ₀ = 0.585. In the most general case, accounting for the DC mode requires modifications to the simulation code. However, Sirko (2005) showed that for a cosmology that includes only matter, the cosmological constant and, optionally, curvature, the DC mode can be accounted for by a simple rescaling of the cosmological parameters. In this paper we use such a rescaling to include the DC mode in P-Gadget that does not support the DC mode explicitly.

4.1. Distribution of Shapes and Misalignment Angles

In Figure 1 we show a comparison between the cumulative distribution functions for the shapes, b/a and c/a in two mass bins for the original 100 h⁻¹ Mpc MBII-100 run and our three 25 h⁻¹ Mpc MBII-25 simulations with different DC modes. The mean values for the shapes are tabulated in Table 1. Because of the limited size of our simulation volumes, we are only able to consider two mass bins. However, this may be sufficient to notice a really strong trend with halo mass; more subtle trends are missed by us and will have to be explored in the future with more precise simulations. Throughout this paper, the galaxy shapes and alignments are analyzed at z = 0.3.

Zoom In Zoom Out Reset image size

Table 1.  Mean 3D Shapes b/a and c/a of the Stellar Component of the Fiducial MBII-100 Simulation and Three MBII-25 Simulations with Different DC Modes of 0 and ±1σ

	MBII-100		MBII-25, ΔDC = 0		MBII-25, ΔDC = +1σ		MBII-25, ΔDC = −1σ
Msubhalo (${h}^{-1}\,{M}_{\odot }$)	$\langle b/a\rangle$	$\langle c/a\rangle$	$\langle b/a\rangle$	$\langle c/a\rangle$	$\langle b/a\rangle$	$\langle c/a\rangle$	$\langle b/a\rangle$	$\langle c/a\rangle$
109.5–1012.0	0.79 ± 0.0	0.61 ± 0.0	0.78 ± 0.0	0.60 ± 0.0	0.75 ± 0.0	0.56 ± 0.0	0.795 ± 0.003	0.62 ± 0.0
1012.0–1015.0	0.74 ± 0.0	0.525 ± 0.002	0.76 ± 0.02	0.524 ± 0.015	0.73 ± 0.02	0.515 ± 0.013	0.77 ± 0.02	0.51 ± 0.015

Download table as:  ASCIITypeset image

For the three smaller-volume simulations we can compute both the mean over the three realizations, and the error in that mean, which we show in these and all subsequent figures with lines and bands respectively. Since the small box simulations may be biased and/or insufficiently accurate, we use the error in the mean as the estimate of our theoretical error due to the limited box size. For example, from Figure 1 it is clear that the differences between the mean of three MBII-25 runs and the original MBII-100 run are comparable to the error on MBII-25, and that error is reasonably modest, about 2%. Hence, by using smaller boxes we do introduce a bias, but this is modest and is comparable to the statistical error of the simulation results.

The distributions of misalignment angles in the same two mass bins are shown in Figure 2, and their mean values are given in Table 2. We find that the galaxies in the lower-mass bin of smaller-volume simulations are more aligned, at about the 2σ level, than in the fiducial MBII-100 run, and in the high-mass bins the low abundance of halos becomes appreciable. In both cases, however, the bias in using smaller boxes is still sufficiently modest (less than 3°) to justify our use of smaller boxes in this first, exploratory work.

Zoom In Zoom Out Reset image size

Table 2.  Mean 3D Misalignment Angles, $\langle \theta \rangle$ (Degrees), between the Major Axis of Galaxies and Their Host Dark Matter Subhalos in the MBII Simulation of 100 h⁻¹ Mpc Size Box and Simulations of 25 h⁻¹ Mpc Box with DC-modes: 0, ±1σ

Msubhalo (${h}^{-1}\,{M}_{\odot }$)	MBII-100	MBII-25, ΔDC = 0	MBII-25, ΔDC = +1σ	MBII-25, ΔDC = −1σ
109.5–1012.0	33259 ± 0078	31151 ± 0603	31761 ± 0662	28153 ± 0545
1012.0–1015.0	27157 ± 0409	31785 ± 3525	25505 ± 2701	2649 ± 435

Download table as:  ASCIITypeset image

The two-point statistics ED and ${w}_{\delta +}$ are shown in Figure 3. The ED and ${w}_{\delta +}$ correlation functions in the MBII-100 simulation and in the mean of MBII-25 runs are in good agreement on small scales and in the high-mass bin. The agreement is worse at large scales in the low-mass bin, but the measurements there are also noisy. The formal error on the mean of three MBII-25 runs is smaller than the difference between the two box sizes, but since the error is estimated from just three runs, it may itself be inaccurate. For the MBII-100 run, we do not calculate the error from independent runs with DC mode incorporated. We can estimate the error due to missing power on large separations to be roughly half of the error in the smaller-volume simulation.

Zoom In Zoom Out Reset image size

Overall, we find that our 25 h⁻¹ Mpc boxes are a suitable, albeit not ideal and moderately biased, tool for exploring the sensitivity of the simulation predictions to the parameters of the star formation and feedback model.

In this section, we explore the effects of modifying the feedback parameters in the simulation on the galaxy shapes and two-point statistics. We follow the methodology of the previous section, and use the three MBII-25 runs with different DC models as our new fiducial simulation set against which we compare runs with varied physics. The details about which parameter is varied in a given model are provided in Section 2.4.

The cumulative shape distributions are plotted in Figure 4 in two mass bins of 10^9.5–12.0 h⁻¹ M_⊙ and 10^12.0–15.0 h⁻¹ M_⊙. Comparing the distributions and the mean values shown in Table 3, we that in the lower-mass bin the axis ratio b/a is larger for the simulation with weaker wind feedback, although the effect is not large, within 2σ of the fiducial model—a deviation comparable to the difference between the mean of three MBII-25 runs and the original MBII-100 run.

Zoom In Zoom Out Reset image size

Table 3.  Mean of b/a, c/a, and θ of the Stellar Shape of Galaxies for Simulations with Varying Star Formation Feedback

	109.5–1012.0(${h}^{-1}\,{M}_{\odot }$)			1012.0–1015.0(${h}^{-1}\,{M}_{\odot }$)
Simulation	b/a	c/a	θ	b/a	c/a	θ
MBII-25	0.779 ± 0.004	0.603 ± 0.003	31.151 ± 0.603	0.76 ± 0.0195	0.524 ± 0.015	31.785 ± 3.525
tsfr-High	0.783 ± 0.003	0.618 ± 0.003	31.883 ± 0.557	0.777 ± 0.02	0.532 ± 0.014	32.512 ± 3.490
tsfr-Low	0.775 ± 0.004	0.582 ± 0.003	33.306 ± 0.648	0.767 ± 0.017	0.533 ± 0.013	30.747 ± 2.973
AGN-High	0.773 ± 0.004	0.604 ± 0.003	30.265 ± 0.600	0.73 ± 0.02	0.513 ± 0.013	29.725 ± 3.18
Wind-High	0.793 ± 0.002	0.636 ± 0.002	37.783 ± 0.427	0.747 ± 0.02	0.535 ± 0.014	32.512 ± 3.49

Download table as:  ASCIITypeset image

Despite all deviations being moderate and not highly significant, some trends are nevertheless intriguing. For example, in the low-mass bin weaker feedback makes galaxies rounder, while in the high-mass bin the (mild) deviation is in the opposite direction.

The cumulative distributions of misalignment angles are shown in Figure 5. The effect of the lower wind loading factor is larger on the angles than on the shapes for lower-mass galaxies—since the wind carries away linear and angular momenta, it can directly affect the orientation of the stellar distribution without affecting the shape too much. For more massive galaxies the effect disappears, however, as in that mass bin the feedback is dominated by AGNs. One can hypothesize that AGNs, being centrally located, are not able to eject large amounts of angular momentum.

Zoom In Zoom Out Reset image size

If this interpretation of our findings is valid, then the critical quantity that controls the distributions of shapes and angles is the angular momentum of the wind; once simulations get it right, their predictions for intrinsic alignment become robust and accurate.

In a previous study, Velliscig et al. (2015) explored the variation in the ellipticities and misalignments, compared with their fiducial model of EAGLE simulation for three different feedback implementations. They investigated models with weaker and stronger stellar feedback and no AGN feedback. They also found that shapes are affected much less than angles, consistent with our hypothesis above. However, they found a larger effect of the stellar feedback on misalignment angles in more massive (>10¹² M_⊙) galaxies, while their measurements for lower mass galaxies were too noisy to be conclusive.

Overall, however, we find a good agreement with the EAGLE simulations, which is encouraging, but not particularly surprising—modern simulations reproduce many observed properties of galaxies fairly.

Two-point statistics for simulations with varied physics are shown in Figure 6. Differences between various models are similar to the level of difference between the MBII-100 and MBII-25 runs: correlation functions agree well on small scales and in the high-mass bin, but exhibit significant variations on scales between 0.1 h⁻¹ Mpc and 1 h⁻¹ Mpc. These variations are non-monotonic and unsystematic, and are likely caused by the lack of statistics in our small box runs. However, just as in previous statistics, we find the largest difference in the run with the low wind mass loading factor.

Zoom In Zoom Out Reset image size

In particular, the dip in the w_δ correlation function at ∼0.3 h⁻¹ Mpc appears to be real—it is insensitive to the numerical details of computing the correlation function such as binning, sample selection, etc. The dip is located close to the radius where the one-halo term transitions to the two-halo term, and may reflect physical processes occurring at the halo–IGM interface. Unfortunately, our simulation volumes are too small to make any strongly statistically significant claim. Using large-volume simulations with similar resolution, we can get a better understanding of these scale-dependent differences and their impact on modeling the intrinsic alignment contamination.

Our primary goal in this paper is to explore the effects of model parameters in the star formation and feedback models on galaxy shapes and alignments using small-volume simulations of size 25 h⁻¹ Mpc on a side. As our fiducial model for the simulation, we adopted the same star formation and feedback model as in the MassiveBlack-II hydrodynamic simulation of galaxy formation (Khandai et al. 2015), which is performed in a box of volume (100 h⁻¹ Mpc)³.

Simulations with significantly (by factors of of 3–10) varying feedback show remarkable consistency with the fiducial run. Within the statistical precision we are able to achieve in our small-volume runs, most of the observational probes are insensitive to the details of subgrid physical modeling, with the exception of misalignment angles. We hypothesize that the angular momentum ejected by galactic winds is the most crucial physical quantity that determines the alignment of stellar shapes, and it remains one of the least robust quantities predicted in modern simulations of galaxy formation. It is likely that the differences in numerical implementations of hydrodynamics can also cause differences in shapes and alignments. While it is known that variations in feedback are more important than numerical implementations when comparing the properties of the simulated galaxies with observations (Scannapieco et al. 2012), a similar analysis has not been performed for galaxy shapes and alignments. This will be explored in a future study.

Our conclusions are also in good agreement with a similar exploration of the role of subgrid physics on intrinsic alignments by the EAGLE simulation team.

Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the United States Department of Energy. AT is supported by the Fermilab Graduate Student Research Program in Theoretical Physics. Simulations have been performed on National Energy Research Supercomputing Center (NERSC) supercomputers "Cori" and "Edison." A.T. thanks Nishikanta Khandai and Tiziana DiMatteo also for sharing their P-Gadget code used to run the simulations in this work.