UNIONS: Red and blue alignments of DESI DR1 galaxies

Introduction

Cosmic shear, the coherent distortion of the images of background galaxies by the gravitational lensing effect of large-scale structure, has emerged as one of the most powerful probes of cosmology. Recent stage-III surveys — the Dark Energy Survey (DES; [2] (DES Collaboration et al. 2026)), the Kilo-Degree Survey (KiDS; [59, 135]), and the Hyper Suprime-Cam survey (HSC; [24] (Li et al. 2023)), and the Ultraviolet Near-Infrared Optical Northern Survey (UNIONS; [49, 44, 46]) — have demonstrated that weak lensing measurements can constrain the amplitude of matter fluctuations to a precision of a few per cent. The upcoming stage-IV surveys, Euclid ([36]), the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; [64]), and the Nancy Grace Roman Space Telescope ([34]), will observe billions of galaxies, making previously sub-dominant systematic effects significant.

Among the most important systematics for cosmic shear is the intrinsic alignment (IA) of galaxies. As galaxies form and evolve within the cosmic web, tidal forces from the surrounding large-scale structure induce coherent orientations in their shapes ([16, 62, 23]). These intrinsic shape correlations mimic or contaminate the cosmic shear signal in two ways: through correlations between the intrinsic shapes of physically nearby galaxies (II correlations) and through correlations between the intrinsic shape of a foreground galaxy and the gravitational shear of a background galaxy (GI correlations). If not properly accounted for, intrinsic alignments can bias cosmological parameter inference from weak lensing surveys (for reviews, see [69, 126, 77, 20]).

Considerable observational effort has been devoted to characterising intrinsic alignments through direct measurements. A strong, radial IA signal has been detected for red, pressuresupported galaxies, whose shapes tend to stretch toward matter overdensities ([70, 117, 71, 43, 110, 57]). The amplitude of this signal depends on luminosity, following a broken power law with a steep rise for bright galaxies ([43, 57]), while Fortuna et al. [42] suggested that this broken power law behaviour was a consequence of the break in the stellar-to-halo mass relation; with halo mass being a more robust proxy of the IA signal. Blue, rotationally-supported galaxies, on the other hand, are expected to acquire spins correlated with the local tidal field via tidal torquing ([112, 22]). Observationally, blue galaxies show alignment signals consistent with zero at low and intermediate redshift ([86, 71, 108, 110, 57, 98]), though hydrodynamical simulations suggest a non-negligible signal may exist, especially at higher redshift ([19, 124, 22, 109, 138, 27]).

Several perturbative models have been developed to describe intrinsic alignments in cosmological analyses. The non-linear linear alignment model (NLA; [14, 62]) extends the linear alignment model of Catelan et al. [16] by replacing the linear matter power spectrum with its non-linear counterpart, an empirical prescription rather than a result of perturbation theory. The tidal alignment and tidal torque model (TATT; [12, 13]) further incorporates higher-order terms, including the quadratic tidal-torquing contribution that theory predicts for disc galaxies. The halo-mass picture of Fortuna et al. [42] has motivated mass-dependent variants of both models: Wright et al. [135] adopted NLA-M in the KiDS-Legacy cosmic shear analysis, parametrising the IA amplitude as a function of halo mass, and Herle et al. [54] introduced TATT-M, calibrated on the FLAMINGO hydrodynamical simulations, which is strongly preferred over the standard NLA model. More recently, effective field theory (EFT) approaches have been proposed to extend the model validity to smaller scales ([131, 132, 6, 85]). NLA and TATT are already marginalised over in cosmological analyses, while the EFT models have so far mostly been validated on simulations, with a first application to cosmic-shear data by Chen et al. [18].

To date, however, most direct IA measurements have focused on bright, red galaxies selected from spectroscopic surveys. This focus on luminous red objects introduces a demographic bias in our understanding of intrinsic alignments, leaving the IA signal of fainter and bluer populations — which dominate the source samples of cosmic shear surveys — poorly constrained. The colour dependence of IA matters for cosmic shear itself: early-type galaxies dominate the IA contamination ([58]), colour-split source samples have served as internal-consistency tests of cosmic-shear analyses ([123]), and restricting the sources to blue galaxies, as proposed by Krause et al. [73] and done for DES Y3 (McCullough et al. [90]), reduces the sensitivity of the cosmological constraints to the IA model. These analyses still marginalise over the alignment of the blue galaxies, so direct measurements of it provide the priors needed to apply such selections in Euclid and other stage-IV surveys. The combination of spectroscopy from the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration [30]) with stage-III lensing data from UNIONS, and ultimately with stage-IV data from Euclid, offers a unique path toward this goal.

Most direct IA measurements to date have used projected two-point correlation functions ([87, 70, 117, 110]). Singh et al. [119], hereafter Si23, recently demonstrated that a multipole-based estimator applied to the three-dimensional correlation function in redshift space provides significantly tighter constraints on IA model parameters than projected estimators, by a factor of 1.5–2 depending on the sample. The multipoles keep the full three-dimensional statistic, whereas $w_{g+}$ integrates it along the line of sight up to a fixed $\Pi_{\max}$ and loses information in the projection. Three-dimensional, redshift-space treatments of alignments have also been developed by Kurita & Takada [74] and Laman et al. (2025).

In this work, we cross-correlate UNIONS galaxy shapes with DESI DR1 spectroscopic redshifts, as Siegel et al. [116], hereafter S26, did with four other imaging surveys, but using the multipole estimators of Si23. We measure the clustering monopole and the IA quadrupole of the DESI Bright Galaxy Survey (BGS), luminous red galaxy (LRG) and emission-line galaxy (ELG) samples, split by galaxy colour and stellar mass, and fit them with the NLA model, and with TATT for the blue samples. Following Fortuna et al. [42], we then ask whether the alignment amplitude of red galaxies is set by halo mass: we measure the galaxy–galaxy lensing signal of UNIONS sources around the BGS red galaxies, fit it with a conditional stellar mass function, and use the resulting stellar-to-halo mass relation to assign an effective halo mass to every alignment sample. This tests the halo-mass scaling assumed by NLA-M and TATT-M down to the low-mass, faint end of the red population. This paper is organised as follows. In Sect. 2, we describe the UNIONS and DESI DR1 datasets and the alignment samples, and in Sect. 3 the correlation-function, covariance and lensing measurements. § 3 presents the IA models and the likelihood. We show our results for red and blue galaxies, including the halo-mass dependence of the alignment amplitude, in Sect. 5, discuss them in Sect. 6, and conclude in Sect. 7.

Data

The overlap between the UNIONS shape catalogue and DESI DR1 is shown in Figure 1: measured from the DESI random catalogues, the DESI DR1 samples overlap the UNIONS shape catalogue over 2397 deg$^{2}$ (BGS), 1204 deg$^{2}$ (LRG) and 1255 deg$^{2}$ (ELG), the largest northern overlap of a stage-III lensing survey with DESI.

UNIONS shape catalogue

The Ultraviolet Near-Infrared Optical Northern Survey (UNIONS; Gwyn et al. [49]) is a collaboration of wide-field imaging surveys combining deep observations from the Canada-France-Hawaii Telescope (CFHT), Pan-STARRS, and Subaru in the Northern hemisphere. Its multi-band ($u griz$) photometry serves a broad range of science, from weak lensing ([106, 89] (Li et al. 2024; Mpetha et al. 2025)) and strong lensing ([17, 111]) to galaxy evolution ([35, 11, 41]) and the Milky Way ([68, 122]), and supports photometric redshift estimation for the Euclid mission. In this work, we use galaxy shape measurements from the CFHT $r$-band data, obtained through the MegaPipe image processing pipeline (Gwyn [50]).

We use an extended version of the UNIONS shape catalogue of Hervas-Peters et al. [56].¹ Galaxy shapes are measured using the ShapePipe pipeline (Guinot et al. [48]; Farrens et al. [39]), which employs the ngmix software (Sheldon [113]) for model fitting and the Metacalibration framework (Sheldon & Huff [114]) for shear calibration. Metacalibration provides a self-calibrating approach to shear estimation by artificially shearing galaxy images and measuring the response of the shape estimator, yielding a shear response $\mathcal{R}$ that accounts for model bias, noise bias, and selection effects. Objects are detected and deblended with SExtractor (Bertin & Arnouts [10]) and the point spread function (PSF) is modelled using PSFEx (Bertin [9]).

Restricted to Dec $\ge22^\circ$, the version of the shape catalogue we use covers 5400 deg$^2$, of which 4800 deg$^2$ lie in the Northern Galactic Cap and the rest in the Southern Galactic Cap, with a median seeing of 0.69″ in the $r$-band. This is 1900 deg$^2$ more than the 3500 deg$^2$ of the catalogue used for the UNIONS cosmic-shear analyses (Guerrini et al. [46]; Goh et al. [44]) and for the previous UNIONS alignment measurement with BOSS/eBOSS (Hervas Peters et al. [57]). Each galaxy is assigned an inverse-variance shape weight $w_{\mathrm{SP}} = (2\sigma_{\mathrm{SN}}^2 + \sigma_{e_i}^2)^{-1}$, where $\sigma_{\mathrm{SN}}$ is the galaxy shape noise and $\sigma_{e_i}$ is the standard deviation of the ellipticity component $i = 1, 2$ estimated by the model fitting procedure.² For a detailed description of an earlier version (v1.4.6.3) of the catalogue and its systematics validation, we refer to Hervas-Peters et al. [56] and Daley et al. [25].

DESI DR1 spectroscopic catalogues

The Dark Energy Spectroscopic Instrument (DESI Collaboration [30]) is a multi-object fibre spectrograph mounted on the 4-metre Mayall telescope at Kitt Peak National Observatory, capable of simultaneously observing 5000 targets. The first data release (DR1; DESI Collaboration et al. 2026) covers the first year of the main survey and provides spectroscopic redshifts for millions of galaxies and quasars over 9500 deg$^2$ in each of the bright and dark programmes.

We use the large-scale structure (LSS) catalogues constructed for the DESI DR1 cosmological analyses (version 1.5; Adame et al. [1]; Ross et al. [107]), which include matched random catalogues and the standard clustering weights. We weight every density sample with the clustering weight provided in these catalogues, $w_{\mathrm{DESI}} = w_{\mathrm{comp}} w_{\mathrm{zfail}} w_{\mathrm{sys}}$, which combines the completeness weight $w_{\mathrm{comp}}$, correcting for fibre-assignment incompleteness, the redshift-failure weight $w_{\mathrm{zfail}}$, and the imaging-systematics weight $w_{\mathrm{sys}}$ (see Adame et al. [1]; Ross et al. [107], for details). The total weight assigned to each galaxy in the density sample is

where $w_{\mathrm{FKP}}(z) = (1+\bar{n}(z)P_0)^{-1}$ is the FKP weight (Feldman et al. [40]) with $P_0$ the power spectrum amplitude at $k \approx0.15\,h\,\mathrm{Mpc}^{-1}$, which optimises the signal-to-noise of clustering measurements. The random catalogues carry the same weights as the data. DESI DR1 provides no imaging-systematics weights for BGS ($w_{\mathrm{sys}} = 1$ in its LSS catalogues), so for BGS we recompute $w_{\mathrm{DESI}}$ in exactly the same way as Ross et al. [107], but with our own $w_{\mathrm{sys}}$, derived following the same procedure as Adame et al. [1]: a linear regression of the galaxy density against the imaging-property maps. The BGS random catalogues are rebuilt accordingly. The correction modulates the density by 2%, and we verified that removing it changes the galaxy–shape quadrupole by less than $0.1\sigma$ in every separation bin. We consider the following tracer populations:

Bright Galaxy Survey (BGS). We use the BGS_ANY sample (Hahn et al. [51]), which extends beyond the flux-limited BGS_BRIGHT selection ($r < 19.5$) to include fainter targets observed during bright-time conditions.

Luminous Red Galaxies (LRGs). The LRG sample (Zhou et al. [137]) is selected based on photometric colour and magnitude cuts designed to isolate massive, red, early-type galaxies at intermediate redshift, $0.4 < z < 1.1$. They have been extensively studied in previous IA measurements (Hirata et al. [61]; Joachimi et al. [70]; Singh et al. [117]; Samuroff et al. [110]).

Emission Line Galaxies (ELGs). The ELG sample (Raichoor et al. [105]) selects star-forming galaxies whose redshifts are measured from their [O II] emission; the DR1 clustering catalogues retain $0.8 < z < 1.6$ [1, 107]; we use the ELG_LOPnotqso LSS catalogue, i.e. the main (LOP) ELG selection with the objects also targeted as quasars removed. These are blue, disc-dominated galaxies for which intrinsic alignments are expected to be weak. We use the ELG sample primarily as a null test for the IA signal.

Alignment samples

To construct our analysis samples, we cross-match the DESI DR1 spectroscopic catalogues with the UNIONS shape catalogue. The DESI galaxies serve as the density (position) tracers, while the subset of DESI galaxies with UNIONS shape measurements constitutes the shape sample. We use a $k$-d tree matching in astropy with a tolerance of $1^{\prime\prime}$ to identify common objects, following the procedure of Hervas Peters et al. [57]. No shape-quality selection is applied to this IA shape sample beyond the cross-match itself: the size, magnitude and signal-to-noise cuts of the lensing selection mainly serve to reject stars and spurious detections, which the match to spectroscopically confirmed galaxies already excludes, and no PSF-leakage correction, which Hervas Peters et al. [57] showed to have a negligible impact on the IA statistics; the ellipticities of each sample are calibrated by its mean Metacalibration response $R$ (Tables B.2 and 2).

We further split the BGS sample by colour to study the colour dependence of the IA signal. Red and blue sub-samples are defined using a Gaussian mixture model (GMM) in the $(r, g-r)$ magnitude–colour plane. Galaxies are binned in eight redshift slices with edges at $z = 0.01, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40$ and $0.50$ (the GMM is fitted on the full BGS catalogue, including the $z < 0.1$ galaxies that are not used in the IA samples), and in each slice a three-component GMM is fitted. The components with the reddest and bluest mean colour define the red and blue sequences; their mean colours and widths are interpolated in redshift with a cubic spline, held constant beyond the first and last slice centres, and each galaxy is assigned to the sequence it is closer to, in units of that sequence’s width, provided it lies within $2\sigma$ of it. 1% of the galaxies satisfy neither condition and are left out of both sub-samples. The magnitudes used to perform this step are taken from the full DESI catalogue, and are corrected from galactic extinction. The resulting colour–redshift distribution of the BGS red, BGS blue, and LRG samples is shown in Figure 2. BGS red galaxies and LRGs are both red but are selected differently: the former are drawn from a flux-limited sample at $z < 0.5$ by their position on the red sequence, the latter by optical–infrared colour cuts designed to yield an approximately constant comoving number density of massive galaxies (Zhou et al. [137]). The two samples overlap in stellar mass only above $\log(M_\star/h^{-2}M_\odot) = 10.9$ (Table B.2).

Stellar masses are taken from the FastSpecFit value-added catalogue (Moustakas et al. [95]), which fits the DESI spectra jointly with the Legacy Imaging Surveys photometry using stellar-population-synthesis templates at the measured spectroscopic redshift. We use its LOGMSTAR estimate throughout. FastSpecFit adopts $H_0 = 100\ {\rm km\,s^{-1}\,Mpc^{-1}}$ internally, so these masses are in units of $h^{-2}M_\odot$, and all stellar masses quoted in this work carry that convention.

Measurements

Correlation functions

We estimate the three-dimensional correlation functions on two grids of pair separations: $(r_p,\Pi)$, where $r_p$ is the transverse (projected) separation and $\Pi$ the line-of-sight separation; and $(s,\mu)$, where $s = \sqrt{r_p^2+\Pi^2}$ is the three-dimensional separation and $\mu= \Pi/s$ is the cosine of the angle between the separation vector and the line of sight. The clustering and intrinsic-alignment statistics are measured on two different logarithmic grids, both spanning 6 to 100 $h^{-1}$ Mpc in $s$ and in $r_p$: a fine grid of 19 bins for the clustering statistics ($\tilde{\xi}_{gg,0}$ and $w_{gg}$), and a coarser grid of 10 bins for the shape cross-correlations ($\tilde{\xi}_{g+,2}$ and $w_{g+}$). In the $(r_p,\Pi)$ frame, $\Pi$ is binned into 20 linear bins of width 10 $h^{-1}$ Mpc over $\Pi\in[-100,100]\ h^{-1}$ Mpc; in the $(s,\mu)$ frame, $\mu$ is binned into 20 linear bins of width 0.01 over $\mu\in[-1,1]$.

The galaxy–galaxy correlation function is estimated using the Landy & Szalay [78] estimator:

where D and R denote the data and random catalogues, respectively. All pair counts are normalised.

The galaxy–shape cross-correlation function is estimated with the estimator of Mandelbaum et al. [87]:

where $S_+\mathrm{D} = \sum_{i\in\mathrm{S},j\in\mathrm{D}}w_iw_je_+(i|j)$ is the weighted sum of ellipticity component $e_+(i|j)$ of shape galaxy $i$ measured relative to the direction towards density galaxy $j$, $S_+\mathrm{R}_{\mathrm{D}}$ is the analogous sum using random points drawn from the density catalogue, and $R_{\mathrm{S}}R_{\mathrm{D}}$ is the random–random normalisation.³ The subtraction of the shape–random term removes additive systematic biases arising from correlations between the shape measurement and the survey geometry. The component $e_+$ is positive when the major axis of the shape galaxy points towards the density galaxy (radial alignment) and negative when it is perpendicular to that direction (tangential alignment).

The multipole pair counts are performed with the GPU-accelerated cucount code⁴, and the projected statistics with pycorr (galaxy–galaxy) and TreeCorr (galaxy–shape). We use the DESI random catalogues, which match the angular and radial selection of each tracer sample (rebuilt for BGS as described in Sect. 2.2). For the density sample, we use the total weight $w_{\mathrm{tot}}$ (Eq. (1)). For the shape sample, each galaxy is weighted by the product of its shape weight $w_{\mathrm{SP}}$ (Sect. 2.1) and the density completeness, redshift-failure and systematic weights $w_{\mathrm{comp}}$, $w_{\mathrm{zfail}}$, $w_{\mathrm{sys}}$, excluding the FKP weight.

Following Si23, we decompose the anisotropic correlation functions into multipole moments. For a general correlation between fields $a$ and $b$ with spins $s_{ab}$, the multipole of order $\ell$ is

where $\mu_r = \Pi/r$ is the cosine of the angle between the separation vector and the line of sight, $\Pi$ is the line-of-sight component, and $L_{\ell}^{\ell,s_{ab}}$ are associated Legendre polynomials. We denote the measured multipoles by $\tilde{\xi}_{ab,\ell}$, and use two of them: the clustering monopole $\tilde{\xi}_{gg,0} \equiv\xi_{gg}$ and the galaxy–shape quadrupole $\tilde{\xi}_{g+,2} \equiv\xi_{g+}$, both as functions of the redshift-space separation $s$. The projected correlation functions $w_{gg}(r_p)$ and $w_{g+}(r_p)$ are obtained by summing the $(r_p,\Pi)$ measurements over $|\Pi| < 100\,h^{-1}\,\mathrm{Mpc}$.

We do not apply an additional transverse cut on $r_p$ when building the multipoles. Si23 discarded pairs with $r_p < 5\,h^{-1}\,\mathrm{Mpc}$ from their multipole estimators, in order to suppress the contribution of redshift-space distortions along the line of sight. Applying the same cut to the low-redshift LRG sample changes the amplitude of $\tilde{\xi}_{g+,2}$ over the fitted range ($s > 10\,h^{-1}\,\mathrm{Mpc}$) by 7%, i.e. $0.2\sigma$ of $A_{\mathrm{IA}}$, and its jackknife uncertainties by less than 3%. We therefore retain all pairs.

Covariance

The covariance of all statistics is estimated by jackknife resampling: the overlapping survey area is divided into $N_{\mathrm{jk}} = 70$ equal-area regions by running the $K$-means algorithm on the random catalogue of each tracer, and $N_{\mathrm{jk}}$ correlation functions are computed, each with one region removed. This number is set so that a BGS region, $6^\circ$ across, spans the largest fitted separation, $100\,h^{-1}\,\mathrm{Mpc}$, at the effective redshift of the BGS density sample; the LRG and ELG regions are larger, as these samples lie at higher redshift. The resulting covariance is corrected for the cross-region pairs (pairs whose two members lie in different regions) following Mohammad & Percival [94]. We do not use the cross-covariance between $\tilde{\xi}_{gg,0}$ and $\tilde{\xi}_{g+,2}$ in the likelihood, as the limited number of jackknife regions makes the off-diagonal blocks noisy; Si23 noted that neglecting it should not have a large impact, since these covariances are dominated by noise.

Galaxy–galaxy lensing

To assign halo masses to the alignment samples (Sect. 5.1.3), we measure the excess surface mass density $\Delta\Sigma(r_p)$ around the BGS red lens samples with the dsigma package (Lange & Huang [81]). The sources are UNIONS galaxies (Hervas-Peters et al. [56]) selected in photometric redshift well behind the $z < 0.5$ lenses, and their redshift distribution is calibrated with clustering redshifts against the DESI spectroscopic samples (d’Assignies et al. 2025). The source selection, the estimator and the calibration are described in Appendix C, together with the stellar-to-halo mass relation derived from this measurement.

Modelling and inference

Our baseline is the non-linear linear alignment (NLA) model (Catelan et al. [16]; Hirata & Seljak [62]; Bridle & King [14]), in which the intrinsic shape field responds linearly to the large-scale tidal tensor $s_{ij}$, computed from the non-linear rather than the linear matter power spectrum, with a single free amplitude $A_{\mathrm{IA}}$, complemented by a linear galaxy bias $b_1$ and the linear-order Kaiser prescription for redshift-space distortions (Kaiser [72]). In our convention a positive $A_{\mathrm{IA}}$ corresponds to shapes elongated along the stretching direction of the tidal field, i.e. pointing towards overdensities (radial alignment), as expected for pressure-supported elliptical galaxies, and a negative $A_{\mathrm{IA}}$ to shapes perpendicular to it (tangential alignment). For the blue samples, for which the tidal-torquing term is expected to be most relevant, we also fit the tidal alignment and tidal torque (TATT) model (Blazek et al. [13]), which expands the intrinsic shape field to second order in $s_{ij}$ and in the matter overdensity $\delta$:

where $C_1$ describes the linear response to the tidal field (tidal alignment), $C_{1\delta}$ the modulation of that response by the local galaxy density, and $C_2$ the quadratic tidal-torquing contribution, which arises when the angular momentum of disc galaxies is acquired from tidal torques and their shapes follow their spin. Their relation to the amplitudes $A_{\mathrm{IA}}$, $b_{\mathrm{TA}}$ and $A_2$ quoted in this paper is given in Appendix A.§ 1. The NLA model retains only the linear term ($b_{\mathrm{TA}} = A_2 = 0$). In the TATT fits we set $b_{\mathrm{TA}} = 1$, keeping $b_1$, $A_{\mathrm{IA}}$ and $A_2$ free. The corresponding real- and redshift-space power spectra are given in Appendix A. In the following we consider only the linear bias approximation: the higher-order terms arising from non-linear galaxy bias in the power spectra are not included (Blazek et al. [13]; Euclid Collaboration: Paviot et al. [37]).

Model prediction

The model predictions for each multipole (Eq. (4)) are obtained through Hankel transforms of the redshift-space power spectra (Eqs. (A.3) and (A.4) in Appendix A):

where $j_\ell$ are spherical Bessel functions and the coefficients $\alpha_{\ell,ab}$ encode the redshift-space distortion contributions. The projected correlation functions follow from the same decomposition by integrating the multipole expansion along the line of sight (Eq. 19 of Si23):

where $r = \sqrt{r_p^2+\Pi^2}$, which accounts for the residual effects of redshift-space distortions at finite $\Pi_{\max}$. For galaxy clustering ($s_{ab}=0$), the non-vanishing multipoles are the monopole ($\ell= 0$), quadrupole ($\ell= 2$), and hexadecapole ($\ell= 4$); for the galaxy–shape cross-correlation ($s_{ab} = 2$), they are the spin-2 quadrupole ($\ell= 2$) and hexadecapole ($\ell= 4$). The corresponding coefficients $\alpha_{\ell,ab}$ are given in Appendix A.§ 2 (Eqs. (A.5)–(A.9)). We evaluate all correlation functions at point separations (no bin averaging) and use $\Pi_{\max} = 100\,h^{-1}\,\mathrm{Mpc}$, as in the measurements. We neglect the gravitational-lensing contamination of the galaxy–shape correlation, i.e. the shear of the shape galaxies by the matter around the density galaxies and the magnification of either sample (the $gG$, $\mu I$ and $\mu G$ terms of Samuroff et al. [110]): adding these terms to the model, with the magnification coefficients measured for our samples with the method of Wenzl et al. [133] (Appendix C), shifts $A_{\rm IA}$ by less than $0.1\sigma$ in every sample, because both samples are spectroscopic and the contamination only becomes appreciable at the largest separations, where the quadrupole carries little weight.

Fitted statistics and scale cuts

We fit the clustering monopole $\tilde{\xi}_{gg,0}$ and the galaxy–shape quadrupole $\tilde{\xi}_{g+,2}$ jointly, and, as a cross-check, the projected $w_{gg}$ and $w_{g+}$ (Eq. (7)). We adopt conservative scale cuts, summarised in Table 1, to exclude the small scales where non-linear clustering and Fingers-of-God effects are not captured by our linear-order model. The clustering monopole, which is the more sensitive to these effects, is cut at $30\,h^{-1}\,\mathrm{Mpc}$ on its finer grid (Sect. 3.1), well inside the linear regime, and the shape quadrupole at $10\,h^{-1}\,\mathrm{Mpc}$. For the TATT fits the quadrupole cut is relaxed to $6\,h^{-1}\,\mathrm{Mpc}$, since the tidal-torquing term contributes predominantly on small scales. The projected cross-check uses the cut commonly adopted in previous analyses, $r_p > 6\,h^{-1}\,\mathrm{Mpc}$ for both $w_{gg}$ and $w_{g+}$ ([70, 117, 71]), which admits pairs down to a three-dimensional separation $s = 6\,h^{-1}\,\mathrm{Mpc}$ (since $s \ge r_p$); Si23 fitted their multipoles from $s = 5\,h^{-1}\,\mathrm{Mpc}$. Our baseline multipole cuts are therefore more conservative than both.

Statistic	Model	$r_{\min}$ [$h^{-1}\,\mathrm{Mpc}$]	$N_{\mathrm{pts}}$
$\tilde{\xi}_{gg,0}$	NLA, TATT	30	8
$\tilde{\xi}_{g+,2}$	NLA	10	8
$\tilde{\xi}_{g+,2}$	TATT	6	10
$w_{gg}$	NLA	6	19
$w_{g+}$	NLA	6	10

Likelihood and parameter inference

The cosmology is fixed to the Planck Collaboration et al. [104] best-fit values ($h = 0.6766$, $\Omega_{\mathrm{c}}h^2 = 0.11933$, $\Omega_{\mathrm{b}}h^2 = 0.02242$, $A_{\mathrm{s}} = 2.105 \times10^{-9}$, $n_{\mathrm{s}} = 0.9665$).

We assume a Gaussian likelihood

where $d_\alpha$ and $m_\alpha$ are the data and model vectors for the monopole ($\alpha= 0$) and quadrupole ($\alpha= 2$), respectively, and $C_\alpha$ is the corresponding jackknife covariance block (Sect. 3.2). We do not apply the Hartlap et al. [53] correction to the inverse covariance, which is derived for covariances estimated from independent realisations, such as mocks, and does not strictly apply to jackknife resamples of a single data set.

The IA model, the likelihood and the fits are implemented in the public IA2PT package.⁵ We sample the posterior with the Nautilus nested sampler (Lange [79]), using flat priors $b_1 \in[0.1,4]$ and $A_{\mathrm{IA}} \in[-10,10]$ and $A_2 \in[-6,6]$ for the TATT fits. The sampler is run until the effective sample size reaches $n_{\mathrm{eff}} = 50\,000$. Quoted uncertainties are the $1\sigma$ marginalised posterior widths, and significances are the posterior median in units of this width. The theoretical predictions are computed using PyCCL (Chisari et al. [21]) with CAMB (Lewis et al. [82]) and HM-Code2020 (Mead et al. [92]) for the non-linear matter power spectrum, and Fast-PT (McEwen et al. [91]; Fang et al. [38]) for the perturbation-theory integrals. Hankel transforms from Fourier to configuration space are performed using the routines in Fast-PT. The clustering monopole is evaluated at the effective redshift of the density sample, and the galaxy–shape quadrupole (and $w_{g+}$) at the pair-weighted effective redshift of the density–shape cross-correlation, $z_{\mathrm{IA}}$, listed in the tables of Appendix § 8. The same $b_1$ is used at both redshifts, without evolution. Evolving it from the density to the pair redshift at constant $b_1D(z)$ instead changes $A_{\mathrm{IA}}$ by at most $6\%$ ($0.6\sigma$, for the full BGS red sample; median shift $0.08\sigma$) and leaves the $\chi^2$ unchanged.

Results

All fits use the NLA model applied to $\tilde{\xi}_{gg,0}$ and $\tilde{\xi}_{g+,2}$, except the TATT fits of Sect. 5.2.3; we interpret the results in Sect. 6.

Red galaxies

As expected, red galaxies are radially aligned. Their best-fit parameters are collected in Table B.2, and their amplitudes are shown as functions of stellar mass and halo mass in Figure 3. The correlation functions of the BGS red mass bins are shown in Fig. B.3 (Appendix § 8).

BGS red

The density sample, common to the BGS red and blue fits, is the full BGS clustering catalogue, whose linear bias is $b_1 = 1.39 \pm0.04$.

We split the BGS red sample into twelve stellar-mass bins built to yield comparable expected signal-to-noise in $\tilde{\xi}_{g+,2}$, which makes them equipopulated at low mass and progressively sparser at high mass.⁶ The mass dependence is non-monotonic (Figure 3): the amplitude decreases from the lowest-mass bins to a dip at $10.6 \le\log(M_\star/h^{-2}M_\odot) \le11.2$, and rises steeply in the two highest-mass bins, up to $A_{\mathrm{IA}} = 5.7 \pm0.6$ at $\log(M_\star/h^{-2}M_\odot) = 11.4$.

LRG

The LRG sample is split into two redshift bins: a low-redshift bin at $0.4 < z < 0.75$ and a high-redshift bin at $0.75 < z < 1.1$. We split each LRG redshift bin into five stellar-mass bins, in a similar way as for the BGS sample. The amplitude increases with mass in both redshift bins, up to $A_{\rm IA}=3$–$4$ in the most massive bin ($\log(M_\star/h^{-2}M_\odot)=11.4$). At the same median stellar mass, the high-redshift LRGs are more strongly aligned than the low-redshift ones: the global amplitudes ($A_{\rm IA}=1.17\pm0.35$ and $2.95\pm0.38$) differ at $3.4\sigma$, and every high-redshift mass bin lies above its low-redshift counterpart. The low-redshift bin at $\log(M_\star/h^{-2}M_\odot)=11.1$ is consistent with zero, while the low-redshift LRG bins otherwise agree with the BGS red bins of similar mass (Figure 3). We discuss this redshift difference in Sect. 6.2.

Halo-mass dependence

Figure 3 also shows the amplitudes of the BGS red and LRG bins against the effective halo mass $M_{\rm eff}$ of each bin. Following Fortuna et al. [42], $M_{\rm eff}$ is the mean mass of the haloes hosting a central galaxy,

where $n(M,z)$ is the halo mass function at redshift $z$ and $\langle N_{\rm c}|M,z\rangle$ the central occupation of the best-fit conditional stellar mass function (Appendix C), integrated over the stellar-mass range of the bin. Satellites therefore do not enter $M_{\rm eff}$; halo masses are defined at 200 times the mean matter density and quoted in $h^{-1}M_\odot$. The relation between median stellar mass and $M_{\rm eff}$ is established on the volume-limited BGS red lens bins and applied to the alignment samples through their median stellar mass; for the LRG samples this is an extrapolation in both stellar mass and redshift.

We fit power laws to the measured $A_{\rm IA}$ and their uncertainties with iminuit (Dembinski et al. [28]). A single power law,

the form adopted by Fortuna et al. [42], provides a poor description of our measurements over the full mass range, with $\chi^2=96.8$ for 20 degrees of freedom (hereafter $\chi^2/\mathrm{dof}=96.8/20$). Restricting it to $M_{\rm eff}>10^{12}\,h^{-1}M_\odot$ gives $\beta_{\rm M}=0.66\pm0.18$, but the fit remains poor ($\chi^2/\mathrm{dof}=80.2/19$), mostly because the three lowest-mass BGS red bins above $10^{12}\,h^{-1}M_\odot$ lie well above the power law and the low-redshift LRG bin that is consistent with zero lies well below it. A broken power law,

continuous at a free pivot mass $M_{\rm b}$, improves the fit substantially ($\chi^2/\mathrm{dof}=32.3/18$), with a nearly flat faint end ($\beta_1=-0.24\pm0.09$) below $M_{\rm b}=1.3\times10^{13}\,h^{-1}M_\odot$ and a steep rise above it ($\beta_2=1.74\pm0.29$; Figure 3).

The low-mass behaviour does not depend on the extrapolation to the LRGs. Restricting the fits to the BGS red samples alone, for which the mapping is direct⁷, the single power law remains strongly disfavoured ($\chi^2/\mathrm{dof}=69.5/10$) while the broken power law gives a good description of the data ($\chi^2/\mathrm{dof}=8.2/8$), with the same faint-end slope and pivot as the joint fit.

Blue galaxies

Unlike the red galaxies, the blue samples show no radial alignment, and the massive BGS blue galaxies are tangentially aligned. The NLA fits are collected in Table 2, the amplitudes are shown as a function of stellar mass in Figure 4, and the correlation functions of the BGS blue mass bins in Figure 5. The multipole and projected ($w_{g+}$) fits agree in every sample, and the multipoles are consistently the more precise: the uncertainty on $A_{\rm IA}$ from $w_{g+}$ is 43–78% larger than from the multipoles on the same data (last column of Table 2).

Sample	$z$ range	$z_{\mathrm{IA}}$	$\log M_\star$	$N_{\mathrm{shapes}}$	$\sigma_e$	$\mathcal{R}$	$b_1$	$A_{\mathrm{IA}}^{\widetilde{\xi}_{g+,2}}(\chi_\nu^2)$	$A_{\mathrm{IA}}^{w_{g+}}(\chi_\nu^2)$	Gain
BGS blue
Global	0.1–0.5	0.25	$10.34 \pm0.53$	827 884	0.226	0.92	$1.390 \pm0.040$	$-0.291 \pm0.113$ (1.80)	$-0.328 \pm0.166$ (1.21)	+47%
$M_\star$ bin 1 (low)	0.1–0.5	0.20	$9.94 \pm0.33$	405 662	0.232	0.93	–	$-0.049 \pm0.111$ (1.89)	$0.069 \pm0.197$ (1.40)	+78%
$M_\star$ bin 2 (high)	0.1–0.5	0.30	$10.68 \pm0.23$	405 663	0.218	0.90	–	$-0.592 \pm0.180$ (1.52)	$-0.585 \pm0.266$ (1.15)	+48%
ELG
Global	0.8–1.6	1.11	$9.43 \pm0.65$	414 759	0.299	0.72	$1.215 \pm0.042$	$-0.375 \pm0.533$ (1.76)	$-0.535 \pm0.760$ (1.30)	+43%

BGS blue

For the full BGS blue sample we measure $A_{\rm IA}=-0.29\pm0.11$ from the multipole analysis, a mild tangential alignment at $2.6\sigma$, consistent with the projected statistic ($A_{\rm IA}=-0.33\pm0.17$).

We split the BGS blue sample into two stellar-mass bins. The low-mass bin is consistent with zero ($A_{\rm IA}=-0.05\pm0.11$). The high-mass bin is tangentially aligned, with $A_{\rm IA}=-0.59\pm0.18$ ($3.3\sigma$), which the projected statistic confirms at lower significance ($2.2\sigma$); to our knowledge, this is the most significant detection of a non-zero alignment of blue galaxies to date. The mild tangential signal of the global sample is thus driven entirely by the high-mass half of the blue population. The detection does not rely on the NLA model. Without any model, the eight fitted quadrupole points of the high-mass bin are inconsistent with zero at $2.2\sigma$ ($\chi^2=16.9$, $p=0.03$). The amplitude is Notes. $z_{\mathrm{IA}}$ as in Table B.2; both $\widetilde{\xi}_{g+,2}$ and $w_{g+}$ are modelled at $z_{\mathrm{IA}}$, the monopole and $w_{gg}$ at the density-sample redshift (0.35 (BGS blue), 1.18 (ELG)). The last column is the gain in precision on $A_{\mathrm{IA}}$ delivered by the multipole analysis relative to the projected one, $\sigma_{w_{g+}}/\sigma_{\widetilde{\xi}_{g+,2}} - 1$.

also stable against the quadrupole scale cut: for $r_{\min}$ between 6 and 20 $h^{-1}\,\mathrm{Mpc}$ it stays tangential and within the uncertainty of the fiducial value (Fig. B.1).

The reduced $\chi^2$ of the multipole fits (Table 2) are high, but the excess comes from the BGS blue monopole that the three samples share rather than from the alignment signal (Appendix § 8).

ELG

The ELG sample at $0.8 < z < 1.6$ yields $A_{\mathrm{IA}} = -0.37 \pm0.53$ from the multipole analysis, consistent with zero and with the projected statistic, and constrains the ELG alignment amplitude to $-1.4 < A_{\mathrm{IA}} < 0.7$ at $2\sigma$. The multipole fit is somewhat poor ($\chi_\nu^2 = 1.76$), with the excess shared equally between the monopole and the quadrupole (Appendix § 8). The ELG correlation functions are shown in Fig. B.2 (Appendix § 8).

TATT fits

To test whether the tangential signal requires a tidal-torquing contribution, we repeat the NLA multipole fits of the blue samples with the full TATT model (Sect. 4), relaxing the quadrupole scale cut to $r_{\min} = 6\,h^{-1}\,\mathrm{Mpc}$ (Table 1). For the global BGS blue sample we obtain $A_{\mathrm{IA}} = -0.85 \pm0.33$ and $A_2 = 1.73 \pm0.75$, and the same pattern, a negative $A_{\mathrm{IA}}$ and a positive $A_2$, holds in both stellar-mass bins (Table B.1). The galaxy bias is unchanged from the NLA fits.

These marginal constraints should not be read as two independent measurements. The two amplitudes are almost perfectly anti-correlated ($\rho= -0.99$; Figure 6), because with the normalisation of Eq. (A.1) a positive $A_2$ produces a radial $\widetilde{\xi}_{g+,2}$, like a positive $A_{\mathrm{IA}}$: around a spherical overdensity, the quadratic term aligns the shapes with the direction of largest tidal eigenvalue in absolute value, which is the radial one. The well-measured quantity is therefore $A_{\mathrm{IA}} + c A_2$, with $c = 0.4$ from the posteriors, which is determined six to eight times more precisely than $A_{\mathrm{IA}}$ alone. This combination reproduces the NLA results: null at low mass ($0.00 \pm0.05$) and tangential at high mass ($-0.19 \pm0.07$).

In our fits the tangential signal is carried by the negative tidal-alignment amplitude, while the positive $A_2$ partly cancels it on small scales; the positive marginal $A_2$ is therefore not evidence for a tidal-torquing origin of the tangential alignment. Fitted to the same scales ($r_{\min} = 6\,h^{-1}\,\mathrm{Mpc}$), TATT lowers the $\chi^2$ of the NLA model by only 2.8–3.9 for one additional parameter, and the Bayesian evidence does not favour it ($\ln B_{\mathrm{TATT/NLA}} = -0.4$ to $-0.7$): the data do not require the tidal-torquing term. Conversely, a torque-only fit ($A_{\mathrm{IA}} = 0$, $A_2$ free) recovers the tangential signal of the high-mass bin as a negative $A_2 = -0.38 \pm0.15$, but describes it less well than NLA on the same scales ($\Delta\chi^2 = 4.6$ for the same number of parameters, $\ln B_{\mathrm{TT/NLA}} = -1.6$), because the quadratic term falls off faster with separation than the measured signal. Setting $b_{\mathrm{TA}} = 0$ instead of unity leaves this picture unchanged: the high-mass sample remains tangentially aligned at a similar significance, the marginal $A_2$ drops by a factor of 2–4, and the evidence favours NLA slightly more ($\ln B_{\mathrm{TATT/NLA}} = -1$). The ELG TATT fit is unconstrained: the $A_2$ posterior piles up against the lower edge of its prior ($A_2 = -6$), so Table B.1 quotes only an upper limit on $A_2$, and the $A_{\mathrm{IA}}$ posterior, which follows $A_2$ along the degeneracy, reflects the prior range rather than the data. Its fit is also poor ($\chi_\nu^2 = 2.21$; Appendix § 8).

Discussion

Halo-mass dependence of red-galaxy alignments

Fortuna et al. [42] found that a single power law with $\beta_M = 0.44 \pm0.04$ describes the halo-mass dependence of the alignment amplitude over $1.9 \times10^{12}$–$3.7 \times10^{14}\ h^{-1}M_\odot$, the range of the compilation they fitted, their own KiDS LRG samples spanning $2.7 \times10^{12}$–$2.6 \times10^{13}\ h^{-1}M_\odot$, and Herle et al. [54] reached the same conclusion in the FLAMINGO simulations above $10^{12}\ h^{-1}M_\odot$, with $\beta_M = 0.77 \pm0.015$. Above $10^{12}\ h^{-1}M_\odot$ our single-power-law slope, $\beta_M = 0.66 \pm0.18$, is consistent with both values, but the single power law is a poor fit, and the broken power law is strongly preferred (Sect. 5.1.3). Our measurements do not extend to the cluster masses ($> 10^{14}\ h^{-1}M_\odot$) that anchor the high-mass end of their sample. The flattening is not confined to masses below those probed by Fortuna et al. [42]: the break

and our BGS red bins at $12.4 \le\log(M_{\mathrm{eff}}/h^{-1}M_\odot) \le12.9$ lie within the range of their KiDS LRG samples, whose amplitudes there ($A_{\mathrm{IA}} = 1$–4, with large uncertainties) lie mostly above ours ($A_{\mathrm{IA}} = 0.9$–1.5). Within this common mass range, the difference may reflect the photometric versus spectroscopic selections or the different definitions of the red samples.

Halo mass may also not be the only variable at play: using the FLAMINGO simulations, Herle et al. [55] showed that the IA amplitude also depends on halo formation redshift, with earlier-forming haloes being more strongly aligned. If the faint red BGS galaxies preferentially occupied early-formed haloes, their alignment would be enhanced relative to the mass-only expectation and the $A_{\mathrm{IA}}$–$M_{\mathrm{eff}}$ relation would flatten at the faint end. In FLAMINGO, however, the formation time mainly shifts the amplitude while leaving its slope with mass unchanged, so this mechanism flattens the relation only if the mix of formation times changes with mass, early-forming haloes being increasingly represented at low mass. Testing the formation-time interpretation requires an observational proxy, such as the Sérsic index, the specific star formation rate or the stellar age (Herle et al. [55]); we leave it to future work.

Redshift dependence of LRG alignments

The high-redshift LRGs are more strongly aligned than the low-redshift ones at the same stellar mass (Sect. 5.1.2), and we have not identified the origin of this difference. Since the alignment amplitude of red galaxies increases with luminosity (Fortuna et al. [43], S26), a difference could arise if the high-redshift LRGs, selected from a flux-limited sample, were brighter at fixed stellar mass. They are, but only slightly: in common 0.1 dex stellar-mass slices, their FastSpecFit rest-frame $r$-band absolute magnitudes (k-corrected to $z = 0.1$) are brighter by 0.1–0.24 mag, as expected from passive evolution. With the bright-end luminosity slope $\beta= 1.2$ of Fortuna et al. [43], this raises the expected amplitude by at most 30%, far short of the factor of 2.5 between the two global amplitudes, so luminosity at fixed stellar mass does not explain the difference. The orientation-dependent target selection of DESI LRGs ([76]) and stellar-mass systematics at $z = 0.9$ could also contribute.

Tangential alignment of massive blue galaxies

The tangential alignment of the massive BGS blue galaxies is three times weaker than the radial alignment of the BGS red sample ($|A_{\mathrm{IA}}| = 0.6$ against $1.7$). The massive blue galaxies are as strongly star-forming as the low-mass ones: their median specific star formation rate from FastSpecFit, $\log(\mathrm{sSFR}/\mathrm{yr}^{-1}) = -11.07$, is the same as that of the low-mass bin ($-11.03$) and 2.7 dex above that of the BGS red sample ($-13.80$). In the linear alignment picture, a tangential alignment requires a negative linear response of the shapes to the tidal field, for which no standard mechanism exists: tidal stretching of pressure-supported galaxies yields a radial alignment. Tidal torquing, in which disc galaxies acquire their angular momentum from the tidal field and their shapes follow their spin, is the natural candidate. Within TATT, Blazek et al. [13] associated a tangential alignment of spiral galaxies with a negative $A_2$, the opposite of the sign we recover, but also showed that tidal torquing generates, on large scales, a contribution with the same scale dependence as the linear tidal alignment. A tangential alignment produced by torquing could therefore be absorbed by a negative $A_{\mathrm{IA}}$ in our fits, without any signature in $A_2$. Since TATT does not improve on NLA, and a torque-only model fits the signal less well than a linear one (Sect. 5.2.3), the data favour the scale dependence of a negative linear alignment, but our measurements cannot determine whether this linear term is itself the large-scale imprint of tidal torquing.

Hydrodynamical simulations do not offer a consensus either. In Horizon-AGN, disc (blue) galaxies are oriented tangentially around spheroidal (red) ones in three dimensions, an effect that is suppressed in projection (Chisari et al. [19]). Comparing TNG300, Horizon-AGN and EAGLE with common shape and morphology definitions, van Heukelum & Chisari [128] found that disc galaxies align radially or not at all; the tangential alignment survives only in Horizon-AGN, for kinematically selected discs around ellipticals at $z = 1$ with reduced shapes, which down-weight the galaxy outskirts. In IllustrisTNG (Pillepich et al. [103]), by contrast, the shapes of massive disc galaxies align with the tidal field (Shi et al. [115]), and blue and spiral galaxies show a radial linear alignment (Samuroff et al. [109]; Zjupa et al. [138]). A tidal-torquing alignment appears only for the most massive spirals at $z = 1$ in TNG100 (Zjupa et al. 2022), and between $z = 0$ and $z = 1$ in MillenniumTNG (Pakmor et al. [101]; Delgado et al. [27]). The sign of the alignment of disc galaxies in simulations thus depends on how the discs are selected, by colour or by kinematics, on how their shapes are measured, and on their stellar mass (van Heukelum & Chisari [128]). Our blue samples are selected by colour, and their shapes are measured by fitting a single Gaussian profile, which, like the reduced shapes, gives most weight to the inner parts of the galaxy, so a quantitative comparison requires simulated galaxies selected and measured in the same way, at the stellar masses of our high-mass bin. Whatever its origin, a non-zero alignment of massive blue galaxies confirms that restricting cosmic-shear sources to blue galaxies reduces, but does not remove, the IA contamination, so that such analyses must still marginalise over it (McCullough et al. [90]).

Comparison with S26

S26 measured IA with the same DESI DR1 tracers cross-matched with four imaging surveys (DES, KiDS, HSC, SDSS), using the projected statistic $w_{g+}$. For red galaxies the two analyses agree: both detect IA at $>10\sigma$ for the BGS red sample and in the LRGs. S26 described their amplitudes by luminosity and colour alone, with no residual stellar-mass dependence, but the dip of our amplitude at intermediate stellar mass is also visible in the stellar-mass panel of their fig:8. For blue galaxies, neither analysis finds a radial alignment ($A_{\mathrm{IA}} = 0.1 \pm0.3$ for their bluest BGS galaxies), and their bright blue BGS galaxies show a tangential alignment, $A_{\mathrm{IA}} = -0.4 \pm0.2$ (read from their figure), close to our high-mass value. The blue amplitudes cannot be compared one to one, since they depend strongly on where the colour boundary is drawn: relaxing the S26 cut from $M_r - M_z < 0.5$ to $0.65$ raises their amplitude to $1.4 \pm0.2$. For ELGs, S26 reported $A_{\mathrm{IA}} = -3.0 \pm1.0$ at $1.15 < z < 1.55$ and an amplitude consistent with zero below; our constraint over the combined range sits between the two. Their TATT fits, from $2\,h^{-1}\,\mathrm{Mpc}$ with $b_{\mathrm{TA}}$ and a second-order bias free, gave $A_2$ consistent with zero, which our $A_2$, not constrained independently of $A_{\mathrm{IA}}$, does not contradict.

Conclusions

We have measured the intrinsic alignments of DESI DR1 galaxies by cross-correlating the BGS, LRG, and ELG spectroscopic samples with the UNIONS weak-lensing shape catalogue, over the area common to the two surveys, 2400 deg$^2$ for BGS and 1200–1250 deg$^2$ for the LRG and ELG samples. We fitted the clustering monopole $\xi_{gg,0}$ and the galaxy–shape quadrupole $\xi_{g+}$ of the three-dimensional correlation function jointly with the NLA model. Each red sample was assigned an effective halo mass through a conditional stellar mass function fitted to the galaxy–galaxy lensing signal measured around the BGS red galaxies, so that the alignment amplitude could be followed as a function of halo mass rather than of a luminosity or stellar-mass proxy alone.

Red galaxies are strongly and radially aligned. In stellar-mass bins, the amplitude dips at intermediate masses and rises steeply at the massive end. As a function of effective halo mass, it is not described by a single power law over $6 \times10^{11}$–$3 \times10^{13}\,h^{-1}M_\odot$, even when the fit is restricted to $M_{\mathrm{eff}} > 10^{12}\,h^{-1}M_\odot$. A broken power law is strongly preferred, with a nearly flat faint-end below $M_b = 1.3 \times10^{13}\,h^{-1}M_\odot$ and a steep rise above it. The same break is recovered from the BGS red samples alone, for which the stellar-to-halo mass relation is measured rather than extrapolated. The flattening lies roughly within the mass range of the KiDS LRG samples of Fortuna et al. [42], so it is a tension with their single power law within the common mass range, possibly due to the different selections of the red samples.

Blue galaxies are not aligned in the same sense. The mild tangential alignment of the full BGS blue sample comes entirely from its high-mass half, which shows a tangential alignment at $3.3\sigma$ ($2.2\sigma$ without any model of its scale dependence), while the low-mass half is consistent with zero. We also fitted TATT to the blue samples. Although the tidal-alignment and tidal-torquing amplitudes are strongly degenerate, the tangential signal is carried by a negative tidal-alignment amplitude, while the tidal-torquing amplitude comes out slightly positive, i.e. of the sign that produces a radial alignment. A torque-only fit describes the data less well than NLA, and the Bayesian evidence does not favour TATT over NLA. The tangential alignment is therefore not produced by the quadratic tidal-torquing term of TATT, although tidal torquing could still contribute through the effective linear term it generates on large scales (Blazek et al. [13]). The ELG sample yields no detection and constrains the alignment amplitude to $-1.4 < A_{\rm IA} < 0.7$ at $2\sigma$ over $0.8 < z < 1.6$.

On identical data, the multipole estimators recover the amplitudes of $w_{g+}$ with uncertainties smaller by 20–43% (Appendix A.3), less than the factor of 1.5–2 reported by Si23, in part because of our more conservative scale cuts.

By reaching the faint end of the red population and blue galaxies down to $\log(M_\star/h^{-2}M_\odot) = 10$, these measurements extend direct IA constraints towards the populations that dominate cosmic-shear source samples. These results carry a direct implication for cosmological analyses. IA priors are usually built by extrapolating a luminosity- or mass-dependent power law calibrated on bright, massive samples down to the faint populations that dominate cosmic-shear source catalogues. Our measurements indicate that a single power law does not capture the alignment of low-mass red galaxies: below the break the amplitude is nearly flat, at $A_{\rm IA} = 1$–$2$, so that a single power law fitted across the break mispredicts it in both directions, too high just below the break and too low at the lowest masses. Mass-dependent IA models of the NLA-M or TATT-M family (Wright et al. [135]; Herle et al. [54]) should allow for a break rather than a single slope. That the relation flattens rather than continues to fall may also point to a second variable, such as halo formation time (Herle et al. [55]), which our data cannot isolate.

Three limitations should temper these conclusions. First, the halo-mass scale of the $A_{\rm IA}$–$M_{\rm eff}$ relation rests on the abundance anchor of the stellar-to-halo mass relation. That anchor depends on how the fibre-assignment incompleteness of BGS is corrected: with individual-inverse-probability weights, the effective halo masses shift coherently, while the shape of the relation and the break are unchanged (Appendix C). A fibre-complete BGS sample is therefore needed to pin the mass scale. Second, the galaxy bias is modelled as linear and constant, with no scale dependence and no evolution between the density and pair redshifts. This is safe for the clustering monopole, fitted above $30\,h^{-1}\,\mathrm{Mpc}$, but the shape quadrupole, in which $b_1$ and $A_{\rm IA}$ are degenerate, extends down to $10\,h^{-1}\,\mathrm{Mpc}$ ($6\,h^{-1}\,\mathrm{Mpc}$ for TATT), where a scale-dependent bias would be absorbed into the alignment amplitudes. Finally, NLA is an empirical description, so the amplitudes quoted here are effective amplitudes for the fitted scales. This matters most for the TATT fits of the blue samples, which rely on scales down to $6\,h^{-1}\,\mathrm{Mpc}$, where perturbative IA models describe the alignment signal poorly; the TATT amplitudes and the model comparison should therefore be regarded as indicative rather than established. DESI DR2, combined with the full UNIONS footprint, will enlarge both the overlapping area and the samples, raise the fibre-assignment completeness of BGS through its additional passes, and allow the measurements to be repeated on large scales alone: it will make it possible to confirm the tangential alignment of massive blue galaxies and to test whether the tidal-torquing signal survives on large scales, to split the ELG sample in redshift, and to calibrate the halo masses of the LRG lenses directly instead of by extrapolation.

Acknowledgements. RP, AC and FHP acknowledge financial support from the Centre national d’études spatiales (CNES), France (ROR: https://ror.org/04h1h0y33). MJH and LVW acknowledge support from NSERC through Discovery Grants. GSHP acknowledges support from the Pan-STARRS project. We are honoured and grateful for the opportunity of observing the Universe from Maunakea and Haleakala, which both have cultural, historical and natural significance in Hawai’i. 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-Hawai’i Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawai’i. 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 Hawai’i, and is supported by the NASA SSO Near Earth Observation Program under grants 80NSSC18K0971, NNX14AM74G, NNX12AR65G, NNX13AQ47G, NNX08AR22G, 80NSSC21K1572 and by the State of Hawai’i. This work was made possible by utilising the CANDIDE cluster at the Institut d’Astrophysique de Paris. The cluster was funded through grants from the PNCG, CNES, DIM-ACAV, the Euclid Consortium, and the Danish National Research Foundation Cosmic Dawn Center (DNRF140). It is maintained by Stéphane Rouberol. We gratefully acknowledge support from the CNRS/IN2P3 Computing Center (Lyon - France) for providing computing and data-processing resources needed for this work. The authors acknowledge the use of the Canadian Advanced Network for Astronomy Research (CANFAR) Science Platform operated by the Canadian Astronomy Data Centre (CADC) and the Digital Research Alliance of Canada (DRAC), with support from the National Research Council of Canada (NRC), the Canadian Space Agency (CSA), CANARIE, and the Canada Foundation for Innovation (CFI). WLKG thanks the University of Edinburgh School of Physics and Astronomy for a postdoctoral Fellowship. We also used NUMPY (Harris et al. [52]), SCIPY (Virtanen et al. [130]), MATPLOTLIB (Hunter [63]), TREECORR (Jarvis et al. [67]; Jarvis [66]), PYCORR with CORRFUNC (Sinha & Garrison [121]), HEALPY (Zonca et al. [139]), HEALPIX (Górski et al. [45]), SKYPROJ, ASTROPY (Astropy Collaboration et al. [5, 4, 3]), IMINUIT (Dembinski et al. [28]), and SCIKIT-LEARN (Pedregosa et al. [102]). The analysis is specified and executed with LIGHTCURVE-CN (https://github.com/LightconeResearch/lightcone-cli) following the ASTRA specification; the complete pipeline is available at https://github.com/rpaviot/unions-desi-ia-paper. The analysis pipeline, the reproducibility directory and the manuscript were developed and edited with the assistance of Claude (Anthropic), used through Claude Code; the author has verified all results and takes full responsibility for the content of this paper.

Intrinsic-alignment model and estimators

Alignment model

Following the pyccl implementation, the coefficients of the TATT expansion (Eq. (5)) are related to the dimensionless amplitudes quoted throughout this paper by

where $\bar{C}_1=5\times10^{-14}\,h^{-2}M_{\odot}^{-1}\,\mathrm{Mpc}^3$ ([14]), $\rho_{\mathrm{crit}}$ is the critical density at $z=0$, so that $\bar{C}_1\rho_{\mathrm{crit}}=0.0139$ is a dimensionless constant, and $D(z)$ is the linear growth factor normalised to unity at $z=0$; all the redshift and cosmology dependence is carried by $\Omega_{\mathrm{m}}/D(z)$.

The NLA model ([16, 62, 14]) retains only the linear term of the TATT expansion, with $C_1$ given by Eq. (A.1). The matter–intrinsic power spectrum in the NLA model, and the galaxy–intrinsic cross-power spectrum in real space that follows from linear galaxy bias, read

where $P_{\delta}$ is the non-linear matter power spectrum, which we compute with HMCode2020, using its CAMB implementation as wrapped in pyccl. The same non-linear $P_{\delta}$ enters the clustering spectrum $P_{\mathrm{gg}}$ (Eq. (A.3)) and the galaxy–intrinsic spectrum $P_{\mathrm{gI}}$. Galaxy shapes are spin-2 quantities; we restrict our analysis to $E$-modes throughout.

The higher-order terms of the TATT expansion (Eq. (5)) generate additional contributions to $P_{\mathrm{gI}}$ involving one-loop perturbation-theory integrals. We do not reproduce them here: the full expressions for the tidal-alignment, density-weighting and tidal-torquing contributions to the intrinsic-alignment power spectra are given in Blazek et al. [13], whose notation and conventions we follow, and we evaluate them with FAST-PT as implemented in pyccl.

Redshift-space distortions

Peculiar velocities along the line of sight induce anisotropies in the estimated three-dimensional separations, known as redshift-space distortions (RSDs). At linear order, the Kaiser effect ([72]) modifies the power spectra by factors of $(1+\beta\mu_k^2)$, where $\beta=f/b_1$ is the ratio of the linear growth rate $f$ to the linear galaxy bias $b_1$, and $\mu_k$ is the cosine of the angle between $\boldsymbol{k}$ and the line of sight. The galaxy clustering and galaxy–shape cross-power spectra in redshift space read (Si23)

where $P_{\mathrm{gI}}$ is the real-space galaxy–intrinsic power spectrum defined in Eq. (A.2). These angular-dependent factors give rise to quadrupole and hexadecapole contributions, which we model via the multipole formalism of Sect. 4.1. They enter the Hankel transforms of Eq. (6) through the coefficients $\alpha_{\ell,ab}$. For galaxy clustering ($\zeta_{ab}=0$), these read ([7, 120])

and for the galaxy–shape cross-correlation ($\zeta_{ab}=2$),

Consistency of the projected and multipole estimators

Our baseline constraints are derived from the multipoles $\tilde{\xi}_{\mathrm{gg},0}$ and $\tilde{\xi}_{\mathrm{g+},2}$, while most IA studies in the literature use the projected statistic $w_{\mathrm{g+}}$. Figure A.1 compares the posterior distributions of the galaxy bias $b_1$ and the IA amplitude $A_{\mathrm{IA}}$ obtained from the two estimators for the low-redshift LRG sample ($0.4<z<0.75$). The two posteriors are in excellent agreement, with $A_{\mathrm{IA}}=1.17\pm0.35$ from the multipoles and $A_{\mathrm{IA}}=1.18\pm0.44$ from $w_{\mathrm{g+}}$: the two estimators recover the same amplitude, and the multipole uncertainty is 20% smaller. The gain is larger for the blue samples, with $w_{\mathrm{g+}}$ uncertainties larger by factors of 1.43 to 1.78 (last column of Table 2), but remains below the factor of 1.5–2 reported by Si23. Part of this difference comes from the scale cuts. Si23 removed pairs with $r_p<5\,h^{-1}\,\mathrm{Mpc}$ (Sect. 4.1) and fitted both multipoles down to $5\,h^{-1}\,\mathrm{Mpc}$, whereas we fit $\tilde{\xi}_{\mathrm{gg},0}$ from 30 and $\tilde{\xi}_{\mathrm{g+},2}$ from 10 $h^{-1}\,\mathrm{Mpc}$ (Sect. 4.2), while our projected statistics use $r_{p,\min}=6\,h^{-1}\,\mathrm{Mpc}$ for both $w_{\mathrm{gg}}$ and $w_{\mathrm{g+}}$. The multipole fits thus rely on fewer and larger scales than the projected ones, so the improvement quoted here should be regarded as a lower bound on what the estimator can deliver.

Table B.1: Best-fit TATT model parameters for the blue samples from the multipole analysis. The density-weighting bias is fixed to $b_{\rm TA}=1$, leaving $b_1$, $A_{\rm IA}$ and $A_2$ free; the fits retain linear galaxy bias only.

Sample	$z_{\rm IA}$	$A_{\rm IA}$	$A_2$	$\chi_\nu^2$
BGS blue, global	0.25	$-0.85 \pm0.33$	$1.73 \pm0.75$	1.58
BGS blue, low $M_\star$	0.20	$-0.69 \pm0.36$	$1.64 \pm0.85$	1.53
BGS blue, high $M_\star$	0.30	$-1.46 \pm0.54$	$2.76 \pm1.16$	1.39
ELG, global	1.11	$-0.07 \pm0.87$	$< 0.96$	2.21

Notes. The sample properties ($z$ range, $\log M_\star$, $N_{\rm shapes}$, $\sigma_e$, $\mathcal{R}$) are listed in Table 2 and are not repeated here; the recovered $b_1$ agrees with the NLA value of Table 2 to better than 0.5% in every case. $A_{\rm IA}$ is the tidal-alignment amplitude (the same parameter as in the NLA fits of Table 2) and $A_2$ the tidal-torquing amplitude, normalised as in Eq. (A.1). The quadrupole scale cut is relaxed to $r_{\min}=6\,h^{-1}\,\mathrm{Mpc}$ for these fits, retaining all ten $\tilde{\xi}_{g+,2}$ points, while the monopole cut is unchanged at $30\,h^{-1}\,\mathrm{Mpc}$; these constraints therefore rely on smaller scales than the NLA baseline, so $\chi_\nu^2$ is not directly comparable with the NLA values of Table 2. The ELG $A_2$ posterior piles up against the lower edge of its prior ($-6 < A_2 < 6$), so we quote its 95% upper limit.

Fit results and correlation functions

This appendix gathers the best-fit parameters and the correlation-function measurements underlying the fits of Sect. 5 that are not shown in the main text.

The $\chi^2$ of the BGS blue multipole fits (Table 2) is driven by the clustering monopole that the three samples share ($\chi^2=15.3$ for eight points); the alignment quadrupole gives $\chi^2=11.2$, 6.1 and 10.0 for eight points in the low-mass, high-mass and global samples. For the ELG, the monopole and the quadrupole contribute about equally ($\chi^2=11.9$ and 12.8, eight points each). The ELG TATT fit ($\chi_\nu^2=2.21$) has no single outlier: its largest quadrupole residual ($2.5\sigma$) lies at $s=21\,h^{-1}\,\mathrm{Mpc}$, well inside the NLA range.

Sample	$z$ range	$z_{\mathrm{IA}}$	$\log M_\star$	$N_{\mathrm{shapes}}$	$\sigma_e$	$\mathcal{R}$	$b_1$	$A_{\mathrm{IA}}^{\tilde{\xi}_{g+,2}}(\chi^2)$
BGS red
Global	0.1–0.5	0.25	$10.73 \pm0.40$	583 617	0.236	0.80	$1.390 \pm0.040$	$1.663 \pm0.147$ (1.61)
$M_\star$ bin 1	0.1–0.5	0.17	$10.03 \pm0.13$	55 434	0.254	0.82	–	$2.211 \pm0.424$ (1.87)
$M_\star$ bin 2	0.1–0.5	0.20	$10.29 \pm0.05$	55 434	0.258	0.81	–	$1.764 \pm0.449$ (1.91)
$M_\star$ bin 3	0.1–0.5	0.21	$10.44 \pm0.04$	55 433	0.257	0.80	–	$1.597 \pm0.492$ (1.39)
$M_\star$ bin 4	0.1–0.5	0.23	$10.56 \pm0.03$	55 435	0.253	0.80	–	$2.070 \pm0.493$ (1.89)
$M_\star$ bin 5	0.1–0.5	0.24	$10.67 \pm0.03$	55 434	0.247	0.80	–	$1.113 \pm0.482$ (1.17)
$M_\star$ bin 6	0.1–0.5	0.25	$10.77 \pm0.03$	55 434	0.240	0.79	–	$0.857 \pm0.471$ (1.20)
$M_\star$ bin 7	0.1–0.5	0.27	$10.86 \pm0.03$	55 433	0.232	0.79	–	$1.474 \pm0.448$ (2.13)
$M_\star$ bin 8	0.1–0.5	0.28	$10.96 \pm0.03$	55 435	0.222	0.79	–	$1.321 \pm0.410$ (1.40)
$M_\star$ bin 9	0.1–0.5	0.29	$11.07 \pm0.03$	47 565	0.213	0.79	–	$1.076 \pm0.451$ (1.64)
$M_\star$ bin 10	0.1–0.5	0.30	$11.17 \pm0.03$	37 374	0.204	0.79	–	$0.219 \pm0.492$ (1.52)
$M_\star$ bin 11	0.1–0.5	0.32	$11.29 \pm0.04$	28 783	0.192	0.79	–	$2.872 \pm0.476$ (1.26)
$M_\star$ bin 12	0.1–0.5	0.33	$11.44 \pm0.08$	20 001	0.170	0.80	–	$5.724 \pm0.548$ (1.64)
LRG
Global (low-$z$)	0.4–0.75	0.59	$11.15 \pm0.16$	201 654	0.217	0.81	$1.952 \pm0.054$	$1.168 \pm0.353$ (1.45)
Global (high-$z$)	0.75–1.1	0.86	$11.13 \pm0.19$	208 369	0.213	0.74	$2.114 \pm0.064$	$2.945 \pm0.382$ (0.64)
$M_\star$ bin 1	0.4–0.75	0.60	$11.03 \pm0.05$	55 044	0.227	0.81	–	$1.006 \pm0.674$ (1.48)
$M_\star$ bin 2	0.4–0.75	0.60	$11.12 \pm0.02$	44 841	0.223	0.81	–	$-0.082 \pm0.579$ (1.51)
$M_\star$ bin 3	0.4–0.75	0.59	$11.19 \pm0.02$	38 050	0.219	0.82	–	$0.982 \pm0.652$ (1.55)
$M_\star$ bin 4	0.4–0.75	0.59	$11.27 \pm0.03$	31 500	0.210	0.81	–	$1.304 \pm0.795$ (1.90)
$M_\star$ bin 5	0.4–0.75	0.59	$11.41 \pm0.10$	30 000	0.191	0.80	–	$2.921 \pm0.722$ (1.65)
$M_\star$ bin 1	0.75–1.1	0.86	$10.95 \pm0.07$	58 721	0.245	0.69	–	$2.560 \pm0.731$ (0.56)
$M_\star$ bin 2	0.75–1.1	0.85	$11.09 \pm0.03$	47 709	0.213	0.75	–	$1.941 \pm0.688$ (0.74)
$M_\star$ bin 3	0.75–1.1	0.86	$11.18 \pm0.03$	38 548	0.204	0.76	–	$2.604 \pm0.737$ (0.67)
$M_\star$ bin 4	0.75–1.1	0.87	$11.28 \pm0.03$	31 098	0.197	0.77	–	$3.915 \pm0.938$ (0.96)
$M_\star$ bin 5	0.75–1.1	0.88	$11.42 \pm0.11$	30 000	0.195	0.77	–	$3.740 \pm0.808$ (0.83)

Notes. $\log M_\star$ is the median stellar mass in units of $\log(M/h^{-2}M_\odot)$, with the standard deviation of the distribution within each bin. $N_{\mathrm{shapes}}$ is the number of galaxies in the shape sample. $\sigma_e$ is the per-component calibrated shape noise, computed following the definition of Heymans et al. [60] as in Hervas Peters et al. [57]. $\mathcal{R}$ is the effective shear response, the weighted average of the diagonal elements of the per-galaxy Metacalibration response matrix, following Siegel et al. [116]. $z_{\mathrm{IA}}$ is the pair-weighted effective redshift of the density–shape cross-correlation, at which the galaxy–shape multipole is modelled; the clustering monopole is modelled at the effective redshift of the density sample, 0.35 (BGS red), 0.60 (LRG, low-$z$), 0.91 (LRG, high-$z$). The density sample of each tracer and redshift bin is its full clustering catalogue (for BGS, all BGS galaxies, shared by the red and blue samples), with which the global shape sample and its stellar-mass sub-samples are cross-correlated. Its linear bias $b_1$ is constrained by the clustering monopole, common to all sub-samples, and is therefore quoted only for the global samples. $A_{\mathrm{IA}}^{\tilde{\xi}_{g+,2}}$ is the IA amplitude from the multipole fit, with the reduced $\chi^2$ of the joint ($\tilde{\xi}_{gg,0},\tilde{\xi}_{g+,2}$) fit in parentheses.

Stellar-to-halo mass relation

To interpret the mass dependence of the IA amplitude, we estimate the stellar-to-halo mass relation (SHMR) of the BGS red sample from the galaxy–galaxy lensing signal, using the conditional stellar mass function (CSMF) halo-model formalism [136, 127, 15]; our modelling pipeline is similar to that of Dvornik et al. [33], and is implemented in the public NRVpy halo-model code⁸. In particular, we use the halo mass function and halo bias of Tinker et al. [125], NFW density profiles [97] with the concentration–mass relation of Duffy et al. [32], and halo masses defined within the radius enclosing 200 times the mean matter density ($M_{200\mathrm{m}}$), which is also the mass definition of the Dark Emulator [99] used for the non-linear halo bias correction described below. The model is evaluated at the effective redshift of each lens bin, in the cosmology of Sect. 4.3, with halo masses in $h^{-1}M_\odot$ and stellar masses in $h^{-2}M_\odot$.

Since the CSMF ties the galaxy abundance to halo mass, the lens samples must be complete in stellar mass. We therefore construct volume-limited samples: the limiting stellar mass of the red sample is estimated as a function of redshift with the automated turn-over method of Wright et al. [134] (see also Dvornik et al. [33]), applied here along the line of sight: in thin stellar-mass slices, the completeness limit is the comoving distance at which the number of galaxies per unit distance turns over, and the sample is split into ten stellar-mass bins, each restricted to the redshift range over which it is complete (Fig. C.1). The bin edges are chosen to yield comparable forecast $\Delta\Sigma$ signal-to-noise per bin, which narrows the bins at the high-mass end where the SHMR constraint resides.

We measure the excess surface mass density $\Delta\Sigma(r_{\mathrm{p}})$ around these ten lens bins (Sect. 3.3). The measurement uses the UNIONS shape catalogue [56] as the source sample. Unlike the IA shape sample, the sources are subject to the ShapePipe quality selection, $0.5 < T/T_{\mathrm{PSF}} < 3$ (ratio of the galaxy to the PSF size; this lower bound is less conservative than the 0.707 adopted for cosmic shear, Hervas-Peters et al. [56]), a SExtractor windowed $r$-band magnitude MAG_WIN $>$ 20, and $10 < \mathrm{S/N} < 200$, which mainly removes stars and spurious detections. We also require the sources to lie outside the bright-star and stellar-halo masks, a manual mask around large nearby galaxies (Messier and NGC objects), and the MAX-MASK mask of image defects and cosmic rays [100], and correct them for PSF leakage [47, 56]. They are further selected in photometric redshift, $0.70 < z_{\mathrm{B}} < 1.25$, using the BPZ [8] point estimate $z_{\mathrm{B}}$ from the UNIONS DR6 ugriz photometry, so that the sources lie well behind the $z < 0.5$ lenses. After masking, the source sample covers 4500 deg$^{2}$ with an effective density of 2.7 arcmin$^{-2}$. These photometric redshifts are used here for the first time, and only to select the sources: their individual values never enter the measurement, whose source redshift distribution is calibrated with clustering redshifts, as described below. The excess surface mass density is estimated as a weighted sum over lens–source pairs,

where the pair weight is $w_{\mathrm{ls}} = w_{\mathrm{l}} w_{\mathrm{s}} \Sigma_{\mathrm{crit,eff}}^{-2}(z_{\mathrm{l}})$, $e_{\mathrm{t,ls}}$ is the tangential ellipticity of the source relative to the lens⁹, $\mathcal{R}$ is the mean Metacalibration shear response, $w_{\mathrm{l}}$ is the lens weight ($w_{\mathrm{DESI}}$, i.e. Eq. (1) without $w_{\mathrm{FKP}}$), and $w_{\mathrm{s}} = w_{\mathrm{SP}}$ is the source shape weight (Sect. 2.1). The second term is the same weighted sum evaluated around the random-lens catalogue (index rs, with $w_{\mathrm{rs}}$ and $z_{\mathrm{r}}$ defined analogously). The subtraction of the random signal decreases the variance at large scales [118]. Rather than relying on individual source photometric redshifts, we compute an effective critical surface density for each lens redshift $z_{\mathrm{l}}$ as a weighted integral over the source redshift distribution, following Fortuna et al. [42]:

where $n(z_{\mathrm{s}})$ is the weighted redshift distribution of the sources. The latter is calibrated with clustering redshifts, obtained by cross-correlating the UNIONS sources with the DESI spectroscopic samples following the methodology of d’Assignies et al. (2025), using the DESI DR1 LRG and ELG LSS catalogues as reference samples: the cross-correlations are measured at $1.5 < r_{\mathrm{p}} < 5\,h^{-1}\,\mathrm{Mpc}$ in slices of $\Delta z = 0.1$ against the LRGs ($0.4 < z < 0.8$), the combined LRG and ELG sample ($0.8 < z < 1.1$) and the ELGs ($1.1 < z < 1.6$), and divided by the bias of each reference slice, obtained from its auto-correlation. The resulting distribution (Fig. C.2) has a mean redshift of 0.97, with a jackknife uncertainty of 0.04, and only 2% of its weight lies below $z = 0.6$. The covariance of $\Delta\Sigma$ is estimated by jackknife over 100 $K$-means regions of the lens footprint. We verified on the brightest lens bins that the boost factor is consistent with unity, as expected from the redshift separation between lenses and sources, and apply no boost or intrinsic-alignment correction; the lens magnification bias, estimated with the finite-difference method of Wenzl et al. [133], which we adapt to the DESI target selections by perturbing the total and fibre fluxes and the stellar masses of the lenses and re-applying the full sample selection, is negligible for these lenses and is not subtracted. Distances are computed in the same Planck cosmology as the IA fits (Sect. 4.3), and $\Delta\Sigma$ is expressed in $h\,M_{\odot}\,\mathrm{pc}^{-2}$ as a function of $r_{\mathrm{p}}$ in $h^{-1}\,\mathrm{Mpc}$.

In the model, the CSMF specifies the occupation of central and satellite galaxies as a function of halo mass, and the halo-model prediction for $\Delta\Sigma$ includes the baryonic contribution of the lens galaxies themselves, the non-linear halo bias correction $\beta^{\mathrm{NL}}$ (Mead & Verde [93]) computed with Dark Emulator, and baryonic feedback through a rescaling of the halo density-profile concentration (Eq. 15 of Dvornik et al. [33]). The free parameters are the four parameters of the central SHMR, the scatter of the central relation, the slope and the two normalisation parameters of the satellite CSMF, and the two profile rescaling amplitudes for matter and satellites. All parameters have flat priors, except the low-mass slope of the SHMR, $\gamma_1$, which is not constrained by our data (the lens bins lie above the knee of the relation) and for which we adopt a Gaussian prior $\gamma_1 \sim\mathcal{N}(7.10,\,2.0)$ centred on the value of Dvornik et al. [33]; the $\gamma_1$ of Table C.1 is therefore prior-dominated. The two profile rescaling amplitudes are bounded above by unity, i.e. the matter and satellite profiles cannot be more concentrated than the dark matter. In addition to the ten lensing profiles, the likelihood includes the measured comoving number density of each bin, which anchors the abundance side of the relation. The number density is measured on the lens sample itself, with the weighting that reproduces the DR1 $n(z)$ tables released with Adame et al. [1] and

dial alignment), whereas the lensing convention takes $e_{\mathrm{t}}$ positive for a tangential alignment around the lens.

Ross et al. [107]: the galaxies are counted with their completeness, redshift-failure and imaging-systematic weights (without the FKP weight), and the count is divided by the comoving volume of the bin over the effective footprint area, i.e. the area traced by the random catalogues weighted by the fraction of tiles observed at each position. We verified that this procedure applied to the DESI LRG sample reproduces the published DESI DR1 LRG $n(z)$. The uncertainty on the number density is estimated by jackknife over the same 100 regions as $\Delta\Sigma$.

The quantity minimised is

where, for each lens bin $i$, $\Delta\Sigma_i$ and $m_i$ are the measured and model lensing profiles over the fitted scales, $C_i$ the corresponding jackknife covariance, $n_{\mathrm{gal},i}$ and $\hat{n}_{\mathrm{gal},i}$ the measured and model number densities. The covariance between different lens bins is neglected. The fit is performed on scales $0.5 < r_{\mathrm{p}} < 60\,h^{-1}\,\mathrm{Mpc}$, i.e. ten $\Delta\Sigma$ points and one number density per bin, 110 data points for 10 free parameters. We sample the posterior with NAUTILUS, with 2000 live points and until the effective sample size reaches 10 000, and quote in Table C.1 the marginal posterior medians with their 16–84% intervals, together with the maximum a posteriori (MAP) sample of the chain, which we use as the best-fit model in what follows. The MAP model yields $\chi^2 = 142.2$ for 100 degrees of freedom: 116.4 from the lensing profiles (100 points, i.e. $\chi^2_{\Delta\Sigma} = 116.4$ for 90 degrees of freedom when the profiles are considered alone) and 25.8 from the number densities (10 points). A direct minimisation with MINUIT (James & Roos [65]) reaches nearly the same minimum ($\chi^2 = 143.2$, every parameter within $0.5\sigma$ of the MAP sample) but pushes both profile rescaling amplitudes to or close to their upper bound and returns parabolic errors that misrepresent the posterior, which is skewed along the $M_0$–$\gamma_2$ and $f_{\mathrm{h}}$–$f_{\mathrm{s}}$ degeneracies: the marginal medians sit about $1\sigma$ away from the mode, and $f_{\mathrm{h}} = 0.83^{+0.12}_{-0.17}$ and $f_{\mathrm{s}} = 0.87^{+0.10}_{-0.17}$ are constrained below unity.

In the MAP model, the satellite fraction of the lens bins decreases monotonically with stellar mass, from 26% in the lowest-mass bin to 2% in the most massive one. Evaluated over the stellar-mass ranges of the twelve BGS red alignment bins, it falls from 31% to below 1%, the last bin being an extrapolation beyond the most massive lens bin. Since $M_{\mathrm{eff}}$ is defined from the central occupation alone (Eq. (9)), it describes the haloes of the central galaxies; the satellites, which reside in more massive host haloes, make up a quarter to a third of the lowest-mass alignment bins and less than 10% above $\log(M_{\star}/h^{-2}M_{\odot}) \simeq11$.

We have also tested the fibre-assignment correction. Instead of the standard DESI completeness weights, we used individual-inverse-probability (IIP) weights, which up-weight each target that receives a fibre by the inverse of its fibre-assignment probability, which Lange et al. [80] showed to give an unbiased galaxy–galaxy lensing signal. We repeated the entire measurement and fit on the DESI DR1 LSS catalogues processed with these weights, keeping the lens binning, the source sample, the redshift distributions and the fitting configuration fixed, so that the completeness treatment is the only difference. The lensing profiles are unchanged within their uncertainties. The number densities, however, are lower, so the abundance anchor shifts the stellar-to-halo mass relation coherently towards higher halo

Parameter	Description	Median	MAP
$\log_{10} M_0$	SHMR normalisation	$10.72^{+0.07}_{-0.06}$	$10.65$
$\log_{10} M_1$	Characteristic halo mass	$11.18^{+0.11}_{-0.09}$	$11.11$
$\gamma_1$	SHMR low-mass slope	$9.1^{+1.4}_{-1.5}$	$9.6$
$\gamma_2$	SHMR high-mass slope	$0.188^{+0.020}_{-0.025}$	$0.208$
$\sigma_c$	Central SHMR scatter	$0.248^{+0.007}_{-0.008}$	$0.249$
$\alpha_s$	Satellite CSMF slope	$-1.21^{+0.11}_{-0.09}$	$-1.27$
$b_0^s$	Satellite norm. intercept	$-0.41^{+0.17}_{-0.17}$	$-0.55$
$b_1^s$	Satellite norm. slope	$1.10^{+0.10}_{-0.11}$	$1.18$
$f_h$	Matter profile rescaling	$0.83^{+0.12}_{-0.17}$	$0.92$
$f_s$	Satellite profile rescaling	$0.87^{+0.10}_{-0.17}$	$1.00$

Notes. $M_0$ is in $h^{-2}M_\odot$ and $M_1$ in $h^{-1}M_\odot$ ($M_{200\mathrm{m}}$). MAP: the maximum a posteriori sample of the chain, the best-fit model of Fig. C.4; $\chi^2 = 142.2$ there (Eq. (C.3), 100 $\Delta\Sigma$ points and 10 number densities for 10 free parameters).

masses. This moves the halo-mass scale of the $A_{\mathrm{IA}}$–$M_{\mathrm{eff}}$ relation of Sect. 5, but not its shape, and leaves our conclusions unchanged.

The resulting central SHMR is shown in Fig. C.3, together with the relation of Dvornik et al. [33] and the weak-lensing constraints compiled by Zu & Mandelbaum (2015). Over the stellar-mass range of our lens bins, $10.3 \le\log(M_\star/h^{-2}M_\odot) \le11.3$, our $\langle M_h|M_\star\rangle$ agrees with the relation of Zu & Mandelbaum (2015) to within 0.1 dex for $10.6 \le\log(M_\star/h^{-2}M_\odot) \le11.15$, and lies up to 0.2 dex above it at the low-mass end and 0.2 dex below it at the high-mass end. It agrees with the points of Velander et al. [129] to within 0.1 dex up to $\log(M_\star/h^{-2}M_\odot) = 11.0$ and lies 0.14 dex below their most massive point, while the early-type points of Mandelbaum et al. [88] lie 0.25–0.55 dex above our relation above $\log(M_\star/h^{-2}M_\odot) = 11$. The measured $\Delta\Sigma$ profiles and the MAP model are shown in Fig. C.4.