AN EMPIRICAL MODEL FOR INTRINSIC ALIGNMENTS: INSIGHTS FROM COSMOLOGICAL SIMULATIONS

,


INTRODUCTION
Weak lensing is an important probe for cosmological galaxy imaging surveys, and will be of major interest for the Rubin Observatory's Legacy Survey of Space and Time (LSST, Ivezić et al. 2019) as well as other Stage IV surveys such as Euclid (Scaramella et al. 2022) and Roman (Akeson et al. 2019).These new surveys will be much more powerful than those in the past, and will require correspondingly more precise analyses.The intrinsic alignments (IA) exhibited by galaxies as a result of their alignment with large-scale structure and other nearby galaxies can contaminate cosmic shear measurements, driving the need for accurate IA modeling (Heymans et al. 2006;Blazek et al. 2011;Joachimi et al. 2013;Krause et al. 2015).
IA has typically been modeled using analytic approaches, usually motivated through cosmological perturbation theory (e.g.Hirata & Seljak 2004;Bridle & King 2007;Blazek et al. 2015Blazek et al. , 2019;;Vlah et al. 2020Vlah et al. , 2021;;Maion et al. 2023;Bakx et al. 2023; Chen & Kok-⋆ E-mail: vanalfen.n@northeastern.eduron 2023).However, these approaches are unable to describe alignments in the fully nonlinear regime.To address this issue, in recent years, interest in halo-based models of the intrinsic alignment of galaxies has taken off, with Fortuna et al. (2020) updating previous work by Schneider & Bridle (2010) to provide a robust theoretical description of intrinsic alignment within dark matter halos.At the same time, efforts have emerged to infuse dark-matter-only N-body simulations with mock galaxies which have correlated alignments.These methods include halo-based modeling (e.g.Joachimi et al. 2013;Hoffmann et al. 2022), tidal field descriptions (e.g.Harnois-Déraps et al. 2022), and machine learning approaches (e.g.Jagvaral et al. 2022).These simulation capabilities are important for validating our approaches to mitigate IA effects in weak lensing surveys, providing realistically complex mock data on which to test our analysis pipelines.These simulations can also be used to directly model IA, including in situations where analytic modeling is not feasible, such as higher-order or map-based analyses.Finally, the resulting mock galaxy catalogs can also help us in understanding and testing galaxy formation and evolution models which could impact the large-scale alignment of galaxies; either directly through comparison and calibration to hydro sims, or through interpreting the resulting behavior of different galaxy formation or evolution models in terms of the HOD IA parameterization.This mock generation process provides the connection between small scale galaxy formation physics either as calculated by a hydro sim or through some other analytic method and IA observables on a wide range of scales.In other words, the parameters of this HOD model are informed by the small-scale galaxy formation physics and allow us to make predictions, including statistical uncertainty, for a wide range of IA observables (Hoffmann et al. 2022;Samuroff et al. 2023).
In this work, we provide a new, flexible, and efficient halo-based method for including IA in mock galaxy catalogs.This method is built on the halotools1 package for generating mock galaxy catalogs (Hearin et al. 2017) and can be applied within larger simulation and analysis frameworks to produce mock galaxy catalogs or for IA modeling (Van Alfen et al. 2024).
Models of the Halo Occupation Distribution (HOD) typically have components to capture the number of galaxies within halos (occupation), as well as the distribution of their positions and velocities within halos (phase space).In this work, we add two new components to this framework: one component to model the statistical distribution of galaxy orientations within their halos, and a second component, the IA strength model, that characterizes the probabilistic alignment between the galaxy orientation and some reference vector associated with the parent halo.These new ingredients to HOD modeling are quite flexible and allow for multivariate dependence on galaxy properties.
This paper is organised as follows: Section 2 presents background on two-point clustering, giving the forms of the correlation functions used to analyze our models.Section 3 describes the cosmologies of the simulations used in this work and discusses how we determine our covariance.In Section 4 we introduce the Dimroth-Watson distribution that we use to characterize the distribution of galaxy orientations within their halos, and we describe the probabilistic models of central and satellite galaxy alignments that will be used throughout the rest of the paper.In Section 5, we examine the impact of altering central and satellite alignment strengths.This provides a baseline for understanding the simplest effect our alignment strengths have on our model.In Section 6, we use high-resolution N-body simulations to investigate the accuracy with which subhalo-orientation correlation functions can be captured by our HOD-type models using radial alignments, as subhalo orientations are often not available.We use the IllustrisTNG hydrodynamical simulation to test the flexibility of our model of alignment strength in Section 7. Section 8 summarizes our results and future work.

TWO-POINT CORRELATION FUNCTIONS AND THE HALO OCCUPATION DISTRIBUTION
In this section, we introduce the basic equations of the analytical model of halo occupation statistics.His-

Theory: The Halo Model and correlation functions
In this paper, we consider the halo model (Cooray & Sheth 2002;Asgari et al. 2023), which provides a flexible approach for us to use in our creation of mock galaxy catalogs.Here we overview the theoretical equations for the two-point correlation functions.The functions are given to provide a conceptual background, but our analysis does not use these formulae to calculate correlation functions.Instead, we measure the correlations from mock data using the estimators described in Section 2.2.
Below, we see the theoretical formulation of the position-position correlations.
ξ gg (r) = ξ 1h gg (r) + ξ 2h gg (r) (1) ξ 2h gg (r) = ξ 2h cc (r) + ξ 2h cs (r) + ξ 2h ss (r) (3) where dn dM is the halo mass function (the differential distribution of number density of halos n as a function of halo mass M ), ⟨N g |M ⟩ is the expectation value of the number of galaxies in a halo of mass M , and ng is the overall number density of the population of galaxies being considered, and we have ignored the central and satellite fractions, which are often absorbed into the definition of the correlation function.In the above equations, subscripts indicate the two populations being correlated.The subscript "g" refers to all galaxies, "h" refers to halo, "c" refers to central galaxy (the most massive galaxy at/near the center of the host dark matter halo), and "s" refers to satellite galaxy (a smaller galaxy in a host dark matter halo besides the central galaxy).The superscripts "1h" and "2h" denote the scale contributing to the correlations.Namely, "1h" refers to the 1-halo regime (galaxies within the same halo), "2h" refers to the 2-halo regime (galaxies residing in separate halos), and no superscript indicates correlations from all scales.In this paper, we use 1 Mpc/h as the approximate delineation between the 1-halo and 2-halo regimes.
The assumptions underlying Eq. ( 4) are: • r ≫ R max vir • On very large scales, ξ hh (r) varies slowly across the length scale of halos where R max vir refers to the largest virial radius of all halos being considered for the 2-halo contribution to the correlation.These assumptions together allow equation ( 4) to treat halos as a single point with a given number density of galaxies concentrated at that point.Clearly, while these assumptions work at large scales, they break down as separations approach the 1-halo regime.At such scales, equation (4) no longer accurately represents the correlation function.The breakdown of this assumption is an example of where simulation-based methods are useful.
In Eq. ( 5), λ(x|M ) is the galaxy profile, i.e., the unitnormalized spatial distribution of galaxies within a halo of mass M .This expression includes both centralsatellite and satellite-satellite pairs.Above, we looked at the position-position correlation functions.Here, we overview the orientation-orientation and position-orientation correlation functions.The equations presented here are analogous to those described in more detail in Fortuna et al. (2020).Following typical naming conventions (e.g. Lee et al. 2008;Tenneti et al. 2015), we use the subscripts "e" and "d" to refer to ellipticity and direction respectively, where direction refers to the separation between galaxies, which in some contexts we call "position."As discussed below, we do not use full shape information in this work, so our use of "ellipticity" refers in practice to orientation.
In Eq. ( 6), e is the three-dimensional ellipticity or orientation (i.e. the unit-normalized ellipticity) of the galaxy.Similarly, the 1-halo part of the ee correlation function and both parts of the ed correlation function are and For ξ ee (r), we can decompose the final term in the integrand into central and satellite contributions: where f cen,i and f sat,i are the fraction of galaxies in sample i that are centrals and satellites, so that f cen,i + f sat,i = 1.
For ξ ed (r), this central-satellite decomposition takes the form: where r = r • r refers to the position vector between two galaxies.

Two-point Correlation Functions Estimators
Below are the estimators for the position-position, position-orientation, and orientation-orientation correlation functions that we use to analyze our model.While earlier we showed the theoretical formulation, here we focus on the estimators we use on mock galaxy catalogs.As stated earlier, to remain consistent with the names used in the literature, we have described these functions using "ellipticity" rather than "orientation."However, while Eqs.( 6)-( 9) can use full ellipticity, the rest of our paper and the estimators below only use the galaxy orientation.Our method, as presented here, does not include full galaxy shape information, which is the topic of ongoing work.We also note that the functions utilized here are convenient for analysing simulations where three-dimensional position and orientation information is available.However, for analysis of observed data, where only projected shape information is present, other correlation functions are typically used (e.g.Mandelbaum et al. 2006;Singh et al. 2023).
The galaxy-galaxy two-point correlation function is given by: where DD(r) is the number of galaxy pairs separated by r and RR(r) is the expected number of pairs for a random distribution.We use this estimator rather than the Landy-Szalay estimator (Landy & Szalay 1993) even though it is somewhat sub-optimal in some cases (Singh et al. 2017).This estimator is faster, and the periodic nature of our box allows us to use analytic randoms, negating much of the sub-optimality.
The ellipticity-direction (ED) correlation function is defined as: The ellipticity-ellipticity (EE) correlation function is defined as: where x is the position vector of a given galaxy, r is the separation vector between two galaxies, and r is the unit vector of the separation vector r.In both of these equations, subtracting the factor of 1/3 accounts for the fact that integrating the dot product of the two vectors (i.e.integrating the cos2 θ between them) over a sphere results in a factor of 1/3 when the vectors are randomly distributed (i.e.there is no correlation).These equations are examined in more depth in Chisari et al. (2016) with discussion about the factor of 1/3 in Hopkins et al. (2005).
Beyond this section, the symbols for equations ( 12)-( 14) will be used.
Namely we use ξ for the position-position correlation function, ω for the position-orientation correlation function, and η for the orientation-orientation correlation function.We do this to separate the rest of our work from the pedagogical discussion here as well as to connect the measurements discussed later with the estimators used to obtain them.

SIMULATIONS AND COVARIANCE
We construct our HOD models using halotools, a library for simulation-based predictions of the galaxyhalo connection.Catalogs of dark matter halos are the root data product underlying halotools-based predictions, and throughout the paper we make use of catalogs provided by the library, although we note that halotools provides full support for user-supplied simulation data.Upon ingesting a catalog of simulated halos, halotools populates each halo with a synthetic galaxy population according to an occupation model, distributing galaxies within each halo according to an assumed phase space model.Orientations are bestowed upon each synthetic galaxy according to the alignment models outlined in Section 4.
Throughout the paper, we will distinguish between a host halo, which is a distinct gravitationally self-bound collection of simulated dark matter particles, and a subhalo, a gravitationally self-bound object that is located inside the boundary of a larger host halo.We will use the term (sub)halo whenever referring to a collection of host halos and subhalos together and/or when referring to an individual object that may be either a subhalo or a host halo.
We have previously defined central galaxies as the most massive galaxy in a dark matter halo, located at or near the center of the halo.Satellite galaxies are all galaxies in a halo other than the central galaxy.The ratio of central galaxies to satellite galaxies depends on the occupation parameters of the HOD model.As a rule of thumb, the galaxy tables produced with halotools using occupation parameters typical of those this work tend to contain roughly three central galaxies to every one satellite galaxy.

Simulations
In this paper, unless otherwise stated, whenever we create a catalog using halotools, we do so using the available halo catalogs from the Bolshoi-Planck 2 (BolPlanck) simulation output at z = 0, the cosmology and simulation parameters of which are given in Table 1 (Klypin et al. 2011).
While most of the paper uses BolPlanck, Figures 3 and  4 use the Small MultiDark Planck (SMDPL) simulation with cosmology described in Table 1.The simulation is evolved by solving for gravitational interactions only using the L-GADGET-2 code, a version of the publicly available cosmological code GADGET-2 (Springel 2005) with a force resolution of 1.5 h −1 kpc.This simulation belongs to the series of MultiDark simulations with Planck cosmology.More details for this simulation are described in Klypin et al. (2016).
For both BolPlanck and SMDPL, (sub)halos are found using the phase-space halo finder ROCKSTAR (Behroozi et al. 2013a), which uses adaptive, hierarchical refinement of friends-of-friends groups in six phase-space dimensions and one time dimension, and tracked over time using the Consistent Trees algorithm (Behroozi et al. 2013b).As demonstrated in Knebe et al. (2011Knebe et al. ( , 2013)), this results in a very robust tracking of (sub)halos (also see Jiang & van den Bosch 2016).
The essential orientation properties used in this paper (the x, y, and z values of the unit vector for the major axis of the ellipsoid) are also generated with the ROCKSTAR halo finder, shown in Appendix B of Rodríguez-Puebla et al. (2016).For a more detailed look at shape finding, see Appendix A.
ROCKSTAR determines virial spherical overdensity (SO hereafter) volumes for each halo, centered on a local density peak (possible for even very ellipsoidal halos), such that the average density inside the sphere is ρh (z) = ∆ vir (z)ρ m (z).Here ρ m (z) = Ω m (z)ρ crit (z), where ρ crit (z) = 3H(z)2/8πG is the critical energy density of the Universe, and ∆ vir (z) is given by a fitting function (Bryan & Norman 1998): For the Planck cosmology, ∆ vir (z) ≃ 360.The radius of each such sphere defines the virial radius R vir of the halo, which is related to the mass of the halo via M vir = (4/3)πR3 vir ρh .Additionally, subhalos in this catalog are distinct, self-bound structures whose centers are found within the virial radius of a more massive host halo.
In Section 7, we compare results from our HOD models (built using the BolPlanck catalogs described above) to the Next Generation Illustris Simulations 3 (Il- lustrisTNG), a suite of hydrodynamical simulations of galaxy formation in cosmological volumes (Nelson et al. 2017;Pillepich et al. 2017;Naiman et al. 2018;Springel et al. 2017;Marinacci et al. 2018).Illustris resolves galaxy shapes, therefore the IA seen in Illustris is modeled directly by hydrodynamic interactions.The IllustrisTNG runs used in this work are two uniform mass resolution cosmological volume simulations with side lengths 205 h −1 Mpc, one "full physics" run, TNG300, including all of the complex physics of galaxy formation, and a gravity-only counterpart, TNG300-Dark.The initial conditions of the simulations were set at z = 127 using the Zeldovich approximation.The adopted cosmological parameters are given in Table 1.
It is unclear exactly to what extent the cosmologies of the underlying halo catalogs shown in Table 1 will affect resulting galaxy tables.However, the methods described here are highly tunable and can produce a wide range of galaxy populations and alignments from a single halo catalog as shown in Section 7.

Covariance
Building our HOD model for IA is inherently stochastic, as both the selection of halos being populated and the alignment of galaxies (discussed below) are drawn from distributions.When fitting this model to "data" (which in this work include both measurements of the simulated halos being populated as well as galaxies in another simulation such as IllustrisTNG), we must account for the fact that both the data and model have stochasticity.When fitting to the halo measurements (Section 6), we find that calculating the covariance using multiple realizations of the galaxy alignments, without changing the halo occupation, is sufficient to allow a test of the model performance while avoiding double counting the variance from the occupation step.When fitting to IllustrisTNG (Section 7), we calculate the covariance from IllustrisTNG and average over multiple model realizations to remove model variance.Further details concerning our covariance are given in Appendix D.

MODELING ALIGNMENTS
In this paper, we make the assumption that both galaxies and (sub)halos can be modeled as triaxial homologous ellipsoids (see Appendix A for details).The orientation of galaxies/halos can then be entirely described by specifying the vectors associated with the minor, intermediate, and semi-major axes of ellipsoids.Here we focus on modeling the orientation of a single axis of the ellipsoid describing a galaxy within its host halo, but we note that the following framework could be extended to model a three-dimensional galaxy orientation and shape, albeit with increased complexity.Our strategy will be to statistically model the misalignment angle between a vector specifying the orientation of each galaxy's (sub)halo, and a vector specifying the orientation of the galaxy, as shown in Figure 1.Mathematically, if we assume that the vector specifying the orientation of each galaxy's host (sub)halo is the unit vector, ẑ, we can model the orientation vector of the galaxy by specifying the angular coordinates on the unit-sphere.
For this purpose, we utilize the Dimroth-Watson distribution, which specifies the spherical polar coordinates of a unit vector.We chose the Dimroth-Watson distribution as it provides a maximum entropy distribution on a sphere while accounting for the 180 degree symmetry we see in our galaxy orientations (Watson 1965).These criteria are difficult to meet with other distributions, but the Dimroth-Watson distribution naturally fits the role for directional, axially symmetric systems on a sphere (Sra & Karp 2013).The probability distribution for the polar angle, θ, and the azimuthal angle, ϕ, are given by where the normalization factor is given by: Note that the azimuthal angle, ϕ, is modeled as a uniform distribution, so that the angle θ is then by definition the misalignment angle between galaxy and host (sub)halo orientation, θ MA .The degree to which galaxies and (sub)halos align is then controlled by the κ parameter.
It is convenient to re-parameterize the strength of alignment κ by: such that µ = 1, 0, −1 corresponds to perfect alignment, random alignment, and perfect anti-alignment, respectively.As illustrated in Figure 2, a positive value of µ corresponds to a preferential parallel alignment between a galaxy and its host (sub)halo.A negative value of µ corresponds to a preferential perpendicular alignment between a galaxy and its host halo.The distribution of misalignment angles is symmetric about cos(θ MA ) = 0, in accord with the symmetries we have assumed of the shapes of galaxies and halos.There is thus no difference in the morphological orientations between a galaxy rotated by 180 • and its original orientation.In principle, for disk galaxies whose orientation is determined by angular momentum, it may be interesting to distinguish between parallel and anti-parallel angular momentum vectors.Although the difference should  (bottom).In both cases, the misalignment angle is determined with respect to some reference axis, which differs for satellite galaxies depending on the model used, either using subhalo alignment (top), or radial alignment (bottom).
not impact the observable galaxy ellipticities, correlations between the angular momentum direction and the density field can allow us to probe tidal torque theory (e.g. Lee 2019) -we leave exploration of this and related effects to future work.

Galaxy Alignment Models
For all investigations, we align the major axes of central galaxies with the major axes of their host halos.Dimroth-Watson distribution model, equation ( 16), for the distribution of galaxy-halo orientation misalignment angles, θ.Different color lines show distributions with different alignment strengths, µ, given by equation ( 18).We see that as µ approaches 1, the Dimroth-Watson distribution approaches delta functions at cos θ = 1 and −1 (perfect alignment).Meanwhile, as µ approaches −1, the distribution approaches a delta function at cos θ = 0 (perfect antialignment).Notably, the value of µ does not deterministically assign misalignment angle (except in the cases of perfect alignment and perfect anti-alignment), but rather it determines the shape of the distribution from which misalignment angles will be drawn.
For satellite galaxies, we explore two models in this paper for determining their orientations: 1. satellites are oriented relative to that of their host dark matter subhalo, called subhalo alignment (as seen in the upper image of Figure 1); 2. satellites are oriented relative to the host halo centric radial vector, called radial alignment (illustrated in the lower image of Figure 1).In this case, there are two different alignment strengths explored.
(a) Constant alignment strength, where the alignment strength is the same for all satellite galaxies.
(b) Distance-dependent alignment strength, where the alignment strength of a satellite depends on the distance to its central galaxy.
In the case of subhalo alignment, the model needs information about subhalo positions and orientations, while the radial alignment model may be used with or without subhalo positions and needs no information about subhalo orientations (though we typically use subhalo positions unless otherwise stated).

Building an HOD Mock
To build an HOD-based mock catalog, we use the HODModelFactory in halotools.Central and satellite galaxies have separately-defined occupation models that determine the number of galaxies of each population present in a given halo.Additionally, each population has its own model for the intra-halo phase space distribution.Throughout the paper, central galaxies are placed at the center of the host halo with the same peculiar velocity as the halo.For satellite galaxies, we consider two different classes of phase space models.The first is defined by the NFWPhaseSpace class: the spatial distribution of satellites is spherically symmetric and has a radial density profile given by the NFW distribution (Navarro et al. 1997); satellite velocities are isotropic and exhibit a velocity dispersion profile determined by the Jeans equation.The second satellite model we consider is given by the SubhaloPhaseSpace class, in which the position and velocity of each satellite is defined by a subhalo that has been randomly selected from the simulated halo, preferentially selecting the most massive subhalos first; when N sat is smaller than the number of subhalos in the host halo, the lowest-mass subhalos are left unoccupied, and when N sat exceeds the number of available subhalos, the remaining satellites are distributed according to the NFWPhaseSpace class.The SubhaloPhaseSpace class provides a way for the synthetic satellite population to inherit the complexity of the intra-halo distributions of simulated subhalos, while at the same time retaining the flexibility of the HOD to occupy host halos with a variable number of satellites.

IMPACT OF GALAXY-HALO MISALIGNMENT ON ORIENTATION CORRELATION FUNCTIONS
In this section, we use our halo model of IA to provide a pedagogical investigation of how galaxy-halo misalignment manifests in orientation correlation functions.We will separately study effects arising from central and satellite galaxies, and we will consider the full range of spatial scales, spanning both 1-halo and 2-halo regimes.Our general strategy for this investigation is to begin by orienting galaxies to be perfectly aligned with the reference vector of the (sub)halo; we then leverage the flexibility of our probabilistic model by programmatically weakening the alignment strength, and studying how the orientation correlation function changes as a result.

Central and Satellite Galaxy Misalignment
We begin by examining the effect of varying the central galaxy -host halo alignment strength on the orientation correlation functions.In Figure 3, we show the positionorientation (ω) on the left and orientation-orientation (η) correlation functions on the right.In both plots, we show subhalos and host halos with a peak mass greater than 10 12 h −1 M ⊙ as black points with error bars.Errors on these measurements are estimated using jackknife resampling of the simulation box, by splitting into 4 3 equal volume sub-samples.Also shown in Figure 3 are color-coded lines, showing ω g gg and η g gg (position-orientation and orientationorientation correlations for all galaxies) for model galaxies in the same sample of (sub)halos with varying levels of central-host halo major axis alignments.Here we use the superscripts "g" and "h" to specify the correlation is being done on the galaxies or halos respectively (e.g.ω h cs refers to the position-orientation correlation between the positions of halos hosting central galaxies and the orientations of subhalos hosting satellite galaxies).The satellite-subhalo alignment strength is kept fixed at full strength, µ sat = 1.Because there is still significant uncertainty in the η measurement in a 400 Mpc box, we use a fitting function to smooth the results (see Appendix B for details).For maximum alignment, µ cen = 1, the correlation functions reproduce the (sub)halo measurement by construction.
Next we show the effect of varying the satellite-subhalo alignment strength on the orientation correlation functions.Similar to previous section, in Figure 4 we show the ω on the left and η on the right, with correlation functions for subhalos and host halos given as black points with error bars.In Figure 4, the color-coded lines correspond to varying levels of satellite-subhalo alignments, while the central-halo alignment strength is kept fixed at full strength, µ cen = 1.Again, when µ sat = 1, the (sub)halo orientation correlation functions are reproduced by construction.
Broadly speaking, we can see that central galaxy alignments produce a much larger effect on the correlation functions relative to satellites, especially at larger scales.The satellite correlations in Figure 4 show noticeable effects at smaller scales, but become much less important on scales larger than ∼ 1 Mpc/h.Meanwhile, Figure 3 shows that changing the alignment strength of centrals has a noticeable effect at all scales.

MODELING SATELLITE ORIENTATIONS WITHOUT SUBHALO ORIENTATIONS
Our analyses in the previous sections have relied upon alignments in which the orientation of a synthetic satellite was variably-correlated with the major axis of its parent subhalo.But relying on subhalo catalogs as the core data product of the model introduces a significant source of systematic error arising from uncertainty in subhalo-finding and numerical/resolution issues on small scales (e.g., van den Bosch et al. 2018;Campbell et al. 2018).The resolution requirements on subhalo finding are steep, especially for predictions involving subhalo internal structure (e.g., Rodríguez-Puebla et al. 2016;Benson 2017;Mansfield & Avestruz 2021), and so making theoretical predictions with models that require resolved subhalo orientations leads to costly demands on the Gpcsized simulations that are required to analyze cosmological surveys.The positions of the subhalos themselves are fairly well-resolved and do not introduce notable uncertainty.
These considerations motivate an alternative formulation of satellite orientations with various levels of subhalo information.We start by looking at radial alignments for cases where no subhalo orientation is present, while still using subhalo positions to place satellites.From there, we move on to examine the ability of this simple model to match halo correlations in the 1-halo and 2-halo regimes, justifying our choice to fit our parameters to the 1-halo regime.Following that discussion, we examine the effect of varying levels of satellite anisotropy, including a fully isotropic model that shows the effect of building an HOD with no subhalo information at all.We then conclude this section with a discussion on the goodness of fit of the best radial alignment models as well as by introducing an alternate method of determining a radial alignment strength in the case where subhalo orientations are present.The colored lines are for model galaxies with varying levels of alignment strength between central galaxies and their host halo, from random alignments (light blue) to perfect alignments (pink).In each model, satellite galaxies take on the same orientation as their subhalo.On the right, a fitting function is used to smooth the results, letting us focus on the overall trend from lowering the central alignments which tends to be noisy for ω correlations.We also see a downward bump around the 1-halo to 2-halo scale transition in the right panel.This suggests to us an effect related to the satellite-satellite correlations, with the sudden upturn going into the1-halo regime being related to satellite pairs within a halo becoming more significant.This effect also only becomes apparent once the central alignment strength has become weak enough that the satellite alignments become more important.The colored lines are for model galaxies with varying levels of alignment strength between satellite galaxies and their subhalo, from random alignments (light blue) to perfect alignments (pink).In each model, central galaxies take on the same orientation as their host halo.On the right, a fitting function is used to smooth the results.
In this section we carry out a targeted study of the flexibility of IA models of satellites that do not rely upon subhalo orientations.Our primary assumption will be that the orientation of a satellite galaxy tends to align radially towards the center of the host halo.We examine this assumption in the simplest alignment model in §6.1, and in subsequent subsections we proceed to incorporate additional complexity and study the impact on orientation correlation functions.In order to assess the flexibility required by an HOD-type model of satellite orientations, we begin by focusing on recovering the orientation correlation functions exhibited by (sub)halos at z = 0 in the gravity-only Bolshoi-Planck simulation.Note that we do not claim that satellite galaxies always perfectly align with their respective subhalos, but we use this as a test of how much subhalo orientation information can be captured when such information is not available.In this way, we probe the full flexibility provided by our model.Following this preliminary investigation, in Section 7, we will proceed to apply our IA model to full-physics hydrodynamical simulations (described in Section 3).

Radial Alignment Modeling
Constant radial misalignment uses a single, constant alignment strength for all galaxies, giving them all the same Dimroth-Watson distribution for their radial alignments.A distance-dependent alignment strength changes the shape of the Dimroth-Watson distribution for each satellite galaxy depending on its distance from its central galaxy.We choose a form of this relationship given by a power law -equation ( 19).Because µ ∈ [−1, 1], this equation is limited to that range as well, assigning -1 to all µ < −1, and 1 to all µ > 1.
Note that equation ( 19) is not the only possible model for radial alignment, but it is one possible simple way to model distance-dependent radial alignments.In particular, this model requires that all satellite galaxies have either preferential alignment, or preferential antialignment (or random in the case of a = 0).We have chosen this particular model for its simplicity, but the same concept of assigning individual alignment strengths based on such a model is general.
Since µ refers to how strongly galaxies tend to align with respect to a given reference vector, the µ values in this section refer to how strongly satellite galaxies align with the radial vector between the satellite and the center of its host halo.In previous sections, we used (sub)halo major axes as reference vectors, so matching the halo correlation in that context would mean perfectly aligning with the halo major axes.Here, however, matching the halo correlations is less straightforward.We would only expect a value of µ = 1 in this case if subhalo shapes are perfectly aligned with their radial vectors.
In any given mock galaxy catalog that we generate, there are roughly three times as many central galaxies as there are satellite galaxies (i.e., the satellite fraction is roughly 25% for the galaxy populations we consider), so looking at the correlation functions between all galaxies dilutes the effect of the satellites, and is heavily influenced by the central correlations.In order to compare the effects of these satellite alignments strengths, we look at ω g cs (the correlation functions for the central galaxy position and satellite galaxy orientation) rather than the full ω correlation function for all galaxies to all galaxies.In this way, we can better see the effect of only the satellite orientations.All fits discussed in this section are fits to central position, satellite orientation correlation functions.We do note that in the 2-halo regime, it is expected that the overall alignment correlation functions are dominated by the impact of central-central alignment only (Schneider & Bridle 2010;Fortuna et al. 2020); this is an important caveat to bear in mind in the following discussions in regards to the impact of our model's ability or inability to fit central-satellite correlations well in the 2-halo regime.
As shown in central plot of Figure 5, and discussed in more detail in Section 6.4, the distance-dependent and constant radial alignment strengths perform comparably.The correlation function from the distance-dependent radial alignment may match that of its halos slightly better, but adds an extra degree of complexity.We also see that the correlation functions match much better in the 1-halo regime (here, defined as anything within 1 Mpc/h) than overall.
For the fits shown in Figure 5, we determined parameters by fitting the 1-halo portion of the ω cs correlation functions from the galaxies to those of the halos using MCMC.We discuss our decision to fit solely to the 1-halo regime in Section 6.2, and we note that we use this type of fitting solely to show the flexibility of the model, as the radial alignment model is intended for situations where subhalo orientations are not known.We find the best constant radial alignment strength to be µ sat = 0.826, and the best distance-dependent alignment strength, following equation ( 19), to have a = 0.804 and γ = −0.04.
While radial alignment is intended for a situation in which no subhalo orientations are available, we do have subhalo orienations and can use them for comparison purposes now.Comparing to the subhalo correlation functions, we see in Figure 5 that the distance-dependent alignment strength produces correlation functions that match within 40% of the halo correlation functions overall (much better on average), and well within 15% in the 1-halo regime (much better on average: about 6%).Meanwhile, the constant alignment strength produces results with very similar values.However, in the 1-halo regime, the average percent difference between the constant model and the halos is slightly worse: about 8% off on average.A closer examination of the relative error can be seen in Section 6.2 in the context of fitting to the 1-halo and 2-halo regimes; both the distance-dependent and constant models here closely match those seen in the 1-halo fit model described there.While there may be a slight preference for a distance-dependent radial alignment, the two methods perform comparably.Certainly, constant radial alignment provides a simpler model and may be sufficient for many needs.
In Figure 5 we restricted our view to central-satellite correlations.We have done this because the radial alignment model inherently contains central-position, satellite-orientation information, but we noted earlier that the 2-halo correlation function can be decomposed into central-central, central-satellite, and satellitesatellite portions.In Figure 6 we see a similar de- Here, we see a slight difference between ω g cs (the galaxy correlation functions) and ω h cs (the corresponding halo correlation function) with the distancedependent alignment (solid red line) performing slightly better at small scales than the constant alignment strength (shown as a dashed blue line).Right: ηcs correlation function (central galaxy orientation and satellite galaxy orientation).In all panels (and in figures below), the shaded regions show the 1-sigma confidence interval of each model.Here we see a generally good agreement except at small scales.Right: ωgs correlation function.We see that, although ωss performs poorly at small scales, the ωcs portion dominates making the overall correlation a good fit.
composition of position-orientation correlation functions into central-satellite (ω cs ), satellite-satellite (ω ss ), and all galaxy-satellite (ω gs ) pairs.The (ω cs ) correlation is identical to the constant radial alignment portion of the central panel of Figure 5, but the other panels look at other contributions.As we have seen, (ω cs ) fits the 1halo regime well and slightly underestimates the 2-halo regime.Meanwhile, the (ω ss ) portion performs reasonably well until we reach small scales.In the (ω gs ) panel, we see that, though (ω ss ) fails at small scales, the (ω cs ) portion dominates, making the overall agreement good.
6.2.Fitting to the 1-halo and 2-halo Regimes We have compared the galaxy correlation functions to those of their respective halos for both the full range, as well for the 1-halo regime alone, generally finding that our models fared better on small scales.Because of this, we restrict ourselves to the 1-halo regime when finding the best-fit parameters for the constant and distancedependent radial alignment models.In this section, we further justify that decision by showing the limits of our model when considering satellite alignments.For both the constant radial alignment strength and the distance-dependent radial alignment strength cases, we ran MCMC chains and compared the galaxy correlation functions to the halo correlation functions in the 1-halo regime.To run these chains, and any time we use MCMC, we use the emcee python package (Foreman-Mackey et al. 2013).
Our attempts to fit to the correlation function across all scales led to a wide spread of a and γ values for equation (19), producing extreme values for µ, and resulting in poor fits and no clear best-fit parameter values.Our attempts to fit to large-scale correlations alone (i.e., the 2-halo regime) resulted in similar issues.Meanwhile, when restricting attention to the small-scale correlations alone, our fitter enjoyed higher-accuracy success.These results, and the fact that we underestimate the 2-halo correlations when fitting to the 1-halo regime only, imply that there is extra environmental information needed to reproduce correlation functions between central position and satellite orientation in the 2-halo regime, and that a simple radial alignment does not fully capture subhalo orientations.This situation is not surprising, since our model only directly includes information from the 1-halo regime (e.g.satellite positions and (sub)halo orientations).
In Figure 7, we see a direct comparison of ω g cs (the central galaxy position, satellite galaxy orientation correlation function), as measured directly from simulated (sub)halos, and two different HOD models: one with alignment strengths found by fitting to the 1-halo regime only (labeled as the 1-halo fit model), and one defined by fitting to the 2-halo regime only (which we will call the 2-halo fit model).Overall, the 1-halo fit model performs better; it matches the 1-halo regime well and performs comparably to the 2-halo fit model in the 2-halo regime.Neither fit suggests that our model can accurately capture the full range of spatial scales, but looking at the 1-halo and 2-halo regimes separately, we see the strengths and weaknesses of each model in more detail.
In the 1-halo range of scales, the 1-halo fit model achieves ∼ 12% accuracy at worst, while the 2-halo fit model only matches within 90%.In the 2-halo regime, the 1-halo fit model is only accurate at the 40% level, but even the 2-halo fit model performs no better, differing by almost as much.
For a clearer picture, the χ 2 per degree of freedom (χ 2 dof ) (describing how well the galaxy correlation functions match those of their halos) of the 1-halo fit model is 5.2, 8.6, and 3.8 for the full, 1-halo, and 2-halo regimes respectively, while the 2-halo fit model had values of 26160, 55792, and 10.8.At this point, the only clear meaning of the 2-halo fit χ 2 dof values is that fitting this model to the 2-halo regime does not work without extra information not yet captured by our model.We remind the reader that at this point, we still have made no attempt to refine the radial alignment model, but rather we have thus far only illustrated why we have chosen to fit our parameters to the 1-halo regime alone.In Section 6.4, we turn attention to fitting distance-dependent radial alignment parameters and refining the goodness of fit.
We conclude this section by noting that our results imply a need for an additional modeling ingredient that captures the effect of environment beyond what is contained in the radial vector, similar to such ingredients that are now commonly used to capture assembly bias in HODs (e.g., Hearin et al. 2016;Yuan et al. 2018).The nature of this correction may require appeal to some measure of large-scale density smoothed on some scale, and/or tidal information.We acknowledge this limitation, however, as demonstrated here, satellite alignments based on radial vectors are a simple, flexible, and effective model when simulated subhalos are unavailable.The development of a model with additional environmental dependence is beyond our current scope and is left as a topic for future work.

Satellite Anisotropy
In §6.1-6.2, we restricted consideration to HOD-style models of satellite alignments in which the spatial positions of satellites were assigned based on the positions of their respective subhalos.We know that such satellite anisotropy affects correlations (Samuroff et al. 2020), and for IA models of satellites based on the radial vector, anisotropy in satellite locations has an influence on the orientation correlation functions.In this section we study , with that of the galaxies, ω g cs (colored lines).Two galaxy series are shown, each one having been generated by an HOD model where the galaxy alignment strengths were fit using MCMC.Bottom: The ratio of the galaxy correlation function over the halo correlation function (ω g /ω h ).The series where the parameters were fit to the 1-halo part of the halo correlation functions fits closely in the 1-halo regime (agreeing within 15%), as expected, while failing to accurately represent the 2-halo regime (differing by up to 40%).Meanwhile, the HOD with alignment strengths obtained by fitting to the 2-halo regime fail to model the 1-halo signal (differing by up to 90%), and don't follow the 2-halo correlation function well (only agreeing within 40%) despite having been fit to that area.
the effect of satellite isotropy and anisotropy on measurements of IA by comparing an isotropic phase space model with anisotropic satellite placement.
We model the anisotropic number density profile of satellites in halos as a triaxial NFW profile following Jing & Suto (2002) and Schneider et al. (2012).
where the relation between the Cartesian and elliptical coordinates is given by: where a, b, and c are the principle axis lengths of the host halo normalized such that a = 1.The anisotropy parameter, β, controls the magnitude of the anisotropy; when β = 0, satellites follow a spherically symmetric NFW profile; when β = 1, satellites follow a triaxial NFW profile with axis ratios equal to that of the underlying dark matter halo; when β > 1, the satellite distribution exhibits a greater degree of anisotropy than the underlying halo.
We compare three levels of anisotropy in the spatial distribution of satellites within their halos: • A subhalo position model, which places satellites at the locations of subhalos.
• A semi-isotropic model, which initially places galaxies using subhalo positions, but then randomly rotates each host halo's satellite population about the halo major axis.This model therefore preserves correlations involving the major axis but erases other anisotropies in the satellite distribution.
• An isotropic model, which places galaxies following a spherically symmetric NFW profile (as in §6.1-6.2.
We stress that in all three models sketched above, the only difference between the models is the level of anisotropy in the intra-halo positions of satellites.For the model of orientations explored in this section, here we assume that centrals are perfectly aligned with the halo shape, and for satellites we adopt a radial alignment model with constant strength of 0.85.As shown in Figure 8, both the subhalo position and semi-isotropic models allow the galaxy correlation functions to match those of the halos within 10% for the ω correlation functions and within 30% for the η correlation function.Meanwhile, models with isotropic satellites predict much weaker correlations, especially in the 1-halo regime, resulting in a 30% and a 60% difference from the subhalo ω and η correlation functions, respectively.The results illustrated in Figure 8 demonstrate that the level of anisotropy of satellite positions within the halo comprises an important degree of freedom for orientation correlation functions and that the major axis of the anisotropy contains most of the relevant information.

Capturing Subhalo IA in Gravity-only Simulations
So far in this section, we have been primarily focused on the position-orientation correlation function; in this section we consider additional orientation correlation functions to study the flexibility of our satellite IA models.Figure 5 shows how three different correlation functions are impacted by satellite IA models with either constant or distance-dependent alignment strength.Each prediction was computed from 50 realizations of the model, and in all cases galaxy positions are placed at the center of (sub)halos.We remind the reader that the satellite fraction is ∼ 25% in the predicted galaxy population, and so our plots show centralto-satellite correlation functions of each type: the central position to satellite position correlation function, ξ, the central position to satellite orientation correlation function, ω, and the central orientation to satellite orientation correlation function, η.We plot central-satellite cross-correlation functions to isolate the effect of satellite alignments on the IA signal.
The orientation correlation function η is better fit in the 1-halo regime, again indicating that additional environmental information beyond simple radial alignment is needed in order to capture physically plausible satellite effects on IA in the 2-halo regime.As mentioned in Section 6.1, by fitting to the 1-halo regime, we found the best constant radial alignment strength to be 0.826, and the best distance-dependent alignment strength to have a = 0.804 and γ = −0.04.
Previously, we derived these values using MCMC, but as an alternate method, we can directly fit the Dimroth-Watson distribution to the PDF of misalignment angles in the simulation.The results of this fitting process are shown in Figure 9, which displays the bestfitting Dimroth-Watson curve to the distribution of misalignment angles of simulated subhalos.To study the distance-dependence of alignment strength, we bin the subhalos (or satellite galaxies in the case of comparing to results obtained via MCMC) by r/r vir , and perform separate fits on each bin.In the right panel of Figure 9, we show a power law fit to the resulting relationship between µ and r/r vir .
At first glance, both the model with distancedependent alignment strength and that with the constant alignment strength appear to give equally valid results for the ω correlation function, as they match within 20% overall and within 6% and 8%, respectively, in the 1-halo regime.However, there are noticeable differences when the full χ 2 dof is considered.Looking at the ω (ED) correlation function, the model with distancedependent alignment strength has a χ 2 dof of 4.3 overall, 3.9 in the 1-halo regime, and 4.0 in the 2-halo regime while the model with constant alignment strength shows a χ 2 dof of 6.3, 8.2, and 4.3 for the overall, 1-halo, and 2-halo regimes respectively.Similarly, if we look at the η (EE) correlation functions, we see that the distancedependent alignment strength matches the halo correlation functions within 10% overall and within 13% in the 1-halo regime, while the constant alignment strength model only matches within 24% overall and 15% in the 1-halo regime.From this, we see that there does seem to be minor improvement from using a distance-dependent alignment model, especially when we limit our investigation to the 1-halo regime where our model works best.However, we should note the small value for γ found above (γ = −0.04),indicating this power law dependence is very close to a constant radial alignment.
In the 2-halo regime, the radial alignment model fails to capture the environmental dependence of subhalo orientations.In order to capture physically plausible satellite contributions to the IA signal on large scales, improving the physical realism beyond a simple radial alignment model will be necessary.We will explore such modeling improvements in future work.

CAPTURING IA WITH REALISTIC COMPLEXITY: TESTING THE MODEL WITH ILLUSTRIS
In the previous sections, we studied how models of satellite IA are able to capture the orientation correlation functions of subhalos in gravity-only N-body simulations.We now turn attention in this section to the IA signal exhibited by galaxies in the IllustrisTNG hydrodynamical simulation.In so doing, we will study IA models of both central and satellite galaxies acting together in concert.Here we compare the correlations of all of the halos themselves (black points) to those of all of the galaxies (black lines), where the isotropy of galaxy placement is the only difference between the three galaxy correlation functions shown.The galaxy correlations shown are the average of fifty iterations (assigning position and alignment fifty times).While all three models align galaxies with respect to the radial vector, the subhalo position model uses subhalo positions to place galaxies while the isotropic model places the galaxies isotropically, and the semi-isotropic model uses the subhalo positions then rotates them around the halo major axis to preserve the radial distribution.Here we see that using the subhalo positions, even if we rotate them, produces a signal remarkably close to that of the subhalos.Placing galaxies isotropically appears to lose some information required to match the halo correlation functions with the parameters given.
For central IA, we use the shape of the parent halo as the reference vector with respect to which central galaxy orientations are (mis)aligned.For satellite IA, we consider two different types of model: one based on the orientation of subhalos, the second HOD-type model based on radial alignments of constant strength.The free parameters of our model modulate the strength of central (satellite) galaxies with respect to the parent (sub)halo, and for all our results we use emcee to fit the parameters regulating the alignment strength.
Prior to fitting the galaxy alignment strengths, we first adjust the occupation model parameters such that the halo occupation of our model more closely matches that of IllustrisTNG300.That is, we ensure that halos of a given mass are populated, on average, with the same number of galaxies as halos of comparable mass in IllustrisTNG300.We used the Zheng07Cens and Zheng07Sats for the central and satellite galaxy occupation models available in halotools, which calculate the mean central and satellite occupations following equations (2) and (5) of Zheng et al. (2007).Table C1 shows the occupation parameters found by fitting the central and satellite occupation functions used by Zheng07Cens and Zheng07Sats to the mean occupations seen in Illustris.Specifically, the parameters were fit on the region where the central occupation transitions from zero central galaxies, to one central galaxy.For simplicity, in our fits σ log M (the width of the cutoff profile) and α (asymptotic slope at high halo mass) were held fixed at 0.26 and 1.0 respectively and not allowed to vary.As noted previously, the models built using halotools tend to have roughly three central galaxies to every one satellite, giv-ing a satellite fraction of roughly 25% with the occupation parameters used.The Illustris satellite fraction varies between 30% and 40% satellite galaxies, depending on the stellar mass threshold chosen (with higher stellar mass thresholds yielding a lower percentage of satellite galaxies).
The corner plots in Figure 10 display the posteriors from our fits to the orientation correlation functions in Illustris.As indicated in the figure labels, we show results for fits to three separate galaxy samples, each defined by a different stellar mass cut.Sample 1 corresponds to log (M * ) ≥ 10.5, sample 2 to log (M * ) ≥ 10, and sample 3 to log (M * ) ≥ 9.5.Two broad trends are apparent in Figure 10.First, the central galaxies have a noticeably larger alignment strength than the satellite galaxies.Second, as the stellar mass threshold increases, so does the alignment strength of both populations.Both of the above trends can be more clearly seen in Figure 11, which shows how the best-fitting alignment strength µ changes with the stellar mass cut.
As seen in Figure 12, our best-fitting HOD model using constant radial satellite alignment generates correlation functions that agree reasonably well with IllustrisTNG.It is by now well established that an HOD model more complex than the one we use here is required for precision fits to the spatial clustering of hydro-simulated galaxies (e.g., Beltz-Mohrmann et al. 2020;Hadzhiyska et al. 2020).For the ξ correlations functions, the HOD models for samples one, two, and three produce correlations within 14.0%, 9.5%, and 12.0% respectively to their corresponding counterpart from the Illustris data.This level of agreement is a sufficient baseline for purposes of Left: Differential distribution of (sub)halo radial misalignment angles for extant (sub)halos with a peak mass greater than 10 12 h −1 M⊙.The best fit Dimroth-Watson distribution (µ = 0.83) is shown as a red line.The shaded region shows the variation in the alignment strength for 80% of subhalos (i.e. the alignment strength assigned to 80% of the subhalos based on where they fall on the power law fit in the right panel).We see here that a single Dimroth-Watson distribution provides a good fit to all (sub)halo misalignment angles.Right: Dependence of the alignment strength for subhalos as a function of radial position, scaled by the virial radius of the host halo.The black points are alignment strengths for bins of r/rvir (with the alignment strength of each bin using the process shown in the left panel) with error bars given by bootstrapping with 50 bootstrap realizations.The red line is a clipped power law fit to these points.As opposed to the left plot where we fit a single Dimroth-Watson distribution to all (sub)halos, we see here that if we bin (sub)halos by distance and fit a Dimroth-Watson distribution to the misalignment angles in each bin, the behavior is approximately described by a power law.In the models used to generate these corner plots, we used a central alignment for central galaxies and subhalo alignment for satellite galaxies.We did the same process using constant radial alignment for satellite galaxies, yielding similar corner plots not shown here.As a general trend, we see the alignment strengths increase with higher mass cutoffs.Left: Best parameters found were µcen = 0.802 and µsat = 0.449.Middle: Best parameters found were µcen = 0.683 and µsat = 0.274.Right: Best parameters found were µcen = 0.528 and µsat = 0.112.dof values for the central-satellite ω (positionorientation) and η (orientation-orientation) correlation functions comparing the functions from our HOD models with those from Illustris.Also included are the alignment strengths that best fit our HOD model to IllustrisTNG.As before, these are the same mass cutoffs seen in Figure 12.
broadly assessing the accuracy of our alignment strength models.
The center and right panels of Figure 12 display the level of success of our best-fitting IA model to capture the position-orientation and orientation-orientation correlation functions, respectively.Table 2 shows the χ 2 dof values for each of these fits.For the satellite IA, we find that a model based on radial alignment performs comparably well to a model based on subhalo alignment, which is encouraging for the prospects of HOD-type approaches to capturing orientation correlation functions.

DISCUSSION & CONCLUSIONS
In this work, we have studied forward modeling techniques that enable predictions of large-scale structure that include intrinsic alignments in galaxy shape.Our methodology is fully simulation-based: when predicting galaxy correlation functions, we first create a synthetic universe of galaxies by populating simulated halos, and then we compute summary statistics with the same point estimators used in the corresponding observational measurements.Our prediction pipeline is computationally efficient, and it has the practical capability to run likelihood analyses in the deeply nonlinear regime with MCMC sampling algorithms.We have demonstrated the effectiveness of our model by reproducing IA signals exhibited by subhalos in gravity-only simulations, as well as by galaxies in hydrodynamical simulations such as Il-lustrisTNG.
Our simulation-based methodology is flexible, and it is able to leverage information from subhalos when available in a high-resolution simulation, or alternatively to populate catalogs of host halos in lower-resolution simulations.When subhalo information is available in the underlying simulation, we can generate more physically realistic predictions by using the orientation of subhalos as the reference vector for satellite galaxies.When using simulations that do not include subhalo orientations (as is often the case due to the computational demands of survey-scale simulations), a radial satellite alignment model (in which satellites preferentially align with respect to the radial vector between central and satellite) can approximate reasonably well the IA signal exhibited by hydro-simulated galaxies in the 1-halo regime.
In detail, we find that a distance-dependent alignment strength, following equation ( 19), may provide a slightly closer match to the 1-halo term signal seen in the correlation functions of either gravity-only subhalos or hydro- Alignment using the central alignment and subhalo alignment models for the central and satellite galaxies respectively.The three points shown for each type of galaxy correspond to samples one, two, and three.The samples are separated based on the minimum mass threshold for the galaxies included in the sample, shown on the x-axis.simulated galaxies.However, even a satellite model with a simple constant radial alignment strength produces a reasonable correlation function on small scales; looking specifically at ω g cs , we see that the distancedependent alignment strength with the radial alignment model gives a χ 2 dof of 4.3 overall, while the constant alignment strength gives 6.3.
As shown in Figure 5, no form of radial alignment model is able to faithfully reproduce the central-satellite correlation functions in the 2-halo regime; this implies the need for the development of a new model of galaxy alignments that more flexibly incorporates environmental information, analogous to ongoing efforts to capture assembly-biased halo occupation statistics (e.g., Contreras et al. 2021;Yuan et al. 2022;Lange et al. 2023).Reference vectors based in part on the tidal field are a promising way to incorporate such environmental effects (Harnois-Déraps et al. 2022).
We find that the anisotropic spatial distribution of satellites within their host halos plays an important role in the strength of orientation correlation functions.Figure 8 shows that models based on the anisotropic positions of subhalos substantially improve upon the realism of the predictions.Meanwhile, models based on a spherically symmetric NFW profile perform poorly in the 1-halo regime.These results imply the need for IA models with improved sophistication in the intrahalo phase space distributions of satellites, echoing the well-established need for such ingredients to capture ordinary two-point clustering (e.g., Orsi & Angulo 2018;Hadzhiyska et al. 2022).
One of our main goals was to test how well the IA signal in IllustrisTNG can be mimicked with a simple HOD- type model of galaxy orientations.In our tests, we used MCMC to fit our models of central and satellite alignment strength to the orientation correlation functions of IllustrisTNG galaxies, repeating this analysis separately for samples defined by three different stellar mass cuts.Figure 12 summarizes the results of this analysis, and Table 2 records the alignment strengths of our bestfitting models.Subhalo-based models are particularly effective in their ability to generate a mock galaxy catalog with correlation functions similar to IllustrisTNG, attaining χ 2 dof close to unity for every sample and correlation function considered.Satellite models based on constant radial alignment perform comparably well to the subhalo models, despite their relative simplicity.In future work, we will generalize the HOD-type methodology used here to incorporate joint correlations between IA strength and galaxy brightness, color, and redshift as seen in simulations (Zjupa et al. 2022) and observations (Samuroff et al. 2023).We will also work towards adding more sophisticated alignment models that take into account more detailed sets of halo properties and other information, including redshift, in order to produce more robust results (Xu et al. 2023).
In a full forward model of galaxy correlation functions ξ, ω, and η, additional ingredients for the probabilistic distribution of galaxy shapes will be needed in addition to the alignment distribution model studied here.Once the model is able to accommodate full shape information, we expect that the results here would show much better agreement corresponding to the expanded functionality of our alignment models.Our mock-population framework lays the groundwork for carrying out Bayesian inference with simulation-based forward models of galaxy shape correlation functions.Our simulation-based approach simplifies the effort required to model systematic effects that can be challenging to capture analytically, such as galaxy assembly bias, halo exclusion in the trans-ition between 1-and 2-halo regimes (van den Bosch et al. 2013;García & Rozo 2019), and beyond-linear halo bias (Mahony et al. 2022;Mead & Verde 2021).Forward models based on mock-population methods also extend naturally to predict higher-order summary statistics, which can provide additional information not contained in the two-point functions studied here (Harnois-Déraps et al. 2021).With the additional effort outlined above, our models could be used to derive posteriors on the true intrinsic alignment strength of galaxies from observations on highly nonlinear scales, and to supply priors on the IA nuisance parameters used in cosmological analyses of larger scales.
where the half-lengths of the principle axis of the ellipsoid are given by a = √ λ a , b = √ λ b , c = √ λ c , (a ≥ b ≥ c).The orientation of this ellipsoid is then specified by the eigenvectors, êa , êb , and êc .
The radial weighting applied in equation (A1) tends to make (sub)halos/galaxies appear more spherical.To alleviate this, we iteratively calculate the reduced inertia tensor for (sub)halos/galaxies.After the initial calculation of the Ĩ, the particle distribution is rotated so that êa , êb , and êc are aligned with the x-, y-, x-axis.The radial distance, r 2 n in equation (A1), is replaced with the elliptical radial distance: This process is repeated until the eigenvalues change by less than 1% between iterations.We implicitly model (sub)halos/galaxies as homologous triaxial ellipsoids.Note the caveats here.halos contain substructure.galaxies may contain multiple components.galaxies are observed to have twisting isophotes.By modeling galaxies as homologous ellipsoids, we can specify a single shape and orientation, vastly simplifying the modeling.

B. FITTING FUNCTIONS
In Figures 3 and 4, we fit a smooth function to our measurements of simulated (subhalo) η (EE) measurements.We primarily do this to reduce the noise in the measurement at any given radius in order to make the figures more clear.We present our fitting function here for completeness.As shown in §7, the alignment correlation functions simultaneously depend on the misalignment distributions as well as the underlying halo occupation statistics.Thus in order to test the efficacy of our misalignment modeling using IllustrisTNG, it's important to first establish an accurate baseline HOD model.In Table C1, we see the occupation model parameters used.These values were found by fitting the occupation parameters in the Zheng07Cens and Zheng07Sats occupation models to fit the occupation distribution of our HOD models with that of Illustris.Specifically, we fit the parameters in the region where the central occupation transitions from zero to one.Note that the values σ log M and α have been fixed for simplicity, only allowing the other three parameters to vary.

D. COVARIANCE
Due to the probabilistic nature of populating an HOD model and aligning with a Dimroth-Watson distribution (discussed later in Section 4), there is stochasticity in both the data vector and the model.We have two different ways in which we select our covariance, depending on our data vector.
In Section 6, we fit our satellite galaxy alignment strengths to find an alignment strength for the radial alignment such that the central position, satellite orientation correlation functions for the galaxies match as closely as possible to those of the halos (keeping central galaxies perfectly aligned with their host halos for simplicity).However, because there is an element of stochasticity both in the population of the dark matter halos with galaxies, and in the realignment of the galaxies, both our data vector and our model have covariance.To account for this, we look at the contribution to the covariance of repopulating halos versus that of assigning galaxy alignments.As seen in Figure D1, we find that most of the covariance comes from aligning galaxies.As such, we use a fixed seed for the galaxy population step in our MCMC, meaning the same dark matter halos are populated on each run, keeping the galaxy alignment step the only stochastic process.For our covariance, we use the covariance of ω g cs (the central galaxy position, satellite galaxy orientation correlation function) from 100 instances of realigning galaxies at that seed, for a given alignment strength (see Section 6 for more detail).We use this covariance in the MCMC, after which a new alignment strength is determined, and the process repeats (generating a new covariance from 100 new iterations at the new alignment strength) until the alignment strength resulting from the MCMC does not change appreciably.
In Section 7, rather than fitting to the halo orientations, we found both central and satellite galaxy alignment strengths by fitting ω g gg (the gg subscript now referring to all galaxies) correlation functions to those of IllustrisTNG300.In this case, since we had a fixed data vector, we used jackknife resampling of the Illustris data, splitting the simulation box into 5x5x5 jackknife regions to determine our covariance.It is important to note that we multiplied this jackknife covariance by the Hartlap factor, shown in equation (D1) (taken from equation 17 in Hartlap et al. (2007) where Ĉ−1 is the adjusted inverse covariance, Ĉ−1 * is the inverse covariance directly from jackknife resampling, n is the number of independent observations (i.e. the number of jackknife regions), and p is the number of elements in the data vector (i.e. the number of bins).This helps approximately debias our estimation of the original inverse jackknife covariance.The diagonals of the covariance matrices generated from measuring the ed correlations of 100 instances of mock galaxy catalogs.The curve for total covariance comes from catalogs where we performed a full repopulation (i.e.no random seed was used and each time, halos were selected and populated randomly before assigning galaxy orientation).
The orientation curve was generated in a similar fashion, but instead of allowing a full repopulation, we used a fixed seed for the halo population (i.e. the same halos were chosen each time with the same galaxy populations) meaning the only source of variance was the galaxy alignments.The final curve shows the difference between these two, Total -Orientation, representing the covariance that comes from just the repopulation of the halos, ignoring the galaxy alignments.Here we see that the majority of the covariance comes from realigning the galaxies.
data we used is publicly available at https://www.tngproject.org/.

Figure 1 .
Figure1.Cartoon of alignment model when subhalos with accurate shapes are available (top) and without subhalos (bottom).In both cases, the misalignment angle is determined with respect to some reference axis, which differs for satellite galaxies depending on the model used, either using subhalo alignment (top), or radial alignment (bottom).

Figure 3 .
Figure3.The effect of degrading the central galaxy-host halo alignment strength on the ω (left) and η (right) orientation correlation functions.The black points with error bars (which are often smaller than the points) are measurements made directly on (sub)halos.The colored lines are for model galaxies with varying levels of alignment strength between central galaxies and their host halo, from random alignments (light blue) to perfect alignments (pink).In each model, satellite galaxies take on the same orientation as their subhalo.On the right, a fitting function is used to smooth the results, letting us focus on the overall trend from lowering the central alignments which tends to be noisy for ω correlations.We also see a downward bump around the 1-halo to 2-halo scale transition in the right panel.This suggests to us an effect related to the satellite-satellite correlations, with the sudden upturn going into the1-halo regime being related to satellite pairs within a halo becoming more significant.This effect also only becomes apparent once the central alignment strength has become weak enough that the satellite alignments become more important.

Figure 4 .
Figure4.The effect of degrading the satellite galaxy-subhalo alignment strength on the ω (left) and η (right) orientation correlation functions.The black points with error bars are measurements made directly on (sub)halos.The colored lines are for model galaxies with varying levels of alignment strength between satellite galaxies and their subhalo, from random alignments (light blue) to perfect alignments (pink).In each model, central galaxies take on the same orientation as their host halo.On the right, a fitting function is used to smooth the results.

Figure 5 .
Figure 5. Left: ξcs correlation function (central galaxy position correlated with satellite galaxy position).Galaxy correlation functions are shown in the solid red and dashed blue lines, with the corresponding halo correlation function shown in black dots.Here, all three share the same ξ (position-position) correlation function as the only difference in the model is the alignment.Middle: ωcs correlation function (central galaxy position and satellite galaxy orientation).Here, we see a slight difference between ω g cs (the galaxy correlation functions) and ω h cs (the corresponding halo correlation function) with the distancedependent alignment (solid red line) performing slightly better at small scales than the constant alignment strength (shown as a dashed blue line).Right: ηcs correlation function (central galaxy orientation and satellite galaxy orientation).In all panels (and in figures below), the shaded regions show the 1-sigma confidence interval of each model.

Figure 6 .
Figure6.The position-orientation correlations for central-satellite (left; ωcs), satellite-satellite (middle; ωss), and all galaxysatellite (right; ωgs) pairs.All correlations shown are using a constant radial alignment model.Left: ωcs function; an identical configuration to the constant radial alignment series in the middle panel of Figure5.Here we see good agreement to the 1-halo regime with slight underestimation in the 2-halo regime.Middle: ωss correlation function.Here we see a generally good agreement except at small scales.Right: ωgs correlation function.We see that, although ωss performs poorly at small scales, the ωcs portion dominates making the overall correlation a good fit.

Figure 7 .
Figure 7. Top: Comparison of the central-satellite positionorientation correlation function for halos, ω hcs (black points), with that of the galaxies, ω g cs (colored lines).Two galaxy series are shown, each one having been generated by an HOD model where the galaxy alignment strengths were fit using MCMC.Bottom: The ratio of the galaxy correlation function over the halo correlation function (ω g /ω h ).The series where the parameters were fit to the 1-halo part of the halo correlation functions fits closely in the 1-halo regime (agreeing within 15%), as expected, while failing to accurately represent the 2-halo regime (differing by up to 40%).Meanwhile, the HOD with alignment strengths obtained by fitting to the 2-halo regime fail to model the 1-halo signal (differing by up to 90%), and don't follow the 2-halo correlation function well (only agreeing within 40%) despite having been fit to that area.

Figure 8 .
Figure8.The effect of subhalo anisotropy on η and ω.Here we compare the correlations of all of the halos themselves (black points) to those of all of the galaxies (black lines), where the isotropy of galaxy placement is the only difference between the three galaxy correlation functions shown.The galaxy correlations shown are the average of fifty iterations (assigning position and alignment fifty times).While all three models align galaxies with respect to the radial vector, the subhalo position model uses subhalo positions to place galaxies while the isotropic model places the galaxies isotropically, and the semi-isotropic model uses the subhalo positions then rotates them around the halo major axis to preserve the radial distribution.Here we see that using the subhalo positions, even if we rotate them, produces a signal remarkably close to that of the subhalos.Placing galaxies isotropically appears to lose some information required to match the halo correlation functions with the parameters given.
Figure9.Left: Differential distribution of (sub)halo radial misalignment angles for extant (sub)halos with a peak mass greater than 10 12 h −1 M⊙.The best fit Dimroth-Watson distribution (µ = 0.83) is shown as a red line.The shaded region shows the variation in the alignment strength for 80% of subhalos (i.e. the alignment strength assigned to 80% of the subhalos based on where they fall on the power law fit in the right panel).We see here that a single Dimroth-Watson distribution provides a good fit to all (sub)halo misalignment angles.Right: Dependence of the alignment strength for subhalos as a function of radial position, scaled by the virial radius of the host halo.The black points are alignment strengths for bins of r/rvir (with the alignment strength of each bin using the process shown in the left panel) with error bars given by bootstrapping with 50 bootstrap realizations.The red line is a clipped power law fit to these points.As opposed to the left plot where we fit a single Dimroth-Watson distribution to all (sub)halos, we see here that if we bin (sub)halos by distance and fit a Dimroth-Watson distribution to the misalignment angles in each bin, the behavior is approximately described by a power law.

Figure 10 .
Figure10.The parameter values found from MCMC for central alignment strength (µcen) and satellite alignment strength (µsat) to let HOD models best fit Illustris correlation functions using a subhalo alignment model for the satellite galaxies.Each sample corresponds to a different mass cutoff for the halos included in the HOD model.These are the same mass cutoffs seen in Figure12.In the models used to generate these corner plots, we used a central alignment for central galaxies and subhalo alignment for satellite galaxies.We did the same process using constant radial alignment for satellite galaxies, yielding similar corner plots not shown here.As a general trend, we see the alignment strengths increase with higher mass cutoffs.Left: Best parameters found were µcen = 0.802 and µsat = 0.449.Middle: Best parameters found were µcen = 0.683 and µsat = 0.274.Right: Best parameters found were µcen = 0.528 and µsat = 0.112.

Figure 12 .
Figure12.Left: Two-point correlation functions offset by 1-dex for clarity Middle: position-orientation correlation function Right: orientation-orientation correlation function.In each panel, the points with error bars are measurements made on the IllustrisTNG300-1 simulation with error bars estimated using jackknife re-sampling of the box.The lines with shaded regions are halo model predictions made by populating a DMO simulation with mock galaxies where the shaded region shows the variation from random realizations of the model.We used a constant radial satellite alignment model in each case to showcase the flexibility of the model in cases where subhalo orientations are not available.The three colors are for three stellar mass threshold samples.
η(r) = A 1 exp (−(r/R 1 ) γ ) + A 2 exp (−(R 2 /r)) 1.0 + (r/R 2 ) α (B1) ω(r) =A 1 exp (−(r/R 1 ) γ ) + A 2 exp (−(R 2 /r)) 1.0 + (r/R 2 ) α Figure D1.The diagonals of the covariance matrices generated from measuring the ed correlations of 100 instances of mock galaxy catalogs.The curve for total covariance comes from catalogs where we performed a full repopulation (i.e.no random seed was used and each time, halos were selected and populated randomly before assigning galaxy orientation).The orientation curve was generated in a similar fashion, but instead of allowing a full repopulation, we used a fixed seed for the halo population (i.e. the same halos were chosen each time with the same galaxy populations) meaning the only source of variance was the galaxy alignments.The final curve shows the difference between these two, Total -Orientation, representing the covariance that comes from just the repopulation of the halos, ignoring the galaxy alignments.Here we see that the majority of the covariance comes from realigning the galaxies.

Table 1 .
Cosmology and other parameters for the three simulations described in Section 3. Additional information about each simulation can be found in said section.