# Strong orientation dependence of surface mass density profiles of dark haloes at large scales

###### Abstract

We study the dependence of surface mass density profiles, which can be directly measured by weak gravitational lensing, on the orientation of haloes with respect to the line-of-sight direction, using a suite of -body simulations. We find that, when major axes of haloes are aligned with the line-of-sight direction, surface mass density profiles have higher amplitudes than those averaged over all halo orientations, over all scales from to we studied. While the orientation dependence at small scales is ascribed to the halo triaxiality, our results indicate even stronger orientation dependence in the so-called two-halo regime, up to . The orientation dependence for the two-halo term is well approximated by a constant shift of the amplitude and therefore a shift in the halo bias parameter value. The halo bias from the two-halo term can be overestimated or underestimated by up to depending on the viewing angle, which translates into the bias in estimated halo masses by up to a factor of two from halo bias measurements. The orientation dependence at large scales originates from the anisotropic halo-matter correlation function, which has an elliptical shape with the axis ratio of up to . We discuss potential impacts of halo orientation bias on other observables such as optically selected cluster samples and a clustering analysis of large-scale structure tracers such as quasars.

###### keywords:

cosmology: theory – methods: numerical – large-scale structure of Universe^{†}

^{†}pubyear: 2017

^{†}

^{†}pagerange: Strong orientation dependence of surface mass density profiles of dark haloes at large scales–Strong orientation dependence of surface mass density profiles of dark haloes at large scales

## 1 Introduction

Weak gravitational lensing provides an important means of measuring matter distributions around galaxies and clusters (e.g., Bartelmann & Schneider, 2001, for a review). In particular, cross-correlating galaxy shapes with positions of foreground galaxies or clusters, which is referred to as stacked weak lensing or galaxy-galaxy lensing, allows one to measure the average matter distribution around the foreground tracers, which can be in turn used to constrain the relation between luminous and matter distributions in the large-scale structure (e.g., Brainerd et al., 1996; Sheldon et al., 2004; Mandelbaum et al., 2006; Okabe et al., 2010; Leauthaud et al., 2012; Okabe et al., 2013; Umetsu et al., 2014; Miyatake et al., 2015; Murata et al., 2017; Prat et al., 2017). Combining the stacked lensing with the auto-correlation function of the same foreground tracers can be used to constrain cosmological parameters to calibrate the bias uncertainty simultaneously (e.g., Baldauf et al., 2010; Oguri & Takada, 2011; More et al., 2015).

The stacked lensing profile of galaxies or clusters consists of two distinct regimes, the “one-halo” and “two-halo” terms in the halo model picture (e.g., Cooray & Sheth, 2002, for a review). Here the one-halo term arises from matter distribution in the same halo, and the lensing signal can be used to constrain the interior mass, i.e. the total mass of gravitationally-bound matter in the haloes. On the other hand, the two-halo term arises from matter in different haloes, more generally matter distributed in the large-scale structure. The two-halo term contains cleaner information on the matter distribution that is relatively easier to theoretically model and is less affected by baryonic physics, which can be used to infer the halo bias including the assembly bias (e.g., Covone et al., 2014; Lin et al., 2016; Miyatake et al., 2016; More et al., 2016; Dvornik et al., 2017; Busch & White, 2017) and/or the underlying matter power spectrum (e.g., Lombriser et al., 2012; Mandelbaum et al., 2013; Sereno et al., 2015; More et al., 2015; van Uitert et al., 2017a; Dark Energy Survey Collaboration, 2017). If the foreground tracers are selected in an unbiased manner such that the matter distribution around the tracers is statistically isotropic, the stacked weak lensing profile is computed from a projection of the spherically averaged three-dimensional halo-matter cross-correlation function at the redshift of foreground tracers (Mandelbaum et al., 2005; Hikage et al., 2013; More et al., 2015).

However, if the sample of foreground tracers has selection effects depending on the line-of-sight direction, referred to as “projection effects” in this paper, it violates the spherical symmetry and could cause a bias in parameters estimated from the weak lensing measurements. For instance, it is often advocated that strong lensing clusters (clusters found via observations of strong lensing phenomena) are affected by selection effects depending on an alignment of the major axis of mass distribution to the line-of-sight direction (Hennawi et al., 2007; Oguri & Blandford, 2009; Meneghetti et al., 2010), because projected surface mass density profiles are sensitive to the orientation of dark haloes (Clowe et al., 2004; Oguri et al., 2005; Gavazzi, 2005; Corless & King, 2007; Limousin et al., 2013). This originates from the fact that dark matter haloes are highly non-spherical, but rather triaxial with a typical major-to-minor axis ratio of for cluster-scale haloes, reflecting the collision-less nature of dark matter (Jing & Suto, 2002; Allgood et al., 2006; Schneider et al., 2012; Vera-Ciro et al., 2014; Vega-Ferrero et al., 2017), which has also been confirmed directly by weak lensing observations (Evans & Bridle, 2009; Oguri et al., 2010; Oguri et al., 2012; Clampitt & Jain, 2016; van Uitert et al., 2017b; Shin et al., 2017). Hence, the line-of-sight projection of the matter distribution for such a sample under the selection effects complicates an interpretation of the lensing observables compared to the case of spherical symmetry. These effects need to be properly taken into account in order to derive an unbiased estimation of the model parameters as well as to use the weak lensing observables to constrain cosmological parameters.

In contrast, it is not well understood how projection effects change surface mass density profiles of dark haloes, especially at large scales in the two-halo term regime. It is known that three-dimensional mass distributions around dark haloes correlate with their shapes, the so-called intrinsic alignment (e.g., Troxel & Ishak, 2015; Joachimi et al., 2015, for a review). This intrinsic alignment implies that there is a larger amount of matter along the major axis of a halo and a smaller amount of matter along the minor axis of a halo than the spherical average, and the correlation extends out to very large scales. However, it is not obvious how the intrinsic alignment affects projected mass distributions, i.e., surface mass density profiles, as a function of the viewing angle of haloes, because the larger amount of matter along the major axis would be compensated by the smaller amount of matter along the minor axis when the matter distribution is projected along the line-of-sight.

In this paper, we study the dependence of surface mass density profiles around dark haloes on viewing angle, focusing on cluster-scale dark haloes. Specifically, we use high resolution -body simulations to quantify the orientation effect on surface mass density profiles from small to large scales up to . We discuss its impact on measurements of halo biases and halo masses from large scale surface mass density profile measurements. Our results will have an important implication for the analysis of ongoing surveys such as Dark Energy Survey (Dark Energy Survey Collaboration, 2016) and Hyper Suprime-Cam survey (Aihara et al., 2017)

This paper is organized as follows. We first review the halo model calculation in Section 2. In Section 3, the details of -body simulations and numerical methods are explained. We present the measurements of surface mass density profiles and the fitting results in Section 4, and discuss the origin of the orientation effect in Section 5. In Section 6, we discuss implications of our results from -body simulations on the analysis of actual observations. Finally, we give our conclusion in Section 7. Throughout this paper, we adopt the flat cold dark matter model with , , , , and at a pivot wavenumber based on Planck measurements of the anisotropy of temperature and polarization (TT, TE, EE+lowP) of cosmic microwave background (Planck Collaboration, 2016).

## 2 Halo model

In this Section, we will briefly review the halo model calculation of surface mass density profiles around haloes, which will be used to estimate possible biases in estimating halo masses from the surface mass density profiles simulated from -body simulations.

### 2.1 Halo density profile

The density profile characterizes the structure of dark haloes. Using -body simulations, Navarro et al. (1996, 1997) proposed a universal form of the spherically averaged radial density profile, the so-called Navarro–Frenk–White (NFW) profile

(1) |

where is the scale density and is the scale radius. The virial mass is defined as

(2) |

where is the virial overdensity, is the critical density, and is the virial radius. For , we use the following formula based on spherical collapse model (Bryan & Norman, 1998)

(3) |

where is matter density at redshift

(4) |

The scale density can be determined from the virial mass as

(5) |

where

(6) |

The parameter is the concentration parameter, which is the ratio of the virial radius to the scale radius, i.e., . Previous studies (e.g., Duffy et al., 2008; Diemer & Kravtsov, 2015) based on -body simulations have shown that the mean relation of the concentration parameter is tightly correlated with halo mass and redshift, and there is a sizable scatter around the mean relation, typically by an amount of for cluster-scale haloes.

The NFW profile has an asymptotic form as for , which indicates that the total enclosed halo mass diverges at . Thus, in the halo model calculations the prescription to truncate the NFW density profile at the virial radius is commonly used (Takada & Jain, 2003a, b). However, this truncation produces an unphysical discontinuity in the weak lensing profiles that are obtained from the line-of-sight projection of an isolated NFW profile. In order to avoid this problem, Baltz et al. (2009) proposed another density profile (hereafter the BMO profile) which has the steeper asymptotic form at outer radii

(7) |

where is the truncation radius. Based on the comparison of halo model calculations with ray-tracing simulations, Oguri & Hamana (2011) proposed to adopt the following value for the truncation radius for ;

(8) |

which we also adopt in this paper. The BMO profile with has a steep outskirts (), and as a result the total enclosed mass converges quickly. Throughout this paper, we use the BMO profile as a fiducial density profile of dark haloes.

### 2.2 Surface mass density profile

Under the assumption that all matter in the Universe is associated with haloes, we can analytically compute the surface mass density profile from small to large scales using the halo model (Guzik & Seljak, 2002; Mandelbaum et al., 2005). In this framework, the surface mass density profile as a function of projected radius , , can be decomposed into the one-halo term and the two-halo term as

(9) | |||||

(10) | |||||

(11) |

where is the mean matter density in the present Universe, is the linear halo bias, is the zeroth order Bessel function, and is the linear matter power spectrum. We compute the linear matter power spectrum using the CAMB code (Lewis et al., 2000). Also see Takada & Jain (2003a) and Oguri & Takada (2011) for the halo model formulation of lensing profiles based on an use of the two- and three-dimensional Fourier transforms of the halo density profile.

In weak lensing observations, we usually measure tangential shear profiles with respect to the center of each foreground halo, which measure “excess” surface mass density profiles rather than . The excess surface mass density profile is related to the surface mass density profile as

(12) |

where is the mean surface density within a circular aperture of radius . Since and carry the same information, in this paper we focus only on . We note that our results on can easily be converted to those on via equation (12).

## 3 Simulations

### 3.1 -body simulations

In this Section, we describe the details of simulations used in the analysis. In order to obtain the distributions of matter and haloes in the Universe, we run -body simulations. For this purpose, we use Tree-PM code Gadget-2 (Springel, 2005). The number of particles is and the length of the simulation box on a side is . The corresponding particle mass is . The initial conditions at the redshift are generated with the parallel code developed in Nishimichi et al. (2009); Nishimichi et al. (2010) and Valageas & Nishimichi (2011), which employs the second-order Lagrangian perturbation theory. We generated 24 initial conditions with different random seeds and ran the -body simulations to obtain the matter distribution at the present Universe, i.e., .

### 3.2 Halo identification and halo shape measurement

Next, we run the halo finding algorithm, the Rockstar (Behroozi et al., 2013), in each -body simulation output. The minimum halo mass which we use in the analysis is . Hence, haloes used in the analysis contain at least particles, which is sufficient for reliable shape measurements of individual haloes (Jing & Suto, 2002). Hereafter, the halo samples are divided into three bins, , , and , according to their virial masses, respectively.

We then compute inertia tensor for each halo. While there are various definitions of the inertia tensor (Bett, 2012, and references therein), in this paper we employ the reduced inertia tensor determined by an iterative scheme to characterize the shape of dark haloes. The reduced tensor for a halo with particles is defined as

(13) |

where is the triaxial radius of the -th member particle in the halo at the -th step. At the first step (), the triaxial radius is the same as the one in Cartesian coordinates. At the -th () step, the triaxial radius is measured in the principal coordinates at the last (-th) step, i.e.,

(14) |

where , and are coordinates along the minor, intermediate, and major axis directions, respectively, and () is the minor-to-major (intermediate-to-major) axis ratio. Since the Cartesian coordinate system is adopted in the first step, the axis ratios are unity, i.e.,

(15) |

We iteratively compute the tensor until the fractional difference of axial ratios between two steps becomes less than . The sum in equation (13) runs over all the particles whose triaxial radii are less than . The eigenvalue (eigenvector) of the inertia tensor corresponds to the length (direction) of the principal axis.

Hereafter, we focus on the direction of the major axis with respect to
the line-of-sight direction, which is chosen to be the
-axis in the simulation box. For each mass bin, we divide the halo
sample based on the directional cosine between the major axis and the
line-of-sight direction, which we denote as .
We divide the halo sample into five bins that are equally spaced with
respect to ^{1}^{1}1The area element of the unit sphere in
polar coordinates is .
Therefore, the number of haloes which belong to each bin is almost
the same when the direction of the major axis is random..

## 4 Results

### 4.1 Halo-matter cross-spectrum

Since the Universe has no preferred direction, spherically averaged mass density profiles around dark haloes should not depend on their orientations. As a sanity check, we calculate three-dimensional halo-matter cross-spectra in the simulations for the five orientation bins to make sure that these cross power spectra, which correspond to spherically averaged mass profiles in real space, do not depend on the halo orientation. Our results in Figure 1 indeed confirm that the cross-spectra are consistent with each other among different orientation bins and the variation is well below statistical scatters over 24 realizations.

### 4.2 Surface mass density profile

In Figure 2, we show surface mass density profiles for different orientation bins, as well as those averaged over all halo orientations. For all mass bins, clear dependence on the orientation can be seen for a wide range of radii, from up to . When we observe haloes whose major axes are aligned with the line-of-sight direction, the amplitudes of the surface mass density profiles are always larger than the average surface mass density profile. The cross-correlation between the intra-halo structure at a few and the large-scale matter distribution might be somewhat surprising. The existence of strong cross-correlation implies that the matter distributions at the totally different scales co-evolve in nonlinear structure formation. We will later discuss the nature of the cross-correlation in more detail.

To study the orientation bias in the surface mass density profiles more quantitatively, we fit the surface mass density profiles using the halo model described in Section 2. We fit surface mass density profiles in simulations with equation (9), leaving , , and as free parameters. For an estimation of the model parameters, we use the standard chi-square fitting:

(16) |

where the fitting range is with 68 equally log-spaced bins, is the surface mass density measured in the simulations, is the halo model prediction given parameters and is the variance of the surface mass density over 24 simulations. In equation (16), we ignore covariance, i.e. correlation between different bins, because the number of realization is too small to accurately estimate covariance and we focus on the parameter bias rather than errors on best-fit parameter values.

The best-fit halo model results shown in Figure 2 indicate that the radius where the fractional differences between haloes with different orientations become minimum () roughly corresponds to the transition between the one-halo and two-halo terms. This implies that the orientation dependences at small and large radii have different origin.

In fact it has been known that the halo triaxiality is the main cause of the orientation dependence of the surface mass density profile in the one-halo regime (Oguri et al., 2005; Gavazzi, 2005; Corless & King, 2007; Limousin et al., 2013). We explicitly check this by computing surface mass density profiles for different halo orientations expected by the triaxial halo model of Jing & Suto (2002). The triaxial density profile is parametrized as

(17) |

where is the characteristic overdensity that is related to the virial overdensity (see Jing & Suto, 2002), is the scale radius, and is the triaxial radius. The triaxial radius is different from the ordinary radius in terms of reference coordinate system and is defined as

(18) |

where , and are Cartesian coordinates in the principal coordinate system of the triaxial halo, and , , and is the length of minor, intermediate, and major axis (i.e., ), respectively. We derive surface mass density profiles of the triaxial halo model following Oguri et al. (2003) and Oguri & Blandford (2009). In particular, for each mass bin we compute surface mass density profiles as a function of the halo orientation using the Monte-Carlo method developed in Oguri & Blandford (2009). In this method, in each mass bin, we randomly generate haloes following the halo mass function of Bhattacharya et al. (2011), and assign their axis ratios and concentration parameters following the probability distribution functions derived in Jing & Suto (2002). We also assign orientations of these haloes with respect to the line-of-sight direction assuming they are randomly oriented. We then connect project triaxial density profiles with corresponding two-dimensional mass () and concentration parameters (), and use the for the spherical NFW halo (equation 10) to predict surface mass density profiles of individual triaxial haloes generated by the Monte-Carlo method. This result is used to compute surface mass density profiles for each bin for the triaxial halo model.

In Figure 3, we compare the results from -body
simulations with the triaxial halo model predictions. We find that the
triaxial halo model reproduces the orientation dependence of surface
mass density profiles at small scales () reasonably well.
While the overall amplitudes of the ratios are slightly
different between the simulations and the triaxial
halo model presumably due to the inaccuracy of the triaxial halo model
especially at high mass end^{2}^{2}2Bonamigo et al. (2015) presented
the probability distribution function of axis-ratios with high resolution
-body simulation. Such extension from Jing &
Suto (2002)
may resolve the discrepancy.,
the model well reproduces the general
trend of the ratios in the one-halo regime, including the radius
dependence of the ratios and the relative behavior of the ratios as a
function of . This comparison supports the idea that the
orientation dependence of surface mass density profiles at small
scales () can mostly be explained by the halo
triaxiality, and the orientation dependence of surface mass density
profiles at large scales () has different origin,
which we will discuss in detail in Section 5.

### 4.3 Bias in estimating halo masses

The halo orientation changes the amplitude of the surface mass density profile at all scales. This shift of the amplitude inevitably causes a bias in cluster properties inferred from an observation of the surface mass density profile, if the cluster properties (parameters) such as halo mass are estimated from the model fitting assuming the spherical symmetry. We ddress such a bias by making a hypothetical model fitting of the surface mass density profiles in simulations with the spherically symmetric halo model presented in Section 2. In doing so, we adopt, as model parameters, the halo mass and the concentration parameter for the one-halo term, and the halo bias for the two-halo term for simplicity. Figure 4 shows the estimated parameters for individual orientation bins relative to the best-fit parameter values for all haloes denoted as , , and . For the halo mass, in addition to the halo mass estimated directly from the one-halo term, we also show the halo mass converted from the estimated halo bias using the fitting formula in Tinker et al. (2010), which we call the two-halo mass.

We find that all the parameters show clear dependence on the halo orientation. For the parameters of the one-halo term, we find that the halo mass and concentration parameter can be overestimated or underestimated up to depending on the viewing angle, which is consistent with the previous analysis based on the triaxial halo model (Oguri et al., 2005). On the other hand, we find that the best-fit halo bias from the two-halo term strongly depends on the halo orientation, which has not been recognized before. We note that the effect of the halo orientation on surface mass density profiles in the two-halo regime can well be approximated by a constant shift of the amplitude over a wide range of radii (see Figure 2), which indicates that the effect can be modeled very well by a change of the value of the halo bias which is assumed to be scale-independent in our halo model. We find that the halo bias can be overestimated or underestimated by up to depending on the viewing angle. This orientation dependence of the halo bias translates into the orientation dependence of the two-halo mass by up to nearly a factor of two, which is much larger than that of the halo mass derived directly from the one-halo term. This large effect of the halo orientation on the two-halo term highlights the importance of the orientation effect on the analysis of surface mass density profiles.

The different dependences of the halo mass and two-halo mass on the halo orientation suggests that it may be possible to break the degeneracy between the halo mass and orientation based on detailed measurements of surface mass density profiles in both the one-halo and two-halo regimes. We leave the exploration of this possibility for future work.

## 5 Origin of the strong orientation dependence at large scales

### 5.1 Effect of projection thickness

The intrinsic alignment of dark haloes implies that there is a larger amount of matter along the major axis of the dark halo, and a smaller amount of matter along the minor axis for compensation. These two effects counteract with each other when we consider the projection along the line-of-sight. This argument also suggests that the orientation dependence of surface mass density profiles should evolve considerably as a function of the projection thickness. For instance, for the halo whose major axis is aligned with the line-of-sight direction (), while they have larger two-halo term amplitudes than the average when the projection thickness is sufficiently larger, the two-halo term amplitudes should become smaller than the average for the smaller projection thickness because only the matter near the direction of the minor axis, which is underdense than the average, is accounted in the projected profile.

We investigate the effect of the projection thickness as follows. In deriving our results, we project matter over the whole simulation box, i.e., the projection thickness is . Here we measure the surface mass density profiles with six different projection thickness, , , , , , and . First, we divide the whole simulation box into 100 slices with thickness, and construct a halo number density field for each slice. Then we compute projected matter density field for each slice projected over the given thickness listed above. We can obtain the surface mass density profile by calculating the cross-correlation between projected halo and matter field with the fast Fourier transform (FFT). We have 100 pairs of projected halo and matter density field, leading to measurements of 100 surface mass density profiles. Finally, we average these 100 measurements.

In Figure 5, we show surface mass density profiles with different projection thickness. As expected, for the case of a thin projection thickness, e.g., , the orientation dependence of surface mass density profiles at large scales show the opposite trend as compared with the case of the full projection. The matter distribution at large separations is lower than the average, supporting a naive picture that the matter distribution in the direction perpendicular to the major axis is indeed in an underdense region such as a void. Kawaharada et al. (2010) has shown that, from X-ray observations, the galaxy cluster A1689, whose major axis is aligned with line-of-sight (Oguri et al., 2005; Umetsu et al., 2015), is likely to be surrounded by void regions in the sky plane (perpendicular to the line-of-sight). As the projection thickness increases, surface mass density profiles of aligned haloes () at large scales gradually become larger because matter clustered along the major axis starts to be projected. For the large projection thickness, e.g., , the results are quite similar to the one with full projection, suggesting that the projection thickness of used for our main analysis is sufficiently large. The results in Figure 5 also imply that an opening angle of the matter enhancement along the major axis at large three-dimensional separations from the haloes is quite wide, because the matter enhancement remains evident even at such large projected radii up to , almost the baryon acoustic oscillation scales.

### 5.2 Ellipse model

To gain a better understanding of the results we have shown so far, we here introduce a phenomenological model to explain the orientation dependence of surface mass density profiles. First recall that, without any loss of generality, the surface mass density profile around haloes is obtained from a line-of-sight projection of the halo-matter cross-correlation function, (Mandelbaum et al., 2006; Miyatake et al., 2015)

(19) |

and the cross-correlation function is originally given as a function of the three-dimensional separation , , where is the coordinate in the line-of-sight direction, and are coordinates in the two-dimensional plane perpendicular to the line-of-sight direction (-direction in our setting). In the above equation we assumed the axial symmetry around the light-of-sight direction, i.e., the cross-correlation is obtained from an angle -average of .

When stacking the matter distribution around haloes whose major axes are aligned with the line-of-sight direction, it violates the spherical symmetry even in a statistical sense. That is, the cross-correlation depends on the separation length and its direction:

(20) |

The persistence of the halo-orientation dependence even at large projection radii would suggest that the cross-correlation can be approximated by a multipole expansion up to the quadrupole, in analogous to the redshift-space distortion effect

(21) |

where . Our results suggest because we found an enhancement of the matter distribution along the major axis (i.e. ). The assumption that the large-scale matter distribution around the aligned haloes is described by a sum of the monopole (the first term) and quadrupole (the second term) might be too simplistic. However, a possible justification is as follows. We argue that such a large-scale cross-correlation between the halo scales, i.e the halo ellipticity at small scales of , and the large-scale structures up to arises from a mode-coupling between the small- and large-scale Fourier modes in nonlinear structure formation. As shown in the several work (e.g., Takada & Hu, 2013; Akitsu & Takada, 2017), the effect of large-scale modes on small-scale structures in the deeply nonlinear regime is well modeled by the second-derivatives of the long-wavelength gravitational potential due to the equivalence principle. If this is the case, the second-derivative tensor of the scalar potential causes only up to a quadrupole pattern in the small-scale matter distribution. We below test this hypothesis.

Motivated by the above picture and also extending the naive picture, equation (21), to a more general case, we assume that the anisotropic cross-correlation function is approximated by

(22) |

where is the spherically-symmetric, i.e., averaged over all directions, correlation function and is chosen to satisfy the following condition:

(23) |

Note that the right hand side of equation (23) is , not . This is because we need the mean radius of the ellipse as the argument in equation (22). Here, we define the radius such that the area of the circle is equal to that of the ellipse defined in equation (23). In order to satisfy this condition, the additional factor of is introduced in the right hand side of equation (23). The correlation function is given by two variables , and we will below treat the ellipticity as a free parameter at each separation, i.e., we allow to vary as a function of circular radius, . This is analogous to the elliptical halo profile in equation (17).

Once the above model (equation 22) is adopted, the projected surface mass density profile is given as

(24) |

where is the cosine angle between the major axes of haloes and the line-of-sight direction. Note that we ignored the mean background density contribution in equation (19) as it is not relevant for the following discussion. We test how this simplified model can reproduce the results we have shown.

In Figure 6, we show a contour map of the cross-correlation function in the polar coordinates converted from , measured from the -body simulations. We selected haloes whose virial masses are in the range of and the directional cosine with respect to line-of-sight satisfies . The line-of-sight (-direction) is taken to the angular direction . The fractional difference of the cross-correlation with respect to the angular averaged correlation function, i.e. , is shown in Figure 7. Figure 6 clearly shows the anisotropic nature of the cross-correlation such that the contour is elongated along the major axis of the halo, up to . Filaments and/or other haloes connected with filaments are more likely to exist along the major axis via the tidal field, which may lead to the elongation. It appears that the shape of the contour can be modeled by an ellipse up to the large scales (see also Schneider et al., 2012). In addition, the fractional difference shown in Figure 7 also displays a clear quadrupolar feature. In other words, this figure does not show a clear signature of higher-order multipole pattern beyond the quadrupole, supporting the argument we gave above.

In order to quantify this behavior, we fit the isoamplitude contour of by the elliptical model equation (22) for each polar radius with varying the parameter . In Figure 8, the best-fit axis ratio is shown as a function of the radius, for radial bin . The best-fit axis ratios slightly depend on the radius, but for a wide radius range the best-fit axis ratio is . For radial bins smaller than , the least chi-square method does not converge. In Figure 9, we show the comparison between the measured cross-correlation functions and the best-fit ellipse model predictions. We find that the ellipse model can explain the angular dependence well except the case of the largest radial bin, . That inconsistency may be caused by the inaccurate averaged cross-correlation measurements in simulations due to the limited size of the simulation box. In any case, we expect that this ellipse model is useful for analyzing the orientation dependence of surface mass density profiles.

## 6 Discussions

In this Section, we discuss potential impacts of our results on the
analysis of observational data. If the halo sample is selected without
any orientation bias, i.e., if the distribution of major axes is
completely random with respect to the line-of-sight direction, the
orientation dependence of surface density mass profiles is averaged out
so that we can simply ignore the orientation dependence derived in
this paper. For example, it is natural to expect that a sample of
normal galaxies is unbiased with respect to their orientation with
respect to the line-of-sight direction^{3}^{3}3The orientation bias
might be relevant even for a galaxy sample, e.g., if galaxies are
selected based on the aperture photometry with the aperture size
smaller than galaxy sizes. Our results presented in this paper
indicates that it is important to explicitly check the possible
orientation bias for any sample used for the analysis of surface
mass density profiles..
However, in some cases sample selections
inevitably introduce the orientation bias, and our results indicate
that correlation analysis of these samples can be significantly
affected by the strong orientation dependence of surface mass density
profiles. Here we discuss several cases where the orientation bias
can be significant.

### 6.1 Lensing selected clusters and galaxies

It has been known that cluster samples selected based on gravitational lensing exhibits the strong orientation bias. Hennawi et al. (2007) showed that clusters that exhibit strong lensing features tend to have their major axes aligned with the the line-of-sight direction, with the median angle of . Oguri & Blandford (2009) argued that the orientation bias is even stronger for strong lensing clusters with large Einstein radii (see also Meneghetti et al., 2010). Hamana et al. (2012) showed that the similar orientation bias exists also for weak lensing selected clusters, i.e., clusters selected from peaks in weak lensing mass maps. Our results predict that surface mass density profiles of these strong and weak lensing selected clusters have larger amplitudes at large radii as compared with those for other cluster samples. We note that the similar orientation bias is also known to exist for galaxy-scale strong lensing (Rozo et al., 2007).

### 6.2 Optically selected clusters

Recent wide-field photometric surveys with multiple passbands allow us to identify clusters of galaxies efficiently based on clustering of red-sequence galaxies. However, due to the projection of red-sequence galaxies along the line-of-sight, this method may preferentially select clusters whose major axes are aligned with the line-of-sight direction. Dietrich et al. (2014) explicitly studied the orientation bias of optically selected clusters using mock galaxy catalogs, and argued the existence of the orientation bias particularly for relatively less massive clusters.

Given the large number of optically selected clusters, high signal-to-noise measurements of surface mass density profiles with stacked weak lensing are possible, particularly in ongoing and future optical imaging surveys. Thus it is of great importance to accurately quantify the orientation bias of optically selected clusters in order to take proper account of the orientation dependence of the surface mass density profile for the accurate mass calibration.

Miyatake et al. (2016) claimed that clusters with member galaxies located closer to the center on average have large halo bias, which was ascribed to the detection of the assembly bias. However, if the average member galaxy separation correlates with the halo orientation, such a correlation can mimic the assembly bias signal and can explain the behavior observed by Miyatake et al. (2016). We check this possible correlation using our simulations. We focus only on haloes with masses in the range of . We derive the distribution of member galaxies in simulations using subhaloes with virial masses larger than . For each dark halo, we use member galaxies within from the halo center in the two dimensional space perpendicular to the line-of-sight direction, and along the line-of-sight direction. The former and latter ranges correspond to the projected cluster radius and the photometric redshift error, respectively. Thus, we compute the average distance of member galaxies, , in the same manner as in Miyatake et al. (2016). Figure 10 indicates that the average member galaxy separation depends very little on the halo orientation, which is insufficient to explain the large variation of the halo bias found in Miyatake et al. (2016). Hence, if the cluster sample used in Miyatake et al. (2016) suffers from the orientation bias, it would be likely from details in the cluster finding algorithm, rather than the simplified method using a catalog of subhaloes in each halo region.

### 6.3 Quasars

Clustering of quasars has been used to infer their host halo masses. However, given the small number of quasars, projected statistics have been used in most cases. Our results indicate that the orientation effect can be important for the clustering analysis of quasars if the jet axis of quasar is aligned with the halo orientation. For instance, Zhang et al. (2009) used -body simulations to show that the spin and minor axis of haloes are well aligned. Thus surface mass density profiles around quasars may also exhibit strong orientation dependence, if the spin axis of haloes is aligned with the jet axis.

If the argument above holds, it implies that type 2 (obscured) quasars have higher apparent clustering signals than type 1 (unobscured) quasars, even if their underlying host halo masses are similar. In fact, recent studies, using projected correlation functions as well as cosmic microwave background lensing measurements, have shown that type 2 quasars are more strongly clustered than type 1 quasars (e.g., Hickox et al., 2011; Donoso et al., 2014; DiPompeo et al., 2016). These observations have been used as a counterargument to the so-called unified model in which the difference between type 1 and type 2 quasars are explained solely by their viewing angles. Our results may imply that this issue needs to be revisited taking account of the possible orientation bias of quasars.

## 7 Conclusion

In this paper we have studied how the surface mass density profiles of haloes vary with the orientation of halo major axis with respect to the projection direction, using a suite of -body simulations. Dark haloes are known to be highly non-spherical by nature, and this non-sphericity leads to variations in the surface mass density profiles depending on the projection direction with respect to the halo shapes. We have found that, if projecting the matter distribution around haloes whose major axes are aligned with the line-of-sight (projection) direction, the average surface mass density profiles have larger amplitudes over all scales we studied, from to , than surface mass density profiles averaged over all orientations. The orientation bias at small projected separation in the one-halo regime is easily understood by the triaxiality of dark matter halo profile. However, we have shown that the orientation bias at large separations up to remains significant. This is opposite to what one would naively expect; since the mass distribution in the direction perpendicular to the major axis (i.e. along the minor or intermediate axis) tends to be in an underdense region such as void, we naively expect that the surface mass density projected at such a large separation is mainly from the underdense region, and therefore has smaller, instead of larger, amplitude than the average. We showed that the mass distribution in the direction of the major axis still affects the projected mass density profile even at large projection separations, and the opening angle of the enhanced mass distribution, viewed from the halo center, is quite wide.

Using the halo model, we have quantified how the halo-orientation dependence of surface mass density profiles causes a bias in parameters such as halo mass, concentration parameter and halo bias if those are estimated from the surface mass density profiles as in the weak lensing observables. We have found that the halo bias can be overestimated or underestimated by up to depending on the viewing angle with respect to the major axis. The orientation dependence of the halo bias from the two-halo term translates into the difference of the two-halo mass by up to a factor of two, which is much larger than the orientation dependence of the inferred halo mass from the one-halo term, which is at most.

In order to understand the origin of the orientation dependence at large scales, we have studied the average matter distribution around haloes, i.e. the halo-matter cross-correlation function, computed with selecting haloes whose major axes are aligned with line-of-sight, because the surface mass density profile is obtained from a projection of the halo-matter cross-correlation function at a fixed projected radius. The cross-correlation displays a clear quadrupole pattern in the mass distribution up to , with a significant opening angle of the mass enhancement along the major axis direction. The isoamplitude contour of the halo-matter cross-correlation at large scales is approximated by an ellipse with a axis ratio of . The strong cross-correlation up to large scales implies that the large-scale, anisotropic matter distribution, which is referred to as the large-scale “tidal” field, is the origin of halo shapes at small scales due to the mode coupling in nonlinear structure formation. It would be interesting to study whether the cross-correlation exists out to even larger separations such as baryon acoustic oscillation scales and even beyond remains or disappears, e.g. due to the mass conservation at very large scales. This is our future study, which will be presented elsewhere.

Our results have various implications. If a sample of large-scale structure tracers such as galaxies, clusters and quasars is affected by selection effects depending on the orientation of host haloes, the analysis of the sample assuming spherical symmetry can lead to a significant bias in estimating their halo masses or halo biases. We have discussed lensing and optically selected clusters as well as quasars for possible examples of samples with the significant orientation bias.

There are other observables that might be affected by the halo orientation bias. First, if a selection of spectroscopic galaxies is somehow affected by halo shapes or filaments in large-scale structure, the redshift-space clustering of galaxies would be modified by the halo orientation bias. It is recently shown that the large-scale anisotropic matter distribution, i.e. the tidal field, causes an anisotropic clustering pattern in the redshift-space galaxy distribution due to the nonlinear mode coupling (Akitsu et al., 2017; Akitsu & Takada, 2017) in addition to the redshift-space distortion and the cosmological distortion such as the Alcock-Paczynski effect. This could cause a bias in cosmological parameters estimated from the measured redshift-space clustering, if the orientation bias is ignored. Our results also imply that the peculiar velocity field in large-scale structure, which causes the redshift-space distortion, would have a strong correlation with shapes of haloes (e.g. Okumura et al., 2017). Another effect is the intrinsic alignment, which is one of the major systematic effects in weak lensing measurements. Our results show that the large-scale matter distribution has a strong correlation with shapes of haloes. If a sample of galaxies used in the shape measurements for weak lensing has a correlation with shapes of haloes, the sample would cause an intrinsic alignment contamination to the weak lensing measurement. These would be interesting to study, and our methods presented in this paper will be useful for such a study.

## Acknowledgements

We thank H. Miyatake, M. Shirasaki, R. Takahashi, and T. Oogi for useful discussions and K. Umestu for comments on an earlier version of the manuscript. KO is supported by Research Fellowships of the Japan Society for the Promotion of Science (JSPS) for Young Scientists and Advanced Leading Graduate Course for Photon Science. TO acknowledges support from the Ministry of Science and Technology of Taiwan under the grant MOST 106-2119-M-001-031-MY3. TN acknowledges financial support from Japan Science and Technology Agency (JST) CREST Grant Number JPMJCR1414. This work was supported by JSPS Grant-in-Aid for JSPS Research Fellow Grant Number JP16J01512. This work was supported in part by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan, and JSPS KAKENHI Grant Number JP26800093, JP23340061, JP26610058, JP15H03654, JP15H05887, JP15H05892, JP15H05893, JP15K21733, and 17K14273. Numerical simulations were carried out on Cray XC30 at the Center for Computational Astrophysics, National Astronomical Observatory of Japan.

## References

- Aihara et al. (2017) Aihara H., et al., 2017, preprint, (arXiv:1704.05858)
- Akitsu & Takada (2017) Akitsu K., Takada M., 2017, preprint, (arXiv:1711.00012)
- Akitsu et al. (2017) Akitsu K., Takada M., Li Y., 2017, Phys. Rev. D, 95, 083522
- Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, MNRAS, 367, 1781
- Baldauf et al. (2010) Baldauf T., Smith R. E., Seljak U., Mandelbaum R., 2010, Phys. Rev. D, 81, 063531
- Baltz et al. (2009) Baltz E. A., Marshall P., Oguri M., 2009, J. Cosmology Astropart. Phys., 1, 015
- Bartelmann & Schneider (2001) Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291
- Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109
- Bett (2012) Bett P., 2012, MNRAS, 420, 3303
- Bhattacharya et al. (2011) Bhattacharya S., Heitmann K., White M., Lukić Z., Wagner C., Habib S., 2011, ApJ, 732, 122
- Bonamigo et al. (2015) Bonamigo M., Despali G., Limousin M., Angulo R., Giocoli C., Soucail G., 2015, MNRAS, 449, 3171
- Brainerd et al. (1996) Brainerd T. G., Blandford R. D., Smail I., 1996, ApJ, 466, 623
- Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
- Busch & White (2017) Busch P., White S. D. M., 2017, MNRAS, 470, 4767
- Clampitt & Jain (2016) Clampitt J., Jain B., 2016, MNRAS, 457, 4135
- Clowe et al. (2004) Clowe D., De Lucia G., King L., 2004, MNRAS, 350, 1038
- Cooray & Sheth (2002) Cooray A., Sheth R., 2002, Phys. Rep., 372, 1
- Corless & King (2007) Corless V. L., King L. J., 2007, MNRAS, 380, 149
- Covone et al. (2014) Covone G., Sereno M., Kilbinger M., Cardone V. F., 2014, ApJ, 784, L25
- Dark Energy Survey Collaboration (2016) Dark Energy Survey Collaboration 2016, MNRAS, 460, 1270
- Dark Energy Survey Collaboration (2017) Dark Energy Survey Collaboration 2017, preprint, (arXiv:1708.01530)
- DiPompeo et al. (2016) DiPompeo M. A., Hickox R. C., Myers A. D., 2016, MNRAS, 456, 924
- Diemer & Kravtsov (2015) Diemer B., Kravtsov A. V., 2015, ApJ, 799, 108
- Dietrich et al. (2014) Dietrich J. P., et al., 2014, MNRAS, 443, 1713
- Donoso et al. (2014) Donoso E., Yan L., Stern D., Assef R. J., 2014, ApJ, 789, 44
- Duffy et al. (2008) Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64
- Dvornik et al. (2017) Dvornik A., et al., 2017, MNRAS, 468, 3251
- Evans & Bridle (2009) Evans A. K. D., Bridle S., 2009, ApJ, 695, 1446
- Gavazzi (2005) Gavazzi R., 2005, A&A, 443, 793
- Guzik & Seljak (2002) Guzik J., Seljak U., 2002, MNRAS, 335, 311
- Hamana et al. (2012) Hamana T., Oguri M., Shirasaki M., Sato M., 2012, MNRAS, 425, 2287
- Hennawi et al. (2007) Hennawi J. F., Dalal N., Bode P., Ostriker J. P., 2007, ApJ, 654, 714
- Hickox et al. (2011) Hickox R. C., et al., 2011, ApJ, 731, 117
- Hikage et al. (2013) Hikage C., Mandelbaum R., Takada M., Spergel D. N., 2013, MNRAS, 435, 2345
- Jing & Suto (2002) Jing Y. P., Suto Y., 2002, ApJ, 574, 538
- Joachimi et al. (2015) Joachimi B., et al., 2015, Space Sci. Rev., 193, 1
- Kawaharada et al. (2010) Kawaharada M., et al., 2010, ApJ, 714, 423
- Leauthaud et al. (2012) Leauthaud A., et al., 2012, ApJ, 744, 159
- Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
- Limousin et al. (2013) Limousin M., Morandi A., Sereno M., Meneghetti M., Ettori S., Bartelmann M., Verdugo T., 2013, Space Sci. Rev., 177, 155
- Lin et al. (2016) Lin Y.-T., Mandelbaum R., Huang Y.-H., Huang H.-J., Dalal N., Diemer B., Jian H.-Y., Kravtsov A., 2016, ApJ, 819, 119
- Lombriser et al. (2012) Lombriser L., Schmidt F., Baldauf T., Mandelbaum R., Seljak U., Smith R. E., 2012, Phys. Rev. D, 85, 102001
- Mandelbaum et al. (2005) Mandelbaum R., Tasitsiomi A., Seljak U., Kravtsov A. V., Wechsler R. H., 2005, MNRAS, 362, 1451
- Mandelbaum et al. (2006) Mandelbaum R., Seljak U., Kauffmann G., Hirata C. M., Brinkmann J., 2006, MNRAS, 368, 715
- Mandelbaum et al. (2013) Mandelbaum R., Slosar A., Baldauf T., Seljak U., Hirata C. M., Nakajima R., Reyes R., Smith R. E., 2013, MNRAS, 432, 1544
- Meneghetti et al. (2010) Meneghetti M., Fedeli C., Pace F., Gottlöber S., Yepes G., 2010, A&A, 519, A90
- Miyatake et al. (2015) Miyatake H., et al., 2015, ApJ, 806, 1
- Miyatake et al. (2016) Miyatake H., More S., Takada M., Spergel D. N., Mandelbaum R., Rykoff E. S., Rozo E., 2016, Physical Review Letters, 116, 041301
- More et al. (2015) More S., Miyatake H., Mandelbaum R., Takada M., Spergel D. N., Brownstein J. R., Schneider D. P., 2015, ApJ, 806, 2
- More et al. (2016) More S., et al., 2016, ApJ, 825, 39
- Murata et al. (2017) Murata R., Nishimichi T., Takada M., Miyatake H., Shirasaki M., More S., Takahashi R., Osato K., 2017, preprint, (arXiv:1707.01907)
- Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
- Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
- Nishimichi et al. (2009) Nishimichi T., et al., 2009, PASJ, 61, 321
- Nishimichi et al. (2010) Nishimichi T., Taruya A., Koyama K., Sabiu C., 2010, J. Cosmology Astropart. Phys., 7, 002
- Oguri & Blandford (2009) Oguri M., Blandford R. D., 2009, MNRAS, 392, 930
- Oguri & Hamana (2011) Oguri M., Hamana T., 2011, MNRAS, 414, 1851
- Oguri & Takada (2011) Oguri M., Takada M., 2011, Phys. Rev. D, 83, 023008
- Oguri et al. (2003) Oguri M., Lee J., Suto Y., 2003, ApJ, 599, 7
- Oguri et al. (2005) Oguri M., Takada M., Umetsu K., Broadhurst T., 2005, ApJ, 632, 841
- Oguri et al. (2010) Oguri M., Takada M., Okabe N., Smith G. P., 2010, MNRAS, 405, 2215
- Oguri et al. (2012) Oguri M., Bayliss M. B., Dahle H., Sharon K., Gladders M. D., Natarajan P., Hennawi J. F., Koester B. P., 2012, MNRAS, 420, 3213
- Okabe et al. (2010) Okabe N., Takada M., Umetsu K., Futamase T., Smith G. P., 2010, PASJ, 62, 811
- Okabe et al. (2013) Okabe N., Smith G. P., Umetsu K., Takada M., Futamase T., 2013, ApJ, 769, L35
- Okumura et al. (2017) Okumura T., Nishimichi T., Umetsu K., Osato K., 2017, preprint, (arXiv:1706.08860)
- Planck Collaboration (2016) Planck Collaboration 2016, A&A, 594, A13
- Prat et al. (2017) Prat J., et al., 2017, preprint, (arXiv:1708.01537)
- Rozo et al. (2007) Rozo E., Chen J., Zentner A. R., 2007, preprint, (arXiv:0710.1683)
- Schneider et al. (2012) Schneider M. D., Frenk C. S., Cole S., 2012, J. Cosmology Astropart. Phys., 5, 030
- Sereno et al. (2015) Sereno M., Veropalumbo A., Marulli F., Covone G., Moscardini L., Cimatti A., 2015, MNRAS, 449, 4147
- Sheldon et al. (2004) Sheldon E. S., et al., 2004, AJ, 127, 2544
- Shin et al. (2017) Shin T.-h., Clampitt J., Jain B., Bernstein G., Neil A., Rozo E., Rykoff E., 2017, preprint, (arXiv:1705.11167)
- Springel (2005) Springel V., 2005, MNRAS, 364, 1105
- Takada & Hu (2013) Takada M., Hu W., 2013, Phys. Rev. D, 87, 123504
- Takada & Jain (2003a) Takada M., Jain B., 2003a, MNRAS, 340, 580
- Takada & Jain (2003b) Takada M., Jain B., 2003b, MNRAS, 344, 857
- Tinker et al. (2010) Tinker J. L., Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010, ApJ, 724, 878
- Troxel & Ishak (2015) Troxel M. A., Ishak M., 2015, Phys. Rep., 558, 1
- Umetsu et al. (2014) Umetsu K., et al., 2014, ApJ, 795, 163
- Umetsu et al. (2015) Umetsu K., et al., 2015, ApJ, 806, 207
- Valageas & Nishimichi (2011) Valageas P., Nishimichi T., 2011, A&A, 527, A87
- Vega-Ferrero et al. (2017) Vega-Ferrero J., Yepes G., Gottlöber S., 2017, MNRAS, 467, 3226
- Vera-Ciro et al. (2014) Vera-Ciro C. A., Sales L. V., Helmi A., Navarro J. F., 2014, MNRAS, 439, 2863
- Zhang et al. (2009) Zhang Y., Yang X., Faltenbacher A., Springel V., Lin W., Wang H., 2009, ApJ, 706, 747
- van Uitert et al. (2017a) van Uitert E., et al., 2017a, preprint, (arXiv:1706.05004)
- van Uitert et al. (2017b) van Uitert E., et al., 2017b, MNRAS, 467, 4131