# Impact of the environment on the angular momentum of
z ~ 0.25 – 1.5 galaxies in the MAGIC survey*

W. Mercier<fnref n="1" />,**, B. Epinat<fnref n="2" />,<fnref n="3" />, T. Contini<fnref n="1" />, and MUSE collaboration

Received ; accepted

## ABSTRACT

**Aims.** Group and cluster environments are thought to impact galaxy properties at intermediate redshift. Among these properties, the angular momentum is a useful proxy to trace environmental effects when the galaxies are infalling into their host structure.

**Methods.** We derive robust estimates of the angular momentum using Hubble Space Telescope images combined with spatially resolved kinematics from the Multi-Unit Spectroscopic Explorer for a sample of roughly 200 galaxies in groups, clusters, and in the field drawn from the MAGIC survey. Using various environmental tracers, we study the shape of the angular momentum-stellar mass (Fall) relation as a function of environment. We also study the impact of bulges on the angular momentum estimate.

**Results.** We find a significant impact of the environment on the zero-point of the Fall relation (0.12 dex without cut and 0.08 dex with mass and redshift cuts) consistent with galaxies in structures being depleted in angular momentum at fixed stellar mass. This effect seems driven by massive galaxies that are found in the densest regions of the structures and does not seem to correlate with their bulge fraction.

**Conclusions.** We find that galaxy groups and clusters at z ~ 0.7 seem on average to deplete their galaxy member's angular momentum with respect to the field. However, the lack of a current comparable sample of massive field galaxies prevent us from drawing clear conclusions. At the moment, we cannot eliminate the possibility that the observed zero point offsets are not mass-driven.

**Key words.** Galaxies: evolution – Galaxies: kinematics and dynamics – Galaxies: clusters: general – Galaxies: groups – Galaxies: high-redshift

## 1. Introduction

In the current paradigm of galaxy evolution, galaxies are expected to form through the condensation of baryons in the centres of DM haloes where, because of external tidal torques, the gas in the proto-galaxy acquires angular momentum before condensing into a disk and forming stars (e.g. Peebles 1969; Fall & Efstathiou 1980). The same is true for Dark Matter (DM) haloes that also acquire angular momentum during their linear phase of structure growth (e.g. Doroshkevich 1970; White 1984). Initially, it was thought that the angular momentum of the baryons (e.g. stars or gas) traces that of the DM component. However, recent simulations have shown that this picture is not entirely correct with processes that can either add or remove angular momentum independently of the halo and proto-galaxy early formation phase (e.g. Genel et al. 2015). For instance, there is strong evidence supporting the fact that galaxies must have smoothly accreted large amounts of cold gas from their circumgalactic medium (CGM) to sustain high values of star formation rate (SFR) across cosmic time (e.g. Bouché et al. 2013, 2016; Zabl et al. 2019). This accretion of fresh gas takes place predominantly in the disk plane of late-type galaxies and thus can

not only drive their star formation history but also the angular momentum growth of their baryonic component (e.g. Danovich et al. 2015; Cadiou et al. 2022). Similarly, feedback processes such as galactic winds have also been found as potential mechanisms to increase the angular momentum of the baryons in a galaxy (e.g. DeFelippis et al. 2017). In addition, galaxy mergers can also substantially redistribute angular momentum either by increasing that of the DM halo (e.g. Hetznecker & Burkert 2006), or by redistributing the mass and thus either increasing the angular momentum of the baryons (e.g. Genel et al. 2015; Jiang et al. 2019) or increasing it if the spins of the two galaxies are aligned prior to the merger (e.g. Lagos et al. 2018).

Numerous studies have tried to constrain the angular momentum of the stellar and gas components in galaxies at various redshifts since the seminal work of Fall (1983) by studying the angular momentum-stellar mass relation, also known as the Fall relation. In the local Universe, some of these studies have focussed on the shape of the relation as a function of galaxy type (e.g. Fall 1983; Romanowsky & Fall 2012; Cortese et al. 2016; Rizzo et al. 2018), stellar mass (e.g. Posti et al. 2018; Di Teodoro et al. 2021), cold gas fraction (e.g. Mancera Piña et al. 2021b), or bulge fraction (e.g. Obreschkow & Glazebrook 2014; Fall & Romanowsky 2018). The general conclusion that can be drawn is that, in the local Universe, galaxies populate the Fall relation that can be described as a single linear relation in log-log space (i.e. no deviation at higher or lower stellar masses). The scatter of the Fall relation seems to be mainly driven by their

<footnote n="1">¹ Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 31400 Toulouse, France</footnote>

<footnote n="2">² Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="3">³ Canada-France-Hawaii Telescope, CNRS, 96743 Kamuela, Hawaii, USA</footnote>

<footnote>* Based on observations made with ESO telescopes at the Paranal Observatory under programs 094.A-0247, 095.A-0118, 096.A-0596, 097.A-0254, 099.A-0246, 100.A-0607, 101.A-0282, 102.A-0327, and 103.A-0563.</footnote>

<footnote>** e-mail: wilfried.mercier@irap.omp.eu</footnote>
morphology, with ellipticals and disk galaxies with prominent bulges located below the relation followed by disk galaxies. The Fall relation has also been studied at higher redshift ($z \sim 1-2$, e.g. Burkert et al. 2016; Swinbank et al. 2017; Harrison et al. 2017; Bouché et al. 2021) where it is found to hold but with a lower zero-point that is consistent with galaxies gaining angular momentum with cosmic time. Interestingly, a recent analysis by Bouché et al. (2021) pointed out that there is a correlation between the zero-point of the Fall relation of the ionised gas and the dynamical state of the galaxies in the sense that dispersion dominated systems are located below the Fall relation of rotationally supported galaxies. This result was actually observed on a few other studies before (e.g. Contini et al. 2016; Burkert et al. 2016) and is consistent with the separation observed in the local Universe between spirals and ellipticals (e.g. Romanowsky & Fall 2012).

From a methodological perspective, simulations put aside, measuring the angular momentum in galaxies is not simple since it requires to know the full 3D mass distribution (stars, gas, and DM) as well as the amplitude and orientation of the velocity vector at each position in the galaxy. Hence, given the impossibility to derive the "true" angular momentum, proxies have been used in the literature that rely on various assumptions. The most common is that the gas and the stars are located in a disk that is dynamically stable because of its own rotation. With a few additional assumptions on the stars or gas kinematics and on the mass distribution in the disk, it is possible to derive the Romanowsky & Fall (2012, hereafter RF12) approximation that is widely used in the literature (e.g. Romanowsky & Fall 2012; Contini et al. 2016; Burkert et al. 2016; Swinbank et al. 2017; Rizzo et al. 2018). Though already mentioned in Romanowsky & Fall (2012) and in subsequent studies, a recent discussion was given in Bouché et al. (2021) where it is shown that, depending on the location where the angular momentum is measured, the RF12 estimator can overestimate the angular momentum by nearly 20%. Hence, because of that effect and of the increasing number of robust kinematics measurements throughout the last ten years, authors have started using more general estimates. Some have chosen to consider an axisymmetric disk model (typically exponential) and either to numerically integrate the expression of the angular momentum for a rotating disk galaxy (e.g. Posti et al. 2018; Bouché et al. 2021; Mancera Piña et al. 2021a,b) or by substituting the integral with a sum (e.g. Rizzo et al. 2018). Cortese et al. (2016) also proposed an interesting method slightly different from those previously mentioned in the sense that the numerical integration of the angular momentum is performed by summing the contribution of each spatial pixel (spaxel) in a data cube obtained through 3D spectroscopy. The advantage of such a method is that it alleviates the assumption of axial symmetry since the observed flux distribution is directly taken into account. However, its main drawback is that it suffers from the usually poor spatial resolution of 3D spectroscopic observations compared to high spatial resolution images obtained for instance by the Hubble Space Telescope (HST).

Finally, a few studies have tried to probe the impact of the environment on the galaxies' angular momentum (e.g. Pelliccia et al. 2019; Pérez-Martínez et al. 2021). The current picture that seems to emerge is that galaxies found in galaxy groups and galaxy clusters at $z \sim 1$ seem to have a deficit of angular momentum with respect to field galaxies located at the same redshift. Because physical mechanisms can affect the angular momentum of the baryons in different ways (i.e. increase or decrease it), an interpretation given in Pelliccia et al. (2019) is that this reduction of angular momentum at fixed stellar mass

could be due to galaxy mergers, in line with their prevalence in dense environments and with recent simulations (e.g. Lagos et al. 2018). However, these results were obtained on moderately large galaxy samples and they required to compare galaxies from different surveys/instruments (e.g. Pérez-Martínez et al. 2021) to reach these conclusions. As in Abril-Melgarejo et al. (2021) and Mercier et al. (2022), we argue that there might be systematic effects when doing so, mainly driven by different selection functions between different datasets. In this analysis, we propose to study the impact of the environment on the stellar angular momentum of galaxies in the MUSE-gAlaxy Groups In COSMOS (MAGIC) survey (Epinat et al., in prep.). As shown in Mercier et al. (2022), this survey is ideal to probe the environmental impact on galaxy dynamics at $z \sim 1$ because of the fact that we observe at the same time galaxies located in structures of varying density (i.e. galaxy groups and galaxy clusters) and foreground and background galaxies located in the field.

This paper is structured as follows. In Sect. 2, we give a brief description of the MAGIC sample, the HST and MUSE observations, how we performed the morphological and kinematics modellings, and how we characterised the galaxies' environment. More importantly, we also describe the sample selection used to study the angular momentum. In Sect. 3, we describe the method used to derive the angular momentum using HST images. Afterwards, we perform the analysis of the Fall relation in Sect. 4, where we first assess the reliability of our method in Sect. 4.1, then we study the impact of the environment without taking into account the contribution of bulges in Sect. 4.2, and finally we discuss the impact of bulges on the shape of the Fall relation and their link to the galaxies' environment in Sect. 4.3. Throughout this paper, we assume a Λ cold dark matter cosmology with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_M = 0.3$, and $\Omega_\Lambda = 0.7$.

## 2. Sample selection and main properties

The sample used in this analysis is part of the MAGIC survey (Epinat et al., in prep.), a deep Multi-Unit Spectroscopic Explorer (MUSE) Guaranteed Time Observation (GTO) survey targeting 15 structures in the Cosmic Evolution Survey (Cosmos) area (Scoville et al. 2007) with 17 different MUSE pointings. The main goal is to study the impact of the environment on the properties of galaxies at intermediate redshift by combining multi-band photometry, in particular Hubble Space Telescope Advanced Camera for Surveys (HST-ACS) observations (Koekemoer et al. 2007; Massey et al. 2010), with resolved spectroscopic properties from MUSE. The analysis performed in this paper is the continuation of two previous ones. The first was on the impact of the environment on the Tully-Fisher Relation (TFR) for a subsample of group galaxies (Abril-Melgarejo et al. 2021). The second focussed on the impact of the environment on various scaling relations such as the size-mass relation, the Main Sequence relation (MS), and the TFR for the full MAGIC sample by comparing between galaxies located in the field and those found in small and large groups (Mercier et al. 2022). We summarise the main observations and their characteristics below and we refer the reader to Epinat et al. (in prep.) for a complete description of the MAGIC survey.

### 2.1. Observations and physical parameters

The basis of the observations are 17 different MUSE data and variance cubes with a spatial sampling of 0.2" and a spectral sampling of 1.25 Å over a spectral range ranging from 4750 Å
to 9350 Å. Observing strategy and data reduction are fully described in Mercier et al. (2022) and Epinat et al. (in prep.). The MUSE Line Spread Function (LSF) was modelled with a second order polynomial function as in Bacon et al. (2017) and Guérou et al. (2017) and the MUSE Point Spread Function (PSF) was modelled by extracting 100 Å wide narrow-band images of stars in each MUSE field and by fitting them with a Moffat profile with free $\beta$ index. The wavelength dependence of the PSF Full Width at Half Maximum (FWHM) was then derived field-by-field by fitting a linear relation to the median curve and the $\beta$ parameter as the mean value of the various measurements weighted by their uncertainties.

The median value of the MUSE PSF FWHM for the 17 fields is 0.67" at 4000 Å and 0.53" at 8000 Å which correspond respectively to 4.8 kpc and 3.8 kpc at $z = 0.7$. Apart from the MUSE data we also have HST images that were used to model the galaxies' morphology. To that end, we used 4" × 4" stamps in the F814W filter since it provides the best spatial resolution possible (PSF FWHM below 0.1", pixel scale of 0.03"). Similarly to what was done with MUSE, the HST PSF was measured by fitting a Moffat profile onto 27 non saturated stars found in the MUSE fields and by deriving the median profile. The final parameters used for the morphological modelling are $FWHM_{HST} = 0.0852$" and $\beta = 1.9$.

Additionally, the two main physical parameters of interest for the previous analyses and for this one are the galaxies' stellar mass and star formation rate (SFR). They were derived from Spectral Energy Distribution (SED) fitting performed with the Code Investigating GALaxy Emission (CIGALE, see Boquien et al. 2019) using Bruzual & Charlot (2003) single stellar populations (SSPs) with a Salpeter (1955) initial mass function (IMF) and a single metallicity value of 0.02 dex, a truncated delayed exponential star formation history (SFH, described in Ciesla et al. 2018, 2021), and a Charlot & Fall (2000) attenuation law with a total to selective extinction ratio $R_V = 3.1$. More details about the grid of parameters and the choice of models will be given in Epinat et al. (in prep.). We note that these two parameters were derived differently in the previous papers. In Abril-Melgarejo et al. (2021), the stellar mass and SFR were derived from SED fitting performed with FAST. In Mercier et al. (2022), only the stellar mass from FAST was used and the SFR was derived from the Kennicutt (1992, 1998) law using the [O II] flux measured in the MUSE cubes in 3" apertures after correcting for Galactic and intrinsic extinction using the Gilbank et al. (2010, 2011) prescription. Overall, we find consistent stellar masses within 0.5 dex when above $10^8 M_\odot$ and slightly larger differences below this value. Additionally, CIGALE tends to find on average larger stellar masses than FAST by roughly 0.05 dex, probably because of differences in the IMF or the SFH used.

## 2.2. Dynamical modelling

A complete dynamical modelling for the entire [O II] emitters sample in the MAGIC survey was done in Mercier et al. (2022). In this analysis, we define as the [O II] emitters sample galaxies that are found in the range $0.2 < z < 1.5$. Furthermore, those for which the [O II] doublet is indeed detected in the MUSE data form the kinematics sample, since we can use their [O II] doublet to extract the ionised gas kinematics. In what follows, we quickly summarise the main steps and in particular we emphasize what has changed since the previous analysis performed in Mercier et al. (2022). A complete description for the morpholog-

ical and kinematics modelling can be found in Sect. 4.1 and 5.1 of Mercier et al. (2022).

Starting from an initial sample of 1142 [O II] emitters, the first step was to model the galaxies' morphology. To that hand, bulge-disk decompositions were performed with GALFIT (Peng et al. 2002) using an exponential disk model for the disk component and a circular de Vaucouleur profile for the bulge. Initially, among the 1142 [O II] emitters 890 could be modelled with GALFIT. For each galaxy various morphological parameters were derived: (i) the effective radius $R_{eff}^1$, (ii) the bulge-to-total flux ratio evaluated at one global effective radius (B/D), (iii) the disk major axis position angle (PA), and (iv) the apparent disk axis ratio $q = b/a$, with *a* and *b* the major and minor axes, respectively. In Mercier et al. (2022), we introduced two expressions to correct the central surface brightness of the disk, $\Sigma_{d,bos}(0)$, and the observed axis ratio $q$ of the galaxy for the finite thickness of the disk. The average redshift dependence of the disk thickness $q_0$ was evaluated thanks to the method presented in Heidmann et al. (1972) and Bottinelli et al. (1983). The central surface brightness of the disk was then corrected assuming a double exponential 3D mass distribution and the disk inclination was derived using the correction given in Bottinelli et al. (1983). In addition, because the rotation velocity that was used for the TFR in Abril-Melgarejo et al. (2021) and Mercier et al. (2022) was measured at the radius $R_{22} = 2.2R_d$ in the plane of the disk, with $R_d$ the disk scale length, we also introduced a corrected stellar mass ($M_{\star,corr}$) so that it is also evaluated at $R_{22}$. The correction assumes a constant mass-to-light ratio throughout the galaxy and takes into account the effect of sky projection and PSF convolution. We refer the reader to Sects. 4.3 and 4.4 and Appendix D.4, of Mercier et al. (2022) for more details. For this analysis and after inspecting the HST images, we decided to re-model the morphology for the 17 following galaxies: 85_CGR23, 82_CGR30, 97_CGR32, 35_CGR35, 74_CGR35, 127_CGR35, 94_CGR51, 36_CGR61, 39_CGR61, 52_CGR61, 121_CGR61, 41_CGR79, 23_CGR84, 122_CGR84, 79_CGR114, 31_CGR172, and 68_CGR172. Among these galaxies, we included a logarithmic spiral pattern to the disk component, available since version three of GALFIT, for the four following galaxies: 85_CGR23, 121_CGR61, 36_CGR61, and 79_CGR114. We did so because these galaxies showed strong spiral and bar features so that their previous disk morphological model used to fit the bar rather than the extended disk. Even if the modelling of the spiral arms and of the bar is not trivial and far from perfect, we did find more robust disk parameters once these features were included, especially for the disk axis ratio and PA. Overall, the remaining 13 galaxies had a center that was slightly misaligned with the bulge component which had the effect to mildly underestimate its contribution. With the updated models, we now have better constraints on the galaxy's centre and on the parameters of the bulge component. Furthermore, after trying to redo their morphology without more success than before, we decided to remove the 22 following galaxies: (i) 12_CGR61, 24_CGR32, 26_CGR61, 29_CGR23, 54_CGR79, 76_CGR30, 76_CGR84, 76_CGR172, 89_CGR79, 240_CGR30, and 240_CGR84 because they are strongly bulge-dominated (probably elliptical) whose disks are too faint to be efficiently constrained, (ii) 55_CGR61, 66_CGR30, 67_CGR87, 73_CGR61, 75_CGR84, 76_CGR172, 79_CGR61, 82_CGR87, and 97_CGR34 because their disks are detected but with too low S/N to be properly constrained, and (iii) 79_CGR30 and

<footnote n="1">¹ This corresponds to the global effective radius that takes into account both the disk and the bulge components given by Eq. 6 of Mercier et al. (2022).</footnote>
84_CGr35 because they both show complex edge-on morphologies that do not match their velocity fields. Indeed, 79_CGr30 looks like two edge-on galaxies whose major axes are aligned but with a velocity field that spans the two objects and that shows undisturbed rotation with a PA rotated by 90° with respect to the stellar components. Similarly, galaxy 84_CGr35 seems to show multiple diffuse components that might be the result of a past merger or a stripping event even though its velocity field also shows an extended rotation which is not aligned with the brightest component visible in the HST image. Therefore, the morphological and mass models for these two objects would be too uncertain to compute their angular momentum.

The second step was to extract the ionised gas kinematics using the [O II] doublet as kinematics tracer. This was performed with CAMEL<fnref n="2" /> by fitting the [O II] doublet spaxel per spaxel and then cleaning the maps to remove isolated spaxels and those with large velocity discontinuities. The detailed procedure of how this was done is described in Sect. 5.1 of Mercier et al. (2022). At the end of this step 271 galaxies were removed from the morphological sample because they had no remaining spaxels with a S/N in the [O II] doublet larger than five.

The last step of the modelling was to fit the velocity fields with MOCING<fnref n="3" /> using mass models as in Sect. 5.1 of Mercier et al. (2022) where it is discussed in detail. The main difference with what we did in our previous analysis is that we re-modelled the ionised gas kinematics for the entire kinematics sample using the updated Moffat MUSE PSF profiles in each field. As in Mercier et al. (2022), we used prior information from the morphological modelling to constrain the contribution of the stellar disk and bulge to the total rotation curve. The remaining component was modelled with a Navarro-Frenk-White (NFW) profile (Navarro et al. 1996). This last component mostly describes the Dark Matter (DM) content of the galaxies but, since we have no constraints on their amount of gas, technically speaking it also includes the contribution of the gas. Given a rotation curve, we used MOCING to model the ionised gas velocity field using the method of line moments (Epinat et al. 2010a) assuming the gas is located in a razor-thin disk. During the fitting process we fixed the parameters related to the disk and bulge models, and we fixed the centre and the inclination in order to remove degeneracies. Thus, the only free parameters are: (i) the DM halo parameters, its scale length $r_s$ and maximum velocity $V_{h,max}$, (ii) the kinematics position angle, and (iii) the systemic redshift $z_s$. Once the best-fit velocity field model was found, MOCING computed a correction map to account for the beam-smearing, an instrumental effect that increases the measured velocity dispersion where the velocity shear is large because of the smearing effect of the PSF, and removed it quadratically from the velocity dispersion map. At the end of this step, we removed four additional galaxies because they lied on the edge of the MUSE field-of-view or had signs of merger in their kinematics maps. Thus, we ended up with a kinematics sample of 571 galaxies.

## 2.3. Environment characterisation

The characterisation of the environment can be done in various ways, each method having its advantages and its drawbacks. For instance, the friends-of-friends (FoF) algorithm used in Knobel et al. (2012) and Iovino et al. (2016) is easy to implement but it only uses spectroscopic redshifts while the Voronoi tessellation Monte-Carlo mapping (VMC) technique discussed in

Lemaux et al. (2017, 2022) and Hung et al. (2020, 2021) combine both photometric and spectroscopic redshifts. Nevertheless, most techniques remain sensitive to the completeness of redshift measurements and the size of the studied fields.

The first step to characterise the galaxies’ environment is to define the structures they belong to. This was done using a 3D FoF (as in Mercier et al. 2022) with a maximum sky projected separation of 450 kpc and a maximum line-of-sight velocity separation of 500 km s⁻¹. The location of the centre of the structures, their systemic redshift, and their dispersion $\sigma_v$ are determined from the distribution of their galaxy members, and their radius is determined from their dispersion (see the survey paper Epinat et al., in prep. for more details). In this analysis we consider four different environmental tracers. The first is the simplest of all and was already used in Mercier et al. (2022): it corresponds to the richness of the structures. However, it is a quite crude tracer of the environment because it does not take into account the distribution of galaxies in the structures. The second tracer is called a global density estimator and is defined as (e.g. Noble et al. 2013; Pelliccia et al. 2019)

$$ \eta = \frac{R_{\text{proj}}}{R_{200}} \times \frac{|\Delta v|}{\sigma_v}, \qquad (1) $$

where $R_{\text{proj}}$ is the projected distance of a galaxy with respect to the structure’s centre, $R_{200}$ is the radius where the density of the structure is equal to 200 times the critical density of the Universe, and $\Delta v$ is the systemic velocity of the galaxy with respect to the structure. This tracer is global in the sense that it does not take into account the distribution of all the galaxies in the structures but only that of the galaxy of interest. However, it has the advantage to combine in a single parameter the distance of the galaxy and its velocity. Thus, it is interesting to probe the dynamics of the galaxy and how bound it is to its host structure.

The two last tracers are different from the others in the sense that they directly probe the density in the structures at the location of the galaxies. Both density estimates rely on the concept of Voronoi tessellations where a single cell is associated to each galaxy. A galaxy located in a low-density environment will have few neighbours and will therefore have a cell spanning a large area whereas the opposite is true for a galaxy located in a high-density environment. The density can then be calculated for each galaxy by taking the inverse of its associated cell’s surface. The reason why we have two density estimates using Voronoi tessellations is because they do not take into account the same data, which can lead to slight variations in the way it is derived. The first density estimate uses the VMC approach by combining photometric redshifts from the Cosmos2015 catalogue of Laigle et al. (2016) and spectroscopic redshifts from zCOSMOS (Lilly et al. 2007) and VUDS (Le Fèvre et al. 2015). Using this method to derive the density of the structures has the advantage that it combines into a probabilistic model both photometric and spectroscopic redshifts that are sampled quite homogeneously throughout the COSMOS field. The method works by producing 100 Monte-Carlo realisations of the redshift distribution where the redshift of each galaxy is randomly assigned given their photometric or spectroscopic redshift and the associated uncertainty. For each realisation, 750 km s⁻¹ wide redshift slices are produced per step of 150 km s⁻¹ in the redshift range 0.2 < z < 1.5. For each slice, Voronoi tessellations are produced and from them we reconstruct a density grid of 75 kpc wide pixels across the field. The final density map is then constructed as the mean of the 100 Monte-Carlo realisations. Finally, because the sampling might differ from one redshift slice to another (e.g.

<footnote n="2">² https://gitlab.lam.fr/bepinat/CAMEL</footnote>

<footnote n="3">³ https://gitlab.lam.fr/bepinat/MocKinG</footnote>
<figure><img src="image_62.png" /><figcaption>Fig. 1. Criteria used for the selection of the kinematics sample. The black points represent galaxies selected according to the surface, S/N, and inclination (removing face-on galaxies only) criteria. Removed galaxies are represented as follows and in this specific order: (i) those removed by the inclination criterion (pink squares), (ii) those removed by the S/N criterion among the remaining galaxies (orange downward pointing triangles), (iii) those removed by the surface criterion among the remaining galaxies but that would have been kept by our previous size criterion used in Mercier et al. (2022, blue upward pointing triangles), and (iv) those removed by the full selection and by our previous size criterion as well (grey dots). We also show galaxies flagged with peculiar kinematics with red contours. The vertical black line shows the surface selection criterion given in Eq. 3 and the other black lines represent the S/N selection for different disk axis ratios (see Eq.4).</figcaption></figure>

due to completeness variations in spectroscopic redshifts), we rather use the overdensity parameter $\delta(z)$ defined as

$$1 + \delta(z) = \frac{\Sigma_{\text{VMC}}(z)}{\Sigma_{\text{med}}(z)}, \qquad (2)$$

where $\Sigma_{\text{VMC}}(z)$ is the mean density map from the VMC method and $\Sigma_{\text{med}}(z)$ is the median of the pixels over the map, both at redshift $z$. The only caveat with this method is that we do not use neither the MUSE spectroscopic, nor the output of the FoF algorithm. To compensate for this problem and to take advantage of the high completeness of MUSE spectroscopic redshifts at $z < 1.5$ in the MUSE fields, we compute a second density estimate using a similar approach. For each structure detected by the FoF algorithm, we perform Voronoi tessellations using MUSE spectroscopic redshifts only and from that we derive the local density at the position of each galaxy. However, this method suffers from edge effects because galaxies on the edges of the MUSE fields or of the structures have undefined boundaries. Thus, to get more precise density estimates near the edges, we also included spectroscopic redshifts from the zCOSMOS catalogue that fall within the redshift range of the structures.

## 2.4. Sample selection

In Mercier et al. (2022) we showed that, depending on the scaling relation studied, different selection criteria had to be used in order to not be biased when fitting these relations. For the analysis of the TFR, we therefore had to remove bulge-dominated, small, and low S/N galaxies. Additionally, we applied two other

selections: the first in inclination and the second in stellar mass fraction. Indeed, we had to apply the three first criteria because such galaxies usually have poor constraints on their disk model, especially their half-light radius and axis ratio. Similarly, the inclination criterion was added to remove face-on galaxies whose line-of-sight kinematics is difficult to measure because of projection effects and edge-on galaxies because their mass models are much less constrained. The last selection criterion in Mercier et al. (2022) was used to select galaxies for which we might have overestimated the contribution of the stellar components to the total rotation curve. These galaxies were visually inspected and those that had consistent morphological and kinematics models (i.e. with an intrinsically high stellar mass fraction) were kept.

For this analysis, we have decided to update these selection criteria. Because each morphological model was re-inspected as part of the derivation of the angular momentum from HST maps (see Sect. 3.2) we manually picked bulge-dominated galaxies in the kinematics sample whose disk model was loosely constrained and removed them. Such galaxies will typically have very-low surface brightness disk components (below the noise level) that can either be due to substantially large magnitudes (i.e. low fluxes), or very large disk sizes that in both cases are inconsistent with the residuals image obtained when removing the bulge model from the initial HST image. These five galaxies are 24_CGR32, 29_CGR23, 240_CGR30, 240_CGR84, and 76_CGR172. Thus, in what follows, we do not apply the bulge-to-total flux ratio selection criterion as in Mercier et al. (2022). Furthermore, we also decided to update the size selection criterion. Indeed, the size selection was used to make sure that we had enough resolution elements in the MUSE kinematics maps to properly constrain the galaxies’ kinematics. Thus, a criterion on the surface of the galaxy rather than on its size would be more appropriate, especially for edge-on galaxies. Therefore, in this analysis, galaxies are selected if their surface on the plane of the sky within an ellipse with major axis equal to one disk effective radius ($R_{eff,d}$) and apparent axis ratio ($q$) equal to that of the disk is larger than the surface of the MUSE PSF within its FWHM. Its simplified form writes

$$2 \sqrt{q} \times R_{eff,d}/\text{FWHM}_{[O\text{II}]}(z) > 1, \qquad (3)$$

where FWHM$_{[O\text{II}]}(z)$ is the MUSE PSF FWHM at the redshift of the galaxy and where $R_{eff,d}$ and FWHM$_{[O\text{II}]}$ must have the same unit. We remind that the [O\text{II}] PSF FWHM in Eq. 3 was measured in the MUSE cubes whereas the disk effective radius and apparent axis ratio were both measured in high resolution HST images. This criterion naturally excludes galaxies with small sizes with respect to the MUSE PSF and edge-on galaxies because their axis ratio is small. Thus, we also updated the inclination criterion and only kept the part that removes face-on systems, that is we ask that the inclination follows $i \ge 25^{\circ}$, as in Mercier et al. (2022). This value seems a good compromise between removing galaxies that are too face-on to be well constrained and keeping a sufficiently large sample to perform the analysis. As an indication a more conservative criterion of $i > 30^{\circ}$ would only remove eight additional galaxies. Following the same argument, we also update the S/N criterion. In Mercier et al. (2022), we only kept galaxies that have an average S/N per spaxel of at least eight across one observed effective radius (i.e. taking into account the PSF). In this analysis, we update this criterion to take into account the ellipticity of the galaxy as follows

$$\frac{(S/N)_{\text{tot}}}{\text{FWHM}_{[O\text{II}]}(z)} \ge 20 \sqrt{\pi} \left[ (x^2 + 1)(qx^2 + 1) \right]^{1/4}, \qquad (4)$$
where $x = 2R_{eff,d}/FWHM_{[OII]}(z)$, with $FWHM_{[OII]}(z)$ and $R_{eff,d}$ in arcsecond. We note that this new criterion, on top of being more physically motivated, actually adds 30 new galaxies that would have been removed by the previous S/N criterion. This number reduces to 15 when combined with the surface and inclination criteria. The new sample selection is shown in Fig. 1. Galaxies selected in surface, S/N, and inclination are shown in black and those removed by the selection with other symbols (one symbol per criterion). Because different criteria can remove the same galaxies, we show them in a specific order as follows: (i) galaxies removed by the inclination criterion with pink squares, (ii) among the remaining galaxies, those removed by the new S/N criterion with orange downward pointing triangles, (iii) among the remaining galaxies, those removed by the surface criterion that would have been kept by our previous size selection used in Mercier et al. (2022) with blue upward pointing triangles, and (iv) the remaining that are removed by our previous size selection as well (grey dots). The black vertical line shows the surface criterion defined in Eq. 3 and the other black lines show the S/N criterion defined in Eq. 4 for different disk axis ratios. The combination of the surface, S/N and inclination criteria yields a sample of 186 galaxies. As an indication, the previous selection used in Mercier et al. (2022) would have yielded 211 galaxies instead. Among them, those that would have been added by the size criterion are mostly located at the limit of the selection (thus quite small) and are significantly inclined ($i \ge 60^{\circ}$).

Finally, rather than removing galaxies with large stellar mass fraction uncertainties, as was done in Mercier et al. (2022) we decided to visually inspect the kinematics maps and the rotation curves and to flag the 30 following galaxies: (i) 38_CGR172, 54_CGR51, 70_CGR79, 74_CGR172, 90_CGR23, 93_CGR114, 101_CGR32, 104_CGR28, 104_CGR172, 105_CGR114, 113_CGR23, 148_CGR30, 185_CGR30, 226_CGR84, 257_CGR84, 267_CGR84, 313_CGR84, 442_CGR32, and 454_CGR32 because there is no visible velocity gradient in their velocity fields so that their best-fit rotation curves probably do not correctly describe the intrinsic velocity of the gas, (ii) 28_CGR26, 96_CGR28, and 172_CGR32 because there is not central peak observed in the HST image as would be expected for an exponential disk even without a bulge so that the contribution of their disk component is overestimated in the inner parts, (iii) 23_CGR84 because it has a very massive disk whose contribution is overestimated in the inner parts so that its contribution to the total rotation curve is too high and therefore produces large velocity residuals, (iv) 87_CGR35 and 345_CGR32 because their velocity field and [OII] emission are quite off-centred from their morphological center so that their velocity gradient is not correctly fitted in the inner parts, (v) 37_CGR84, 85_CGR23, and 278_CGR84 because, even though they do show some signs of rotation, their velocity fields are too perturbed for a rotating razor-thin gas disk model to be appropriate enough, (vi) 106_CGR84 because there are two kinematically distinct components in its velocity field whose centres match the locations of two morphologically distinct component in its HST image (perhaps following a merging event), therefore rendering the kinematics model uncertain, and (vii) 85_CGR35 because it is a double-peaked galaxy in its HST image with little rotation and whose velocity field residuals show that a rotating razor-thin gas disk model is not sufficient. Nevertheless, contrary to Mercier et al. (2022) where we removed such galaxies, we do not do so in this analysis. Instead, we will keep these objects and we will identify them clearly during the analysis of the angular momentum to see how they scale with the other galaxies. This way, our sample

is comparable to others found in the literature where such objects are usually not excluded and we avoid adding a complex selection function that could bias our sample.

The choice of whether a galaxy has enough rotation in its velocity field to be flagged or not is ultimately a question of perspective. Thus, we have tried to remain as conservative as possible. Therefore, in the 156 galaxies that are not flagged there might still be a few objects that are at the edge of what could be called a rotating galaxy. To illustrate the difference between a galaxy that does not have clear signs of rotation and one that does, as well as galaxies at the limit, we show examples of dynamical models for such cases in Fig. 2. For each sub-figure the morphological model is shown in the leftmost column, the velocity field model in the middle one, and the velocity dispersion in the rightmost one. We show on the top left galaxy 18_CGR114 whose morphology and kinematics are sufficiently well fitted for the analysis of the angular momentum and on the top right 70_CGR79 that does not have any velocity gradient in its velocity field and was therefore flagged. On the bottom row, we also show on the left galaxy 17_CGR34 that has a small and slightly disturbed velocity field that was not flagged and on the right 278_CGR84, a massive galaxy with a large but quite disturbed velocity field as can be seen in its residuals. Because the rotating razor-thin gas disk model does not seem appropriate to model the kinematics of this galaxy, it is unlikely that we correctly constrain its intrinsic gas kinematics and we therefore also flagged it.

# 3. Deriving the angular momentum

## 3.1. General considerations

Deriving an accurate estimate of the stellar angular momentum $J_{\star}$ in galaxies is not straightforward given that it involves knowing a priori the 3D stellar mass density $\rho_M$ as well as the 3D velocity vector $V$ of the stars. Indeed, the most general form it can take when integrated in a volume $V$ is

$$J_{\star}(\mathbf{r}) = \int_{V} d^3\mathbf{r} \rho_{M}(\mathbf{r}) \mathbf{r} \times \mathbf{V}(\mathbf{r}), \qquad (5)$$

where $\times$ represents the vector product operation. However, Eq. 5 is hardly usable in practice unless one is working with simulations (e.g. Brook et al. 2012; Stewart et al. 2013; Zavala et al. 2016; Cadiou et al. 2022) or in the vicinity of the Milky Way where 6D phase-space positions are available from Gaia (e.g. del Pino et al. 2021). Therefore, assumptions must be made on both the stellar mass density and the velocity of the stars in order to constrain their angular momentum from morphological and kinematics data. One of the most widely used expression is that of Romanowsky & Fall (2012) (hereafter RF12) which assumes that the stellar mass distribution can be described by a razor-thin exponential disk with a constant rotation curve. More recently, Posti et al. (2018) and Bouché et al. (2021) showed that the RF12 approximation is usually not precise enough to correctly compute the angular momentum of galaxies at intermediate redshift, especially in the case of low mass galaxies (typically $\log_{10} M_{\star}/M_{\odot} \le 9-9.5$) since they tend to have shallower inner slopes in their rotation curve than their high mass counterparts. A more general expression than RF12 can be derived from Eq. 5 assuming a razor-thin disk under rotation only, that is neglecting radial and vertical motions (see Appendix B for a derivation), without any assumption on the rotation curve. It writes
<figure><img src="image_63.png" /><figcaption>Fig. 2. Examples of dynamical models obtained in Mercier et al. (2022) and used in this analysis. For each sub-figure and from top to bottom we have the HST image, the GALFIT model, and the residuals in the leftmost column, the velocity field map extracted with CAMEL, the best-fit MOCING velocity field, and the residuals in the middle column, and the raw velocity dispersion map extracted with CAMEL, the beam-smearing model from MOCING, and the beam-smearing corrected velocity dispersion map from MOCING in the rightmost column. The top left sub-figure shows a galaxy with a seemingly correct dynamical model, the top right one shows a low S/N galaxy without any velocity gradient in its velocity field map that was flagged, the bottom left one shows a low S/N galaxy with a weak but non-negligible velocity gradient that was not flagged, and the bottom right one shows a massive galaxy with a kinematically disturbed velocity field that was flagged as well.</figcaption></figure>

$$J_{\star} = \int_{S} d\theta dR R^{2} \Sigma_{M}(R, \theta) V_{\theta}(R, \theta) \hat{z}, \qquad (6)$$

where S is the surface in the plane of the disk over which the angular momentum is integrated, R is the distance to the cen-

ter of the disk, $\theta$ is the azimuthal angle, $\Sigma_M$ is the surface mass density given by the integral of the 3D mass distribution along the vertical direction, $V_\theta$ is the rotation velocity, and $\hat{z}$ represents the direction perpendicular to the plane of the disk. Usually, the component of $J_{\star}$ orthogonal to the disk is normalised by the
galaxy or disk mass. It is called the specific angular momentum and it writes

$$j_{\star} = \frac{\int_{S_1} d\theta dR R^2 \Sigma_M(R, \theta) V_{\theta}(R, \theta)}{\int_{S_2} d\theta dR R \Sigma_M(R, \theta)}. \qquad (7)$$

In Eq. 7, we have written two different surfaces $S_1$ and $S_2$. A common practice is to integrate the normalisation factor on the whole disk surface (i.e. $\theta \in [0, 2\pi)$ and $R \in \mathbb{R}$) to get the total stellar mass or in a given aperture. However, we stress that the choice of the normalisation remains free and that different values of the stellar mass may be used depending on the context and the objectives of the analysis. Similarly, $S_1$ and $S_2$ do not have to be the same surface as long as it is relevant for the discussion. Besides, Eqs. 6 and 7 are technically speaking also valid for thick disks as long as they are (i) symmetric with respect to the plane of the disk and that (ii) it is possible to separate the vertical profile from the surface brightness distribution in the 3D stellar mass density (see Appendix B.1.1 for a more thorough discussion on that topic). If one further assumes a constant mass-to-light ratio throughout the galaxy and a rotation velocity that is independent of the azimuthal angle then the angular momentum integrated in a circular aperture of radius $R_1$ within the plane of the disk is given by

$$j_{\star}(R_1) = \frac{\int_{0}^{R_1} dR R^2 \Sigma(R, \theta) V_{\theta}(R)}{\int dR R \Sigma(R, \theta)}, \qquad (8)$$

where $\Sigma$ is the intrinsic surface brightness distribution (i.e. not sky projected) and where we have made implicit the limits of the integral for the normalisation factor since they are free to choose. We note that Eq. 8 is close to the expression used for instance in Bouché et al. (2021) and Mancera Piña et al. (2021a) that both rely on the assumption of a surface brightness distribution that only depends on $R$. On the other hand, our expression is slightly more general because it can account for both axisymmetric and non-axisymmetric disks (e.g. galaxies with spiral arms, bars or clumps). One potential caveat of using Eq. 8 is that it assumes a rotation velocity independent of $\theta$ which is unlikely for galaxies with non-axisymmetric features that will produce non-axisymmetric gravitational potentials. However, modelling such features is not straightforward even though technically possible (e.g. see Chemin et al. 2016) and is beyond the scope of this analysis. Besides, very few galaxies in our sample actually feature spiral arm and bar structures (less than 3%) and clumps will only locally affect the gravitational potential of their host galaxy so that their impact should be hardly visible in the velocity fields extracted from the MUSE cubes.

## 3.2. Angular momentum from HST maps

In the previous section we have defined the main equations used to derive the angular momentum. To extend beyond what was done in Bouché et al. (2021) or Mancera Piña et al. (2021a), we will use Eq. 8 in combination with HST maps to compute the angular momentum of galaxies in the MAGIC survey. There are two important quantities that must be derived beforehand: (i) the intrinsic surface brightness $\Sigma$ and (ii) the rotation velocity $V_{\theta}$. For the former we use the HST F814W images that were used in Mercier et al. (2022) to derive the galaxies’ morphological parameters and for the latter the best-fit rotation curves also

derived in Mercier et al. (2022) from the mass modelling performed on the kinematics maps extracted from the MUSE cubes (see Sect. 2.2 for a summary). To do so, we approximate Eq. 8 by discretising it along the spatial dimensions $x'$ and $y'$ of the HST image. Thus, each pixel at position $(x', y')$ has a specific angular momentum that writes

$$j_{\text{pix}}(x', y') = \frac{R F(x', y') V_{\theta}(R)}{\sum_{x', y'} F(x', y')}, \qquad (9)$$

where $R = \sqrt{x^2 + y^2}$ is the distance of the pixel to the centre of the galaxy in the plane of the disk, with $(x, y)$ the coordinates in the plane of the disk, and $F(x', y')$ is the flux in the pixel. Similarly to Eq. 8 we have let the bounds of the denominator implicit. The conversion from the apparent position $(x', y')$ of a pixel on the plane of the sky to its location $(x, y)$ in the plane of the disk is done by taking into account the position of the centre, the position angle (PA), and the inclination of the disk, all derived in Mercier et al. (2022) from the bulge-disk decomposition performed with GALFIT, assuming a razor-thin disk geometry (i.e. elliptical isophotes). The specific angular momentum measured in a circular aperture of radius $R_1$ in the plane of the disk is therefore just the sum of the specific angular momentum of each pixel that falls within the aperture:

$$j_{\star}(R_1) = \sum_{\{x', y' | R < R_1\}} j_{\text{pix}}(x', y'). \qquad (10)$$

Contrary to the surface brightness, we cannot use in Eq. 8, 9 and 10 the velocity field extracted from the MUSE cubes because (i) the MUSE observations are much less spatially resolved than the HST data (0.2" per spaxel for a PSF FWHM of 0.5" on average for MUSE versus 0.03" per pixel for a PSF FWHM of roughly 0.1" for HST), (ii) the velocity fields extracted from the cubes are too severely affected by beam smearing, especially in the inner parts where the velocity gradient and the ionised gas flux are large, and (iii) the velocity fields are projected onto the sky and there is no trivial way to invert the projection, in particular along the minor axis. Thus, we use instead the rotation curves obtained from the forward mass models performed on the velocity fields in Mercier et al. (2022). These have the advantage of being unaffected by beam-smearing or projection effects and they can be interpolated to any radius $R$. However, the main drawback is that these rotation curves are that of the ionised gas whereas $V_{\theta}$ in Eqs. 8, 9, and 10 should be that of the stars. We provide in Sect. 3.3 a discussion regarding the reliability of such an estimate of the galaxies’ stellar angular momentum.

## 3.3. Application of the HST formalism to MAGIC

Our goal for this analysis is to estimate the stellar angular momentum for intermediate redshift galaxies in the MAGIC survey. Given Eqs. 5, 8, 9, and 10, this means that we should measure the mass distribution and the kinematics of the stellar component. However, our sample is comprised of a majority of star-forming galaxies with strong emission lines (e.g. [O II]) so that we would be able to derive the stellar kinematics only for a small fraction of the entire sample. However, we do have a proxy for the stellar mass distribution through the galaxies’ HST F814W images. Hence, our method is to estimate the stellar angular momentum using Eqs. 9 and 10 with the stellar flux distribution given by the galaxies’ HST images but using the rotation curves derived from the [O II] doublet in Mercier et al. (2022). The underlying
assumption of this method is that there is co-rotation between the gas and the stars in the galaxies. We argue that this estimate should be close to the true stellar angular momentum for a couple of reasons. First, if the underlying assumptions behind the kinematics modelling performed in Mercier et al. (2022) are appropriate (e.g. rotating razor-thin disk) then the dynamics of the stars and the gas should be driven by the gravitational potential of the galaxy, which is the same for both components. Hence, their bulk dynamics (i.e. not taking into account local variations, for instance because of out-of-equilibrium motions) should match. The only effect that can differ between the gas and the stars is the impact of asymmetric drift. However, the asymmetric drift-corrected circular velocity is supposed to trace the total mass distribution and therefore, as long as we take that effect into account, using the ionised gas kinematics as a proxy for that of the stars should be appropriate enough. As an indication, a similar argument can be found in Guérou et al. (2017) where it is shown that the asymmetric drift-corrected stellar and ionised gas kinematics are consistent with each other for galaxies in the same redshift interval as ours.

We note that similar methods have already been used in a few previous studies (e.g. Cortese et al. 2016; Di Teodoro et al. 2021). However, they usually take into account the contribution of all the components (disk, bulge, spiral arms, bar, etc.) to the observed surface brightness distribution. Doing so may not be entirely appropriate because (i) the aforementioned equations are technically valid if the stellar mass is located in a disk, which is not the case for the bulge component, and (ii) bulges may be dispersion dominated systems in which case they angular momentum should be negligible. Thus, to not be biased by the bulge component that can contribute significantly to the flux in the inner parts, we remove from the HST image its best-fit GALFIT bulge model. In the few cases where the flux becomes negative near the centre in the new bulge-removed HST image, for instance because its contribution was slightly overestimated by GALFIT, we replace the pixels with negative values by the flux of the disk model at the same location. In what follows, we will nevertheless also investigate the impact of adding the bulge to the flux distribution to the Fall relation.

Finally, we must provide a last discussion regarding the impact of noise on our estimate. In Eqs. 9 and 10 we use the flux of the pixels to approximate Eq. 8 when using the HST images to derive the stellar angular momentum. However, HST images are affected by noise which can impact the angular momentum estimate. Because the background signal has already been removed from the images, either beforehand or by using an additional sky background component during the morphological modelling performed with GALFIT, we expect the noise to have an almost null mean value. Thus, if the integral in Eq. 8 or the sum in Eq. 10 are taken along a circle at fixed $R$ and if there are a sufficiently large number of pixels along the path of integration/summation, then we expect the contribution of the noise to the angular momentum to be close to null as well. This argument does not hold near the centre where we do not have enough pixels to properly sample the noise distribution and because systematics may arise from the removal of bulges. However, because the rotation velocity is negligible in the inner parts, the angular momentum should be small at $R \approx 0$ and the noise contribution should still be negligible.

## 4. Analysis

For the analysis, we use the 186 galaxies from the kinematics sample, 30 of which have been flagged in Sect. 2.4 as having peculiar kinematics.

<figure><img src="image_64.png" /><figcaption>Fig. 3. Fall relation for the kinematics sample in the MAGIC survey using bulge-removed HST maps. The stellar mass, angular momentum, and its normalisation were all measured within $R_{22}$, and the normalisation corresponds to the mass of the disk component only. Galaxies selected from the kinematics sample are represented with black points and other symbols represent galaxies removed by different selection criteria (see Fig. 1 for more details). Black circles with red contours correspond to galaxies flagged with peculiar kinematics.</figcaption></figure>

As discussed in the previous section, measuring the angular momentum of galaxies can mean multiple things depending on the way it is performed and the convention used. There are at least three important aspects that need to be discussed. First, where to measure the angular momentum? Should it be at infinity, thus at the cost of extrapolating the morphological and kinematics measurements? In one or two effective radii? At some kinematics radius such as $R_{22}$, the distance where the velocity induced by the disk component is the highest? Aspect number two is: which normalisation should we use? In the same radius where the angular momentum is measured or rather at infinity, once again at the cost of extrapolation? And finally, the last aspect is: how should we compute it? Should we use HST maps or analytical approximations? If the former, should we take into account the flux of the bulge or not? All these questions are important because they affect the value of the angular momentum and thus we address them in what follows. By default, we use the formalism based on HST maps discussed in Sect. 3.2, we measure the angular momentum within $R_{22}$, and similarly for the normalisation$^4$. We also remove the contribution of the bulge to the HST map by default before deriving the angular momentum. Later on, we investigate whether changing the radius, the normalisation, or taking into account the bulge have an impact on our results and if so how.

To study the shape of the Fall relation, we use the same fitting methodology as in Mercier et al. (2022). We always fit the relation using the expression

$$ \log_{10} j_{\star} [\text{kpc km s}^{-1}] = \beta + \alpha(\log_{10} M_{\star} [\text{M}_{\odot}] - p), \quad (11) $$

<footnote>$^4$ We note that photometric and kinematics measurements beyond this radius become less robust because of the significant drop in S/N.</footnote>
where $\beta$ acts as the zero point of the relation and $\alpha$ as the slope, $p$ is the pivot point taken to be the decimal logarithm of the median stellar mass in the kinematics sample, and $\log_{10} j_*$ is the dependent variable in the fit. Thus, we always use $p = 9.8$. The point of the pivot is to reduce the correlation between the slope and the zero point during the fit. Indeed, the further the data points are from the origin, the larger the zero point is impacted when there is a small variation in the slope. With a pivot, the data points are closer to the origin and the effect of the pivot is thus reduced. In what follows, we mainly consider the change in zero point between different Fall relations when fixing the slope unless there is a significant variation between subsamples, in which case we take it into account. To perform the fits we use a slightly updated version of LtsFrr<fnref n="5" /> (Cappellari et al. 2013) that implements the Least Trimmed Squares (LTS) technique described in Rousseeuw & Driessen (2006) to find and remove outliers from the fit. The only feature that our updated version adds is the possibility to fit while fixing the slope at the same time. When fixing the slope, we use the following procedure: first, we fit the entire sample under consideration with a free slope and second, we fit each subsample separately using this slope. Another advantage of LtsFrr is that it can take into account uncertainties on both the dependent (angular momentum) and independent (stellar mass) variables during the fit. For the stellar mass, we use the uncertainty from the SED fitting and, according to Mercier et al. (2022), we add a systematic uncertainty of 0.2 dex. For the angular momentum, we did not find in the literature (e.g. Bouché et al. 2021; Mancera Piña et al. 2021a) a typical systematic uncertainty that the authors add to perform the fit, though it is mentioned in Romanowsky & Fall (2012) that they use bootstrap resampling to estimate their uncertainties. After some tests, we found that adding no systematic uncertainty on the angular momentum can bias the fit because we treat asymmetrically the dependent and independent variables. Thus, we decided to add a systematic uncertainty of 0.1 dex on the angular momentum when fitting. This value is quite close to the typical uncertainty on the angular momentum evaluated by computing the difference between the HST-derived angular momentum estimates when using a rotation curve derived from the mass models and when using a simpler flat model rotation curve (see Sect. 4.1 for more details). Additionally, we also add quadratically for each galaxy an estimate of the uncertainty on $j_*$ due to the propagation of uncertainties taken as the standard deviation of 100 Monte-Carlo realisations of their HST map and kinematics parameters assuming Gaussian perturbations. This method nevertheless tends to under-predict the uncertainties on $j_*$ which therefore justifies the addition of a systematic uncertainty.

The Fall relation is shown in Fig. 3 for the kinematics sample. Similar symbols to Fig. 1 are used, with selected galaxies represented with black circles and flagged galaxies with red contours. The stellar mass (bulge plus disk), angular momentum, and its normalisation are all measured within $R_{22} = 2.2R_d$, with $R_d$ the disk scale length, and we normalised the latter by the disk mass only. We recover the typical Fall relation with the angular momentum that seems to scale linearly with stellar mass. Our range of angular momentum values and the scatter in the relation are also consistent with recent results obtained with MUSE (e.g. Bouché et al. 2021). The main difference between our two studies is regarding the component that is traced by the angular momentum estimate since, in Bouché et al. (2021), it is that of the ionised gas that is traced, whereas in this analysis we try to estimate the stellar angular momentum. Still, we find a similar

Fall relation. Among the galaxies that would have been selected by the size criterion, the edge-on objects that would have been removed by the inclination selection (blue upward pointing triangles) tend to be spread out further than the selected galaxies above and below the relation. Those that would have been retained by the size and inclination selection criteria (orange downward pointing triangles) tend to be located throughout the Fall relation with more or less the same scatter, but populating more the low-mass end. Other galaxies that would have been removed anyhow by the size selection criterion (grey dots) mainly populate the low-mass and low-angular momentum part of the relation and even seem to be mostly found below the average relation. The fact that these galaxies are mostly found at low-mass is expected since, as was shown in Mercier et al. (2022), they fall on the size-mass relation. Finally, we note that flagged galaxies mostly populate the bottom part of the Fall relation which makes sense since the majority do not exhibit any visible velocity gradient. Thus, their rotation velocity is loosely constrained so that a DM halo with a nearly null maximum velocity is fitted to minimise the residuals. In other terms, the contribution of the stellar disk and bulge to the total rotation curve is likely overestimated for these galaxies.

## 4.1. Assessing the reliability of the method

First, we must ensure that we correctly measure the angular momentum when using HST maps. There are two main aspects that can change the value of the angular momentum with respect to more traditional estimates: (i) the fact that we derive the ionised gas kinematics and thus the total rotation curve using mass models. In Mercier et al. (2022), we showed that this could produce a limiting line on the TFR acting as a lower bound that had to be taken into account. And (ii) the fact that we use HST maps instead of assuming an axisymmetric surface brightness profile. Thus, in what follows we discuss how these two aspects impact the shape of the Fall relation. In our case, there are four different ways we can compute the angular momentum: (i) using the HST maps and the rotation curve derived from the mass modelling, (ii) using the HST maps and a simpler flat model rotation curve described in Wright et al. (2007); Epinat et al. (2009, 2010b); Abril-Melgarejo et al. (2021), see also Eq. B.9, which is not physically motivated but that can nevertheless robustly recover the rotation velocity for slowly declining or constant rotation curves at large radii, (iii) numerically integrating Eq. 8 using the best-fit axisymmetric exponential disk model and the rotation curve from the mass models, and (iv) using the analytical expression for the angular momentum of an exponential disk with a flat rotation curve that is given by Eqs. B.10 and B.13. The Fall relations for these four estimates are shown in Fig. 4. All are derived within $R_{22}$ using the bulge-removed HST maps and similarly for the normalisation of $j_*$. The same symbols have been used as for Fig. 1. This figure can be read along two directions. First, horizontally where the first row shows the Fall relations obtained with HST maps when using the rotation curve derived from mass models (top left) and when using a flat model for the rotation curve (top right). The second row shows instead the Fall relations when using an axisymmetric exponential disk surface brightness model with a rotation curve from mass models (bottom left) and with a flat rotation curve model (bottom right). Plots located on the same row therefore allow to see the effect of using different rotation curve models. The figure can also be read vertically, in which case the rotation curve model remains the same between the two plots and only the surface brightness map changes (HST on the top, exponential disk on the bottom).

<footnote n="5">⁵ https://pypi.org/project/ltsfit/</footnote>
<figure><img src="image_65.png" /><figcaption>**Fig. 4.** Fall relation for the kinematics sample using various formalisms. On the top row is shown the angular momentum derived using HST maps and on the bottom row that derived using an exponential disk model. The leftmost column represents the Fall relation when using the rotation curve derived from the mass modelling and the rightmost column the Fall relation when using a flat model for the rotation curve. Stellar masses and angular momenta are all measured within $R_{22}$, as is the normalisation used for $j_\star$. Galaxies selected for the analysis are shown as black points and those flagged with peculiar kinematics with red contours. Other symbols correspond to galaxies removed by the selection (see Fig. 1 for more details). The two black lines in each plot show the limit of the Fall relation for an exponential disk without DM halo and with two different disk scale lengths of 5 kpc (dashed line) and 15 kpc (plain line, see Eq. B.18).</figcaption></figure>

Hence, this allows to see the effect of using non-axisymmetric HST maps compared to a more traditional axisymmetric exponential disk model

To begin with, we consider the effect of using mass models for the rotation curve. By comparing horizontally the plots in a same row in Fig. 4, we can see that there are large discrepancies found below roughly $10^9 M_\odot$ which correspond for the vast majority to galaxies removed by the sample selection, and in particular to galaxies with small sizes (grey dots). Above this stellar mass threshold, the two relations are more alike though slightly different. The selected galaxies (black points) are consistently spread along a line with a somewhat similar slope and zero point, though we can clearly see that the Fall relation obtained with mass models seems to have a lower limiting line similar to what could be seen for the TFR in Mercier et al. (2022). This effect

is striking for the galaxies flagged with peculiar kinematics (red contours) that are found for the majority along this limiting line when using mass models but that are spread across it when using the flat model. This result actually makes sense because most of the flagged galaxies are so because they lack a clear velocity gradient in their velocity field (e.g. 70-CGR79 in Fig. 2). This means that the contribution of the stellar components to the rotation curve of these galaxies is certainly overestimated in which case the kinematics model can only add a small fraction of DM. In turn, this produces a correlation between the galaxies’ stellar mass and their stellar angular momentum. The only aspects that render this correlation not perfect are the non-axisymmetric features in the HST maps such as spiral arms or clumps and the relative contribution of the bulge component with respect to the disk.
Then, we investigate whether using HST maps instead of an exponential disk model has an important impact on the estimate of the angular momentum. This is shown by comparing vertically the plots in a same column in Fig. 4. Contrary to the previous comparison, we see that the HST maps have only a small impact on the Fall relation, with an average absolute difference in angular momentum below 0.1 dex which is about the same order of magnitude as the average uncertainty in the kinematics sample. Larger differences can be seen at the lowest stellar masses but they mostly concern small galaxies (grey dots) that certainly have loosely constrained disks. This result holds independently of whether we consider a rotation curve from mass models or from a flat model, which strengthens the conclusion that the choice of rotation curve is the leading factor that will affect the value of the derived angular momentum. We can also see the same kind of limiting lower line as for the HST maps when using mass models. We have previously mentioned that this line is due to galaxies whose baryon contribution was overestimated, which in turn produces a correlation between their stellar mass and their stellar angular momentum. For HST maps it is not trivial to derive an expression for this limiting line because it will depend on the flux distribution in the image. However, in the case where the angular momentum is derived through numerical integration, thus assuming an axisymmetric exponential disk, it is possible to derive an expression for this limiting line. More details are given in Appendix B.5, but basically the idea is to derive the expression of the angular momentum for a disk galaxy that has no DM halo. We do not take into account the impact of the bulge because the expression would already become much more complex. Besides, adding a bulge component will increase the flux and the rotation velocity which will in turn increase the angular momentum. For an exponential disk alone, the velocity of the ionised gas must compensate the disk's gravity and so there is a link between its mass and its angular momentum. In the specific case where we measure the angular momentum within $R_{22}$, there is an analytical expression for the limiting line (see Eq. B.18 that depends on the disk's mass and on its radius. The latter only affects the zero point of the line (smaller disks have less angular momentum at fixed stellar mass) but not its slope which is always equal to one-half. Two examples of such limiting lines are shown on the bottom right plot of Fig. 4 for galaxies with disk a scale length of 5 kpc or 15 kpc.

Hence, these two comparisons suggest that the choice of the rotation curve is what has the largest impact on the estimate of the angular momentum. In particular, we must be careful when using the mass models to not be biased by baryon-dominated galaxies that lie on the bottom of the Fall relation along the kind of limiting line that has been previously discussed. The same selection and kinematics modelling were applied to the entire sample independently of the characteristics of the galaxies' environment so we expect the effect of this lower limiting line to be independent of environment. Therefore it should not bias this analysis unless perhaps if the fraction of baryon-dominated galaxies changes with environment. Besides, as shown in Fig. 4, selected galaxies from the kinematics sample are the least affected compared to very low-mass galaxies as they seem to populate similar Fall relations when using mass models or a flat one for the rotation curve.

## 4.2. Impact of environment

As discussed in Sect. 2.3, we have different proxies for the environment that we can use. Each proxy does not probe the environment in the same way and therefore it is interesting to consider all of them for this analysis.

### 4.2.1. Environment traced by structure richness

<figure><img src="image_66.png" /></figure>

**Fig. 5.** Fall relation between galaxies in the field (black points) and those found in the richest structures (more than ten members - red circles). The stellar mass, angular momentum, and its normalisation are all measured within $R_{22}$ and the normalisation corresponds to the mass of the disk component only. Galaxies beyond the mass cut at $M_{\star} = 10^{10} M_{\odot}$ are shown with semi-transparent symbols and detected outliers are shown with yellow stars. Two best-lines are shown for each sub-sample (black line for the field and red for the structures). The dashed lines represent the best-fit Fall relations using the entire mass range whereas the plain lines correspond to the aforementioned stellar mass cut and an additional redshift cut $0.5 < z < 0.9$.

To begin with, we consider the simplest environment estimator that was used in Mercier et al. (2022) and that defines the environment in terms of the structures' richness (i.e. number of galaxy members). We consider the two following sub-samples: (i) field galaxies, that is galaxies that were not associated to any structure by the FoF algorithm or that were associated to a structure with less than three members, and (ii) galaxies located in large structures that have at least ten galaxy members. As in Mercier et al. (2022), we do not consider smaller structures with less than ten members because we lack statistics for this sub-sample to properly constrain it. The Fall relation for these two sub-samples is illustrated in Fig. 5 with field galaxies shown with black points and those in structures with red circles. As discussed in our previous analysis of the TFR, there is an excess of massive galaxies found in large structures with respect to the field. This could be problematic because if the slope of the Fall relation is different for the high-mass end independently of environment, then it will mostly impact galaxies in the structures and we may wrongly interpret this result as an effect of environment. Hence, we fit the Fall relations with (plain lines) and without (dashed lines) a mass cut at $M_{\star} < 10^{10} M_{\odot}$, which is the stellar mass below which the field galaxies and those in large structures have the same stellar mass distribution, to ensure we are not biased by the most massive galaxies. As in Mercier et al. (2022), we also applied a redshift cut $0.5 < z < 0.9$ to ensure
<figure><img src="image_67.png" /><figcaption>**Fig. 6.** Dispersion (data minus best-fit line, neither mass nor redshift cuts) along $j_*$ as a function of stellar mass. Black points represent field galaxies, red circles those in structures, and yellow stars outliers. Galaxies flagged with peculiar kinematics are also shown with cyan crosses. The median, 16th, and 84th percentiles around the best-fit lines are shown as crosses and error bars for various 1 dex wide and partially overlapping stellar mass bins separated by 0.5 dex. We note that we used the best-fit line of each sub-sample separately (field or structures) to derive the galaxies’ dispersion. There is no significant difference in scatter between both subsamples, except in the last mass bin. However, we note that there are too few galaxies in the field around $10^{11} M_{\odot}$ to properly constrain this last bin. This result also holds when applying mass and redshift cuts.</figcaption></figure>

we are not probing a potential redshift evolution of the relation given that galaxies in structures are mainly located at $z \approx 0.7$ whereas field galaxies are spread throughout the entire redshift range $0.2 \le z \le 1.5$. In Fig. 5, we directly show the best-fit linear relation (plain line) that takes into account both the stellar mass and the redshift cuts. However, we note that, even though the zero-points are higher when we only apply the mass cut, the offset between the two subsamples remains the same (0.08 dex), modulo a slightly steeper slope. As an indication, we also show the galaxies removed by the selection in the background with grey symbols and outliers detected by LrsFrr as yellow stars.

Interestingly, this first result may suggest that there is a visible impact of the environment on the Fall relation where field galaxies have a higher angular momentum (0.12 dex without any cut, 0.08 dex with mass and redshift cuts applied) than those in large structures at fixed stellar mass. We note that this result holds even when we apply the stellar mass and redshift cuts even though the slopes differ significantly. The zero-points are actually quite well constrained with an uncertainty of at most 0.02 dex, though their offset is significantly below the typical scatter in angular momentum of 0.3 dex measured in stellar mass bins of 1 dex. We also note that the result holds even when using the flat model rotation curve (see Fig. A.1). As discussed in Sect. 4.1, the advantage of the flat model is that it is not constrained by the mass modelling and thus it can probe galaxies that might be located below the Fall relation. Without mass an redshift cuts, we actually find a much larger difference of 0.22 dex (instead of 0.12 dex for mass models) between subsamples. However, when applying both cuts we recover the same 0.08 dex zero-points difference, modulo a slightly steeper slope. There may be two interpretations for this difference in zero point. Either there is a global shift of the Fall relation in the largest

structures with respect to the field or there is a larger scatter in a specific direction (seemingly towards lower angular momentum at fixed stellar mass for galaxies in structures). To investigate these interpretations, we computed the dispersion along $j_*$ around the best-fit line for each sub-sample separately, as is shown in Fig. 6. On the left-hand side, we represent the dispersion around the best-fit line (without any cut) when using the rotation curve derived from the mass models, similarly to Fig. 5, and on the right-hand side the dispersion (without cuts as well) when using the flat model rotation curve (as in Fig. A.1). To help visualise the scatter around the Fall relation, we represent the median with crosses and the 16th and 84th percentiles with error bars for each sub-sample, computed in 1 dex wide stellar mass bins. Each bin is separated by 0.5 dex, starting at $10^9 M_{\odot}$ and ending at $10^{10.5} M_{\odot}$. We also performed the same exercise using the best-fit lines obtained when applying the mass and redshift cuts and we arrive at a similar conclusion. The only notable difference, which is visible in Fig. 5, is that nearly every massive galaxy is located below this line, but there is not visible trend with environment. We note that because we used four different best-fit Fall relations to derive the scatter (two for the two subsamples with the mass models and two for the two subsamples with the flat model), Fig. 6 only allows us to directly appreciate the dispersion around the relations and not offsets between different subsamples though, when combined with Figs. 5 and A.1, it is possible to draw some conclusions on that point. Looking at Fig. 6, we can notice that there is no significant difference between field galaxies and those in large structures in terms of scatter. For most of the stellar mass bins, the median value represented by the cross symbol is close to zero which shows that the slope is well constrained. Only for the highest mass bin can we see a significant difference. However, we note that in this mass
range there are very few galaxies found in the field and therefore we cannot draw clear conclusions for this bin. The same applies when using the flat model rotation curve and when applying a mass cut or mass and redshift cuts together. Similarly, removing galaxies flagged with peculiar kinematics before fitting the Fall relation does not change our conclusions in terms of zero-point difference and scatter around the best-fit relation. However, we note that there is a significant change in slope when doing so in the sense that the new value is much closer to the slope obtained when fitting the relation with a stellar mass cut. The effect is especially important for the Fall relation when using the flat model rotation curve. This indicates that flagged galaxies, because they mostly populate the bottom part of the Fall relation (mass models) or are located below it (flat models), tend to drive the slope towards a lower value.

Thus, the assumption that the difference in zero-points observed between subsamples is due to a larger scatter for the galaxies in the large structures does not seem to hold. Hence, this would point towards the other explanation that there is a global negative offset (in angular momentum at fixed stellar mass) of the relation for galaxies located in the largest structures with respect to the field. We note that there is a significant change in slope when we remove the most massive galaxies beyond $10^{10} M_{\odot}$ that make the best-fit Fall relation under-predict the angular momentum of low-mass galaxies below $10^9 M_{\odot}$ (see Fig. 5). The reason is that the best-fit line is actually driven by intermediate mass galaxies in the range $10^9 M_{\odot} < M_{\star} < 10^{10} M_{\odot}$. Hence, this global offset mostly applies for galaxies beyond $10^{10} M_{\odot}$.

### 4.2.2. Environment traced by other estimators

Using the structures' richness to probe the effect of the environment is an interesting first step but, for galaxies located in structures, it lacks some details. Indeed, depending on the type of structure (e.g. compact versus diffuse group) and the location of the galaxies in this structure, the environment can affect them differently. Thus, it is interesting to weigh our previous findings by looking at the Fall relation as a function of other environmental tracers.

First, we consider the global environment estimator $\eta$ defined in Sect. 2.3. This estimator does not take into account the distribution of galaxies in the structures, but it allows to separate galaxies depending on their location and on their dynamics with respect to the structure they belong to. Low values of $\eta$ means a galaxy is either located in the inner parts of the structure or that it has a low velocity relative to the structure's velocity dispersion (or both). Hence, by classifying galaxies using this parameter we can separate those that have certainly been largely affected by their environment from those that have been barely affected. In this analysis, we follow a similar classification to what was done in Pelliccia et al. (2019). We consider the three following classes:

(i) $0.1 < \eta$ corresponding to "backsplash" galaxies, that is galaxies that have already passed the pericentre of their orbit once (e.g. Balogh et al. 2000; Gill et al. 2005) and that are thus expected to be highly affected by their environment,

(ii) $0.1 < \eta < 0.4$ corresponding to recently accreted galaxies that may have started to be affected, and

(iii) $\eta > 0.4$ corresponding to galaxies located in the outer parts and/or loosely bound to the structure (e.g. because of a large systemic velocity).

<figure><img src="image_68.png" /><figcaption>**Fig. 7.** Fall relation colour coded according to the global environment tracer $\eta$ defined in Sect. 2.3. The stellar mass, angular momentum, and its normalisation are all measured within $R_{22}$, and the normalisation corresponds to the mass of the disk component only. We have used a similar classification as Pelliccia et al. (2019) with "backsplash" galaxies in red, recently accreted galaxies in blue, and galaxies in the outer regions of the structures in orange. Field galaxies are represented with black points. As an indication we also show galaxies removed by the selection in grey in the background. The best-fit line for each class of galaxies are shown with dashed lines when there is neither mass, nor redshift cut and with plain lines when there are both cuts applied. The stellar mass and angular momentum histograms are also shown on the sides of the figure.</figcaption></figure>

The Fall relation for these three classes is shown in Fig. 7 with "backsplash" galaxies in red, those recently accreted in blue, and those loosely bound in orange. Field galaxies are also shown as black points and the stellar mass and angular momentum histograms for each class are represented on the sides of the main plot. The first striking feature is the difference in the stellar mass distributions between the field and the three aforementioned classes. We have already discussed in previous sections that galaxies located in structures populate the high-mass end of the Fall relation whereas only a handful of field galaxies are found beyond $10^{10} M_{\odot}$. In this figure, we actually see a clear transition from galaxies located in the outer parts of the structures, to recently accreted galaxies that populate mostly the $10^{9.5} M_{\odot} < M_{\star} < 10^{10.5} M_{\odot}$ range, and then to "backsplash" galaxies that are found at the high-mass end. This transition is also visible in the angular momentum histograms, though perhaps not as clearly. We note that there are only a few galaxies with $\eta > 0.4$ which makes it difficult to constrain their Fall relation. We have decided to keep them as a separate class rather than integrating them into the "recently accreted" class to avoid being potentially biased by the fact that these two populations of galaxies may be impacted by the environment in different ways.

To probe the impact of the environment, we fitted each class of galaxies separately, using the same slope found in Sect. 4.2.1. The best-fit lines are shown with (plain lines) and without (dashed lines) mass and redshift cuts applied, as in Fig. 5. We note that the best-fit lines with the mass and redshift cuts applied are difficult to interpret for various reasons. For "backsplash" galaxies, this is because most are located beyond $10^{10} M_{\odot}$. For recently accreted galaxies, there is a signif-
<figure><img src="image_69.png" /><figcaption>**Fig. 8.** Fall relation using the mass models rotation curve colour coded according to two local density estimates. On the left-hand side is shown the projected density (galaxies per unit area) when using MUSE and zCOSMOS spectroscopic redshifts only and on the right-hand side the overdensity using photometric and spectroscopic redshifts described in Sect. 2.3. Galaxies are split between low- (orange), medium- (blue), and high-density (red) environments. For each plot, the best-fit line for each class of galaxies is shown with a dashed line (plain line when mass and redshift cuts are applied). The stellar mass, angular momentum, and its normalisation are measured within $R_{22}$. The normalisation is taken as the mass of the disk component only.</figcaption></figure>

icant fraction below $10^{10}$ M$_\odot$, but the vast majority are clus-
tered above $10^{9.5}$ M$_\odot$ which makes it particularly difficult to fit
a straight line without having large uncertainties on the zero-
point. And for the last class, it is because, as discussed before,
there is already few galaxies without any cut, thus even less re-
main after the cuts are applied. When no cut is applied, we find
a visible impact of the environment on the zero-point of the rela-
tion with field galaxies that lie on the top, followed by galaxies
in the outer parts of the structures (difference in zero-point of
$\Delta\beta = -0.04$ dex with respect to the field), then by recently ac-
creted galaxies ($\Delta\beta = -0.09$ dex), and finally by "backsplash"
galaxies ($\Delta\beta = -0.19$ dex). However, as soon as we take into ac-
count the mass cut or both the mass and redshift cuts, these dif-
ferences disappear. Only for the second class (recently accreted
galaxies) do we find a large offset from the best-fit lines of the
field galaxies and of those found in the two other classes, but
we remind that, as discussed before, this class is by far the least
constrained of the three, and so this result is not particularly sig-
nificant. The same result applies when using the flat model ro-
tation curve with a larger zero point difference without any cut
($\Delta\beta = -0.25$ dex) and a similar one with mass and redshift cuts
($\Delta\beta = -0.08$ dex), as illustrated in Fig. A.2. Combined with the
previously mentioned fact that "backsplash" galaxies are clearly
separated in stellar mass from the other galaxies, but not as much
in angular momentum (see the two histograms in Fig. 7), this
may be taken as an indication that massive "backsplash" galax-
ies deviate from the Fall relation found at lower stellar masses.
We note that this interpretation is consistent with the variations
of slope seen in Sect. 4.2.1 when applying or not the mass cut.
However, this effect may not be necessarily environmentally-

driven. Because such galaxies are the only ones to occupy the
high-mass end, we cannot eliminate the possibility that this de-
viation is solely caused by their mass, that is that if we had a
larger sample of massive field galaxies we would find the same
trend.

The parameter $\eta$ provides us information regarding the loca-
tion of the galaxies and/or how tightly bound they are to their
structure. This can give us indications regarding the impact of
the environment on the galaxies’ angular momentum because
we may expect a galaxy that has already completed an orbit
through the structure to be more affected than a galaxy that has
just started being accreted. However, structures such as galaxy
groups are not completely isotropic and, depending on their ex-
act location, galaxies may reside in more or less dense regions.
In essence, the density of a region could mean that of the intra-
cluster medium (ICM) or that of the galaxies (average number of
galaxies per unit area). In our case, we are only able to probe the
latter using one of the two local estimates discussed in Sect. 2.3.
These two estimates are close but they do not probe exactly the
same scales and they do not rely on the same data, so it is in-
teresting to treat them in parallel. The VMC estimate uses a
combination of photometric and spectroscopic redshifts (but not
from MUSE) and derived the density within redshift slices with a
fixed width whereas the other uses a combination of MUSE and
zCOSMOS spectroscopic redshifts only and is derived structure-
per-structure. Hence, the latter provides more precise density es-
estimates but the former spans uniformly the entire field. In par-
ticular, because of depth variations from one redshift slice to
another, we will consider the overdensity rather than the den-
sity for the estimate using the combination of photometric and spectroscopic redshifts (see Sect. 2.3 for more details).

We show the Fall relation for these two density estimates in Fig. 8. On the left-hand side the data are colour coded according to the estimator that uses the spectroscopic redshifts only (hence density estimate, noted $\Sigma_{\text{proj}}$) and on the right-hand side the data are colour coded according to the estimate using photometric and spectroscopic redshifts (hence overdensity, noted $\delta$). For both plots, the Fall relation is split into three classes:

(i) low-density environment where $\log_{10} \Sigma_{\text{proj}} < 0$ or $\log_{10}(1 + \delta) < 0.5$ (this also includes field galaxies),

(ii) high-density environment where $\log_{10} \Sigma_{\text{proj}} > 2$ or $\log_{10}(1 + \delta) > 1$, and

(iii) medium-density environment in-between the two previous classes.

The following results do not depend strongly on the criteria used to separate the galaxies into classes. For instance, only a handful of galaxies are added to the low-density environment when using $\log_{10} \Sigma_{\text{proj}} < 1$ instead of the aforementioned value. Nevertheless, we keep the criterion $\log_{10} \Sigma_{\text{proj}} < 0$ because it seems physically appropriate to separate low- from medium-density environments. We note that because of the size of the MUSE fields, this population of galaxies will be dominated by field galaxies. It is interesting to see that we recover similar results to what was previously seen in other figures. In particular, we find that massive galaxies are mostly located in high-density (or overdensity) environments, though the trend is not as clearly visible as in Fig. 7. By comparing Figs. 7 and 8, we can note that some galaxies identified as "backsplash" are found in low-density environments. This may be expected given that such galaxies can be classified as "backsplash" even if they are located far from the centre where the density is lower (e.g. because they have passed the pericentre and have now reached the apocentre of their orbit). When fitting each class separately, we find similar results to those obtained with the previously discussed environment tracers. Without any mass cut, we recover the transition from the field/low-density environment to the medium-/high-density ones. When using the density estimator (left-hand plot in Fig. 8), the zero-point of the medium- and high-density galaxy classes are switched compared to the result obtained with the parameter $\eta$. However, when using the overdensity parameter (right-hand side plot), the same zero-point is found for both classes. When applying the mass and redshift cuts, the differences between the field, the small-, and the large-density classes vanish. However, we note that for both estimators there is a significant zero-point offset ($\Delta\beta = -0.12$ dex with respect to the field) for the class probing medium-density environments. This offset is mainly driven by $M_* \sim 10^9 - 10^{9.5} M_{\odot}$ galaxies which tend to populate the bottom part of the Fall relation when using both estimators.

### 4.2.3. Summary of the impact of the environment

It seems that, independently of the environmental tracer used, galaxies located in richer and denser structures and/or that are more dynamically bound to their host structure have less angular momentum at fixed stellar mass. Nevertheless, this effect seems to be mainly driven by the most massive galaxies. Indeed, those galaxies tend to be located below the Fall relation that is found for lower mass galaxies. This effect is even more striking when using an unconstrained rotation curve (flat model) to derive the angular momentum. Hence, this effect seems real and not induced by our method.

Because we lack a statistically significant sample of field galaxies in the highest stellar mass bin ($M_* > 10^{10} M_{\odot}$), we cannot definitely reject the hypothesis that the observed difference in the zero-point of the Fall relation as a function of environment is not an effect of the galaxies’ mass rather than environment. Still, we can note that there is a strong correlation with environment: these massive galaxies that show a deficit of angular momentum tend to be on average more dynamically bound to their host structure than their lower-mass counterparts and are located in denser environments. It is also interesting to note that a similar zero point offset of 0.1 dex between a sample of field galaxies and another of galaxies in structures was also obtained in Pelliccia et al. (2019) where the two samples have comparable stellar mass distributions at the high-mass end. Their interpretation of this result is that this loss of angular momentum is mainly due to gas-poor galaxy mergers happening in the structures. However, they also clearly state that other processes might be responsible for this observed angular momentum loss (e.g. tidal interactions). Galaxy mergers might be a viable explanation in our case as well because the effect seems more predominant in high-density environments where they are more likely to happen (e.g. Lin et al. 2010). But once again, without a comparable sample of field galaxies at the high-mass end, we cannot conclude for certain this effect is environmentally driven.

### 4.3. Impact of bulges

Some authors have argued that the presence of a bulge can affect the location of the galaxies on the Fall relation (e.g. Obreschkow & Glazebrook 2014; Fall & Romanowsky 2018). As plainly discussed by Obreschkow & Glazebrook (2014), galactic bulges can have some rotation but it is nevertheless usually much smaller than that of the disk component and therefore their contribution to the unnormalised angular momentum (i.e. not divided by the mass) is negligible. However, because a significant fraction of the stellar mass can be found in the bulge, especially when evaluated relatively close to the centre (e.g. at one effective radius), they can nevertheless impact the position of the galaxies on the Fall relation and they may even explain to some extent the dispersion seen in the relation (e.g. Fall & Romanowsky 2018). Hence, there are various ways we can look at the impact of bulges: (i) following Obreschkow & Glazebrook (2014) and Fall & Romanowsky (2018), we can explore the location of the galaxies on the Fall relation as a function of the bulge fraction and the environment, (ii) we can check the impact of the bulge when taking it into account for the normalisation of the angular momentum, and (iii) we can check the effect of keeping the bulge contribution in the HST maps when estimating the angular momentum, that is assuming the bulge co-rotates with the disk and that its flux distribution can be approximated as a thick disk. Each of these three approaches can potentially give us complementary information regarding the link between angular momentum, stellar mass, and bulge fraction, ideally as a function of the galaxies’ environment.

#### 4.3.1. Impact of normalisation

To begin with, we look at the position of the galaxies on the same Fall relation as before (i.e. normalising by the disk’s mass) as a function of their bulge-to-disk flux ratio (B/D) and we check how this changes when we consider instead the total stellar mass for the normalisation. Indeed, so far we have always normalised the Fall relation using the disk’s mass but what is usually done
<figure><img src="image_70.png" /></figure>

**Fig. 9.** Fall relation colour coded according to the bulge-to-disk flux ratio (B/D) evaluated at the global effective radius for the entire kinematics sample (selected galaxies with plain symbols and those removed by the selection criteria with semi-transparent symbols). The stellar mass, angular momentum, and its normalisation are all measured within $R_{22}$, and the normalisation corresponds to the mass of the disk only (left-hand side) or to the total mass (disk plus bulge, right-hand side). Galaxies are separated between those in the field (dots) and those in structures (circles). We note that, contrary to Figs. 5 and 6, we also include galaxies located in structures with less than ten members.

<figure><img src="image_71.png" /></figure>

**Fig. 10.** Fall relation for the entire kinematics sample using the mass models rotation curve and binned in 0.5 dex wide stellar mass bins separated by 0.25 dex. The stellar mass, angular momentum, and its normalisation are all computed within $R_{22}$. Galaxies are split between those with a bulge-to-disk ratio (B/D) evaluated at one effective radius larger than one (red upward pointing triangles) and those with B/D below unity (blue downward pointing triangles). The different line styles represent various ways to derive the angular momentum. The plain and dashed correspond to the angular momentum derived using the bulge-removed HST map, whereas the dotted and dash-dotted lines represent that derived using the full HST map (i.e. including the bulge). Furthermore, the plain and dotted lines show the angular momentum derived with a normalisation taken as the disk mass only, whereas the dashed and dash-dotted lines show that derived with the total mass (i.e. disk plus bulge).

in similar studies is to rather use the total stellar mass (e.g. Bouché et al. 2021; Mancera Piña et al. 2021b). Hence, we check whether this has an impact on our results or not. In this section, we still derive the angular momentum from the bulge-removed HST maps and only in Sect. 4.3.2 do we analyse the effect of including the flux of the bulge when deriving the angular momentum.

The Fall relation using the mass models rotation curve colour coded as a function of B/D evaluated at one global effective radius is shown in Fig. 9. On the left-hand side is shown the relation when using the disk mass within $R_{22}$ as normalisation and on the right-hand side is shown the relation when using the total stellar mass within $R_{22}$ as normalisation. Furthermore, galaxies are split between those in the field (points) and those in structures (circles). The full kinematics sample is shown but, for indication, we represent galaxies removed by the selection with semi-transparent symbols. Overall, galaxies are mainly disk dominated with B/D < 1 and bulge-dominated galaxies (defined here with B/D > 1) are mainly found at the highest stellar masses. There is nevertheless a non-negligible fraction of bulge-dominated galaxies throughout the entire stellar mass range. As already observed in Mercier et al. (2022), there is no obvious trend with environment, except for the highest stellar mass bin which is dominated by bulge-dominated galaxies. Yet, from visual inspection, it seems that lower mass bulge-dominated galaxies are mostly found on the bottom part of the Fall relation compared to disk-dominated galaxies. This is particularly clear when considering the Fall relation with the total stellar mass used as normalisation (right-hand side plot) since the inclusion of the bulge mass necessarily reduces more significantly the angular momentum for bulge-dominated galaxies than for the disk-dominated ones. This effect is also clearly visible in the Fall relation that uses the flat model rotation curve (see Fig. A.3).

To highlight this observation, we compute the median value along the Fall relation in 0.5 dex wide stellar mass bins sepa-
rated by steps of 0.25 dex, as shown in Fig. 10. We separate disk-dominated galaxies represented with downward pointing triangles from bulge-dominated galaxies represented with upward pointing triangles. The lines of interest for the current discussion are the two plain lines (disk mass as normalisation) and the dashed line (total stellar mass as normalisation), which for the latter is only shown for bulge-dominated galaxies since disk-dominated galaxies nearly unaffected. For each subsample, the 16th and 84th percentiles are also shown with error bars. It is clear that bulge-dominated galaxies are located on average well below the disk-dominated galaxies subsample and that the same applies for the scatter which is preferentially found below the median value. As already inferred from Fig. 9, the effect is amplified when the normalisation uses the total stellar mass. This result is actually similar to the dependence of the Fall relation with the bulge fraction presented in Fall & Romanowsky (2018). We also note that, based on previous findings in Sect. 4.2, the location of the most bulge-dominated galaxies is consistent with objects found in high-density environment and with the so-called "backsplash" galaxies discussed in Sect. 4.2.2, though without a comparable sample of field galaxies in this mass range it is difficult to assess whether this is environmentally driven.

### 4.3.2. Effect of co-rotating bulges

Including bulges in the analysis of the angular momentum can be done by integrating them in the normalisation, as discussed in Sect. 4.3.1, and/or by including their contribution to the HST flux map used to compute the angular momentum. For the former, the net effect is to increase the normalisation and thus to reduce the angular momentum (for instance, see Fig. 9), whereas for the latter the expected net effect is to increase the flux (mostly in the inner parts) and thus to increase the angular momentum. If we consider that bulges are dispersion dominated systems, then it is more appropriate to remove their contribution to the HST image before deriving the galaxies’ stellar angular momentum, as has been done so far. However, some simulations seem to suggest that the stellar mass may be partially redistributed from the disk to the bulge component after a wet merger (e.g. Lagos et al. 2018), producing a rotating bulge. If this scenario is correct, then it is interesting to compare the location of the bulge-dominated galaxies when including the bulge component in our derivation of the angular momentum, assuming it is co-rotating with the disk. To do so, we derive the angular momentum using the complete HST maps (i.e. without removing the bulge) and using the mass models rotation curve. When doing so, bulge-dominated galaxies are actually brought back along the Fall relation of disk-dominated galaxies when using the disk mass as normalisation, as illustrated by the dotted line in Fig. 10, and closer, though mostly still below, when using the total stellar mass as normalisation (dashed-dotted line in Fig. 10). Overall, there is a good agreement between disk- and bulge-dominated galaxies when using either normalisation, except at intermediate stellar mass where there is still a deficit of angular momentum in bulge-dominated galaxies compared to the disk-dominated ones. The difference observed in the lowest stellar mass bins when using the disk mass as normalisation is driven by a few outlier bulge-dominated galaxies located above the main Fall relation because the statistics are lower for these bins.

To interpret these results, let us consider two galaxies that have similar disks and DM haloes but the second galaxy has a more massive bulge component. The second galaxy will always move towards higher stellar masses in the Fall relation because we consider the total stellar mass in the relation. Then,

the second galaxy might further move depending on the relation that is considered: (i) in the Fall relation that uses the disk-only flux map and the disk mass as normalisation, the second galaxy will not move along $j_*$ because the bulge does not intervene. For bulge-dominated galaxies (B/D > 1), this means a stellar mass offset of at least 0.3 dex which brings back the most massive bulge-dominated galaxies to a similar Fall relation as the disk-dominated ones. (ii) The second galaxy will have a lower $j_*$ because the bulge contributes only to the normalisation when using the disk-only flux map and the total stellar mass for the normalisation of the Fall relation. For bulge-dominated galaxies, this will produce an offset in $j_*$ of at least -0.3 dex which is what is seen in Fig. 10. (iii) For the Fall relation that uses the total flux map (i.e. including the bulge flux) but the disk mass as normalisation, the second galaxy will move towards higher values of $j_*$ because the flux will be higher. The expected offset cannot be derived as simply as before because the impact on the angular momentum will depend on the shapes and amplitudes of the rotation curve and of the bulge component. Still, as a crude estimate, if we double the mass of the disk component, which should give us a sort of upper limit, it will produce the same offset of 0.3 dex. (iv) Finally, when we take the contribution of the bulge into account both in the normalisation and in the flux map, then the two effects should, at least partially, compensate each other. This is also what we see in Fig. 10. Hence, depending on the mass range probed, the variation in zero point between disk- and bulge-dominated galaxies when varying the way the Fall relation is derived can be explained. For massive bulge-dominated galaxies, the picture seems consistent with galaxies having their angular momentum nearly entirely in their disk component since, based on the arguments presented above, including the bulge as a co-rotating component over-predicts the angular momentum compared to disk-dominated galaxies at a similar disk’s mass. For lower mass galaxies, the picture seems different. This is particularly true at intermediate stellar masses ($M_* \sim 10^{9.5} - 10^{10} M_\odot$) where even co-rotating bulges seem to under-predict compared to disk-dominated galaxies. However, contrary to high-mass galaxies, these objects would have been removed for the majority by our selection and so necessarily have a higher intrinsic uncertainty on their angular momentum. Hence, this difference may not be completely significant for this population of galaxies.

## 5. Conclusions

In this analysis, we have studied the environmental dependence of the angular momentum-stellar mass (Fall) relation for a sample of nearly 200 intermediate redshift ($0.2 \le z \le 1.5$) galaxies from the MAGIC survey that are found either in the field or in structures of varying density. These galaxies were chosen from a larger sample of [O II] emitters in MAGIC using robust selection criteria. We have modelled the galaxies morphology using their HST F814W images and their ionised gas kinematics from their MUSE cubes using the [O II] doublet as kinematics tracer. We develop a new method to estimate the galaxies’ stellar angular momentum that combines our robust rotation curves obtained from MUSE with high-resolution HST images. This approach allows us to alleviate the assumption of axisymmetry for the disk component and therefore to derive more precise values of the angular momentum.

After discussing the impact of the sample selection and assessing the reliability of this method compared to simpler approaches to measure the angular momentum, we discuss the impact of the environment and of the galaxies’ bulge component on
the Fall relation. First, we fit the Fall relation, separating galaxies between those in the field and those in rich structures, and we find a significant 0.12 dex offset that reduces to 0.08 dex when applying mass and redshift cuts. This result indicates that galaxies in the richest structures are depleted in angular momentum by roughly 20% compared to those in the field at the same redshift and in the same stellar mass range, consistent with other studies (e.g. Pelliccia et al. 2019). There is no clear indication that this result is correlated to the scatter found in the relation. However, it is more significant when massive galaxies are taken into account. Because such galaxies are predominantly found in structures, we argue that it may be an indication that they are more impacted by the environment than their lower mass counterparts. This result is further strengthened when comparing the location of the galaxies on the Fall relation as a function of various environmental tracers. Using global and local density estimates, we find that massive galaxies are indeed found below the Fall relation observed at lower stellar masses and that these objects are associated to the population of so-called "backsplash" galaxies and to the highest density parts of their host structure. We also show that the shape of the Fall relation and the location of bulge-dominated galaxies can significantly change when normalising by either the disk or the total stellar mass and when taking the contribution of the bulge component in the flux distribution into account for the derivation of the angular momentum. We find that high mass bulge dominated galaxies are consistent with a scenario where their angular momentum is nearly entirely found in their disk component. For lower mass galaxies, the lack of well constrained statistics prevent us from drawing any conclusions.

Thus, this analysis hints at the fact that there might be a visible impact of the environment on the Fall relation. However, our current lack of a complete sample of massive galaxies found at the same redshift prevents us from drawing clear conclusions. It may be that massive galaxies in structures suffer from a depletion of angular momentum, but without comparable field galaxies, we cannot eliminate the other possibility that this effect is mass-driven. To provide a definite answer, we will need to assemble a large sample of massive field galaxies from similar MUSE surveys that probe an even much larger volume than MAGIC currently does.

*Acknowledgements.* We dedicate this article in memory of Hayley Finley. This work was supported by the Programme National Cosmology et Galaxies (PNCG) of CNRS/INSU with INP and IN2P3, co-funded by CEA and CNES. This work has been carried out through the support of the ANR FOGHAR (ANR-13-BS05-0010-02), the OCEVU Labex (ANR-11-LABX-0060), and the A*MIDEX project (ANR-11-IDEX-0001-02), which are funded by the "Investissements d'avenir" French government program managed by the ANR. This work has been carried out thanks to the support of the Ministry of Science, Technology and Innovation of Colombia (MINCIENCIAS) PhD fellowship program No. 756-2016. This research has made use of MATPLOTLIB (Hunter 2007), scipy (Virtanen et al. 2020), NUMPY (Harris et al. 2020) and ASTROPY (Astropy Collaboration et al. 2013, 2018). We acknowledge David Carton for his investment in the build-up of the project.

## References

* Abril-Melgarejo, V., Epinat, B., Mercier, W., et al. 2021, A&A, 647, A152

* Astropy Collaboration, Price-Whelan, A. M., Sipócz, B. M., et al. 2018, AJ, 156, 123

* Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

* Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1

* Balogh, M. L., Navarro, J. F., & Morris, S. L. 2000, ApJ, 540, 113

* Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, Astronomy & Astrophysics, 622, A103

* Bottinelli, L., Gouguenheim, L., Paturel, G., & de Vaucouleurs, G. 1983, A&A, 118, 4

* Bouché, N., Finley, H., Schroetter, I., et al. 2016, The Astrophysical Journal, 820, 121

* Bouché, N., Murphy, M. T., Kacprzak, G. G., et al. 2013, Science, 341, 50

* Bouché, N. F., Genel, S., Pellissier, A., et al. 2021, Astronomy & Astrophysics, 654, A49

* Brook, C. B., Stinson, G., Gibson, B. K., et al. 2012, Monthly Notices of the Royal Astronomical Society, 419, 771

* Bruzual, G. & Charlot, S. 2003, Monthly Notices of the Royal Astronomical Society, 344, 1000

* Burkert, A., Förster Schreiber, N. M., Genzel, R., et al. 2016, The Astrophysical Journal, 826, 214

* Cadiou, C., Dubois, Y., & Pichon, C. 2022, Monthly Notices of the Royal Astronomical Society, 514, 5429

* Cadiou, C., Pontzen, A., & Peiris, H. V. 2022, Stellar angular momentum can be controlled from cosmological initial conditions, arXiv:2206.11913 [astro-ph]

* Cappellari, M., Scott, N., Alatalo, K., et al. 2013, MNRAS, 432, 1709

* Charlot, S. & Fall, S. M. 2000, The Astrophysical Journal, 539, 718

* Chemin, L., Huré, J.-M., Soubiran, C., et al. 2016, Astronomy & Astrophysics, 588, A48

* Ciesla, L., Buat, V., Boquien, M., et al. 2021, Astronomy & Astrophysics, 653, A6

* Ciesla, L., Elbaz, D., Schreiber, C., Daddi, E., & Wang, T. 2018, Astronomy & Astrophysics, 615, A61

* Contini, T., Epinat, B., Bouché, N., et al. 2016, A&A, 591, A49

* Cortese, L., Fogarty, L. M. R., Bekki, K., et al. 2016, MNRAS, 463, 170

* Cortese, L., Fogarty, L. M. R., Bekki, K., et al. 2016, Monthly Notices of the Royal Astronomical Society, 463, 170

* Danovich, M., Dekel, A., Hahn, O., Ceverino, D., & Primack, J. 2015, Monthly Notices of the Royal Astronomical Society, 449, 2087

* DeFelippis, D., Genel, S., Bryan, G. L., & Fall, S. M. 2017, The Astrophysical Journal, 841, 16

* del Pino, A., Fardal, M. A., van der Marel, R. P., et al. 2021, The Astrophysical Journal, 908, 244

* Di Teodoro, E. M., Posti, L., Ogle, P. M., Fall, S. M., & Jarrett, T. 2021, Monthly Notices of the Royal Astronomical Society, 507, 5820

* Doroshkevich, A. G. 1970, Astrophizika, 6, 581

* Epinat, B., Amram, P., Balkowski, C., & Marcelin, M. 2010a, MNRAS, 401, 2113

* Epinat, B., Amram, P., Balkowski, C., & Marcelin, M. 2010b, Monthly Notices of the Royal Astronomical Society, 401, 2113

* Epinat, B., Contini, T., Le Fèvre, O., et al. 2009, Astronomy & Astrophysics, 504, 789

* Fall, S. M. 1983, in Internal Kinematics and Dynamics of Galaxies, ed. E. Athanassoula, Vol. 100, 391-398

* Fall, S. M. & Efstathiou, G. 1980, Monthly Notices of the Royal Astronomical Society, 193, 189

* Fall, S. M. & Romanowsky, A. J. 2018, The Astrophysical Journal, 868, 133

* Genel, S., Fall, S. M., Hernquist, L., et al. 2015, The Astrophysical Journal Letters, 804, L40

* Gilbank, D. G., Baldry, I. K., Balogh, M. L., Glazebrook, K., & Bower, R. G. 2010, MNRAS, no

* Gilbank, D. G., Baldry, I. K., Balogh, M. L., Glazebrook, K., & Bower, R. G. 2011, MNRAS, 412, 2111

* Gill, S. P. D., Knebe, A., & Gibson, B. K. 2005, Monthly Notices of the Royal Astronomical Society, 356, 1327

* Graham, A. W., Driver, S. P., Petrosian, V., et al. 2005, AJ, 130, 1535

* Guérou, A., Krajnović, D., Epinat, B., et al. 2017, A&A, 608, A5

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

* Harrison, C. M., Johnson, H. L., Swinbank, A. M., et al. 2017, MNRAS, 467, 1965

* Heidmann, J., Heidmann, N., & de Vaucouleurs, G. 1972, Memoirs of the Royal Astronomical Society, 75, 85

* Hetznecker, H. & Burkert, A. 2006, Monthly Notices of the Royal Astronomical Society, 370, 1905

* Hung, D., Lemaux, B. C., Gal, R. R., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 5524

* Hung, D., Lemaux, B. C., Gal, R. R., et al. 2021, Monthly Notices of the Royal Astronomical Society, 502, 3942

* Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90

* Iovino, A., Petropoulou, V., Scodeggio, M., et al. 2016, Astronomy & Astrophysics, 592, A78

* Jiang, F., Dekel, A., Kneller, O., et al. 2019, Monthly Notices of the Royal Astronomical Society, 488, 4801

* Kennicutt, Jr., R. C. 1992, ApJ, 388, 310

* Kennicutt, Jr., R. C. 1998, ARA&A, 36, 189

* Knobel, C., Lilly, S. J., Iovino, A., et al. 2012, The Astrophysical Journal, 753, 121
Koekemoer, A. M., Aussel, H., Calzetti, D., et al. 2007, ApJS, 172, 196

Lagos, C. d. P., Stevens, A. R. H., Bower, R. G., et al. 2018, MNRAS, 473, 4956

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, The Astrophysical Journal Supplement Series, 224, 24

Le Fèvre, O., Tasca, L. A. M., Cassata, P., et al. 2015, Astronomy & Astrophysics, 576, A79

Lemaux, B. C., Cucciai, O., Le Fèvre, O., et al. 2022, Astronomy & Astrophysics, 662, A33

Lemaux, B. C., Tomczak, A. R., Lubin, L. M., et al. 2017, Monthly Notices of the Royal Astronomical Society, 472, 419

Lilly, S. J., Fevre, O. L., Renzini, A., et al. 2007, The Astrophysical Journal Supplement Series, 172, 70

Lin, L., Cooper, M. C., Jian, H.-Y., et al. 2010, ApJ, 718, 1158

Mancera Piña, P. E., Posti, L., Fraternali, F., Adams, E. A. K., & Oosterloo, T. 2021a, Astronomy & Astrophysics, 647, A76

Mancera Piña, P. E., Posti, L., Pezzulli, G., et al. 2021b, Astronomy & Astrophysics, 651, L15

Massey, R., Stoughton, C., Leauthaud, A., et al. 2010, MNRAS, 401, 371

Mercier, W., Epinat, B., Contini, T., et al. 2022, arXiv:2204.08724 [astro-ph], arXiv: 2204.08724

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Noble, A. G., Webb, T. M. A., Muzzin, A., et al. 2013, The Astrophysical Journal, 768, 118

Obreschkow, D. & Glazebrook, K. 2014, The Astrophysical Journal, 784, 26

Peebles, P. J. E. 1969, Astrophysical Journal, 155, 393

Pelliccia, D., Lemaux, B. C., Tomczak, A. R., et al. 2019, Monthly Notices of the Royal Astronomical Society, 482, 3514

Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, The Astronomical Journal, 124, 266

Pérez-Martínez, J. M., Ziegler, B., Dannerbauer, H., et al. 2021, A&A, 646, A53

Posti, L., Fraternali, F., Di Teodoro, E. M., & Pezzulli, G. 2018, Astronomy & Astrophysics, 612, L6

Rizzo, F., Fraternali, F., & Iorio, G. 2018, MNRAS, 476, 2137

Romanowsky, A. J. & Fall, S. M. 2012, The Astrophysical Journal Supplement, 203, 17

Romanowsky, A. J. & Fall, S. M. 2012, The Astrophysical Journal Supplement Series, 203, 17

Rousseeuw, P. & Driessen, K. 2006, Data Min. Knowl. Discov., 12, 29

Salpeter, E. E. 1955, The Astrophysical Journal, 121, 161

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Sérsic, J. L. 1963, Boletín de la Asociación Argentina de Astronomía La Plata Argentina, 6, 41

Stewart, K. R., Brooks, A. M., Bullock, J. S., et al. 2013, The Astrophysical Journal, 769, 74

Swinbank, A. M., Harrison, C. M., Trayford, J., et al. 2017, MNRAS, 467, 3140

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261

White, S. D. M. 1984, Astrophysical Journal, 286, 38

Wright, S. A., Larkin, J. E., Barczys, M., et al. 2007, The Astrophysical Journal, 658, 78

Zabl, J., Bouché, N. F., Schroetter, I., et al. 2019, Monthly Notices of the Royal Astronomical Society, 485, 1961

Zavala, J., Frenk, C. S., Bower, R., et al. 2016, Monthly Notices of the Royal Astronomical Society, 460, 4466
## Appendix A: Additional figures

<figure><img src="image_72.png" /><figcaption>**Fig. A.1.** Fall relation between galaxies in the field (black points) and those found in the richest structures (more than ten members - red circles) using the flat model rotation curve. See Fig. 5 for details regarding the legend.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>**Fig. A.3.** Fall relation using the flat model rotation curve colour coded according to the bulge-to-disk ratio (B/D) evaluated at the global effective radius for the entire kinematics sample. The stellar mass, angular momentum, and its normalisation are all measured within $R_{22}$, and the normalisation corresponds to the mass of the disk only. See Fig. 9 for the details of the legend.</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>**Fig. A.2.** Fall relation colour coded according to the global environment tracer $\eta$ using the flat model rotation curve. See Fig. 7 regarding the details of the legend.</figcaption></figure>
## Appendix B: Deriving angular momentum

Our goal is to derive the angular momentum of the stars in the MAGIC survey. We can split distribution of stars into two main components: the stellar disk and the stellar bulge. Among the two, the disk is the simplest to derive its angular momentum because of the assumptions that are usually made on its geometry and on its kinematics. Besides, bulges are usually thought of as stellar systems dynamically supported by random motions. In our case, because we model the bulge component as a spherical system, this means that its angular momentum should be nearly null. Thus, we will make the assumption that the bulge component does not possess any angular momentum and we will focus on deriving that of the disk component only.

### Appendix B.1: General derivation

The angular momentum of a system is defined from its 3D mass distribution $\rho(r)$ and 3D velocity field $V(r)$ as

$$J = \int_V d\tau \rho(r) r \times V(r), \qquad (B.1)$$

where $V$ represents the volume over which the integral is computed. In cylindrical coordinates, the cross product writes

$$r \times V(r) = RV_{\theta}\hat{z} + (zV_R - RV_{\theta})\hat{\theta} - zV_{\theta}\hat{R}, \qquad (B.2)$$

where $V(r) = V_R \hat{R} + V_{\theta} \hat{\theta} + V_z \hat{z}^6$. Assuming we want to derive the angular momentum of a disk-like distribution where only the tangential component of the velocity field is non-zero and does not depend on $z$, that is $V_R = V_z = 0$ and $V_{\theta} = V_{\theta}(R, \theta)$, then the angular momentum simplifies to

$$J = \int_V d\theta \, dR \, dz \, \rho(r) [R^2 V_{\theta}(R, \theta) \hat{z} - Rz V_{\theta}(R, \theta) \hat{R}]. \qquad (B.3)$$

### Appendix B.1.1: Angular momentum of thick disks

The angular momentum of a disk galaxy with a non-zero thickness profile $f(z)$ is similar to that of a razor-thin disk galaxy with the same surface mass density $\Sigma_M(R)$ if the three following criteria are met:

(i) $\rho(r) = \Sigma_M(R, \theta)f(z)$,

(ii) $\int_R dz \, f(z) = 1$,

(iii) $f(z) = f(-z)$.

Criterion (i) means that it is possible to separate the surface mass density profile from the thickness profile. For instance, this is the case for a double exponential profile where $f(z) = \exp\{-|z/h_z|/h_z\}$. Among the three conditions, criterion (ii) is always met since it is required in order to recover the surface mass density when integrating the 3D mass distribution along the vertical axis with respect to the disk plane, that is $\Sigma_M(R, \theta) = \int_R dz \, \rho(r)$. Finally, criterion (iii) is usually assumed to simplify the shape of the vertical profile. If these three conditions are met, one can see that the rightmost integral in Eq. B.3 vanishes when integrating along the vertical axis, whereas only the $R$ and $\theta$ integrations remain in the leftmost one. Thus, we

get an angular momentum integral similar to that of a razor-thin disk<fnref n="7" />

$$J = \int_S d\theta \, dR \, R^2 \, \Sigma_M(R, \theta) V_{\theta}(R, \theta) \hat{z}, \qquad (B.4)$$

where $S$ represents the surface over which the $R$ and $\theta$ integrations are carried out. It is common practice to normalise the angular momentum by the total stellar mass to define a specific angular momentum $j = |J|/M_\star$ (e.g. Fall 1983; Romanowsky & Fall 2012; Bouché et al. 2021) but it is also possible to normalise it by the mass of the disk only. For instance, if we normalise the angular momentum measured over some surface $S$ by the total disk mass, we get

$$j = \frac{\int_S d\theta \, dR \, R^2 \, \Sigma_M(R, \theta) V_{\theta}(R, \theta)}{\int_{[0,2\pi) \times \mathbb{R}_+} d\theta \, dR \, R \, \Sigma_M(R, \theta)}, \qquad (B.5)$$

or in the case where there is no dependence on $\theta$ for both the surface mass density and the velocity:

$$j = \frac{\int_R dR \, R^2 \, \Sigma_M(R) V_{\theta}(R)}{\int_{\mathbb{R}_+} dR \, R \, \Sigma_M(R)}, \qquad (B.6)$$

### Appendix B.2: Angular momentum for special rotation curves

In most general terms, if one wants to compute the specific angular momentum for any surface mass density and velocity profile, one needs to numerically integrate Eq. B.5 or B.6 depending on whether there is a $\theta$ dependence or not. However, there exist for a few specific velocity profiles simplified formulae which can be computed either numerically or analytically in certain cases. We present in the following simplified expressions for some rotation curves, assuming no $\theta$ dependence on both the surface mass density and the rotation curve.

### Appendix B.2.1: Constant rotation curve

If one assumes a constant rotation curve $V_{\theta}(R) = V_{\theta}$, then the angular momentum up to radius $r$ is given by

$$j(r) = V_{\theta} \times M_2 \{\Sigma_M\}_{0}^{r} / M_1 \{\Sigma_M\}_{0}^{\infty}, \qquad (B.7)$$

where $M_k\{\Sigma_M\}_\alpha^\beta$ is the $k$th radial moment of the surface brightness distribution $\Sigma_M$ integrated between $\alpha$ and $\beta$, that is

$$M_k\{\Sigma_M\}_\alpha^\beta = \int_\alpha^\beta dR \, R^k \, \Sigma_M(R). \qquad (B.8)$$

### Appendix B.2.2: Flat model

If one assumes the rotation curve to be a flat model as in Abril-Melgarejo et al. (2021), that is described by the following formula:

$$V_{\theta}(R) = V_t \times \begin{cases} R/R_t, & \text{if } R < R_t \\ 1, & \text{otherwise} \end{cases} \qquad (B.9)$$

<footnote n="6">⁶ In the following, when the velocity field is fully described by $V_{\theta}$, we will refer to it as the galaxy rotation curve.</footnote>

<footnote n="7">⁷ For a razor-thin disk we have $f(z) = \delta(z)$, with $\delta$ a Dirac distribution, so criteria (i) to (iii) are also met.</footnote>
then the angular momentum is given by

$$j(r) = \frac{V_t}{R_t M_1 \{\Sigma_M\}_0^\infty} \times \begin{cases} M_3 \{\Sigma_M\}_0^R, & \text{if } R < R_t \\ M_3 \{\Sigma_M\}_0^{R_t} + R_t M_2 \{\Sigma\}_R^R, & \text{otherwise} \end{cases} \qquad (B.10)$$

## Appendix B.3: Angular momentum for a Sérsic profile

If the surface mass density can be further described by a Sérsic profile (Sérsic 1963) with Sérsic index $n$, effective radius $R_{eff}$ and central surface mass density $\Sigma_{M,0}$, then we have

$$J = \Sigma_{M,0} \int_R dR R^2 e^{-b_n(R/R_{eff})^{1/n}} \int_{\theta} d\theta V_{\theta}(R, \theta) \hat{z}, \qquad (B.11)$$

where $b_n$ is the solution of the equation $\Gamma(2n) = 2\gamma(2n, b_n)$ (Graham et al. 2005), with $\gamma$ and $\Gamma$ the lower incomplete and complete gamma functions, respectively. Without any $\theta$ dependence, the specific angular momentum writes

$$j = \frac{2 b_n^{2n}}{(2n)!} \int_R dR (R/R_{eff})^2 V_{\theta}(R) e^{-b_n(R/R_{eff})^{1/n}}, \qquad (B.12)$$

where the normalisation factor comes from the analytical expression for the total mass of a Sérsic profile (e.g. Mercier et al. 2022). In the case of a flat rotation curve as in Eq. B.9, it is possible to compute analytically Eq. B.8 and therefore Eq. B.10. Indeed, we have

$$M_k \{\Sigma_M\}_\alpha^\beta = \frac{n \Sigma_{M,0} R_{eff}^{k+1}}{b_n^{n(k+1)}} \left[ \begin{aligned} & \gamma \left( n(k+1), b_n \left( \frac{\beta}{R_{eff}} \right)^{1/n} \right) \\ & - \gamma \left( n(k+1), b_n \left( \frac{\alpha}{R_{eff}} \right)^{1/n} \right) \end{aligned} \right]. \qquad (B.13)$$

## Appendix B.4: Recovering the RF12 approximation

The Romanowsky & Fall (2012) approximation, hereafter referred as RF12, assumes an exponential disk with a constant rotation curve and computes the angular momentum at infinity. Thus, we can use Eq. B.10 and B.13 to check whether we recover it. First, we can see that for a flat model and in the limit where the rotation curve becomes constant, that is $R_t \to 0$, Eq. B.10 reduces to B.7. Furthermore, evaluated at infinity and for an exponential disk model (i.e. $n = 1$) Eq. B.13 simplifies to

$$M_k \{\Sigma_M\}_0^\infty = \Sigma_{M,0} R_d^{k+1} \Gamma(k+1), \qquad (B.14)$$

where $\Gamma$ is the complete gamma function and $R_d = R_{eff,d}/b_1$ is the disk scale length, with $b_1 \approx 1.6783$. Thus, given Eq. B.7, we do recover the Romanowsky & Fall (2012) approximation

$$j = 2V_{\theta}R_{d}. \qquad (B.15)$$

## Appendix B.5: Angular momentum for an exponential disk in equilibrium

It can be interesting to derive the angular momentum for an exponential disk model that is stable against its own gravity, that is without any DM halo or bulge component. We consider the case where the normalisation is taken at the same radius as the angular momentum. Given the expression of the integrated mass and circular velocity for such a model (see Eqs. C.1, D.8, and D.9 of Mercier et al. (2022), its integral writes

$$j(R) = \frac{\sqrt{\pi G R_d \Sigma_{M,d}(0)}}{2\pi R_d^3 \gamma(2, R/R_d)} \int_0^R dr \, r^3 e^{-r/R_d} f\left(\frac{r}{2R_d}\right), \qquad (B.16)$$

where $G$ is the gravitational constant, $R_d$ is the disk scale length, $\Sigma_{M,d}(0)$ is the central surface mass density of the disk, $\gamma$ is the lower incomplete gamma function, and $f$ is a function of Bessel functions defined in Appendix D of Mercier et al. (2022). If we define $y = x/R_d$, then it simplifies to

$$j(R) = \frac{\sqrt{GM_d(R)R_d}}{2\sqrt{2\pi}\gamma^2(2, R/R_d)} \times \int_0^{R/R_d} dy \, y^3 e^{-y} f(y/2), \qquad (B.17)$$

where $M_d(R)$ is the mass of the disk within the radius $R$. At $R_{22} = 2.2R_d$, the radius of maximum velocity for the disk, both $\gamma$ and the integral on the right-hand side become constants and Eq. B.17 further reduces to

$$j(R_{22}) [\text{kpc km s}^{-1}] \approx 3.373 \times 10^{-4} \times \sqrt{M_d(R_{22})R_d}, \qquad (B.18)$$

where $M_d(R_{22})$ is in $M_\odot$ and $R_d$ in kpc.
Chapter 8

A new method for a detailed study
of the dark matter content in
galaxies

At the end of the first year of this Thesis was discussed the possibility to produce resolved stellar mass maps from SED fitting with the aim to better constrain the distribution of DM within galaxies. Indeed, as discussed in Chapter 5, we can use MUSE cubes to extract the galaxies' kinematics maps which can then be used to model the rotation of their ionised gas. If mass models are used, as was done in Mercier et al. (2022) and in Chapter 6, we can then constrain the fraction of DM found in these galaxies. However, deriving precise constraints on the DM fraction is difficult. Among the different sources of uncertainty that can propagate, an important one comes from the distribution of the stellar mass between the disk and the bulge components. Indeed, based on the implementation of the mass models for the analyses of the scaling relations and of the angular momentum presented in Chapter 7, there are three components that contribute to the total rotation curve: (i) the stellar bulge, (ii) the stellar disk, and (iii) the DM halo. To be complete we should also take into account the contribution of the gas, especially given that at intermediate to high redshift its content in galaxies is larger than in the local Universe (e.g. Padmanabhan & Kulkarni 2017; Wang et al. 2022). However, we do not have constraints on the amount of gas at these redshifts. An important point is that, when modelling the ionised gas kinematics, the DM halo is what is truly fitted<fnref n="1" />, being unknown, and the stellar components are fixed based on the morphological models as discussed in Sect. 5.3.3. Thus, to properly constrain the kinematics parameters of the DM halo and therefore its mass distribution we need first to have robust constraints on the distribution of mass in the stellar disk and bulge.

In Mercier et al. (2022), the disk and bulge morphological parameters were derived from a single HST map in the F814W band. These parameters (more specifically each component’s central surface brightness) were then converted in units of mass assuming a constant M/L throughout the galaxy. This M/L was derived by dividing the SED-based stellar mass with the total flux integrated in the same aperture as where the mass was evaluated. Thus, there are two potentially

<footnote n="1">¹Because we do not constrain the gas, it technically corresponds to the combination of the DM halo and the gas distribution.</footnote>
strong assumptions that had to be made: (i) that the F814W band traces well enough the distribution of mass and (ii) that the M/L is constant throughout the galaxies. The first assumption is not completely correct because at the redshifts we are concerned with (0.2 $\le$ z $\le$ 1.5) the central wavelength of this band corresponds to a rest-frame wavelength that falls between roughly 3300 Å and 6800 Å. Hence, the F814W band misses the red part of the spectrum that traces best the stellar mass. This is certainly one of the most important limitations of the modelling that I performed in Mercier et al. (2022) and it was a strong motivation to develop a framework to produce resolved stellar mass maps. The second assumption is certainly not correct as well but its effect will rather be on the relative amount of stellar mass found in the stellar disk with respect to the bulge. Indeed, with a constant M/L we implicitly assume that both components can be described by the same stellar population. However, we expect the bulge component to be redder than the disk because it is supposed to host a majority of older low-mass stars (e.g. Quilley & de Lapparent 2022) whose M/L is higher. Thus, with a constant M/L and if the disk dominates the total flux, in particular the massive young stars if there is star formation in the disk, we might underestimate the amount of mass located in the bulge. As an indication, this is typically what is seen when producing resolved stellar mass maps of spiral galaxies as described in Sect. 8.1.2.3, with M/L values that can vary between the stellar bulge and disk components by roughly a factor of five on average.

A first solution that could not be applied in Mercier et al. (2022) because we lacked high-resolution infrared bands but that could be done with future data would be to use two different bands to estimate two M/L values, one for the disk and one for the bulge (e.g. Junais et al. 2020). Though, this would still assume that the two bands would trace sufficiently well each component separately. But, if there are not just two but multiple bands all the way through from the UV to the IR then it would be even more efficient to use the full information contained in each of these bands to derive the stellar mass in each pixel separately rather than trying to estimate the M/L directly. This is the option that we decided to investigate and it will be referred to in the next parts as resolved or pixel-per-pixel SED fitting. Besides, such a method does not reduce to the sole estimation of the stellar mass in each pixel. Indeed, since it involves SED fitting, it is technically speaking possible to produce resolved maps for any physical parameters that comes out of the SED modelling. Notably, this applies for the SFR for which it might also be interesting to know its spatial distribution to improve the kinematics modelling. Producing resolved stellar mass maps is not completely new. For instance, early works did so on local galaxies by parametrising the M/L in a given band as a function of multiple colours (e.g. Zibetti et al. 2009) or by directly building SED templates and then applying them on the pixel level (e.g. Welikala et al. 2008). Follow-up studies also applied the same technique to higher redshift galaxies (e.g. Welikala et al. 2011; Wuyts et al. 2012; Lang et al. 2014; Tacchella et al. 2018), while other works rather focussed on constraining from multi-band photometry 1D M/L radial profiles which can then be used to reconstruct resolved stellar mass maps (e.g. Suess et al. 2019).

## 8.1 Resolved maps from multi-band HST photometry

The first step of this project was to select a small sample of galaxies that would be used to test the method and compare the results of the morphological and kinematics models when using the resolved stellar mass maps with the ones derived from a single band and a constant M/L value. The two important constraints to select this sample were (i) to have multiple high-resolution bands from the UV to the IR and (ii) to have well resolved morphological and
kinematics maps for these galaxies, in the sense that we should have enough spaxels in both cases to test the impact of the mass maps on the morphological and kinematics model outputs.
We decided to use a sample of 12 galaxies, nine of which are from the MXDF survey (Bacon et al. 2021) and the remaining three are located in the MOSAIC of the MUSE-HUDF survey (Bacon et al. 2017). The nine galaxies of the MXDF were selected in Bouché et al. (2022) from the [OII] emitters catalogue of Inami et al. (2017) keeping those with the highest S/N and eliminating face-on objects to properly constrain their kinematics. In Bouché et al. (2022), these galaxies were used to model their morphology with GALFIT and their kinematics with the GALPAK<sup>3D</sup> code in order to constrain the shape of their DM halo profile in the inner parts (cusp-core problem). Thus, they were specifically chosen to have well resolved morphology and kinematics. Besides, these galaxies are all located in the MXDF which is itself found within the HUDF for which the nine following HST photometric bands are available (from UV to NIR): F435W, F606W, F775W, F814W, F850LP, 105W, 125W, 140W, and 160W. Additionally, since we had contributed to the analysis performed in Bouché et al. (2022), especially on the morphological modelling of the nine MXDF galaxies with GALFIT, it further made sense to select these galaxies as part of the test sample. The main characteristics of this sample are presented in Table 8.1 with their (1) ID, (2) MUSE spectroscopic redshift, (3) total exposure time, (4) [OII] S/N, (5) total stellar mass, (6) SFR derived from SED fitting (Bouché et al. 2022) using the code Multi-wavelength Analysis of Galaxy Physical Properties (MAGPHYS)<fnref n="2" />, (7) stellar disk inclination derived from GALFIT, (8) B/T evaluated at two effective radii (based on the GALFIT best-fit model), and (9) global effective radius derived from a single Sérsic profile fit on their F160W HST map. The HST images for the 12 galaxies in the different HST bands are also represented in Figs. 8.1 and 8.2. The largest galaxy is ID1 and it has a quite massive bulge in the central parts (roughly 1.5 × 10¹⁰ M<sub>☉</sub> at two effective radii, that is around 35% of the total stellar mass; see Table 8.1 for more details). Thus, it will be used as the main test galaxy in what follows, especially in Sect. 8.2 for the machine learning application. We see that there is a non-negligible fraction of edge-on systems in the test sample and galaxies with various bulge contributions, as is also shown in Table 8.1. Among the galaxies there are two special cases with nearby companions: ID3 and ID943. For both galaxies the companions have been masked using segmentation maps (see Sect. 8.1.2.1 for more details). Most of the galaxies in the sample are quite red with HST images in the reddest band (F160W) brighter by at least 1 mag (galaxy ID37) and by at most 4 mag (galaxy ID 15) than in the bluest band (F435W). Looking at galaxy ID1 we can also see that the bluest bands (435W and 606W) tend to trace the spiral arms and clumps where star formation is most likely to happen. On the other hand, these features tend to fade in the IR bands where the bulk of the stellar mass is probed.

In order to perform the SED fitting and extract maps of physical parameters, we relied on two codes: LEPHARE and CIGALE. Among the two, CIGALE is the most recent and probably the most versatile. Initially, LEPHARE was designed to derive photometric redshifts but it is also possible to use it with a fixed redshift to estimate physical parameters. Broadly speaking, the two codes work in the same way and include similar models. The only difference is that CIGALE is an energy-balance code, meaning that the energy of the stars and gas absorbed by dust in the UV and in the optical is re-emitted entirely by the dust model in the mid-infrared (MIR) and far-infrared (FIR). The development of the framework to produce resolved stellar mass and SFR maps was started in 2021 as part of H. Plombat's Master project that I supervised. In addition, previous data processing of the 12 selected galaxies to align images and match their pixel scale

<footnote n="2">²http://www.iap.fr/magphys/</footnote>
<table><caption>Table 8.1: Main characteristics of the test sample used for resolved SED fitting. The table is adapted from Bouché et al. (2022) and H. Plombat’s Master project.</caption><thead><tr><th>ID</th><th>z</th><th>t<sub>exp</sub></th><th>S/N</th><th>log<sub>10</sub> M<sub>*</sub></th><th>log<sub>10</sub> SFR</th><th>i</th><th>B/T</th><th>R<sub>eff</sub></th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th><th>(4)</th><th>M<sub>☉</sub></th><th>M<sub>☉</sub> yr<sup>-1</sup></th><th>°</th><th>%</th><th>"</th></tr><tr><th></th><th></th><th></th><th></th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th></tr></thead><tbody><tr><td colspan="9" style="text-align: center; font-weight: bold;">MXDF</td></tr><tr><td>3</td><td>0.62</td><td>42</td><td>15</td><td>10.08<sup>+0.02</sup><sub>-0.06</sub></td><td>2.62<sup>+0.80</sup><sub>-1.20</sub></td><td>64</td><td>4</td><td>0.83</td></tr><tr><td>15</td><td>0.67</td><td>136</td><td>9</td><td>10.23<sup>+0.02</sup><sub>-0.10</sub></td><td>0.99<sup>+2.53</sup><sub>-0.99</sub></td><td>86</td><td>0</td><td>0.82</td></tr><tr><td>37</td><td>0.98</td><td>136</td><td>61</td><td>8.87<sup>+0.06</sup><sub>-0.17</sub></td><td>0.80<sup>+1.22</sup><sub>-0.18</sub></td><td>56</td><td>0</td><td>0.45</td></tr><tr><td>912</td><td>0.62</td><td>136</td><td>124</td><td>9.19<sup>+0.19</sup><sub>-0.06</sub></td><td>1.29<sup>+1.60</sup><sub>-0.30</sub></td><td>27</td><td>0</td><td>0.28</td></tr><tr><td>919</td><td>1.10</td><td>136</td><td>98</td><td>9.83<sup>+0.02</sup><sub>-0.06</sub></td><td>4.28<sup>+0.52</sup><sub>-0.55</sub></td><td>28</td><td>21</td><td>0.40</td></tr><tr><td>937</td><td>0.73</td><td>44</td><td>18</td><td>9.35<sup>+0.08</sup><sub>-0.06</sub></td><td>1.50<sup>+0.34</sup><sub>-0.36</sub></td><td>82</td><td>16</td><td>0.74</td></tr><tr><td>943</td><td>0.66</td><td>136</td><td>32</td><td>9.28<sup>+0.06</sup><sub>-0.06</sub></td><td>1.00<sup>+0.57</sup><sub>-0.40</sub></td><td>77</td><td>0</td><td>0.73</td></tr><tr><td>982</td><td>1.10</td><td>92</td><td>22</td><td>9.75<sup>+0.09</sup><sub>-0.07</sub></td><td>3.25<sup>+0.84</sup><sub>-0.97</sub></td><td>70</td><td>43</td><td>0.59</td></tr><tr><td>1002</td><td>0.99</td><td>47</td><td>35</td><td>9.43<sup>+0.02</sup><sub>-0.07</sub></td><td>1.18<sup>+0.02</sup><sub>-0.26</sub></td><td>73</td><td>21</td><td>0.52</td></tr><tr><td colspan="9" style="text-align: center; font-weight: bold;">mosaic</td></tr><tr><td>1</td><td>0.62</td><td></td><td></td><td>10.65<sup>+0.06</sup><sub>-0.08</sub></td><td>8.19<sup>+3.07</sup><sub>-2.34</sub></td><td>36</td><td>35</td><td>1.45</td></tr><tr><td>7</td><td>0.62</td><td></td><td></td><td>10.05<sup>+0.05</sup><sub>-0.04</sub></td><td>5.51<sup>+1.17</sup><sub>-1.08</sub></td><td>40</td><td>26</td><td>0.55</td></tr><tr><td>13</td><td>1.00</td><td></td><td></td><td>9.92<sup>+0.06</sup><sub>-0.06</sub></td><td>2.63<sup>+0.90</sup><sub>-0.69</sub></td><td>32</td><td>0</td><td>0.39</td></tr></tbody></table>

**Notes:** (1) ID from Inami et al. (2017), (2) MUSE spectroscopic redshift, (3) exposure time, (4) signal-to-noise ratio for [O II], (5) total stellar mass, (6) star formation rate, (7) disk inclination, (8) bulge-to-total flux ratio in the F160W HST band evaluated at two effective radii, and (9) global effective radius. The MUSE data for galaxy ID3 are from the MOSAIC of the MUSE-HUDF because it lies on the side of the MXDF FoV so that the deepest data are in the MOSAIC.

and their resolution was performed by J. C. B. Pineda. Afterwards, I kept working on the project and developed a proper python implementation.

In what follows, I will discuss the framework that we developed, starting in Sect. 8.1.1 by mentioning the fundamental principles that are shared between SED fitting codes. The properties of multi-band photometric observations and how they are handled by such codes will be discussed and, in Sect. 8.1.2, I will mention the pre-processing steps performed by J. C. B. Pineda, as well as how the HST maps errors were handled for the SED fitting part. Afterwards, a comparison will be made between the stellar mass maps reconstructed from resolved SED fitting and the distribution of mass inferred from single-band observations assuming a constant M/L throughout the galaxies. Finally, the python implementation of the framework will be discussed in Sect. 8.1.3. In particular,
examples of automatic pixel-per-pixel SED fitting modellings performed with both LEPHARE and CIGALE will be shown and differences between the two implementations will be mentioned.
<figure><img src="image_75.png" /><figcaption>Figure 8.1: HST images in the nine available bands for the first six galaxies in the test sample used for resolved SED fitting. For each galaxy the same logarithmic scale is used for the nine bands. All of the images are shown in their coarsest pixel scale (0.06" per pixel) but in their original resolution and have already been masked using segmentation maps from SExTRACTOR.</figcaption></figure>
<figure><img src="image_76.png" /><figcaption>Figure 8.2: HST images for the remaining six galaxies in the test sample used for resolved SED fitting. See Fig. 8.1 for the description.</figcaption></figure>
8.1.1 Fundamental principles of SED fitting

Before delving into the implementation we used to perform pixel-per-pixel SED fitting, I must
first discuss the basic principles behind the concept of SED fitting and how it is possible to derive
physical parameters such as stellar masses or SFRs from multi-band photometric observations.
Because the two codes we are interested in that are LEPHARE and CIGALE only handle
photometric data, that is fluxes measured in filters, and not spectra nor a combination of both, I
will solely focus on the implementation of such codes.

**8.1.1.1 Multi-band photometric observations**

<figure><img src="image_77.png" /></figure>

Figure 8.3: Illustration of the total throughput (filter, mirrors, and instrument) for the nine HST filters used to derive resolved stellar mass maps. The lines are approximately coloured by that of their central wavelength, except for the infrared filters which are coloured by shades of grey (the lighter the longer the central wavelength). Filters corresponding to the ACS WFC instrument are shown with solid lines whereas those from the WFC3 instrument in its IR channel are shown with dashed lines.

The input data that handle SED fitting codes are multi-band photometric observations, that is fluxes or magnitudes measured in a given band. By nature, the constituents of a galaxy (e.g. stars, gas, dust, etc.) emit and absorb (e.g. gas, dust) light in a continuum of wavelengths which is what is referred to as a spectrum. The spectrum of a galaxy is therefore the combination of the spectra of each of its constituents, for instance of all its stars of different ages and metallicity, of the different phases of gas (atomic, molecular, or ionised), or of the dust. In addition, these components can interact with each other. Typically, the dust will partially obscure the optical part of the spectrum emitted by the stars and the gas and will re-emit it in the IR. Thus, a photometric observation in a single band is no more than integrating the redshifted spectrum of the galaxy in the spectral window of the band.

This picture would be correct if there was no loss of light along the optical path all the way
through to the detector. However, because of the optics of the telescope, the instrument, the
detector, etc., some light is lost along the way. This is what is called throughput: the amount of
light that is preserved when passing through some part or the entirety of the system. In addition, the filters themselves are not perfect boxcar filters that let through 100% of the light in their targeted wavelength range. Thus, depending on the filter and the instrument used the spectrum of the galaxy will be more or less impacted as a function of wavelength which will change its shape and therefore the flux measured in a given band. This is illustrated in Fig. 8.3 where I show the total throughput, that is including the filter, the mirrors, and the instrument, for the nine HST filters that we used to produce resolved maps. Those from the Hubble Space Telescope Advanced Camera for Surveys (HST-ACS) that are located in the optical are represented with a color that roughly matches their wavelength whereas the IR filters from the Hubble Space Telescope Wide Field Camera 3 (HST-WFC3) are shown with varying shades of grey and with a dashed line. Hence, a first key ingredient in SED fitting codes is to provide, if not already present, the total throughput as a function of wavelength for each band that is going to be used in the modelling.

<figure><img src="image_78.png" /><figcaption>Figure 8.4: Illustration of the impact of redshifting and integrating a spectrum in various HST filters on a galaxy’s SED. The grey line represents a fictitious rest-frame spectrum at z = 0.7 adapted from the MUSE spectrum of galaxy 328_CGR84 which, for readability, has been offset by 2800 Å. The red line represents the redshifted spectrum as it would be observed and the blue data points show the integrated flux in the five HST ACS WFC filters that we used (flux scale is shown in blue on the right hand-side). For illustration, the various filters’ total throughput curves are shown as faint dotted line in the background.</figcaption></figure>

So, to convert a high-resolution spectrum into a SED, that is a combination of observed fluxes in multiple bands, at least two steps must be taken: (i) take into account the impact of cosmology by redshifting the spectrum to the galaxy’s cosmological redshift and (ii) integrate it into each band taking into account the total throughput. Based on Sect. 1.1, there are two aspects to take into account to perform the first step. First, we need to scale appropriately the amplitude of the spectrum according to whether it is measured in units of frequency or wavelength. Second, we
must redshift its wavelength range, the effect of which will be to shift it redward and to dilate it. The fact that it is dilated is important because if we consider a filter with a given spectral window then it will probe a smaller rest-frame spectral range for higher redshift galaxies. I show an illustration of how a fictitious spectrum is affected by these steps in Fig. 8.4. The spectrum has been adapted from that of galaxy 328_CGr84. In grey is represented the rest-frame spectrum offset by 2800 Å so that it is visible in the wavelength range of the redshifted spectrum (z = 0.7) shown in red. The integrated fluxes in the five HST ACS WFC bands that we used (from shorter to longer wavelengths: F435W, F606W, F775W, F814W, and F850LP) are shown with blue circle. The filters' throughput curves are also shown as faint dotted lines in the background. As discussed previously we see the effects of redshifting the spectrum, in particular how it spreads it into a larger wavelength range. We can also see that, even if the global shape of the spectrum is retained in the SED, we lose substantial details compared to the original spectrum which therefore limits the complexity of the modelling that can be applied.

Transforming a high-resolution spectrum into a SED technically requires a few more steps than those I have discussed above. Indeed, if we consider a SED that is made of space-based and ground-based observations then the latter will be affected in different ways than the former because of the impact of the atmosphere that will both absorb and emit, and that as a function of wavelength. For instance, the atmosphere will emit thermal radiation in the IR that will change throughout the night as the atmosphere cools down. Furthermore, independently of whether the observations are space-based or ground-based there are further aspects that need to be taken into account such as, for instance, the attenuation by dust in the Milky Way or the zodiacal light. The former can be modelled using a Galactic attenuation curve and an amplitude, typically taken as $A_V$ (see discussion in Sect. 8.1.1.2). Since dust is primarily located in the disk of the Galaxy, it will more severely affect observations made therein and thus a map of the spatial distribution of dust and its properties must also be taken into account (e.g. Planck Collaboration et al. 2014). The latter effect is rather due to dust located in the ecliptic plane of the Solar System that scatters light emitted from the Sun and produces a faint glow. Such additional effects may or may not be included in SED fitting codes depending on the implementation. As an example, CIGALE assumes that fluxes have already been pre-treated beforehand and therefore does not take them into account whereas LEPHARE does include a mean atmospheric extinction when computing the throughput of the filters.

### 8.1.1.2 Key concepts: IMF, SSPs, SFH, and dust attenuation

Now that we know how to transform a high-resolution spectrum into a SED, we need to actually produce models of spectra to fit the observed SEDs. This is the core concept of SED fitting codes: produce high-resolution spectra from models, degrade them to match the different bands, and fit them onto the observed SED to find the best solution. To produce such models there are at least four key ingredients that are shared between every SED fitting code: the IMF, SSP, SFH, and dust absorption. I will quickly go through each one to explain how they are used in combination with the others to produce composite models of galaxies' spectra.

SED fitting codes such as CIGALE were initially designed to model SEDs from integrated photometry, that is the SEDs combining the light coming from the entire galaxies or from most of their extension. Thus, the light that is received is the combination of a very large number of stars, dust, and gas particles whose spatial distribution is not very well constrained and it would therefore be vain to try to model each star or dust and gas particles separately. Instead, the
concept of SSPs was developed. A SSP is a population of stars that was born at the same time
from the same gas cloud. Thus, a SSP is composed of a large number of stars with different
masses but with the same age and metallicity that will evolve along an isochrone. In practice, a
SSP is stored as a spectrum that changes with time as stars evolve along and out of the stellar
main sequence. Massive blue stars will burn their fuel faster and therefore will leave the main
sequence and “die” quicker than older red stars. Thus, the SSP spectrum will tend to get redder
with cosmic time. Usually, different SSP spectra with a different evolution are provided for
discrete values of metallicity, though implementations may vary between SED fitting codes. For
instance, in CIGALE the following metallicities are available: log₁₀(O/H) = 10⁻<fnref n="4" />, 4 × 10⁻<fnref n="4" />,
4 × 10⁻<fnref n="3" />, 8 × 10⁻<fnref n="3" />, 2 × 10⁻², and 5 × 10⁻² dex. Furthermore, SSP spectra are normalised so that
their amplitude with wavelength only depends on the proportion of stars in each mass bin, which
is also an important parameter for a SSP. Indeed, depending on how many stars there are in
low-mass or high-mass bins, the shape of the spectrum will necessarily be affected. Thus, we also
need to have a function that tells us how likely it is to sample a star of a given mass from a
collapsing gas cloud. This function is called the IMF. Commonly used IMFs are those of Salpeter
(1955), Kroupa (2002), or Chabrier (2003). Once an IMF is chosen the SSP can be derived for a
given metallicity using a stellar evolution code that will sample it, evolve the stellar population
along its isochrone, and compute the total spectrum at each time step. Typical SSPs used in SED
fitting codes are those of Bruzual & Charlot (2003) and Maraston (2005). Discussing the origin of
the differences between these SSPs is beyond the scope of this discussion but I can just note that
they can be affected by lots of different treatments from, for instance, the modelling of the IMF,
the treatment of the post and pre-stellar main sequence phases, to whether they integrate the
contribution of multiple-stars systems.

However, SSPs are not sufficient in themselves to create synthetic spectra because the spectrum of
a galaxy is not just due to a single population of stars but to the sum of the spectra of all the
stellar populations in the galaxy with different ages. Because these populations have different ages
their spectrum will not evolve by the same amount and it will therefore differ in shape from one
age to another. Furthermore, this sum must be weighted by the amount of stars created at each
time step. The function that describes how many stars are born as a function of cosmic time is
referred to as the SFH. When measured at a given time step or across a certain period it is rather
called the SFR. There is a multitude of SFHs that can be used such as an exponentially declining
star-formation law, a delayed exponential law that includes a rising part before the exponential
decline, periodic, or even empirically calibrated laws (e.g. see Boquien et al. 2019). These
examples pertain to the class of parametric models, that is models for which we assume a certain
parametrised shape of the SFH. Another class of SFH models not implemented in LEPHARE
and CIGALE but that is for instance in PROSPECTOR<fnref n="3" /> (Leja et al. 2019; Johnson et al. 2021) or
STARLIGHT<fnref n="4" /> (Cid Fernandes et al. 2005) is non-parametric models. These models are
non-parametric in the sense that they do not assume a shape of the SFH or rather than they
assume a SFH that is piecewise constant. The amplitude of each piece being a free parameter,
these models are therefore usually as parametric if not more than parametric models.

Nevertheless, they can be useful to better constrain the shape of the SFH though their use can be
limited with only broad-band photometry (Leja et al. 2019). The production of a full galaxy’s
spectrum is then quite simple once we know the SFH. Indeed, let us assume that a galaxy with an
age t<sub>age</sub> can be modelled with a SSP that we will note SSP(t) (t = 0 corresponds to the beginning

<footnote n="3">³https://prospect.readthedocs.io/en/latest/index.html</footnote>

<footnote n="4">⁴http://www.starlight.ufsc.br/</footnote>
of the galaxy and $t = t_{age}$ corresponds to the time when the light that reaches us was emitted), and with a SFH that can be split in three parts: SFH(0), SFH($t_1$), and SFH($t_{age}$), with $0 < t_1 < t_{age}$. At $t = 0$, the spectrum of the galaxy will be given by $SSP(0) \times SFH(0)$. At $t = t_1$, the SSP emitted at $t = 0$ will evolve into $SSP(t_1)$ and a new SSP will emit (thus its spectrum will be that given by $SSP(0)$) with an amount of stars given by $SFH(t_1)$. Hence, the spectrum at $t = t_1$ writes $SSP(t_1)SFH(0) + SSP(0)SFH(t_1)$. Finally, at $t = t_{age}$, the first SSP will evolve to $SSP(t_{age})$, the second to $SSP(t_1)$, and a new one will emit with an amount of stars given by $SFH(t_{age})$. Hence, the final spectrum is $SSP(t_{age})SFH(0) + SSP(t_1)SFH(t_1) + SSP(0)SFH(t_{age})$. If we generalise to continuous SSPs and SFHs, then we get the general expression to compute the model galaxy's spectrum:

$$F_{\lambda} = \int_{0}^{t_{age}} dt \, \text{SSP}_{\lambda}(t)\text{SFH}(t_{age} - t), \qquad (8.1)$$

where I have used the subscript $\lambda$ for the SSP to emphasize that it is per unit wavelength as for $F_{\lambda}$ on the left-hand side of the equation.

To be complete, we need at the very least one more aspect to take into account. We need to consider the fact that the spectrum of galaxies are affected by the absorption and the scattering by dust that reduces the amount of flux we receive in the UV and in the optical and increases it in the FIR. The sheer effect of absorption and scattering by dust irrespective of the geometry of the stars, gas, and dust is called dust absorption. In essence, the more dust along the LOS, that is the higher its column density, the larger the dust absorption. On the other hand, the complex modelling of the interplay between the light of the obscured and unobscured stars, the ISM, and the dust, taking into account their spatial distribution, and potential scattering back onto the LOS is referred to as dust attenuation. This is what SED fitting codes implement. To do so, they must know the amplitude of the attenuation which is related to the column density of the dust, and its shape, that is by how much is attenuated the emitted light as a function of wavelength. The former is parametrised by the attenuation in the V-band, noted $A_V$, that was already discussed in Sect. 3.2.5.2 and the latter by the attenuation curve, noted $A_{\lambda}/A_V$. The modelling of the attenuation curve is complex and depends on the properties of the ISM, on the distribution of the stars and the grains, and on the chemistry of the grains among other things. Thus, I will not extend further and I will refer to Salim & Narayanan (2020) for a review on the topic. Two typical attenuation laws that are used are those of Calzetti et al. (2000) and Charlot & Fall (2000). These four components (IMF, SSP, SFH, and dust attenuation) are the key ingredients of SED fitting but, depending on the code and its objectives, there can be other components that might be in some cases as important. For instance, in both LEPHARE and CIGALE the contribution of emission lines can be modelled and integrated into the galaxy's spectrum model using recipes that change from one code to another. Another example in CIGALE is its treatment of the dust emission in the FIR with various dust emission models and using the principle of energy balance where the absorbed energy from the stars and the gas is used as input for the energy emitted by the dust.

### 8.1.1.3 Fitting galaxies' SEDs with LePHARE and Cigale

Both codes will find the best-fit solution in the same way. Given the models and the list of values for their parameters they will generate a multi-dimensional grid of normalised SED model. The best-fit solution is therefore the one that minimises the following $\chi^2$ expression:
$$ \chi^2 = \sum_i \left( \frac{f_i - \alpha \times m_i}{\sigma_i^2} \right)^2, \qquad (8.2) $$

where $f_i$, $m_i$, and $\sigma_i$ are the observed flux, normalised SED model, and photometric uncertainty in the band $i$, respectively, and $\alpha$ is the de-normalisation factor for the model that writes (Boquien et al. 2019)

$$ \alpha = \frac{\sum_i f_i m_i / \sigma_i^2}{\sum_i m_i^2 / \sigma_i^2}. \qquad (8.3) $$

CIGALE uses a slightly more general expression but that reduces to Eqs. 8.2 and 8.3 in our case. The scale factor $\alpha$ is then used to scale appropriately the normalised extensive physical properties (e.g. stellar mass or SFR) associated to the best-fit model. Therefore, in both codes both the stellar mass and the SFR are products of the SED fitting rather than parameters that are fitted. Finally, CIGALE implements a Bayesian method to estimate more robustly some physical parameters. For each model in the grid CIGALE computes a likelihood as $\exp(-\chi^2/2)$ and then likelihood weighted mean and standard deviation estimates for each parameter.

## 8.1.2 Application to resolved SED fitting

### 8.1.2.1 Pre-processing steps

A few pre-processing steps were applied by J. C. B. Pineda on the HST images of the 12 selected galaxies. Originally, each band had its own resolution and pixel scale. The images in the F435W, F606W, F775W, F814W, and F850LP ACS WFC bands had the finest pixel scale of 0.03" per pixel, though they were also available at 0.06" per pixel, whereas the remaining four F105W, F125W, F140W, and F160W WFC3 IR bands had a coarser scale of 0.06" per pixel. Additionally, the images' world coordinate systems (WCS) were aligned between the bands of the same instrument (ACS WFC or WFC3) but not between the two instruments. Finally, some objects had companions in their HST images that needed to be masked. Therefore there was a need to equate the resolutions and the pixel scales of the images and align them before producing resolved maps. The way this was achieved is as follows: (i) the sky background was estimated in empty regions of the images and subtracted for each image in each band, (ii) the HST PSF FWHM was derived in each band by fitting as many stars as possible, the retained value being the median FWHM of all the fits (see Table 8.2 for the full list of values), (iii) for each galaxy, East-West and North-South offsets were measured in each band using the images at a scale of 0.06" per pixel and a correction was applied to align the images, (iv) the images were convolved by a Gaussian kernel<fnref n="5" /> to match the worst resolution in our data, that is in the F160W band with a FWHM of 0.22", and (v) the images were masked using segmentation maps from Inami et al. (2017) derived with the Source-Extractor (SEXTRACTOR)<fnref n="6" /> software (Bertin & Arnouts 1996) using multiple HST bands (see also Rafelski et al. 2015).

The images were originally given in units of electrons per second so they needed to be converted in physical units before being handled by LEPHARE or CIGALE. The former accepts magnitudes in the AB system (Oke 1974) whereas the latter accepts flux densities in Jy, that is per unit

<footnote n="5">⁵A Moffat profile might have been more appropriate but a Gaussian kernel was simpler to match the resolution between each band.</footnote>

<footnote n="6">⁶https://www.astromatic.net/software/sextractor/</footnote>
<table><caption>Table 8.2: Zero point and PSF FWHM values for the nine HST bands used for the derivation of the resolved maps. The table is adapted from H. Plombat's Master project.</caption><thead><tr><th>HST band</th><th>Zero point<br>AB magnitude</th><th>PSF FWHM<br>"</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th></tr></thead><tbody><tr><td>F435W</td><td>25.68</td><td>0.110</td></tr><tr><td>F606W</td><td>26.51</td><td>0.134</td></tr><tr><td>F775W</td><td>25.69</td><td>0.123</td></tr><tr><td>F814W</td><td>25.94</td><td>0.129</td></tr><tr><td>F850LP</td><td>24.87</td><td>0.120</td></tr><tr><td>F105W</td><td>26.27</td><td>0.208</td></tr><tr><td>F125W</td><td>26.23</td><td>0.211</td></tr><tr><td>F140W</td><td>26.45</td><td>0.219</td></tr><tr><td>F160W</td><td>26.94</td><td>0.220</td></tr></tbody></table>

**Notes:** (1) HST band, (2) zero point in the AB magnitude system, and (3) HST PSF FWHM. The zero point values are those used in Eq. 8.5 to convert the HST maps from electron counts to AB magnitudes and/or flux densities.

frequency. Formally, AB magnitudes are defined as the magnitude $m_{AB}$ corresponding to an equivalent flux density $F_{\nu}$ that would produce the same observed count rate but with a flat spectrum per unit of frequency. Its expression is

$$m_{AB} = -2.5 \log_{10} F_{\nu} - 48.6, \qquad (8.4)$$

where $F_{\nu}$ must be in erg s$^{-1}$ Hz$^{-1}$ cm$^{-2}$. To convert images from electrons per second to AB magnitudes one can use the following formula:

$$m_{AB} = -2.5 \log_{10} \text{CR} + \text{Zpt}_{AB}, \qquad (8.5)$$

where CR is the count rate measured in electrons per second and Zpt<sub>AB</sub> is the zero point value that must be used to get AB magnitudes. This zero point can be derived using two keywords given in HST image headers: an inverse sensitivity noted PHOTFLAM and a pivot wavelength noted PHOTPLAM. Its expression is given by

$$\text{Zpt}_{AB} = -2.5 \log_{10} \text{PHOTFLAM} - 5 \log_{10} \text{PHOTPLAM} - 2.408, \qquad (8.6)$$

where PHOTFLAM must be in erg Å$^{-1}$ e$^{-1}$ cm$^{-2}$ and PHOTPLAM in Å. The different zero points are listed in Table 8.2 for the nine HST bands. With Eqs. 8.4, 8.5, and 8.6 it is possible to convert from count rates to fluxes and magnitudes and vice versa.

### 8.1.2.2 Handling flux uncertainties

An important aspect for any kind of fit, including SED fitting, is to properly estimate the uncertainty on the observations. This is particularly true for SED fitting because we do not want
a band with a large uncertainty to contribute too much to the estimation of the galaxy's physical parameters. This is even more true in the case of pixel-per-pixel SED fitting. Indeed, if we assume that the photons that hit the detector follow Poisson's statistics, then we know that the standard deviation for a signal S will be given by $N = \sqrt{S}$. This gives us an estimate of the uncertainty on the signal and therefore on the noise level in the absence of additional noise sources. The S/N is then equal to $\sqrt{S}$ which means the more signal the relatively more precise it becomes. Incidentally, this also means that the S/N in each individual pixel will be lower than that of the integrated signal coming from the entire galaxy and therefore that we need more precise estimates of the noise than we would with an integrated SED to have precise constraints on the stellar mass and SFR in each pixel.

To do so, we started from the weight maps given alongside the HST images that correspond to the inverse of the variance of the images, excluding the contribution of the noise due to the signal (i.e. Poisson's noise). However, these are not standard images since they have been treated through a drizzling program (Illingworth et al. 2013) called MULTIDRIZZLE and described in Koekemoer et al. (2003). This drizzling step allows to combine multiple dithered images of the same FoV to produce a final image that has a smaller pixel scale than the original images but at the cost of introducing correlations between the noise in each of the new pixels (Casertano et al. 2000). This renders the modelling of the propagation of the signal-induced noise difficult. Besides, as argued in Casertano et al. (2000), the relevant noise, for instance to decide on the significance of a detection, is rather the background noise. Hence, the noise induced by the signal was not implemented in their weight maps, which in turn produces error maps that are relatively flat. Nevertheless, even if the background noise is sufficiently high, we would expect the signal-induced noise to be larger where the flux of the galaxy peaks. Therefore, we decided to include an estimate of Poisson's noise into these maps. If the detector measures a count rate S over an exposure time $T_e$, then the signal corresponding to the detected number of electrons is $S_e = S \times T_e$. If we assume that the electrons obey Poisson's statistics, then this signal must also be equal to its variance. In addition, because HST images are given in units of electrons per second, the HST weight maps provide us with variance maps (noted $\sigma^2$) for the background signal in units of electrons squared per second squared. Thus, the total variance due to the background signal over the entire exposure time is $\sigma_b^2 = T_e^2 \times \sigma^2$. The total variance $\sigma_{tot}^2$ will be the sum of the two and the corresponding variance $\sigma_{e/s}^2$ in electrons per second will therefore write

$$ \begin{aligned} \sigma_{e/s}^2 &= \frac{\sigma_{tot}^2}{T_e^2} \\ &= \sigma^2 + \frac{S}{T_e}, \end{aligned} \qquad (8.7) $$

where S is the signal in electrons per second and $\sigma^2$ is the variance in electrons squared per second squared coming from the HST weight maps. As already discussed above, the drizzling step changes the pixels' dimensions and introduces correlations into the noise so that Eq. 8.7 is only an approximation of the real variance. A last important aspect is to properly take into account the effect of convolution. Indeed, the HST images were convolved to the worst resolution available in the F160W band (see Sect. 8.1.2.1). As long as we can neglect the covariance between pixels, the variance map must also be convolved but using the square of the kernel that was used for the images (see in Appendix D for a demonstration):
$$ \sigma^2(x_0, y_0) = \sum_{x,y} K^2(x_0 - x, y_0 - y)\sigma^2(x, y) + \sum_{\substack{x,x',y,y' \\ x \neq x' \land y \neq y'}} K(x_0 - x, y_0 - y)K(x_0 - x', y_0 - y')Cov((x, y), (x', y')) \quad (8.8) $$

where $\sigma^2(x_0, y_0)$ is the variance of the convolved signal and $Cov((x, y), (x', y'))$ is the covariance between two pixels at position $(x, y)$ and $(x', y')$. If the covariance is null or negligible, then the variance of the convolved signal is given by the variance map convolved by the square of the kernel. In theory, neglecting the covariance is not entirely appropriate since we know there must be some, though we do not know how much there is. Still, neglecting it is at least consistent with the previous assumption that the signal follows Poisson's statistics.

### 8.1.2.3 Pixel-per-pixel SED fitting with LEPHARE and Cigale

Once the pre-processing steps were done and the variance maps correctly estimated, the last part was to extract the flux in each pixel and in each band, produce tables according to the input format required by LEPHARE and CIGALE, define the models and their parameters for the SED fitting, run the fit, and finally recover the output tables, reconstruct the resolved maps, and analyse them. Discussing the variation of the results when using different SED fitting models is important and shall be done in a future analysis. Instead, in what follows, the goal is to see how much impact resolved stellar mass maps have on the galaxies' morphology and on the derived fraction of mass found in their DM component when propagating their best-fit morphological model to their kinematics modelling compared to the previous assumption of a constant M/L throughout the galaxies. Thus, this constant M/L used in the following parts will be derived with the same SED fitting models as those used to produce the resolved maps. For LEPHARE, we used the Bruzual & Charlot (2003) SSP library available as an additional package with a Calzetti et al. (2000) attenuation law. In LEPHARE, the galaxy spectra are automatically computed for a given SSP library with fixed metallicity values, IMF, and SFH. Thus, for the Bruzual & Charlot (2003) library it comes with a Chabrier (2003) IMF, three metallicity values, and an exponentially declining SFH. For CIGALE we also used a Bruzual & Charlot (2003) SSP library with a Chabrier (2003) IMF, a fixed metallicity of 0.02 dex, a delayed exponential SFH with a final burst-quenching phase (see Sect. 3.2.5.2), and a dust attenuation that is modelled with a single powerlaw for both young and old stars with an adjustable attenuation amplitude between both. Finally, we also included the nebular module to take into account the contribution of emission lines. Three examples of resolved stellar mass maps for galaxies ID1, ID13, and ID1002 are shown in Fig. 8.5. On the leftmost column is shown the map derived using LEPHARE, on the second column from the left the map derived with CIGALE, on the third column the stellar mass map obtained when multiplying the HST map in the F160W band by the M/L derived using CIGALE's total mass, and on the rightmost column the same as the F160W HST map but using the F435W band. The same mass scale has been used for each image and each galaxy. As already discussed in Sect. 8.1, we can see that the F160W band seems to trace much more robustly the stellar mass distribution than the F435W band. Overall, LEPHARE and CIGALE find consistent stellar mass maps between each other and their total mass derived by summing up the contribution of each unmasked pixel is nearly the same between the two codes. Furthermore, we checked whether the
<figure><img src="image_79.png" /><figcaption>Figure 8.5: Examples of resolved stellar mass maps (mass per pixel) produced with LEPHARE (leftmost column) and CIGALE (second column from the left) for galaxies ID1, ID13, and ID1002. Also represented are the stellar mass maps obtained with a constant M/L value when using the F160W (third column from the left) and F435W (rightmost column) band images. For these two latter maps the total stellar mass used to compute the M/L was that of CIGALE. The value of the total stellar mass integrated over all pixels is given in the maps derived with LEPHARE and CIGALE.</figcaption></figure>

total stellar mass derived in this manner was consistent with the stellar mass we would get if we
performed the SED fitting on the integrated SED rather than on each pixel individually and we
do find consistent results. On average CIGALE finds slightly larger pixel-per-pixel stellar mass
uncertainties (around 25%) than LEPHARE (around 15%) but this difference may not be
entirely significant given that LEPHARE adds by default lower systematic uncertainties on the
fluxes than CIGALE does<fnref n="7" /> before performing the SED fitting. These systematic uncertainties are
ubiquitous in SED fitting codes and are important to not underestimate the uncertainty on the
fluxes. In some way, they can be seen as the uncertainty on the models rather than on the input
data. The difference between the two codes is mostly found in the details. Looking at galaxy ID1,
we can see that LEPHARE finds much more structure than CIGALE or said otherwise that the
stellar mass map derived with CIGALE is smoother than that derived with LEPHARE. The reason

<footnote n="7">⁷Note however that these are the typical values used but we could have chosen the same value for both.</footnote>
behind that is not entirely clear yet and it is still being investigated at the moment. The fact that we use different SFHs might play a role. It would not be so much about the final burst-quenching phase that happens so late that it effectively only affects the SFR and not the stellar mass, but rather that the SFH used in CIGALE is delayed. What might be more likely though is that LEPHARE allows for three metallicity values whereas we fixed the metallicity in CIGALE to one value only. If the metallicity is different then it may also affect the age because of the age-metallicity degeneracy which in turn will affect the amount of stellar mass that is created. A more thorough analysis of the reliability of the stellar mass maps derived using different codes and models shall be carried out in the near future. For now, we will keep using these maps and we will investigate their effect on the morphological and kinematics parameters in Sect. 8.1.4. But before doing so, I will quickly discuss in the following section the current implementation of the code.

### 8.1.3 A python library to easily produce resolved maps from pixel-per-pixel SED fitting

Once we had converged on how to implement the framework to produce resolved maps of physical parameters, I implemented it into a python module<fnref n="8" />. This module has the advantage of producing resolved maps in a simple way by automatising as much as possible the entire pipeline while still retaining the flexibility that is required to deal with different SED fitting codes and potentially with different input data. Originally, the code only worked with LEPHARE but it has been updated since then to also work with CIGALE. No plan has been made yet to integrate other SED fitting codes but the modular design of the module should allow any one who is interested to do so in a relatively simple way. Thus I will focus on CIGALE and LEPHARE implementations below<fnref n="9" />. In essence, there are five key python objects that must be used to produce and handle resolved maps with this module: `Filter`, `FilterList`, `Catalogue`, `SED`, and `Output`. The `Filter` object corresponds to a single band and contains the data as well as the variance map and zero point to convert input data into fluxes or magnitudes. All the `Filter` objects must be assembled into a `FilterList` with a common mask, the code that is going to be used for the SED fitting, the redshift of the galaxy, and a few more optional keywords. The `FilterList` must then be converted into a `Catalogue` object, either `CigaleCat` for CIGALE or `LePhareCat` for LEPHARE. It essentially contains the data that has already been pre-formatted to be correctly handled by either SED fitting codes. Once the data have been prepared accordingly, the SED fitting can be run by first creating a `SED` object and then calling it. Each code has its own `SED` object: `CigaleSED` or `LePhareSED`. These two objects do not expect the same input parameters since the SED fitting codes do not work exactly in the same way. I will quickly discuss each separately below. Nevertheless, a `SED` object contains all the modules and physical parameters that are needed to perform the SED fitting. Once created, either can be run and it will launch the corresponding SED fitting code on the data contained in the `Catalogue`. At the end, they will both return an `Output` object (`CigaleOutput` or `LePhareOutput`) that will read the result and store it. If the SED fitting has already been done beforehand, it is also possible to directly load the result without having to redo it. For LEPHARE, the code will perform some transforms to

<footnote n="8">⁸https://github.com/WilfriedMercier/SED</footnote>

<footnote n="9">⁹This is not a SED fitting code since it relies on LEPHARE and CIGALE. Thus, it cannot work adequately if neither codes are properly set up first. Please check that the SED fitting codes work in stand-alone with the intended modules on a test case first before using this module. As an example, if one wants to use LEPHARE with the external Bruzual & Charlot (2003) SSP library that can be recovered on LEPHARE's website, then please check that the external library is properly installed before trying to use it.</footnote>
the data before creating the **Catalogue** so it is recommended to use the **SED.link** method to directly transfer the transforms so that it can correctly produce resolved maps with the right unit. It is also possible to provide them manually if necessary. For **CIGALE**, the only thing the **SED.link** method does is to provide the shape of the input/output data to produce resolved maps. Once again, the shape can also be given manually if required. More details about the specific input for each object and each method can be found in the online documentation<fnref n="10" />. An example of the code that has been used to produce the stellar mass map for galaxy ID1 represented in Fig. 8.5 is shown below. It is split in three parts as follows: (i) a common part that spans the creation of the **Filter** objects all the way through to the **Catalogue**, (ii) a part specific to LEPHARE, and (iii) another part specific to CIGALE. Given that LEPHARE pre-computes the models in multiple steps before running the SED fitting, it is possible to skip them to speed up the process if the correct models have already been generated (e.g. in a previous fit). This can be done so using the keywords `skipSEDgen`, `skipFiltergen`, `skipMagQSO`, `skipMagStar`, and `skipMagGal` when calling the `LePhareSED` object. Multiple properties, such as the libraries to use, the filters to consider, or the systematic error to add can be set using the `properties` keyword when generating the `LePhareSED` object. For CIGALE, it works slightly differently. Because CIGALE is built using a modular design, the same has been used for **CigaleSED**. The **SFH**, **SSP**, **nebular**, and **attenuation** modules must be provided when creating the **CigaleSED** object using the corresponding **SED.cigmod** objects. A list of at least one module must be given for each. When calling the **CigaleSED** object additional parameters can be given such as the bands to use for the fit or the number of cores for parallelisation. All the details with specific examples can be found in the documentation.

<footnote n="10">¹⁰https://wilfriedmercier.github.io/SED/</footnote>
from astropy.io import fits
import os.path as opath
import SED
# HST zero points
zeropoints = [25.68, 26.51, 25.69, 25.94, 24.87, 26.27, 26.23, 26.45, 25.94]
# HST bands (same order as zero points)
bands = ['435', '606', '775', '814', '850', '105', '125', '140', '160']
# Bands names in Cigale
band_names = [f'F{band}LP' if band == '850' else f'F{band}W' for band in bands]
# Band names in LePhare
band_names = [f'F{band}LP' if band == '850' else f'F{band}W' for band in bands]
dataFiles = [] # Flux maps
data2Files = [] # Flux maps convolved by the PSF squared
varFiles = [] # Variance maps
# Get input file names
for band in bands:
    file = opath abspath(f'1_{band}_smooth.fits')
    dataFiles.append(file)
    file2 = opath abspath(f'1_{band}_kernel Smooth.fits')
    data2Files.append(file2)
    vfile = opath abspath(f'1_{band}_variance sm.fits')
    varFiles.append(vfile)
# Get mask file
mfile = opath abspath(f'mask 1.fits')
with fits.open(mfile) as hdul:
    mask = hdul[0].data == 0 # True where data must be masked
# Generate a list of Filter objects
filt s = []
for band, data, data2, var, zpt in zip(band_names, dataFiles, data2Files, varFiles, zeropoints):
    filt s.append(SED.Filter(band, data, data2, var, zpt))
# Create the FilterList object
scaleFactor is a parameter important for LePhare only. Cigale does not use it.
cleanMethod handles pixels with negative values in different ways
flist = SED.FilterList(filt s, mask,
    code = SED.SEDcode.CIGALE,
    redshift = 0.622,
    cleanMethod = SED.CleanMethod. ZERO,
    scaleFactor = 1,
    texpFac = 4
)
# Convert to CigaleCat or LePhareCat. The conversion to the correct catalogue is
# done based on the code given in FilterList
catalogue = flist.toCatalogue(galName)

Code 8.1: First part of the example of a python script to perform SED fitting.
```python
1 # Custom properties
2 filt_used = [opath.join('hst_filters', i) for i in ['HST_ACS_WFC.F435W', 'HST_ACS_WFC.F606W',
3                                         'HST_ACS_WFC.F775W', 'HST_ACS_WFC.F814W',
4                                         'HST_ACS_WFC.F850LP', 'HST_WFC3_IR.F105W',
5                                         'HST_WFC3_IR.F125W', 'HST_WFC3_IR.F140W',
6                                         'HST_WFC3_IR.F160W']]
7
8 # Systematic error set to 0.03 in this case
9 # See LePhareSED docstring for other parameters
10 properties = {'FILTER_LIST': filt_used,
11              'ERR_SCALE': [0.03]*len(filt_used)}
12
13 # Create LePhare SED object
14 sedobj = sed.LePhareSED(f'{gal}', properties=properties)
15
16 # Output parameters to keep
17 oparams = ['K_COR()', 'MASS BEST', 'MASS INF', 'MASS MED', 'MASS SUP', 'SFR BEST', 'SFR INF',
18            'SFR MED', 'SFR SUP', 'SSFR BEST', 'SSFR INF', 'SSFR MED', 'SSFR SUP', 'AGE BEST',
19            'AGE INF', 'AGE MED', 'AGE SUP']
20
21 # Run SED fitting
22 # skipXXX parameters allow to skip the models generation if already done in a previous
23 # run to save time
24 # Comment if the SED fitting has already been done
25 output = sedobj(catalogue,
26                 outputParams = oparams,
27                 skipSEDgen = False,
28                 skipFilterGen = False,
29                 skipMagQSO = False,
30                 skipMagStar = False,
31                 skipMagGal = False
32             )
33
34 # Uncomment to load the result if SED fitting has already been done
35 #output = sed.LePhareOutput('1/1.out')
36
37 # Recommended with LePhare
38 output.link(filterList)
39
40 # Generate a resolved stellar mass map
41 mstar = output.toImage('mass_med')
```

Code 8.2: Second part of the example of a python script to perform SED fitting that is specific to LEPHARE.
```python
1 # Modules are given as lists because Cigale accepts multiple modules
2 # e.g. multiple SSPs or SFHs
3
4 SFH = [SED.cigmod.SFHDELAYEDBQmodule(tau_main = [250, 500, 1000, 2000, 4000, 6000, 8000],
5                                  age_main = [2500, 5000, 7500, 10000, 12500],
6                                  age_bq    = [10, 25, 50, 75, 100, 150, 200],
7                                  r_sfr     = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.25, 1.5,
8                                     1.75, 2.0, 5.0, 10.0],
9                                  sfr_A     = [1.0],
10                                 normalise = True
11 )]
12
13 SSP = [SED.cigmod.BCO3module(imf              = SED.IMF.CHABRIER,
14                             separation_age = [8],
15                             metallicity  = [0.02]
16 )]
17
18 nebular = [SED.cigmod.NEBULARmodule(logU         = [-2.0],
19                                     f Esc          = [0.0],
20                                     f Dust          = [0.0],
21                                     lines width   = [300.0],
22                                     include emission = True
23 )]
24
25 attenuation = [SED.cigmod.DUSTATT_POWERLAWmodule(Av_young = [0.0, 0.25, 0.5, .75, 1.0,
26                                                     1.25, 1.5, 1.75, 2.0, 2.25,
27                                                     2.5, 2.75, 3.0],
28                                                     Av_old_factor     = [0.44],
29                                                     uv_bump_wavelength = [217.5],
30                                                     uv_bump_width    = [35.0],
31                                                     uv_bump_amplitude  = [0.0, 1.5, 3.0],
32                                                     powerlaw_slope   = [-0.7],
33                                                     filters         = ' & '.join(band_names)
34 )]
35
36 # Create the Cigale SED object
37 sedobj = SED.CigaleSED(galName, band_names,
38                       uncertainties = [True]*len(band_names),
39                       SFH            = SFH,
40                       SSP            = SSP,
41                       nebular        = nebular,
42                       attenuation    = attenuation,
43                       flux_uncertainty = 0.1
44 )
45
46 # Run SED fitting
47 # Comment if the SED fitting has already been done
48 output = sedobj(catalogue,
49                  ncores          = 16,          # Number of threads to use
50                  physical_properties = None,   # None means all properties will be computed
51                  bands          = band_names,  # Estimate the flux for all the bands
52                  save_best_sed   = False,
53                  save_chi2       = False,
54                  lim_flag        = False,
55                  redshift_decimals = 2,
56                  blocks          = 1,
57              )
58
59 # Generate a resolved stellar mass map
60 mass_star = output.toImage('bayes.stellar.m_star', shape=(100, 100))
```

Code 8.3: Third part of the example of a python script to perform SED fitting that is specific to CIGALE.
8.1.4 A quick look at the impact of resolved maps on morphological
and kinematics models

<figure><img src="image_80.png" /><figcaption>Figure 8.6: Comparison of the best-fit GALFIT results for galaxy ID1 when using CIGALE’s stellar mass map (top row) and when using the F160W HST band (bottom row). For comparison, the images in the F160W band have been transformed into mass maps by multiplying by the M/L derived using CIGALE’s total mass. The disk size, the disk axis ratio, and the B/T computed at one effective radius are given for each fit and the ellipsis that encircles one effective radius is shown on each model as a blue line. The residuals (observations - model) are shown on the rightmost column.</figcaption></figure>

To conclude this section, a quick analysis of the impact of resolved maps on the morphological
and kinematics parameters derived with the formalisms discussed in Chapters 4 and 5 is
presented. A more statistically significant analysis shall be done first on the entire test sample
and then on a larger sample of galaxies in the MUSE-HUDF in the future. In the meantime,
below is shown a comparison of the morphology and the kinematics of galaxy ID1, that is the
largest galaxy in the test sample which has, based on its F160W map, a non-negligible bulge
contribution. Besides, this galaxy will also be used as benchmark for the machine learning
application in Sect. 8.2, thus it makes sense to perform its full dynamical modelling with this
framework. The best-fit GALFIT bulge-disk decompositions for galaxy ID1 using the resolved
stellar mass map derived with CIGALE (top row) and the F160W map (bottom row) are shown in
Fig. 8.6. The images in the F160W band have been converted into stellar masses using the M/L
derived using the total stellar mass from CIGALE. The rightmost column represents the residuals
in absolute units, hence observations minus model. As an indication, the ellipsis encircling half of
the disk’s total mass is shown with a blue line. As can be seen, the model from CIGALE finds a
smaller disk and a weaker bulge contribution than the M/L model at the galaxy’s half-light
radius. The disk ellipticity is also smaller, meaning that the galaxy is found to be more face-on,
which should have an impact on the kinematics model when de-projecting the circular velocity.
The reason why the model of the stellar mass map finds a smaller but more prominent disk is
because it fits a denser disk model in the inner parts that forces the bulge component to be
smaller to fit the remaining residuals near the centre. Thus, at one global effective radius the disk
already significantly dominates the flux budget. On the other hand, the central peak visible in the
F160W band is rather fitted by the bulge component, certainly because it is steeper so that the
disk alone is not sufficient to fit the inner parts. Thus, this leaves the disk with the possibility to
have a shallower core and rather fit the outer parts, hence increasing its radius. The residuals are
quite similar between the two models. The stellar mass in the arms is underestimated and that in
the inter-arms region is overestimated which is expected given that we do not try to model the
arms or the clumps but only a smooth disk. The spirals arms are more visible in the residuals of
the F160W band modelling but that is only because there is more structure visible in this band
than in the stellar mass map.

GALFIT was originally designed to fit photometric observations in unit of electrons or flux and not
in unit of mass. Nevertheless, by tweaking some header keywords in the input file, it is possible to
use it to fit stellar mass maps. However, there is one more aspect that must be handled correctly:
the variance map that is used to weight the value of each pixel during the fit. By default, GALFIT
will convert the input map in electrons and then it will assume that the pixels follow Poisson’s
statistics unless a variance map is provided by the user. However, it is unlikely that our maps
simply follow Poisson’s statistics and it is not yet clear whether the uncertainty on the mass in
each pixel that comes out of the SED fitting is sufficient for the morphological modelling given
that it is mostly dominated by the systematic uncertainty added by CIGALE before the fit is done.
In turn, this will generate correlations between the stellar mass and its uncertainty in each pixel
which will therefore produce relatively flat stellar mass S/N maps. Such maps are quite different
from what would be expected for Poisson’s statistics and it will therefore have an impact on the
fit. As an indication, Poisson’s statistics was assumed in Fig. 8.6. Instead, if we provide the
uncertainty from the SED fitting, we get the following morphological parameters: $R_{eff,d} = 1.29''$,
$B/T(R_{eff}) = 52\%$, and $q = 0.92$ which are closer to the values found using the F160W map,
except for the disk axis ratio that remains nearly unchanged. This example illustrates the
importance of the weight map on the best-fit parameters, which is one of the last aspects that we
are currently working on.

Finally, we can also check the impact of using the stellar mass map, compared to the F160W
HST band, on the best-fit velocity field model and on the derived DM fraction when propagating
the bulge and disk components to the kinematics modelling of the ionised gas component. This
comparison performed on galaxy ID1 is shown in Fig. 8.7. The two velocity fields and their
residuals are shown on the second and third columns with respect to the left-hand side. The