# Intrinsic correlations of galaxy sizes in a hydrodynamical cosmological simulation

Harry Johnston, Dana Sophia Westbeek, Sjoerd Weide, Nora Elisa Chisari, Yohan Dubois, Julien Devriendt, Christophe Pichon

► To cite this version:

Harry Johnston, Dana Sophia Westbeek, Sjoerd Weide, Nora Elisa Chisari, Yohan Dubois, et al. Intrinsic correlations of galaxy sizes in a hydrodynamical cosmological simulation. Monthly Notices of the Royal Astronomical Society, 2023, 520 (1), pp.1541-1566. 10.1093/mnras/stad201. hal-03808454

HAL Id: hal-03808454

https://hal.science/hal-03808454

Submitted on 16 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Intrinsic correlations of galaxy sizes in a hydrodynamical cosmological simulation

Harry Johnston,<fnref n="1" /> Dana Sophia Westbeek,<fnref n="1" /> Sjoerd Weide,<fnref n="1" /> Nora Elisa Chisari,<fnref n="1" />★ Yohan Dubois,<fnref n="2" />
Julien Devriendt<fnref n="3" /> and Christophe Pichon<fnref n="2" />,<fnref n="4" />

Accepted 2023 January 17. Received 2023 January 17; in original form 2022 September 22

## ABSTRACT

Residuals between measured galactic radii and those predicted by the Fundamental Plane (FP) are possible tracers of weak lensing magnification. However, observations have shown these to be systematically correlated with the large-scale structure. We use the Horizon-AGN hydrodynamical cosmological simulation to analyse these intrinsic size correlations (ISCs) for both elliptical (early-type) and spiral (late-type) galaxies at $z = 0.06$. We fit separate FPs to each sample, finding similarly distributed radius residuals, $\lambda$, in each case. We find persistent $\lambda\lambda$ correlations over three-dimensional separations $0.5–17 h^{-1}$ Mpc in the case of spiral galaxies, at $>3\sigma$ significance. When relaxing a mass-selection, applied for better agreement with galaxy clustering constraints, the spiral $\lambda\lambda$ detection strengthens to $9\sigma$; we detect a $5\sigma$ density-$\lambda$ correlation; and we observe intrinsically-large spirals to cluster more strongly than small spirals over scales $\lesssim 10 h^{-1}$ Mpc at $>5\sigma$ significance. Conversely, and in agreement with the literature, we observe lower-mass, intrinsically-small ellipticals to cluster more strongly than their large counterparts over scales $0.5–17 h^{-1}$ Mpc at $>5\sigma$ significance. We model $\lambda\lambda$ correlations using a phenomenological non-linear size model, and predict the level of contamination for cosmic convergence analyses. We find the systematic contribution to be of similar order to, or dominant over the cosmological signal. We make a mock measurement of an intrinsic, systematic contribution to the projected surface mass density $\Sigma(r)$, and find statistically significant low-amplitude, positive (negative) contributions from lower-mass spirals (ellipticals), which may be of concern for large-scale ($\gtrsim 7 h^{-1}$ Mpc) measurements.

**Key words:** Galaxies: elliptical and lenticular, cD–Galaxies: general–Galaxies: spiral–large-scale structure of Universe.

## 1 INTRODUCTION

Weak gravitational lensing is typically studied through estimation of the complex gravitational shear field $\gamma$, traced by the measurable ellipticities of source galaxies over a range of redshifts. Cosmic shear correlators defined with this evolving field, and the evolving galaxy density field, yield information describing the amount and clustering of matter in the late-time Universe (Hikage et al. 2019; Asgari et al. 2020; Heymans et al. 2021; Abbott et al. 2022; Secco et al. 2022).

The other component of weak lensing is the magnification sourced by the scalar convergence field $\kappa$, often probed via lensing of the cosmic microwave background (CMB; Aghanim et al. 2020; Fang et al. 2022). In photometric galaxy weak lensing surveys, the projected surface mass density $\Sigma \propto \kappa$ at a given redshift magnifies sources behind it. Effectively, it changes the observed flux coming from a distant source and the solid angle it subtends, while preserving its surface brightness. The increase in flux promotes galaxies across the survey detection threshold. Simultaneously, the background number density is diluted by the increase in solid angle. The relative strength of these competing effects is determined by

the faint-end slope of the survey luminosity function (Bartelmann & Schneider 2001); how many faint objects are ‘waiting’ to be promoted across the flux limit.

Magnification can induce galaxy number density fluctuations and correlations that must be correctly modelled in order to avoid catastrophic biases in cosmological parameter inference (Cardona et al. 2016; Hoekstra, Viola & Herbonnet 2017; Thiele, Duncan & Alonso 2020; Unruh et al. 2020; Duncan et al. 2022; Mahony et al. 2022). However, it also has potential as an aide to weak lensing halo-mass calibration (Rozo & Schmidt 2010; Hildebrandt et al. 2011; Umetsu et al. 2011; Hildebrandt et al. 2013; Ford et al. 2014; Duncan et al. 2016), and as a cosmological probe; various techniques have been developed to detect magnification correlations via measured galaxy sizes, magnitudes, redshifts, number densities, and even shears (Myers et al. 2005; Scranton et al. 2005; Hildebrandt, van Waaerbeke & Erben 2009; Ménard et al. 2010; Morrison et al. 2012; Schmidt et al. 2012; Alsing et al. 2015; Garcia-Fernandez et al. 2018; Liu et al. 2021).

Whilst recent work suggests that gains in the precision of cosmological parameter inference from the inclusion of magnification in cosmic shear analyses are modest when galaxy clustering is also included (Duncan et al. 2014; Lorenz, Alonso & Ferreira

<footnote n="1">¹Institute for Theoretical Physics, Utrecht University, Princetonplein, NL-3584 CC, Utrecht, The Netherlands</footnote>

<footnote n="2">²CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis, Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3">³Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="4">⁴IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, F-91191 Gif-sur-Yvette, France</footnote>

<footnote>*E-mail: n.e.chisari@uu.nl</footnote>
<figure><img src="image_2.png" /></figure>

# Intrinsic correlations of galaxy sizes in a hydrodynamical cosmological simulation

Harry Johnston,<fnref n="1" /> Dana Sophia Westbeek,<fnref n="1" /> Sjoerd Weide,<fnref n="1" /> Nora Elisa Chisari,<fnref n="1" />★ Yohan Dubois,<fnref n="2" />
Julien Devriendt<fnref n="3" /> and Christophe Pichon<fnref n="2" />,<fnref n="4" />

Accepted 2023 January 17. Received 2023 January 17; in original form 2022 September 22

## ABSTRACT

Residuals between measured galactic radii and those predicted by the Fundamental Plane (FP) are possible tracers of weak lensing magnification. However, observations have shown these to be systematically correlated with the large-scale structure. We use the Horizon-AGN hydrodynamical cosmological simulation to analyse these intrinsic size correlations (ISCs) for both elliptical (early-type) and spiral (late-type) galaxies at $z = 0.06$. We fit separate FPs to each sample, finding similarly distributed radius residuals, $\lambda$, in each case. We find persistent $\lambda\lambda$ correlations over three-dimensional separations $0.5–17 h^{-1}$ Mpc in the case of spiral galaxies, at $>3\sigma$ significance. When relaxing a mass-selection, applied for better agreement with galaxy clustering constraints, the spiral $\lambda\lambda$ detection strengthens to $9\sigma$; we detect a $5\sigma$ density-$\lambda$ correlation; and we observe intrinsically-large spirals to cluster more strongly than small spirals over scales $\lesssim 10 h^{-1}$ Mpc at $>5\sigma$ significance. Conversely, and in agreement with the literature, we observe lower-mass, intrinsically-small ellipticals to cluster more strongly than their large counterparts over scales $0.5–17 h^{-1}$ Mpc at $>5\sigma$ significance. We model $\lambda\lambda$ correlations using a phenomenological non-linear size model, and predict the level of contamination for cosmic convergence analyses. We find the systematic contribution to be of similar order to, or dominant over the cosmological signal. We make a mock measurement of an intrinsic, systematic contribution to the projected surface mass density $\Sigma(r)$, and find statistically significant low-amplitude, positive (negative) contributions from lower-mass spirals (ellipticals), which may be of concern for large-scale ($\gtrsim 7 h^{-1}$ Mpc) measurements.

**Key words:** Galaxies: elliptical and lenticular, cD–Galaxies: general–Galaxies: spiral–large-scale structure of Universe.

## 1 INTRODUCTION

Weak gravitational lensing is typically studied through estimation of the complex gravitational shear field $\gamma$, traced by the measurable ellipticities of source galaxies over a range of redshifts. Cosmic shear correlators defined with this evolving field, and the evolving galaxy density field, yield information describing the amount and clustering of matter in the late-time Universe (Hikage et al. 2019; Asgari et al. 2020; Heymans et al. 2021; Abbott et al. 2022; Secco et al. 2022).

The other component of weak lensing is the magnification sourced by the scalar convergence field $\kappa$, often probed via lensing of the cosmic microwave background (CMB; Aghanim et al. 2020; Fang et al. 2022). In photometric galaxy weak lensing surveys, the projected surface mass density $\Sigma \propto \kappa$ at a given redshift magnifies sources behind it. Effectively, it changes the observed flux coming from a distant source and the solid angle it subtends, while preserving its surface brightness. The increase in flux promotes galaxies across the survey detection threshold. Simultaneously, the background number density is diluted by the increase in solid angle. The relative strength of these competing effects is determined by

the faint-end slope of the survey luminosity function (Bartelmann & Schneider 2001); how many faint objects are ‘waiting’ to be promoted across the flux limit.

Magnification can induce galaxy number density fluctuations and correlations that must be correctly modelled in order to avoid catastrophic biases in cosmological parameter inference (Cardona et al. 2016; Hoekstra, Viola & Herbonnet 2017; Thiele, Duncan & Alonso 2020; Unruh et al. 2020; Duncan et al. 2022; Mahony et al. 2022). However, it also has potential as an aide to weak lensing halo-mass calibration (Rozo & Schmidt 2010; Hildebrandt et al. 2011; Umetsu et al. 2011; Hildebrandt et al. 2013; Ford et al. 2014; Duncan et al. 2016), and as a cosmological probe; various techniques have been developed to detect magnification correlations via measured galaxy sizes, magnitudes, redshifts, number densities, and even shears (Myers et al. 2005; Scranton et al. 2005; Hildebrandt, van Waaerbeke & Erben 2009; Ménard et al. 2010; Morrison et al. 2012; Schmidt et al. 2012; Alsing et al. 2015; Garcia-Fernandez et al. 2018; Liu et al. 2021).

Whilst recent work suggests that gains in the precision of cosmological parameter inference from the inclusion of magnification in cosmic shear analyses are modest when galaxy clustering is also included (Duncan et al. 2014; Lorenz, Alonso & Ferreira

<footnote n="1">¹Institute for Theoretical Physics, Utrecht University, Princetonplein, NL-3584 CC, Utrecht, The Netherlands</footnote>

<footnote n="2">²CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis, Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3">³Sub-department of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK</footnote>

<footnote n="4">⁴IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, F-91191 Gif-sur-Yvette, France</footnote>

<footnote>*E-mail: n.e.chisari@uu.nl</footnote>
2018; Duncan et al. 2022; Mahony et al. 2022), the potential for magnification analyses independent from cosmic shear systematics offers a valuable consistency test (Hildebrandt et al. 2009; Van Waerbeke 2010; Alsing et al. 2015; Ghosh, Durrer & Schäfer 2021), and future space surveys could enable magnification to approach the statistical power of cosmic shear (Casaponsa et al. 2013; Heavens, Alsing & Jaffe 2013).

Some works have explored the possibility of measuring weak lensing magnification by correlating the sizes of galaxies as characterized by residuals with respect to the Fundamental Plane (Bertin & Lombardi 2006; Huff & Graves 2014; Freudenburg, Huff & Hirata 2020). The Fundamental Plane (Djorgovski & Davis 1987) is a tight set of scaling relations, originating from the virial theorem. Assuming homologous mass-to-light ratios and constant surface brightness, one derives the Faber–Jackson relation between luminosity and velocity dispersion, and the Fundamental Plane (FP), relating the radius, velocity dispersion, and surface brightness of elliptical (early-type) galaxies. Breakdowns of these assumptions result in the well-studied ‘tilting’ of the FP (Djorgovski & Davis 1987; Bernardi et al. 2003; Hyde & Bernardi 2009; Saglia et al. 2010; Cappellari et al. 2013a,b; Saulder et al. 2013).

La Barbera et al. (2010) had previously shown that the best-fitting FP for a given sample of elliptical galaxies is sensitive to the local density contrast, and Joachimi, Singh & Mandelbaum (2015) and Singh, Yu & Seljak (2021) made detections of auto and cross-correlations between FP radius residuals and the large-scale structure for elliptical galaxies from SDSS DR8 (Saulder et al. 2013), and BOSS CMASS & LOWZ (Alam et al. 2015), respectively.

Intrinsic galaxy alignments – wherein tidal forces orient galaxies toward local density peaks in three dimensions – are thought to induce an orientation-dependent scatter in the galaxy size distribution, after shapes are projected onto the two-dimensional surface of the celestial sphere (Hirata 2009; Martens et al. 2018). However, Singh et al. (2020) showed that these were not sufficient to explain observed density-FP residual correlations, which must therefore feature some other physical or systematic contributions.

In analogy to the intrinsic alignments (IA; Catelan, Kamionkowski & Blandford 2001; Hirata & Seljak 2004) of galaxies as a contaminant to cosmic shear, intrinsic spatial correlations between galaxy sizes, or between sizes and the density field, could be mistakenly attributed to lensing magnification, and thus bias the cosmological interpretations of measured observables (Ciarlariello, Crittenden & Pace 2015; Ciarlariello & Crittenden 2016).

Moreover, intrinsic size correlations (ISCs) could further mimic magnification in promoting galaxies across detection thresholds. A hybrid of the lensing-induced size bias (Schmidt et al. 2009a), and the lensing-independent tidal alignment bias (Hirata 2009; Martens et al. 2018), an *intrinsic size bias* would be induced by aperture selections – a density-dependent selection effect with the potential to contaminate all measurable galaxy statistics (Schmidt et al. 2009b). To our knowledge, such an effect has not been explicitly studied in the literature, though developments in effective field theories offer possible avenues to do so (see e.g. Agarwal et al. 2021).

Besides the possibility of intrinsic contamination of lensing statistics derived from galaxy sizes, intrinsic galaxy size correlations are themselves of interest for astrophysics and cosmology. The scaling of galaxy sizes with their environments and other properties, as it relates to the divergent dynamics of early- and late-type objects, and their interplay with galaxy merger events, are all promising laboratories for studies of galaxy formation and evolution (Shen et al. 2003; Governato et al. 2007; Vale & Ostriker 2008; Naab, Johansson & Ostriker 2010; Oser et al. 2010; Newman et al. 2012; Dubois et al.

2013, 2016; Kravtsov 2013; Cappellari et al. 2013b; Welker et al. 2014, 2017).

This work uses data from the hydrodynamical Horizon-AGN simulation (Dubois et al. 2014) to explore the landscape of intrinsic size correlations with the benefit of precise determination of galaxy properties, and in the absence of gravitational lensing. We characterize intrinsic size distributions according to deviations from the fitted Fundamental Plane, similarly to Joachimi et al. (2015) and Singh et al. (2021), and we extend this concept to include an FP for spiral (late-type) galaxies, inspired by Shen, Mo & Shu (2002) who investigated spiral FPs motivated by the Tully–Fisher relation. This is a significant extension, as one can assume that spiral galaxies will dominate the deep samples utilized for studies of lensing magnification (Huff & Graves 2014).

We describe our simulated data in Section 2. Our Fundamental Planes and definitions of intrinsic sizes are detailed in Section 3. In Section 4, we outline our estimators for intrinsic size and galaxy density correlations in the simulation box, and our methods for predicting intrinsic contamination of magnification signals. Section 5 discusses our measured correlations, contamination predictions, and their implications, and our concluding remarks are presented in Section 6.

Throughout, we work with the $\Lambda$CDM cosmology of the Horizon-AGN simulation: $\Omega_m = 0.272$, $\Omega_b = 0.045$, $\Omega_\Lambda = 0.728$, $\sigma_8 = 0.81$, $H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_s = 0.967$. All quoted distances are comoving distances, though applications of a factor $h$ (the dimensionless Hubble parameter) can differ, and are specified by the quoted units.

## 2 SIMULATION

Horizon-AGN is a cosmological hydrodynamical simulation run with a standard flat $\Lambda$CDM cosmology {$\Omega_m = 0.272, \Omega_b = 0.045, \Omega_\Lambda = 0.728, \sigma_8 = 0.81, H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}, n_s = 0.967$} compatible with the constraints from WMAP7 (Komatsu et al. 2011). The box has a width of $L_{\text{box}} = 100 h^{-1} \text{ Mpc}$ and 1024$^3$ dark matter particles, each having a mass of $8 \times 10^7 \text{ M}_\odot$. Run with the adaptive mesh refinement code RAMSES (Teyssier 2002), and an initial gas mass resolution of $10^7 \text{ M}_\odot$, the mesh is adaptively refined down to $\Delta x = 1$ proper kpc. For more details on the simulation of gas cooling, star formation, and stellar and black hole feedback processes, see Dubois et al. (2014) and Dubois et al. (2016).

At each redshift snapshot, galaxies are identified in the simulation by running the ADAPTAHOP subclump finder (Aubert, Pichon & Colombi 2004) against the stellar particle distribution. Galactic stellar masses are computed as the sum over all ($\ge$50) stellar particles attributed to a given galaxy by ADAPTAHOP. Galaxy luminosities are computed using single stellar population models (Bruzual & Charlot 2003) with a Salpeter initial mass function. The flux contributed by each star depends on its mass, metallicity, and age. Total fluxes are convolved with Sloan Digital Sky Survey filters (Gunn et al. 2006) without dust extinction to obtain absolute AB magnitudes and rest-frame galaxy colours (e.g. $g - r$).

Horizon-AGN has been shown to reproduce a host of observed stellar and galactic observations: luminosity and stellar mass functions, the main sequence of star formation, the rest-frame colour distribution from UV to infrared, and the cosmic star formation history (Kaviraj et al. 2017). Dubois et al. (2016) showed the simulation to agree with observations of galaxy size-stellar mass relations from 3D-HST and CANDELS (Van Der Wel et al. 2014), and cited feedback from active galactic nuclei as important in achieving extended galactic profiles over cosmic evolution. Hatfield
2018; Duncan et al. 2022; Mahony et al. 2022), the potential for magnification analyses independent from cosmic shear systematics offers a valuable consistency test (Hildebrandt et al. 2009; Van Waerbeke 2010; Alsing et al. 2015; Ghosh, Durrer & Schäfer 2021), and future space surveys could enable magnification to approach the statistical power of cosmic shear (Casaponsa et al. 2013; Heavens, Alsing & Jaffe 2013).

Some works have explored the possibility of measuring weak lensing magnification by correlating the sizes of galaxies as characterized by residuals with respect to the Fundamental Plane (Bertin & Lombardi 2006; Huff & Graves 2014; Freudenburg, Huff & Hirata 2020). The Fundamental Plane (Djorgovski & Davis 1987) is a tight set of scaling relations, originating from the virial theorem. Assuming homologous mass-to-light ratios and constant surface brightness, one derives the Faber–Jackson relation between luminosity and velocity dispersion, and the Fundamental Plane (FP), relating the radius, velocity dispersion, and surface brightness of elliptical (early-type) galaxies. Breakdowns of these assumptions result in the well-studied ‘tilting’ of the FP (Djorgovski & Davis 1987; Bernardi et al. 2003; Hyde & Bernardi 2009; Saglia et al. 2010; Cappellari et al. 2013a,b; Saulder et al. 2013).

La Barbera et al. (2010) had previously shown that the best-fitting FP for a given sample of elliptical galaxies is sensitive to the local density contrast, and Joachimi, Singh & Mandelbaum (2015) and Singh, Yu & Seljak (2021) made detections of auto and cross-correlations between FP radius residuals and the large-scale structure for elliptical galaxies from SDSS DR8 (Saulder et al. 2013), and BOSS CMASS & LOWZ (Alam et al. 2015), respectively.

Intrinsic galaxy alignments – wherein tidal forces orient galaxies toward local density peaks in three dimensions – are thought to induce an orientation-dependent scatter in the galaxy size distribution, after shapes are projected onto the two-dimensional surface of the celestial sphere (Hirata 2009; Martens et al. 2018). However, Singh et al. (2020) showed that these were not sufficient to explain observed density-FP residual correlations, which must therefore feature some other physical or systematic contributions.

In analogy to the intrinsic alignments (IA; Catelan, Kamionkowski & Blandford 2001; Hirata & Seljak 2004) of galaxies as a contaminant to cosmic shear, intrinsic spatial correlations between galaxy sizes, or between sizes and the density field, could be mistakenly attributed to lensing magnification, and thus bias the cosmological interpretations of measured observables (Ciarlariello, Crittenden & Pace 2015; Ciarlariello & Crittenden 2016).

Moreover, intrinsic size correlations (ISCs) could further mimic magnification in promoting galaxies across detection thresholds. A hybrid of the lensing-induced size bias (Schmidt et al. 2009a), and the lensing-independent tidal alignment bias (Hirata 2009; Martens et al. 2018), an *intrinsic size bias* would be induced by aperture selections – a density-dependent selection effect with the potential to contaminate all measurable galaxy statistics (Schmidt et al. 2009b). To our knowledge, such an effect has not been explicitly studied in the literature, though developments in effective field theories offer possible avenues to do so (see e.g. Agarwal et al. 2021).

Besides the possibility of intrinsic contamination of lensing statistics derived from galaxy sizes, intrinsic galaxy size correlations are themselves of interest for astrophysics and cosmology. The scaling of galaxy sizes with their environments and other properties, as it relates to the divergent dynamics of early- and late-type objects, and their interplay with galaxy merger events, are all promising laboratories for studies of galaxy formation and evolution (Shen et al. 2003; Governato et al. 2007; Vale & Ostriker 2008; Naab, Johansson & Ostriker 2010; Oser et al. 2010; Newman et al. 2012; Dubois et al.

2013, 2016; Kravtsov 2013; Cappellari et al. 2013b; Welker et al. 2014, 2017).

This work uses data from the hydrodynamical Horizon-AGN simulation (Dubois et al. 2014) to explore the landscape of intrinsic size correlations with the benefit of precise determination of galaxy properties, and in the absence of gravitational lensing. We characterize intrinsic size distributions according to deviations from the fitted Fundamental Plane, similarly to Joachimi et al. (2015) and Singh et al. (2021), and we extend this concept to include an FP for spiral (late-type) galaxies, inspired by Shen, Mo & Shu (2002) who investigated spiral FPs motivated by the Tully–Fisher relation. This is a significant extension, as one can assume that spiral galaxies will dominate the deep samples utilized for studies of lensing magnification (Huff & Graves 2014).

We describe our simulated data in Section 2. Our Fundamental Planes and definitions of intrinsic sizes are detailed in Section 3. In Section 4, we outline our estimators for intrinsic size and galaxy density correlations in the simulation box, and our methods for predicting intrinsic contamination of magnification signals. Section 5 discusses our measured correlations, contamination predictions, and their implications, and our concluding remarks are presented in Section 6.

Throughout, we work with the ΛCDM cosmology of the Horizon-AGN simulation: $\Omega_m = 0.272$, $\Omega_b = 0.045$, $\Omega_\Lambda = 0.728$, $\sigma_8 = 0.81$, $H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $n_s = 0.967$. All quoted distances are comoving distances, though applications of a factor *h* (the dimensionless Hubble parameter) can differ, and are specified by the quoted units.

## 2 SIMULATION

Horizon-AGN is a cosmological hydrodynamical simulation run with a standard flat ΛCDM cosmology {$\Omega_m = 0.272, \Omega_b = 0.045, \Omega_\Lambda = 0.728, \sigma_8 = 0.81, H_0 = 70.4 \text{ km s}^{-1} \text{ Mpc}^{-1}, n_s = 0.967$} compatible with the constraints from WMAP7 (Komatsu et al. 2011). The box has a width of $L_{\text{box}} = 100 \text{ h}^{-1} \text{ Mpc}$ and 1024³ dark matter particles, each having a mass of $8 \times 10^7 \text{ M}_\odot$. Run with the adaptive mesh refinement code RAMSES (Teyssier 2002), and an initial gas mass resolution of $10^7 \text{ M}_\odot$, the mesh is adaptively refined down to $\Delta x = 1$ proper kpc. For more details on the simulation of gas cooling, star formation, and stellar and black hole feedback processes, see Dubois et al. (2014) and Dubois et al. (2016).

At each redshift snapshot, galaxies are identified in the simulation by running the ADAPTAHOP subclump finder (Aubert, Pichon & Colombi 2004) against the stellar particle distribution. Galactic stellar masses are computed as the sum over all (≥50) stellar particles attributed to a given galaxy by ADAPTAHOP. Galaxy luminosities are computed using single stellar population models (Bruzual & Charlot 2003) with a Salpeter initial mass function. The flux contributed by each star depends on its mass, metallicity, and age. Total fluxes are convolved with Sloan Digital Sky Survey filters (Gunn et al. 2006) without dust extinction to obtain absolute AB magnitudes and rest-frame galaxy colours (e.g. *g* − *r*).

Horizon-AGN has been shown to reproduce a host of observed stellar and galactic observations: luminosity and stellar mass functions, the main sequence of star formation, the rest-frame colour distribution from UV to infrared, and the cosmic star formation history (Kaviraj et al. 2017). Dubois et al. (2016) showed the simulation to agree with observations of galaxy size-stellar mass relations from 3D-HST and CANDELS (Van Der Wel et al. 2014), and cited feedback from active galactic nuclei as important in achieving extended galactic profiles over cosmic evolution. Hatfield
et al. (2019) measured galaxy clustering in Horizon-AGN light-cones designed to mimic VIDEO photometric observations (Jarvis et al. 2013), finding agreement for galaxies with stellar masses $M_* \gtrsim 10^{10.5} M_\odot$.

The successful reproduction of trends in galaxy size and clustering is crucial to our work here, as we seek to take advantage of precisely-known, lensing-free simulated quantities to gather meaningful predictions for intrinsic size correlations in the real Universe. As such, we consider galaxy samples both with and without a cut to the stellar mass $M_* > 10^{10.5} M_\odot$, and advise that the lower-mass results be taken with circumspection. Throughout this work, sample two-point correlations will feature the mass selection on (i) both, or (ii) neither of the samples.

For all galaxies, the object's semimajor *a* and semiminor *b* axes are computed in projection along the z-axis of the box (Chisari et al. 2015), and two-dimensional, anisotropic surface areas are defined as $\pi ab$.

The velocity dispersion per galaxy $\sigma = \frac{1}{3}\sqrt{\sigma_r^2 + \sigma_\theta^2 + \sigma_z^2}$ is computed via the radial $\sigma_r$, tangential $\sigma_\theta$, and vertical $\sigma_z$ dispersion components, defined with respect to the angular momentum vector of each galaxy, and the rotational velocity $V = \bar{v}_\theta$ is the average of stellar tangential velocity components $v_\theta$ (Dubois et al. 2016). A ‘circular’ velocity $V_c$ is also defined, as the rotational velocity at the virial radius $R_{vir*}$ – the radius within which the virial theorem is satisfied by the stellar particles assigned to the object.

Spiral galaxies are rotationally supported, having coherent stellar motions and low-stellar velocity dispersion, whilst pressure-supported ellipticals have high-stellar velocity dispersion due to the random motions of stars. We thus follow previous works in defining the boundary between elliptical and spiral galaxies according to the ratio of rotational velocity $V$ to velocity dispersion $\sigma$. We divide the sample at $V/\sigma = 0.6$ – though this choice is fairly arbitrary, it is roughly where Dubois et al. (2014) saw the ‘spin-flip’ occurring in Horizon-AGN, whereby the alignments of galaxies’ spin axes with nearby filamentary orientations transitions from (spiral/low-mass) parallel to (elliptical/high-mass) perpendicular (Welker et al. 2014; Codis et al. 2018; Bate et al. 2019; Veena et al. 2019; Kraljic et al. 2021; Lee et al. 2021b). As a marker of a morphological transition brought about by mergers, and having dramatic implications for correlations between galaxy spins/shapes and local tidal fields, this seems a sensible place to separate our putative intrinsic size correlations into spiral and elliptical contributions. We note, however, that a more detailed morphological classification would be of interest for studies seeking to bridge the gap between the elliptical and spiral FPs, and for studies of ISCs and IA in general.

A unified Fundamental Plane for elliptical and spiral galactic radii may prove elusive, given the divergent dynamics of such systems. However, Ferrero et al. (2021) show with hydrodynamical simulations and observations that the Tully–Fisher and Faber–Jackson scaling relations can be unified upon consideration of the ratio of stellar-to-dark halo mass enclosed within an effective (i.e. stellar half-mass) radius, going on to demonstrate a unified stellar mass plane. This, along with constraints upon galaxy stellar mass-size scaling relations (Kawinwanichakij et al. 2021; Nedkova et al. 2021; Rodriguez et al. 2021), and the identification of other tight scaling relations for spiral galaxies (Lagos et al. 2016; Matthee & Schaye 2018; Mancera Piña et al. 2021), should place a generalized FP for the radii of galaxies, agnostic of morphological type, within the realm of possibility.

If a unified plane were to be identified besides offering further insights into galaxy formation and evolution via FPs and ISCs, one might achieve a reduction in the number of parameters required to

model ISCs for magnification studies, and open up the possibility for size-based convergence analyses utilizing both morphological types, with colour-split FPs then offering a cross-checking mechanism. In this work, for simplicity, we define separate FPs for ellipticals and for spirals.

## 2.1 Elliptical Fundamental Plane properties

The effective radius $R_{eff}$ in this work is equal to the geometric mean of three half-mass radii, each computed after projection of an object's stellar particle distribution along a Cartesian axis of the simulation box (Dubois et al. 2016). Rosito et al. (2021) use a slightly different definition in their study of FPs in Horizon-AGN, where elliptical radii are taken as the three-dimensional radii containing half of each object's stellar mass – these estimates are none the less comparable, inhabiting a similar dynamic range.

The surface brightness $I$ is computed as the object luminosity (the sum of stellar particle luminosities) per unit area ($A = \pi ab$). We note that our combination of projected and three-dimensional quantities is likely to introduce some additional scatter and discrepancies with respect to the FPs of Rosito et al. (2021), who avoided projections. They also explored FPs replacing the surface brightness $I$ with the surface mass density $\Sigma$ – we shall refer to these planes as L-FP and $\Sigma$-FP, respectively, in forthcoming sections.

The final parameter of the elliptical FP is the stellar velocity dispersion $\sigma$ defined above.

## 2.2 Spiral Fundamental Plane properties

Shen et al. (2002) worked from the Tully–Fisher (TF) relation, which describes the positive scaling of spiral galaxies’ luminosity $L$ with the circular velocity $V_c$. Monte Carlo sampling their detailed disc dynamics model in search of a third variable for the FP, they found that the shape of the rotation curve and the disc scale-length $R_d$ were correlated with the scatter around the TF relation. We neglect to make specific estimations of $R_d$ for our simulated galaxies, finding the virial radius of stellar particles $R_{vir*}$ to be similarly correlated with the TF scatter. We thus make use of $R_{vir*}$ for our spiral FP, as a simple addition to the TF relation, and advise that follow-up work use estimates of more observationally tractable quantities such as the disc-scale length.

We note that the two radii under consideration have different dynamic ranges with a typical empirical ratio $R_{eff} \sim 0.1 R_{vir*}$. Both radii are measured using only the stellar particles assigned to an object by ADAPTAHOP, but the effective radius is an averaged half-mass radius (Section 2.1), whilst the virial radius is defined such that the motions of enclosed particles satisfy the virial theorem. We thus emphasize that the intrinsic sizes we are to define (Section 3.1) are probing inner-galactic radii for elliptical galaxies, ~3–50 $h^{-1}$ kpc, and outer-galactic radii, ~10–300 $h^{-1}$ kpc, for spirals.

Using abundance matching techniques and observational galaxy data, Kravtsov (2013) found galaxy half-mass and total (i.e. stellar plus dark matter) virial radii to scale almost linearly, with a power-law slope of ~0.95, across two decades of radius and eight decades of stellar mass – thus including all morphological types – and with a scatter of ~0.2 dex.

Our Horizon-AGN samples display similar scatters in the $R_{eff}-R_{vir*}$ relation, but shallower power-law slopes of ~0.5–0.6 for spiral samples, as well as on-average larger $R_{eff}$ and expectedly smaller $R_{vir*}$ (given that our virial radii are estimated only from the stellar particle distribution), resulting in a normalization ~10 × larger than seen by Kravtsov (2013). Concurrently, we find that $R_{vir*}$ enables reasonable
et al. (2019) measured galaxy clustering in Horizon-AGN light-cones designed to mimic VIDEO photometric observations (Jarvis et al. 2013), finding agreement for galaxies with stellar masses $M_* \gtrsim 10^{10.5} M_\odot$.

The successful reproduction of trends in galaxy size and clustering is crucial to our work here, as we seek to take advantage of precisely-known, lensing-free simulated quantities to gather meaningful predictions for intrinsic size correlations in the real Universe. As such, we consider galaxy samples both with and without a cut to the stellar mass $M_* > 10^{10.5} M_\odot$, and advise that the lower-mass results be taken with circumspection. Throughout this work, sample two-point correlations will feature the mass selection on (i) both, or (ii) neither of the samples.

For all galaxies, the object's semimajor *a* and semiminor *b* axes are computed in projection along the z-axis of the box (Chisari et al. 2015), and two-dimensional, anisotropic surface areas are defined as $\pi ab$.

The velocity dispersion per galaxy $\sigma = \frac{1}{3}\sqrt{\sigma_r^2 + \sigma_\theta^2 + \sigma_z^2}$ is computed via the radial $\sigma_r$, tangential $\sigma_\theta$, and vertical $\sigma_z$ dispersion components, defined with respect to the angular momentum vector of each galaxy, and the rotational velocity $V = \bar{v}_\theta$ is the average of stellar tangential velocity components $v_\theta$ (Dubois et al. 2016). A ‘circular’ velocity $V_c$ is also defined, as the rotational velocity at the virial radius $R_{vir*}$ – the radius within which the virial theorem is satisfied by the stellar particles assigned to the object.

Spiral galaxies are rotationally supported, having coherent stellar motions and low-stellar velocity dispersion, whilst pressure-supported ellipticals have high-stellar velocity dispersion due to the random motions of stars. We thus follow previous works in defining the boundary between elliptical and spiral galaxies according to the ratio of rotational velocity $V$ to velocity dispersion $\sigma$. We divide the sample at $V/\sigma = 0.6$ – though this choice is fairly arbitrary, it is roughly where Dubois et al. (2014) saw the ‘spin-flip’ occurring in Horizon-AGN, whereby the alignments of galaxies’ spin axes with nearby filamentary orientations transitions from (spiral/low-mass) parallel to (elliptical/high-mass) perpendicular (Welker et al. 2014; Codis et al. 2018; Bate et al. 2019; Veena et al. 2019; Kraljic et al. 2021; Lee et al. 2021b). As a marker of a morphological transition brought about by mergers, and having dramatic implications for correlations between galaxy spins/shapes and local tidal fields, this seems a sensible place to separate our putative intrinsic size correlations into spiral and elliptical contributions. We note, however, that a more detailed morphological classification would be of interest for studies seeking to bridge the gap between the elliptical and spiral FPs, and for studies of ISCs and IA in general.

A unified Fundamental Plane for elliptical and spiral galactic radii may prove elusive, given the divergent dynamics of such systems. However, Ferrero et al. (2021) show with hydrodynamical simulations and observations that the Tully–Fisher and Faber–Jackson scaling relations can be unified upon consideration of the ratio of stellar-to-dark halo mass enclosed within an effective (i.e. stellar half-mass) radius, going on to demonstrate a unified stellar mass plane. This, along with constraints upon galaxy stellar mass-size scaling relations (Kawinwanichakij et al. 2021; Nedkova et al. 2021; Rodriguez et al. 2021), and the identification of other tight scaling relations for spiral galaxies (Lagos et al. 2016; Matthee & Schaye 2018; Mancera Piña et al. 2021), should place a generalized FP for the radii of galaxies, agnostic of morphological type, within the realm of possibility.

If a unified plane were to be identified besides offering further insights into galaxy formation and evolution via FPs and ISCs, one might achieve a reduction in the number of parameters required to

model ISCs for magnification studies, and open up the possibility for size-based convergence analyses utilizing both morphological types, with colour-split FPs then offering a cross-checking mechanism. In this work, for simplicity, we define separate FPs for ellipticals and for spirals.

## 2.1 Elliptical Fundamental Plane properties

The effective radius $R_{eff}$ in this work is equal to the geometric mean of three half-mass radii, each computed after projection of an object's stellar particle distribution along a Cartesian axis of the simulation box (Dubois et al. 2016). Rosito et al. (2021) use a slightly different definition in their study of FPs in Horizon-AGN, where elliptical radii are taken as the three-dimensional radii containing half of each object's stellar mass – these estimates are none the less comparable, inhabiting a similar dynamic range.

The surface brightness *I* is computed as the object luminosity (the sum of stellar particle luminosities) per unit area ($A = \pi ab$). We note that our combination of projected and three-dimensional quantities is likely to introduce some additional scatter and discrepancies with respect to the FPs of Rosito et al. (2021), who avoided projections. They also explored FPs replacing the surface brightness *I* with the surface mass density $\Sigma$ – we shall refer to these planes as L-FP and $\Sigma$-FP, respectively, in forthcoming sections.

The final parameter of the elliptical FP is the stellar velocity dispersion $\sigma$ defined above.

## 2.2 Spiral Fundamental Plane properties

Shen et al. (2002) worked from the Tully–Fisher (TF) relation, which describes the positive scaling of spiral galaxies’ luminosity $L$ with the circular velocity $V_c$. Monte Carlo sampling their detailed disc dynamics model in search of a third variable for the FP, they found that the shape of the rotation curve and the disc scale-length $R_d$ were correlated with the scatter around the TF relation. We neglect to make specific estimations of $R_d$ for our simulated galaxies, finding the virial radius of stellar particles $R_{vir*}$ to be similarly correlated with the TF scatter. We thus make use of $R_{vir*}$ for our spiral FP, as a simple addition to the TF relation, and advise that follow-up work use estimates of more observationally tractable quantities such as the disc-scale length.

We note that the two radii under consideration have different dynamic ranges with a typical empirical ratio $R_{eff} \sim 0.1 R_{vir*}$. Both radii are measured using only the stellar particles assigned to an object by ADAPTAHOP, but the effective radius is an averaged half-mass radius (Section 2.1), whilst the virial radius is defined such that the motions of enclosed particles satisfy the virial theorem. We thus emphasize that the intrinsic sizes we are to define (Section 3.1) are probing inner-galactic radii for elliptical galaxies, ~3–50 $h^{-1}$ kpc, and outer-galactic radii, ~10–300 $h^{-1}$ kpc, for spirals.

Using abundance matching techniques and observational galaxy data, Kravtsov (2013) found galaxy half-mass and total (i.e. stellar plus dark matter) virial radii to scale almost linearly, with a power-law slope of ~0.95, across two decades of radius and eight decades of stellar mass – thus including all morphological types – and with a scatter of ~0.2 dex.

Our Horizon-AGN samples display similar scatters in the $R_{eff}-R_{vir*}$ relation, but shallower power-law slopes of ~0.5–0.6 for spiral samples, as well as on-average larger $R_{eff}$ and expectedly smaller $R_{vir*}$ (given that our virial radii are estimated only from the stellar particle distribution), resulting in a normalization ~10 × larger than seen by Kravtsov (2013). Concurrently, we find that $R_{vir*}$ enables reasonable
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Distributions of some spiral (top) and elliptical (bottom) galaxy properties in the Horizon-AGN $z = 0.06$ snapshot, selected according to a $V/\sigma = 0.6$ boundary. Solid-line (dashed-line) histograms denote Masscut (Full) samples to be fitted with Fundamental Planes, after application of all selections (shown here as grey shading). Black histograms display the Masscut/Full galaxy populations before any additional selections (see Section 2) on virial radii (**left**), effective radii (**middle-left**), spin magnitudes (**middle-right**), and stellar masses (**right**). The mass-selection works to remove small-radius spirals and ellipticals from the samples, affecting 70 and 40 per cent of selected objects, respectively.</figcaption></figure>

fits of a spiral Fundamental Plane (Section 3). This complicates the interpretation of our measured intrinsic size correlations, which we shall attempt to make clear in the coming sections.

The remaining parameters for the spiral FP are the circular velocity $V_c$ defined above, and the absolute r-band magnitude $M_r$.

## 2.3 Additional selections

In refining our FPs for spiral and elliptical galaxies, some objects were found to be poorly described. These included:

(i) Elliptical galaxies with effective radii $R_{\text{eff}} < 3 h^{-1} \text{kpc}$, approaching the resolution limit of the simulation;

(ii) Elliptical galaxies with stellar masses $M_* < 10^{10} M_{\odot}$;

(iii) Spiral galaxies with stellar virial radii $R_{\text{vir}} < 10 h^{-1} \text{kpc}$;

(iv) Galaxies with very large spin magnitudes (see Aubert et al. 2004, equation 2), which were found on inspection to be recently merged/merging structures that contaminate lower-mass samples.

We exclude these objects from our FP samples, as each contributes to heavily increased scatter, tilting, curving, or other irregular structures on the fitted planes. We also cut away the aforementioned lower-mass galaxies with $M_* < 10^{10.5} M_{\odot}$ (Hatfield et al. 2019) prior to fitting a more conservative set of FPs (see Section 2), which we denote as ‘Masscut FP’, as opposed to the ‘Full FP’ samples which include those lower-mass objects.

Fig. 1 shows the resulting spiral (top) and elliptical (bottom) galaxy property distributions, with solid-line histograms giving the Masscut samples, and dashed-line histograms the Full samples. Colours denote the galaxies that we fit with FPs after all selections (shown as grey shading) are applied, whilst black-lined histograms show the total Masscut/Full sample populations before any additional selections. One sees that the confluence of stellar mass, radius, and

<figure><img src="image_2.png" /><figcaption>**Figure 2.** Colour-magnitude diagrams for the full sample of spirals (blue) and ellipticals (red) in the Horizon-AGN $z = 0.06$ snapshot, selected according to a $V/\sigma = 0.6$ boundary. The colour-magnitude diagram clearly indicates the presence of a ‘blue cloud’ and a ‘red sequence’. The $V/\sigma$ cut serves to separate the populations efficiently, though some of the ellipticals still have significantly blue colours.</figcaption></figure>

spin selections serves to exclude populations of small, low-mass, fast-spinning galaxies, which are known to be overproduced in the simulation (Kaviraj et al. 2017).

Fig. 2 shows the colour-magnitude diagram for the two full samples. Spirals occupy the ‘blue cloud’ and ellipticals have mostly redder colours, clearly defining a ‘red sequence’. There is a significant elliptical fraction that overlaps with the blue cloud given our selection cuts. This is a known issue in cosmological hydrodynamical simulations. Despite colours largely matching observations, as demonstrated in Kaviraj et al. (2017), it is common to find bluer colours than in observations. Nevertheless, it is clear that despite some overlap, the $V/\sigma$ cut works to separate the two populations efficiently. This is also evidenced in terms of the quality of the fits
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Distributions of some spiral (top) and elliptical (bottom) galaxy properties in the Horizon-AGN $z = 0.06$ snapshot, selected according to a $V/\sigma = 0.6$ boundary. Solid-line (dashed-line) histograms denote Masscut (Full) samples to be fitted with Fundamental Planes, after application of all selections (shown here as grey shading). Black histograms display the Masscut/Full galaxy populations before any additional selections (see Section 2) on virial radii (**left**), effective radii (**middle-left**), spin magnitudes (**middle-right**), and stellar masses (**right**). The mass-selection works to remove small-radius spirals and ellipticals from the samples, affecting 70 and 40 per cent of selected objects, respectively.</figcaption></figure>

fits of a spiral Fundamental Plane (Section 3). This complicates the interpretation of our measured intrinsic size correlations, which we shall attempt to make clear in the coming sections.

The remaining parameters for the spiral FP are the circular velocity $V_c$ defined above, and the absolute r-band magnitude $M_r$.

## 2.3 Additional selections

In refining our FPs for spiral and elliptical galaxies, some objects were found to be poorly described. These included:

(i) Elliptical galaxies with effective radii $R_{\text{eff}} < 3 h^{-1} \text{kpc}$, approaching the resolution limit of the simulation;

(ii) Elliptical galaxies with stellar masses $M_* < 10^{10} M_{\odot}$;

(iii) Spiral galaxies with stellar virial radii $R_{\text{vir}} < 10 h^{-1} \text{kpc}$;

(iv) Galaxies with very large spin magnitudes (see Aubert et al. 2004, equation 2), which were found on inspection to be recently merged/merging structures that contaminate lower-mass samples.

We exclude these objects from our FP samples, as each contributes to heavily increased scatter, tilting, curving, or other irregular structures on the fitted planes. We also cut away the aforementioned lower-mass galaxies with $M_* < 10^{10.5} M_{\odot}$ (Hatfield et al. 2019) prior to fitting a more conservative set of FPs (see Section 2), which we denote as ‘Masscut FP’, as opposed to the ‘Full FP’ samples which include those lower-mass objects.

Fig. 1 shows the resulting spiral (top) and elliptical (bottom) galaxy property distributions, with solid-line histograms giving the Masscut samples, and dashed-line histograms the Full samples. Colours denote the galaxies that we fit with FPs after all selections (shown as grey shading) are applied, whilst black-lined histograms show the total Masscut/Full sample populations before any additional selections. One sees that the confluence of stellar mass, radius, and

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Colour-magnitude diagrams for the full sample of spirals (blue) and ellipticals (red) in the Horizon-AGN $z = 0.06$ snapshot, selected according to a $V/\sigma = 0.6$ boundary. The colour-magnitude diagram clearly indicates the presence of a ‘blue cloud’ and a ‘red sequence’. The $V/\sigma$ cut serves to separate the populations efficiently, though some of the ellipticals still have significantly blue colours.</figcaption></figure>

spin selections serves to exclude populations of small, low-mass, fast-spinning galaxies, which are known to be overproduced in the simulation (Kaviraj et al. 2017).

Fig. 2 shows the colour-magnitude diagram for the two full samples. Spirals occupy the ‘blue cloud’ and ellipticals have mostly redder colours, clearly defining a ‘red sequence’. There is a significant elliptical fraction that overlaps with the blue cloud given our selection cuts. This is a known issue in cosmological hydrodynamical simulations. Despite colours largely matching observations, as demonstrated in Kaviraj et al. (2017), it is common to find bluer colours than in observations. Nevertheless, it is clear that despite some overlap, the $V/\sigma$ cut works to separate the two populations efficiently. This is also evidenced in terms of the quality of the fits
of the fundamental relations of each sample, which we describe in Section 3.

## 3 FUNDAMENTAL PLANES

Rosito et al. (2021) studied the morphology, kinematics, and scaling relations of elliptical (early-type) galaxies in Horizon-AGN. They found tight FPs with no clear redshift evolution for $z \le 3$, and that energetic feedback from active galactic nuclei (omitted from the initially identical companion simulation, Horizon-noAGN) is necessary for the reproduction of several observations, including the Fundamental Plane. Our elliptical FPs are not directly comparable to those of Rosito et al. (2021), who estimate a three-dimensional half-mass radius (different to our $R_{eff}$, defined by Dubois et al. 2016) and use this to compute circularized surface areas. However, we do see a qualitatively similar plane, with equivalent signs in the coefficients, a comparable root-mean-square (rms) of residuals, and some minor asymmetry about the 1:1 relation (see Section 3.2). Given that Rosito et al. (2021) observed no redshift evolution of the FP, we limit our analysis here to a single simulation snapshot at $z = 0.06$, though we note that follow-up work incorporating higher-$z$ snapshots would be desirable.

We define our FPs similarly to Joachimi et al. (2015) and Rosito et al. (2021), but omit the redshift dependence employed by the former, since we consider only a single simulated redshift snapshot. We retain the surface brightness $I$ for our FP, in contrast to Rosito et al. (2021), who showed that the ‘L-FP’ (surface brightness FP) for their Horizon-AGN sample had a less symmetric scatter about the 1:1 relation, as compared with an FP using the surface stellar mass density $\Sigma$ ($\Sigma$-FP).

We find the opposite to be true for our samples, most likely as a consequence of the discrepant projected/three-dimensional quantities already discussed, though we also have minor differences in selections (e.g. Rosito et al. 2021 selected central galaxies according to numbers of stellar particles, whilst we make various cuts against stellar mass). Whilst the L-FP/$\Sigma$-FP difference is small at the level of the planes (hence we neglect to investigate in great detail), we shall see in Section 5.1 that the $\Sigma$-FP erases intrinsic size correlations as seen by the L-FP.

For elliptical galaxies, the FP is then given by

$$ \log R_{\text{eff}} = a \log \sigma + b \log I + c, \quad (1) $$

where $R_{eff}$ is the object’s effective radius, $\sigma$ is the stellar velocity dispersion, and $I$ is the surface brightness (for which Rosito et al. 2021 substituted surface stellar mass density $\Sigma$). Throughout this work, we shall denote the base-10 logarithm as ‘log’, and the natural logarithm as ‘ln’.

Following Shen et al. (2002), equation 13, we use a radius estimate to tighten the TF relation between luminosity $L$ (or absolute magnitude) and circular velocity $V_c$. Our spiral galaxy FP is thus given as

$$ \log R_{\text{vir*}} = \alpha M_r + \beta \log V_c + \gamma, \quad (2) $$

where $M_r$ is the simulated r-band absolute magnitude, and the circular velocity $V_c$ is estimated as the rotational velocity at the virial radius $R_{vir*}$, which we use in lieu of the disc scale length $R_d$ employed by Shen et al. (2002) (see Section 2.2).

We normalize each FP parameter to its median value before fitting the coefficients $a, b, c, \alpha, \beta, \gamma$ to each defined galaxy sample

via ordinary linear regression,<fnref n="1" /> and discuss the resulting FPs in Section 3.2.

### 3.1 Intrinsic sizes

Following Joachimi et al. (2015), we characterize the ‘intrinsic sizes’ of galaxies according to residuals between their measured radii and corresponding predictions from fitted Fundamental Planes. The dimensionless intrinsic size parameter $\lambda$ is given as

$$ \lambda \equiv \ln \left( \frac{R_I}{R_{FP}} \right) \approx \frac{R_I}{R_{FP}} - 1, \quad (3) $$

for the measured intrinsic radius $R_I$, which we take as $R_{eff}$ or $R_{vir*}$, for ellipticals and spirals, respectively, and the predicted radius $R_{FP}$ from the fitted FP. We will label the set of galaxies whose radii are predicted to be larger than they are in the simulation as $\lambda_+$, and those which are smaller as $\lambda_-$.

In this work, we are primarily interested in the intrinsic variability of galaxy sizes as a possible contaminant to the size fluctuations that one might attribute to the weak lensing convergence field $\kappa$. If we consider that lensing magnification operates on the intrinsic radius $R_I$ to produce the observed radius of a galaxy $R_O$, then we can explicitly define the latter as

$$ R_O = R_I(1 + \kappa) \approx R_{FP}(1 + \lambda)(1 + \kappa), \quad (4) $$

thus motivating our search for spatially-correlated, lensing-independent size fluctuations $\lambda$ that could bias size-based estimates of the convergence $\kappa$ (see Section 4).

### 3.2 Fundamental Plane results

Our fitted FPs (shown in Fig. 3) yield residuals $\lambda$ (inset axes) that are closely comparable with Gaussian distributions (orange curves), with means $\langle \lambda \rangle$ of $\mathcal{O}[10^{-10}]$ or smaller. Viewing the planes edge-on, we observe tilting in each, such that (small) large measured radii are (over-) underpredicted by the FPs; we characterize the degree of tilting by fitting a linear coefficient $m$ and a constant offset to each two-dimensional distribution $\log R_{FP} - \log R_{\text{vir*}/\text{eff}}$, showing the results in green in Fig. 3, for comparison with 1:1 relations, given as cyan lines.

The subsets of lower-mass objects in each Full FP are shown in Fig. 3 as coloured contours, encompassing 68 and 95 per cent of the supplemental objects, and labelled by the mass range of the subset: $\log(M_*/M_\odot) < 10.5$. The minimum elliptical stellar mass is $10^{10} M_\odot$, whilst spirals in the Full FP sample go down to $\sim 10^{8.5} M_\odot$ (see Fig. 1). The FP root-mean-square residuals $\sigma_{FP}$ are given as black and coloured numbers<fnref n="2" /> per-panel, corresponding to the Full FP and the low-mass subset, respectively.

Fitted FP coefficients are displayed in Table 1. The signs of the coefficients show that Masscut ellipticals with large radii have higher velocity dispersion and lower surface mass density. These trends are consistent with those seen for the Full FP sample, excepting a smaller intercept $c$, which is reflected in Figs 1 & 3 as the low-mass subset can be seen to inhabit smaller radii.

<footnote n="1">¹See e.g. Magoulas et al. (2012), Said et al. (2020), Howlett et al. (2022) for more complex, censored 3D Gaussian models for FP fitting, which are useful in the context of measurement errors. For our simulated, effectively noiseless quantities, we assume that linear regression will suffice.</footnote>

<footnote n="2">²These correspond to the rms of $\lambda$ divided by $\ln 10$, for comparison with previous work (Joachimi et al. 2015; Rosito et al. 2021; Singh et al. 2021).</footnote>
of the fundamental relations of each sample, which we describe in Section 3.

## 3 FUNDAMENTAL PLANES

Rosito et al. (2021) studied the morphology, kinematics, and scaling relations of elliptical (early-type) galaxies in Horizon-AGN. They found tight FPs with no clear redshift evolution for $z \le 3$, and that energetic feedback from active galactic nuclei (omitted from the initially identical companion simulation, Horizon-noAGN) is necessary for the reproduction of several observations, including the Fundamental Plane. Our elliptical FPs are not directly comparable to those of Rosito et al. (2021), who estimate a three-dimensional half-mass radius (different to our $R_{eff}$, defined by Dubois et al. 2016) and use this to compute circularized surface areas. However, we do see a qualitatively similar plane, with equivalent signs in the coefficients, a comparable root-mean-square (rms) of residuals, and some minor asymmetry about the 1:1 relation (see Section 3.2). Given that Rosito et al. (2021) observed no redshift evolution of the FP, we limit our analysis here to a single simulation snapshot at $z = 0.06$, though we note that follow-up work incorporating higher-$z$ snapshots would be desirable.

We define our FPs similarly to Joachimi et al. (2015) and Rosito et al. (2021), but omit the redshift dependence employed by the former, since we consider only a single simulated redshift snapshot. We retain the surface brightness $I$ for our FP, in contrast to Rosito et al. (2021), who showed that the ‘L-FP’ (surface brightness FP) for their Horizon-AGN sample had a less symmetric scatter about the 1:1 relation, as compared with an FP using the surface stellar mass density $\Sigma$ ($\Sigma$-FP).

We find the opposite to be true for our samples, most likely as a consequence of the discrepant projected/three-dimensional quantities already discussed, though we also have minor differences in selections (e.g. Rosito et al. 2021 selected central galaxies according to numbers of stellar particles, whilst we make various cuts against stellar mass). Whilst the L-FP/$\Sigma$-FP difference is small at the level of the planes (hence we neglect to investigate in great detail), we shall see in Section 5.1 that the $\Sigma$-FP erases intrinsic size correlations as seen by the L-FP.

For elliptical galaxies, the FP is then given by

$$ \log R_{\text{eff}} = a \log \sigma + b \log I + c, \quad (1) $$

where $R_{eff}$ is the object’s effective radius, $\sigma$ is the stellar velocity dispersion, and $I$ is the surface brightness (for which Rosito et al. 2021 substituted surface stellar mass density $\Sigma$). Throughout this work, we shall denote the base-10 logarithm as ‘log’, and the natural logarithm as ‘ln’.

Following Shen et al. (2002), equation 13, we use a radius estimate to tighten the TF relation between luminosity $L$ (or absolute magnitude) and circular velocity $V_c$. Our spiral galaxy FP is thus given as

$$ \log R_{\text{vir}*} = \alpha M_r + \beta \log V_c + \gamma, \quad (2) $$

where $M_r$ is the simulated r-band absolute magnitude, and the circular velocity $V_c$ is estimated as the rotational velocity at the virial radius $R_{vir*}$, which we use in lieu of the disc scale length $R_d$ employed by Shen et al. (2002) (see Section 2.2).

We normalize each FP parameter to its median value before fitting the coefficients $a, b, c, \alpha, \beta, \gamma$ to each defined galaxy sample

via ordinary linear regression,<fnref n="1" /> and discuss the resulting FPs in Section 3.2.

### 3.1 Intrinsic sizes

Following Joachimi et al. (2015), we characterize the ‘intrinsic sizes’ of galaxies according to residuals between their measured radii and corresponding predictions from fitted Fundamental Planes. The dimensionless intrinsic size parameter $\lambda$ is given as

$$ \lambda \equiv \ln \left( \frac{R_I}{R_{FP}} \right) \approx \frac{R_I}{R_{FP}} - 1, \quad (3) $$

for the measured intrinsic radius $R_I$, which we take as $R_{eff}$ or $R_{vir*}$, for ellipticals and spirals, respectively, and the predicted radius $R_{FP}$ from the fitted FP. We will label the set of galaxies whose radii are predicted to be larger than they are in the simulation as $\lambda_+$, and those which are smaller as $\lambda_-$.

In this work, we are primarily interested in the intrinsic variability of galaxy sizes as a possible contaminant to the size fluctuations that one might attribute to the weak lensing convergence field $\kappa$. If we consider that lensing magnification operates on the intrinsic radius $R_I$ to produce the observed radius of a galaxy $R_O$, then we can explicitly define the latter as

$$ R_O = R_I(1 + \kappa) \approx R_{FP}(1 + \lambda)(1 + \kappa), \quad (4) $$

thus motivating our search for spatially-correlated, lensing-independent size fluctuations $\lambda$ that could bias size-based estimates of the convergence $\kappa$ (see Section 4).

### 3.2 Fundamental Plane results

Our fitted FPs (shown in Fig. 3) yield residuals $\lambda$ (inset axes) that are closely comparable with Gaussian distributions (orange curves), with means $\langle \lambda \rangle$ of $\mathcal{O}[10^{-10}]$ or smaller. Viewing the planes edge-on, we observe tilting in each, such that (small) large measured radii are (over-) underpredicted by the FPs; we characterize the degree of tilting by fitting a linear coefficient $m$ and a constant offset to each two-dimensional distribution $\log R_{FP} - \log R_{vir*}/\text{eff}$, showing the results in green in Fig. 3, for comparison with 1:1 relations, given as cyan lines.

The subsets of lower-mass objects in each Full FP are shown in Fig. 3 as coloured contours, encompassing 68 and 95 per cent of the supplemental objects, and labelled by the mass range of the subset: $\log(M_*/M_\odot) < 10.5$. The minimum elliptical stellar mass is $10^{10} M_\odot$, whilst spirals in the Full FP sample go down to $\sim 10^{8.5} M_\odot$ (see Fig. 1). The FP root-mean-square residuals $\sigma_{FP}$ are given as black and coloured numbers<fnref n="2" /> per-panel, corresponding to the Full FP and the low-mass subset, respectively.

Fitted FP coefficients are displayed in Table 1. The signs of the coefficients show that Masscut ellipticals with large radii have higher velocity dispersion and lower surface mass density. These trends are consistent with those seen for the Full FP sample, excepting a smaller intercept $c$, which is reflected in Figs 1 & 3 as the low-mass subset can be seen to inhabit smaller radii.

<footnote n="1">¹See e.g. Magoulas et al. (2012), Said et al. (2020), Howlett et al. (2022) for more complex, censored 3D Gaussian models for FP fitting, which are useful in the context of measurement errors. For our simulated, effectively noiseless quantities, we assume that linear regression will suffice.</footnote>

<footnote n="2">²These correspond to the rms of $\lambda$ divided by $\ln 10$, for comparison with previous work (Joachimi et al. 2015; Rosito et al. 2021; Singh et al. 2021).</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Fundamental Planes fitted to elliptical (left) and spiral (right) galaxies selected from the Horizon-AGN simulation (Section 2). Hexagonally-binned two-dimensional histograms show the (log) $R_{FP}-R_{eff/vir}$ distributions for the Full Fundamental Planes (top), containing all galaxies that passed our selections, with the coloured contours marking 68 and 95 per cent of the lower-mass galaxies that are lost from the Masscut Fundamental Plane samples (bottom). The 1:1 relations are shown as solid cyan lines, and can be compared with linear fits to all points on each plane, shown as dashed green lines, and accompanied by the fitted gradients. Inset figures show the Fundamental Plane radius residuals $\lambda \equiv \ln(R_{FP}/R_{eff})$, overlain with Gaussian distributions (orange). The root-mean-square deviations $\sigma_{FP} = \text{rms}(\lambda)/\ln 10$ are given in black for the respective plane, and in colour for the subsets of lower-mass objects on the Full Fundamental Planes. Each of the fitted planes exhibits some degree of tilting, with spiral and Masscut planes most affected.</figcaption></figure>

<table><caption>**Table 1.** Sample details and Fundamental Plane constraints (if applicable) for elliptical, spiral (size tracers), and lens (density tracer) galaxy samples defined in the Horizon-AGN simulation. Columns give the sample, the number of galaxies, the ratio of the mean sample luminosity to a pivot luminosity corresponding to an absolute magnitude $M = -22$, the FP coefficients for elliptical (equation 1) and spiral (equation 2) FPs, and the root-mean-square deviation from the FP, respectively.</caption><thead><tr><th>FP/lens sample</th><th><i>N</i></th><th>⟨<i>L</i>⟩/<i>L</i><sub>piv</sub></th><th><i>a</i>|<i>α</i></th><th><i>b</i>|<i>β</i></th><th><i>c</i>|<i>γ</i></th><th>σ<sub>FP</sub></th></tr></thead><tbody><tr><td>Elliptical FP</td><td>6254</td><td>0.55</td><td>1.210 ± 0.093</td><td>-0.340 ± 0.055</td><td>0.729 ± 0.013</td><td>0.0737</td></tr><tr><td>Elliptical FP (&gt;10<sup>10.5</sup> M<sub>☉</sub>)</td><td>3684</td><td>0.84</td><td>1.041 ± 0.142</td><td>-0.380 ± 0.072</td><td>0.837 ± 0.017</td><td>0.0708</td></tr><tr><td>Spiral FP</td><td>26215</td><td>0.24</td><td>-0.217 ± 0.016</td><td>-0.587 ± 0.105</td><td>-1.488 ± 0.006</td><td>0.0866</td></tr><tr><td>Spiral FP (&gt;10<sup>10.5</sup> M<sub>☉</sub>)</td><td>6394</td><td>0.59</td><td>-0.286 ± 0.029</td><td>-1.289 ± 0.217</td><td>-1.295 ± 0.013</td><td>0.0877</td></tr><tr><td>Lens</td><td>7479</td><td>0.49</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Lens (&gt;10<sup>10.5</sup> M<sub>☉</sub>)</td><td>3741</td><td>0.84</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Our spiral FP coefficients and intercept exhibit agreeable signs with those found by Shen et al. (2002),<fnref n="3" /> who suggested negative values for each (if radii are taken in kiloparsecs). Both the

magnitude and circular velocity slopes are significantly steeper for the Masscut FP sample than for the Full FP, whilst the intercept is slightly less negative. These differences reflect the large space of $M_r-V_c$ opened-up by the addition of lower-mass galaxies.

We observe more of a tilt in the spiral plane, and a marginally larger rms scatter $\sigma_{FP}$, though the spiral FP is qualitatively similar

<footnote n="3">³Considering those of their FP definitions that are comparable with our own; their Eqs. 13 & 14, Tables 1–3.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Fundamental Planes fitted to elliptical (left) and spiral (right) galaxies selected from the Horizon-AGN simulation (Section 2). Hexagonally-binned two-dimensional histograms show the (log) $R_{FP}-R_{eff/vir}$ distributions for the Full Fundamental Planes (top), containing all galaxies that passed our selections, with the coloured contours marking 68 and 95 per cent of the lower-mass galaxies that are lost from the Masscut Fundamental Plane samples (bottom). The 1:1 relations are shown as solid cyan lines, and can be compared with linear fits to all points on each plane, shown as dashed green lines, and accompanied by the fitted gradients. Inset figures show the Fundamental Plane radius residuals $\lambda \equiv \ln(R_{FP}/R_{eff})$, overlain with Gaussian distributions (orange). The root-mean-square deviations $\sigma_{FP} = \text{rms}(\lambda)/\ln 10$ are given in black for the respective plane, and in colour for the subsets of lower-mass objects on the Full Fundamental Planes. Each of the fitted planes exhibits some degree of tilting, with spiral and Masscut planes most affected.</figcaption></figure>

<table><caption>**Table 1.** Sample details and Fundamental Plane constraints (if applicable) for elliptical, spiral (size tracers), and lens (density tracer) galaxy samples defined in the Horizon-AGN simulation. Columns give the sample, the number of galaxies, the ratio of the mean sample luminosity to a pivot luminosity corresponding to an absolute magnitude $M = -22$, the FP coefficients for elliptical (equation 1) and spiral (equation 2) FPs, and the root-mean-square deviation from the FP, respectively.</caption><thead><tr><th>FP/lens sample</th><th><i>N</i></th><th>⟨<i>L</i>⟩/<i>L</i><sub>piv</sub></th><th><i>a</i>|<i>α</i></th><th><i>b</i>|<i>β</i></th><th><i>c</i>|<i>γ</i></th><th>σ<sub>FP</sub></th></tr></thead><tbody><tr><td>Elliptical FP</td><td>6254</td><td>0.55</td><td>1.210 ± 0.093</td><td>-0.340 ± 0.055</td><td>0.729 ± 0.013</td><td>0.0737</td></tr><tr><td>Elliptical FP (&gt;10<sup>10.5</sup> M<sub>☉</sub>)</td><td>3684</td><td>0.84</td><td>1.041 ± 0.142</td><td>-0.380 ± 0.072</td><td>0.837 ± 0.017</td><td>0.0708</td></tr><tr><td>Spiral FP</td><td>26215</td><td>0.24</td><td>-0.217 ± 0.016</td><td>-0.587 ± 0.105</td><td>-1.488 ± 0.006</td><td>0.0866</td></tr><tr><td>Spiral FP (&gt;10<sup>10.5</sup> M<sub>☉</sub>)</td><td>6394</td><td>0.59</td><td>-0.286 ± 0.029</td><td>-1.289 ± 0.217</td><td>-1.295 ± 0.013</td><td>0.0877</td></tr><tr><td>Lens</td><td>7479</td><td>0.49</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Lens (&gt;10<sup>10.5</sup> M<sub>☉</sub>)</td><td>3741</td><td>0.84</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Our spiral FP coefficients and intercept exhibit agreeable signs with those found by Shen et al. (2002),<fnref n="3" /> who suggested negative values for each (if radii are taken in kiloparsecs). Both the

magnitude and circular velocity slopes are significantly steeper for the Masscut FP sample than for the Full FP, whilst the intercept is slightly less negative. These differences reflect the large space of $M_r-V_c$ opened-up by the addition of lower-mass galaxies.

We observe more of a tilt in the spiral plane, and a marginally larger rms scatter $\sigma_{FP}$, though the spiral FP is qualitatively similar

<footnote n="3">³Considering those of their FP definitions that are comparable with our own; their Eqs. 13 & 14, Tables 1–3.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 4.** The intrinsic sizes $\lambda$ of spiral (blues) and elliptical (reds) galaxies, as estimated from Fundamental Plane residuals (Section 3), versus various galaxy properties (see Section 2 for details on these and other properties). $\lambda_{\text{spin}}$ is the spin parameter (Aubert et al. 2004, equation 2), not to be confused with Fundamental Plane residuals $\lambda$. Shown are the relations both for conservative Masscut Fundamental Plane samples ($M_* > 10^{10.5} M_{\odot}$; darker colours), and for Full Fundamental Plane samples (lighter colours; see Table 1 for sample details). Points and error bars show the means and 16$^{th}$, 84$^{th}$ percentiles of $\lambda$ for six equipopulated bins defined on each of the x-axes. Several galaxy property variables exhibit correlations with the Fundamental Plane residuals $\lambda$, which must be taken into account when considering their spatial correlations. Panel headings indicate in which cases the quantity in the x-axis is part of the Fundamental Plane fits. Deviations from the horizontal in this case indicate a residual dependency on the parameter that could not be captured by the best-fit Fundamental Plane. For panels without headings, deviations from the horizontal suggest that there is a correlation between sizes and the parameter considered, and that including it as a variable in the Fundamental Plane fit could improve the predictability of galaxy sizes.

to that of ellipticals. We note also that the Full FP is less tilted than the Masscut FP for both ellipticals and spirals. These findings are promising for FP studies wishing to explore the evolution and statistics of deep, spiral-dominated galaxy samples.

A comparison with the previous study of the elliptical galaxy FP in Horizon-AGN (Rosito et al. 2021) is not direct; as discussed in Section 2, that analysis used differently defined three-dimensional/projected radii and surface areas. These, as well as minor selection differences, conspire to yield incompatible coefficients in the fit – though they deviate in the same directions from the virial-theorem-constant-$M/L$ prediction ($a=2, b=-1$), and the rms residuals are comparable (~0.071–0.074 in this work, versus 0.067 in Rosito et al. 2021).

Fig. 4 displays relationships between our Masscut and Full FP residuals $\lambda$ and galaxy properties. Elliptical $\lambda$ exhibit a weakly positive correlation with the effective radius, as expected given the slight tilting of the plane. The softly positive trend of $\lambda$ with measured radii is similar to that seen by Singh et al. (2021), who also saw strong negative correlations between $\lambda$ and absolute magnitude, which we do not observe. We observe slightly stronger positive correlations of spiral $\lambda$ with virial (due to tilting) and effective radii, and with stellar mass. Correlations between $\lambda$ and other FP variables are weak,

barring some correlations of spiral $\lambda$ with projected surface quantities – none of which feature on the spiral FP – and an anticorrelation with the morphological metric $V/\sigma$. Thus faster-spinning spirals tend to be intrinsically smaller.

The Masscut FPs are qualitatively similar to the Full FPs, but featuring slightly stronger tilts in both cases (though within $1\sigma$ of the gradients quoted in Fig. 3), and losing significant numbers of low-radius objects. Differences between $\lambda$ estimates with and without the mass selection are uncorrelated with stellar mass.

We note that the mass-selection tends to exacerbate correlations between FP residuals and other galaxy properties (dark versus light points in Fig. 4) in all cases – thus Full FP residuals are less contaminated by object characteristics other than intrinsic size fluctuations. We take these various correlations into consideration when interpreting measured ISCs in the coming sections.

Our Horizon-AGN samples thus yield good fits of the FP, with elliptical coefficients and rms scatter that are comparable with those found for observational data (Hyde & Bernardi 2009; Saulder et al. 2013; Joachimi et al. 2015; Saulder et al. 2019; Singh et al. 2020) and hydrodynamical simulations (Horizon-AGN; Rosito et al. 2021, IllustrisTNG; Lu et al. 2020). For simulated spiral galaxies, the planes are comparable with the theoretical planes of Shen et al. (2002),
<figure><img src="image_6.png" /></figure>

**Figure 4.** The intrinsic sizes $\lambda$ of spiral (blues) and elliptical (reds) galaxies, as estimated from Fundamental Plane residuals (Section 3), versus various galaxy properties (see Section 2 for details on these and other properties). $\lambda_{\text{spin}}$ is the spin parameter (Aubert et al. 2004, equation 2), not to be confused with Fundamental Plane residuals $\lambda$. Shown are the relations both for conservative Masscut Fundamental Plane samples ($M_* > 10^{10.5} M_{\odot}$; darker colours), and for Full Fundamental Plane samples (lighter colours; see Table 1 for sample details). Points and error bars show the means and 16$^{th}$, 84$^{th}$ percentiles of $\lambda$ for six equipopulated bins defined on each of the x-axes. Several galaxy property variables exhibit correlations with the Fundamental Plane residuals $\lambda$, which must be taken into account when considering their spatial correlations. Panel headings indicate in which cases the quantity in the x-axis is part of the Fundamental Plane fits. Deviations from the horizontal in this case indicate a residual dependency on the parameter that could not be captured by the best-fit Fundamental Plane. For panels without headings, deviations from the horizontal suggest that there is a correlation between sizes and the parameter considered, and that including it as a variable in the Fundamental Plane fit could improve the predictability of galaxy sizes.

to that of ellipticals. We note also that the Full FP is less tilted than the Masscut FP for both ellipticals and spirals. These findings are promising for FP studies wishing to explore the evolution and statistics of deep, spiral-dominated galaxy samples.

A comparison with the previous study of the elliptical galaxy FP in Horizon-AGN (Rosito et al. 2021) is not direct; as discussed in Section 2, that analysis used differently defined three-dimensional/projected radii and surface areas. These, as well as minor selection differences, conspire to yield incompatible coefficients in the fit – though they deviate in the same directions from the virial-theorem-constant-$M/L$ prediction ($a=2, b=-1$), and the rms residuals are comparable (~0.071–0.074 in this work, versus 0.067 in Rosito et al. 2021).

Fig. 4 displays relationships between our Masscut and Full FP residuals $\lambda$ and galaxy properties. Elliptical $\lambda$ exhibit a weakly positive correlation with the effective radius, as expected given the slight tilting of the plane. The softly positive trend of $\lambda$ with measured radii is similar to that seen by Singh et al. (2021), who also saw strong negative correlations between $\lambda$ and absolute magnitude, which we do not observe. We observe slightly stronger positive correlations of spiral $\lambda$ with virial (due to tilting) and effective radii, and with stellar mass. Correlations between $\lambda$ and other FP variables are weak,

barring some correlations of spiral $\lambda$ with projected surface quantities – none of which feature on the spiral FP – and an anticorrelation with the morphological metric $V/\sigma$. Thus faster-spinning spirals tend to be intrinsically smaller.

The Masscut FPs are qualitatively similar to the Full FPs, but featuring slightly stronger tilts in both cases (though within $1\sigma$ of the gradients quoted in Fig. 3), and losing significant numbers of low-radius objects. Differences between $\lambda$ estimates with and without the mass selection are uncorrelated with stellar mass.

We note that the mass-selection tends to exacerbate correlations between FP residuals and other galaxy properties (dark versus light points in Fig. 4) in all cases – thus Full FP residuals are less contaminated by object characteristics other than intrinsic size fluctuations. We take these various correlations into consideration when interpreting measured ISCs in the coming sections.

Our Horizon-AGN samples thus yield good fits of the FP, with elliptical coefficients and rms scatter that are comparable with those found for observational data (Hyde & Bernardi 2009; Saulder et al. 2013; Joachimi et al. 2015; Saulder et al. 2019; Singh et al. 2020) and hydrodynamical simulations (Horizon-AGN; Rosito et al. 2021, IllustrisTNG; Lu et al. 2020). For simulated spiral galaxies, the planes are comparable with the theoretical planes of Shen et al. (2002),
and feature stronger tilts than the elliptical planes. Still, correlations between FP residuals and key properties, such as the virial radius and stellar mass, are not drastic – particularly in the case of the Full FP. We thus advance to the measurement of intrinsic size correlations so-derived, applying caution as regards the impacts of tilting, λ-property correlations, and of our usage of virial radii to characterize spiral sizes.

## 4 INTRINSIC SIZE CORRELATIONS

We follow Joachimi et al. (2015) and Singh et al. (2021) in measuring galaxy clustering and intrinsic size correlations, though we do not integrate signals along the line-of-sight (the simulation box z-axis in our work), preferring to leverage our exact knowledge of object positions to assess correlations as functions of three-dimensional pair separations. The relevant correlation functions are estimated via

$$ \hat{\xi}_{gg}(r) = \frac{DD - 2DR_D + R_D R_D}{R_D R_D}, \qquad (5) $$

$$ \hat{\xi}_{g\lambda}(r) = \frac{D\lambda - R_D\lambda}{R_\lambda R_D}, \qquad (6) $$

$$ \hat{\xi}_{\lambda_1\lambda_2}(r) = \frac{\lambda_1\lambda_2}{R_{\lambda_1}R_{\lambda_2}}, \qquad (7) $$

where DD, RD, RR are unweighted pair-counts between galaxy D and random point R samples (with subscripts denoting galaxy samples to which the randoms correspond), for measurements of galaxy clustering. These, and the weighted pair-counts Dλ, R_Dλ, λ_1λ_2, are given by

$$ XY = \sum_{i,j} \Theta_{ij}, \qquad (8) $$

$$ X\lambda = \sum_{i,j} \lambda_i \Theta_{ij}, \qquad (9) $$

$$ \lambda_1\lambda_2 = \sum_{i,j} \lambda_{1,i}\lambda_{2,j} \Theta_{ij}, \qquad (10) $$

where samples X, Y denote the galaxy D or randoms R, R_D, R_λ samples, i, j are indices running over all objects in the two samples, and Θ_{ij} are binary filters applying a log-spaced binning in three-dimensional separation r. Notice equation (5) is the usual Landy-Szalay estimator (Landy & Szalay 1993). All estimators presented in equations (5), (6), and (7) are normalized by random-random pair counts (versus data-data and data-random pair counts). This is to allow direct comparison of our results to previous works (Joachimi et al. 2015; Singh et al. 2021), and to ensure that they can be interpreted and modelled as correlation functions in the sections that follow.

We measure the autoclustering correlations $\hat{\xi}_{gg}$ for each of our defined FP samples (Table 1), and also split those samples into intrinsically large λ_+ and small λ_- subsets for assessment of the clustering variation across the FP.

We define two additional ‘lens’ samples of elliptical galaxies, with masses >10^10 and >10^10.5 M_⊙ – the latter, more conservative mass selection matching that applied to our Masscut FP samples (see Section 2) – for use as density tracers in the density-size cross-correlations $\hat{\xi}_{g\lambda}$. The lens samples are similar to the elliptical FP samples, but retain galaxies lost to poor fitting of the FPs, thus offering slight reductions in shot-noise. Intrinsic size autocorrelations $\hat{\xi}_{\lambda\lambda}$ are measured only within FP samples, since λ probes different radii between ellipticals and spirals (Section 2), and is derived from differently fitted and tilted planes between the Masscut and Full set-ups.

We note that, whilst our fitted λ distributions are extremely close to symmetric by construction, any $\langle \lambda \rangle \neq 0$ will result in scale-invariant additive contributions to λ-weighted correlation functions, most visible at large separations where correlations are weak, and stemming from inherent clustering contributions (Singh et al. 2021). We ensure that this contribution is compensated by subtracting $\langle \lambda \rangle$ from each λ distribution under consideration, though in practice, for these samples, this subtraction is small and not necessary for the measurement of stable signals in the simulation box.

We work with seven three-dimensional separation bins in the range 0.5–30 h⁻¹ Mpc, and estimate covariances according to a delete-one jackknife resampling of 64 subvolumes defined on the simulation box. The covariance estimates are given by

$$ \widehat{\text{Cov}} = \frac{63}{64} \sum_{i=1}^{64} (\xi_i - \bar{\xi})^T (\xi_i - \bar{\xi}), \qquad (11) $$

where $\xi_i$ is a correlation function measured upon removal of the i'th jackknife subsample, $\bar{\xi}$ is the average of all jackknife measurements, and T denotes the conjugate transpose of the difference vector.

We exclude the final data point from assessments of signal significance, since the variance on these scales is unlikely to be captured by subvolumes of scales ~L-box/4 = 25 h⁻¹ Mpc. The resulting Hartlap correction factor (Hartlap, Simon & Schneider 2007) for estimates of the 6-bin signal inverse covariance is then ~0.9, and the maximum scale under consideration is ~17 h⁻¹ Mpc.

FP size residuals λ are not re-normalized, so as to preserve the amplitude of the signal, which otherwise would be inconsistent with the estimation of contamination to magnification signals. To explore the symmetry of intrinsic size correlation signals, we also measured various correlations for absolute sizes |λ|, and for sample selections λ_+ and λ_-, finding that most such splittings yielded signals corresponding to the intrinsic clustering of objects, merely down-weighted by powers of $\langle \lambda \rangle$.<fnref n="4" />

The exception to this statement is for galaxy clustering measured within λ_+ and λ_− samples, which shows variably significant differences for spiral and elliptical samples; meaning that, depending on the galaxy morphology, intrinsically large/small objects cluster differently in the simulation. We follow Joachimi et al. (2015) in defining a statistic to capture this variability, given as

$$ \Delta_{gg}(r) = \frac{\xi_{gg}(r, \lambda > 0)}{\xi_{gg}(r, \lambda < 0)} - 1. \qquad (12) $$

Retaining individual clustering measurements from jackknife resampling, we are able to estimate the covariance of this quantity directly, thus achieving a partial cancellation of shot-noise and cosmic variance.

Random catalogues are oversampled (20 × relative to their corresponding galaxy sample, given by subscripts on R) sets of points, uniform-randomly distributed in the simulation box with unit weights. All correlations and jackknife covariances are implemented within TREECORR (Jarvis et al. 2013), taking advantage of the functionality designed for scalar convergence fields κ, and observing periodic boundaries of the simulation box during pair-counting of three-dimensional correlations.

<footnote n="4">⁴One could re-subtract $\langle \lambda \rangle$ from the split distributions, but the meaning of the FP residual would then be lost; we neglect to explore the split-λ signals further, noting that any interesting asymmetries in the behaviour of λ_+, λ_− galaxies should manifest as non-zero signals measured on the original λ distributions.</footnote>
and feature stronger tilts than the elliptical planes. Still, correlations between FP residuals and key properties, such as the virial radius and stellar mass, are not drastic – particularly in the case of the Full FP. We thus advance to the measurement of intrinsic size correlations so-derived, applying caution as regards the impacts of tilting, λ-property correlations, and of our usage of virial radii to characterize spiral sizes.

## 4 INTRINSIC SIZE CORRELATIONS

We follow Joachimi et al. (2015) and Singh et al. (2021) in measuring galaxy clustering and intrinsic size correlations, though we do not integrate signals along the line-of-sight (the simulation box z-axis in our work), preferring to leverage our exact knowledge of object positions to assess correlations as functions of three-dimensional pair separations. The relevant correlation functions are estimated via

$$ \hat{\xi}_{gg}(r) = \frac{DD - 2DR_D + R_D R_D}{R_D R_D}, \qquad (5) $$

$$ \hat{\xi}_{g\lambda}(r) = \frac{D\lambda - R_D\lambda}{R_\lambda R_D}, \qquad (6) $$

$$ \hat{\xi}_{\lambda_1\lambda_2}(r) = \frac{\lambda_1\lambda_2}{R_{\lambda_1}R_{\lambda_2}}, \qquad (7) $$

where DD, RD, RR are unweighted pair-counts between galaxy D and random point R samples (with subscripts denoting galaxy samples to which the randoms correspond), for measurements of galaxy clustering. These, and the weighted pair-counts Dλ, R_Dλ, λ_1λ_2, are given by

$$ XY = \sum_{i,j} \Theta_{ij}, \qquad (8) $$

$$ X\lambda = \sum_{i,j} \lambda_i \Theta_{ij}, \qquad (9) $$

$$ \lambda_1\lambda_2 = \sum_{i,j} \lambda_{1,i}\lambda_{2,j} \Theta_{ij}, \qquad (10) $$

where samples X, Y denote the galaxy D or randoms R, R_D, R_λ samples, i, j are indices running over all objects in the two samples, and Θ_{ij} are binary filters applying a log-spaced binning in three-dimensional separation r. Notice equation (5) is the usual Landy-Szalay estimator (Landy & Szalay 1993). All estimators presented in equations (5), (6), and (7) are normalized by random-random pair counts (versus data-data and data-random pair counts). This is to allow direct comparison of our results to previous works (Joachimi et al. 2015; Singh et al. 2021), and to ensure that they can be interpreted and modelled as correlation functions in the sections that follow.

We measure the autoclustering correlations $\hat{\xi}_{gg}$ for each of our defined FP samples (Table 1), and also split those samples into intrinsically large λ_+ and small λ_- subsets for assessment of the clustering variation across the FP.

We define two additional ‘lens’ samples of elliptical galaxies, with masses >10^10 and >10^10.5 M_⊙ – the latter, more conservative mass selection matching that applied to our Masscut FP samples (see Section 2) – for use as density tracers in the density-size cross-correlations $\hat{\xi}_{g\lambda}$. The lens samples are similar to the elliptical FP samples, but retain galaxies lost to poor fitting of the FPs, thus offering slight reductions in shot-noise. Intrinsic size autocorrelations $\hat{\xi}_{\lambda\lambda}$ are measured only within FP samples, since λ probes different radii between ellipticals and spirals (Section 2), and is derived from differently fitted and tilted planes between the Masscut and Full set-ups.

We note that, whilst our fitted λ distributions are extremely close to symmetric by construction, any $\langle \lambda \rangle \neq 0$ will result in scale-invariant additive contributions to λ-weighted correlation functions, most visible at large separations where correlations are weak, and stemming from inherent clustering contributions (Singh et al. 2021). We ensure that this contribution is compensated by subtracting $\langle \lambda \rangle$ from each λ distribution under consideration, though in practice, for these samples, this subtraction is small and not necessary for the measurement of stable signals in the simulation box.

We work with seven three-dimensional separation bins in the range 0.5–30 h⁻¹ Mpc, and estimate covariances according to a delete-one jackknife resampling of 64 subvolumes defined on the simulation box. The covariance estimates are given by

$$ \widehat{\text{Cov}} = \frac{63}{64} \sum_{i=1}^{64} (\xi_i - \bar{\xi})^T (\xi_i - \bar{\xi}), \qquad (11) $$

where $\xi_i$ is a correlation function measured upon removal of the i'th jackknife subsample, $\bar{\xi}$ is the average of all jackknife measurements, and T denotes the conjugate transpose of the difference vector.

We exclude the final data point from assessments of signal significance, since the variance on these scales is unlikely to be captured by subvolumes of scales ~L-box/4 = 25 h⁻¹ Mpc. The resulting Hartlap correction factor (Hartlap, Simon & Schneider 2007) for estimates of the 6-bin signal inverse covariance is then ~0.9, and the maximum scale under consideration is ~17 h⁻¹ Mpc.

FP size residuals λ are not re-normalized, so as to preserve the amplitude of the signal, which otherwise would be inconsistent with the estimation of contamination to magnification signals. To explore the symmetry of intrinsic size correlation signals, we also measured various correlations for absolute sizes |λ|, and for sample selections λ_+ and λ_-, finding that most such splittings yielded signals corresponding to the intrinsic clustering of objects, merely down-weighted by powers of $\langle \lambda \rangle$.<fnref n="4" />

The exception to this statement is for galaxy clustering measured within λ_+ and λ_− samples, which shows variably significant differences for spiral and elliptical samples; meaning that, depending on the galaxy morphology, intrinsically large/small objects cluster differently in the simulation. We follow Joachimi et al. (2015) in defining a statistic to capture this variability, given as

$$ \Delta_{gg}(r) = \frac{\xi_{gg}(r, \lambda > 0)}{\xi_{gg}(r, \lambda < 0)} - 1. \qquad (12) $$

Retaining individual clustering measurements from jackknife resampling, we are able to estimate the covariance of this quantity directly, thus achieving a partial cancellation of shot-noise and cosmic variance.

Random catalogues are oversampled (20 × relative to their corresponding galaxy sample, given by subscripts on R) sets of points, uniform-randomly distributed in the simulation box with unit weights. All correlations and jackknife covariances are implemented within TREECORR (Jarvis et al. 2013), taking advantage of the functionality designed for scalar convergence fields κ, and observing periodic boundaries of the simulation box during pair-counting of three-dimensional correlations.

<footnote n="4">⁴One could re-subtract $\langle \lambda \rangle$ from the split distributions, but the meaning of the FP residual would then be lost; we neglect to explore the split-λ signals further, noting that any interesting asymmetries in the behaviour of λ_+, λ_− galaxies should manifest as non-zero signals measured on the original λ distributions.</footnote>
## 4.1 Lensing contamination

We focus now on ways in which ISCs might contaminate galaxy statistics by mimicking the size variations induced by weak lensing magnification, in a near-exact analogy to the intrinsic alignment phenomenon that contaminates studies of cosmic shear.

### 4.1.1 Density-magnification

Huff & Graves (2014) give an estimator for the weak lensing convergence $\hat{\kappa}$ as

$$ \log(1 + \hat{\kappa}) \equiv \Delta \log R, \qquad (13) $$

where $\Delta \log R = \log R_O - \log R_{FP}$ is the residual between an observed galaxy radius $\log R_O$ and an FP prediction thereof, given by $\log R_{FP}$. Cross-correlating this quantity with foreground lenses, as demonstrated by Huff & Graves (2014), one can construct a two-point estimate of the projected surface density $\Sigma$ at lens plane $z_l$, via

$$ \Sigma_{\text{crit}} \kappa = \Sigma(d_l \bar{\theta}), \qquad (14) $$

for the comoving distance to the lens plane $d_l$, and angular separation vector $\bar{\theta}$, with critical surface density

$$ \Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{d_a(z_s)}{d_a(z_l)(d_a(z_s) - d_a(z_l))(1 + z_l)^2}, \qquad (15) $$

where $d_a$ denote angular diameter distances, and subscripts s, l denote source (background) and lens (foreground) samples.

Intrinsic size fluctuations will contribute to the estimate $\hat{\kappa}$ (Alsing et al. 2015; Ciarlariello et al. 2015; Ciarlariello & Crittenden 2016). If these fluctuations do not vanish under spatial averaging, e.g. if the distribution is asymmetric, or if they are correlated with one another or with structure, then this may introduce biases into estimators such as equation (14) via some induced deviation from the true mean size of objects at a given redshift.

We return to equation (4), which re-defines the observed radius $R_O$ in terms of the lensing convergence $\kappa$, the FP-predicted radius $R_{FP}$, and the lensing-independent intrinsic size fluctuation $\lambda$, with $\langle \lambda \rangle \to 0$, i.e. the global distribution of intrinsic galaxy sizes is symmetric. This is an assumption, but one supported by the fitted low-z (and therefore low-lensing) FP of Joachimi et al. (2015), and borne out by our own simulation data.<fnref n="5" /> We then expand equation (13) as

$$ \log\{(1 + \kappa)(1 + \lambda)\} = \log \frac{R_O}{R_{FP}}. \qquad (16) $$

We note that $\lambda$ here (and in equation 4) is approximately equivalent to our definition in equation (3), modulo terms of $O(\lambda^2)$ (Joachimi et al. 2015). The estimator of equation (13) then becomes

$$ \Delta \log R \equiv \log(1 + \kappa + \lambda), \qquad (17) $$

where $\kappa + \lambda$ signifies that the Huff & Graves (2014) estimator is actually estimating the combination of intrinsic and lensing-derived size fluctuations. If intrinsic sizes are uncorrelated with structure, then one expects $\lambda$ to vanish under spatial averaging, and the unbiased convergence estimate of equation (13) is recovered. This

would not, however, preclude the possibility of autocorrelations in $\lambda$, which would have the potential to contaminate any magnification–magnification correlation estimates based on measured galaxy sizes (Section 4.1.2).

An ISC contamination of estimates for $\Sigma(r)$ may be of lesser concern, as one expects that lens and source samples should not feature ISCs if the constituent galaxies are not co-located in space; background intrinsic sizes know nothing about distant, foreground lenses. In analogy to the contamination of galaxy–galaxy lensing (GGL) by intrinsic alignments, one would only expect ISC contamination of $\Sigma(r)$ in the presence of large photometric redshift errors – however, such errors cannot be discounted, and are known to be more prevalent among the spiral galaxies that dominate weak lensing samples (Rozo et al. 2016).

In Huff & Graves (2014), angular cross-correlations $w_{ls}(\theta)$ between lens and source samples were used to estimate the fraction $f_l$ of foreground sources that had been scattered out to higher redshifts by photo-z errors. These galaxies acquire systematically biased FP residuals, as they are thought to be far more distant than they are in reality. The interloper fraction $f_l(z_l, z_s, \theta) \in [0, 1]$, per lens and source redshift $z_l$ and $z_s$, and angular $\theta$ bin, was thus used to down-weight the convergence estimate as $(1 - f_l)\kappa$; to weight a term in the estimator $f_l \Delta \log R_{err}$ that corrects the FP residual of a redshift interloper, and to outright discard lens-source bin pairs and angular scales for which the interloper fraction is very high.

In our simulated data, the convergence $\kappa = 0$. Thus we can estimate the amplitude of an ISC contribution $\Sigma_{\lambda}(r)$ by cross-correlating FP residuals $\Delta \log R$ with galaxy positions, assuming some distribution of photometric redshift errors. The hypothetical set-up is that galaxies situated at some low redshift have been mistaken for higher-redshift sources, and thus made their way into an estimate of the projected surface density. Assuming some value or functional form for $f_l$, which would in reality be estimated/modelled using data/simulations, we can also approximate the mitigation strategy of Huff & Graves (2014).

We note that $f_l$ could in principle vary widely between zero and unity, though all efforts will seek to minimize it. For a single $\theta$-bin $\in [0.6, 6]$ arcmin,<fnref n="6" /> Stölzner et al. (2022) estimate the lens-source versus lens–lens angular clustering ratio $w_{ls}/w_{ll}$ ($\neq f_l$) to increase with redshift, and to be as large as ~0.4 (~0.1), for spectroscopic lenses, without (with) an outlier mitigation strategy for the LSST DESC cosmoDC2 mock catalogue (Korytov et al. 2019). One expects larger ratios for photometric lens samples, as redshift errors should deprecate the autocorrelation $w_{ll}$ more severely than the cross-correlation $w_{ls}$. These ratios are not equivalent to $f_l$ – the fraction of ‘source’ objects that are in fact situated at the lens plane – but they do indicate that the strength of spurious, outlier-driven clustering is expected to remain relatively high for the foreseeable future. Noting also that Huff & Graves (2014) excluded any redshift-angular separation bins for which $f_l > 0.5$ (and found higher thresholds to change little), we explore a simple grid formed of three points in $f_l$ at 0.0,0.2,0.4, and a standalone case of $f_l = 1$, which is the true case for our hypothetical set-up; all of our ‘source’ objects are redshift interlopers, co-located with the lenses.

For simplicity, we henceforth assume insignificant redshift evolution in $f_l$, and in the intrinsic size field, and its correlations with itself and the density field. Thus we treat the box as if it were

<footnote n="5">⁵We might expect a violation of this symmetry in deep lensing data; at high redshifts, shape measurements necessitate resolution cuts that will exclude the small-R tail of objects, even if the true $\lambda$ distribution is symmetric. This is the source of the lensing and size biases (Schmidt et al. 2009a,b; Ciarlariello et al. 2015; Ciarlariello & Crittenden 2016), as well as an intrinsic size selection bias in the presence of ISCs.</footnote>

<footnote n="6">⁶Corresponding to scales of order a few hundred $h^{-1}$ kpc at $z \sim 0.1$, and a few/several $h^{-1}$ Mpc at $z \sim 2$, which are coincident with those observed by Huff & Graves (2014).</footnote>
## 4.1 Lensing contamination

We focus now on ways in which ISCs might contaminate galaxy statistics by mimicking the size variations induced by weak lensing magnification, in a near-exact analogy to the intrinsic alignment phenomenon that contaminates studies of cosmic shear.

### 4.1.1 Density-magnification

Huff & Graves (2014) give an estimator for the weak lensing convergence $\hat{\kappa}$ as

$$ \log(1 + \hat{\kappa}) \equiv \Delta \log R, \qquad (13) $$

where $\Delta \log R = \log R_O - \log R_{FP}$ is the residual between an observed galaxy radius $\log R_O$ and an FP prediction thereof, given by $\log R_{FP}$. Cross-correlating this quantity with foreground lenses, as demonstrated by Huff & Graves (2014), one can construct a two-point estimate of the projected surface density $\Sigma$ at lens plane $z_l$, via

$$ \Sigma_{\text{crit}} \kappa = \Sigma(d_l \bar{\theta}), \qquad (14) $$

for the comoving distance to the lens plane $d_l$, and angular separation vector $\bar{\theta}$, with critical surface density

$$ \Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{d_a(z_s)}{d_a(z_l)(d_a(z_s) - d_a(z_l))(1 + z_l)^2}, \qquad (15) $$

where $d_a$ denote angular diameter distances, and subscripts s, l denote source (background) and lens (foreground) samples.

Intrinsic size fluctuations will contribute to the estimate $\hat{\kappa}$ (Alsing et al. 2015; Ciarlariello et al. 2015; Ciarlariello & Crittenden 2016). If these fluctuations do not vanish under spatial averaging, e.g. if the distribution is asymmetric, or if they are correlated with one another or with structure, then this may introduce biases into estimators such as equation (14) via some induced deviation from the true mean size of objects at a given redshift.

We return to equation (4), which re-defines the observed radius $R_O$ in terms of the lensing convergence $\kappa$, the FP-predicted radius $R_{FP}$, and the lensing-independent intrinsic size fluctuation $\lambda$, with $\langle \lambda \rangle \to 0$, i.e. the global distribution of intrinsic galaxy sizes is symmetric. This is an assumption, but one supported by the fitted low-z (and therefore low-lensing) FP of Joachimi et al. (2015), and borne out by our own simulation data.<fnref n="5" /> We then expand equation (13) as

$$ \log\{(1 + \kappa)(1 + \lambda)\} = \log \frac{R_O}{R_{FP}}. \qquad (16) $$

We note that $\lambda$ here (and in equation 4) is approximately equivalent to our definition in equation (3), modulo terms of $O(\lambda^2)$ (Joachimi et al. 2015). The estimator of equation (13) then becomes

$$ \Delta \log R \equiv \log(1 + \kappa + \lambda), \qquad (17) $$

where $\kappa + \lambda$ signifies that the Huff & Graves (2014) estimator is actually estimating the combination of intrinsic and lensing-derived size fluctuations. If intrinsic sizes are uncorrelated with structure, then one expects $\lambda$ to vanish under spatial averaging, and the unbiased convergence estimate of equation (13) is recovered. This

would not, however, preclude the possibility of autocorrelations in $\lambda$, which would have the potential to contaminate any magnification–magnification correlation estimates based on measured galaxy sizes (Section 4.1.2).

An ISC contamination of estimates for $\Sigma(r)$ may be of lesser concern, as one expects that lens and source samples should not feature ISCs if the constituent galaxies are not co-located in space; background intrinsic sizes know nothing about distant, foreground lenses. In analogy to the contamination of galaxy–galaxy lensing (GGL) by intrinsic alignments, one would only expect ISC contamination of $\Sigma(r)$ in the presence of large photometric redshift errors – however, such errors cannot be discounted, and are known to be more prevalent among the spiral galaxies that dominate weak lensing samples (Rozo et al. 2016).

In Huff & Graves (2014), angular cross-correlations $w_{ls}(\theta)$ between lens and source samples were used to estimate the fraction $f_l$ of foreground sources that had been scattered out to higher redshifts by photo-z errors. These galaxies acquire systematically biased FP residuals, as they are thought to be far more distant than they are in reality. The interloper fraction $f_l(z_l, z_s, \theta) \in [0, 1]$, per lens and source redshift $z_l$ and $z_s$, and angular $\theta$ bin, was thus used to down-weight the convergence estimate as $(1 - f_l)\kappa$; to weight a term in the estimator $f_l \Delta \log R_{\text{err}}$ that corrects the FP residual of a redshift interloper, and to outright discard lens-source bin pairs and angular scales for which the interloper fraction is very high.

In our simulated data, the convergence $\kappa = 0$. Thus we can estimate the amplitude of an ISC contribution $\Sigma_{\lambda}(r)$ by cross-correlating FP residuals $\Delta \log R$ with galaxy positions, assuming some distribution of photometric redshift errors. The hypothetical set-up is that galaxies situated at some low redshift have been mistaken for higher-redshift sources, and thus made their way into an estimate of the projected surface density. Assuming some value or functional form for $f_l$, which would in reality be estimated/modelled using data/simulations, we can also approximate the mitigation strategy of Huff & Graves (2014).

We note that $f_l$ could in principle vary widely between zero and unity, though all efforts will seek to minimize it. For a single $\theta$-bin $\in [0.6, 6]$ arcmin,<fnref n="6" /> Stölzner et al. (2022) estimate the lens-source versus lens–lens angular clustering ratio $w_{ls}/w_{ll}$ ($\neq f_l$) to increase with redshift, and to be as large as ~0.4 (~0.1), for spectroscopic lenses, without (with) an outlier mitigation strategy for the LSST DESC cosmoDC2 mock catalogue (Korytov et al. 2019). One expects larger ratios for photometric lens samples, as redshift errors should deprecate the autocorrelation $w_{ll}$ more severely than the cross-correlation $w_{ls}$. These ratios are not equivalent to $f_l$ – the fraction of ‘source’ objects that are in fact situated at the lens plane – but they do indicate that the strength of spurious, outlier-driven clustering is expected to remain relatively high for the foreseeable future. Noting also that Huff & Graves (2014) excluded any redshift-angular separation bins for which $f_l > 0.5$ (and found higher thresholds to change little), we explore a simple grid formed of three points in $f_l$ at 0.0,0.2,0.4, and a standalone case of $f_l = 1$, which is the true case for our hypothetical set-up; all of our ‘source’ objects are redshift interlopers, co-located with the lenses.

For simplicity, we henceforth assume insignificant redshift evolution in $f_l$, and in the intrinsic size field, and its correlations with itself and the density field. Thus we treat the box as if it were

<footnote n="5">⁵We might expect a violation of this symmetry in deep lensing data; at high redshifts, shape measurements necessitate resolution cuts that will exclude the small-R tail of objects, even if the true $\lambda$ distribution is symmetric. This is the source of the lensing and size biases (Schmidt et al. 2009a,b; Ciarlariello et al. 2015; Ciarlariello & Crittenden 2016), as well as an intrinsic size selection bias in the presence of ISCs.</footnote>

<footnote n="6">⁶Corresponding to scales of order a few hundred $h^{-1}$ kpc at $z \sim 0.1$, and a few/several $h^{-1}$ Mpc at $z \sim 2$, which are coincident with those observed by Huff & Graves (2014).</footnote>
situated at redshifts up to $z \sim 0.2$ (and extrapolate model fits at $z = 0.06$ to $\sim 0.8$ in Section 4.1.2), and make predictions for ISC lensing contamination with the caveat that evolution of the intrinsic size field, cross-talking with $z$-dependent photo-$z$ quality, will pose further complications for a magnification analysis.<fnref n="7" /> Since we retain lower-redshift FP residuals when scattering objects out to $z_s$, the part of the FP error term $\Delta\log R_{\text{err}}$ that corrects $R_{\text{FP}}(z_s) \rightarrow R_{\text{FP}}(z_1)$ is effectively already applied, and without any error. We thus mimic the correction term simply according to the ratio of angular diameter distances $d_a$ at the mistaken source redshift and true lens redshift (Huff & Graves 2014)

$$ \Delta \log R_{\text{err}}(z_1, z_s) = \log \left( \frac{d_a(z_s) R_{\text{FP}}(z_1)}{d_a(z_1) R_{\text{FP}}(z_s)} \right) \\ \rightarrow \log \left( \frac{d_a(z_s)}{d_a(z_1)} \right). \qquad (18) $$

The final intrinsic size field estimator, featuring the photo-$z$ correction term, then becomes

$$ \Delta \log R = \log(1 + \hat{\lambda}) + f_i \Delta \log R_{\text{err}}(z_1, z_s). \qquad (19) $$

In order to estimate the intrinsic contribution, we translate the centre of the simulation box to a fixed comoving distance $\chi(z_1 = 0.2)$, and convert galaxy three-dimensional Cartesian coordinates into RA, Dec, and $\chi$. Size tracer objects are thus situated at $z \sim 0.2$, with corresponding RA, Dec, but we shall mistake them to be at some higher $z_s$, whilst all lenses are at $z_1 \sim 0.2$. This means that $\Sigma_{\lambda}(r)$ is most likely to feature at large-$r$, since small-angle lens-source pairs will be thought to probe highly non-linear physical scales on the lens plane that are difficult to model, and often excluded from GGL analyses (see Singh et al. 2020).

As mock lenses, we employ the Lens samples described in Section 4 (Table 1). For source galaxies, Huff & Graves (2014) selected 8.4 million photometric SDSS (York et al. 2000) elliptical galaxies between $z = 0.1-0.55$. We fit an analytic $n(z)$ to the reconstructed redshift distribution of the SDSS DR8 photometric sample (Sheldon et al. 2012), given by

$$ n(z) = z^{\alpha} \exp \left\{ - \left( \frac{z\sqrt{2}}{z_0} \right)^{\beta} \right\}, \qquad (20) $$

with $\alpha = 1.58$, $z_0 = 0.34$, $\beta = 1.86$. Assuming that our measurable density-FP residual correlations are representative across the redshift range of interest, we bootstrap our FP samples to estimate the uncertainty in $\Sigma_{\lambda}(r)$, neglecting cosmic variance; our predictions might thus be considered pessimistic. Each FP sample (Table 1) is re-sampled with replacement 1000 times, and the bootstrapped galaxies<fnref n="8" /> are assigned redshifts $z_s$ drawn from the fitted $n(z)$

(equation 20). We note that the hard $f_i < 0.5$ cut imposed by Huff & Graves (2014) functions to exclude regimes dominated by source galaxies that are physically coincident with lenses; we thus exclude galaxies assigned to $0.1 < z_s < 0.3$, i.e. we impose that no $z_s$ sits within $z_1 \pm 0.1$. We deliberately allow a small number of $z_s < 0.1$ objects, but find that they make no discernible difference to the measured correlation.

We estimate $\Sigma_{\lambda}(r)$ itself by cross-correlating the intrinsic size field estimates $\hat{\lambda}(f_i, z_s)$ (equation 19), weighted by critical surface densities $\Sigma_{\text{crit}}(z_1, z_s)$ (equation 13), with angular lens sample positions, collecting the weighted-average intrinsic size residual $\langle \hat{\lambda} \rangle(\vartheta)$ that would contaminate<fnref n="9" /> the convergence estimate $\langle \hat{k} \rangle(\vartheta)$. The contamination estimator is thus given as

$$ \hat{\Sigma}_{\lambda}(\vartheta) = \frac{\sum_{i,j} \Sigma_{\text{crit}}(z_i, z_j) \hat{\lambda}_j \langle i|j \rangle}{\sum_{i,j} \Sigma_{\text{crit}}(z_i, z_j)}, \qquad (21) $$

where $i, j$ index lens and source samples, respectively, $\hat{\lambda}_j$ gives the intrinsic size ($\hat{k}$-contaminant) estimate for galaxy $j$ (equation 19), and the bin-filter $\langle i|j \rangle$ here applies the angular binning $\theta$. Angular separations $\theta$ are then converted into transverse comoving separations<fnref n="10" /> $r$, all evaluated at $z_1 = 0.2$. We further measure the average size residual around randomly distributed points, for subtraction from the galaxy signal; this removes spurious large-scale signals, and improves the covariance properties of the estimator (Singh et al. 2017).

As for the box correlations (Section 4), here one must consider the mean FP residual $\langle \lambda \rangle$, which can contribute spurious signals derived from the clustering of galaxies. Observational galaxy FP analyses incorporate redshift dependence via polynomial fitting (Joachimi et al. 2015; Singh et al. 2021), or else fit the FP in narrow redshift bins (Huff & Graves 2014; Singh et al. 2021), thus achieving $\langle \lambda \rangle \sim 0$ over the redshift range, or enforcing it via subtraction. However, selections and effective weighting in $f_i$ could reintroduce a non-zero mean residual, if imposed after fitting of the FP. We assess the impact of the mean residual for our set-up by re-measuring $\Sigma_{\lambda}(r)$ after enforcing $\langle \lambda \rangle = 0$, and find that the difference in signal is negligible due to our subtraction of the signal measured around random points. The random signal subtraction is itself necessary to remove a spurious excess signal at large $r$, even if $\langle \lambda \rangle = 0$ is enforced.

For these angular correlations, we avoid a consideration of periodic boundaries by limiting the lens sample to be at least $5 h^{-1}$ Mpc away from the inner-edge of the box on the X-Y plane. At the lens redshift of $z_1 = 0.2$, this corresponds to an angle of $\sim 1.7$ arcmin, beyond which we discard measured correlations. Huff & Graves (2014) estimated the effect of photo-$z$-induced misestimation of $\Sigma_{\text{crit}}$ upon their signal to be less than 10 per cent, which should hold for our work (modulo noise) as we adopt their redshift distribution and consequent lensing geometry – as we shall see, the $f_i$-weighted correction term, and in particular our mass-selections, wield far larger influence over the signals in any case.

### 4.1.2 Magnification–magnification

We also consider the potential for ISCs to contaminate estimates of magnification autocorrelations $\kappa\kappa$. Again in analogy to the

<footnote n="7">⁷We note also that, if ISCs are in fact stronger at higher-$z$, then our estimates of contamination are conservative. If they are weaker at high-$z$, then their amplitude must have grown over time. Concurrently, the Universal star formation rate, and elliptical galaxy fraction, are evolving, particularly as observed in flux-limited samples contending with the Malmquist bias. The relative strength of spiral/elliptical ISCs would then be of interest with respect to the evolution of the contamination with redshift, similarly to cosmic shear contamination by IA (Fortuna et al. 2020).</footnote>

<footnote n="8">⁸We note that the elliptical samples are subsets of the lens samples given in Table 1, and would thus result in unrealistic, duplicated galaxy positions between foreground lenses and background sources. We verify that discarding lens objects that are duplicated in the source sample prior to measuring each of the bootstrap correlations has a negligible impact upon the signal prediction other than to raise the noise level.</footnote>

<footnote n="9">⁹This contamination will not function via simple addition, since the denominator in our estimator includes only the interloper objects; $\Sigma_{\lambda}(r)$ will contribute to any measured $\Sigma(r)$ only after a re-weighting according to the relative numbers of lens-source and lens-interloper pairs in each angular bin.</footnote>

<footnote n="10">¹⁰Presented in kpc, and not $h^{-1}$ kpc, for more direct comparison with Huff & Graves (2014).</footnote>
situated at redshifts up to $z \sim 0.2$ (and extrapolate model fits at $z = 0.06$ to $\sim 0.8$ in Section 4.1.2), and make predictions for ISC lensing contamination with the caveat that evolution of the intrinsic size field, cross-talking with $z$-dependent photo-$z$ quality, will pose further complications for a magnification analysis.<fnref n="7" /> Since we retain lower-redshift FP residuals when scattering objects out to $z_s$, the part of the FP error term $\Delta \log R_{\text{err}}$ that corrects $R_{\text{FP}}(z_s) \rightarrow R_{\text{FP}}(z_1)$ is effectively already applied, and without any error. We thus mimic the correction term simply according to the ratio of angular diameter distances $d_a$ at the mistaken source redshift and true lens redshift (Huff & Graves 2014)

$$ \Delta \log R_{\text{err}}(z_1, z_s) = \log \left( \frac{d_a(z_s) R_{\text{FP}}(z_1)}{d_a(z_1) R_{\text{FP}}(z_s)} \right) \\ \rightarrow \log \left( \frac{d_a(z_s)}{d_a(z_1)} \right). \qquad (18) $$

The final intrinsic size field estimator, featuring the photo-$z$ correction term, then becomes

$$ \Delta \log R = \log(1 + \hat{\lambda}) + f_i \Delta \log R_{\text{err}}(z_1, z_s). \qquad (19) $$

In order to estimate the intrinsic contribution, we translate the centre of the simulation box to a fixed comoving distance $\chi(z_1 = 0.2)$, and convert galaxy three-dimensional Cartesian coordinates into RA, Dec, and $\chi$. Size tracer objects are thus situated at $z \sim 0.2$, with corresponding RA, Dec, but we shall mistake them to be at some higher $z_s$, whilst all lenses are at $z_1 \sim 0.2$. This means that $\Sigma_{\lambda}(r)$ is most likely to feature at large-$r$, since small-angle lens-source pairs will be thought to probe highly non-linear physical scales on the lens plane that are difficult to model, and often excluded from GGL analyses (see Singh et al. 2020).

As mock lenses, we employ the Lens samples described in Section 4 (Table 1). For source galaxies, Huff & Graves (2014) selected 8.4 million photometric SDSS (York et al. 2000) elliptical galaxies between $z = 0.1-0.55$. We fit an analytic $n(z)$ to the reconstructed redshift distribution of the SDSS DR8 photometric sample (Sheldon et al. 2012), given by

$$ n(z) = z^{\alpha} \exp \left\{ - \left( \frac{z\sqrt{2}}{z_0} \right)^{\beta} \right\}, \qquad (20) $$

with $\alpha = 1.58$, $z_0 = 0.34$, $\beta = 1.86$. Assuming that our measurable density-FP residual correlations are representative across the redshift range of interest, we bootstrap our FP samples to estimate the uncertainty in $\Sigma_{\lambda}(r)$, neglecting cosmic variance; our predictions might thus be considered pessimistic. Each FP sample (Table 1) is re-sampled with replacement 1000 times, and the bootstrapped galaxies<fnref n="8" /> are assigned redshifts $z_s$ drawn from the fitted $n(z)$

(equation 20). We note that the hard $f_i < 0.5$ cut imposed by Huff & Graves (2014) functions to exclude regimes dominated by source galaxies that are physically coincident with lenses; we thus exclude galaxies assigned to $0.1 < z_s < 0.3$, i.e. we impose that no $z_s$ sits within $z_1 \pm 0.1$. We deliberately allow a small number of $z_s < 0.1$ objects, but find that they make no discernible difference to the measured correlation.

We estimate $\Sigma_{\lambda}(r)$ itself by cross-correlating the intrinsic size field estimates $\hat{\lambda}(f_i, z_s)$ (equation 19), weighted by critical surface densities $\Sigma_{\text{crit}}(z_1, z_s)$ (equation 13), with angular lens sample positions, collecting the weighted-average intrinsic size residual $\langle \hat{\lambda} \rangle(\vartheta)$ that would contaminate<fnref n="9" /> the convergence estimate $\langle \hat{k} \rangle(\vartheta)$. The contamination estimator is thus given as

$$ \hat{\Sigma}_{\lambda}(\vartheta) = \frac{\sum_{i,j} \Sigma_{\text{crit}}(z_i, z_j) \hat{\lambda}_j \langle i|j \rangle}{\sum_{i,j} \Sigma_{\text{crit}}(z_i, z_j)}, \qquad (21) $$

where $i, j$ index lens and source samples, respectively, $\hat{\lambda}_j$ gives the intrinsic size ($\hat{k}$-contaminant) estimate for galaxy $j$ (equation 19), and the bin-filter $\langle i|j \rangle$ here applies the angular binning $\theta$. Angular separations $\theta$ are then converted into transverse comoving separations<fnref n="10" /> $r$, all evaluated at $z_1 = 0.2$. We further measure the average size residual around randomly distributed points, for subtraction from the galaxy signal; this removes spurious large-scale signals, and improves the covariance properties of the estimator (Singh et al. 2017).

As for the box correlations (Section 4), here one must consider the mean FP residual $\langle \lambda \rangle$, which can contribute spurious signals derived from the clustering of galaxies. Observational galaxy FP analyses incorporate redshift dependence via polynomial fitting (Joachimi et al. 2015; Singh et al. 2021), or else fit the FP in narrow redshift bins (Huff & Graves 2014; Singh et al. 2021), thus achieving $\langle \lambda \rangle \sim 0$ over the redshift range, or enforcing it via subtraction. However, selections and effective weighting in $f_i$ could reintroduce a non-zero mean residual, if imposed after fitting of the FP. We assess the impact of the mean residual for our set-up by re-measuring $\Sigma_{\lambda}(r)$ after enforcing $\langle \lambda \rangle = 0$, and find that the difference in signal is negligible due to our subtraction of the signal measured around random points. The random signal subtraction is itself necessary to remove a spurious excess signal at large $r$, even if $\langle \lambda \rangle = 0$ is enforced.

For these angular correlations, we avoid a consideration of periodic boundaries by limiting the lens sample to be at least $5 h^{-1}$ Mpc away from the inner-edge of the box on the X-Y plane. At the lens redshift of $z_1 = 0.2$, this corresponds to an angle of $\sim 1.7$ arcmin, beyond which we discard measured correlations. Huff & Graves (2014) estimated the effect of photo-$z$-induced misestimation of $\Sigma_{\text{crit}}$ upon their signal to be less than 10 per cent, which should hold for our work (modulo noise) as we adopt their redshift distribution and consequent lensing geometry – as we shall see, the $f_i$-weighted correction term, and in particular our mass-selections, wield far larger influence over the signals in any case.

### 4.1.2 Magnification–magnification

We also consider the potential for ISCs to contaminate estimates of magnification autocorrelations $\kappa\kappa$. Again in analogy to the

<footnote n="7">⁷We note also that, if ISCs are in fact stronger at higher-$z$, then our estimates of contamination are conservative. If they are weaker at high-$z$, then their amplitude must have grown over time. Concurrently, the Universal star formation rate, and elliptical galaxy fraction, are evolving, particularly as observed in flux-limited samples contending with the Malmquist bias. The relative strength of spiral/elliptical ISCs would then be of interest with respect to the evolution of the contamination with redshift, similarly to cosmic shear contamination by IA (Fortuna et al. 2020).</footnote>

<footnote n="8">⁸We note that the elliptical samples are subsets of the lens samples given in Table 1, and would thus result in unrealistic, duplicated galaxy positions between foreground lenses and background sources. We verify that discarding lens objects that are duplicated in the source sample prior to measuring each of the bootstrap correlations has a negligible impact upon the signal prediction other than to raise the noise level.</footnote>

<footnote n="9">⁹This contamination will not function via simple addition, since the denominator in our estimator includes only the interloper objects; $\Sigma_{\lambda}(r)$ will contribute to any measured $\Sigma(r)$ only after a re-weighting according to the relative numbers of lens-source and lens-interloper pairs in each angular bin.</footnote>

<footnote n="10">¹⁰Presented in kpc, and not $h^{-1}$ kpc, for more direct comparison with Huff & Graves (2014).</footnote>
phenomena of IA, a density-intrinsic size correlation $\delta\lambda$ could induce spurious contributions to size-based $\kappa\kappa$ estimates over large separations in redshift, as the same lenses source both background convergence and foreground ISCs – the gravitational-intrinsic (GI) correlation, in IA parlance (Hirata & Seljak 2004) – whilst the impact of autocorrelations $\lambda\lambda$ would be largely limited to tomographic autocorrelations which include many closely-associated galaxy pairs – the intrinsic–intrinsic (II) correlation.

In tidal torque theory (Schafer 2009), intrinsic spiral galaxy ellipticities and alignments are expected to be pure orientation effects that will not yield a GI correlation (Tugendhat, Reischke & Schafer 2020; Ghosh et al. 2021). However, this null prediction will not necessarily hold for intrinsic spiral sizes, because the size relates to the isotropic part (i.e. the trace) of the tidal shear tensor; itself proportional to the local density via Poisson’s equation (Tugendhat & Schäfer 2018). Thus one should expect a GI magnification correlation from spirals, if indeed spiral sizes are correlated with the local density contrast $\delta$.

Moreover, one might expect a significant spiral intrinsic size autocorrelation to contribute strongly to redshift bin autocorrelations $\kappa\kappa$ as an II term, since spiral galaxies dominate deep lensing data; a fact that also makes exclusion of such objects difficult, given the deleterious impact upon the signal-to-noise of the desired measurement.

Lacking measurements of $\kappa\kappa$ correlations in the literature with which to compare, we elect instead to make a similar ansatz to those made by Joachimi et al. (2015) and Alsing et al. (2015)$^{11}$; that the intrinsic size field $\lambda(\vec{x})$ can be characterized by some linear scaling of the matter density contrast $\delta(\vec{x})$, parametrized separately for spiral and elliptical sizes as $B_{\text{spi}}$ and $B_{\text{ell}}$, respectively.

We fit this model to our measured $\lambda\lambda$ correlations in the range $3 < r < 17 h^{-1}$ Mpc, where the lower limit excludes highly non-linear scales, whilst the upper is dictated by our subsample covariance estimation procedure (see Section 4). We thus constrain the absolute values of the B parameters, the squares of which modulate the amplitude under our ansatz: $\xi_{\lambda\lambda} = B^2 \xi_{\delta\delta}$, where $\xi_{\delta\delta}$ is the matter autocorrelation function, computed for the Horizon-AGN cosmology using NBODYKIT$^{12}$ (Hand et al. 2018).

For comparison with Joachimi et al. (2015), we assume two simple Gaussian redshift distributions centred on $z = 0.4$, $0.8$ with widths $\sigma = 0.1$ – we thus assume the phenomenological model from $z = 0.06$ to hold out to $z \sim 0.8$, which is clearly optimistic. However, our predictions should be considered conservative in the case that ISCs weaken over cosmic time, and as we shall see, the landscape of spiral/elliptical predictions makes for interesting conclusions in any case. We convert the model into $C_{\ell}$ expectations for the auto $\lambda\lambda$ and cross-correlations $\lambda\kappa$, to be compared with the expected convergence signal $\kappa\kappa$ via (Joachimi et al. 2015)

$$ 
\begin{align}
C_{\kappa\kappa}^{ij}(\ell) &= \int_{0}^{\chi_{\text{hor}}} d\chi \frac{q^{i}(\chi) q^{j}(\chi)}{\chi^2} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, \chi \right), \\
C_{\lambda\kappa}^{ij}(\ell) &= B \int_{0}^{\chi_{\text{hor}}} d\chi \frac{p^{i}(\chi) q^{j}(\chi)}{\chi^2} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, \chi \right), \\
C_{\lambda\lambda}^{ij}(\ell) &= B^2 \int_{0}^{\chi_{\text{hor}}} d\chi \frac{p^{i}(\chi) p^{j}(\chi)}{\chi^2} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, \chi \right), \tag{22}
\end{align}
$$

where we take $P_{\delta}$ as the matter power spectrum with non-linear corrections$^{13}$ (Smith et al. 2003; Takahashi et al. 2012), $p^i(\chi)$ is the probability distribution of comoving distances $\chi(z)$ in the $i$th tomographic sample, $\chi_{\text{hor}}$ is the comoving distance to the horizon, $\ell$ is the angular multipole, $B$ is the relevant linear intrinsic size field parameter, and $q^i(\chi)$ is the lensing efficiency of sample $i$, given by

$$ q(\chi) = \frac{3H_0^2 \Omega_m}{2c^2} \frac{\chi}{a(\chi)} \int_{\chi}^{\chi_{\text{hor}}} d\chi' p(\chi') \frac{\chi' - \chi}{\chi'}, \quad (23) $$

assuming a flat universe, with present-day Hubble parameter $H_0$ and matter energy-density fraction $\Omega_m$, scale factor $a$, and speed of light $c$.

We note here that a physically-motivated model for these correlations is highly desirable both for prospective studies of cosmic convergence, and for gaining insight into galaxy evolution from intrinsic sizes. For example, Ghosh et al. (2021) construct a unified linear model of elliptical galaxy intrinsic alignments and size correlations, asserting that the ‘elasticity’ of ellipticals – the constant of proportionality $D$ between observed shapes/sizes and the magnitude and orientation of the tidal shear (Tugendhat & Schäfer 2018) – is responsible for both.

A more complex and promising avenue is the effective field theory (EFT) model of Vlah, Chisari & Schmidt (2020), which is capable of modelling scalar biased tracers – such as galaxy sizes, and other properties – in principle accounting for all contributions up to a chosen order. Yielding coefficients describing the strength of contributions from all possible field combinations, insights into the evolution of galaxy sizes could be readily derived. An application of this model to ISCs, in an extended simulation analysis, could thus be of great value.

# 5 RESULTS AND DISCUSSION

Here, we detail the results of our measurements of intrinsic size correlations, and our so-derived predictions for contamination of cosmic convergence statistics. We emphasize here the difference between absolute size and intrinsic size as we have defined it (Section 3.1), noting that an absolutely large object can be intrinsically small (towards the right of a panel but above the cyan line, in Fig. 3) in comparison with its fellows of similar morphology and characteristics. We also note that many previous works explore the variability of galaxy sizes at fixed stellar mass. As Fig. 5 shows, this is similar to considering the variability of FP residuals within an interval centred on some radius; that is, examining samples at fixed stellar mass is roughly akin to examining distinct regions of the FP.

## 5.1 Intrinsic size correlations

Fig. 6 displays the three-dimensional two-point correlations measured for our various samples (Table 1) defined in the Horizon-AGN simulation box. The top (bottom) row shows correlations measured for spiral (elliptical) FP samples, and columns show the galaxy clustering $\xi_{gg}$, the $\lambda_+$ versus $\lambda_-$ clustering difference $\Delta_{gg}$, the lens density-FP sample size correlations $\xi_{g\lambda}$, and the size autocorrelation $\xi_{\lambda\lambda}$, respectively. Large points and solid lines depict conservative

<footnote>$^{11}$Though Alsing et al. (2015) modelled intrinsic and magnification-induced size variations according to joint distributions of absolute size and flux, as opposed to deviations from the FP.</footnote>

<footnote>$^{12}$https://nbodykit.readthedocs.io</footnote>

<footnote>$^{13}$We note that Alsing et al. (2015) used the linear power spectrum for II and the geometric mean of non-/linear spectra for GI ISCs, as suggested by Kirk et al. (2012) in the context of IA. Recent years have seen declining usage of this model in IA contexts, and we elect to use the full non-linear power spectrum here.</footnote>
phenomena of IA, a density-intrinsic size correlation $\delta\lambda$ could induce spurious contributions to size-based $\kappa\kappa$ estimates over large separations in redshift, as the same lenses source both background convergence and foreground ISCs – the gravitational-intrinsic (GI) correlation, in IA parlance (Hirata & Seljak 2004) – whilst the impact of autocorrelations $\lambda\lambda$ would be largely limited to tomographic autocorrelations which include many closely-associated galaxy pairs – the intrinsic–intrinsic (II) correlation.

In tidal torque theory (Schafer 2009), intrinsic spiral galaxy ellipticities and alignments are expected to be pure orientation effects that will not yield a GI correlation (Tugendhat, Reischke & Schafer 2020; Ghosh et al. 2021). However, this null prediction will not necessarily hold for intrinsic spiral sizes, because the size relates to the isotropic part (i.e. the trace) of the tidal shear tensor; itself proportional to the local density via Poisson’s equation (Tugendhat & Schäfer 2018). Thus one should expect a GI magnification correlation from spirals, if indeed spiral sizes are correlated with the local density contrast $\delta$.

Moreover, one might expect a significant spiral intrinsic size autocorrelation to contribute strongly to redshift bin autocorrelations $\kappa\kappa$ as an II term, since spiral galaxies dominate deep lensing data; a fact that also makes exclusion of such objects difficult, given the deleterious impact upon the signal-to-noise of the desired measurement.

Lacking measurements of $\kappa\kappa$ correlations in the literature with which to compare, we elect instead to make a similar ansatz to those made by Joachimi et al. (2015) and Alsing et al. (2015)$^{11}$; that the intrinsic size field $\lambda(\vec{x})$ can be characterized by some linear scaling of the matter density contrast $\delta(\vec{x})$, parametrized separately for spiral and elliptical sizes as $B_{\text{spi}}$ and $B_{\text{ell}}$, respectively.

We fit this model to our measured $\lambda\lambda$ correlations in the range $3 < r < 17 h^{-1}$ Mpc, where the lower limit excludes highly non-linear scales, whilst the upper is dictated by our subsample covariance estimation procedure (see Section 4). We thus constrain the absolute values of the B parameters, the squares of which modulate the amplitude under our ansatz: $\xi_{\lambda\lambda} = B^2 \xi_{\delta\delta}$, where $\xi_{\delta\delta}$ is the matter autocorrelation function, computed for the Horizon-AGN cosmology using NBODYKIT$^{12}$ (Hand et al. 2018).

For comparison with Joachimi et al. (2015), we assume two simple Gaussian redshift distributions centred on $z = 0.4$, $0.8$ with widths $\sigma = 0.1$ – we thus assume the phenomenological model from $z = 0.06$ to hold out to $z \sim 0.8$, which is clearly optimistic. However, our predictions should be considered conservative in the case that ISCs weaken over cosmic time, and as we shall see, the landscape of spiral/elliptical predictions makes for interesting conclusions in any case. We convert the model into $C_{\ell}$ expectations for the auto $\lambda\lambda$ and cross-correlations $\lambda\kappa$, to be compared with the expected convergence signal $\kappa\kappa$ via (Joachimi et al. 2015)

$$ 
\begin{align} 
C_{\kappa\kappa}^{ij}(\ell) &= \int_{0}^{\chi_{\text{hor}}} d\chi \frac{q^i(\chi) q^j(\chi)}{\chi^2} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, \chi \right), \nonumber \\ 
C_{\lambda\kappa}^{ij}(\ell) &= B \int_{0}^{\chi_{\text{hor}}} d\chi \frac{p^i(\chi) q^j(\chi)}{\chi^2} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, \chi \right), \nonumber \\ 
C_{\lambda\lambda}^{ij}(\ell) &= B^2 \int_{0}^{\chi_{\text{hor}}} d\chi \frac{p^i(\chi) p^j(\chi)}{\chi^2} P_{\delta} \left( \frac{\ell + 1/2}{\chi}, \chi \right), \tag{22}
\end{align} 
$$

where we take $P_{\delta}$ as the matter power spectrum with non-linear corrections$^{13}$ (Smith et al. 2003; Takahashi et al. 2012), $p^i(\chi)$ is the probability distribution of comoving distances $\chi(z)$ in the $i$th tomographic sample, $\chi_{\text{hor}}$ is the comoving distance to the horizon, $\ell$ is the angular multipole, $B$ is the relevant linear intrinsic size field parameter, and $q^i(\chi)$ is the lensing efficiency of sample $i$, given by

$$ q(\chi) = \frac{3H_0^2 \Omega_m}{2c^2} \frac{\chi}{a(\chi)} \int_{\chi}^{\chi_{\text{hor}}} d\chi' p(\chi') \frac{\chi' - \chi}{\chi'}, \quad (23) $$

assuming a flat universe, with present-day Hubble parameter $H_0$ and matter energy-density fraction $\Omega_m$, scale factor $a$, and speed of light $c$.

We note here that a physically-motivated model for these correlations is highly desirable both for prospective studies of cosmic convergence, and for gaining insight into galaxy evolution from intrinsic sizes. For example, Ghosh et al. (2021) construct a unified linear model of elliptical galaxy intrinsic alignments and size correlations, asserting that the ‘elasticity’ of ellipticals – the constant of proportionality $D$ between observed shapes/sizes and the magnitude and orientation of the tidal shear (Tugendhat & Schäfer 2018) – is responsible for both.

A more complex and promising avenue is the effective field theory (EFT) model of Vlah, Chisari & Schmidt (2020), which is capable of modelling scalar biased tracers – such as galaxy sizes, and other properties – in principle accounting for all contributions up to a chosen order. Yielding coefficients describing the strength of contributions from all possible field combinations, insights into the evolution of galaxy sizes could be readily derived. An application of this model to ISCs, in an extended simulation analysis, could thus be of great value.

# 5 RESULTS AND DISCUSSION

Here, we detail the results of our measurements of intrinsic size correlations, and our so-derived predictions for contamination of cosmic convergence statistics. We emphasize here the difference between absolute size and intrinsic size as we have defined it (Section 3.1), noting that an absolutely large object can be intrinsically small (towards the right of a panel but above the cyan line, in Fig. 3) in comparison with its fellows of similar morphology and characteristics. We also note that many previous works explore the variability of galaxy sizes at fixed stellar mass. As Fig. 5 shows, this is similar to considering the variability of FP residuals within an interval centred on some radius; that is, examining samples at fixed stellar mass is roughly akin to examining distinct regions of the FP.

## 5.1 Intrinsic size correlations

Fig. 6 displays the three-dimensional two-point correlations measured for our various samples (Table 1) defined in the Horizon-AGN simulation box. The top (bottom) row shows correlations measured for spiral (elliptical) FP samples, and columns show the galaxy clustering $\xi_{gg}$, the $\lambda_+$ versus $\lambda_-$ clustering difference $\Delta_{gg}$, the lens density-FP sample size correlations $\xi_{g\lambda}$, and the size autocorrelation $\xi_{\lambda\lambda}$, respectively. Large points and solid lines depict conservative

<footnote>$^{11}$Though Alsing et al. (2015) modelled intrinsic and magnification-induced size variations according to joint distributions of absolute size and flux, as opposed to deviations from the FP.</footnote>

<footnote>$^{12}$https://nbodykit.readthedocs.io</footnote>

<footnote>$^{13}$We note that Alsing et al. (2015) used the linear power spectrum for II and the geometric mean of non-/linear spectra for GI ISCs, as suggested by Kirk et al. (2012) in the context of IA. Recent years have seen declining usage of this model in IA contexts, and we elect to use the full non-linear power spectrum here.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** The Full elliptical sample Fundamental Plane (as in top-left-hand panel of Fig. 3), coloured according to six uniform bins in log-stellar mass log ($M_*/M_\odot$). The 1:1 relation, a linear fit to all points on the plane (Section 3.2), and contours illustrating $M_* < 10^{10.5} M_\odot$ objects (Section 2), are reproduced from Fig. 3, here in black. Annotations $\lambda_-$ and $\lambda_+$ denote the sides of the plane corresponding to intrinsically small and large objects, respectively. Selections in stellar mass are seen to isolate regions of the Fundamental Plane, and similar trends are seen for spiral planes.</figcaption></figure>

estimates, measured for FP and lens samples with a mass-selection $M_* > 10^{10.5} M_\odot$ (Section 2; Hatfield et al. 2019), whilst small points and dotted lines include all FP galaxies, and elliptical lenses down to $10^{10} M_\odot$. Grey shading marks scales $> 17 h^{-1} \text{ Mpc}$, beyond which uncertainties are likely to be underestimated by our jackknife procedure (Section 4) – we thus exclude these data points from reports of statistical significance.

We observe strong non-linearity in the galaxy clustering correlations (see also Lee et al. 2021a), such that we are unable to make good fits of linear bias models to any of the measured $\xi_{gg}$. As such, we do not follow Joachimi et al. (2015) and Singh et al. (2021) in fitting jointly to the clustering $\xi_{gg} \propto b_g^2$ and density-size cross-correlation $\xi_{g\lambda} \propto b_g B$, thereby breaking a degeneracy between the galaxy bias $b_g$ and the ISC parameter $B$; we prefer instead to model the $\lambda\lambda$ correlations (see Section 5.5).

### 5.1.1 Spiral intrinsic size correlations

Our measured spiral galaxy clustering reveals differences (Fig. 6; top-left-hand panel) when the samples are split into intrinsically large $\lambda_+$ and small $\lambda_-$ objects for the measurement of the $\Delta_{gg}$ function (top-middle-left-hand panel). In the Masscut case, low signal-to-noise and inter-scale correlations result in a $\Delta_{gg}$ signal that is not significantly non-zero ($0.6\sigma$; large points). However, in the Full FP case (small points) a $\Delta_{gg}$ signal of similar form is detected at a significance of $5.9\sigma$, suggesting that, when lower-mass spirals ($M_* < 10^{10.5} M_\odot$) are included, spiral galaxies with intrinsically large virial radii $\lambda_+$ are more strongly clustered than their small-radius counterparts $\lambda_-$ over scales $\lesssim 10 h^{-1} \text{ Mpc}$.

This detection of a difference in clustering across the spiral FP is supported by a $5.2\sigma$ detection of a positive density-intrinsic size correlation $\xi_{g\lambda}$ (top-middle-right-hand panel), the significance of which is similarly reduced (to $1.6\sigma$) upon removal of the lower-mass spirals for the Masscut FP; thus larger spirals are more apt to be found at the peaks of the density distribution, as traced by the lens sample, with the significant detection again relying upon the inclusion of lower-mass spiral galaxies.

Cebrián & Trujillo (2014) found ~10 per cent larger $R_{\text{eff}}$ for low-z late-type galaxies in the field, and Matharu et al. (2019) found z ~1 star-forming cluster galaxies to be 17 per cent smaller than field spirals, both at fixed stellar mass. We consider the stellar particle virial radius $R_{\text{vir*}}$ here, and find $R_{\text{eff}} \sim (0.2-0.3) R_{\text{vir*}}^{0.5-0.6}$ for spiral galaxies (scatter ~0.2 dex; see Section 2.2), such that the variation of absolute radius with density is unlikely to differ between virial/effective radii. Thus, unless the intrinsic variations in $R_{\text{vir*}}$ and $R_{\text{eff}}$ with density (akin to those found by fixing stellar mass; Fig. 5) have opposite signs, the simulation disagrees with the findings of Cebrián & Trujillo (2014) and Matharu et al. (2019) – modulo differences in sample selection, and in the different controlling variables ($M_*$ versus expected radius).

We make a $9\sigma$ detection of an intrinsic size autocorrelation $\xi_{\lambda\lambda}$ (Fig. 6, top-right-hand panel) for the Full spiral FP, and in contrast to $\Delta_{gg}$ and $\xi_{g\lambda}$, this signal remains significant at ~3.4$\sigma$ even after removing lower-mass spirals. Thus the simulation strongly suggests that intrinsically large or small spiral galaxies tend to have similarly-sized neighbours, as reckoned by the stellar particle virial radius (Section 2.2).

As Fig. 4 shows, any correlation between spiral $\lambda$ and stellar mass is weaker for the Full than for the Masscut spiral FP. It therefore seems unlikely that these trends are simple manifestations of mass-selection and changing galaxy bias, since the Full FP sample is far more evenly distributed on the $M_*-\lambda$ plane. The tilt in the spiral FP might offer some explanation, as it results in a tendency of intrinsically large spirals to also be absolutely large, and absolute size does correlate with stellar mass (the likely source of stellar mass-$\lambda$ correlations in Fig. 4). However, the galaxies responsible for these significant signals have typically small-to-intermediate virial radii in absolute terms, and are distributed far more symmetrically about the spiral Fundamental Plane (see Fig. 3; blue contours in top-right-hand panel), such that the tilt alone cannot satisfactorily explain the array of significant signals from the Full spiral FP.

We thus conclude that spiral ISCs are in fact present in the simulation. Here, we must reiterate that the lower-mass objects ($M_* < 10^{10.5} M_\odot$; Section 2) most strongly revealing these size correlations are contributors to a general underestimation of the galaxy clustering amplitude by Horizon-AGN, as compared with observations (Hatfield et al. 2019). However, their symmetric distribution around the FP, the lack of $\lambda$ versus property correlations, and the persistence of the $\lambda\lambda$ signal after limiting to higher-mass spirals, prevent us from discounting these correlations.

### 5.1.2 Elliptical intrinsic size correlations

The signals we measure for elliptical FP samples are more sensitive to our analysis choices with respect to sample selections and the Fundamental Plane. We find that the usage of galaxy surface mass density $\Sigma$ in the FP (as in Rosito et al. 2021), as opposed to the surface brightness $I$, effectively erases measurable ISCs when holding the galaxy sample constant. We also see great variability of FP fits and measured ISCs with respect to the inclusion of smaller, lower-mass elliptical galaxies in the simulation.

Briefly, relaxation of the global elliptical mass-selection $M_* > 10^{10} M_\odot$ rapidly drives steep, positive $\Delta_{gg}$ correlations on small scales, becoming negative at intermediate scales, whilst simultaneously erasing any $\xi_{g\lambda}$ and $\xi_{\lambda\lambda}$. Moreover, the total clustering signal $\xi_{gg}$ outstrips both the $\lambda_+$ and $\lambda_-$ clustering profiles on small scales, with $\lambda_+$ rising to dominate on large scales. Given the confused correlation picture which aligns poorly with literature findings;
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The Full elliptical sample Fundamental Plane (as in top-left-hand panel of Fig. 3), coloured according to six uniform bins in log-stellar mass log ($M_*/M_\odot$). The 1:1 relation, a linear fit to all points on the plane (Section 3.2), and contours illustrating $M_* < 10^{10.5} M_\odot$ objects (Section 2), are reproduced from Fig. 3, here in black. Annotations $\lambda_-$ and $\lambda_+$ denote the sides of the plane corresponding to intrinsically small and large objects, respectively. Selections in stellar mass are seen to isolate regions of the Fundamental Plane, and similar trends are seen for spiral planes.</figcaption></figure>

estimates, measured for FP and lens samples with a mass-selection $M_* > 10^{10.5} M_\odot$ (Section 2; Hatfield et al. 2019), whilst small points and dotted lines include all FP galaxies, and elliptical lenses down to $10^{10} M_\odot$. Grey shading marks scales $> 17 h^{-1} \text{ Mpc}$, beyond which uncertainties are likely to be underestimated by our jackknife procedure (Section 4) – we thus exclude these data points from reports of statistical significance.

We observe strong non-linearity in the galaxy clustering correlations (see also Lee et al. 2021a), such that we are unable to make good fits of linear bias models to any of the measured $\xi_{gg}$. As such, we do not follow Joachimi et al. (2015) and Singh et al. (2021) in fitting jointly to the clustering $\xi_{gg} \propto b_g^2$ and density-size cross-correlation $\xi_{g\lambda} \propto b_g B$, thereby breaking a degeneracy between the galaxy bias $b_g$ and the ISC parameter $B$; we prefer instead to model the $\lambda\lambda$ correlations (see Section 5.5).

### 5.1.1 Spiral intrinsic size correlations

Our measured spiral galaxy clustering reveals differences (Fig. 6; top-left-hand panel) when the samples are split into intrinsically large $\lambda_+$ and small $\lambda_-$ objects for the measurement of the $\Delta_{gg}$ function (top-middle-left-hand panel). In the Masscut case, low signal-to-noise and inter-scale correlations result in a $\Delta_{gg}$ signal that is not significantly non-zero ($0.6\sigma$; large points). However, in the Full FP case (small points) a $\Delta_{gg}$ signal of similar form is detected at a significance of $5.9\sigma$, suggesting that, when lower-mass spirals ($M_* < 10^{10.5} M_\odot$) are included, spiral galaxies with intrinsically large virial radii $\lambda_+$ are more strongly clustered than their small-radius counterparts $\lambda_-$ over scales $\lesssim 10 h^{-1} \text{ Mpc}$.

This detection of a difference in clustering across the spiral FP is supported by a $5.2\sigma$ detection of a positive density-intrinsic size correlation $\xi_{g\lambda}$ (top-middle-right-hand panel), the significance of which is similarly reduced (to $1.6\sigma$) upon removal of the lower-mass spirals for the Masscut FP; thus larger spirals are more apt to be found at the peaks of the density distribution, as traced by the lens sample, with the significant detection again relying upon the inclusion of lower-mass spiral galaxies.

Cebrián & Trujillo (2014) found ~10 per cent larger $R_{\text{eff}}$ for low-z late-type galaxies in the field, and Matharu et al. (2019) found z ~1 star-forming cluster galaxies to be 17 per cent smaller than field spirals, both at fixed stellar mass. We consider the stellar particle virial radius $R_{\text{vir*}}$ here, and find $R_{\text{eff}} \sim (0.2-0.3) R_{\text{vir*}}^{0.5-0.6}$ for spiral galaxies (scatter ~0.2 dex; see Section 2.2), such that the variation of absolute radius with density is unlikely to differ between virial/effective radii. Thus, unless the intrinsic variations in $R_{\text{vir*}}$ and $R_{\text{eff}}$ with density (akin to those found by fixing stellar mass; Fig. 5) have opposite signs, the simulation disagrees with the findings of Cebrián & Trujillo (2014) and Matharu et al. (2019) – modulo differences in sample selection, and in the different controlling variables ($M_*$ versus expected radius).

We make a $9\sigma$ detection of an intrinsic size autocorrelation $\xi_{\lambda\lambda}$ (Fig. 6, top-right-hand panel) for the Full spiral FP, and in contrast to $\Delta_{gg}$ and $\xi_{g\lambda}$, this signal remains significant at ~3.4$\sigma$ even after removing lower-mass spirals. Thus the simulation strongly suggests that intrinsically large or small spiral galaxies tend to have similarly-sized neighbours, as reckoned by the stellar particle virial radius (Section 2.2).

As Fig. 4 shows, any correlation between spiral $\lambda$ and stellar mass is weaker for the Full than for the Masscut spiral FP. It therefore seems unlikely that these trends are simple manifestations of mass-selection and changing galaxy bias, since the Full FP sample is far more evenly distributed on the $M_*-\lambda$ plane. The tilt in the spiral FP might offer some explanation, as it results in a tendency of intrinsically large spirals to also be absolutely large, and absolute size does correlate with stellar mass (the likely source of stellar mass-$\lambda$ correlations in Fig. 4). However, the galaxies responsible for these significant signals have typically small-to-intermediate virial radii in absolute terms, and are distributed far more symmetrically about the spiral Fundamental Plane (see Fig. 3; blue contours in top-right-hand panel), such that the tilt alone cannot satisfactorily explain the array of significant signals from the Full spiral FP.

We thus conclude that spiral ISCs are in fact present in the simulation. Here, we must reiterate that the lower-mass objects ($M_* < 10^{10.5} M_\odot$; Section 2) most strongly revealing these size correlations are contributors to a general underestimation of the galaxy clustering amplitude by Horizon-AGN, as compared with observations (Hatfield et al. 2019). However, their symmetric distribution around the FP, the lack of $\lambda$ versus property correlations, and the persistence of the $\lambda\lambda$ signal after limiting to higher-mass spirals, prevent us from discounting these correlations.

### 5.1.2 Elliptical intrinsic size correlations

The signals we measure for elliptical FP samples are more sensitive to our analysis choices with respect to sample selections and the Fundamental Plane. We find that the usage of galaxy surface mass density $\Sigma$ in the FP (as in Rosito et al. 2021), as opposed to the surface brightness $I$, effectively erases measurable ISCs when holding the galaxy sample constant. We also see great variability of FP fits and measured ISCs with respect to the inclusion of smaller, lower-mass elliptical galaxies in the simulation.

Briefly, relaxation of the global elliptical mass-selection $M_* > 10^{10} M_\odot$ rapidly drives steep, positive $\Delta_{gg}$ correlations on small scales, becoming negative at intermediate scales, whilst simultaneously erasing any $\xi_{g\lambda}$ and $\xi_{\lambda\lambda}$. Moreover, the total clustering signal $\xi_{gg}$ outstrips both the $\lambda_+$ and $\lambda_-$ clustering profiles on small scales, with $\lambda_+$ rising to dominate on large scales. Given the confused correlation picture which aligns poorly with literature findings;
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Galaxy two-point correlations, measured as functions of three-dimensional separations *r*, for Fundamental Plane (Section 3) samples selected from the Horizon-AGN simulation (Section 2). For spiral (top) and elliptical (bottom) galaxies, columns give the galaxy clustering $\xi_{gg}(r)$, the fractional difference in clustering $\Delta_{gg}(r)$ between intrinsically large/small galaxies, the density-intrinsic size cross-correlation $\xi_{g\lambda}(r)$, and the intrinsic size autocorrelation $\xi_{\lambda\lambda}(r)$, respectively. The clustering panels also display the large (small) galaxy clustering signals $\xi_{g\lambda}+g_{\lambda+}(\xi_{g\lambda}-g_{\lambda-})$ from which $\Delta_{gg}$ is derived (equation 12). Larger data points denote more conservative signals measured for Masscut Fundamental Plane samples ($M_* > 10^{10.5} M_{\odot}$; Section 2), whilst smaller points are measured using the Full Fundamental Plane samples featuring lower-mass galaxies. Pink solid (dotted) curves in the right-most panels give the results of fitting the 1-parameter phenomenological model (described in Section 4.1.2) separately to elliptical/spiral $\lambda\lambda$ correlations measured in the Masscut (Full) Fundamental Plane samples. Intrinsic size correlations are highly sensitive to galaxies' morphology and stellar mass range.</figcaption></figure>

poorer FP fits and induced asymmetry in $\lambda$ (Section 2.3); the steeply increasing stellar mass-function (as $M_*$ decreases from $10^{10} M_{\odot}$) of simulated ellipticals (Fig. 1), which are known to be overproduced at the low-mass end (Kaviraj et al. 2017); and the aforementioned concerns around low-mass clustering in Horizon-AGN, we opt to exclude these objects from our Full FP sample (and Lens samples; Table 1), imposing that elliptical stellar masses $M_* > 10^{10} M_{\odot}$.

For our fiducial choices, and upon inclusion of lower-mass ellipticals with $10^{10} < M_*/M_{\odot} < 10^{10.5}$, we are in agreement with the observational findings of Joachimi et al. (2015). That work studied fainter elliptical galaxies at low redshifts ($z < 0.2$), and made significant detections of negative $\Delta_{gg}$ and $g\lambda$, as well as a tentative detection of positive $\lambda\lambda$. We reproduce the negative $\Delta_{gg}$ detection at $\sim 5.6\sigma$ when considering the Full FP sample, and this is accompanied by weak hints of negative $\xi_{g\lambda}$ ($1.7\sigma$) and positive $\xi_{\lambda\lambda}$ ($1.4\sigma$), wherein scales are highly correlated. Thus we observe ellipticals with intrinsically small effective radii (Section 2.1) to be more clustered than their large counterparts, in opposition to the virial radius trend for spirals, and over a slightly broader range of scales. However, none of these signals are detected at $>0.7\sigma$ when the mass-selection $M_* > 10^{10.5} M_{\odot}$ is imposed.

Singh et al. (2021) detected positive $g\lambda$ correlations in higher-redshift ($0.16 < z < 0.7$) luminous red galaxy (LRG) data. However, their findings are compatible with ours and those of Joachimi et al. (2015), given that their sample is comparatively dominated by brighter, more massive objects; they were able to reproduce the negative signal seen by Joachimi et al. (2015) by making selections in galaxy luminosity and colour. Interestingly, our Masscut FP signals weaken towards zero as compared with the Full FP signals in each case. This is agreeable with the findings of Singh et al. (2021), since many of their low-luminosity selections – bringing their sample more into parity with our own – yield signals consistent with zero, or weakly negative ($\sim 2\sigma$) in the faintest cases.

Both Joachimi et al. (2015) and Singh et al. (2021) argued that the consistent assignment of central and satellite galaxies to $\lambda_+$ and $\lambda_-$, respectively, drove the $g\lambda$ trend; that is, elliptical FPs and their residuals are dependent upon galaxies' environments (Bernardi et al. 2003; D'Onofrio et al. 2008; La Barbera et al. 2010; Magoulas et al. 2012; Cappellari et al. 2013a; Hou & Wang 2015; Singh et al. 2021; Howlett et al. 2022). However, Saglia et al. (2010), Saulder et al. (2019), and Singh et al. (2021) reported correlations between elliptical FP residuals and stellar mass/luminosity, which could partially account for a perceived environmental dependence.

In Horizon-AGN, we see no strong correlations between elliptical $\lambda$ and luminosity in Fig. 4. Further splitting of the elliptical galaxy samples in search of high-mass/luminosity signals is complicated by the already-small sample size – we are unable to successfully fit an FP to more than $\sim 6200$ ellipticals (Table 1, Section 3.2). We thus explore measurable correlations between the properties of both elliptical and spiral galaxies in our FP samples with estimates of the local galaxy density contrast.

## 5.2 Environmental dependence

We investigate the correlations of properties of our selected samples with their environments, using the local galaxy density contrast $\delta_{g,s}$, in spheres of radius $S h^{-1} M_{\odot}$, as a proxy. We note that the impact of tides imposed by the anisotropic cosmic web structure will not be well-characterized by this isotropic density metric. Future analyses of simulated ISCS should thus consider environmental geometry and top-down scale coupling in galaxy formation and evolution, which has been shown to influence galaxy properties that are related to sizes (Pichon et al. 2011; Codis et al. 2012; Codis, Pichon & Pogosyan 2015; Laigle et al. 2015; Kraljic et al. 2018; Musso et al. 2018).

We display property versus $\delta_{g,s}$ correlations in Fig. 7, with blue (red) curves denoting spiral (elliptical) samples, and solid (dashed) curves denoting the Masscut (Full) FPs. Galaxies are geometrically-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Galaxy two-point correlations, measured as functions of three-dimensional separations *r*, for Fundamental Plane (Section 3) samples selected from the Horizon-AGN simulation (Section 2). For spiral (top) and elliptical (bottom) galaxies, columns give the galaxy clustering $\xi_{gg}(r)$, the fractional difference in clustering $\Delta_{gg}(r)$ between intrinsically large/small galaxies, the density-intrinsic size cross-correlation $\xi_{g\lambda}(r)$, and the intrinsic size autocorrelation $\xi_{\lambda\lambda}(r)$, respectively. The clustering panels also display the large (small) galaxy clustering signals $\xi_{g\lambda}+g_{\lambda+}(\xi_{g\lambda}-g_{\lambda-})$ from which $\Delta_{gg}$ is derived (equation 12). Larger data points denote more conservative signals measured for Masscut Fundamental Plane samples ($M_* > 10^{10.5} M_{\odot}$; Section 2), whilst smaller points are measured using the Full Fundamental Plane samples featuring lower-mass galaxies. Pink solid (dotted) curves in the right-most panels give the results of fitting the 1-parameter phenomenological model (described in Section 4.1.2) separately to elliptical/spiral $\lambda\lambda$ correlations measured in the Masscut (Full) Fundamental Plane samples. Intrinsic size correlations are highly sensitive to galaxies' morphology and stellar mass range.</figcaption></figure>

poorer FP fits and induced asymmetry in $\lambda$ (Section 2.3); the steeply increasing stellar mass-function (as $M_*$ decreases from $10^{10} M_{\odot}$) of simulated ellipticals (Fig. 1), which are known to be overproduced at the low-mass end (Kaviraj et al. 2017); and the aforementioned concerns around low-mass clustering in Horizon-AGN, we opt to exclude these objects from our Full FP sample (and Lens samples; Table 1), imposing that elliptical stellar masses $M_* > 10^{10} M_{\odot}$.

For our fiducial choices, and upon inclusion of lower-mass ellipticals with $10^{10} < M_*/M_{\odot} < 10^{10.5}$, we are in agreement with the observational findings of Joachimi et al. (2015). That work studied fainter elliptical galaxies at low redshifts ($z < 0.2$), and made significant detections of negative $\Delta_{gg}$ and $g\lambda$, as well as a tentative detection of positive $\lambda\lambda$. We reproduce the negative $\Delta_{gg}$ detection at $\sim 5.6\sigma$ when considering the Full FP sample, and this is accompanied by weak hints of negative $\xi_{g\lambda}$ ($1.7\sigma$) and positive $\xi_{\lambda\lambda}$ ($1.4\sigma$), wherein scales are highly correlated. Thus we observe ellipticals with intrinsically small effective radii (Section 2.1) to be more clustered than their large counterparts, in opposition to the virial radius trend for spirals, and over a slightly broader range of scales. However, none of these signals are detected at $>0.7\sigma$ when the mass-selection $M_* > 10^{10.5} M_{\odot}$ is imposed.

Singh et al. (2021) detected positive $g\lambda$ correlations in higher-redshift ($0.16 < z < 0.7$) luminous red galaxy (LRG) data. However, their findings are compatible with ours and those of Joachimi et al. (2015), given that their sample is comparatively dominated by brighter, more massive objects; they were able to reproduce the negative signal seen by Joachimi et al. (2015) by making selections in galaxy luminosity and colour. Interestingly, our Masscut FP signals weaken towards zero as compared with the Full FP signals in each case. This is agreeable with the findings of Singh et al. (2021), since many of their low-luminosity selections – bringing their sample more into parity with our own – yield signals consistent with zero, or weakly negative ($\sim 2\sigma$) in the faintest cases.

Both Joachimi et al. (2015) and Singh et al. (2021) argued that the consistent assignment of central and satellite galaxies to $\lambda_+$ and $\lambda_-$, respectively, drove the $g\lambda$ trend; that is, elliptical FPs and their residuals are dependent upon galaxies' environments (Bernardi et al. 2003; D'Onofrio et al. 2008; La Barbera et al. 2010; Magoulas et al. 2012; Cappellari et al. 2013a; Hou & Wang 2015; Singh et al. 2021; Howlett et al. 2022). However, Saglia et al. (2010), Saulder et al. (2019), and Singh et al. (2021) reported correlations between elliptical FP residuals and stellar mass/luminosity, which could partially account for a perceived environmental dependence.

In Horizon-AGN, we see no strong correlations between elliptical $\lambda$ and luminosity in Fig. 4. Further splitting of the elliptical galaxy samples in search of high-mass/luminosity signals is complicated by the already-small sample size – we are unable to successfully fit an FP to more than $\sim 6200$ ellipticals (Table 1, Section 3.2). We thus explore measurable correlations between the properties of both elliptical and spiral galaxies in our FP samples with estimates of the local galaxy density contrast.

## 5.2 Environmental dependence

We investigate the correlations of properties of our selected samples with their environments, using the local galaxy density contrast $\delta_{g,s}$, in spheres of radius $S h^{-1} M_{\odot}$, as a proxy. We note that the impact of tides imposed by the anisotropic cosmic web structure will not be well-characterized by this isotropic density metric. Future analyses of simulated ISCS should thus consider environmental geometry and top-down scale coupling in galaxy formation and evolution, which has been shown to influence galaxy properties that are related to sizes (Pichon et al. 2011; Codis et al. 2012; Codis, Pichon & Pogosyan 2015; Laigle et al. 2015; Kraljic et al. 2018; Musso et al. 2018).

We display property versus $\delta_{g,s}$ correlations in Fig. 7, with blue (red) curves denoting spiral (elliptical) samples, and solid (dashed) curves denoting the Masscut (Full) FPs. Galaxies are geometrically-
<figure><img src="image_7.png" /><figcaption>**Figure 7.** Correlations between galaxy sample (Table 1) average properties (rows) and the local galaxy density contrast $\delta_{g,S}$ within spheres of radius $S$ (columns) – thus columns increase the smoothing scale for the local density estimate, from left to right. Solid (dashed) curves give Masscut (Full) Fundametal Plane sample objects, with spiral (elliptical) samples shown in blue (red) – see Table 1 for sample details. Excepting the top row, which gives the density-binned mean FP residual $\langle\lambda\rangle$, points give $\Delta X = \langle X\rangle_{\delta_{g,S}} / |\langle X\rangle| - 1$, the fractional deviation of the mean of quantity $X$ in a local density bin from the mean of $X$ over the entire type/mass-selected subsample. Errors give the standard error on the mean. Hatching indicates linear axis spacing in the range $-1 < \delta_{g,S} < 1$, with logarithmic spacing outside. Rows are, respectively, the Fundamental Plane radius residual $\lambda$, the measured log-effective radius, the velocity dispersion $\sigma$ for ellipticals, the surface brightness $I$ for ellipticals, the absolute magnitude, the circular velocity for spirals, the surface mass density $\Sigma$ for ellipticals, the stellar mass, the ellipticity $e$, the stellar mass-to-light ratio $M_*/L$, and the surface area (see Section 2 for property definitions). Elliptical and spiral galaxy samples exhibit distinct correlations of their properties with the local density, which are sensitive to the inclusion of lower-mass objects, and to the smoothing scale $S$ applied to the density field.</figcaption></figure>

binned according to the number density (computed using all galaxies from Table 1) of spheres centred on their locations, and then normalized by the mean density in the box to give the galaxy density contrast $\delta_{g,S}$. Thus from left to right in each panel, curves report statistics for

low- to high-density regions, and the columns increase the smoothing scale from left to right. Points then give the fractional deviation of the density-binned mean property, relative to the unbinned mean property, with error bars equal to the appropriately propagated
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Correlations between galaxy sample (Table 1) average properties (rows) and the local galaxy density contrast $\delta_{g,S}$ within spheres of radius $S$ (columns) – thus columns increase the smoothing scale for the local density estimate, from left to right. Solid (dashed) curves give Masscut (Full) Fundametal Plane sample objects, with spiral (elliptical) samples shown in blue (red) – see Table 1 for sample details. Excepting the top row, which gives the density-binned mean FP residual $\langle\lambda\rangle$, points give $\Delta X = \langle X\rangle_{\delta_{g,S}} / |\langle X\rangle| - 1$, the fractional deviation of the mean of quantity $X$ in a local density bin from the mean of $X$ over the entire type/mass-selected subsample. Errors give the standard error on the mean. Hatching indicates linear axis spacing in the range $-1 < \delta_{g,S} < 1$, with logarithmic spacing outside. Rows are, respectively, the Fundamental Plane radius residual $\lambda$, the measured log-effective radius, the velocity dispersion $\sigma$ for ellipticals, the surface brightness $I$ for ellipticals, the absolute magnitude, the circular velocity for spirals, the surface mass density $\Sigma$ for ellipticals, the stellar mass, the ellipticity $e$, the stellar mass-to-light ratio $M_*/L$, and the surface area (see Section 2 for property definitions). Elliptical and spiral galaxy samples exhibit distinct correlations of their properties with the local density, which are sensitive to the inclusion of lower-mass objects, and to the smoothing scale $S$ applied to the density field.</figcaption></figure>

binned according to the number density (computed using all galaxies from Table 1) of spheres centred on their locations, and then normalized by the mean density in the box to give the galaxy density contrast $\delta_{g,S}$. Thus from left to right in each panel, curves report statistics for

low- to high-density regions, and the columns increase the smoothing scale from left to right. Points then give the fractional deviation of the density-binned mean property, relative to the unbinned mean property, with error bars equal to the appropriately propagated
standard error on the mean. The x-axes are log-spaced, except where hatching indicates linear spacing between $-1 < \delta_{g,S} < 1$.

The observed relations between $\langle\lambda\rangle$ and $\delta_{g,S}$ (Fig. 7, top row) are consistent with our measured ISCs; high-density regions host spirals (ellipticals) which tend to have positive (negative) $\lambda$, hence we see opposite signs for spiral and elliptical $\Delta_{gg}$, $\xi_{g\lambda}$ in Fig. 6. The deprecation of signals on application of the mass-selection is also explained, as we see that the inclusion of lower-mass objects creates or steepens density-$\langle\lambda\rangle$ correlations in most panels (dashed curves). We start by discussing elliptical galaxies, which are more sensitive to the mass-selection at the level of ISCS.

### 5.2.1 Elliptical environmental dependence

For elliptical galaxies, correlations are pronounced between $\delta_{g,S}$ and effective radius $R_{\text{eff}}$, velocity dispersion $\sigma$, virial radius $R_{\text{vir*}}$, absolute magnitude $M_r$, stellar mass $M_*$, and to slightly lesser extents, surface brightness $I$, and ellipticity $e$, suggesting systematic trends in the early-type galaxy population between variably-dense environments in the simulation.

The qualitative forms of $\delta_{g,S} - \{R_{\text{eff}}, \sigma, I, M_*\}$ are each in agreement with the two-point correlations measured by Singh et al. (2021). We note that the surface brightness $I$ correlates negatively with the local density, whilst the surface stellar mass density relation $\delta_{g,S} - \Sigma$ is noisier, displaying no clear correlation. This is suggestive of variable elliptical mass-to-light ratios across different environments (La Barbera et al. 2010; Cappellari et al. 2013a,b; Suess et al. 2019a,b), and may be related to the erasure of ISCs measured on the $\Sigma$-FP. Indeed, we see positive correlations of the stellar mass-to-light ratio $M_*/L$ with local density, in agreement with La Barbera et al. (2010) (though we note that our stellar masses and luminosities are computed for all stellar particles assigned to a galaxy, and not at any principled radius; Cappellari et al. 2013a). Coupled with the positive $\delta_{g,S}-M_*$ correlation, we infer that the surface area of our ellipticals generally increases with local density, and this is confirmed by the measured $\delta_{g,S}$-surface area correlation (bottom row).

However, the relaxation of the mass-selection $M_* > 10^{10.5} M_\odot$ in the elliptical FP sample can be seen in Fig. 7 to slash many of these trends (solid versus dashed red curves). The additional lower-mass galaxies primarily influence the mean properties of high-density environments by dragging them back towards the global means (given by zero-lines in the figure), with the notable exceptions of $\langle\lambda\rangle$, for which a negative $\langle\lambda\rangle$ at high-$\delta_{g,S}$ is induced by the lower-mass objects, and $M_*/L$ and $I$, which remain consistent.

Thus, relative to the Masscut sample, regions of high density in the Full FP sample are seen to contain more compact, fainter, lower stellar mass, lower velocity dispersion elliptical galaxies (as seen by Poggianti et al. 2013; Cappellari et al. 2013a; Baldry et al. 2020). These changes in density-mean property gradients could be sourced by the addition of e.g. larger ellipticals to intermediate-density environments, but given that the lower-mass galaxies are known to be of lower absolute radius (Fig. 1), we can infer that they are indeed preferentially located in high-density environments. The expanded population of ellipticals in dense environments then induces the gradient in $\delta_{g,S} - \langle\lambda\rangle$, which manifests as non-zero, negative $\Delta_{gg}$ and $\xi_{g\lambda}$ signals for the Full FP sample, where the former is statistically significant.

### 5.2.2 Spiral environment dependence

Our Masscut spiral FP samples show positive correlations between the local density and the circular velocity $V_c$, the stellar mass $M_*$,

and the stellar mass-to-light ratio $M_*/L$, which greatly strengthen upon relaxation of the mass-selection (dashed versus solid blue lines). Correlations between $\delta_{g,S}$ and effective/virial radii (absolute magnitudes) meanwhile are weak for the Masscut sample, and strongly positive (negative) for the Full sample, suggesting that the Masscut sample is more homogenized, in terms of absolute size, across density regimes.

The lack of environmental variation in the radii of high-mass spirals $M_* > 10^{10.5} M_\odot$ is in agreement with the observational findings of Maltby et al. (2010) and Lani et al. (2013), who saw no strong evidence for any such scaling. However, Maltby et al. (2010), Cebrián & Trujillo (2014), and Matharu et al. (2019) also saw in observations that lower-mass spirals tended to be larger in the field than in dense environments, which may disagree with what we see when relaxing the mass-selection in our simulated samples (dashed blue lines). Whilst the correlation of absolute size with density may be explained by the matching stellar mass-density correlation, a weaker $\langle\lambda\rangle - \delta_{g,S}$ correlation persists, suggesting that (modulo $\langle\lambda\rangle - R_{\text{eff/vir*}}$ correlations; Fig. 4) spirals are larger in high-density environments at fixed stellar mass (Fig. 5). Trends showing increased luminosity at higher densities for the Full sample are agreeable with some studies of the environmental dependence of the late-type luminosity function (Mo et al. 2004; Croton et al. 2005; Zucca et al. 2009; Eardley et al. 2015), whilst the lack of correlation seen for the Masscut sample agrees with others (Tempel et al. 2011; Zandivarez & Martínez 2011; Poudel et al. 2016). More detailed investigations of these differences are beyond the scope of this work.

Noting that the radii of lower-mass spirals are smaller than those existing only in the Masscut sample (Fig. 1; low-mass also have also lower $V_c$, not shown), the changing density-radius trends must be driven by the addition of small-radius objects to the field. This offers some explanation for the loss of signal significance upon application of the mass-selection, seen in Section 5.1.1; whilst the $\delta_{g,S}-\langle\lambda\rangle$ trend remains for the Masscut spiral sample, the measured correlations $\Delta_{gg}$ and $\xi_{g\lambda}$ are rendered statistically insignificant by the loss of signal-to-noise due to the removal of ~75 per cent (Table 1) of the source objects, preferentially from the low-density field. Meanwhile the loss of signal-to-noise is less severe for the $\xi_{g\lambda}$ correlation, which relies less upon a broad sampling of environments, and is thus maintained at $>3\sigma$ significance through the mass-selection (though there are some indications that the Masscut spiral $\xi_{g\lambda}$ may be contaminated by $\lambda$-property correlations; see Appendix A).

## 5.3 Intrinsic size correlations summarized

We further elucidate the intrinsic/absolute size-density trends by considering the FPs again, now coloured by the local density estimates from Section 5.2, shown in Fig. 8. For the smoothing scale $S = 1 h^{-1}$ Mpc, one clearly sees gradients in $\delta_{g,S}$ in each panel, in directions almost orthogonal to the 1:1 relation, here given in black, with $\lambda_-$, $\lambda_+$ annotations to denote the intrinsically small and large radii. The dashed black lines and contours illustrate the tilts of the planes, and the distribution of lower-mass objects removed from the Masscut FPs, respectively (as described in Section 3.2).

The tendency of absolutely small, intrinsically small (large) ellipticals to inhabit more (less) dense environments is seen in both of the left-hand panels, and is more pronounced when the smoothing scale is reduced to $S = 0.5 h^{-1}$ Mpc, but quickly weakens for $S \ge 2 h^{-1}$ Mpc (not shown) – signifying that more localized structures are the drivers of this environmental dependence.

Similarly, one begins to see hints of an opposite gradient in the top-right-hand panel, for Full FP spirals, and even more so in the
standard error on the mean. The x-axes are log-spaced, except where
hatching indicates linear spacing between −1 < δg,S < 1.

The observed relations between ⟨λ⟩ and δg,S (Fig. 7, top row) are consistent with our measured ISCs; high-density regions host spirals (ellipticals) which tend to have positive (negative) λ, hence we see opposite signs for spiral and elliptical Δgg, ξgλ in Fig. 6. The deprecation of signals on application of the mass-selection is also explained, as we see that the inclusion of lower-mass objects creates or steepens density-⟨λ⟩ correlations in most panels (dashed curves). We start by discussing elliptical galaxies, which are more sensitive to the mass-selection at the level of ISCS.

5.2.1 Elliptical environmental dependence

For elliptical galaxies, correlations are pronounced between δg,S and effective radius R<sub>eff</sub>, velocity dispersion σ, virial radius R<sub>vir*</sub>, absolute magnitude M<sub>r</sub>, stellar mass M*, and to slightly lesser extents, surface brightness I, and ellipticity *e*, suggesting systematic trends in the early-type galaxy population between variably-dense environments in the simulation.

The qualitative forms of δ<sub>g,S</sub> − {R<sub>eff</sub>, σ, I, M<sub>r</sub>} are each in
agreement with the two-point correlations measured by Singh et al.
(2021). We note that the surface brightness *I* correlates negatively
with the local density, whilst the surface stellar mass density relation
δ<sub>g,S</sub>Σ is noisier, displaying no clear correlation. This is suggestive of
variable elliptical mass-to-light ratios across different environments
(La Barbera et al. 2010; Cappellari et al. 2013a,b; Suess et al.
2019a,b), and may be related to the erasure of ISCs measured on
the Σ-FP. Indeed, we see positive correlations of the stellar mass-to-
light ratio M<sub>*</sub>/L with local density, in agreement with La Barbera et al.
(2010) (though we note that our stellar masses and luminosities are
computed for all stellar particles assigned to a galaxy, and not at any
principled radius; Cappellari et al. 2013a). Coupled with the positive
δ<sub>g,S</sub>–M<sub>*</sub> correlation, we infer that the surface area of our ellipticals
generally increases with local density, and this is confirmed by the
measured δ<sub>g,S</sub>–surface area correlation (bottom row).

However, the relaxation of the mass-selection $M_{*} > 10^{10.5} M_{\odot}$ in the elliptical FP sample can be seen in Fig. 7 to slash many of these trends (solid versus dashed red curves). The additional lower-mass galaxies primarily influence the mean properties of high-density environments by dragging them back towards the global means (given by zero-lines in the figure), with the notable exceptions of $\langle\lambda\rangle$, for which a negative $\langle\lambda\rangle$ at high-$\delta_{g,S}$ is induced by the lower-mass objects, and $M_{*}/L$ and $I$, which remain consistent.

Thus, relative to the Masscut sample, regions of high density in the Full FP sample are seen to contain more compact, fainter, lower stellar mass, lower velocity dispersion elliptical galaxies (as seen by Poggianti et al. 2013; Cappellari et al. 2013a; Baldry et al. 2020). These changes in density-mean property gradients could be sourced by the addition of e.g. larger ellipticals to intermediate-density environments, but given that the lower-mass galaxies are known to be of lower absolute radius (Fig. 1), we can infer that they are indeed preferentially located in high-density environments. The expanded population of ellipticals in dense environments then induces the gradient in δ<sub>g,S</sub> − ⟨λ⟩, which manifests as non-zero, negative Δ<sub>gg</sub> and ξ<sub>gλ</sub> signals for the Full FP sample, where the former is statistically significant.

5.2.2 Spiral environment dependence

Our Masscut spiral FP samples show positive correlations between
the local density and the circular velocity $V_c$, the stellar mass $M_*$,

and the stellar mass-to-light ratio M*/L, which greatly strengthen
upon relaxation of the mass-selection (dashed versus solid blue
lines). Correlations between δg,S and effective/virial radii (absolute
magnitudes) meanwhile are weak for the Masscut sample, and
strongly positive (negative) for the Full sample, suggesting that the
Masscut sample is more homogenized, in terms of absolute size,
across density regimes.

The lack of environmental variation in the radii of high-mass
spirals M* > 10<sup>10.5</sup> M<sub>☉</sub> is in agreement with the observational
findings of Maltby et al. (2010) and Lani et al. (2013), who saw
no strong evidence for any such scaling. However, Maltby et al.
(2010), Cebrián & Trujillo (2014), and Matharu et al. (2019) also
saw in observations that lower-mass spirals tended to be larger
in the field than in dense environments, which may disagree with
what we see when relaxing the mass-selection in our simulated
samples (dashed blue lines). Whilst the correlation of absolute size
with density may be explained by the matching stellar mass-density
correlation, a weaker ⟨λ⟩ − δg,S correlation persists, suggesting
that (modulo ⟨λ⟩ − R<sub>eff/vir*</sub> correlations; Fig. 4) spirals are larger
in high-density environments at fixed stellar mass (Fig. 5). Trends
showing increased luminosity at higher densities for the Full sample
are agreeable with some studies of the environmental dependence of
the late-type luminosity function (Mo et al. 2004; Croton et al. 2005;
Zucca et al. 2009; Eardley et al. 2015), whilst the lack of correlation
seen for the Masscut sample agrees with others (Tempel et al. 2011;
Zandivarez & Martínez 2011; Poudel et al. 2016). More detailed
investigations of these differences are beyond the scope of this work.

Noting that the radii of lower-mass spirals are smaller than those
existing only in the Masscut sample (Fig. 1; low-mass also have also
lower $V_c$, not shown), the changing density-radius trends must be
driven by the addition of small-radius objects to the field. This offers
some explanation for the loss of signal significance upon application
of the mass-selection, seen in Section 5.1.1; whilst the $\delta_{g,S}-\langle\lambda\rangle$ trend
remains for the Masscut spiral sample, the measured correlations $\Delta_{gg}$
and $\xi_{g\lambda}$ are rendered statistically insignificant by the loss of signal-
to-noise due to the removal of ~75 per cent (Table 1) of the source
objects, preferentially from the low-density field. Meanwhile the loss
of signal-to-noise is less severe for the $\xi_{g\lambda}$ correlation, which relies
less upon a broad sampling of environments, and is thus maintained
at >3$\sigma$ significance through the mass-selection (though there are
some indications that the Masscut spiral $\xi_{g\lambda}$ may be contaminated
by $\lambda$-property correlations; see Appendix A).

**5.3 Intrinsic size correlations summarized**

We further elucidate the intrinsic/absolute size-density trends by
considering the FPs again, now coloured by the local density
estimates from Section 5.2, shown in Fig. 8. For the smoothing scale
$S = 1 h^{-1}$ Mpc, one clearly sees gradients in $\delta_{g,S}$ in each panel, in
directions almost orthogonal to the 1:1 relation, here given in black,
with $\lambda_-$, $\lambda_+$ annotations to denote the intrinsically small and large
radii. The dashed black lines and contours illustrate the tilts of the
planes, and the distribution of lower-mass objects removed from the
Masscut FPs, respectively (as described in Section 3.2).

The tendency of absolutely small, intrinsically small (large) ellipticals to inhabit more (less) dense environments is seen in both of the left-hand panels, and is more pronounced when the smoothing scale is reduced to $S = 0.5 h^{-1}$ Mpc, but quickly weakens for $S \ge 2 h^{-1}$ Mpc (not shown) – signifying that more localized structures are the drivers of this environmental dependence.

Similarly, one begins to see hints of an opposite gradient in the
top-right-hand panel, for Full FP spirals, and even more so in the
<figure><img src="image_8.png" /><figcaption>**Figure 8.** Elliptical (**left**) and spiral (**right**) Fundamental Planes, as in Fig. 3, but coloured here by the galaxy number density contrast $\delta_{g,S}$ estimated for galaxies in spheres of radius $S = 1\ h^{-1}\ \text{Mpc}$ (Section 5.2). The 1:1 relation, the linear fits to the points on each plane, and the contours showing the lower-mass objects' distributions are the same as those in Fig. 3, now given in black. Annotations $\lambda_-$ and $\lambda_+$ illustrate the sides of the plane corresponding to intrinsically small and large galaxies, respectively. Clear and opposite gradients in local density are seen across the Fundamental Planes of spiral and elliptical galaxies, with or without the application of the conservative mass selection (**bottom**; see Section 2).</figcaption></figure>

Masscut case (bottom-right). For spirals, the gradients are weaker for $S = 0.5\ h^{-1}\ \text{Mpc}$, but stronger for $S = 2, 4\ h^{-1}\ \text{Mpc}$, and still visible for $S = 8\ h^{-1}\ \text{Mpc}$; the structures driving this correlation of local density with FP residuals are thus more extended than those driving the elliptical correlations.

The evolving galaxy size-mass-morphology distribution is likely to have a complex dependence upon the anisotropic cosmic web environment of voids, walls, filaments, saddle points, and nodes, where the relative efficiency and vorticity of gas flows dictate the build-up of galactic angular momenta (Pichon et al. 2011; Codis et al. 2012; Welker et al. 2014; Codis et al. 2015; Laigle et al. 2015); the bulk velocity flows transport galaxies through the different environments over cosmic time (Codis et al. 2015; Kraljic et al. 2018; Laigle et al. 2018); and the variable efficiency of mass accretion and mergers strongly influence galaxies' formation and evolution (Welker et al. 2017; Musso et al. 2018). Each such process has the potential to force galaxies out of kinetic equilibrium, and/or influence the new equilibrium to be reached after relaxation; indeed the intrinsic size distribution may be dependent upon the balance of different kinetic equilibria and dis-equilibria.

We defer a detailed analysis of the intrinsic size distribution across anisotropic cosmic environments to future work. For now we turn to Welker et al. (2017), who studied the impacts of different stellar growth processes upon galaxy size and morphological evolution over $z = 1–5$ in Horizon-AGN. Whilst we are extrapolating their findings from $z = 1–5$ down to $z = 0$ for our interpretations here, Welker et al. (2017) found mergers of any kind to be rare, violent events, occurring on average only twice in the history of a galaxy over the epoch $z = 1–5$. Given that the merger rate decreases with cosmic time (Rodriguez-Gomez et al. 2015), we assume that their conclusions are largely applicable to our work here.

They found that multiple minor, gas-deprived mergers had similar effects to major mergers in terms of destroying disc structures and forming spheroids. Minor, gas-rich mergers were even seen to cause contraction of the remnant objects (compatibly with the gas compaction paradigm of Dekel & Burkert 2014; Inoue et al. 2016). Thus wet mergers are a possible pathway for the compact elliptical galaxies driving our measured correlations, provided that they are not consumed by later mergers. One does see that the absolutely-and-intrinsically small galaxies constitute a minority of the supplementary objects in the Full elliptical FP, as they sit
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Elliptical (**left**) and spiral (**right**) Fundamental Planes, as in Fig. 3, but coloured here by the galaxy number density contrast $\delta_{g,S}$ estimated for galaxies in spheres of radius $S = 1\ h^{-1}\ \text{Mpc}$ (Section 5.2). The 1:1 relation, the linear fits to the points on each plane, and the contours showing the lower-mass objects' distributions are the same as those in Fig. 3, now given in black. Annotations $\lambda_-$ and $\lambda_+$ illustrate the sides of the plane corresponding to intrinsically small and large galaxies, respectively. Clear and opposite gradients in local density are seen across the Fundamental Planes of spiral and elliptical galaxies, with or without the application of the conservative mass selection (**bottom**; see Section 2).</figcaption></figure>

Masscut case (bottom-right). For spirals, the gradients are weaker for $S = 0.5\ h^{-1}\ \text{Mpc}$, but stronger for $S = 2, 4\ h^{-1}\ \text{Mpc}$, and still visible for $S = 8\ h^{-1}\ \text{Mpc}$; the structures driving this correlation of local density with FP residuals are thus more extended than those driving the elliptical correlations.

The evolving galaxy size-mass-morphology distribution is likely to have a complex dependence upon the anisotropic cosmic web environment of voids, walls, filaments, saddle points, and nodes, where the relative efficiency and vorticity of gas flows dictate the build-up of galactic angular momenta (Pichon et al. 2011; Codis et al. 2012; Welker et al. 2014; Codis et al. 2015; Laigle et al. 2015); the bulk velocity flows transport galaxies through the different environments over cosmic time (Codis et al. 2015; Kraljic et al. 2018; Laigle et al. 2018); and the variable efficiency of mass accretion and mergers strongly influence galaxies' formation and evolution (Welker et al. 2017; Musso et al. 2018). Each such process has the potential to force galaxies out of kinetic equilibrium, and/or influence the new equilibrium to be reached after relaxation; indeed the intrinsic size distribution may be dependent upon the balance of different kinetic equilibria and dis-equilibria.

We defer a detailed analysis of the intrinsic size distribution across anisotropic cosmic environments to future work. For now we turn to Welker et al. (2017), who studied the impacts of different stellar growth processes upon galaxy size and morphological evolution over $z = 1–5$ in Horizon-AGN. Whilst we are extrapolating their findings from $z = 1–5$ down to $z = 0$ for our interpretations here, Welker et al. (2017) found mergers of any kind to be rare, violent events, occurring on average only twice in the history of a galaxy over the epoch $z = 1–5$. Given that the merger rate decreases with cosmic time (Rodriguez-Gomez et al. 2015), we assume that their conclusions are largely applicable to our work here.

They found that multiple minor, gas-deprived mergers had similar effects to major mergers in terms of destroying disc structures and forming spheroids. Minor, gas-rich mergers were even seen to cause contraction of the remnant objects (compatibly with the gas compaction paradigm of Dekel & Burkert 2014; Inoue et al. 2016). Thus wet mergers are a possible pathway for the compact elliptical galaxies driving our measured correlations, provided that they are not consumed by later mergers. One does see that the absolutely-and-intrinsically small galaxies constitute a minority of the supplementary objects in the Full elliptical FP, as they sit