© The Authors 2025

1. Introduction

In recent years, cosmic shear measurements have demonstrated a capacity to constrain the dark-matter distribution on cosmological scales (Dark Energy Survey and Kilo-Degree Survey Collaboration 2023; More et al. 2023). Over the coming decade, they represent a key probe in determining our Universe’s cosmological scenario at later times and to gain insights into the nature of dark energy and dark matter at lower and intermediate redshift. Weak lensing is one of the primary cosmological probes for a number of upcoming surveys, including Euclid (Euclid Collaboration: Mellier et al. 2025), Vera Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019), and Nancy Grace Roman Space Telescope (Roman; Eifler et al. 2021). The cosmic shear signal is estimated by averaging over the shape of background galaxies, which are coherently deformed as their light travels through the large scale structures. By studying the correlations between the shapes of neighbouring galaxies, we can infer the foreground distribution of matter, particularly dark matter. In the coming decade, unprecedented statistical power will be reached by modern weak lensing experiments, which are poised to observe billions of galaxies. This increase in statistical power means that systematic contributions, which had previously been of lesser importance, can lead to important biases if they are not accounted for. These systematic contributions can be broadly split into two categories: biases at the measurement level and biases at the modelling level. In the latter, three sources of biases are particularly critical: baryonic feedback, a non-linear effect that suppresses mode fluctuations in the range of k ∼ [0.1, 10] Mpc/h, photometric-redshift biases caused by a poor estimation of the galaxies redshift distribution, and galaxy intrinsic alignment; the latter is due to the fact that galaxy orientations are not uniformly distributed on the sky, but correlated with their environments. For a review, we refer to Mandelbaum (2018).

It was shown in Catelan et al. (2001) and Crittenden et al. (2001) that tidal forces arising from the density field have the potential to produce coherent orientations in neighboring galaxies, which can pose a challenge for cosmic shear. Specifically, the difficulty lies in disentangling the correlations induced directly on the galaxy shapes by the tidal forces from those caused by the shearing of the light profile due to gravitational lensing. The effects of tidal forces on the orientation of galaxies have been described in Catelan et al. (2001) and in Mackey et al. (2002). For a review and visual explanations, we refer to Joachimi et al. (2015), Lamman et al. (2024).

To mitigate the intrinsic alignment, various solutions have been proposed. To remove the intrinsic alignment from the cosmic-shear signal, two main methods exist: ‘down-weighting’ (King & Schneider 2002) and ‘nulling’ (Heavens & Joachimi 2011). Down-weighting consists of suppressing the correlation between pairs of galaxies that are close in redshift and subject to intrinsic alignment. This is achieved by applying a Gaussian weighting factor in the correlation function estimator to suppress the contributions from close pairs. This weighting factor has a variance reflecting the photometric redshift uncertainty, which implies that imprecise photometric redshifts will lead to a large loss of statistical power, as many pairs are suppressed. The concept of ‘nulling’ uses the specific dependency of the lensing kernel on the comoving distance to suppress any combination of bins that might be too sensitive to the interplay between shear and intrinsic alignment. For a review, we refer to Kirk et al. (2015).

These methods, however, have been shown to strongly reduce the statistical power of a weak lensing survey and to be very dependent on precise photometric redshifts. The approach taken in recent stage-III cosmic shear surveys has therefore consisted of the joint modelling of the weak lensing and the intrinsic alignment signal. Two main models have been applied: the non-linear alignment model (NLA) and its higher-order expansion, the tidal alignment and tidal torque model (TATT). These models are described in Sects. 4.2.1 and 4.2.2. They are theoretically motivated, but the physical range of their validity remains an open question. In particular, Bakx et al. (2023) recently showed that the NLA (TATT) prediction of the alignment of subhaloes in the DarkQuest simulation is only able to describe correlations on modes smaller than k ≤ 0.08 h/Mpc (k ≤ 0.15 h/Mpc). This shortcoming has motivated the inclusion of even higher order terms in the effective field theory approach. With eight free parameters (Vlah et al. 2020), this model is harder to constrain than the NLA and TATT models, which have one and three parameters, respectively, when the redshift and luminosity dependencies are fixed. A similar EFT approach can be found in Chen & Kokron (2024). We note that while these theoretically motivated models are essential in the prediction of the weak lensing correlation functions, they are also used to add intrinsic alignment to convergence maps used in simulation-based approaches, which require the full-field information (Harnois-Déraps et al. 2021; Kacprzak et al. 2023; Jeffrey et al. 2024). Similarly, this model has been applied in a direct inference of the galaxies’ alignment with the tidal field, which is itself inferred from SDSS spectroscopic tracers (Tsaprazi et al. 2022). Other models based on the halo model (Schneider & Bridle 2010; Fortuna et al. 2021a) achieve higher accuracy at smaller scales by incorporating an intrinsic ellipticity field within each halo. These models predict distinct effects on central and satellite galaxies, with probability distributions derived from halo-occupation distribution models based on an NFW profile and dependent on halo mass.

To properly characterise the intrinsic alignment, various hydrodynamical simulations have been used (Chisari et al. 2015; Samuroff et al. 2021). While they offer a picture that is broadly consistent with observations, no precise calibration has been reached from these various simulations, as the selection of galaxies and their mapping to phenomenological properties remains challenging. In Tenneti et al. (2016), it was shown that the choice of elliptical galaxies can produce drastic changes, as the ellipticity-density correlation can go from positively to negatively correlated depending on the simulation. A more novel route to calibrating the intrinsic alignment, which makes use of the power of N-body simulations to describe the large-scale structure of the Universe, is to add galaxies empirically. This is performed by drawing misalignment angles from probability distribution functions, which are calibrated on observations or hydrodynamical simulations and draw central and satellite galaxy populations from a halo occupation distribution model (Hoffmann et al. 2022; Van Alfen et al. 2024).

In parallel, and complementary to the simulation efforts, many studies have tried to constrain intrinsic alignment models by direct measurements, which require precise redshift information to disentangle lensing contributions from intrinsic alignment. They have shown a strong intrinsic alignment between red galaxies at scales relevant for the two-halo term (i.e. ∼6 − 150 Mpc) outside the one-halo regime (Joachimi et al. 2013; Singh et al. 2015, hereafter Si15; Johnston et al. 2019; Fortuna et al. 2021b; Samuroff et al. 2018). For red galaxies, the intrinsic alignment was shown to depend on luminosity following a broken power law, (Fortuna et al. 2021a; Samuroff et al. 2023, hereafter Sa23). For blue spiral galaxies, direct measurements tend to prefer an absence of intrinsic alignment (Mandelbaum et al. 2011; Johnston et al. 2019; Tonegawa et al. 2024; Sa23). The alignment between blue and red samples was adequately measured in Johnston et al. (2019), was made possible by the very densely sampled and complete spectroscopic galaxy and mass assembly survey (GAMA, Hopkins et al. 2013), overlapping with the KiDS dataset (Giblin et al. 2021). Very recently, McCullough et al. (2024) showed a strong difference in cosmological parameter estimation from cosmic shear when splitting the DES Y3 in red and blue galaxy samples. This difference was attributed to the varying sensitivities of these samples to the intrinsic alignment models.

This work contributes to the challenge of propagating intrinsic alignment measurements into a sensible prior for stage-IV surveys. This task is complicated by a few subtleties affecting observations and the resulting amplitude of intrinsic alignment. One such effect is the dependence of the intrinsic alignment amplitude on the shape measurement method: as tidal forces act stronger on the outskirts of galaxies, these extended regions are more affected by intrinsic alignment than the central part. This was first observed in Singh & Mandelbaum (2016), who found a 40% amplitude difference between the re-Gaussianisation, isophotal, and de-Vaucouleurs shape measurement methods. While these three shape measurement methods are less frequently used in stage-III surveys, a recent work (MacMahon-Gellér & Leonard 2024) showed how this effect can be used to gain information on intrinsic alignment effects. In this work, the authors applied the multi-estimator method (MEM) described in Leonard et al. (2018), which uses the different deformations in isophotes between intrinsic alignment and cosmic shear to disentangle the two effects. They found that a 40% difference in the intrinsic alignment amplitude is required for a 1σ detection in LSST Year 1 data. In Georgiou et al. (2019a,b), a strong change in IA amplitude was shown by varying the weight function of the shape measurement. In addition Georgiou et al. (2019a) noted a chromatic effect in intrinsic alignment by changing the photometric band in which the galaxy shape is measured.

In this work, we applied the projected two-point correlation function, but other estimators have recently been developed. In Linke et al. (2024), a three-point correlation function with LOWZ data was found to be in good agreement with the two-point correlation functions, promising better constraining power for future surveys. Another important avenue is the three-dimensional measurement in Fourier space, known as the Yamamoto estimator (Yamamoto et al. 2006). It was initially used by the galaxy clustering community, but it has been applied successfully to constrain intrinsic alignment (Kurita et al. 2021; Kurita & Takada 2023). Similarly, Singh et al. (2024) showed that including real-space multipoles can also help in reaching tighter constraints.

While intrinsic alignment is a treated as a source of contamination with respect to to weak lensing measurements, it is also a cosmological probe in its own right. Examples include the detection of baryon acoustic oscillations (BAOs) in the intrinsic-alignment correlation functions (van Dompseler et al. 2023; Xu et al. 2023a) and constraints on primordial non-Gaussianities (Schmidt et al. 2015; Chisari et al. 2016; Kurita & Takada 2023).

In this work, we used the UNIONS imaging survey together with BOSS/eBOSS spectroscopic galaxies to measure the intrinsic alignment. We compared our measurements to other samples by investigating the dependence on luminosity before taking a deeper look at intrinsic alignment measurements from SDSS lensing and DES Y3 lensing on the same set of galaxies. We verified that our measurements are robust against systematic diagnostics, in particular, with respect to the point spread function (PSF) multiplicative and additive biases.

This paper is organised as follows. We start in Sect. 2 by presenting the datasets used to constrain intrinsic alignment models. In Sect. 3, we we present the various estimators correlating shapes and positions. Sect. 4 introduces the NLA and TATT models and the calculations of the predictions for the estimators. We also present the calculation of the analytical covariance matrices, which we validate with jackknife estimates. We give our results in Sect. 5 and preset our conclusions in Sect. 6.

2. Data

2.1. UNIONS catalogues

We used 3500 deg² of imaging data from UNIONS, the Ultraviolet Near-Infrared Optical Northern Survey (Gwyn et al. 2025). In Fig. 1 we show the different footprints. UNIONS groups deep-field observations obtained by the CFHT (Canada-France-Hawai’i Telescope on Mauna Kea), Pan-STARRS (Panoramic Survey Telescope and Rapid Response System on Maui), and Subaru (on Mauna Kea) telescopes in the northern hemisphere. The final observations will contain multi-band (ugriz) information to provide precise photometric redshifts for Euclid. In this work, we use shape measurements realised in the r band observed by CFHT. The astrometry and image reduction were obtained with the Megapipe pipeline (Gwyn 2008). Throughout this work, we use a catalogue created with the ShapePipe¹ pipeline (Farrens et al. 2022). A first realisation of the catalogue covering 1500 deg² and systematics diagnostics are described in Guinot et al. (2022). The Shapepipe pipeline uses the ngmix (Sheldon 2015) software for shape measurement, with a 2D Gaussian profile to fit the galaxy light distribution, and Metacalibration (Sheldon & Huff 2017) for the shear calibration. Metacalibration is a framework in which galaxies are artificially sheared during the shape measurement process to evaluate a finite ∂ε/∂γ correction. This gives a shear response, R, evaluated between positively and negatively sheared postage stamps via finite differences. This artificial shearing makes it possible to quantify the model, noise, and selection biases by evaluating how the estimated ellipticity varies when it is changed by a known shear. Objects were detected and deblended with the SExtractor software (Bertin & Arnouts 1996). In this analysis, we included objects with FLAGS ≤ 2 from SExtractor as we are working with spectroscopically confirmed galaxies only. The catalogue comes with inverse variance weights related to errors on the shape measurements given by ngmix,

$$\begin{array}{c}\hfill w_{\mathrm{SP}}=\frac{1}{2{\sigma }_{\mathrm{SN}}^2+{\sigma }_{e_1}^2+{\sigma }_{e_2}^2},\end{array}$$(1)

Fig. 1. Distribution of the different spectroscopic and lensing samples used in this work. The plot on the left shows the North Galactic Cap (NGC), where most UNIONS galaxies lie. On the right is a portion of the South Galactic Cap (SGC), with the DES Y3 observations in the most southern part. There is no overlap between the UNIONS and DES Y3 catalogues. The apparent variation in shades between the two panels is due to the two projections resulting in different point densities.

where σ_SN is the galaxy shape noise and σ_ei represents the variance of the ellipticity component i = 1, 2 estimated by the model fitting procedure.

The PSF was modelled with MCCD² (Liaudat et al. 2021), a software based on a minimisation scheme analogous to PSFEx, but which includes both a ‘local’ component fitted on a single CCD and a ‘global’ contribution. This coupling allows it to capture effects varying across the entire focal plane.

The objects entering the lensing catalogue are selected with a signal-to-noise ratio of S/N > 10, estimated from the fitted flux divided by the flux error. The UNIONS catalogue reaches an r-magnitude depth of r ≈ 24.7, as shown in Fig. 3. For validation, we compared a correlation function measurement obtained with ShapePipe to one produced with lensfit (shown in Appendix B).

Fig. 3. Un-weighted normalised r-band magnitude distributions of the BOSS/eBOSS spectroscopic samples and UNIONS photometric galaxies. We can see that the spectroscopic galaxies from the BOSS/eBOSS surveys fall into the brighter tail of the UNIONS distribution. The galaxies which dominate the lensing signal are comparatively fainter. Note: we could not use our large overlap with LOWZ galaxies as these bright galaxies (r < 20) are mostly absent from our lensing sample due to size cuts.

2.2. DES Y3 catalogue

For a comparison between UNIONS and DES intrinsic alignment measurements, we reproduced the Dark Energy Survey (DES) measurements close to those presented in Sa23. For this purpose, we used the publicly available DES Y3 Metacalibration shape catalogue. The catalogue is described in Gatti et al. (2021). It contains 100 204 026 galaxies covering 4139 square degrees of sky area. Images were taken in the g, r, i, z, and Y bands and galaxy shapes were fitted on r, i, z band images simultaneously.

The core of the galaxy shear estimation methods is common to DES and UNIONS, which both use the ngmix model fitting software and calibrate the shear response with the Metacalibration scheme. This is critical for our comparison, as it was shown in Georgiou et al. (2019a) and Singh & Mandelbaum (2016) that different shape measurement algorithms can lead to very different intrinsic alignment amplitudes.

To create the shear and redshift galaxy sample, we first matched the BOSS and eBOSS spectroscopic galaxies to the full DES Y3 source galaxies (see Sect. 2.4). We then removed spurious objects by applying cuts based on the size and S/N, as well as the likelihood of being contaminated by binary stars, as described in Gatti et al. (2021). Below we highlight differences between the UNIONS and DES Y3 catalogues and methods which are important in order to compare the measurements from each dataset:

• DES uses the PIFF (Jarvis et al. 2021) PSF model while the UNIONS catalogue is built with MCCD. PIFF in the DES Y3 setup fits the PSF CCD by CCD, while MCCD captures the entire focal plane at once.

• UNIONS measures galaxy shapes by fitting the CFHT r band, while DES provides joint fits using the r, i, and z band images.

• The seeing between the two data sets is very different; UNIONS carried out observations with a median seeing of 0.69″ in r, while DES images have a median seeing of 0.95″, 0.88″, 0.83″ in r, i, and z, respectively (Abbott et al. 2021).

2.3. BOSS/eBOSS catalogues

The BOSS/eBOSS tracer catalogues were selected for galaxy clustering studies. Using these catalogues to measure the intrinsic alignment signal is therefore advantageous, as these galaxies have been selected by phenomenological properties and have a well-characterised galaxy bias. Such homogeneous galaxy samples are not representative of the complete population of a typical weak-lensing galaxy set and the corresponding intrinsic alignment measurements cannot be adopted directly as prior for weak lensing. To measure the intrinsic alignment, we need to model the galaxy bias parameters, over which we perform the marginalisation. Therefore, these parameters need to be well constrained to avoid degeneracies with intrinsic alignment parameters.

The catalogue contains weights for each galaxy described in Reid et al. (2015), Raichoor et al. (2020), which are the product of the following components:

• w_cp, to account for fiber collisions in galaxy pairs;

• w_noz, to account for model redshift failures;

• w_sys, tot, to remove spurious observationally induced fluctuations of non-cosmological origin;

• w_FKP, to minimise the mean square fluctuations of the density power spectrum, P_δ(k), introduced in Feldman et al. (1994). They are defined as ${\mathit{w}}_{\mathrm{FKP}}(z)={(1+\overline{n}(z)P_0)}^{-1}$ with P₀ being the value of the power spectrum at k ≈ 0.15 h⁻¹ Mpc and $\overline{n}(z)$ the mean galaxy density at a redshift, z.

In the following, we briefly describe the different BOSS/eBOSS samples used in this paper. The redshift distribution of the different samples is shown in Fig. 2.

Fig. 2. Un-weighted normalised redshift distribution n(z) of the three BOSS/eBOSS galaxy samples. The UNIONS n(z) is approximate and for reference only, as blinding was implemented as a shift in the redshift distribution. The n(z) has been obtained by matching galaxies in the CFHTLenS W3 field and calibrated through deep spectroscopic reference samples.

2.3.1. CMASS

The eBOSS CMASS sample (Dawson et al. 2013) was selected on stellar masses as a passively evolving galaxy sample. It contains galaxies in the redshift range z ∈ [0.4, 0.8]. The targets were selected using colour cuts for which the philosophy and functioning are presented in Eisenstein et al. (2001). These cuts are i < 19.9, d_⊥ > 0.55 and i < 19.86 + 1.6 × (d_⊥ − 0.8). The quantity d_⊥ = (r − i)−(g − r)/8 quantifies a line in colour space perpendicular to the position of the locus of the desired galaxy population at z > 0.4. For this sample, the combined weight w is

$$\begin{array}{c}\hfill w=(w_{\mathrm{cp}}+w_{\mathrm{noz}}-1)w_{\mathrm{sys},\mathrm{tot}}{w}_{\mathrm{FKP}}.\end{array}$$(2)

2.3.2. LRG

The eBOSS luminous red galaxy (LRG) DR16 sample (Bautista et al. 2018) was selected with the conditions r − i > 0.98, and r − W1 > 2.0 (r − i), where W1 is the 3.4 μm band from the WISE survey, and the r and i bands are from the SDSS imaging survey. This results in a sample of high-redshift red galaxies with z ∈ [0.6, 1.1].

The total weight for eBOSS LRGs is (Bautista et al. 2020):

$$\begin{array}{c}\hfill w=w_{\mathrm{cp}}{w}_{\mathrm{noz}}{w}_{\mathrm{sys},\mathrm{tot}}{w}_{\mathrm{FKP}}.\end{array}$$(3)

2.3.3. ELGs

The eBOSS emission-line galaxy (ELG) sample (Raichoor et al. 2020) is made up of selected ELGs in the range 0.6 < z < 1.1 over an observed sky area of 1170 deg². The target selection was made via a cut in the g band magnitude, which is sensitive to [OII] emitters, and a box selection in the grz colour space. In the NGC, the g cut is 21.825 < g < 22.9. These galaxies form a younger and bluer sample, which means that no intrinsic alignment is expected in these galaxies and that they have a lower clustering bias. The ELG galaxies follows the weight scheme Eq. (3) adopted in eBOSS for the LRGsample.

2.4. Matching procedure

To obtain a sample with precise redshift information, we matched the publicly available SDSS catalogues described above with the UNIONS shape catalogue introduced in Sect. 2.1. We used the k-D tree implementation in astropy with a tolerance of 1″, which is the maximum distance between galaxy pairs to be considered a match. We obtained the matched sub-sample of our respective catalogues with precise ellipticity and redshift measurements as our shape sample. For the density sample, we use all galaxies of the respective SDSS catalogues. The different properties of the surveys are described in Table 1.

3. Correlation and covariance estimations

3.1. Shape-density correlation estimator

To estimate the shape-density correlations of intrinsic galaxy alignment, we use the common projected two-point cross-correlation function w_g, + (Mandelbaum et al. 2006). The intrinsic-alignment measurement is similar to galaxy-galaxy lensing, but differs from the latter as the shape tracer can be either at lower or higher redshift compared to the density tracer, at very limited separations along the line of sight.

Following Mandelbaum et al. (2006), we first computed the 3D correlation function using the modified Landy-Szalay estimator:

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{g}+}(r_{\mathrm{t}},\mathrm{\Pi })=\frac{1}{N_{\mathrm{r},{\mathrm{r}}_{\mathrm{s}}}(r_{\mathrm{t}},\mathrm{\Pi })}& \times ({\displaystyle \sum _{\mathit{sd}}}{\gamma }_{+,sd}{\mathrm{\Delta }}_{r,\mathrm{\Pi }}({\mathit{x}}_s-{\mathit{x}}_d)\hfill \\ \hfill & -{\displaystyle \sum _{\mathit{sr}}}{\gamma }_{+,sr}{\mathrm{\Delta }}_{r_{\mathrm{t}},\mathrm{\Pi }}({\mathit{x}}_s-{\mathit{x}}_r)).\hfill \end{array}$$(4)

The variables s, d, and r indicate indices of objects in the shape, density, and random sample, respectively. The estimator is a sum of the tangential shear, γ_+, sd, of objects in the shape sample around the density sample, corrected by a contribution of γ_+, sr around random points. The normalizing quantity N_r, rs(r_t, Π) is the number of pairs between the random position and random shape samples. The pairs are binned in the tangential separations r_t and line-of-sight distance Π. The subtraction of the random correlations is introduced to mitigate geometrical effects, for example, the shape correlation with the survey mask.

The bin indicator function Δ of the 2D separation vector x = (x_⊥,x_∥) is defined as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{r,\mathrm{\Pi }}(\mathit{x})=\{\begin{array}{c}\begin{array}{cc}1\hfill & \mathrm{if}|log(x_{\perp }/r_{\mathrm{t}})|<\frac{\mathrm{\Delta }logr}{2}\mathrm{and}|x_{\Vert }-\mathrm{\Pi }|<\frac{\mathrm{\Delta }\mathrm{\Pi }}{2},\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\hfill \end{array}\end{array}$$(5)

With this characteristic function, we can separate galaxy pairs into logarithmic perpendicular bins with width Δlog r, and linear radial (line-of-sight) bins with a width, ΔΠ.

Next, we integrate the correlation function Eq. (4) along the line of sight to project out redshift space distortion (RSD) effects stemming from peculiar velocities,

$$\begin{array}{c}\hfill w_{\mathrm{g}+}(r_{\mathrm{t}})={\displaystyle {\int }_{-{\mathrm{\Pi }}_{\max }}^{{\mathrm{\Pi }}_{\max }}}{\xi }_{\mathrm{g}+}(r_{\mathrm{t}},\mathrm{\Pi })\mathrm{d}\mathrm{\Pi }.\end{array}$$(6)

We varied the comoving limit and found a stable value of Π_max = 150 Mpc, which is consistent with previous works (Si15, Sa23). In practice, we use the Treecorr package described in Jarvis (2015) to compute Eq. (4) in slices of ΔΠ = 15 Mpc, and sum up the correlations to obtain w_g+.

3.2. Density-density correlation estimator

To estimate the galaxy bias parameters, as motivated in Sect. 4.1, we computed the two-point galaxy correlation function using the Landy-Szalay estimator:

$$\begin{array}{c}\hfill {\xi }_{\mathrm{gg}}(r_{\mathrm{t}},\mathrm{\Pi })=\frac{N_{\mathit{dd}}(r_{\mathrm{t}},\mathrm{\Pi })-2N_{\mathit{dr}}(r_{\mathrm{t}},\mathrm{\Pi })+N_{\mathit{rr}}(r_{\mathrm{t}},\mathrm{\Pi })}{N_{r,r}(r_{\mathrm{t}},\mathrm{\Pi })}·\end{array}$$(7)

As for the shape-density estimator, the random correlations account for systematic geometrical effects. We integrate along the line of sight to wash out RSD effects:

$$\begin{array}{c}\hfill w_{\mathrm{gg}}(r_{\mathrm{t}})={\displaystyle {\int }_{-{\mathrm{\Pi }}_{\max }}^{{\mathrm{\Pi }}_{\max }}}{\xi }_{\mathrm{gg}}(r_{\mathrm{t}},\mathrm{\Pi })\mathrm{d}\mathrm{\Pi }.\end{array}$$(8)

3.3. Shape-shape correlation estimator

The majority of the intrinsic alignment constraining power is obtained from w_g+, as position-shape measurements have a higher S/N than the shape-shape measurements. While w_g+ is sensitive to gravitational-intrinsic correlations, as explained in Sect. 4.2.1, the projected shape-shape correlation function w₊₊ is proportional to the intrinsic-intrinsic power spectrum. Previous measurements were presented in Si15, Sa23.

The corresponding estimator is defined as

$$\begin{array}{c}\hfill {\xi }_{++}(r_{\mathrm{t}},\mathrm{\Pi })=\frac{1}{N_{{\mathrm{r}}_{\mathrm{s}},{\mathrm{r}}_{\mathrm{s}}}(r_{\mathrm{t}},\mathrm{\Pi })}{\displaystyle \sum _{ss\prime }}{\xi }_{++,ss\prime }{\mathrm{\Delta }}_{r_{\mathrm{t}},\mathrm{\Pi }}({\mathit{x}}_s-{\mathit{x}}_{s\prime }),\end{array}$$(9)

which is then integrated along the line of sight,

$$\begin{array}{c}\hfill w_{++}={\displaystyle {\int }_{-{\mathrm{\Pi }}_{\max }}^{{\mathrm{\Pi }}_{\max }}}{\xi }_{++}(r_{\mathrm{t}},\mathrm{\Pi })\mathrm{d}\mathrm{\Pi }.\end{array}$$(10)

We included this estimator for completeness and for comparison with previous works.

3.4. Computing luminosities

When estimating the intrinsic luminosities, various observational effects that modify the observed galaxy flux have to be corrected. For this we compute k-correction using (Blanton & Roweis 2007) as implemented in the python package k-correct. Furthermore, the e-correction is a quantity that reflects the evolutionary processes of galaxies, such as star formation, aging stellar populations, and mergers of galaxies, which change the galaxy brightness as the universe evolves. We obtained the e-correction from the Ezgal package (Mancone & Gonzalez 2012), for which we assumed the commonly used Bruzual & Charlot (2003) model as the stellar population synthesis model. Lastly, we corrected for the dust extinction with the E-corrections. These values are given by the SDSS catalogue and are based on the Schlegel et al. (1998) maps. To compute the luminosities, we used the DevFlux apertures from the SDSS photometric imaging survey; this allowed us to get an approximation of the SED necessary for the e-correction. The absolute magnitude as a function of the observed SDSS magnitude, r, was computed with the formula

$$\begin{array}{c}\hfill M_r=r-5[{log}_{10}\left(D_{\mathrm{l}}(z)\right)-1]-(k+e)+E,\end{array}$$(11)

where D_l is the luminosity distance in pc/h. This equation and the correction quantities are well established for red galaxy populations and cannot be generically applied, as they require a good understanding of galaxy evolution. Luminosity mean values for different samples are shown in Table 1. We use this information to split our CMASS galaxies in luminosity sub-samples and compare the luminosity scaling to results obtained from the literature in Sect. 5.3. Our luminosity estimates are in good agreement with those from Sa23, even though we do not expect exact matches as their measurements are based on the DES r band filter, which differs slightly from SDSS. To obtain ⟨L/L₀⟩, we used the standard pivot luminosity, L₀, defined to match an absolute magnitude of M_r = −22.

4. Modeling intrinsic alignment and galaxy clustering

4.1. Modelling Galaxy clustering

Galaxies are biased tracers of the underlying dark matter field. At the linear level, the relation is simply δ_g(x) = b₁δ_m(x), which consequently is for the matter power spectrum

$$\begin{array}{c}\hfill P_{\mathrm{gg}}=b_1^2{P}_{\delta },\end{array}$$(12)

where P_δ represents the linear matter power spectrum. We adopted a Eulerian bias model in this work. To go beyond the linear order, we expanded the galaxy overdensity as function of the matter density and tidal field up to second order as

$$\begin{array}{c}\hfill {\delta }_{\mathrm{g}}(x)=b_1{\delta }_{\mathrm{m}}(x)+\frac{b_2}{2}({\delta }_{\mathrm{m}}^2(x)-⟨{\delta }_{\mathrm{m}}^2⟩)+\frac{b_s}{2}(s^2(x)-⟨s^2⟩)b_{3\mathrm{n}1}\psi +\cdots ,\end{array}$$(13)

where s is the tidal tensor, introduced in McDonald & Roy (2009), and expressed in Fourier space as

$$\begin{array}{c}\hfill s_{\mathit{ij}}(k)=({\widehat{k}}_i{\widehat{k}}_j-\frac{1}{3}{\delta }_{\mathit{ij}}){\delta }_{\mathrm{m}}(k).\end{array}$$(14)

The number of bias parameters that can be constrained by observations depends on the precision of the data, and the angular scales used in the analysis. In this study, we only vary b₁ and b₂, which is sufficient to obtain a good fit. We fix the values of the remaining bias parameters to

$$\begin{array}{c}\hfill b_{s^2}=-\frac{4}{7}(b_1-1);b_{3\mathrm{nl}}=(b_1-1).\end{array}$$(15)

The higher order terms from Eq. (13) can be used to establish an expansion of P_gg (Pandey et al. 2020). We used FASTPT (McEwen et al. 2016) to computes the various perturbation contributions in the expansion.

4.2. Intrinsic alignment models

Ellipticities are commonly decomposed into three components: the shear, γ_G, the intrinsic ellipticity, γ_I, and a random contribution, γ_rand,

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{tot}}={\gamma }_{\mathrm{G}}+{\gamma }_{\mathrm{I}}+{\gamma }_{\mathrm{rand}}.\end{array}$$(16)

The first two components do not account for the entire shape noise. Chisari & Dvorkin (2013) showed that the predicted dispersion of ellipticities using these terms is two orders of magnitude smaller than the observed value. An additional component, γ_rand, was therefore introduced to quantify unaccounted processes. While not random per se, these processes are not correlated with the large-scale structure; therefore, they form a noise contribution. In the following sections, we present two models to predict γ_I.

4.2.1. NLA

The non-linear alignment (NLA) model (Bridle & King 2007), derived from the linear alignment model (Hirata & Seljak 2004) quantifies the contribution of intrinsic alignment to the ellipticity. It is based on the study of Catelan et al. (2001), which asserts that the mean intrinsic ellipticities of galaxies forming in a tidal field is related to the Newtonian gravitational potential, Ψ_P, via the relation

$$\begin{array}{c}\hfill {\gamma }_{\mathrm{I}}=-\frac{C_1}{4\pi G}({\mathrm{\nabla }}_x^2-{\mathrm{\nabla }}_y^2,2{\mathrm{\nabla }}_x{\mathrm{\nabla }}_y){\mathrm{\Psi }}_{\mathrm{P}},\end{array}$$(17)

where the derivatives along x and y are taken in the transverse plane. The constant C₁ takes the fiducial value of 5 × 10⁻¹⁴ M_⊙ h⁻² Mpc⁻³ based on the SuperCOSMOS survey from Brown et al. (2002). This ties the tidal deformations to cosmological quantities as the Newtonian potential can be related to the linear matter overdensities δ_m via the relation in Fourier space

$$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\mathrm{p}}(k)=-4\pi G\frac{\overline{\rho }}{\overline{D}(z)}{k}^{-2}{\delta }_{\mathrm{m}}(k),\end{array}$$(18)

where $\overline{D}(z)$ is the rescaled growth factor $\overline{D}(z)\propto (1+z)D(z)$ and $\overline{\rho }$ is the mean density of the Universe.

These relations motivated the expressions of the intrinsic-intrinsic and intrinsic-matter power spectra as

$$\begin{array}{c}\hfill P_{\mathrm{II}}(k,z)=F^2(z)P_{\delta }(k,z);P_{\delta I}(k,z)=F(z)P_{\delta }(k,z).\end{array}$$(19)

The function F(z) accounts for the pre-factors in the Poisson equation (18),

$$\begin{array}{c}\hfill F(z)=-A_{\mathrm{IA}}{C}_1{\rho }_{\mathrm{crit}}\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{\overline{D}(z)}·\end{array}$$(20)

Here, ρ_crit is the critical density of the Universe, and Ω_m the matter density parameter, both expressed at z = 0. The A_IA parameter is usually left free to account for a varying amplitude of the intrinsic alignment contribution; in this study, we aim to obtain constraints on this parameter. While initially P_δ(k, z) in Eq. (19) was understood as the linear power spectrum, Bridle & King (2007) instead used the non-linear power spectrum to better account for small-scale processes. While not motivated theoretically, this modification provided a better fit to the data and gave rise to the name non-linear alignment model.

4.2.2. TATT

The tidal alignment and tidal torque (TATT) model was introduced as a second-order correction to the NLA model in Blazek et al. (2019), analogous to the inclusion of higher-order biases in galaxy clustering.

An expansion in the linear density field and collection of contributing terms allow us to write the unprojected intrinsic ellipticity tensor as

$$\begin{array}{c}\hfill {\gamma }_{\mathit{ij}}^{\mathrm{I}}(x)=C_1{s}_{\mathit{ij}}+C_2(s_{\mathit{ik}}{s}_{\mathit{kj}}-\frac{1}{3}{\delta }_{\mathit{ij}}{s}^2)+C_{1\delta }\left({\delta }_{\mathrm{m}}{s}_{\mathit{ij}}\right).\end{array}$$(21)

Summation over repeated indices is implied. The most common parametrisation established in Blazek et al. (2019) defines the first parameter as

$$\begin{array}{c}\hfill C_1(z)=-A_1{\overline{C}}_1{\rho }_{\mathrm{crit}}{\mathrm{\Omega }}_{\mathrm{m}}{D}^{-1}(z),\end{array}$$(22)

which is the parameter already present in the NLA model.

The amplitude of the quadratic contribution of the tidal tensor is parameterised by C₂, given by

$$\begin{array}{c}\hfill C_2(z)=A_{2}5{\overline{C}}_1{\rho }_{\mathrm{crit}}{\mathrm{\Omega }}_{\mathrm{m}}{D}^{-2}(z).\end{array}$$(23)

The term quantifying the correlation between the tidal tensor and the overdensity is

$$\begin{array}{c}\hfill C_{\delta 1}(z)=b_{\mathrm{TA}}{C}_1=-A_{1\delta }{\overline{C}}_1{\rho }_{\mathrm{crit}}{\mathrm{\Omega }}_{\mathrm{m}}{D}^{-1}(z).\end{array}$$(24)

From a theoretical point of view, b_TA can be interpreted as the linear galaxy bias b₁. In the TATT model, however, it is an additional free parameter quantifying the correlations between the tidal tensor and the galaxy density field. We note that we choose A₁, A₂, and A_1δ to be constant for simplicity, but they can also be specified as a function of redshift. Using these parameters one can establish an expression for P_gI and P_II. We refer to Blazek et al. (2019) for details.

4.3. Redshift-space distortions

Peculiar velocity effects, known as redshift space distortions (RSDs), induce anisotropies along the line of sight, which have to be accounted for in our theoretical model. RSDs are added at the power-spectrum level following the Kaiser (1987) parametrisation,

$$\begin{array}{c}\hfill P_{\mathrm{gg}}^s(\mathit{k},z)=b_1^2{[1+\beta (z)\mu ]}^2{P}_{\mathrm{gg}}(k,z),\end{array}$$(25)

with μ as the cosine of the angle between k and the line of sight. The superscript s specifies that the power spectrum is defined in redshift space. The function β(z) = f(z)/b₁(z) is the ratio of the growth rate f(z) and the linear galaxy bias b₁. In the linear regime, only the first three even multipoles are non-vanishing. They are given by Kaiser (1987) as

$$\begin{array}{cc}\hfill P_0^s(k,z)& =(1+\frac{2}{3}\beta (z)+\frac{1}{5}{\beta }^2(z))P(k,z);\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill P_2^s(k,z)& =(\frac{4}{3}\beta (z)+\frac{4}{7}{\beta }^2(z))P(k,z);\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill P_4^s(k,z)& =\frac{8}{35}{\beta }^2(z)P(k,z).\hfill \end{array}$$(28)

The configuration-space multipoles can then be estimated by the Hankel transform

$$\begin{array}{c}\hfill {\xi }_{\ell }(r)=\frac{i^{\ell }}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}{k}^2{\mathrm{j}}_{\ell }(kr)P_{\ell }(k)\mathrm{d}k,\end{array}$$(29)

where j_ℓ represent the spherical Bessel functions of order ℓ. The 2D anisotropic two-point correlation function can then be reconstructed as

$$\begin{array}{c}\hfill {\xi }_{\mathrm{gg}}^s(r_{\mathrm{t}},\mathrm{\Pi })={\displaystyle \sum _{\ell =0,2,4}}{\mathcal{P}}_{\ell }(\mathrm{\Pi }/d){\xi }_{\ell }(d),\end{array}$$(30)

with $d=\sqrt{r_{\mathrm{p}}^2+{\mathrm{\Pi }}^2}$ representing the distance between the galaxies. The estimation of w_gg is then obtained by integrating ξ_gg^s(r_t,Π) along the line of sight. We note that w_g+ is very mildly affected by RSDs and we did not include these effects here.

4.4. Predicting the projected correlation functions

We modelled the projected correlation functions measured in Eqs. (6), (8) and (10) via a Limber projection, and using a Hankel transform:

$$\begin{array}{cc}\hfill w_{\mathrm{gg}}\left(r_{\mathrm{t}}\right)=& {\displaystyle \int \mathrm{d}}z{\mathcal{W}}^{\mathrm{dd}}(z){\displaystyle \int \frac{\mathrm{d}{k}_{\perp }{k}_{\perp }}{2\pi }}{\mathrm{J}}_0\left(k_{\perp }{r}_{\mathrm{t}}\right)P_{\mathrm{gg}}(k_{\perp },z);\hfill \\ \hfill w_{\mathrm{g}+}\left(r_{\mathrm{t}}\right)=& -{\displaystyle \int \mathrm{d}}z{\mathcal{W}}^{\mathrm{ds}}(z){\displaystyle \int \frac{\mathrm{d}{k}_{\perp }{k}_{\perp }}{2\pi }}{\mathrm{J}}_2\left(k_{\perp }{r}_{\mathrm{t}}\right)P_{\mathrm{gI}}(k_{\perp },z);\hfill \\ \hfill w_{++}\left(r_{\mathrm{t}}\right)=& {\displaystyle \int \mathrm{d}}z{\mathcal{W}}^{\mathrm{ss}}(z)\hfill \\ \hfill & \times {\displaystyle \int \frac{\mathrm{d}{k}_{\perp }{k}_{\perp }}{4\pi }}\{[{\mathrm{J}}_0\left(k_{\perp }{r}_{\mathrm{t}}\right)+{\mathrm{J}}_4\left(k_{\perp }{r}_{\mathrm{t}}\right)]P_{\mathrm{II},\mathrm{EE}}(k_{\perp },z)\hfill \\ \hfill & +[{\mathrm{J}}_0\left(k_{\perp }{r}_{\mathrm{t}}\right)-{\mathrm{J}}_4\left(k_{\perp }{r}_{\mathrm{t}}\right)]P_{\mathrm{II},\mathrm{BB}}\}.\hfill \end{array}$$(31)

The function J_ν is the first-kind Bessel function of an order, ν, and 𝒲^ij is the projection kernel (Mandelbaum et al. 2011) given by:

$$\begin{array}{c}\hfill {\mathcal{W}}^{\mathit{ij}}(z)=\frac{n^i(z)n^j(z)}{{\chi }^2(z)\mathrm{d}\chi /\mathrm{d}z}\times {[{\displaystyle \int \mathrm{d}z}\prime \frac{n^i(z\prime )n^j(z\prime )}{{\chi }^2(z\prime )\mathrm{d}\chi /\mathrm{d}z\prime }]}^{-1}.\end{array}$$(32)

χ(z) is the comoving distance at distance z, and nⁱ(z) is the redshift distribution of the shape (i = s) or density (i = d) sample.

4.5. Modelling the lensing contamination

By construction, the estimator w_g+ defined in Eq. (5) is able to pick out intrinsic-alignment correlations at high S/N values. It is, however, also sensitive to a number of other parameters and effects. Some of these effects include a marginalisation over cosmological parameters, baryonic feedback, weak gravitational lensing magnification and shape correlations. In this section, we focus on the latter source of contamination.

To account for the lensing contribution in our measurement, we follow the method of Joachimi et al. (2011), Tonegawa & Okumura (2022). The galaxy-galaxy lensing contribution is given by the angular power spectrum (Joachimi et al. 2011)

$$\begin{array}{cc}\hfill C_{\mathrm{g},\mathrm{G}}(\ell ;z_1,z_2)=b_1{\displaystyle {\int }_0^{{\chi }_{\lim }}}\mathrm{d}\chi \prime & \frac{p_n(\chi \prime |\chi (z_1))q_{\epsilon }(\chi \prime |\chi (z_2))}{{\chi }^{\prime 2}}\hfill \\ \hfill & \times P_{\delta }(\frac{\ell }{\chi \prime },\chi \prime ).\hfill \end{array}$$(33)

This requires knowledge of the galaxy bias. The density tracer p_x(χ, χ₁) and the lensing kernel q_ε(χ′,χ(z₂)) are given in full generality for a photometric survey. With spectroscopic surveys, the conditional density tracer p_x(χ|χ₁) simplifies to the n(χ₁) as spectroscopic redshift estimates can be considered exact at the precision we work at. This implies that the lensing kernel is expressed as:

$$\begin{array}{c}\hfill q_x(\chi ,{\chi }_1)=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\frac{\chi }{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{hor}}}\mathrm{d}\chi \prime p_x(\chi \prime |{\chi }_1)}\frac{\chi \prime -\chi }{\chi \prime },\end{array}$$(34)

It can be further simplified by the insertion of a Dirac delta as p_x(χ′|χ₁):

$$\begin{array}{c}\hfill q(\chi ,{\chi }_1)=\frac{3H_0^2{\mathrm{\Omega }}_{\mathrm{m}}}{2c^2}\frac{\chi }{a(\chi )}\frac{{\chi }_1-\chi }{{\chi }_1}·\end{array}$$(35)

This allows to rewrite Eq. (33) in the simpler form:

$$\begin{array}{c}\hfill C_{\mathrm{g},\mathrm{G}}(\ell ,|{\chi }_1,{\chi }_2)=b_{\mathrm{g}}\frac{3H_0^2{\mathrm{\Omega }}_m}{2c^2}\frac{n({\chi }_1)}{a({\chi }_1)}\frac{{\chi }_2-{\chi }_1}{{\chi }_1{\chi }_2}{P}_{\delta }(\frac{\ell }{{\chi }_1},{\chi }_1).\end{array}$$(36)

Now we do not express our correlation function in terms of χ₁, χ₂, but we instead use the coordinate system, r_t, Π, in physical space. To transform the coordinates, we need to start by defining new coordinates which relate the two systems:

$$\begin{array}{c}\hfill z_m=\frac{1}{2}(z_2+z_1);r_t=\theta \xi (z_m);\mathrm{\Pi }=\frac{c}{H(z_m)}(z_2-z_1).\end{array}$$(37)

As we need the correlation functions for a variety of z_m, Π, r_t we built a grid of ξ(θ, z₁, z₂) in the z₁, z₂ space. We therefore have a grid of correlation functions that we can then use to evaluate the integrated correlation function:

$$\begin{array}{c}\hfill w_{\mathrm{gG}}={\displaystyle {\int }_0^{{\mathrm{\Pi }}_{\max }}}\mathrm{d}\mathrm{\Pi }{\displaystyle \int \mathrm{d}}{z}_{m}W(z_m){\xi }_{\mathrm{gG}}(r_t,\mathrm{\Pi },z_m).\end{array}$$(38)

The kernel here is described in Eq. (32). The weak-lensing contribution has a scale dependence which is analogous to the intrinsic alignment signal on the scales of interest and contributes 1% for CMASS-UNIONS, 3% for LRG-UNIONS, and 6% for ELG-UNIONS. This is consistent with the findings in Sa23 when we make adjustments for the relative signal amplitudes.

Regarding the other sources of contamination to the integrated correlation function listed above are sub-dominant on the scales of interest, and we choose, therefore, not to model them here. The marginalisation over cosmological parameters has been studied and we refer to Johnston et al. (2019) for more details.

4.6. Covariance matrix

4.6.1. Analytical covariance matrix

To estimate the covariance matrix for our different estimators, we use an analytical model, which was first described in Chisari & Dvorkin (2013), and used in Samuroff et al. (2021) and Sa23. To compute the covariance matrix we follow the procedure laid out in Samuroff et al. (2021). We compute the full covariance matrix of the joint vector containing the integrated correlation functions w_g+, w_gg, w₊₊. Each block of the analytical covariance matrix can be decomposed into three types of contributions: a Gaussian part that accounts for cosmic variance, a shot noise part, and a mixed term.

The full general equation for the covariance matrix between two integrated correlation function estimates w_αβ and w_γε for α, β, γ, ε ∈ [g, +] at angular scales r_t, i and r_t, j is given by

$$\begin{array}{cc}& \mathrm{Cov}\left[w_{\alpha \beta }\left(r_{\mathrm{t},i}\right)w_{\gamma \epsilon }\left(r_{\mathrm{t},j}\right)\right]=\frac{1}{\mathcal{A}\left(z_{\mathrm{eff}}\right)}{\displaystyle {\int }_0^{\infty }}\frac{k\mathrm{d}k}{2\pi }{\overline{\mathrm{\Theta }}}_{\alpha \beta }\left(k{r}_{\mathrm{t},i}\right){\overline{\mathrm{\Theta }}}_{\gamma \epsilon }\left(k{r}_{\mathrm{t},j}\right)\hfill \\ \hfill & \times [(P_{\alpha \gamma }(k)+{\delta }_{\alpha \gamma }{N}^{\alpha })\times (P_{\beta \epsilon }(k)+{\delta }_{\beta \epsilon }{N}^{\beta })\hfill \\ \hfill & +(P_{\alpha \epsilon }(k)+{\delta }_{\alpha \epsilon }{N}^{\alpha })\times (P_{\beta \gamma }(k)+{\delta }_{\beta \gamma }{N}^{\beta })].\hfill \end{array}$$(39)

The pre-factor is the inverse of 𝒜(z_eff), the comoving area at redshift z_eff. To obtain this effective redshift, we take the average of the sample redshift distribution n(z), weighted by the square of the number of galaxy pairs N_gal²(z). This is motivated by the fact that the integrated correlation function w_αβ is sensitive to the number of galaxy pairs. When the shape and density samples differ for w_g+, we use the effective area of the smaller sample, as this is where pairs are formed. The effective redshift approximation is accurate for our samples but might need to be refined for more complex redshift distributions. We will improve upon this approximation in future work.

The $\overline{\mathrm{\Theta }}$ functions are bin-averaged Bessel functions,

$$\begin{array}{c}\hfill {\overline{\mathrm{\Theta }}}_{\alpha \beta }(k,r_{\mathrm{t},i})=\frac{1}{A_i}{\displaystyle {\int }_{r_{\mathrm{t},i}^{\min }}^{r_{\mathrm{t},i}^{\max }}}2\pi r{\mathrm{\Theta }}_{\alpha \beta }(k,r)\mathrm{d}r.\end{array}$$(40)

with A_i the area of bin i, Θ_gg = J₀, Θ_g+ = J₂ and Θ₊₊ = J₀ + J₄.

The two different shot-noise terms are $N^g=1/{\overline{n}}_{\mathrm{density}}(z_{\mathrm{eff}})$ corresponding to Poisson noise from discrete galaxy positions, and $N^{+}={\sigma }_{\epsilon }^2/{\overline{n}}_{\mathrm{shape}}(z_{\mathrm{eff}})$ from intrinsic galaxy shapes³. We use σ_ε as the single-component shape noise term defined in Gatti et al. (2021) as

$$\begin{array}{c}\hfill {\sigma }_e^2=\frac{1}{2}[\frac{{\displaystyle {\sum }_i}{\left(w_i{e}_{i,1}\right)}^2}{{\left({\displaystyle {\sum }_i}{w}_i\right)}^2}+\frac{{\displaystyle {\sum }_i}{\left(w_i{e}_{i,2}\right)}^2}{{\left({\displaystyle {\sum }_i}{w}_i\right)}^2}]\left[\frac{{\left({\displaystyle {\sum }_i}{w}_i\right)}^2}{{\displaystyle {\sum }_i}{w}_i^2}\right].\end{array}$$(41)

The estimator with the highest S/N carrying the bulk of information on intrinsic-alignment parameters is the shape-density integrated correlation function w_g+ defined in Eq. (31). We explicitly write out the covariance Eq. (39) for the special case of (αβ) = (γε) = (g + ) as

$$\begin{array}{cc}\hfill \mathrm{Cov}& \left[w_{\mathrm{g}+}\left(r_{\mathrm{t},i}\right)w_{\mathrm{g}+}\left(r_{\mathrm{t},j}\right)\right]=\frac{1}{\mathcal{A}\left(z_{\mathrm{eff}}\right)}{\displaystyle {\int }_0^{\infty }}\frac{k\mathrm{d}k}{2\pi }{\overline{\mathrm{\Theta }}}_{\alpha \beta }\left(k{r}_{\mathrm{t},i}\right){\overline{\mathrm{\Theta }}}_{\gamma \epsilon }\left(k{r}_{\mathrm{t},j}\right)\hfill \\ \hfill & \times [(P_{\mathrm{gg}}(k)+\frac{1}{{\overline{n}}_{\mathrm{density}}(z_{\mathrm{eff}})})\times (P_{\mathrm{II}}(k)+\frac{{\sigma }_{\epsilon }^2}{{\overline{n}}_{\mathrm{shape}}(z_{\mathrm{eff}})})+P_{\mathrm{gI}}^2(k)].\hfill \end{array}$$(42)

The mixed power spectrum P_gI²(k) does not have a shot-noise contribution as it corresponds to a cross-correlation, resulting in the corresponding Kroenecker deltas in Eq. (39) to vanish.

We plot the diagonal of this covariance in Fig. 4. On small scales, the covariance matrix is shot-noise dominated. On larger scales, the mixed term of the clustering power spectrum and galaxy shape noise becomes the leading term. We note that this mixed term is missing from the covariance matrix of Sa23 (their Eqs. 26 and 27). Since this term is dominant for r_t ⪆ 20 Mpc, this omission resulted in a bias of their results, and an underestimation of the error bars published in Sa23. We confirmed this finding by applying our analytical covariance matrix without this mixed term to DES Y3 data, leading to the same biases as in Sa23. We note that this term was still present in the earlier publication Samuroff et al. (2021). In Sect. 5, we present a re-analysis of the DES Y3 and BOSS correlations and show the corrected intrinsic-alignment constraints.

Fig. 4. Square root of the diagonal of the analytical covariance matrix of wg+ for the CMASS-UNIONS sample (green solid line). The individual contributions are cosmic variance (dashed; the PggPII term is identical to the PgI contribution); mixed clustering – ellipticity shot noise (dotted); mixed intrinsic – density shot noise (long dash-dotted); and pure shot noise (short dash-dotted). See Eq. (42) and text for details. For comparison, we plot the jackknife variance (purple solid line), which is in excellent agreement with the analytical covariance over the relevant scales.

4.6.2. Jackknife covariance matrix

To validate our analytical covariance matrix, we compared it to a data resampling covariance using the jackknife method (Si15, Johnston et al. 2019; Fortuna et al. 2021b; Joachimi et al. 2011). We split the observed area into n_s patches using the K-means algorithm implemented in Treecorr. A number of n_s realisations were built by removing one patch from the data set at a time. The correlation function was computed on the remaining data. The idea is that the patches contain the largest scales probed by the correlation function. By removing a patch, we can estimate the variation in the data without any external analytical formulation or simulation. The jackknife estimation of the covariance matrix is given by:

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{ij}}=\frac{n_{\mathrm{s}}-1}{n_{\mathrm{s}}}{\displaystyle \sum _{k=1}^{n_{\mathrm{s}}}}(x_i^k-\overline{x_i})(x_j^k-\overline{x_j}),\end{array}$$(43)

with $x_i^k$ as the k-th correlation function computed from the n_s jackknife samples. The ij indices are the bins in r_t between which we want to compute the covariance. The average $\overline{x}$ is taken over the n_s jackknife estimations.

Figure 4 shows the excellent agreement between the analytical and jackknife covariance over the angular scales of interest. Since the jackknife covariance is noisy, in particular on the off-diagonal elements (Samuroff et al. 2021), we used the analytical covariance for the parameter inference.

4.7. Scale cuts

As described in Sect. 2.3, the eBOSS samples were originally selected for galaxy clustering on large scales, with the explicit goal of measuring the baryonic acoustic oscillation scale around 100 h⁻¹ Mpc. The sample has therefore not been optimised for small-scale correlations, where the intrinsic alignment becomes most relevant. To include smaller scales would require a more densely sampled survey such as GAMA or DESI.

Furthermore, we did not include small scales in our analysis as baryonic feedback and non-linearities can become dominant on these scales. To be consistent with previous works, we introduced a small-scale cut at 9 Mpc (Joachimi et al. 2011, Si15, Sa23). We also excluded points above 100 Mpc from the fit. This was motivated by the projection effects described in Singh & Mandelbaum (2016): The integration over Π along the line of sight causes a suppression for large Π, as galaxy shapes pointing towards the density tracer at large line-of-sight separation will appear smaller along the radial axis when projected on the 2D observer plane. While this projection effect exists at all transverse separations r_t, it is sub-dominant where the amplitude of w_g+ is large. The sensitivity of w_g+ to P_gI in the [r_T, Π] plane is quantified in Joachimi et al. (2011). To mitigate these projections effect, the 3D Yamamoto estimator (Kurita et al. 2021) or the real-space 3D correlator (Singh et al. 2024) can be used, which circumvent the issue by using the full 3D information.

For consistency with previous works (Sa23), we take 21 log-spaced bins in the range r_t ∈ [0.1, 350] Mpc to compute our projected correlation functions. For the cosmological parameters we assumed h = 0.69, n_S = 0.97, A_S = 2.15 × 10⁻⁹, and Ω_m = 0.3. We used the emcee package to fit our models. All correlation functions are computed with the TreeCorr package (Jarvis 2015).

5. Results

5.1. Integrated correlation functions

We show in Fig. 5 the measurement of w_g+ and w₊₊ for the different UNIONS-BOSS/eBOSS samples as described in Sect. 2. We observe strong intrinsic alignment in our two red galaxy samples, CMASS and LRG, consistent with what has been shown in the literature. In grey, we shade the areas excluded from the fit, as motivated in Sect. 4.7. For CMASS-UNIONS, our measured w_g+ is well fitted by the NLA model down to scales of r_t ≈ 5 Mpc. Below this scale, the CMASS galaxies are sparsely sampled, and we therefore do not see the appearance of a one-halo term. Overall we remain on linear scales where the NLA model performs well. The w₊₊ plot shows an apparent positive correlation, which seems under-fitted by the NLA model. The reduced χ² summarised in Table 2 indicates that the fit is good, but the null hypothesis is also not excluded for this data vector. Our fit of $A_{\mathrm{IA}}=4.{02}_{-0.31}^{+0.31}$ represents one of the most significant detections in the literature.

Fig. 5. Measurements of wg+ (upper panels) and w++ (lower row) for the UNIONS-CMASS (left), UNIONS-LRG (middle), and UNIONS-ELG (right column) samples. The solid lines show the best-fit NLA model, where the shaded areas indicate the 1σ uncertainty associated with the fit. The grey areas are excluded from the fit as explained in Sect. 4.7.

In the LRG sample, the picture is similar, although the measurement is slightly more scattered. The fit follows the data points well down to scales of r_T ≈ 8 Mpc. The large error bars are due to the scarce density of the LRG galaxies, which in turn are responsible for the large uncertainty on the A_IA fit, even though the 80 000 galaxies form a consequential sample for this type of measurement. The w₊₊ correlation function is noisy and has negative points. Its error bars are shot-noise dominated and we do not give any strong interpretation of w₊₊ as the measurement is compatible both with the NLA model and with the null hypothesis. The total fit gives $A_{\mathrm{IA}}=3.3_{-1.0}^{+1.0}$, a value we compare to other measurements in Sect. 5.4.

For the ELG sample, we expected no intrinsic alignment signal as this is a sample mainly composed of blue galaxies. These have repeatedly shown to exhibit intrinsic alignment amplitudes consistent with zero. While we have measured a relatively high value of A_IA, the error bars are still compatible with A_IA = 0 at the 1.1σ level. The w₊₊ vector also appears to exhibit a strong correlation on larger scales, but here again, the points are correlated and the data vector is compatible with the null hypothesis.

All results are summarised in Table 2. The global reduced χ² indicates that our overall fit is good. We note that while we report the individual reduced χ² for each data vector separately, they are for reference only as the fit is performed jointly on all 3 data vectors, and the covariance between them is not negligible. We considered 5 degrees of freedom for w₊₊ and w_g+ (6−1), 6 for w_gg (8−2) and 17 for the total fit (20−3).

5.2. Model comparison

With our high S/N measurement, we test the TATT model (Blazek et al. 2019) to see if the additional degrees of freedom can help in modelling the w_g+ signal. The TATT model is a higher-order expansion, as explained in Sect. 4.2.2. We therefore fit w_g+ and w₊₊ down to 4 Mpc. The constraints on the second-order TATT parameters from CMASS-UNIONS can be seen in Fig. 7. From Fig. 6, the NLA and TATT best-fit models are in very good agreement, within 1σ down to small scales. It appears that the scales outside the fitted range are better captured by the NLA model. Considering the reduced χ² we get,

$$\begin{array}{c}\hfill \mathrm{NLA}:{\chi }^2=1.40;\mathrm{TATT}:{\chi }^2=1.32.\end{array}$$(44)

Fig. 6. Best-fit lines and 1σ contours for the NLA and TATT, respectively. The TATT model is fit down to 4 Mpc but both models appear very consistent. The 1σ TATT contours are obtained by drawing 500 models from the posterior and calculating the 68% confidence interval at eat rt.

Fig. 7. Posterior contours of the NLA and TATT model fitted to the CMASS-UNIONS sample. We see no deviation from the NLA case as A2 and bTA appear consistent with 0.

These reduced χ² have been obtained for the full data vectors meaning, 8 + 6 + 6 − 3 degrees of freedom for NLA and 8 + 8 + 8 − 5 degrees of freedom for TATT. To compare them, we can use the Bayesian information criterion (BIC) defined as

$$\begin{array}{c}\hfill \mathrm{BIC}=klog(N_{\mathrm{p}})+{\chi }^2,\end{array}$$(45)

where k is the number of parameters and N_p is the number of fitted points. The BIC is a commonly employed metric which balances goodness of fit and model complexity by penalizing an excessive number of parameters, which would indicate overfitting. We respectively obtain:

$$\begin{array}{c}\hfill \mathrm{NLA}:\mathrm{BIC}=32.8;\mathrm{TATT}:\mathrm{BIC}=40.7.\end{array}$$(46)

From Δ_BIC(TATT − NLA) = 7.9 we conclude that the additional TATT degrees of freedom are not needed here. We can nevertheless make a few observations, as follows. First, the A₁ and A₂ parameters are strongly degenerate. This is expected as A₁ is the proportionality factor for s_ij, and A₂ for s_ij², displaying similar behaviour for small values of the tidal tensor. To decorrelate the two parameters, the largest possible variations of the tidal field must be probed so that the linearisation breaks down and the quadratic regime can be distinguished from the linear regime. It is not obvious if this will happen on scales relevant for weak lensing, since we are in the perturbation-theory domain where variations are small.

Furthermore, we observe a degeneracy between b_TA and A₂. As b_TA is multiplied with the overdensity δ_ms_ij, a similar argumentation holds for these two parameters as for the (A₁, A₂) pair. From theoretical arguments b_TA can be expected to be comparable to b₁ as it correlates the tidal tensor with the overdensity field, traced here by galaxies. The fact that we constrain b_TA to values compatible with zero, far away from b₁ ≈ 2, can be an indication that the linear assumptions allowing us to take the product of the linearised density field with the tidal tensor can not be naively applied here.

Finally, we see no correlation between b_TA and A₁. This can be interpreted by the model picking up on the density field and decorrelating it from the tidal tensor. As both parameters are multiplied with each other to form the parameter A_d, the fact that raising one is not simply equivalent to lowering the other indicates that the measurement effectively decorrelates matter variations from tidal-field modifications. Upcoming measurements will allow us to go down to smaller scales and properly understand the behaviour of these parameters for different samples.

5.3. Luminosity dependence

We measured the luminosity scaling of our intrinsic-alignment measurements of red galaxies and compared the results to previous studies. As described in Sect. 3.4, multiple observations have shown a correlation between luminosity and intrinsic alignment. Such a correlation is motivated by recent findings in Fortuna et al. (2025), who showed that luminosity scales with halo mass. A similar correlation was already hinted at in Xu et al. (2023b). This explains why luminosity is a good predictor of intrinsic alignment, as halo masses determine the amplitude of the tidal field. In Fig. 8 we plot our UNIONS measurements for the CMASS sample split into three equi-populated bins with limits log(L/L₀)∈[−0.8; −0.14; 0.12; 0.6]. The intrinsic-alignment amplitude shows a clear increase with luminosity. Results are shown in Table 3.

Fig. 8. Integrated correlation function, wg+, as measured on the CMASS-UNIONS sample split in 3 bins of luminosity, showing increasing luminosity from left to right panels. One can see a clear increase in the amplitude of the signal.

Fig. 9 shows the NLA amplitude A_IA from our UNIONS and DES measurements with the CMASS and LRG samples, together with a collection of data points from the literature. Our samples are in agreement with measurements from Joachimi et al. (2011), Si15, Johnston et al. (2019), Fortuna et al. (2021b). The CMASS-UNIONS points are slightly higher than the CMASS-DES points, but both seem to scale in the same way with luminosity. On the contrary, the LRG-UNIONS point is lower than the LRG-DES point. We will come back to this discrepancy in Sect. 5.4.

Fig. 9. Intrinsic alignment amplitude as a function of luminosity as measured by different surveys. The fitted function is α⟨L/L0⟩β. We can see very similar scaling with luminosity when comparing our CMASS-DES and CMASS-UNIONS measurements.

To quantify the luminosity-dependence and evolution of the intrinsic-alignment amplitude, we fitted all data points from Fig. 9, assuming their independence. Various functional forms have been tested; Joachimi et al. (2011) introduced a scaling relation of

$$\begin{array}{c}\hfill P_{\mathrm{gI}}(k,z,L)=A_{\mathrm{IA}}{b}_1{P}_{\delta \mathrm{I}}(k,z){\left(\frac{1+z}{1+z_0}\right)}^{\eta }{\left(\frac{L}{L_0}\right)}^{\beta }.\end{array}$$(47)

Motivated by recent works which disfavour any strong z-dependency (Si15; Johnston et al. 2019; Sa23), we fixed η = 0. In addition to this single-power-law model, a broken power law was introduced in Fortuna et al. (2021b) to better account for the flattening of the A_IA − log L relation at low luminosities.

We test the single and double power law models, defined as:

$$\begin{array}{cc}\hfill A_{\mathrm{IA}}(L)& =A_{\beta }{\langle \frac{L}{L_0}\rangle }^{\beta };\hfill \\ \hfill A_{\mathrm{IA}}(L)& =A_{\beta }{\langle \frac{L}{L_{\mathrm{break}}}\rangle }^{\beta }\{\begin{array}{cc}\beta ={\beta }_1\hfill & \mathrm{if}L<L_{\mathrm{break}}\hfill \\ \beta ={\beta }_2\hfill & \mathrm{if}L>L_{\mathrm{break}}.\hfill \end{array}\hfill \end{array}$$(48)

Figure 9 demonstrates that both models provide decent fits to the data. The reduced χ² for the single and double power-law are χ² = 2.42 with 28 degrees of freedom (30−2) and χ² = 2.19 with 26 degrees of freedom (30−4), respectively. These high values of reduced χ² indicate that both models perform rather poorly which is expected given the scatter in Fig. 9. We can compute the BIC and thus obtain:

$$\begin{array}{c}\hfill \mathrm{Single}\mathrm{PL}:\mathrm{BIC}=74.6;\mathrm{Double}\mathrm{PL}:\mathrm{BIC}=70.6.\end{array}$$(49)

The difference Δ_BIC(SPL − DPL) = 4.0 indicates that the double power law is slightly favoured, a result consistent with the recent literature.

While converging on the A_IA(L) relationship is important, we want to highlight the complexity that Fig. 9 is omitting by using a single predictive parameter. While the redshift dependence has been found to be sufficiently well captured by the naïve NLA model, many observational properties strongly affect the intrinsic alignment amplitude, as explained in the introduction. One can therefore not apply this intrinsic alignment-luminosity relationship to derive priors for a different survey without taking the risk of biasing the measurements. Including additional effects which explain the scatter in Fig. 9 will be the subject of a future work. Figure 3 shows that these spectroscopically selected galaxies have apparent magnitudes which fall inside our UNIONS source sample, although they are within the bright tail of the distribution. To establish a reliable prior, it is essential to map the intrinsic alignment amplitudes across redshift, luminosity, and colour space, for the entire weak-lensing source sample, with sufficient accuracy. Achieving better coverage of this three-dimensional parameter space requires either complete spectroscopic samples which include fainter objects, or deep narrow-band fields that provide reliable photometric redshifts for all galaxies. This study examines only a small subset of this parameter space. However, the agreement between different surveys, as discussed in Sect. 5.4, is a necessary step for further progress.

5.4. Comparing UNIONS to other lensing samples

One of the main motivations of this work is a comparison of the measured intrinsic alignment quantities between UNIONS and DES Y3. The interest lies in the understanding of the universal nature of intrinsic alignment. As stated in the introduction, observational effects can influence this kind of measurement. Important effects discussed in the literature include the galaxy shape measurement method or the colour of the sampled galaxies. Here, we have a unique opportunity to study a homogeneous set of galaxies selected by the SDSS imaging survey, with two distinct shape measurement catalogues. While a number of similarities can be noted in the ShapePipe and DES Y3 pipelines, a few important differences between the two datasets exist, as explained in Sect. 2. A good agreement between the UNIONS and DES Y3 constraints would, therefore, indicate a certain control on the effects affecting intrinsic alignment measurements. This is crucial if one wants to use existing measurements as a prior for future surveys such as Rubin/LSST and Euclid. Tight priors improve the statistical power to constrain cosmological parameters but will induce biases if inaccurate.

To test the compatibility of different shape measurement samples, we can compare our measurements to a third dataset for the CMASS-type galaxies (Kurita & Takada 2023). This dataset has the advantage of covering the entire CMASS sample and, therefore, overlaps with both UNIONS and DES. That study was performed with a newly developed, more accurate, 3D correlation function estimator (Kurita et al. 2021), which decreases the error bars on the estimated parameters. Their measurement relies on the shape catalogue described in Reyes et al. (2012), Mandelbaum et al. (2013), which is based on the ‘re-Gaussianisation’ method described in Hirata & Seljak (2003). This method computes moments of the shape and corrects for a non-Gaussian PSF. We note that we make use of this shape catalogue for LOWZ shapes with which we validate our pipeline in Appendix A.

In Fig. 10, we show the overlap between the different posteriors. For the CMASS density tracer, the UNIONS and DES Y3 posteriors agree at 1.20σ. The samples do not differ in image selections but they differ slightly in redshift. The DES Y3 sample was chosen in the range z ∈ [0.4, 0.6], the UNIONS sample covers the range z ∈ [0.4, 0.7], and the Kurita sample is in the range z ∈ [0.5, 0.7]. From previous studies on redshift dependence, we do not anticipate these discrepancies to result in significant variations. Although Fig. 10 might suggest a coherent north-south difference in intrinsic alignment, we do not consider this significant, given the inversion of the UNIONS and DES amplitudes in the LRG plot. Nonetheless, we report the differing mean galactic extinctions in Table 1, as this is the only spatially correlated property we could identify that might influence shape measurements.

Fig. 10. 1D marginal posterior of the NLA intrinsic-alignment amplitude AIA of the same SDSS density tracer. Left panel: UNIONS, DES Y3 (both our measurements), and the results from Kurita & Takada (2023) in the North (NGC) and South Galactic Cap (SGC) using CMASS as density tracer. Right panel: Our measured UNIONS and DES Y3 results with the LRG sample.

For the LRG sample, we find an agreement between the DES and UNIONS measurements at the 1.51σ level. We conclude that as long as the error bars are properly taken into account the measurements are in sufficiently good agreement to form a prior for stage IV surveys.

5.5. Systematic diagnostics

In Fig. 11, we show the cross-modes, w_g×, computed from our samples. These are expected to be compatible with zero in the absence of systematic contributions (Joachimi et al. 2011), analogously to lensing B-modes in cosmic shear. In Table 4 we show the deviation from zero of the cross-mode w_g×, and compare it to the one of the intrinsic-alignment mode w_g+. We convert the reduced χ² into a Gaussian σ estimate for easier interpretation. The cross-mode w_g× is compatible with zero at 1σ for all samples, which strengthens the confidence in our results. In contrast, the deviations from zero for w_g+ testify of our high S/N measurement. For the ELG sample, both w_g+ and w_g× have remained consistent with zero at the 1σ level.

Fig. 11. Measurements of wg× for the CMASS-UNIONS, LRG-UNIONS and ELG-UNIONS samples. Deviations from 0 are characterised in Table 4 and are all compatible at the 1σ level.

5.6. Quantifying PSF errors

The effects of PSF mispecifications on cosmic shear have been studied in detail (Heymans et al. 2012; Jarvis et al. 2021; Giblin et al. 2021). Even small PSF residuals can produce biases in cosmological analyses and mimic the correlations induced by weak gravitational lensing. Quantities such as PSF leakage (Jarvis et al. 2016), ρ-statistics (Rowe 2010), and τ-statistics (Gatti et al. 2021; Guerrini et al. 2024) have been introduced to measure and marginalise over the biases. For direct measurements of intrinsic alignments, the impact of PSF errors has been less studied, as their effect is comparatively weaker. Singh & Mandelbaum (2016) computed the additive bias from PSF mis-modelling using the quantity

$$\begin{array}{c}\hfill A_{\mathrm{PSF}}(r_{\mathrm{t}})=\frac{1}{r_{\mathrm{t},\max }-r_{\mathrm{t},\min }}{\displaystyle \underset{r_{\mathrm{t},\min }}{\overset{r_{\mathrm{t},\max }}{\int }}}\mathrm{d}{r}_{\mathrm{t}}\frac{⟨{\epsilon }_{\mathrm{obs}}{\epsilon }_{\mathrm{PSF}}⟩(r_{\mathrm{t}})}{⟨{\epsilon }_{\mathrm{PSF}}{\epsilon }_{\mathrm{PSF}}⟩(r_{\mathrm{t}})}·\end{array}$$(50)

This allows us to quantify a systematic contribution to the shape-shape correlation function ξ₊₊ given by

$$\begin{array}{c}\hfill {\xi }_{++}^{\mathrm{sys}}\approx A_{\mathrm{PSF}}^2\times {\xi }_{++}.\end{array}$$(51)

A shape measurement method that does not correct for the PSF, such as the isophotal method used in Singh & Mandelbaum (2016), can display a large amplitude of A_PSF.

Given the progress in imaging and shape measurement in recent years and the increasingly stringent requirements on additive biases, we calculate this quantity. As shown in Fig. 12, we measure values of below 2% for our red galaxies, which is sufficiently low to be ignored within our analysis. Interestingly, it varies between samples, which is the effect of a varying leakage across size and S/N, see also Li et al. (2023). In contrast to red galaxies, we find a 4% additive bias for the ELG sample; such a bias can be calibrated via subtraction as in Li et al. (2023). We will address this in future work, as the current error bars on the ELG-UNIONS sample are large due to the limited number of galaxies available for this study.

Fig. 12. Amplitude expressed by Eq. (50) of the additive bias contamination of the shape-shape correlation function. From left to right, we show our measurements using the Si15 data corresponding to the isophotal, de Vaucouleurs, and re-Gaussianisation shape catalogues, followed by UNIONS-CMASS, UNIONS-LRG, and UNIONS-ELG.

The measurement of A_PSF informs us about additive leakage to the shape-shape correlation. Multiplicative and additive biases due to PSF contamination for shape-density correlations have been recently quantified in Zhang et al. (2024) (hereafter Z24). Additive biases are not expected to contribute significantly, since the UNIONS PSF is not expected to show a dependence on the positions of SDSS galaxies. The so-called λ-statistics developed in Z24 allow us to quantify the presence or absence of such biases. A multiplicative bias induced by the PSF can act as an additional calibration factor. In Z24 λ-statistics were introduced for galaxy-galaxy lensing. Here, we adapt their work for the intrinsic-alignment correlation function, w_g+.

We started with the PSF error propagation model from Paulin-Henriksson et al. (2008),

$$\begin{array}{c}\hfill {\epsilon }^{\mathrm{obs}}={\epsilon }^{\mathrm{int}}+(1+m)\gamma +c+\delta \epsilon +\alpha {\epsilon }^{\mathrm{psf}}.\end{array}$$(52)

In this model, the multiplicative bias term m and additive bias term c do not capture PSF information, which differs from the standard generic parametrisation where m and c represent jointly all systematic contributions. The PSF leakage term αε^psf captures a residual shape contribution from leakage of the PSF, typically attributed to an inaccurate PSF correction during the shape measurement process. The δε term, on the other hand, captures residuals in the PSF measurement at the galaxy position and is due to errors in PSF measurement, modelling, and interpolation and can be understood as a PSF additive bias.

To estimate the PSF error at the galaxy position Paulin-Henriksson et al. (2008) derived a formula using Gaussian error propagation of:

$$\begin{array}{c}\hfill \delta \epsilon =({\epsilon }^{\mathrm{obs}}-{\epsilon }^{\mathrm{psf}})\frac{\delta T^{\mathrm{psf}}}{T}-\frac{T^{\mathrm{psf}}}{T}\delta {\epsilon }^{\mathrm{psf}}.\end{array}$$(53)

The residuals on the ellipticity, δε^psf, and size, δT^psf, were obtained from a set of reserve stars, which are not used in the PSF modelling process. As we did not have access to this information at each position, we treated these residuals statistically as a constant pre-factor.

With w_g+ we can correlate the intrinsic ellipticity with the position of a close-by galaxy (i.e. $d<\sqrt{{\mathrm{\Pi }}_{\max }^2+r_t^2}$) via the estimator

$$\begin{array}{c}\hfill {\xi }_{+}=⟨{\epsilon }_{+}n⟩.\end{array}$$(54)

We use ξ₊^obs to denote the observed correlation function and ξ₊^theo as the underlying correlation of the intrinsic shape we are trying to quantify. We can define the change in correlation induced by a PSF mis-specification as

$$\begin{array}{c}\hfill \delta {\xi }_{+}={\xi }_{+}^{\mathrm{obs}}-{\xi }^{\mathrm{theo}}.\end{array}$$(55)

By plugging Eq. (52) into Eq. (54) we can write:

$$\begin{array}{cc}\hfill \delta {\xi }_{+}=\langle {\epsilon }_{\mathrm{t}}^{\mathrm{s}}n\rangle & +\langle \frac{T^{\mathrm{psf}}}{T}\frac{\delta T^{\mathrm{psf}}}{T^{\mathrm{psf}}}{\epsilon }_{\mathrm{t}}^{\mathrm{s}}n\rangle -\langle \frac{T^{\mathrm{psf}}}{T}\frac{\delta T^{\mathrm{psf}}}{T^{\mathrm{psf}}}{\epsilon }_{\mathrm{t}}^{\mathrm{psf}}n\rangle \hfill \\ \hfill & -\langle \frac{T^{\mathrm{psf}}}{T}\delta {\epsilon }_{\mathrm{t}}^{\mathrm{psf}}n\rangle +\langle \alpha {\epsilon }_{\mathrm{t}}^{\mathrm{psf}}n\rangle .\hfill \end{array}$$(56)

The $\langle {\epsilon }_{\mathrm{t}}^{\mathrm{s}}n\rangle$ contribution vanishes as intrinsic ellipticities do not correlate with random foreground positions. This allows us to define the λ-statistics,

$$\begin{array}{cc}\hfill {\lambda }_1& =\langle {\epsilon }_{\mathrm{t}}^{\mathrm{psf}}n\rangle ;\hfill \\ \hfill {\lambda }_2& =\langle \frac{\delta T^{\mathrm{psf}}}{T^{\mathrm{psf}}}{\epsilon }_{\mathrm{t}}^{\mathrm{psf}}n\rangle ;\hfill \\ \hfill {\lambda }_3& =\langle \delta {\epsilon }_{\mathrm{t}}^{\mathrm{psf}}n\rangle .\hfill \end{array}$$(57)

Since we do not expect the UNIONS PSF to be correlated with SDSS positions as stated above, the density tracer n can here be considered as random positions. The λ-statistics then indicate additive errors due to random PSF errors, and potential PSF correlations with the mask or the survey footprint.

The computation of λ₁, λ₂, λ₃, and α are detailed in Z24. With that, we rewrite the PSF contamination with respect to ξ₊ as

$$\begin{array}{c}\hfill \delta {\xi }_{+}=\langle \frac{T^{\mathrm{psf}}}{T}\rangle \langle \frac{\delta T^{\mathrm{psf}}}{T^{\mathrm{psf}}}\rangle {\xi }_{+}^{\mathrm{theo}}+\alpha {\lambda }_1-\langle \frac{T^{\mathrm{psf}}}{T}\rangle ({\lambda }_2+{\lambda }_3).\end{array}$$(58)

We can remove the dependence on the theoretical prediction ξ₊^theo by combining Eqs. (55) and (58) to write the bias of the measured quantity,

$$\begin{array}{c}\hfill \delta {\xi }_{+}={\xi }_{+}^{\mathrm{obs}}-\frac{{\xi }_{+}^{\mathrm{obs}}-\alpha {\lambda }_1+\langle \frac{T^{\mathrm{psf}}}{T}\rangle ({\lambda }_2+{\lambda }_3)}{1+\langle \frac{T^{\mathrm{psf}}}{T}\rangle \langle \frac{\delta T^{\mathrm{psf}}}{T^{\mathrm{psf}}}\rangle },\end{array}$$(59)

To obtain the contamination with respect to w_g+, we considered ∫_Πmin^Π_maxδξ₊dΠ. We plot δw_g+/w_g+ for the CMASS, LRG, and ELG samples in Fig. 13. The PSF contamination remains under 2% on the scales of interest. This indicates that PSF errors are subdominant compared to our statistical power. In the future, when direct intrinsic alignment measurements reach precisions of a few percent, a significant measurement of PSF-tracer correlations can be used to account and calibrate PSF biases.

Fig. 13. Relative contribution to the wg+ signal of PSF induced correlations. We see that on the scales of interest the contribution is below 2% which we judge satisfactory for current precision. The error bars blow up when wg+ is close to 0.

6. Conclusion

In this paper, we made use of the UNIONS shape catalogue and the broad overlap with BOSS/eBOSS spectroscopy to carry out a direct measurement of the intrinsic alignment of galaxies. We obtain very significant measurements using 201 639 CMASS, 78 134 LRG, and 14 762 ELG galaxies. Our results are consistent with previous studies showing a strong intrinsic alignment for red galaxies and an alignment compatible with zero for ELGs. These results demonstrate that the UNIONS observations, thanks to their excellent image quality, form a competitive dataset with which we are able to measure the intrinsic alignment in the northern sky. In particular, we highlight the following points:

• We obtained detections of intrinsic alignment at 13σ with CMASS, at 3σ with LRG and at 1σ with ELG. The latter result is compatible with the null hypothesis. The NLA model was able to produce a compelling fit for these measurements in the range of r_t ∈ [9, 100] Mpc.

• We tested whether the highest S/N measurement (CMASS-UNIONS) would allow for a detection of the higher-order terms of the TATT. To that end, we fit transverse scales down to 4 Mpc. We found no preference for non-zero values of these parameters. However, in the data, the tidal field contribution, A₁, de-correlated from the density × tidal field contribution, b_TA. The BIC result seems to indicate a preference for the linear NLA model over TATT.

• Our measured NLA intrinsic-alignment amplitude scales with the luminosity, as also observed in the literature. A double power law seems to provide a reasonable fit, although the reduced χ² is large due to the high scatter.

• By adapting a novel framework from Z24, we estimated the bias contributions of the PSF leakage to the intrinsic alignment measurement, which is below 2%.

• The cross-component w_g× was measured to be consistent with zero at 1σ for all samples, indicating an absence of obvious systematic contributions.

• We re-analysed the DES Y3 data, and find intrinsic-alignment amplitude in agreement with UNIONS, at 1.2σ for CMASS and 1.5σ for LRGs. We found further agreement with other measurements obtained in Kurita & Takada (2023) for CMASS galaxies. We conclude that all of these measurements are compatible with each other and, therefore, they can serve as the basis for an informative prior for stage-IV surveys, as long as the uncertainties are properly accounted for.

This work demonstrates that the UNIONS sample is mature and ready for competitive science cases. In the near future, we aim to use the upcoming spectroscopic data in the northern hemisphere to probe the intrinsic alignment on smaller scales and across a wider ranger of colours and luminosities. The goal is to improve our understanding at the dawn of stage-IV weak lensing surveys.

Acknowledgments

We thank Benjamin Joachimi and Elisa Chisari for informing discussions. FHP acknowledges support from the Centre National d’Études Spatiales (CNES). We are honored and grateful for the opportunity of observing the Universe from Maunakea and Haleakala, which both have cultural, historical and natural significance in Hawaii. This work is based on data obtained as part of the Canada-France Imaging Survey, a CFHT large program of the National Research Council of Canada and the French Centre National de la Recherche Scientifique. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA Saclay, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l’Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. This research is based in part on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. Pan-STARRS is a project of the Institute for Astronomy of the University of Hawaii, and is supported by the NASA SSO Near Earth Observation Program under grants 80NSSC18K0971, NNX14AM74G, NNX12AR65G, NNX13AQ47G, NNX08AR22G, 80NSSC21K1572 and by the State of Hawaii. LB is supported by the PRIN 2022 project EMC2 – Euclid Mission Cluster Cosmology: unlock the full cosmological utility of the Euclid photometric cluster catalogue (code no. J53D23001620006). H. Hildebrandt is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative Research Center SFB1491, an ERC Consolidator Grant (No. 770935), and the DLR project 50QE2305.

References

Appendix A: Validating with LOWZ data from Singh 2015

To validate our pipeline, our initial attempt was to use the Sa23 DES measurements. During the course of this validation, we realised that their w_g+ covariance matrix was incomplete as explained in Sect. 4.6. To ensure good agreement with previous data sets we therefore validated our w_g+ measurement and covariance matrix estimate with data from Si15. Due to some uncertainty in the way the extinction mask is applied, we only managed to reach a 2% agreement on the number of galaxies in the shape catalogue. Nevertheless, the agreement seen in Fig. A.1 is visually very good. When fitting the NLA model we recover A₁ = 4.62 ± 0.44, in excellent agreement with the A₁ = 4.6 ± 0.5 stated in the Si15 paper.

Fig. A.1. Reproduction of the LOWZ measurement from Si15 for validation purposes. While we recovered a 2% agreement on the number of galaxies in the shape sample, the signals are in excellent agreement on all scales above 0.5 Mpc. At small scales, the number of pairs becomes very small, which makes the measurement noisier.

To validate our w_gg estimation we were satisfied by reaching agreement on b₁, b₂ with Sa23. We also use this opportunity to validate our analytical covariance matrix against the jackknife matrix estimates from LOWZ, we show the result in Fig. A.2. As can be expected from the correlation function plot, the agreement is good on all scales.

Fig. A.2. Validation of our covariance matrix estimate with the LOWZ sample from Reyes et al. (2012), Nakajima et al. (2012). We observe that our analytical prediction, the jackknife from Singh 2015, and our jackknife all agree very well, in particular on the scales of interest, while they deviate slightly below 1 Mpc.

Appendix B: Comparison to Lensfit

While the catalogue produced by ShapePipe (Farrens et al. 2022) is the main UNIONS catalogue, a second catalogue was produced by lensfit to cross-check any systematics. For example it can be seen in Li et al. (2023) that both agree very well when measuring lensing masses of AGN-black holes. The strategy adopted in the lensfit pipeline is described in Miller et al. (2013). In short, the pipeline adopts a likelihood-based method by fitting a bulge+disk model in a Bayesian framework, which includes nuisance parameters, such as the galaxy position, which are marginalised over.

Although it was not used for the main results (as the catalogue was not fully validated by the time of this work), we nevertheless verified that a good agreement exists as can be seen in Fig. B.1, where the CMASS-UNIONS data vectors generated with lensfit and ShapePipe are in good agreement on all scales. While this is not a sufficient check to diagnose all systematics, it increases our confidence in our results.

Fig. B.1. Comparison of lensfit and ShapePipe estimates of wg+. The good agreement increases confidence in the robustness of the measurement.

We note that it is not necessary for the w_g+ from the two different shape methods to be in agreement, as they are based on very different model-fitting functions. This is known to potentially result in different intrinsic alignment amplitudes, as explained throughout this paper. While we initially believed that we would be able to obtain further constraints from the comparison of these two shape measurement samples, recent work by MacMahon-Gellér & Leonard (2024) showed that a dedicated method was necessary to draw conclusions from a relative change in amplitude.

All Tables

All Figures

Fig. 1. Distribution of the different spectroscopic and lensing samples used in this work. The plot on the left shows the North Galactic Cap (NGC), where most UNIONS galaxies lie. On the right is a portion of the South Galactic Cap (SGC), with the DES Y3 observations in the most southern part. There is no overlap between the UNIONS and DES Y3 catalogues. The apparent variation in shades between the two panels is due to the two projections resulting in different point densities.

In the text

Fig. 3. Un-weighted normalised r-band magnitude distributions of the BOSS/eBOSS spectroscopic samples and UNIONS photometric galaxies. We can see that the spectroscopic galaxies from the BOSS/eBOSS surveys fall into the brighter tail of the UNIONS distribution. The galaxies which dominate the lensing signal are comparatively fainter. Note: we could not use our large overlap with LOWZ galaxies as these bright galaxies (r < 20) are mostly absent from our lensing sample due to size cuts.

In the text

	Fig. 2. Un-weighted normalised redshift distribution n(z) of the three BOSS/eBOSS galaxy samples. The UNIONS n(z) is approximate and for reference only, as blinding was implemented as a shift in the redshift distribution. The n(z) has been obtained by matching galaxies in the CFHTLenS W3 field and calibrated through deep spectroscopic reference samples.
In the text

Fig. 4. Square root of the diagonal of the analytical covariance matrix of wg+ for the CMASS-UNIONS sample (green solid line). The individual contributions are cosmic variance (dashed; the PggPII term is identical to the PgI contribution); mixed clustering – ellipticity shot noise (dotted); mixed intrinsic – density shot noise (long dash-dotted); and pure shot noise (short dash-dotted). See Eq. (42) and text for details. For comparison, we plot the jackknife variance (purple solid line), which is in excellent agreement with the analytical covariance over the relevant scales.

In the text

	Fig. 5. Measurements of wg+ (upper panels) and w++ (lower row) for the UNIONS-CMASS (left), UNIONS-LRG (middle), and UNIONS-ELG (right column) samples. The solid lines show the best-fit NLA model, where the shaded areas indicate the 1σ uncertainty associated with the fit. The grey areas are excluded from the fit as explained in Sect. 4.7.
In the text

	Fig. 6. Best-fit lines and 1σ contours for the NLA and TATT, respectively. The TATT model is fit down to 4 Mpc but both models appear very consistent. The 1σ TATT contours are obtained by drawing 500 models from the posterior and calculating the 68% confidence interval at eat rt.
In the text

	Fig. 7. Posterior contours of the NLA and TATT model fitted to the CMASS-UNIONS sample. We see no deviation from the NLA case as A2 and bTA appear consistent with 0.
In the text

	Fig. 8. Integrated correlation function, wg+, as measured on the CMASS-UNIONS sample split in 3 bins of luminosity, showing increasing luminosity from left to right panels. One can see a clear increase in the amplitude of the signal.
In the text

	Fig. 9. Intrinsic alignment amplitude as a function of luminosity as measured by different surveys. The fitted function is α⟨L/L0⟩β. We can see very similar scaling with luminosity when comparing our CMASS-DES and CMASS-UNIONS measurements.
In the text

	Fig. 10. 1D marginal posterior of the NLA intrinsic-alignment amplitude AIA of the same SDSS density tracer. Left panel: UNIONS, DES Y3 (both our measurements), and the results from Kurita & Takada (2023) in the North (NGC) and South Galactic Cap (SGC) using CMASS as density tracer. Right panel: Our measured UNIONS and DES Y3 results with the LRG sample.
In the text

	Fig. 11. Measurements of wg× for the CMASS-UNIONS, LRG-UNIONS and ELG-UNIONS samples. Deviations from 0 are characterised in Table 4 and are all compatible at the 1σ level.
In the text

	Fig. 12. Amplitude expressed by Eq. (50) of the additive bias contamination of the shape-shape correlation function. From left to right, we show our measurements using the Si15 data corresponding to the isophotal, de Vaucouleurs, and re-Gaussianisation shape catalogues, followed by UNIONS-CMASS, UNIONS-LRG, and UNIONS-ELG.
In the text

	Fig. 13. Relative contribution to the wg+ signal of PSF induced correlations. We see that on the scales of interest the contribution is below 2% which we judge satisfactory for current precision. The error bars blow up when wg+ is close to 0.
In the text

	Fig. A.1. Reproduction of the LOWZ measurement from Si15 for validation purposes. While we recovered a 2% agreement on the number of galaxies in the shape sample, the signals are in excellent agreement on all scales above 0.5 Mpc. At small scales, the number of pairs becomes very small, which makes the measurement noisier.
In the text

	Fig. A.2. Validation of our covariance matrix estimate with the LOWZ sample from Reyes et al. (2012), Nakajima et al. (2012). We observe that our analytical prediction, the jackknife from Singh 2015, and our jackknife all agree very well, in particular on the scales of interest, while they deviate slightly below 1 Mpc.
In the text

	Fig. B.1. Comparison of lensfit and ShapePipe estimates of wg+. The good agreement increases confidence in the robustness of the measurement.
In the text