4 h 20 min, one is slightly above two and a half hours and four fields only have one hour. The location of the FoVs, as well as their rotation angle with respect to the north, and their on-source exposure time are summarised in Table. 3.1. Each OB was processed through the MUSE standard pipeline described in Weilbacher et al. (2020) that includes all the common steps to process raw data from 3D spectrographs such as but not limited to producing bad pixel tables, bias frames, dark frames, flat-fields, flux calibration, sky subtraction. Sky subtraction was applied on each exposure separately before aligning and then combining them using stars in the FoV and was further improved by applying the Zurich Atmosphere Purge (ZAP) software (Soto et al. 2016) on the final data cube. For more details about the data reduction process, see the survey paper (Epinat et al., in prep.).

### 3.2.2 Redshift determination

Contrary to what was done in the MUSE-HDFS and MUSE-HUDF surveys (Bacon et al. 2015, 2017), galaxies in the MAGIC survey were not blindly detected since multi-band photometric observations were required to derive key physical parameters such as their stellar mass or SFR. Instead, any object listed in the COSMOS2015 catalogue of Laigle et al. (2016) that fell in one of the MUSE FoVs was added to the MAGIC catalogue and had its redshift and physical properties derived. Recently, this step was further refined by using the newest COSMOS2020 catalogue of Weaver et al. (2022). Because this last step was done quite recently, these new galaxies do not appear in Abril-Melgarejo et al. (2021), Mercier et al. (2022), and in Chapter. 6 but are taken into account in the analysis of the angular momentum presented in Chapter 7. In both COSMOS2015 and COSMOS2020 catalogues the $3\sigma$ limiting magnitude measured in an aperture of $3''$ is $z_{\text{app}}^{++} < 25.9$ (25.7 in COSMOS2020, see Laigle et al. 2016; Weaver et al. 2022). This means that galaxies brighter than this limit are present in the COSMOS catalogue and are therefore also part of the MAGIC catalogue. Thus, the MAGIC survey is technically speaking almost magnitude-limited in the sense that the flux limit does not come from the MUSE observations themselves but from COSMOS. Nevertheless, as will be quickly discussed in Sect. 3.2.3, this is not particularly problematic as the MAGIC survey is nearly complete up to an upper apparent magnitude of roughly 24.5 in the $z^{++}$ band, that is within the detection limit of the COSMOS catalogue.

Beforehand, galaxies in the MAGIC survey only had photometric redshifts from COSMOS, except those observed with the VIMOS multi-object spectrograph as part of the zCOSMOS and COSMOS-WALL observations. However, MUSE is a 3D spectrograph which means that for any object in the catalogue a spectrum can be extracted and a (potentially tentative) redshift can be measured. Thus, for each object a PSF-weighted spectrum was extracted and a redshift was measured similarly to what was done in Bacon et al. (2015, 2017) and Inami et al. (2017) for the MUSE-HDFS and MUSE-HUDF surveys. To be more precise, this was achieved with the help of the redshift finding algorithm Manual and Automatic Redshifting Software (MARZ) described in Hinton et al. (2016) using absorption and emission lines and continuum features. Depending on the redshift of the galaxy, the main absorption lines are the Balmer lines, CaII H, CaII K, and the G-band from CH molecules. The main emission lines are also the Balmer lines (chiefly $H\alpha$ but also $H\beta$), [OII], [OIII], whereas continuum features are the Balmer break, the Lyman break, and D4000. We refer the reader to Fig. 2.5 for more details as to which lines and features are available depending on the redshift of the galaxies and whether AO is used or not. For each object a redshift was assigned with a given confidence flag (hereafter CONFID) ranging from zero to three:
<figure><img src="image_13.png" /><figcaption>Figure 3.2: Examples of MUSE spectra and their main spectral features extracted in a 3" aperture for six different types of galaxies. From top to bottom: (i), (ii), and (iii) 81\_CGR84, 50\_CGR84, and 129\_CGR35, three galaxies with similar redshift ($z \approx 0.7$) and stellar mass ($M_{*} \approx 10^{9} M_{\odot}$) and with a low, intermediate, and high spectroscopic redshift confidence flag (CONFID), respectively, (iv) 328\_CGR84, a galaxy similar to 129\_CGR35 but with a high stellar mass, (v) 139\_CGR28, a galaxy similar to 129\_CGR35 but found at higher redshift, and (vi) a galaxy similar to 129\_CGR35 but found at lower redshift. To the right of each spectrum are shown three MUSE narrow band images for different absorption and emission lines taken by collapsing the cube in a window of ten spectral pixels (spexels) around the line after removing the continuum and the corresponding HST image in the F814W band. All images have a dimension of 3" $\times$ 3". The gap visible between roughly 5800 Å and 6000 Å for the three first galaxies on the top is due to the sodium laser of the AO system. We can see how much more resolved become the spectral features when going from a low to a high CONFID value as well as how it becomes more and more difficult to detect strong emission lines when a galaxy get close to the MUSE desert that starts at $z \approx 1.5$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 3.3: Redshift distribution of the whole MAGIC survey from Epinat et al. (in prep.) including galaxies and stars (shown in the negative bin), for different confidence flags (CONFID): (1) tentative, (2) confident, and (3) highly confident. An inset zoomed-in on the [O II] emitters redshift range is also shown with smaller redshift bins to highlight the location of the structures. The bulk of the galaxies are located in the redshift range $0.2 \le z \le 1.5$ which corresponds to the redshift at which the [OII] doublet falls in the MUSE wavelength range ([OII] emitters). Beyond $z \approx 1.5$ there is a sheer drop in detection because of the MUSE redshift desert. The detections rise slightly again after $z \approx 2.9$ because Lyα enters the wavelength range. Overall, most [OII] emitters have confident MUSE spectroscopic redshifts ($CONFID \ge 2$).</figcaption></figure>

a value of zero means that no spectroscopic redshift was assigned, a value of one represents a tentative redshift that was chosen without any strong spectroscopic features visible, a value of two represents a confident redshift measurement that is based on a few faint lines or on a single strong feature (e.g. bright emission line), and a value of three corresponds to a highly confident redshift measurement that was based on multiple strong features (usually multiple bright absorption and emission lines). Examples of different spectra are shown in Fig. 3.2. I represent from top to bottom three galaxies with similar redshift $z \approx 0.7$ and stellar mass $M_{\star} \approx 10^9 M_{\odot}$ but with different CONFID values (81\_CGR84, 50\_CGR84, and 129\_CGR35), as well as galaxy 328\_CGR84 that is similar to 129\_CGR35 in terms of redshift and CONFID but with a much higher stellar mass around $10^{11} M_{\odot}$ and two more galaxies, 139\_CGR28 and 232\_CGR32, also similar to 129\_CGR35 in terms of stellar mass and CONFID with the first found at a higher redshift ($z \approx 1.4$) and the second at a lower redshift ($z \approx 0.3$). MUSE narrow-band images are also shown on the right-hand side of each spectrum centred on various absorption or emission lines and computed by collapsing the cube in a ten spexels window around each line after removing the continuum. As an indication, I also show the corresponding HST images of the galaxies in the F814W band. We can clearly see how spectral features become much more visible, and in particular how emission lines become brighter, as we go from CONFID = 1 to CONFID = 3. Additionally, we also see the wealth of emission lines available for relatively low redshifts and how bright emission and absorption lines fall out the MUSE range when we reach the beginning of the redshift desert at $z \approx 1.5$.
The redshift distribution of the entire MAGIC survey is shown in Fig. 3.3 with the highest CONFID shown in red and the lowest in orange, including galaxies but also stars represented in the negative redshift bin. The bulk of the galaxies are [OII] emitters, as expected given that all of the targeted groups fall within this redshift range and because it is one of the brightest emission lines available after Hα. The [OII] emitters represent a key population of emission line galaxies that was targeted by the MAGIC survey and which has been used in Abril-Melgarejo et al. (2021); Mercier et al. (2022) and in Chapters 6 and 7 to derive the galaxies’ kinematics and study the impact of the environment on their dynamical properties. Beyond $z \approx 1.5$ a sheer drop in detection is visible in Fig. 3.3 which is due to the galaxies entering the MUSE redshift desert (see Sect. 2.2.3). Beyond $z \approx 2.9$, the detection rises moderately because of the Lyα line entering the MUSE wavelength range. Nevertheless, the growth is not as significant as for the MUSE-HDFS and MUSE-HUDF surveys. This is because, as was previously stated, no blind detection was performed (except in a single field, CGR30) which significantly reduces the number of observed Lyα emitters. Besides, every field in the MAGIC survey is shallower than both the MUSE-HDFS and the MUSE-HUDF, thus further hindering the ability to detect such high-redshift objects. Overall, most of the derived MUSE spectroscopic redshifts are secure with CONFID values larger than or equal to two. For the entire redshift range, the objects with a CONFID equal to one only amount to roughly 16%. When restricting to [OII] emitters only, this number drops below 9%, showing the strength of MUSE in deriving precise spectroscopic redshifts in this redshift range.

### 3.2.3 Survey completeness

The goal of this short section is not to give a full account of the completeness of the full MAGIC survey but rather to highlight the values that we get for the main population of galaxies that were targeted, that is [OII] emitters. A more thorough description will be given in the survey paper (Epinat et al., in prep.). The completeness is shown in Fig. 3.4 for the entire redshift range as a function of the apparent $z^{++}$ magnitude ($z_{app}^{++}$). In our case, the completeness was defined as the number of galaxies with secure spectroscopic redshifts (CONFID $\ge$ 2) divided by the total number of galaxies with a photometric redshift, and this in bins of redshift and apparent magnitude. Thus, the completeness is a measure of our capacity to securely measure the galaxies’ redshifts. A value of unity means that 100% of the galaxies in the bin found in the photometric catalogue had their spectroscopic redshift precisely measured, whereas a null value means that no spectroscopic redshift could be correctly measured in the bin. We can see that the survey is nearly complete (at least above 80%) for the [OII] emitters redshift range up to an apparent upper magnitude of $z_{app}^{++} \approx 24.5$. Above this value, it quickly drops close to a few percent and to zero beyond the COSMOS limiting magnitude. The very high completeness for the [OII] emitters combined with the fact that the [OII] doublet is a bright emission line therefore make them the perfect galaxy population in the MAGIC survey to study the impact of the environment on their dynamical properties. Note that the completeness in the [OII] emitters range is different from the fraction of star-forming galaxies with the [OII] doublet detected. Indeed, as already stated in footnote 1 of this chapter, the [OII] emitters sample corresponds to galaxies in the redshift range $0.25 \le z \le 1.5$ where the [OII] doublet can be theoretically detected. Thus, among this sample are found star-forming galaxies that have the [OII] doublet in their spectrum as well as passive galaxies that do not have it but whose redshift can be determined from their absorption lines. As an indication, among the 1142 galaxies that are part of the [OII] emitters sample, 20% do not have any [OII] doublet detected at all and this number rises to 40% when including galaxies whose
[OII] doublet detection is below the $3\sigma$ level. As an indication, this corresponds to a flux detection limit of roughly$^2$ $2 \times 10^{-18}$ erg s$^{-1}$ cm$^{-2}$ which can be converted into a [OII]-derived SFR lower limit of roughly 0.1 M$_\odot$ yr$^{-1}$ at the median redshift of the structures in MAGIC (i.e. $z \approx 0.7$, see Sect. 3.2.5.2 for the details of the derivation of the SFR from the [OII] doublet flux).

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

Figure 3.4: Completeness of the MAGIC survey as a function of redshift and apparent ($z_{app}^{++}$) magnitude in the $z^{++}$ band adapted from Epinat et al. (in prep.). The orange dashed line represents the limiting magnitude of COSMOS and the black dotted vertical lines show, from left to right, the minimum redshift for [OII] emitters, the maximum redshift for [OII] emitters, and the minimum redshift for Ly$\alpha$ emitters to be detected in the MUSE wavelength range. Up to a $z_{app}^{++} \approx 24.5$ the completeness is very high for the [OII] emitters range.

Complementarily to what is shown in Fig. 3.3, the completeness drops to zero beyond $z = 1.5$ because of the [OII] doublet leaving the MUSE wavelength range and the galaxies entering the MUSE redshift desert. Past $z \approx 2.9$, it rises significantly again, though only for the intrinsically brightest sources, due to the Ly$\alpha$ line entering the wavelength range.

### 3.2.4 PSF and LSF measurements

An important aspect of dynamical studies such as those carried out as part of the MAGIC survey require a precise knowledge of both the PSF and the LSF. Indeed, in our MUSE data the LSF will affect each spexel by broadening the absorption and emission lines, whereas the PSF will blend in each spexel the signal coming from every other spexel in the cube. When extracting and modelling the stellar or gas kinematics of a galaxy out of a data cube, both the PSF and the LSF impact the velocity field and velocity dispersion map and it is referred to as beam smearing (see Sect. 5.3.2 for a discussion of its impact on kinematics maps and models).

<footnote>$^2$This value was estimated in Mercier et al. (2022) but there is a typo in the text with the minus sign missing. The correct value is the one given in this Thesis.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 3.5: Illustration of the fitting procedure to measure the wavelength dependence of the MUSE PSF FWHM for CGR35. Examples of GALFIT Gaussian fits for the two stars found in the field are shown on the top for four different wavelengths. Below is shown the raw (dotted lines) and the rolling-averaged wavelength variation of the PSF FWHM for the two stars, as well as the median curve used to derive the linear trend (thick black solid line). The FWHM values measured from the examples of the two stars are shown as circles with the same colour as the frames of the images. The raw lines are cut between roughly 5700 Å and 6000 Å since no data was available in this range because of the sodium notch filter (due to AO, shown in light yellow). The rolling-averaged curves were interpolated before deriving the median curve.</figcaption></figure>
### 3.2.4.1 LSF variation with wavelength

The MUSE LSF was the easiest of the two quantities to parametrise. To do so, we used the prescription of Bacon et al. (2015) and Guérou et al. (2017) who fitted with a Gaussian profile sky emission lines spread throughout the MUSE wavelength range in the MUSE-HDFS and UDF-10 (MUSE-HUDF) data cubes. As discussed in Sect. 3.1 of Guérou et al. (2017), for both data cubes the LSF FWHM variation with wavelength is well described by a similar second order polynomial given by

$$ FWHM_{LSF}(\lambda) = \lambda^2 \times 5.866 \times 10^{-8} - \lambda \times 9.187 \times 10^{-4} + 6.040, \quad (3.1) $$

where FWHM$_{LSF}$ and λ are both in angstrom, with λ the observed wavelength. Equation 3.1 was found to be very robust when measured from one spaxel to another, with variations of the order of 0.05 Å, and for two different data cubes obtained from observations carried at different observing periods. Thus, we did not try to model once again the LSF in MAGIC and used Eq. 3.1 throughout.

### 3.2.4.2 PSF modelling

The modelling of the PSF is a bit trickier because it is expected to change from one cube to another and must therefore be modelled independently for each MUSE field. Furthermore, it is also expected to vary with wavelength and might as well change depending on where the galaxies are located in the FoV. However, the lack of bright stars at multiple enough locations in the cubes made it difficult, if not impossible at all, to precisely measure any spatial variation of the PSF. Therefore, only the wavelength dependence was modelled and it was assumed to be constant at a given wavelength throughout the MUSE FoV.

Initially, the PSF was measured in Abril-Melgarejo et al. (2021) for all the MUSE fields except CGR35, CGR87, and CGR172 since the observations for these three fields were not over yet at the time. Because the only galaxies studied in Abril-Melgarejo et al. (2021) were all located in groups and because earlier results from the MUSE-HDFS (Bacon et al. 2015) and MUSE-HUDF (Bacon et al. 2017) surveys showed that the PSF FWHM dependence with wavelength can be approximated by a linearly decreasing function, the way it was first measured was to (i) find all the bright, but non-saturated, stars located in the fields, (ii) select every group studied in the analysis, (iii) for each group extract, for each star in the FoV, narrow-band images around the observed wavelengths of [OII] and [OIII], and then (iv) fit each narrow-band with a Gaussian profile to retrieve the PSF FWHM at two different wavelengths. In Mercier et al. (2022), the value of the MUSE PSF FWHM at the observed wavelength of [OII] (being the emission line used for the kinematics) was computed for each galaxy by interpolating in-between or extrapolating beyond the two values measured in Abril-Melgarejo et al. (2021), assuming a linearly decreasing relation. In addition, CGR35, CGR87, and CGR172 were also incorporated into the analysis and therefore required their MUSE PSF to be measured as well. A quick description of how it was done is given in Sect. 2.1 of the paper and in Chapter 6, but we describe it in slightly more details below.

Since the approach was different from Abril-Melgarejo et al. (2021) in the sense that we were now looking at galaxies in structures and in the field throughout the whole redshift range of [OII] emitters ($0.2 \le z \le 1.5$), I decided to use a more robust method to derive the FWHM dependence with wavelength for these fields. After identifying bright but non-saturated stars, as well as a low
surface brightness unresolved object located at $z \approx 2.67$ for CGR87, 100 sub-data cubes with dimensions $10 \times 10$ pixels<fnref n="3" /> were extracted for each star (and for the unresolved object). The sub-data cubes were then collapsed into narrow-band images with a fixed redshift slice of $\Delta z = 0.01$ that converts to a wavelength slice of $\Delta\lambda = \Delta z \times \lambda$. Each narrow-band was fitted independently of the others using GALFIT with symmetric (i) Gaussian and (ii) Moffat profiles with free $\beta$ parameter. Since results did not differ significantly between both models and to remain consistent with what was done in Abril-Melgarejo et al. (2021) the Gaussian model was kept. For each MUSE field, the wavelength dependence of the PSF FWHM was measured by (i) applying a rolling-average to the FWHM-wavelength curve of each star to remove small scale variations, (ii) median combining the averaged curves of all the stars in each field, and then (iii) fitting each median line with a linearly decreasing function. This procedure is illustrated for CGR35 in Fig. 3.5. Examples of GALFIT fits performed on two stars for four different wavelength values are shown on the top and the PSF wavelength dependence is shown below. In particular, the initial raw dependence is shown with dotted lines, the 5700 Å to 6000 Å range missing because of the sodium notch filter due to AO observations, the rolling-averaged curves as thin solid lines, and the median relation as the thick black solid line.

Even if the linear fit does not capture the full complexity of the wavelength dependence of the MUSE PSF FWHM, the clear advantage of using this method over the one used in Abril-Melgarejo et al. (2021) is that the relation that is derived does not depend on potential fluctuations that would happen at the [OII] and [OIII] wavelengths, therefore rendering it much more robust. Thus, as part of the MAGIC survey paper (Epinat et al., in prep.), it was decided to update the MUSE PSF values using the method described above at the exceptions of (i) extracting narrow-band images in an aperture of $6'' \times 6''$, (ii) collapsing them in a fixed wavelength interval of 100 Å, and (iii) modelling the stars with a flat background plus a Gaussian/Moffat profile, the latter being retained in the end. Thus, at present all MUSE fields have had their PSF measured in a consistent way and these values are the ones used in the analysis of the angular momentum in Chapter 7. The parameters for the parametrisation of the PSF are given in columns (6), (7), and (8) of Table 3.1. Column (6) represents the FWHM at a wavelength of 7000 Å (FWHM$_{7000}$) and column (7) represents the linear gradient of the FWHM with wavelength ($\partial$FWHM$_{PSF}$/$\partial\lambda$). Thus, the parametrisation writes

$$ \text{FWHM}_{\text{PSF}}(\lambda) = \frac{\partial \text{FWHM}_{\text{PSF}}}{\partial \lambda} \times (\lambda - 0.7) + \text{FWHM}_{7000}, \quad (3.2) $$

where $\lambda$ is the observed wavelength in µm. Column (8) provide the $\beta$ parameter of the Moffat profile. For the three latest MUSE fields CGR35, CGR87, and CGR172, the difference between this new modelling and the previous one is negligible since they are nearly identical. For the other fields, it leads to differences in FWHM of at most 0.3" at the blue end of the spectrum and 0.2" at the red end. Furthermore, because the PSF FWHM was updated after the analysis carried out in Mercier et al. (2022), we checked whether this would affect the best-fit kinematics parameters, in particular the circular velocity (see Sect. 5.2.1 and Chapter 6). Apart from a few galaxies that had poorly constrained kinematics parameters in the first place, no significant impact could be found.

<footnote n="3">³MUSE pixels, 0.2" across</footnote>
### 3.2.5 Derivation of galaxies’ physical parameters

Galaxies’ key physical parameters such as their stellar mass or their SFR were derived from SED fitting using the photometric bands available in the COSMOS catalogue of Laigle et al. (2016) in combination with their MUSE spectroscopic redshift. SED fitting is a complicated modelling process that relies on many assumptions and libraries. The point of this section is to briefly describe how the physical parameters were derived before discussing the physical properties of field galaxies and those in structures. More details can be found in Sect. 8.1.1 where the models, their parameters, and their underlying assumptions used in SED fitting, especially in the case of the Code Investigating GALaxy Emission (CIGALE), will be discussed.

#### 3.2.5.1 Early derivations with FAST

In Epinat et al. (2018), Boselli et al. (2019), and Abril-Melgarejo et al. (2021), these parameters were derived using the SED fitting code FAST (Kriek et al. 2009). In total, 29 out of the 32 bands of the COSMOS2015 catalogue were used and each galaxy’s SED was modelled using a synthetic library generated from the Single Stellar population (SSP) models of Conroy & Gunn (2010) with a Chabrier (2003) initial mass function (IMF), a Calzetti et al. (2000) extinction law, and an exponentially declining star formation history (SFH). These physical parameters were also used in Mercier et al. (2022) except for the SFR. Indeed, when compared with that given in the COSMOS catalogue, computed with the Photometric Analysis for Redshift Estimate (LEPHARE)<fnref n="4" /> code, the SFR values had an average difference of about 0.8 dex, which was three to four times larger than the discrepancy observed in stellar mass. Besides, the MS relation obtained with FAST is quite discretised and galaxies seem to be capped to an upper limiting line. Thus, because stellar masses seemed consistent between the two SED modellings but not the SFRs, because of the sparsity of the $M_\star$–SFR grid, and because of the presence of this limiting line, we decided to not rely on the SFR values from FAST.

#### 3.2.5.2 Star formation rates from the [OII] doublet

In Mercier et al. (2022), the SFRs were computed from the [OII] doublet luminosity using the Kennicutt (1992, 1998) relation after correcting for Galactic and intrinsic dust extinctions:

$$ \text{SFR} [\text{M}_\odot \text{yr}^{-1}] = (1.4 \pm 0.4) \times 10^{-41} L_{[\text{OII}]} [\text{erg s}^{-1}], \quad (3.3) $$

where $L_{[\text{OII}]}$ is the [OII] luminosity corrected for extinction. This relation relies on a Salpeter (1955) IMF but it is possible to convert it to other IMFs. For instance, to recover the SFR with a Kroupa (2002) IMF Eq. 3.3 must be multiplied by a factor of 1.5 (Gilbank et al. 2010) and in the case of a Chabrier (2003) IMF Eq. 3.3 must be multiplied by a factor of 0.62 (Boogaard et al. 2018). Ideally, we would have used $H\alpha$ to compute the SFR as it is more robust (Kennicutt 1998) but the line is not available for the entire redshift range for the [OII] emitters. Besides, modulo a larger scatter, the [OII] doublet is also a good tracer of star formation. The luminosity is linked to the flux of the galaxy through the usual equation

$$ L_{[\text{OII}]} = 4\pi D_L^2 F_{[\text{OII}],\text{corr}}, \quad (3.4) $$

<footnote n="4">⁴https://www.cfht.hawaii.edu/~arnouts/LEPHARE/lephare.html</footnote>
where $D_L$ is the luminosity distance (see Sect. 1.1.4) and $F_{[O II],corr}$ is the flux corrected of extinction. Following Kennicutt (1998), the correction of the observed flux $F_{[O II]}$<fnref n="5" /> must be computed by taking into account the Galactic extinction $A_{[O II],MW}$ at the observed [OII] wavelength and the intrinsic extinction of the galaxy computed at the rest-frame Hα wavelength ($A_{Hα}$) using the formula

$$F_{[O II],corr} = F_{[O II]} \times 10^{0.4(A_{H\alpha} + A_{[O II],MW})}. \qquad (3.5)$$

In Eq. 3.5, we must use the intrinsic extinction at the rest-frame wavelength of Hα rather than at that of [OII], as might be expected, because of the way the [OII] fluxes were calibrated in Kennicutt (1998). For both extinctions, one needs a model to go from an observed or derived quantity (e.g. from SED fitting) such as the extinction in the V band to the value at the required wavelength. We decided to use a Cardelli et al. (1989) law for the Galactic extinction. We do not write explicitly its expression because it is a complex polynomial function by segments whose coefficients vary depending on which part of the spectrum it is evaluated. It provides a parametrisation of the ratio between the Galactic extinction $A(\lambda)$ at a given wavelength and that in the V band ($A_V$) as a function of wavelength. Thus, one needs to know $A_V$ to derive the Galactic extinction at the wavelength of the observed [OII] doublet. This was done by using a total-to-selective extinction ratio $R_V = 3.1$, as is usual for the Milky Way, and the colour excess between the B and V bands $E(B-V)$ provided by FAST. By definition, the V band extinction and the colour excess are linked through the following relation:

$$A_V = R_V \times E(B-V), \qquad (3.6)$$

where $E(B-V)$ is defined as

$$ \begin{aligned} E(B-V) &= A_B - A_V \\ &= (m_B - M_B) - (m_V - M_V) \\ &= (m_B - m_V) - (M_B - M_V), \end{aligned} \qquad (3.7) $$

with $A_B$ the extinction in the B band and where the symbol $m$ represents an apparent magnitude and the symbol $M$ an absolute magnitude corrected for extinction. The k-correction does not appear in Eq. 3.7 because the dust in the Milky Way affects the redshifted light and not the intrinsic one. Since the extinction in the B band is higher in our Galaxy than in the V band, the flux in the B band will be more severely affected and the apparent magnitude in this band will therefore increase more significantly producing an excess of colour, hence the name. We could have done a similar derivation for the intrinsic extinction using the $A_V$ value provided by FAST and a model for the extinction law such as a Calzetti et al. (2000) curve but, as for the SFR, the large uncertainties on this parameter dissuaded us from using it. Instead, we relied on a prescription given by Gilbank et al. (2010, 2011) that empirically parametrises the Hα extinction as a function of stellar mass for the [OII] doublet when used to derive the SFR with the Kennicutt law. This parametrisation writes

$$A_{H\alpha} = 51.201 - 11.199 \log_{10} \left( \frac{M_{\star}}{M_{\odot}} \right) + 0.615 \log_{10} \left( \frac{M_{\star}}{M_{\odot}} \right), \qquad (3.8)$$

<footnote n="5">⁵The total flux is measured from the best-fit double Gaussian fit performed on the [OII] doublet. See Sect. 5.3.1.2 for more details.</footnote>
for stellar masses above $10^9$ M$_\odot$ and a constant value below.

### 3.2.5.3 New derivation with Cigale

For the survey paper (Epinat et al., in prep.) and following the release of the new COSMOS2020 version of the catalogue (Weaver et al. 2022), the physical parameters were derived again with the latest photometric measurements using CIGALE$^6$. A description of the main parameters and their range of values is shown in Table 3.2. Briefly, for this modelling the SSPs that were used are those of Bruzual & Charlot (2003) using a Salpeter (1955) IMF with a fixed metallicity of 0.02 dex. The SFH that was used is called a truncated delayed exponential law (see Ciesla et al. 2018, 2021) and it corresponds to a delayed exponential law (i.e. SFH ~ $t \times \exp\{-t/\tau\}$, with $\tau$ the e-folding time of the SFH) followed by a truncation episode of constant star formation that can be lower (quenching) or larger (burst) than the value of the SFH just before the episode. This phenomenological model allows to reduce biases when measuring the stellar mass and the SFR of quenched and starburst galaxies from SED modelling (Ciesla et al. 2018, 2021). Thus, the truncation episode was forced to happen at the very end of the SFH, with an age below 100 Myr (see Table. 3.2 for the full list of values). Finally, we used a modified starburst attenuation law (Boquien et al. 2019) which, for the parameters that were used, reduces to a Charlot & Fall (2000) law with a total-to-selective extinction ratio $R_V = 3.1$ and a colour excess varying from 0.001 to 0.7.

Two important parameters that come out of SED fitting are the galaxies’ stellar mass and their SFR. In CIGALE, the stellar mass is not a free parameter as in for instance LEPHARE but is derived by scaling the SFH to the data (for more details, see Boquien et al. 2019). Similarly, the SFR is not fitted directly but is a product of the modelling. What is truly fitted are the free parameters of the SFH that will determine the value of the SFR. By default, CIGALE provides three different values: (i) the instantaneous SFR (noted SFR$_{inst}$) which corresponds to the SFH evaluated during the last 1 Myr time-step, (ii) an average value computed over a short period of 10 Myr, and (iii) an average value computed over a longer period of 100 Myr. A comparison between the three different values is shown in Fig. 3.6. It is clear that the instantaneous SFR and the average value over 10 Myr are similar for most of the galaxies. The main reason is certainly that below 10 Myr the truncation episode that leads either to a burst or to a quenching of star formation has already happened for the majority of the galaxies. But, the SFH becomes constant afterwards and therefore there is no further changes in SFR. Only for galaxies that are modelled with a late truncation episode does this change, though the effect remains moderate. On the other hand, the discrepancy is much larger when using the average value over 100 Myr. Indeed, this corresponds to the time when the truncation is allowed to take place for the first time. Thus, this value will be sensitive to both the smooth part of the SFH and to the truncation. This is visible in Fig. 3.6 where the 100 Myr average value is larger than the instantaneous one for most of the galaxies given that CIGALE models most of them with a quenching phase. Star formation rates can be derived in various ways (emission line, ultra-violet or blue bands, emission in the infrared, etc.) but none of these methods provides an SFR estimate over a period as short as 1 Myr. The 10 Myr average value can be seen as a good compromise since it roughly corresponds to the timescale over which the SFR is measured when using emission lines. However, as is shown in Fig. 3.6, it is quite sensitive to the late truncation episode so one has somehow to trust that such an episode occurred in order to use it. Lastly, the 100 Myr average value is also a SFR proxy candidate which is much

<footnote>$^6$https://cigale.lam.fr/</footnote>
<figure><img src="image_17.png" /><figcaption>Figure 3.6: Comparison between the instantaneous SFR derived by CIGALE during the last 1 Myr time-step and the average values derived during the last 10 and 100 Myr. We only show [O II] emitters with CONFID ≥ 2. For most galaxies the instantaneous value and that averaged over 10 Myr are similar. However, there is a large discrepancy when using the 100 Myr averaged value. This indicates that for most galaxies the truncation episode happens between the last 100 Myr and 10 Myr of the SFH and is usually modelled as a quenching phase.</figcaption></figure>
less affected by the late truncation phase and more by the smooth part of the SFH. Thus, it can also be seen as a good compromise at the cost of having higher SFR values on average. Therefore, in what follows, we will use the two average values and we will drop the instantaneous SFR.

<table><caption>Table 3.2: Models and grid of parameters for the SED modelling with CIGALE of galaxies in MAGIC</caption><thead><tr><th>Parameter</th><th>Values</th><th>Description</th></tr><tr><th>(1)</th><th>(2)</th><th>(3)</th></tr></thead><tbody><tr><td colspan="3"><strong>Truncated delayed SFH</strong></td></tr><tr><td>Age (Gyr)</td><td>1e-03, 2e-03, 8e-03, 2.3e-02, 6.7e-02, 0.2e-01, 0.6, 1.6, 4.5, 12.9</td><td>Age of the oldest stars that contribute to the SED</td></tr><tr><td>τ (Gyr)</td><td>1e-03, 3e-03, 9e-03, 2.7e-02, 8.1e-02, 0.24, 0.73, 2.2, 6.6, 20</td><td>e-folding time of the exponential part of the SFH</td></tr><tr><td>Age trunc. (Myr)</td><td>1, 12, 23, 34, 45, 56, 67, 78, 89, 100</td><td>Age of the truncation episode</td></tr><tr><td>SFR ratio</td><td>1e-04, 2.2e-03, 4.6e-02, 1, 2.2e+01, 4.6e+02, 1e+04</td><td>Ratio of SFR after and before the truncation episode</td></tr><tr><td colspan="3"><strong>Bruzual et Charlot (2003) SSP</strong></td></tr><tr><td>IMF</td><td>Salpeter</td><td>Initial Mass Function used in the modelling</td></tr><tr><td>Metallicity</td><td>0.02</td><td>Metallicity values</td></tr><tr><td colspan="3"><strong>Modified starburst attenuation law</strong></td></tr><tr><td>E(B-V)</td><td>1e-3, 8e-2, 0.16, 0.23, 0.32, 0.39, 0.47, 0.55, 0.62, 0.7</td><td>Colour excess of the nebular lines and of the stellar continuum</td></tr><tr><td>R<sub>V</sub></td><td>3.1</td><td>Total to selective extinction ratio.</td></tr></tbody></table>

**Notes:** (1) Parameter name, (2) range of values used in the modelling, and (3) short description. Only the SFH, SSP, and attenuation law models and their parameters are shown in this table, the other models (nebular and dust template) having default values. The dust template model used is that of Dale et al. (2014).

When stellar mass is plotted against the SFR it is referred to as the MS. In Fig. 3.7 we show the MS for the entire MAGIC sample ($CONFID \ge 2$), that is including galaxies beyond the [OII] emitters range. Field galaxies and those in structures (see Sect. 3.3) are shown separately as squares and circles with a black contour, respectively. Additionally, I used the MS/red sequence separation line from Abril-Melgarejo et al. (2021) to highlight the difference between the population of star-forming galaxies that lie along the MS (shown in blue) and that of passive/quenched galaxies (shown in red - it also contains a few green valley galaxies). Before applying the separation, I normalised the SFRs to a redshift $z_0 = 0.7$ as in Mercier et al. (2022) to remove the intrinsic redshift dependence of the relation. This normalisation writes
$$ \log_{10} \text{SFR}_z = \log_{10} \text{SFR} - 2.8 \log_{10} \left( \frac{1+z}{1+z_0} \right), \qquad (3.9) $$

where SFR$_z$ is the normalised SFR and where the factor of 2.8 is derived from Speagle et al. (2014) and is more adapted for our stellar mass range than the value of 1.4 used in Abril-Melgarejo et al. (2021). However, the SFR values shown in the plot are those directly coming out of CIGALE. Above the main plot are also shown histograms for four different physical parameters. From left to right we have (i) the main stellar population age, (ii) the ratio of the SFR before and after the episode of burst/quenching of star formation, (iii) the SFR measured in the last CIGALE time step of 10 Myr, and (iv) the total stellar mass of the galaxies. Item number (ii), that is the SFR ratio, measures by how much the smooth exponentially declining SFH is changed when the burst/quenching episode happens. If it is below unity, it means that a quenching phase was modelled by CIGALE at the end of the SFH and, if it is above unity, it rather means that a burst phase was modelled. However, the significance of whether the galaxy experienced or not a burst/quenching phase is not as simple as just looking at the SFR ratio given that it also depends on the value of the SFH just before the episode happened.

A few points can be noted. First, the environment within which galaxies reside (field versus structures) is different between the MS and red sequence populations in the sense that passive galaxies are mostly found in structures whereas star-forming galaxies are an equal mix of galaxies in the field and in structures. Secondly, the red sequence population does not have the same physical parameter distribution than the MS population with on average: (i) older galaxies (peaking at an age of about 2 – 3 Gyr versus 300 Myr for MS galaxies), (ii) lower SFR values (by definition), and (iii) it is mostly populated by massive galaxies (the majority above $10^{10} M_\odot$). In addition, galaxies on the MS exhibit a much wider distribution of ages, dominating the young population and the low SFR ratios with the bulk of the objects modelled with a quenching episode in their SFH (negative values, see histogram in Fig. 3.7). On the other hand, galaxies on the red sequence are mostly modelled with old stellar populations and with a SFR ratio that peaks near unity. These differences arise because passive galaxies are modelled with old SFHs that smoothly increase with time, reaching their peak value and then smoothly decreasing to low SFRs. A non-negligible fraction (of the order of 10%) is also modelled with a quenching phase in their last billion years. On the other hand, star-forming galaxies have a much wider range of SFHs which, depending on their position on the MS, can be modelled with or without a burst/quenching phase. Among the objects that do have this phase, some experience a high exponential growth of their SFH followed by a quick quenching phase, while others experience a weak exponential growth of their SFH, then followed by a very strong burst episode. In both cases, these types of SFH certainly allow CIGALE to produce blue SEDs since the old stellar population required to produce the red part of the spectrum would not have had time to dominate the emission yet.

### 3.2.5.4 Comparison of physical parameters between the different modellings

In Abril-Melgarejo et al. (2021), both the stellar mass and the SFR derived from FAST were used to separate galaxies between star-forming and passive objects. In Mercier et al. (2022), we kept using the stellar mass from FAST but computed a new estimate of the SFR using the [OII] doublet and in the survey paper (Epinat et al., in prep.), as well as in the analysis of the angular momentum performed in Chapter 7, we derived again new values of the stellar mass and the SFR using CIGALE. Among the two parameters, the stellar mass is probably the most important one
<figure><img src="image_18.png" /><figcaption>Figure 3.7: Main sequence relation obtained from the SFR and the stellar mass values derived with CIGALE for the full MAGIC sample ($CONFID \ge 2$) and histograms for, from left to right, (i) the age of the galaxy, (ii) the ratio of the SFR before and after the burst/quenching episode, (iii) the instantaneous SFR, and (iv) the galaxies’ stellar mass. As an indication, I show field galaxies with square symbols and galaxies in structures (see Sect. 3.3) with circles with a black contour. Furthermore, galaxies are separated between star-forming (blue symbols in the MS and a blue solid line in the histograms) and passive systems (red symbols in the MS and a red dashed line in the histograms). To split the sample between star-forming and passive galaxies I used the separation given in Abril-Melgarejo et al. (2021) after normalising the SFRs to $z_0 = 0.7$ using the same expression as in Mercier et al. (2022). This normalisation corrects for the redshift evolution of the MS and is well suited for the stellar mass range that is probed. For illustration purposes, the raw values of the SFR and stellar mass are shown.</figcaption></figure>
<figure><img src="image_19.png" /></figure>

Figure 3.8: Comparison between the stellar mass derived from CIGALE and that from FAST. Only the [OII] emitters with CONFID ≥ 2 are shown. A positive difference means that CIGALE overestimates the stellar mass with respect to FAST. Overall, values are consistent within plus or minus 0.1 dex and with a small median bias of 0.05 dex (median value represented as a plain line and the 25th and 75th quantiles as dashed lines). Only low-mass objects according to CIGALE were modelled as more massive galaxies by FAST which is certainly due to different SFHs. Objects with a difference below −0.5 dex have been identified with downward pointing red triangles.

for the analysis since it explicitly appears in the mass models (see Sect. 5.2.4 for examples) and in
the various scaling relations that have been investigated (see Chapters 6 and 7). On the contrary,
the SFR does not impact the dynamical models. Thus, as long as the stellar mass is well
constrained, a large uncertainty on the SFR should not be too problematic to derive the galaxies'
dynamics<sup>7</sup>. Only when the SFR appears explicitly in a scaling relation it becomes important to
constrain it efficiently. This was the case in Mercier et al. (2022) and in Chapter 6 where the MS
relation, the size-mass relation, and the TFR were investigated. Hence, I show in what follows a
comparison of the stellar mass estimated with CIGALE and with FAST and I discuss how the MS
relation is impacted when using different stellar mass and SFR estimates.

The difference in stellar mass between CIGALE and FAST as a function of that of CIGALE is shown in Fig. 3.8 (a positive difference means that the value from CIGALE is larger), focussing on the [OII] emitters with secure spectroscopic redshifts (CONFID ≥ 2) only. For the bulk of the galaxies, CIGALE and FAST find consistent stellar masses within at most 0.5 dex, except towards the lowest CIGALE-derived stellar masses where FAST derives larger values (downward pointing red triangles in the figure). The median difference between the two codes is only of 0.05 dex when excluding outliers at low CIGALE stellar masses, whereas the 25th and 75th quantiles are found at −0.09 dex and 0.2 dex, meaning that most galaxies are found within plus or minus 0.15 dex from the median value. This shows that on average CIGALE tends to produce slightly larger (by roughly 25%) stellar masses than FAST. These differences can certainly be explained by different SFH and IMF models being used, in particular given that a Salpeter (1955) IMF should produce on average higher stellar masses than a Chabrier (2003) IMF. Furthermore, such a discrepancy between two SED fitting codes is consistent with the literature (e.g. Walcher et al. 2011).

Concerning the larger discrepancy observed at low stellar masses, we can compare it with Fig. 3.9
that shows the MS relation using the different stellar mass and SFR estimates at our disposal. In
particular, we can focus on the top left plot that shows the MS for CIGALE using the 10 Myr

<footnote n="7"><sup>7</sup>Though, a sample selection based on the SFR can have an indirect impact on the analysis.</footnote>
<figure><img src="image_20.png" /></figure>

Figure 3.9: Comparison of four different MS relations that we can build when using the various stellar mass and SFR estimates presented in Sect. 3.2.5. Only the [O II] emitters with secure spectroscopic redshifts are shown (CONFID ≥ 2). On the left is shown the MS relation using the stellar mass and SFR from CIGALE (top: average over 10 Myr, bottom: average over 100 Myr). The top right plot shows the MS using the values from FAST and the bottom right plot shows the relation when using the stellar mass from FAST and the SFR derived from the [O II] doublet. The same thick red line was added to each plot to guide the eye. It corresponds to the best-fit line of CIGALE’s 100 Myr relation obtained by fitting the MS after removing the red sequence while taking into account the uncertainties on both the stellar mass and the SFR during the fit. All SFRs and stellar masses were estimated in the same aperture of 3'' on the plane of the sky. Galaxies with a large stellar mass difference identified in Fig. 3.8 are also shown with the same downward pointing red triangle symbol.
average SFR estimate, as well as the top right plot that shows the MS using the values from FAST. Similarly to Fig. 3.8, the galaxies with the largest discrepancy ($\Delta \log_{10} M_* < -0.5$) are identified with downward pointing red triangles. With CIGALE, these objects are above the MS whereas with FAST they are located throughout. Thus, the difference observed in stellar mass for these objects is most likely the result of using two different parametrisations of the SFH, FAST lacking a burst/quenching phase contrary to CIGALE.

As previously mentioned, Fig. 3.9 represents four different MS relations that can be built using the different stellar mass and SFR estimates. The top right plot shows the relation derived with FAST that was used in Abril-Melgarejo et al. (2021). The bottom right plot represents the relation that uses the SFR estimated from the [OII] doublet and the stellar mass from FAST that was used in Mercier et al. (2022, see also Chapter 6). On the left are shown the MS relations using the parameters from CIGALE (the top plot uses the 10 Myr average SFR and the bottom the 100 Myr average values). These MS relations have not been normalised to correct for their redshift evolution. In each of these relations we recover the shape of the typical MS with low-mass galaxies that tend to have lower SFR values than high-mass galaxies, but both the slope, zero point, and the scatter of the MS significantly change from one relation to another, as discussed in more details below. Focussing first on the two MS relations using the parameters from CIGALE, we see a positive vertical offset of the zero point of about 1 dex when using the 100 Myr average SFR values with respect to the 10 Myr average, consistent with what is visible on the right-hand side of Fig. 3.6. In both cases the red sequence is clearly apparent below the MS, with galaxies spread along 3-4 dex of SFR values. However, the scatter along the vertical direction in the MS is larger by at least 0.3 dex when using the 10 Myr average SFR, certainly because of the burst/quenching phase that spreads galaxies above and below the relation that was built by the smooth part of the SFH. We see on the top right plot the discretisation of the $M_* - \text{SFR}$ grid that produces strong correlation lines. This affects particularly the top of the MS where no galaxies are found above the upper limiting line discussed in Sect. 3.2.5.1, as well as the red sequence. Furthermore, FAST's SFR values were derived using a Chabrier (2003) IMF, whereas we used a Salpeter (1955) IMF for CIGALE. This should affect the shape of the MS and in particular the estimate of the SFR. According to Boogaard et al. (2018), if the number of massive stars (> 10$M_\odot$) is assumed to remain the same between the two IMFs, then converting from a Chabrier (2003) to a Salpeter (1955) IMF can be done by diving the SFR by a factor of 0.62. Hence, the MS from FAST can be compared in an IMF-independent way to the other MS relations in Fig. 3.9 by vertically offsetting it by roughly 0.2 dex. Note that this value is quite consistent with the 0.3 dex zero point difference observed between this MS and that from CIGALE when using the 100 Myr average SFR values. Finally, the relation in the bottom right plot is probably the most different of the four. As for FAST, this MS has a zero point lower by roughly 0.3 dex with respect to CIGALE's 100 Myr average SFR relation, though this difference cannot be accounted for by IMF differences since it uses the same Salpeter (1955) IMF as CIGALE. Furthermore, this relation also has a lower slope of about 0.7 when the other relations rather find a slope of 0.8 and it tends to have a SFR scatter lower by 0.1 dex on average with respect to CIGALE's MS. One of the reasons for such a difference might have been the Gilbank et al. (2010, 2011) parametrisation of $A_{H\alpha}$ that was used to derive the SFR. However, since it provides larger extinction values for higher stellar masses, it should actually produce a larger slope with respect to CIGALE's MS, which is the opposite of what is observed. The two most striking features in this relation are (i) its low scatter (as already discussed above) and (ii) the lack of a red sequence. These two points can be explained by the fact that we are considering a sample of [OII] emitters, thus star-forming galaxies. Hence, even
galaxies that have accumulated a significant fraction of old stars so that their spectrum is sufficiently red to fall on CIGALE's red sequence should still have enough star formation to produce a bright [OII] doublet. Therefore, they should still appear on the MS when using this SFR tracer. Besides, galaxy-per-galaxy variations in dust extinction that might affect their vertical location on the MS and its scatter are smoothed by the Gilbank et al. (2010, 2011) parametrisation since it only provides an average parametrisation of $A_{H\alpha}$ with stellar mass.

## 3.3 Structure identification and density estimation

To properly constrain the impact of the environment on galaxy properties and galaxy evolution, the environment must be accurately quantified first. To this end, many techniques have been proposed that allow to identify which galaxies belong to structures and which galaxies do not (i.e. field galaxies). Such techniques include among others Voronoi tessellation, Delaunay triangulation (e.g. Marinoni et al. 2002), weighted adaptive kernel smoothing, $n^{th}$ nearest neighbours, or Friends of Friends (FoF) algorithms (e.g. Knobel et al. 2012; Iovino et al. 2016). Some methods are purely based on either photometric or spectroscopic observations whereas others such as the Voronoi tessellation Monte-Carlo mapping technique described in Lemaux et al. (2017, 2022) and Hung et al. (2020, 2021) combine both observations. Because, as discussed in Darvish et al. (2015), different estimators yield similar but still slightly different results, we used multiple estimators to define the structures and estimate their density. I mention that this work was not carried out by myself but by B. Epinat and B. Lemaux. However, I have extensively used the results of the FoF algorithm for the analysis of the scaling relations in Mercier et al. (2022) and in Chapter 6 and for the analysis of the angular momentum in the MAGIC survey (see Chapter 7). Hence, it is important that I provide a brief account of how this was achieved.

### 3.3.1 Structure identification with a FoF algorithm

The first step was to determine which galaxies are associated to structures and which galaxies belong to the field. Following Bacon et al. (2015), we used a FoF algorithm to detect structures in our MUSE fields. Classically, FoF algorithms define structures through an iterative process given two tunable parameters that are the transverse and longitudinal (i.e. along the line-of-sight) separations. In our case, the former is given by the physical distance between galaxies on the plane of the sky and the latter by the difference of MUSE spectroscopic redshifts between two galaxies. When galaxies are gravitationally bound to the same structure, then this redshift difference is related to their relative motion along the LOS in the structure. Hence, the two parameters used for the FoF algorithm to detect structures are a transverse proper (i.e. physical) distance and a LOS velocity separation. The algorithm then iterates by finding the “friends” of each galaxy in the input catalogue, that is galaxies that respect the projected distance and velocity separations criteria, then find the “friends” of “friends” using the same criteria, and so on until the structures are determined. Classically, the two parameters of the FoF algorithm have to be tuned from simulations beforehand to account for the survey target selection and success rate (e.g. Knobel et al. 2012; Iovino et al. 2016). In our case, and as discussed previously, we do not have any pre-selection, we have a high completeness in the [OII] emitters range (i.e. below $z \approx 1.5$) where we want to detect the structures, and we have low uncertainties on the velocity difference between two galaxies ($\Delta v \approx 20 \text{ km s}^{-1}$). Besides, our goal is to characterise small groups, where low-mass galaxies may be found, as much as large groups and clusters. Therefore,
we ran the FoF algorithm on the entire MAGIC sample$^8$, selecting galaxies that have a secure spectroscopic redshift ($CONFID \ge 2$), without any magnitude cut. Based on the prescription given in Knobel et al. (2012) that applies for galaxy groups with more than six members, we used a transverse separation of 450 kpc and a velocity separation of 500 km s$^{-1}$.

### 3.3.2 Local and global density estimators

Once the structures have been identified, we can estimate the density at each location in the FoV. The first way to do so is through the derivation of local density estimators that compute a density estimate locally at a given position in the FoV taking into account the spatial distribution in the structure. In our case, we have used two different estimates that do not rely on the same data. First, we have considered the Voronoi tessellation Monte-Carlo mapping technique (Lemaux et al. 2017, 2022; Hung et al. 2020, 2021) that combines photometric and spectroscopic information but that does not take into account the output of the FoF algorithm. Second, we have derived a somewhat similar but finer density estimate that also uses Voronoi tessellation but for each detected group separately and that uses MUSE and zCOSMOS (Lilly et al. 2007) spectroscopic redshifts. Both methods have advantages and caveats that are discussed below and thus they provide complementary information.

The Voronoi tessellation Monte-Carlo mapping technique is interesting because it can be applied on photometric and spectroscopic data and thus over a larger FoV than our MUSE observations, which allows to mitigate the oversampling of spectroscopic redshifts in our MUSE FoVs. To apply the method, we used photometric information from the COSMOS catalogue of Laigle et al. (2016) and spectroscopic information from the zCOSMOS and VIMOS Ultra-Deep Survey (VUDS) spectroscopic catalogues (Lilly et al. 2007; Le Fèvre et al. 2015). Only these catalogues were considered because, contrary to other catalogues (e.g. COSMOS-Wall, Iovino et al. 2016), they provide a homogeneous sampling in redshift space. The way the method works is by assigning to each source a redshift and an associated uncertainty drawn from a random distribution. Then, redshift slices in the range $0.2 \le z \le 1.5$ are produced with a step of 150 km s$^{-1}$ and for each slice Voronoi tessellations are performed. A first density map is computed as the inverse of the cell's area which is then transformed into a grid with rectangular cells with 75 kpc wide pixels in physical distance. This process is repeated 100 times and the final density is computed as the mean map. Because sampling may differ from one redshift slice to another, it is more useful to use the overdensity $\delta(z)$ rather than the density which is defined as $1 + \delta(z) = \Sigma(z)/\Sigma_{\text{med}}(z)$, where $\Sigma(z)$ is the density map centred at redshift $z$ and $\Sigma_{\text{med}}(z)$ is the median value of the map at the same redshift. The clear advantage of this method is that we can assign an overdensity value for each galaxy. This includes field galaxies but also galaxies that are located near the edges of the MUSE FoV. However, its disadvantages are that (i) it does not use MUSE spectroscopic redshifts and (ii) it only provides a density estimate in a redshift slice rather than associated to a given structure.

Hence, to correct for these two caveats we derived a second local density estimator which uses MUSE spectroscopic redshifts. This estimator also relies on Voronoi tessellation but, in this case, there is no more need to perform Monte Carlo realisations because the galaxies associated to structures in MAGIC all have secure spectroscopic redshifts. The clear advantage of this method compared to the previous one is that we can perform the Voronoi tessellation and compute the local density for each structure separately, therefore taking into account the output of the FoF. In

<footnote>$^8$That is on the entire catalogue rather than on each MUSE field separately.</footnote>
particular, this means that we have a density estimate galaxy-per-galaxy whereas the previous method averaged it over 75 kpc wide pixels. Its caveats are that (i) we do not have a density estimate for field galaxies and (ii) we underestimate the density for galaxies near the sides of the MUSE FoVs because we certainly miss a few galaxies that belong to the structure but that are located beyond the FoVs. To help circumvent the latter point, we also have computed the same estimator but including galaxies in zCOSMOS that are located beyond the FoVs.

Finally, we also consider a last estimator that probes global rather than local environment. This estimator combines the projected distance of a galaxy with respect to the structure's centre ($R_{proj}$) with its line-of-sight velocity with respect to the structure's systemic velocity ($\Delta v$, e.g. Noble et al. 2013; Pelliccia et al. 2019). This estimator writes

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

where $R_{200} = \sqrt{3}\sigma_v/[10H(z)]$ is the radius where the density of the structure is equal to 200 times the critical density of the Universe, with $H(z)$ the Hubble parameter at redshift $z$, and $\sigma_v = c \times \Delta z/(1+z)$ is the velocity dispersion of the structure, with $\Delta z$ the difference in redshift between the galaxy and the structure (see Epinat et al. for more details). Infalling galaxies that are still located far away from the structure's centre and that have large velocities with respect to the systemic velocity will therefore have small $\eta$ values whereas already accreted galaxies located in the inner parts and with low velocity will have larger values. Hence, assuming the density is higher in the central parts, $\eta$ can also be used as a proxy for the density at the galaxy's location.

Because it uses the projected distance and the galaxy's velocity, it is also a complementary estimator to the two previous local estimators. However, $\eta$ requires to define a structure and its centre (in our case taken as the barycentre of the galaxies' associated to the structure by the FoF algorithm) which is not possible for field galaxies. Hence, this estimator is solely restricted to galaxies located in structures found by the FoF algorithm.
# Chapter 4

# Morphological modelling

The main goal of this thesis is to perform the analysis of the impact of the environment on the dynamical properties of intermediate redshift galaxies in the MAGIC survey. In Mercier et al. (2022) and in Chapter 6 I studied the impact of the environment on scaling relations and in Chapter 7 I discuss the angular momentum of galaxies in MAGIC. For both analyses, a crucial aspect was to model the galaxies’ morphology to extract key morphological parameters. This step was performed during the first year of this thesis. There were actually two different, though not entirely estranged, reasons why this modelling was important: (i) to disentangle the contribution to the observed surface brightness of the disk and of the bulge components and (ii) to reliably model the morphology to constrain afterwards the contribution of the stars to the total gravitational potential of the galaxy in order to produce mass models for the kinematics modelling (see Chapter 5 for more details). Indeed, the first point is important because even at intermediate redshift ($z \approx 1$) star-forming disk galaxies tend to host two main components: first, the stellar disk in itself and second, an inner spherical or spheroidal bulge component. Thus, deriving accurate morphological parameters for both components was important, especially in Mercier et al. (2022) where we made use of the fraction of flux found in the bulge with respect to the disk to select a sample of star-forming galaxies and of the disk size to study the size-mass relation. The second point is as important as the first, if not even more, because, as is discussed in Sect. 5.3.3, we used the parameters derived from the morphological modelling to constrain the ionised gas kinematics in order to derive more precise DM fractions.

Thus, in this chapter I describe how I performed the morphological modelling for nearly 900 [OII] emitters in the MAGIC survey. I begin with a presentation of the main disk model that I have been using and its main characteristics in Sect. 4.1. I also go through an important aspect that is sky projection and how it can affect the shape of galaxies seen on the plane of the sky. In particular, I focus on its effect on two different disk models: (i) the razor-thin disk and (ii) the double exponential disk. In Sect. 4.2, I quickly argue, illustrated with an example, why taking into account the contribution of the bulge component during the morphological modelling is important and, in Sect. 4.3, I discuss the multi-component decomposition that was done in Abril-Melgarejo et al. (2021) and Mercier et al. (2022). Lastly, I present in Sect. 4.4 the morphological modelling performed on the MAGIC sample during this thesis and I mention a small python wrapper that I developed at the beginning of my PhD to automatise the modelling with GALFIT for a large number of galaxies.
<figure><img src="image_21.png" /><figcaption>Figure 4.1: Example of an azimuthally averaged surface brightness profile for NGC2776 from Pohlen & Trujillo (2006). Left: r' band SDSS image with the dashed blue line representing the noise limit. Right: surface brightness profile in the g' band (blue triangles) and in the r' band (red triangles), with the best-fit exponential disk profile shown in black.</figcaption></figure>

## 4.1 Morphology of disk galaxies

Various models for the surface brightness distribution of disks of galaxies have been suggested in the literature over the years. Some were based on observations of galaxies in the local Universe e.g. Sérsic 1963) whereas others have been theoretically derived so as to find a potential-density pair which solves Poisson equation (see Sect. 5.2). For instance, this is the case for Kuzmin (1956) solution of the eponymous Kuzmin’s disk.

### 4.1.1 Sérsic model

Among all the surface brightness models one in particular, the Sérsic profile (Sérsic 1963), has been used in numerous studies to represent the surface brightness distribution Σ of both disk-like and elliptical galaxies at various redshifts (e.g. Trujillo et al. 2007; Maltby et al. 2010; Mowla et al. 2019; Costantin et al. 2022). This axisymmetric model is usually written as

$$ \Sigma(R) = \Sigma_{\text{eff}} \exp \left\{ -b_n \left[ \left( \frac{R}{R_{\text{eff}}} \right)^{1/n} - 1 \right] \right\}, \quad (4.1) $$

where $R_{\text{eff}}$, called the effective or half-right radius, corresponds to the distance where the integrated flux reaches half of the total flux, and $\Sigma_{\text{eff}}$ is the surface brightness evaluated at $R_{\text{eff}}$. In Eq. 4.1, the term $b_n$ is a function of the Sérsic index $n$ and is the solution of the following equation (Graham et al. 2005):

$$ \Gamma(2n) = 2\gamma(2n, b_n), \quad (4.2) $$

where $\Gamma$ and $\gamma$ are the complete and lower incomplete gamma functions, respectively. The central surface density of such a model is given by
$$ \Sigma(0) = \Sigma_{eff} e^{b_n}, \qquad (4.3) $$

and the flux F integrated up to radius R can be computed as

$$ F(< R) = 2\pi n \Sigma_{eff} R_{eff}^2 e^{b_n} \gamma \left( 2n, b_n \left( \frac{R}{R_{eff}} \right)^{1/n} \right) / b_n^{2n}. \qquad (4.4) $$

This model is therefore entirely described by three parameters only and has the advantage to represent quite accurately the surface brightness profile of various types of galaxies. Indeed, the disk component of galaxies is well represented by an exponential disk model given by setting n to unity. In this case Eq. 4.1 simplifies to

$$ \Sigma(R) = \Sigma_{eff} \exp \left\{ -b_1 \left[ \frac{R}{R_{eff}} - 1 \right] \right\}, \qquad (4.5) $$

where $b_1 \approx 1.6783$. Equation 4.5 is sometimes written with the disk scale length $R_d = R_{eff,d}/b_1$ instead which corresponds to the e-folding length of the profile.

Another model, usually used to represent elliptical galaxies or the bulge component of disk-like galaxies, is the de Vaucouleurs profile (de Vaucouleurs 1948) and corresponds to a Sérsic profile with n = 4, that is

$$ \Sigma(R) = \Sigma_{eff} \exp \left\{ -b_4 \left[ \left( \frac{R}{R_{eff}} \right)^{1/4} - 1 \right] \right\}, \qquad (4.6) $$

where $b_4 \approx 7.6692$. Additionally, the Sérsic profile also reduces to a Gaussian profile when n = 1/2. Thus, the Sérsic profile is really well suited to fit the global (average) morphology of galaxies independently of whether they are more disk-like or more elliptical. For more complicated cases such as galaxies with bulges, bars, spiral arms, rings, tails, and so on either more complex surface brightness models are required or in some cases multiple Sérsic profiles can be combined together (see Sec. 4.3).

### 4.1.2 Sky projection: the curse of astrophysicists

Generally speaking, what really matters to understand the evolution of galaxies is their intrinsic properties, but these are not directly observed. For instance, understanding in detail the dynamical evolution of a galaxy would require at the very least measurements of its 3D stellar and/or gas distribution as well as its velocity vector at each position in 3D space. Unfortunately, there is an effect that affects virtually any observed quantity and inevitably leads to a loss of crucial information: line-of-sight integration also known as sky projection. Effectively, what this means is that any observed quantity is the combination of the contribution of each emitting material element along the line of sight. This effect is technically speaking an utter hindrance because it completely blurs any spatially resolved information along the line of sight and renders extremely difficult, if not impossible, the inversion process. In practice, sky projection affects any quantity derived from flux measurements such as the morphology of galaxies obtained from photometric observations in a single or in multiple bands. Additionally, it also affects the derived galaxies kinematics. Indeed, because we cannot measure proper motions on the plane of the sky and spatially resolve stars in galaxies at intermediate to high redshift, as can be done with GAIA
<figure><img src="image_22.png" /><figcaption>Figure 4.2: Schematics of the geometry used to integrate a razor-thin and thick disk distributions (see Sects. 4.1.2.1 and 4.1.2.2). The axis $x = x'$ (not shown on the figure but going in to have oriented axes) is the axis of rotation, with $x - y$ the plane of the galaxy and $x' - y'$ the plane of the sky, and $z'$ is the direction of the line-of-sight. The angle $i$ represents the inclination of the galaxy. This geometry leads to the coordinate transformations given by Eqs. 4.11 and 4.12.</figcaption></figure>

in the Milky Way and in its vicinity, we are bound to measure the velocity along the line of sight using the Doppler shift. The shift, as well as the shape of the emission or absorption line that is used to derive the velocity, will be affected by both the stellar light distribution of the emitting material and its velocity vector along the line-of-sight.

The impact of sky projection is thus very important and does somewhat hinders our ability to understand in acute details the dynamics of galaxies. In this sense, and because it is not limited to the sole field of extragalactic astrophysics, there is no overstatement in calling sky projection/line-of-sight integration the curse of astrophysicists. In what follows, we discuss how sky projection can affect the modelling of the morphology of galaxies.

Let us assume that the stellar light emitted in 3D space by a galaxy can be represented by its stellar light distribution $\rho(\vec{r})$, and that the galaxy is located sufficiently far away so that its light rays reach us parallel to each other. If $z'$ represents the line of sight distance with respect to the plane of the sky that passes through the center of the galaxy (see Fig. 4.2 for a schematics representing the geometry of the problem), then the observed surface brightness profile at position $(x', y')$ on the plane of the sky writes

$$ \Sigma_{\text{obs}}(x', y') = \int_{z' \in \mathbb{R}} dz' \rho(\vec{r}). \qquad (4.7) $$

Note that Eq. 4.7 is correct as long as the medium is optically thin, otherwise $\rho(\vec{r})$ must be understood as the stellar light distribution already attenuated by the medium the light rays will have to pass through; for instance by taking into account dust absorption and/or scattering. For a spherically symmetric distribution, that is $\rho(\vec{r}) = \rho(r)$, Eq. 4.7 simplifies to Abel integral equation
which can be inverted, most of the time numerically, to derive $\rho$ as a function of $\Sigma_{\text{obs}}$ (for a derivation of the inverse of Abel integral equation, see Appendix B of Binney & Tremaine 2008).

### 4.1.2.1 Projection for a razor-thin disk

For more complicated distributions, solving analytically Eq. 4.7 is usually not possible, and so is inverting it to recover $\rho$ given $\Sigma_{\text{obs}}$. However, there is a family of distributions, used extensively in the literature, for which Eq. 4.7 admits an analytical solution: razor-thin disks. A disk is said to be razor-thin when it has no thickness, in which case its stellar light distribution and surface brightness profile are related through the following equation

$$ \rho(\vec{r}) = \Sigma(x, y)\delta(z), \qquad (4.8) $$

where $\delta$ is a Dirac distribution, $z$ is the vertical distance with respect to the plane $x - y$, and $\Sigma$ is the surface density in the plane of the disk. If $\rho$ represents the 3D stellar light distribution we will call $\Sigma$ the intrinsic surface brightness distribution and if $\rho$ represents the 3D mass distribution we will rather call $\Sigma$ the surface mass density. Since stellar mass is not a direct observable and because I mainly focussed on modelling the surface brightness and stellar light distribution of galaxies during this thesis, I will exclusively use the term surface brightness throughout the following parts.

If the galaxy is seen face-on, that is $z' = z$ (see Fig. 4.2), then Eq. 4.7 can be easily solved and one finds that $\Sigma_{\text{obs}} = \Sigma$. In other terms, for a face-on razor-thin disk galaxy the observed surface brightness is identical to the intrinsic one. On the other hand, when the galaxy is seen edge-on the integral will be carried out perpendicularly to the $z$ direction so that we will get $\Sigma_{\text{obs}} \propto \delta(z)$. Hence an edge-on razor-thin disk will have its flux distributed along a line (that is a null surface) so that its surface brightness should diverge. For a galaxy which is neither face-on, nor edge-on, we need to consider the geometry of the problem to solve Eq. 4.7. For a razor-thin disk, the geometry is shown in Fig. 4.2 with the sky represented as the plane $x' - y'$, where $x'$ is not shown on the figure since it is orthogonal to the plane $y' - z'$. The disk, contained in the plane $x - y$, is inclined by an angle $i$ around the axis $x = x'$. If we define $r$ as the distance of a point in the $(y', z')$ plane and $\theta$ as the oriented angle between the $r$ and $y'$ axes ($-\pi/2 < \theta < \pi/2$), then we have

$$ y' = r \cos \theta, \qquad z' = r \sin \theta, \qquad (4.9) $$

$$ y = r \cos(\theta - i), \qquad z = r \sin(\theta - i). \qquad (4.10) $$

Since the integral is computed along a line of constant $y'$, we can insert Eq. 4.9 into Eq. 4.10 after developing the cosine and sine terms to get

$$ y = z' \sin i + y' \cos i, \qquad (4.11) $$

$$ z = z' \cos i - y' \sin i. \qquad (4.12) $$

For a razor-thin disk, Eq. 4.7 writes

$$ \Sigma_{\text{obs}}(x', y') = \int_{z' \in \mathbb{R}} dz' \Sigma(x, y)\delta(z). \qquad (4.13) $$
From Eq. 4.12, we have $dz = dz' \cos i$, which means the integral in Eq. 4.13 will only be evaluated at $z = 0$ where $y = y' / \cos i$. Thus, we get

$$ \Sigma_{\text{obs}}(x', y') = \frac{\Sigma(x', y' / \cos i)}{\cos i}. \qquad (4.14) $$

The isophotes of $\Sigma_{\text{obs}}$ on the plane of the sky will be given by the set of points $P = (x', y')$ such that $\Sigma_{\text{obs}}(x', y')$ is constant. From Eq 4.14, the isophotes of $\Sigma_{\text{obs}}$ are therefore given by

$$ \{P = (x', y') | \Sigma(x', y' / \cos i) = C\}, \qquad (4.15) $$

where the braces represent a set (in this case of points) and $C = \Sigma_{\text{obs}}(x', y') \cos i$ is a constant term that defines the isophotal value. If the intrinsic surface brightness follows cylindrical symmetry, that is $\Sigma(x, y) = \Sigma(R)$, then the isophotes are given by the equation

$$ x'^2 + \left( \frac{y'}{\cos i} \right)^2 = R^2 \qquad (4.16) $$

where $R$ is a constant term. This corresponds to an ellipse whose axis ratio $q = b/a$, with $a$ and $b$ respectively the semi-major and semi-minor axes, and ellipticity $e = \sqrt{1 - q^2}$ are respectively given by

$$ q = \cos i \quad \text{and} \quad e = \sin i. \qquad (4.17) $$

In Eq. 4.14, we recover the fact that the observed surface brightness must be equal to the intrinsic one when the galaxy is seen face-on and that it diverges for edge-on galaxies. In what follows, we will see that Eq. 4.17 is not valid any more when considering disks with non-zero thickness.

### 4.1.2.2 Projection of thick disks: the case of the double exponential disk

Razor-thin disks are a very common assumption that is made in many morphological and kinematical studies. Nevertheless there has been increasing evidence that disk-like galaxies tend to have significant thickness profiles, especially when going towards higher redshifts. Indeed, studies of the 3D shape of galaxies as a function of cosmic time have been carried out for instance by van der Wel et al. (2014a) and Zhang et al. (2019) using SDSS and HST data in the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) field. By modelling the apparent axis ratio distribution with an ellipsoid model for the shape of their star forming galaxies they found that galaxy disks become thicker with larger stellar masses and at higher redshift up to $z = 2.5$. Thus, the assumption of razor-thin disks may be inappropriate when modelling the morphology of intermediate to high redshift galaxies and may therefore bias the morphological parameters derived from this simple assumption, especially regarding the disk inclination. In the most general case, in order to derive the surface brightness distribution of a thick disk, one has to solve Eq. 4.7. However, it is common practice to make a few assumptions when modelling the thickness profile of disk-like galaxies. One of the main assumptions that is made is that it is possible to separate the intrinsic surface brightness profile $\Sigma$ from the thickness profile $h$ (e.g. van der Kruit & Searle 1981a; Bizyaev & Mitronova 2002; Mosenkov et al. 2015; Bizyaev et al. 2020), that is the stellar light distribution writes

$$ \rho(\vec{r}) = \Sigma(x, y)h(\vec{r}), \qquad (4.18) $$
<figure><img src="image_23.png" /><figcaption>Figure 4.3: Examples of sky projection for a double exponential disk as a function of the disk inclination: 0° (top left), 25° (top right), 50° (bottom left), and 75° (bottom right). Normalised surface brightness contours are shown at half (orange), one tenth (red), and one hundredth (blue) the central value for single exponential (dashed line), double exponential (solid line), and sech law (dotted line) profiles. All profiles share the same disk scale length of ten pixels, scale height of 2 pixels ($q_0 = 0.2$), and central surface brightness.</figcaption></figure>

where the following integral always holds:

$$ \int_{z \in \mathbb{R}} dz \ h(\vec{r}) = 1. \qquad (4.19) $$

Equation 4.19 is just a different way to say that the integral of the stellar light distribution along the line of sight for a face-on galaxy yields the intrinsic surface brightness profile of the galaxy (that is independent of projection effects) independently of the shape of the thickness profile. Because of Eq. 4.19 h must be integrable and normalised, and thus cannot be constant. At the very least, it has to be a function of z, the vertical distance with respect to the disk x - y plane but, in essence, it could also vary with x and y. The question of the variation of the
thickness profile with the distance R to the centre was discussed as early as van der Kruit & Searle (1981a), though at the time they could not find any significant dependence with R. On the other hand, Bizyaev et al. (2014) did find non-negligible radial gradients of the scale height when modelling a sample of edge-on galaxies in the SDSS. However, given the uncertainties on their modelling and the impact of dust extinction and dust scattering, they were only able to measure significant radial gradients for galaxies with massive bulges. Therefore, as they discuss in their analysis, the radial dependence they measure has probably nothing to do with an intrinsic variation of the disk scale height with R, but more with a contamination of the bulge component. Thus, let us simplify the shape of the thickness profile and assume in what follows that it is only a function of z, that is $h(\vec{r}) = h(z)$. The question we need to ask ourselves now is: what is its shape ? A first answer can be given from the theoretical side by assuming that the distribution of light along z at a given position R corresponds to an isothermal sheet: this is the so-called sech law (e.g. Spitzer 1942). Its expression can be derived from the Jean's equations in cylindrical coordinates (see Appendix C.4 for a quick derivation) and is given by

$$h(z) = \frac{1}{4h_z} \operatorname{sech}^2 \left( \frac{z}{2h_z} \right), \qquad (4.20)$$

where we have defined $h_z = z_0/2$ the scale height of the disk, with $z_0$ defined in Appendix C.4, and where we have normalised Eq. C.23 so that it obeys Eq. 4.19. In the case of a thin disk or far enough from the disk plane, that is when $|z| \gg h_z$, we have $\operatorname{sech}(x) \approx 2 \exp\{-|x|\}$ and the sech law simplifies to (van der Kruit & Searle 1981b)

$$h(z) = \frac{1}{2h_z} \exp \left\{ \frac{-|z|}{h_z} \right\}, \qquad (4.21)$$

where we have also normalised it to obey Eq. 4.19. When combined with an exponential disk for the surface brightness distribution we obtain the so-called double exponential disk:

$$\rho(R, z) = \frac{\Sigma(0)}{2h_z} \exp \left\{ -\frac{R}{R_d} - \frac{|z|}{h_z} \right\}, \qquad (4.22)$$

and we call $q_0 = h_z/R_d$ the intrinsic axis ratio of the disk.

The sky projection of a double exponential disk is described in details in Appendix E of Mercier et al. (2022), so we directly provide below the solution. Assuming an optically thin medium, the sky projected surface brightness distribution of a double exponential disk writes

$$\Sigma_{\text{obs}}(x', y') = \frac{\Sigma(0)}{2q_0 \sin i} \int_{\mathbb{R}} dv \exp \left\{ -\sqrt{\alpha^2 + v^2} - |\beta v - \gamma| \right\}, \qquad (4.23)$$

where $v = y/R_d$, with y related to the line-of-sight coordinate $z'$ through Eq. 4.11, and where we have defined the following variables:

$$\alpha = \frac{x'}{R_d}, \quad \beta = \frac{1}{q_0 \tan i}, \quad \gamma = \frac{y'}{h_z} \left( \sin i + \frac{\cos^2 i}{\sin i} \right). \qquad (4.24)$$

Generally speaking, Eq. 4.23 cannot be solved analytically and has to be integrated numerically. Nevertheless there exists an analytical solution along the minor axis $y'$ given by
$$ 
\begin{aligned}
\Sigma_{\text{obs}}(0, y') &= \frac{\Sigma(0)}{2q_0 \sin i} \int_{\mathbb{R}} dv \exp \{-|v| - |\beta v - \gamma|\} \\
&= \frac{\Sigma(0)}{q_0 \sin i} \frac{e^{-\gamma} - \beta e^{-\gamma/\beta}}{1 - \beta^2}.
\end{aligned} 
\quad (4.25) 
$$

We show in Fig. 4.3 an example of sky projection for a double exponential disk with $R_d = 10$ pixels and $h_z = 2$ pixels computed through numerical integration. As a comparison, we show surface brightness contours for this profile (solid line) and for a single exponential disk (dashed line) and a sech law (dotted line) both with similar parameters. When viewed face-on, the three profiles match perfectly as can be expected from Eq. 4.19. However, the more inclined the disk the larger the difference between the single exponential disk model (i.e. razor-thin) and the double exponential disk becomes. In particular, we see that for a given set of intrinsic disk model parameters the projected single exponential disk will appear more concentrated than the double exponential disk. Furthermore, the isophotes of a projected double exponential disk are not elliptical contrary to a razor-thin disk since the intrinsic thickness of the disk will contribute to the shape of the isophotes on top of the sky projection. Nevertheless, even up to quite high inclinations ($i \le 80^{\circ}$) the elliptical isophote approximation remains acceptable. Finally, we can see that for any surface brightness contour and inclination value the sech law profile gives almost the same isophotes as the double exponential disk, which shows that the approximation which was made to derive Eq. 4.21 was appropriate.

### 4.1.2.3 From single to double exponential disk

<figure><img src="image_24.png" /><figcaption>Figure 4.4: Ratio of the intrinsic central surface brightness to the value measured by a single exponential fit as a function of the intrinsic disk inclination $i_0$ and axis ratio $q_0$. This ratio was computed using the approximation of an oblate spheroid system to relate the intrinsic axis ratio to the intrinsic inclination.</figcaption></figure>

Rarely will a surface brightness fit be performed by numerically integrating a 3D model along the line-of-sight since it can be prone to numerical artefacts, might suffer from degeneracies, and
would certainly increase the computation time. For all these reasons, 2D fitting softwares such as GALFIT (see Sect. 4.4.1 for a quick description) only provide surface brightness distributions which can be projected onto the sky assuming a razor-thin disk geometry. Nevertheless, we saw in Sect. 4.1.2.2 that a double exponential profile projected onto the sky has isophotes which are nearly elliptical. Therefore, if one knows or at least assumes a given axis ratio, it might be possible to recover the intrinsic disk parameters of the double exponential profile given a single exponential razor-thin disk fit onto the observed surface brightness distribution.

If one knows $q_0$, then there are three main parameters that need to be recovered: the intrinsic central surface brightness $\Sigma(0)$, the disk scale length, and the intrinsic inclination of the disk $i_0$. Estimating the impact on the disk scale length is not straightforward since it would involve integrating Eq. 4.23 to compute the e-folding length of the projected surface brightness distribution, which can only be done numerically and would depend upon the other parameters.

Nevertheless, for face-on or mildly inclined galaxies, the disk scale length from the single exponential fit should give a fairly reasonable estimate of the intrinsic value. Regarding the central surface brightness, it is possible to compute the analytical correction that is required to go from the single exponential fit value to the intrinsic one. To do so, one must (i) compute the relation between the observed central surface brightness and the intrinsic one ($\Sigma_{\text{exp}}(0)$) for a single exponential fit, which is given by evaluating Eq. 4.14 at $x' = y' = 0$, and (ii) compute the same relation but in the case of a projected double exponential disk, which is given by evaluating Eq. 4.25 at $y' = 0$. The correction writes

$$ \frac{\Sigma(0)}{\Sigma_{\text{exp}}(0)} = \frac{q_0 \sin i_0 + \cos i_0}{q}, \qquad (4.26) $$

where $q$ and $q_0$ are respectively the observed (from the single exponential fit) and intrinsic axis ratio. Equation 4.26 can only be solved if the intrinsic axis ratio and inclination of the disk are known a priori but, if the former can usually be assumed, the latter is usually what one wants to derive. Thus, to compute the central surface brightness correction, one needs to derive the correction for the inclination first. Similarly to the disk scale length, there is no analytical formula that links the observed axis ratio to the intrinsic inclination given the other parameters.

Therefore, in Mercier et al. (2022) we approximated the dependence of the intrinsic inclination to the observed and intrinsic axis ratio by considering that the double exponential disk could be approximated by an oblate spheroid with similar $q_0$. Then, it follows from Bottinelli et al. (1983) that we have the following relation:

$$ \cos^2 i_0 = \frac{q^2 - q_0^2}{1 - q_0^2}, \qquad (4.27) $$

in which case Eq. 4.26 reduces to

$$ \frac{\Sigma(0)}{\Sigma_{\text{exp}}(0)} = \frac{q_0 \sin i_0 + \cos i_0}{\sqrt{q_0^2 \sin^2 i_0 + \cos^2 i_0}}. \qquad (4.28) $$

I show in Fig. 4.4 the ratio from Eq. 4.28 as a function of the intrinsic axis ratio and inclination. For perfectly face-on and edge-on galaxies, the measured central surface brightness from a single exponential fit is not biased. However, for other inclinations and for non-zero intrinsic thickness values the central surface brightness derived from a single Sérsic fit will underestimate the intrinsic value. The correction that needs to be applied will vary with $i_0$ depending on $q_0$ and reach a maximum value equal to $\sqrt{2}$.
## 4.2 Impact of bulges

Assuming the disk component of galaxies can be approximated by a razor-thin disk, their inclination can be derived by either measuring the ellipticity of their isophotes and then using Eq. 4.17 or by directly fitting a sky projected 2D model to the observed surface brightness using Eq. 4.14. In particular, a common practice (e.g. Trujillo et al. 2007; Longhetti et al. 2007; Kormendy et al. 2009; Law et al. 2012; van der Wel et al. 2012, 2014b) is to fit a single Sérsic profile with a free Sérsic index *n* to account for both disk-like (*n* ≈ 1) and elliptical (*n* ≈ 4) galaxies, as well as those with in-between shapes. However, this approximation may break as soon as the contribution of other components to the overall surface brightness distribution can no longer be neglected, in particular that of the central bulge. Therefore, fitting a single Sérsic profile to a galaxy with a non-negligible bulge contribution may bias the measure of the disk parameters. To get a rough estimate of how and how much the bulge component might affect the results of a single Sérsic fit performed to recover the disk parameters, I generated a grid of galaxy models using GALFIT (for more details and for an example of a modelling performed with GALFIT see Sect. 4.4). Each model was generated by combining an exponential disk profile for the disk component with a symmetrical de Vaucouleurs profile to represent the surface brightness distribution of a spherically symmetric bulge, and was convolved with a PSF represented by a circular Moffat profile with a FWHM of 3 pixels and a powerlaw slope of 1.9, similar to the PSF used to model the morphology of galaxies from HST F814W images in the MAGIC survey (Abril-Melgarejo et al. 2021; Mercier et al. 2022). Once the models generated, they were fitted with a single Sérsic profile with all its parameters let free, including its Sérsic index *n*. Because using a single bulge-disk model would not be representative of the type of systematics which may arise when fitting with a single Sérsic profile, I generated instead a grid of models with different bulge and disk parameters. The most important parameters which control the bulge-to-total flux ratio (B/T) are the bulge and disk total magnitudes and effective radii. In particular, only the difference between the bulge total magnitude and that of the disk will affect B/T rather than their exact values. Therefore, without loss of generality, I set the disk total magnitude to 20, with the magnitude photometric zero point set to an arbitrary value of 30, and varied the other aforementioned parameters. The following grid was used:

* Bulge total magnitude: 18, 20, 20.5, 21, 22, and 23 dex.

* Bulge effective radius: 15, 10, 5, 3, and 2 pixels.

* Disk effective radius: 30, 25, 20, 15, 10, and 5 pixels.

* Disk axis ratio: 0.17, 0.5, 0.77, 0.94, and 1, which corresponds to razor-thin disk inclinations of 0, 20, 40, 60, and 80°, respectively.

The reason for choosing these values was twofold: (i) to be consistent with the typical sizes and B/T values found in the morphological analysis of the MAGIC survey (see Chapter 3), (ii) so that there would be approximately 25% of the models falling in one of the following B/T bins: B/T < 15%, 15% ≤ B/T ≤ 50%, 50% ≤ B/T ≤ 75%, and B/T > 75%. In total, this corresponds to 900 bulge-disk models generated and then fitted with GALFIT. We show in Fig. 4.5 the variation of various best-fit parameters as a function of the true disk axis ratio and B/T. The parameters for each individual fit are shown as small markers whereas the median value for each B/T and axis ratio bins are shown as large markers with a dashed line. From left to right and top
<figure><img src="image_25.png" /><figcaption>Figure 4.5: Variation of the best-fit parameters for various composite galaxy models with non negligible bulge contribution as a function of the true disk axis ratio and bulge-to-total flux ratio (B/T). Each model was fitted with GALFIT using a single Sérsic profile with free Sérsic index. The B/T values for each model are shown as small markers whereas the median values for each B/T are shown as large markers with dashed lines. From left to right and top to bottom: (i) difference between the true total magnitude (disk and bulge, noted *mag<sub>tot</sub>*) and that of the best-fit model (mag), (ii) ratio between the model effective radius ($R_{eff}^{GF}$) and the true global effective radius (noted $R_{eff}$, see Sect. 4.3.2 for its definition), (iii) difference between the true disk magnitude ($mag_d$) and that of the best-fit model, (iv) ratio between the model effective radius and the true disk effective radius $R_{eff,d}$, (v) best-fit Sérsic index, and (vi) ratio between the model axis ratio $q$ and the true disk axis ratio $q_d$. Even if the global parameters (total magnitude and global effective radius) are quite well recovered, a single Sérsic fit overestimates the flux contained in the disk and its axis ratio while it underestimates its size, especially for large bulge contributions.</figcaption></figure>
to bottom (i) difference between the true total magnitude (noted mag<sub>tot</sub>), including the disk and the bulge components, and that of the best-fit model (noted mag), (ii) ratio between the single Sérsic model effective radius R<sub>eff</sub><sup>GF</sup> and the true global effective radius R<sub>eff</sub> which corresponds to the radius which contains half of the total light, that is including the disk and the bulge components (see Sect. 4.3.2 for how it is computed in practice), (iii) difference between the true disk magnitude (mag<sub>d</sub>) and that of the best-fit model, (iv) ratio between the model effective radius and the true disk effective radius R<sub>eff,d</sub>, (v) best-fit Sérsic index, and (vi) ratio between the model axis ratio q and the true disk axis ratio q<sub>d</sub>

The first row of Fig. 4.5 shows how well does a single Sérsic fit recover the global properties of a bulge-disk model that are the total magnitude and the global effective radius. On average a single Sérsic fit recovers almost the same total magnitude and finds an effective radius which is slightly biased towards larger values though still close to the correct one. As we go towards larger B/T, then the scatter in these two plots increases significantly even though the median values are not as much affected. The second and third rows focus on the other hand on the efficiency of using the single Sérsic fit parameters as proxies for the disk properties. We can see that, for small bulge contributions, a single Sérsic fit converges towards a near exponential model with small bias and scatter for each parameter. But, the larger the bulge contribution the more biased and the more scattered the parameters become. This is especially relevant for the axis ratio since it is a morphological parameter that needs to be precisely constrained to properly model the dynamics of galaxies (see Sect. 5.3 for a discussion on that topic). Indeed, we see that the more elliptical the disk component (i.e. low q<sub>d</sub>), the more underestimated is its axis ratio estimate q when fitting with a single Sérsic profile.

## 4.3 Multi-component decomposition

### 4.3.1 Bulge-disk decomposition

One way to circumvent the biases which arise when modelling with a single Sérsic profile the surface brightness distribution of galaxies is to use multiple components. When doing so, the most general form the total surface brightness Σ can take is

$$ \Sigma_{\text{obs}}(x, y) = \sum_{i} \Sigma_{\text{obs},i}(x, y), \qquad (4.29) $$

where Σ<sub>obs,i</sub> corresponds to the projected surface brightness distribution of the component i. When the multi-component decomposition simplifies to a combination of a bulge and a disk, a common practice (e.g. Laurikainen et al. 2010; Abril-Melgarejo et al. 2021; Chu et al. 2021; Bouché et al. 2022; Costantin et al. 2022; Mercier et al. 2022; Quilley & de Lapparent 2022) is to model each component as an axisymmetric Sérsic profile. To remove degeneracies the disk is usually chosen to be an exponential disk and the bulge can either be a Sérsic profile with free Sérsic index or a de Vaucouleurs profile. In the latter case, the total surface brightness writes

$$ \Sigma(R) = \Sigma_{\text{eff,d}} \exp \left\{ -b_1 \left[ \frac{R}{R_{\text{eff,d}}} - 1 \right] \right\} + \Sigma_{\text{eff,b}} \exp \left\{ -b_4 \left[ \left( \frac{R}{R_{\text{eff,b}}} \right)^{1/4} - 1 \right] \right\}, \quad (4.30) $$

where (Σ<sub>eff,d</sub>, R<sub>eff,d</sub>) and (Σ<sub>eff,b</sub>, R<sub>eff,b</sub>) are the disk and bulge Sérsic parameters, respectively. For disk-like galaxies, the bulge component will only dominate in the inner parts. Nevertheless, this
<figure><img src="image_26.png" /><figcaption>Figure 4.6: Variation of the flux distribution (black lines) and bulge-to-total flux ratio (B/T, blue lines) with radius for a weak bulge (left-hand plot) and a strong bulge (right-hand plot). For each plot, two cases are shown: (i) a small bulge with $R_{eff,b}/R_{eff,d} = 0.1$ represented as thin lines and (ii) a larger bulge with $R_{eff,b}/R_{eff,d} = 0.5$ represented as thick lines. The bulge flux distribution is shown as dotted lines and the total flux distribution (bulge and disk) as solid lines. The same disk model was used throughout and its flux distribution is represented as a thin dashed line in both plots. The bump visible for the thick blue line on the right-hand plot is due to the fact that the bulge is sufficiently bright and large so that its flux distribution reaches its plateau after that of the disk, as is shown on the inset in the right plot.</figcaption></figure>

does not mean that the bulge contribution becomes null at large radii. Indeed, using Eq. 4.4 we can write the bulge-to-total flux ratio integrated up to radius R as

$$B/T(R) = \frac{B/D(R)}{1 + B/D(R)}, \qquad (4.31)$$

where the bulge-to-disk ratio B/D is given by

$$B/D(R) = \frac{4b_1^2 \Sigma_b(0) R_{eff,b}^2 \gamma \left(8, b_4 \left(\frac{R}{R_{eff,b}}\right)^{1/4}\right)}{b_4^8 \Sigma_d(0) R_{eff,d}^2 \gamma \left(2, b_1 \frac{R}{R_{eff,d}}\right)}. \qquad (4.32)$$

From the series expansion of the incomplete gamma function evaluated for small values of x given by $\gamma(s, x) \approx x^s/s$, we get that B/T near R = 0 is given by

$$B/T(R \to 0) = [1 + \Sigma_d(0)/\Sigma_b(0)]^{-1}. \qquad (4.33)$$
Thus, the bulge-to-total flux ratio computed at the centre is just a different way to write the ratio of the central surface brightnesses between the disk and the bulge components. At infinity, B/T is given by

$$B/T(R \to \infty) = \left[ 1 + \frac{2b_4^8 \Sigma_d(0) R_{eff,d}^2}{(8!)b_1^2 \Sigma_b(0) R_{eff,b}^2} \right]^{-1}, \quad (4.34)$$

where $2b_4^8/[(8!)b_1^2] \approx 210.75$. Therefore, the only way to have a null bulge contribution at large distances is that either the disk contribution is infinite while that of the bulge is finite, or that the central surface brightness and/or effective radius of the bulge component are equal to zero, that is that there is no bulge at all. Similarly, even if the bulge usually dominates in the inner parts (i.e. $\Sigma_b(0) > \Sigma_d(0)$), this does not necessarily mean that B/T will reach a value close to one. However, in practice bulges have quite small radii compared to that of the disk component<fnref n="1" />. This means that, most of the time, galaxies with a non-negligible bulge contribution at large radii will have B/T almost equal to one near the centre.

A few examples of B/T values computed at different distances for various bulge and disk parameters are shown in Fig. 4.6. The B/T values are shown as blue solid lines whereas the black lines represent the flux distribution of the bulge component (dotted), the disk component (dashed), and the sum of the two (solid). The left-hand plot represents the case of a galaxy with a weak bulge contribution and the right plot shows the case of a strong bulge contribution. In both plots, thick lines represent galaxies with a large bulge ($R_{eff,b}/R_{eff,d} = 0.5$) and thin lines galaxies with a small bulge ($R_{eff,b}/R_{eff,d} = 0.1$). As previously stated and as shown in Fig. 4.6, for typical galaxies without massive bulges, B/T gets close to zero far enough from the centre. In the case of the most massive bulge with the largest size (thick blue line, right plot), we can see a bump in B/T after $R \approx 2R_{eff,d}$ which means the contribution of the bulge component slightly increases again after this radius. This bump is actually physical and not an artefact (assuming such massive and large bulges would exist in nature in the first place). The explanation is that, near the centre the bulge dominates, as can be seen in Fig. 4.6, but at the same time the disk contribution increases since its inner slope is larger than that of the bulge and therefore B/T decreases. However, because the bulge is large, it reaches the plateau of its flux distribution after that of the disk component. Thus what happens is that, till the disk has not reached its plateau, B/T decreases but, once it is on the plateau, the flux of the bulge keeps mildly rising while that of the disk remains constant, therefore slightly increasing the value of B/T before the bulge reaches its own plateau and B/T remains nearly constant afterwards.

### 4.3.2 Global effective radius

An interesting quantity to derive is the radius $R_n$ which encloses nth percent of the total flux. If measured on the plane of the sky, then this can easily be done using an isophotal fitting approach (e.g. Jedrzejewski 1987; Stone et al. 2021), albeit taking care of the background signal and of the flux from nearby objects. However, we know that galaxies are composed of multiple components which have different geometries (e.g. spherical bulge or razor-thin disk). Isophotal fitting will therefore bias the measurement of the radius, especially when galaxies have a non-negligible bulge contribution.

<footnote n="1">¹At least in the case of the MAGIC survey, see for instance Fig. A.2 of Mercier et al. (2022) or in Chapter 6.</footnote>
This issue disappears when one fits the surface brightness distribution with a multi-component decomposition because the geometry of each component is directly taken into account into the fit, which means that the recovered individual surface brightness profiles are intrinsic. Assuming a multi-component decomposition as described in Eq. 4.29, cylindrical symmetry for the individual surface brightness profiles (i.e. $\Sigma_i(x, y) = \Sigma_i(R)$), and that all components share the same centre, one can find $R_n$ by solving the following equation

$$ \sum_i \int_0^{R_n} dRR \, \Sigma_i(R) = \frac{n}{100} \times \sum_i \int_{\mathbb{R}} dRR \, \Sigma_i(R). \quad (4.35) $$

When $n = 50$, the solution to this equation is called the global half-light or effective radius (noted $R_{eff}$). Generally speaking, Eq. 4.35 can be solved by computing the two integrals numerically. However, in the case of a bulge-disk decomposition similar to Eq. 4.30, one can use the expression of the flux of a Sérsic profile integrated within a radius $R$ given by Eq. 4.4 to derive the following analytical form:

$$ F_{\text{tot,d}} \left[ \gamma \left( 2, b_1 \frac{R_{\text{eff}}}{R_{\text{eff,d}}} \right) - 0.5 \right] + \frac{F_{\text{tot,b}}}{\Gamma(8)} \left[ \gamma \left( 8, b_4 \left( \frac{R_{\text{eff}}}{R_{\text{eff,b}}} \right)^{1/4} \right) - \frac{\Gamma(8)}{2} \right] = 0, \quad (4.36) $$

where $F_{tot,b}$ and $F_{tot,d}$ are the bulge and disk total fluxes, respectively. If one defines the total magnitude of a component $i$ as $\text{mag}_i = -2.5 \log_{10} F_{tot,i} + \text{zpt}$, where zpt is a zero point that is the same for all the components and normalises by the total flux, then Eq. 4.36 writes

$$ f(R_{\text{eff}})/f(\infty) = 0, \quad (4.37) $$

where the function $f$ is defined as

$$ f(x) = 10^{-\text{mag}_d/2.5} \left[ \gamma \left( 2, b_1 \frac{x}{R_{\text{eff,d}}} \right) - 0.5 \right] + 10^{-\text{mag}_b/2.5} \left[ \gamma \left( 8, b_4 \left( \frac{x}{R_{\text{eff,b}}} \right)^{1/4} \right) / \Gamma(8) - 0.5 \right], \quad (4.38) $$

and Eq. 4.37 can be solved by searching for its unique zero located between $R_{eff,d}$ and $R_{eff,b}$. Indeed, if $R_{eff}$ is larger than both $R_{eff,d}$ and $R_{eff,b}$, then the flux at $R_{eff}$ will be the sum of $F_d(R_{eff}) > F_{d,tot}/2$ and $F_b(R_{eff}) > F_{b,tot}/2$. Thus, its value will be larger than the expected $F_{tot}/2$ value. Therefore, $R_{eff}$ cannot be greater than max($R_{eff,d}$, $R_{eff,b}$), and the same argument can be given when $R_{eff}$ is less than both $R_{eff,d}$ and $R_{eff,b}$. Additionally, it can be proven that there is a unique zero that is a solution of Eq. 4.37 by noticing that $f$ is a monotonously increasing function of $x$ whose normalised form $f(x)/f(\infty)$ is bounded between -1 for $x = 0$ and 1 for $x = \infty$.

## 4.4 Morphological modelling of galaxies in the MAGIC sample

### 4.4.1 Modelling with Galfit

In Mercier et al. (2022), we modelled the surface brightness distribution of more than a thousand galaxies in the MAGIC sample by performing a bulge-disk decomposition using the morphological
modelling tool GALFIT<fnref n="2" /> (Peng et al. 2002). GALFIT is a powerful, fast, and versatile forward
modelling tool that can fit any number of components on images using various models that can be
projected onto the sky. A custom PSF model can be provided by the user and masks can be
applied to reduce the influence of nearby objects on the fit. Running GALFIT only requires to
write an input file (called FEEDME) which contains the models, the initial values of their
parameters, and a few metadata. A second optional file (called CONSTRAINTS) can be given to
provide additional constraints between parameters. We provide below an example of typical
FEEDME and CONSTRAINTS files that were used to model the galaxies in the MAGIC sample.
In addition to the bulge and disk components we also added a sky background model which
produces a sky gradient in both x and y directions. This model can be useful to remove the
contribution to the disk and bulge components of a potential sky background in the image or that
of a nearby source which would induce a sky gradient in its direction. For all the components, two
values must be provided for each parameter, except for the centre position which needs four (two
for the x coordinate and two for the y coordinate). The first value corresponds to the initial guess
and the second one defines whether the parameter must be fixed (0) or free (1). In our case, the
free parameters are (i) the centre positions, (ii) the total magnitudes, (iii) the effective radii,
(iv) the disk axis ratio, (v) the disk PA, and (vi) the sky model parameters. On the other hand,
the Sérsic indices of the disk and the bulge were fixed to *n* = 1 and *n* = 4 and the bulge axis ratio
and PA were fixed to *q* = 1 and PA = 0°, respectively. Furthermore, a few constraints were
placed: (i) the offset between the disk and the bulge components was fixed to its initial value.
Because both models were initialised to the same position, in practice this means that we forced
the disk and the bulge to share the same centre, (ii) to limit potential degeneracies, we forced the
centre to only be able to move in a box whose side length is equal to ten pixels, (iii) we set limits
on the total magnitudes and effective radii of the disk and bulge components. These limits forced
GALFIT to try different bulge-disk combinations before potentially converging to a single disk or
single bulge model. If a galaxy has a negligible disk or bulge component, GALFIT will try to find
the combination of parameters which minimises its contribution to the surface brightness
distribution. In practice, this means its total magnitude will converge to the maximum value and
its size is likely to converge to either its upper or lower bound. With an upper magnitude of
100 dex for a magnitude zero point of 30 dex, this effectively means that for our HST data the
model will have a null contribution to the surface brightness and to the total flux.

<footnote n="2">²https://usersobs.carnegiescience.edu/peng/work/galfit/galfit.html</footnote>
<table><thead><tr><th>Image Component</th><th>Name</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>1</td><td># IMAGE PARAMETERS</td><td></td><td></td></tr><tr><td>2</td><td></td><td></td><td></td></tr><tr><td>3 A)</td><td>input.fits</td><td></td><td># Input data image (FITS file)</td></tr><tr><td>4 B)</td><td>output.fits</td><td></td><td># Name for the output image</td></tr><tr><td>5 D)</td><td>psf.fits</td><td></td><td># Input PSF image for convolution</td></tr><tr><td>6 G)</td><td>example_constraints</td><td></td><td># File with parameter constraints</td></tr><tr><td>7 H)</td><td>0 134 0 134</td><td></td><td># Image region to fit (xmin xmax ymin ymax)</td></tr><tr><td>8 I)</td><td>134 134</td><td></td><td># Size of convolution box (x y)</td></tr><tr><td>9 J)</td><td>30.0</td><td></td><td># Magnitude photometric zeropoint</td></tr><tr><td>10 O)</td><td>regular</td><td></td><td># Display type (regular, curses, both)</td></tr><tr><td>11 P)</td><td>0</td><td></td><td># 0=normal run; 1,2=make model/imgblock;</td></tr><tr><td>12</td><td></td><td></td><td># 3=normal run and separate components</td></tr></tbody></table>

<table><thead><tr><th>Disk Component</th><th>Name</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>15</td><td># Disk component</td><td></td><td></td></tr><tr><td>16 O)</td><td>sersic</td><td></td><td># Object type</td></tr><tr><td>17 1)</td><td>70.0 70.0 1 1</td><td></td><td># position x, y [pixel]</td></tr><tr><td>18 3)</td><td>22.0 1</td><td></td><td># total magnitude</td></tr><tr><td>19 4)</td><td>15.00 1</td><td></td><td># effective radius R_e [pixel]</td></tr><tr><td>20 5)</td><td>1.00 0</td><td></td><td># Sersic exponent</td></tr><tr><td>21 9)</td><td>0.30 1</td><td></td><td># axis ratio (b/a)</td></tr><tr><td>22 10)</td><td>90.0 1</td><td></td><td># position angle [Degrees: Up=0, Left=90]</td></tr><tr><td>23 Z)</td><td>0</td><td></td><td># Skip this model in output image?</td></tr><tr><td>24</td><td># Bulge component</td><td></td><td></td></tr><tr><td>25</td><td></td><td></td><td></td></tr><tr><td>26 O)</td><td>sersic</td><td></td><td># Object type</td></tr><tr><td>27 1)</td><td>70.0 70.0 1 1</td><td></td><td># position x, y [pixel]</td></tr><tr><td>28 3)</td><td>22.0 1</td><td></td><td># total magnitude</td></tr><tr><td>29 4)</td><td>5.00 1</td><td></td><td># effective radius R_e [pixel]</td></tr><tr><td>30 5)</td><td>4.00 0</td><td></td><td># Sersic exponent</td></tr><tr><td>31 9)</td><td>1.00 0</td><td></td><td># axis ratio (b/a)</td></tr><tr><td>32 10)</td><td>0.0 0</td><td></td><td># position angle [Degrees: Up=0, Left=90]</td></tr><tr><td>33 Z)</td><td>0</td><td></td><td># Skip this model in output image?</td></tr><tr><td>34</td><td># sky</td><td></td><td></td></tr><tr><td>35</td><td></td><td></td><td></td></tr><tr><td>36 O)</td><td>sky</td><td></td><td># Object type</td></tr><tr><td>37 1)</td><td>1.00 1</td><td></td><td># sky background [ADU counts]</td></tr><tr><td>38 2)</td><td>0.000 1</td><td></td><td># dsky/dx (sky gradient in x)</td></tr><tr><td>39 3)</td><td>0.000 1</td><td></td><td># dsky/dy (sky gradient in y)</td></tr><tr><td>40 Z)</td><td>0</td><td></td><td># Skip this model in output image?</td></tr></tbody></table>

Code 4.1: Example of GALFIT FEEDME file.

<table><thead><tr><th>Component</th><th>Parameter</th><th>Constraint</th></tr></thead><tbody><tr><td>1</td><td>x</td><td>offset</td></tr><tr><td>2</td><td>1_2</td><td>offset</td></tr><tr><td>3</td><td>1_2</td><td>offset</td></tr><tr><td>4</td><td>1</td><td>-5.0 5.0</td></tr><tr><td>5</td><td>1</td><td>-5.0 5.0</td></tr><tr><td>6</td><td>1</td><td>10.0 to 100.0</td></tr><tr><td>7</td><td>2</td><td>10.0 to 100.0</td></tr><tr><td>8</td><td>1</td><td>0.5 to 100.0</td></tr><tr><td>9</td><td>2</td><td>0.5 to 100.0</td></tr></tbody></table>

Code 4.2: Example of GALFIT CONSTRAINTS file.
4.4.2 Automatic modelling for a large number of galaxies

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

Figure 4.7: Example of a current summary file for galaxy 83_CGr23 generated by my python wrapper for automatic modelling of a large number of galaxies with GALFIT. From left to right: HST image, GALFIT best-fit model, and residuals. Flux units are arbitrary. This type of summary file is automatically generated at the end of the fitting process for each object that was fitted and was used in Mercier et al. (2022) to quickly identify galaxies that needed improved morphological models.

One of my first steps was to model the morphology of a large sample of mainly star forming galaxies from the MAGIC survey. From the 2730 objects in the initial MAGIC catalogue, including nearby stars, intermediate, and high-redshift galaxies ($z > 1.5$), we selected 1142 [O II] emitters with reliable spectroscopic redshifts<fnref n="3" />. The bulge-disk decomposition performed on these 1142 [O II] emitters is similar to what was described in Sect. 4.4.1 and a full description can be found in Mercier et al. (2022).

Since GALFIT was designed to model in details the morphology of galaxies by manually creating configuration files and iterating through the models and their parameters, using it directly on such a large number of galaxies is not entirely well suited. Thus, to accelerate and automatise the modelling process I developed a python wrapper around GALFIT<fnref n="4" /> which allows to easily model any number of galaxies using a common set of models with similar initial parameter values. This wrapper was built to be as versatile as possible and therefore is not restricted to bulge-disk decompositions only, but can be used with any combination of models provided in GALFIT and for any kind of photometric data. At the end of the modelling, summary pdf files are created to quickly visualise the modelling and assert whether a model needs to be improved or not. An example of such a summary file for galaxy 83_CGr23 is shown in Fig. 4.7. Because of the way it was built, this wrapper makes it as easy to model a large number of galaxies as single objects without having to manually tune any configuration file. This code was successfully used in Mercier et al. (2022) where I performed a bulge-disk decomposition on HST ACS F814W images of the 1142 [O II] emitters in the MAGIC sample. After removing low signal-to-noise ratio (S/N) objects and very small objects, we managed to get reliable bulge-disk decompositions for 890 galaxies<fnref n="5" />.

<footnote n="3">³[O II] emitters are defined as intermediate redshift galaxies for which we might be able to detect the [O II]λλ3727, 3729 doublet in the MUSE wavelength range, that is, located in the redshift range $0.25 \le z \le 1.5$.</footnote>

<footnote n="4">⁴https://github.com/WilfriedMercier/wilfried/tree/master/galfit</footnote>

<footnote n="5">⁵The morphological models can be found as supplementary material of Mercier et al. (2022) at the following address: https://www.aanda.org/articles/aa/olm/2022/09/aa43110-22/aa43110-22.html</footnote>
Below is shown an example of a python code which uses this GALFIT wrapper to produce FEEDME and CONSTRAINTS files similar to what is illustrated in the code snippets 4.1 and 4.2, which runs GALFIT onto those files, and which produces a pdf summary file once the fit is done:
```python
1 from galfit import run_galfit
2
3 # Sky model
4 sky = {'name': 'sky'}
5
6 # Header properties:
7 #   - do not include the input/output and constraint file names
8 #   - psf must be the same for all objects
9 header = {'psfImage': 'psf.fits'}
10
11 # Sersic models:
12 #   - objects will have the same initial parameters
13 disk = {'name': 'sersic', 'x': 70, 'y': 70, 'n': 1, 'fixedParams': ['n'],
14        're': 15, 'mag': 22, 'bOvera': 0.3, 'PA': 90}
15
16 bulge = {'name': 'sersic', 'x': 70, 'y': 70, 'n': 4, 'fixedParams': ['n', 'PA'],
17        'bOvera'], 're': 5, 'mag': 22}
18
19 # Constraints
20 absRangeM = {'type':'absoluteRange', 'value':[10, 100]}
21 absRangeRe = {'type':'absoluteRange', 'value':[0.5, 100]}
22 relRange  = {'type':'relativeRange', 'value':[-5, 5]}
23 offset    = {'type':'offset', 'value':0.0}
24
25 constraint = [('components': [1, 2], 'parameter':'x',
26               'constraint': offset),
27               {'components': [1, 2], 'parameter':'y',
28               'constraint': offset},
29               {'components': 1, 'parameter':'x',
30               'constraint': rRange},
31               {'components': 1, 'parameter':'y',
32               'constraint': rRange},
33               {'components': 1, 'parameter':'mag',
34               'constraint': absRangeM},
35               {'components': 2, 'parameter':'mag',
36               'constraint': absRangeM},
37               {'components': 1, 'parameter':'re',
38               'constraint': absRangeRe},
39               {'components': 2, 'parameter':'re',
40               'constraint': absRangeRe},
41           ]
42
43 inputNames = ['input.fits']
44 outputNames = ['output.fits']
45 feedmeNames = ['input feedme']
46 constraintNames = ['input_constraints']
47
48 run_galfitfeedmeNames, header=header, pathIn='/'
49         listProfiles = [disk, bulge, sky],
50         inputNames = inputNames,
51         constraintNames = constraintNames,
52         outputNames = outputNames,
53         constraints = constraint)
```

Code 4.3: GALFIT python wrapper code example.
# Chapter 5

# Dynamical modelling

For the dynamical analysis of the galaxies in the MAGIC survey performed in Abril-Melgarejo et al. (2021) and Mercier et al. (2022) (see also Chapter 6), and for the analysis of the angular momentum in Chatper 7, I had to derive the galaxies’ ionised gas kinematics, more specifically the kinematics that is traced by the [O II] doublet, being one of the brightest emission lines visible in the MUSE cubes in the redshift range $0.2 \le z \le 1.5$. For this, I used the moment maps formalism described in Epinat et al. (2010) that assumes that the gas is located in a rotating razor-thin disk. For the star-forming galaxies at intermediate redshift that I studied, this approximation is appropriate because such galaxies mostly have disk geometries. Besides, the gas being collisional it is expected that, with enough time, it should collapse into a rotating disk through a redistribution of angular momentum. This hypothesis is therefore common in kinematics studies found in the literature (e.g. Epinat et al. 2010; Übler et al. 2017; Bouché et al. 2022). There are at least two major difficulties with such a modelling: (i) the line-of-sight integration (as for the morphology) that removes the 3D information that would be required to have a complete modelling of the 3D gas or stellar kinematics and (ii) the low spatial and/or spectral resolution of our MUSE data around and beyond $z=1$ which limits the complexity of the models we can use. Therefore, to perform such models we must rely on assumptions. The razor-thin disk geometry of the gas is certainly one of them that is appropriate for most galaxies. The fact that we only consider rotation during the modelling is another. Some authors have nevertheless started investigating new models beyond the simple rotating disk assumption (e.g. Bizyaev et al. 2017; Wang & Lilly 2022; Lopez-Coba et al. 2022) but these remain nevertheless fairly limited yet.

Thus, I describe in this chapter the different steps that I took to produce and fit resolved kinematics map models for almost 600 [O II] emitters using a forward modelling method that involves mass models. I start in Sects. 5.1 and 5.2 by describing the very general theoretical background behind the kinematics and mass modelling that I performed. These two sections, even if they have been used to derive the ionised gas kinematics only, are technically speaking applicable to stellar kinematics as well. More specifically, in Sect. 5.1 I derive the kinematics of a test particle under the influence of generic forces, especially in the context of cylindrical coordinates. In Sect. 5.2, I turn towards mass modelling and I show how the kinematics of a test particle can be related to the gravitational potential of a given mass distribution under the assumption of equilibrium. I discuss the effect of the asymmetric drift correction on the rotation velocity and I provide detailed examples of mass models that were used in Mercier et al. (2022). I
also quickly mention whether the choice of parametrisation for the DM halo profile has an important impact or not on our modelling of the velocity fields for our MUSE data of intermediate redshift galaxies and on the derived DM mass fraction. In Sect. 5.3, I explain how the kinematics maps extraction works and how I did it for the MAGIC sample. I also describe the theory behind the moment maps formalism that is used to model the galaxies’ velocity fields, taking into account the instrumental effects on both the velocity fields and the velocity dispersion maps. Finally, I conclude this chapter by mentioning how the modelling was performed in practice for the MAGIC survey and how it changed between the analyses performed in Abril-Melgarejo et al. (2021) and Mercier et al. (2022) and that of the angular momentum presented in Chapter 7.

## 5.1 Kinematics of a test particle

To begin with, we will consider a test particle whose position $\vec{r}(t)$ can vary in time under the influence of various forces. For instance, the motion of the test particle can be impacted by the gravitational potential generated by the mass distribution of the galaxy within which it is embedded, as well as that of nearby objects, but also by other kind of forces such as those generated by outflows from AGN feedback and/or star formation, pressure gradients, and so on. In very general terms, the kinematics of a test particle can be described by Newton’s equation

$$ \frac{\partial \vec{V}}{\partial t}(\vec{r}) = \vec{F}(\vec{r}), \qquad (5.1) $$

where $\partial \vec{V}/\partial t$ is the total acceleration, with $\vec{V}$ the velocity vector, and $\vec{F}$ is the sum of the different forces contributing to the kinematics. In Eq. 5.1 and in what follows we always make the dependence on time implicit. Since Eq. 5.1 is too general to be used directly, in what follows we will simplify it by making some assumptions about the forces at play and the geometry of the problem.

### 5.1.1 Kinematics in cylindrical coordinates

Because of the symmetries found in some galaxies, in particular their disk component, it is common practice, without loss of generality, to write the velocity vector in cylindrical coordinates $(R, \theta, z)$. In these coordinates, $R$ and $\theta$ represent the radial distance and azimuthal angle in a plane at fixed $z$, respectively, and $z$ represents the vertical distance with respect to the plane $z = 0$. In this coordinate system, and with the corresponding basis vectors $(\hat{e}_R, \hat{e}_\theta, \hat{e}_z)$, we can write the velocity vector as

$$ \vec{V}(\vec{r}) = V_R(\vec{r}) \hat{e}_R + V_\theta(\vec{r}) \hat{e}_\theta + V_z(\vec{r}) \hat{e}_z, \qquad (5.2) $$

where $V_R(\vec{r}) = \dot{R}$ corresponds to radial motions, $V_\theta(\vec{r}) = R\dot{\theta}$ to rotation and $V_z(\vec{r}) = \dot{z}$ to vertical motions. Similarly, the total acceleration writes

$$ \frac{\partial \vec{V}}{\partial t}(\vec{r}) = a_R(\vec{r}) \hat{e}_R + a_\theta(\vec{r}) \hat{e}_\theta + a_z(\vec{r}) \hat{e}_z, \qquad (5.3) $$

where the three components $a_R$, $a_\theta$ and $a_z$ can be written as
$$ 
\begin{align} 
a_R &= \dot{V}_R - V_\theta^2/R, \nonumber \\ 
a_\theta &= V_R V_\theta/R + \dot{V}_\theta, \nonumber \\ 
a_z &= \dot{V}_z, \tag{5.4}
\end{align} 
$$

and where we have made implicit the dependence on $\vec{r}$. A few points are worth mentioning. To begin with, we see that vertical forces (i.e. along $\hat{e}_z$) will only affect the vertical velocity component, whereas radial and azimuthal forces can affect both the radial and azimuthal components of the velocity vector. Besides, if all the forces cancel each other out along the vertical axis it only implies that $V_z$ will be constant in time, but not necessarily in space. If we consider the special case where the contribution of forces cancel each other out along $\hat{e}_\theta$, for instance to study the kinematics of a test particle located in an axisymmetric disk, and if we further assume that the radial velocity component is constant (potentially null) then the velocity of the test particle is purely circular and given by

$$ V_{\theta} = \sqrt{-R a_{R}}. \tag{5.5} $$

We note that Eq. 5.5 remains correct even if there are azimuthal forces, that is along $\hat{e}_\theta$, as long as $\dot{V}_R = 0$. However, then there would also be a non-zero velocity component along the radial direction to take into account and/or $V_\theta$ would become time-dependent (see the second line of Eq. 5.4).

## 5.2 Mass modelling

In this section, we assume that we want to derive the kinematics of a test particle embedded in a galaxy. Generally speaking, this galaxy could be isolated or have companion objects such as other galaxies in its vicinity. Each object will have its own mass distribution, made of various constituents (e.g. stars, gas, dust, DM halo, etc.), which can be decomposed into multiple components (e.g. bulge, disk, halo, etc.). Each component of each galaxy will have its own mass distribution $\rho_{M,i}$, which will generate a gravitational potential $\Phi_i$ following Poisson equation

$$ \nabla^2 \Phi_i(\vec{r}) = -4\pi G \rho_{M,i}(\vec{r}), \tag{5.6} $$

where $G$ is the gravitational constant and $\nabla^2$ is the Laplacian operator. For a given mass distribution $\rho_{M,i}$, assuming one can solve Eq. 5.6, then the corresponding acceleration $\vec{a}_i = \vec{F}_i/m$, with $m$ the mass of the test particle, is given by

$$ \vec{a}_i(\vec{r}) = -\nabla \Phi_i(\vec{r}), \tag{5.7} $$

where $\nabla$ is the gradient operator. Thus, to derive the kinematics of a test particle in a given mass distribution, one has to solve Eq. 5.5 in combination with Eq. 5.7. Because the total mass distribution $\rho_M$ is nothing more than the sum of the individual mass distributions, so will be the total gravitational potential $\Phi$, that is we have

$$ \rho_M(\vec{r}) = \sum_i \rho_{M,i}(\vec{r}) \quad \text{and} \quad \Phi(\vec{r}) = \sum_i \Phi_i(\vec{r}). \tag{5.8} $$
### 5.2.1 Mass modelling for a kinematics tracer in a plane

Because we wrote the equation of motion of a test particle in cylindrical coordinates in Sect. 5.1.1, it can be useful to also rewrite Eq. 5.7 in the same way, for instance if the kinematics tracer that is used is located in a plane. This is given by

$$ \vec{\nabla}\Phi(\vec{r}) = \frac{\partial\Phi}{\partial R}(\vec{r}) \hat{e}_R + \frac{1}{R} \frac{\partial\Phi}{\partial\theta}(\vec{r}) \hat{e}_{\theta} + \frac{\partial\Phi}{\partial z}(\vec{r}) \hat{e}_z. \quad (5.9) $$

The radial acceleration will be given by $-\partial\Phi/\partial R$, so that if we assume that $\dot{V}_R = 0$ we find the solution

$$ V_{\theta}(\vec{r}) = \sqrt{R \frac{\partial \Phi}{\partial R}(\vec{r})}. \quad (5.10) $$

Note that $V_{\theta}$ appearing in Eq. 5.10 is sometimes noted $V_c$ or $V_{\text{circ}}$ instead and referred to as the circular velocity. Thus, using Eq. 5.8 we get the law of velocity addition for the circular velocity:

$$ V_{\theta}^{2}(\vec{r}) = \sum_{i} V_{\theta,i}^{2}(\vec{r}), \quad (5.11) $$

where $V_{\theta,i}$ is the circular velocity that is solution of Eq. 5.10 for the mass distribution $\rho_{M,i}$.

### 5.2.2 Asymmetric drift correction

Equation 5.10 is a common starting point for kinematic studies using mass models since it relates the theoretical rotation or circular velocity to a given mass distribution. However, it is also common to find in such studies a correction for the circular velocity called asymmetric drift. This corrective term is included so as to take into account the impact of a radial pressure gradient to the kinematics. The asymmetric drift correction to the circular velocity can be derived from Jean's equations under the main assumption of cylindrical symmetry (e.g. see Chapter 4 of Binney & Tremaine 2008). Nevertheless, its mathematical form can be recovered without invoking Jean's equations. Following Dalcanton & Stilp (2010) and Bouché et al. (2022), the acceleration $\vec{a}_P$ generated by such a pressure gradient can be written as

$$ \vec{a}_P(\vec{r}) = -\frac{1}{\rho_M} \frac{\partial P}{\partial R}(\vec{r}) \hat{e}_R, \quad (5.12) $$

with $P$ the pressure and where the negative sign comes from the fact that a positive pressure gradient will induce an inward force, and a negative gradient will induce an outward force. In this case, Eq. 5.10 is rewritten as

$$ V_{\theta, \text{corr}} = \sqrt{R \frac{\partial \Phi}{\partial R}(\vec{r}) - \frac{R}{\rho} \frac{\partial P}{\partial R}(\vec{r})}, \quad (5.13) $$

where $V_{\theta, \text{corr}}$ is the asymmetric drift corrected circular velocity, that is the rotation velocity that is affected by both the gravitational potential and the pressure gradient. Noting that $d\log R = dR/R$ and assuming that $P = \rho_M \sigma_R^2$, with $\sigma_R$ the velocity dispersion along the radial direction, we get

$$ V_{\theta, \text{corr}}^2 = V_{\theta}^2 + V_{\text{ac}}^2, \quad (5.14) $$
where $V_{ac}$ is the asymmetric drift correction given by

$$V_{ac} = \sqrt{-\sigma_R^2 \left[ 2 \frac{\partial \log \sigma_R}{\partial \log R} + \frac{\partial \log \rho_M}{\partial \log R} \right]}, \quad (5.15)$$

and where we have made implicit the dependence on $\vec{r}$. Equations 5.14 and 5.15 are the most general expressions of the asymmetric drift correction usually found in the literature<fnref n="1" />. To determine the correction, one has to know a priori the mass density of the kinematics tracer used (stars, ionised or cold gas), but also the radial component of the velocity dispersion. Because the dispersion is usually measured along the line of sight from spectroscopic data (noted $\sigma_V$), deriving $\sigma_R$ is not straightforward. A common assumption is isotropy, that is $\sigma(\vec{r}) = \sigma_R = \sigma_V$. We note that in this expression $\sigma$ only corresponds to the velocity dispersion related to the pressure gradient. However, in practice, the velocity dispersion which can be measured from kinematics data is the result of multiple physical processes (e.g. stellar or AGN feedback, accretion, mergers, etc.). Thus, directly using $\sigma_V$ in Eq. 5.14 might lead to an overestimate of the asymmetric drift correction, which in turn could lead to biases in the estimate of the circular velocity.

### 5.2.2.1 Asymmetric drift correction for thick disks

In the case where the mass distribution of the kinematics tracer is located in a disk as given by Eq. 4.18, Eq. 5.15 simplifies to

$$V_{ac} = \sqrt{-\sigma_R^2 \left[ 2 \frac{\partial \log \sigma_R}{\partial \log R} + \frac{\partial \log \Sigma_M}{\partial \log R} + \frac{\partial \log h_z}{\partial \log R} \right]}, \quad (5.16)$$

where $\Sigma_M$ is the mass density of the kinematics tracer in the plane of the disk, that is the equivalent of the surface brightness distribution but in terms of mass. For a double exponential disk model or even a sech law the thickness profile $h_z$ is independent of $R$ so that its derivative in Eq. 5.16 vanishes. To further simplify the expression of the asymmetric drift correction, the assumption of constant velocity dispersion is usually made, in which case Eq. 5.16 becomes

$$V_{ac} = \sigma_R \sqrt{-\frac{\partial \log \Sigma_M}{\partial \log R}}. \quad (5.17)$$

If the mass density is described by a Sérsic profile then Eq. 5.17 is simply given by

$$V_{ac} = \sigma_R \sqrt{\frac{b_n}{n} \left( \frac{R}{R_{eff}} \right)^{1/n}}, \quad (5.18)$$

where $R_{eff}$ is the effective radius of the Sérsic profile. This simplifies even further for an exponential disk as it becomes

$$V_{ac} = \sigma_R \sqrt{\frac{R}{R_d}}, \quad (5.19)$$

<footnote n="1">¹Similar expressions can be derived from Jeans' equations by assuming that $\sigma_\phi = \sigma_R$ and $\overline{v_R v_z} = 0$. See Appendix C for the notation and Sect. 4.8.2 of Binney & Tremaine (2008) for the derivation.</footnote>
where $R_d$ is the disk scale length of the kinematics tracer mass distribution. Other expressions for the asymmetric drift correction that depend on other hypotheses ($h_z$ and/or $\sigma_R$ dependent on $R$, spheroidal mass distribution, etc.) can be found in the literature. These expressions go beyond the scope of the kinematics modelling and analysis performed in Mercier et al. (2022) (see Sect. 6) and in this Thesis and therefore we refer the interested reader to Appendix A of Bouché et al. (2022) for more details and references.

### 5.2.3 Circular velocity for a spherical distribution of mass

Spherical mass distributions are useful first approximations used to describe multiple components in galaxies, whether that be a bulge or a stellar or DM halo. Such distributions are interesting from a kinematical point of view because their circular velocity can be easily derived. Indeed, let us write Poisson’s equation in cylindrical coordinates, assuming the gravitational potential is spherically symmetric as should be appropriate for distributions:

$$ \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial \Phi}{\partial r} \right) = 4\pi G \rho_M. \qquad (5.20) $$

Multiplying by $r^2$ and integrating along $r$ from the center to some radius $R$ we get

$$ R^2 \left. \frac{\partial \Phi}{\partial r} \right|_R = GM(<R), \qquad (5.21) $$

where $M(<R)$ is the integrated mass up to radius $R$. If Eq. 5.21 is evaluated in the plane of the galaxy disk, that is at $z=0$, then we recover on the left hand side of the equation $R$ times the circular velocity so that, for any spherically symmetric mass distribution, it writes

$$ V_{\theta} = \sqrt{\frac{GM(<R)}{R}}. \qquad (5.22) $$

We note that Eq. 5.22 can also be used inversely, that is to compute the integrated mass within some radius given $V_{\theta}$ as long as the mass distribution remains spherically symmetric.

### 5.2.4 Examples of mass models

There exist multiple mass models that can be used to model the dynamics of galaxies. Some can be derived from analytical mass density-gravitational potential pairs (e.g. Kuzmin’s disk) whereas others must be computed numerically for a given mass distribution. In the latter case and for an arbitrary mass density, one has to find first the corresponding gravitational potential and then compute the circular velocity. In Mercier et al. (2022), we used and described four different mass models: (i) a spherically symmetric Hernquist’s model (Hernquist 1990) to describe the bulge component of galaxies, (ii) a razor-thin exponential disk model, also referred to as Freeman’s disk (Freeman 1970), (iii) a double exponential disk, and (iv) a spherically symmetric NFW (Navarro et al. 1996) halo for the DM distribution. For details about the derivation of the circular velocity for each mass distribution we invite the reader to refer to Appendix D of Mercier et al. (2022) which can also be found in Chapter 6 of this Thesis. An example adapted from galaxy 104_CGR79 of the shape of the circular velocity, mass density, and cumulative mass profiles for the Hernquist’s model, double exponential disk, and NFW halo is shown in Fig. 5.1.
### 5.2.4.1 Hernquist bulge

The Hernquist’s model is an interesting mass distribution that is used to describe galaxy bulges since, as it was shown in Hernquist (1990), once projected onto the plane of the sky, and under the assumption of spherical symmetry, it fits quite well the usual circular de Vaucouleurs profile used to model the surface brightness distribution of the bulge component in galaxies. It writes

$$ \rho_M(r) = \frac{M_b a}{2\pi r (r+a)^3}, \qquad (5.23) $$

where $r$ is the radial distance with respect to the center of the profile, $M_b$ the total bulge mass and $a = r_{1/2,M}/(1 + \sqrt{2})$ its scale radius with $r_{1/2,M}$ the half-mass size. Integrating Eq. 5.23 to get the mass and then using Eq. 5.22, the circular velocity in the plane of the disk writes

$$ V_b(R) = 2V_{b,max} \frac{\sqrt{ar}}{a+r}, \qquad (5.24) $$

where $V_{b,max} = 0.5 \times \sqrt{GM_b/a}$ is the maximum circular velocity reached at $r = a$. A caveat of Eqs. 5.23 and 5.24 is that they require $M_b$ and $a$ to be known which are parameters of the mass density rather than parameters of the surface brightness distribution which is what is observed/fitted in practice. Therefore, in Mercier et al. (2022) we computed a mapping between Hernquist and de Vaucouleurs parameters by generating a grid of Hernquist’s models whose parameters are appropriate for the galaxies in the MAGIC sample, sky projected them, and then fitted them with de Vaucouleurs profile. The observed correlation between the two Hernquist parameters $F_b = M_b/\Upsilon$ and $a, \Upsilon$ being the mass-to-light ratio (M/L) assumed to be constant throughout the galaxy, and the two de Vaucouleurs parameters $\Sigma_{eff}$ and $R_{eff}$ were then fitted with the two following scaling relations:

$$ \log_{10} a \, [\text{kpc}] = -0.454 + 0.725 \log_{10} R_{eff,b} \, [\text{kpc}] \qquad (5.25) $$

$$ \log_{10} F_b / \Sigma_{eff,b} \, [\text{cm}^2] = 1.194 + 1.750 \log_{10} R_{eff,b} \, [\text{kpc}]. \qquad (5.26) $$

For typical bulge sizes around 1 – 2 kpc the relative difference between the true maximum bulge induced circular velocity and that inferred using the mapping given in Eqs. 5.25 and 5.26 is, at its lowest, below 10%. On the other hand, bulges with very small and very large sizes should have larger relative differences, however such bulges are typically associated with small fluxes since they are either too large to be physical or too small to be resolved in our HST data.

### 5.2.4.2 Freeman’s disk

A very common disk model used in the literature is the so-called Freeman’s disk. This model assumes a razor-thin disk with an exponential surface brightness distribution (see Eq. 4.8) and has an analytical expression for its circular velocity that was derived by Freeman (1970) based on a previous work by Toomre (1963). The velocity writes

$$ V_{RT}(R) = V_{RT,max} \times \frac{y f(y)}{1.075 f(1.075)}, \qquad (5.27) $$

where the subscript RT stands for razor-thin, $f(y) = \sqrt{I_0(y)K_0(y) - I_1(y)K_1(y)}$, $y = R/(2R_d)$, and with $I_i$ and $K_i$ modified Bessel functions of the first and second kind, respectively. In Eq. 5.27,
the maximum circular velocity is reached at $R = 2.15R_d$ and can be found numerically by finding the zero of the derivative of the rotation velocity as given in Freeman (1970). It is equal to

$$V_{\text{RT,max}} = 2.15f(1.075) \sqrt{\pi G R_d \Sigma_M(0)}, \quad (5.28)$$

with $\Sigma_M(0) = \Upsilon\Sigma(0)$ the central surface mass density of the disk. The factor 1.075 in Eqs. 5.27 and 5.28 just comes from the fact that the function $f(y)$ is evaluated at the radius of maximum velocity.

<figure><img src="image_28.png" /><figcaption>Figure 5.1: Illustration of the various mass models discussed in Sect. 5.2.4. In each figure is represented a Hernquist profile (red line, see Sect. 5.2.4.1), a double exponential disk (blue dashed line, see Sect. 5.2.4.3 and Eq. 4.22), and a NFW profile (black dotted line, see Sect. 5.2.4.4). On the top left are shown the rotation curves, on the bottom left the mass density profiles, and on the right hand-side the cumulative mass profiles. The mass density profile for the double exponential disk is taken in the plane of the disk. The total mass of the bulge and disk component are indicated with red and blue arrows, respectively. The NFW profile reaches a virial mass $M_{200} \approx 2 \times 10^{12} M_{\odot}$ at its virial radius $R_{200} \approx 47R_{22}$, making the galaxy DM-dominated with a DM fraction of about 91% at $R_{200}$.</figcaption></figure>

### 5.2.4.3 Double exponential disk

In the case of a double exponential disk as described by Eq. 4.22, there is no known analytical formula for the circular velocity. However, an analytical approximation was derived in the case of
thin disks, that is with small scale height $h_z$, in Bovy J. online book “Dynamics and Astrophysics of Galaxies. Princeton University Press, Princeton, NJ” (in preparation)<fnref n="2" />:

$$V_d^2(R) = V_{RT}^2(R) - V_{corr,max}^2 \times \frac{R \ e^{1-R/R_d}}{R_d}, \qquad (5.29)$$

The maximum of the correction $V_{corr,max}$ is reached at $R_d$ and is given by

$$V_{corr,max} = \sqrt{2\pi G h_z \Sigma_M(0)/e}. \qquad (5.30)$$

Equations 5.29 and 5.30 give a fairly reasonable approximation to the real rotation curve generated by the potential of a double exponential disk in the plane of the disk for small $h_z$. For most of the radial range this approximation gives a circular velocity that is off from the real value by around 2%, except in the inner parts where the difference rises significantly. More specifically, the main caveat of Eq. 5.29 is that $V_d^2$ becomes negative as soon as the correction factor becomes greater than the razor-thin disk circular velocity squared appearing on the right hand side of the equation, which always happens in the inner parts. This is not too problematic if one just sets the circular velocity to zero once this happens but the issue is now that the circular velocity reaches zero before reaching the center of the galaxy, which is unphysical for axisymmetric models. The straightforward solution to this problem would be to numerically compute the circular velocity without relying on Eq. 5.29 but (i) with a much longer computational time and (ii) at the cost of numerical stability given that it would involve two integrals with Bessel functions (see Eq. D.11 in Mercier et al. 2022).

The solution I adopted was to apply a correction in the inner parts to the approximation given by Eq. 5.29. Indeed, near the center the rotation curve of a flattened mass distribution is expected to rise quickly so that the approximation of a linearly rising rotation curve in the inner parts should be appropriate enough and would certainly be better than assuming a null velocity throughout. This derivation is presented in Appendix D.5 of Mercier et al. (2022) (see also Sect. 6) and was done by finding a scaling relation between the intrinsic axis ratio $q_0 = h_z/R_d$ of the galaxy and the radius $R_0$ at which the tangential line to the rotation curve given by Eq. 5.29 passes through the center of the galaxy. The relation reads

$$\log(R_0/R_d) = 0.76679 + 0.86230x - 0.13703x^2 - 0.02308x^3 + 0.00452x^4 + 0.00102x^5, \quad (5.31)$$

where $x = \log q_0$.<fnref n="3" /> The choice of a fifth order polynomial fit was based mostly on the shape of the observed relation and on the need to keep a low error (typically below 2%) in the range of $q_0$ values we were interested in. Therefore, the rotation curve for the double exponential disk we adopted in Mercier et al. (2022) writes

$$V_{DE} = \begin{cases} V_d(R_0) \times R/R_0 & \text{if } R < R_0, \\ V_d(R) & \text{otherwise,} \end{cases} \qquad (5.32)$$

where the subscript DE stands for double exponential, $R_0$ is solution of Eq. 5.31, and $V_d$ is the circular velocity derived by J. Bovy and given by Eqs. 5.29 and 5.30. On top of what has been previously said, Eq 5.32 also has the advantage of being continuous and differentiable everywhere.

<footnote n="2">²See Eq. 8.73 in Chapter 8 at https://galaxiesbook.org/.</footnote>

<footnote n="3">³Note that there are a few typos in Appendix D.5 of Mercier et al. (2022). The equation to solve to find $R_0$ given $q_0$ is indeed Eq. D.24 of the paper with $y_0 = R_0/(2R_d)$. The correct expression for the solution is the one given in Eq. 5.31 of this Thesis.</footnote>
### 5.2.4.4 NFW dark matter profile

Multiple models of DM exist in the literature either cores or cusps in their central parts, for instance the NFW (Navarro et al. 1996), Einasto (Einasto 1965), Burkert (Burkert 1995), or Dekel-Zhao (Freundlich et al. 2020) profiles just to cite a few. In Mercier et al. (2022) we have decided to use a spherically symmetric NFW profile to model the DM halo component that has strong effects on the dynamics of galaxies. This mass density profile writes

$$ \rho(r) = \delta_c \rho_{\text{crit}} \frac{r_s}{r(1 + r/r_s)^2}, \qquad (5.33) $$

where $r_s = r_{200}/c$ is the scale radius, with $r_{200}$ the virial radius of the halo where the mean overdensity is equal to 200 and $c$ the halo concentration, $\rho_{\text{crit}} = 3H_0^2/(8\pi G)$ the Universe closure density and $\delta_c$ the halo characteristic overdensity (for more details, see Navarro et al. 1996). The corresponding circular velocity writes

$$ V_h(r) = \frac{V_{h,\text{max}}}{0.46499} \sqrt{\frac{\log(1 + r/r_s)}{r/r_s} - \frac{1}{1 + r/r_s}}, \qquad (5.34) $$

where $V_{h,\text{max}}$ is the maximum circular velocity reached at $r \approx 2.163 r_s$ (similarly to the Freeman's disk it has to be derived numerically). This DM halo model has the disadvantage that, if integrated to infinity, its total mass diverges. Thus, its total mass as well as its extent are to be taken at the virial radius of the halo where the structure is supposed to be the outermost bound. Hence, one has to be careful when using an NFW profile to not integrate it beyond its virial radius when computing its dynamical effect. Hopefully for us, DM haloes are found to be sufficiently large so that their extent goes beyond the galaxies' optical radius where the ionised gas kinematics is effectively measured in our MUSE data.

### 5.2.4.5 Of the choice of DM profiles for dynamical models of intermediate redshift galaxies

Throughout the years it has become increasingly obvious that the dynamics and the observable matter content of galaxies cannot be reconciled unless an additional DM component is added. Two main families of DM exist in the literature, differing in the slope in the inner parts. The first, derived from DM-only simulations, is called cuspy and it typically corresponds to models such as the NFW profile that have a non-zero slope in their inner parts so that the DM density is infinite at the center. The second family corresponds to core profiles and is favoured by observations, for instance the pseudo isothermal sphere (Begeman 1987). A third family, widely used in the literature to study the cusp-core problem, are parametric models with varying inner and/or outer slopes. Typical parametric DM profiles are that of Einasto, Zhao (Zhao 1997), or Dekel-Zhao. These profiles behave very differently near the center even though they all share a similar behaviour at large radii with constant or slowly declining rotation curves. Thus, it could be expected that their impact on the dynamics of galaxies is significant, especially in the inner parts. However, we argue that the choice of the DM profile, and therefore the shape of the halo near the center, is not particularly significant to model the kinematics of intermediate redshift galaxies given our MUSE data. This point is particularly relevant in this case and in the study performed in Mercier et al. (2022) since the goal was never to study the shape of the DM haloes but rather to compute their mass far enough from the central parts.
To illustrate that point, one of the most massive disk-dominated galaxies with a non negligible bulge contribution was selected from the MAGIC sample: galaxy ID 104\_CGR79. This galaxy was found to have massive disk and bulge components of $M_d(R_{22}) \approx 7.0 \times 10^{10}$ M$_\odot$ and $M_b(R_{22}) \approx 2.5 \times 10^{10}$ M$_\odot$, respectively, with a dominating DM component$^4$. Its size and mass therefore make it a perfect candidate to check the impact of different DM parametrisations on the goodness of fit of the velocity field. To do so, the velocity field extracted from the MUSE cube was fitted using a mass model combining a double exponential disk with intrinsic axis ratio $q_0 \approx 0.2$ (see Eq. 9 of Mercier et al. 2022 for more details on how to derive the disk thickness), a Hernquist bulge, and four different types of spherically symmetric DM haloes: (i) a cuspy NFW profile, as was done in Mercier et al. (2022), (ii) a parametric Einasto model which can be either cuspy or cored, (iii) a pseudo isothermal sphere which is a cored profile, and (iv) a Hubble modified model (Reynolds 1913; Hubble 1930; Rood et al. 1972; Binney & Tremaine 2008) which is cored as well. I skip the details of how the kinematics modelling was performed or how the kinematic maps were extracted from the MUSE cube since these will be covered in details in Sect. 5.3.

The results of the four different fits are shown in Fig. 5.2. The first row shows the various best-fit DM rotation curves for the NFW (solid line), Einasto (dashed line), pseudo isothermal sphere (dotted line), and Hubble modified (dashed-dotted line) profiles. As an indication, I also show the disk (blue solid line) and bulge (red solid line) rotation curves used in the fit, the radius $R_{22}$ where the dynamical mass was measured (vertical blue dashed line), and the extent of the MUSE data ($R_{last}$ - vertical red solid line). The second and third row represent the best-fit velocity fields and residuals, respectively. Almost all the rotation curves display a similar inner slope, except for NFW. This comes from the fact that only the NFW profile is cuspy (the Einasto profile converges towards a cored profile with $\alpha \approx 0.6$). Beyond $R_{last}$, where the rotation curves are extrapolated, the profiles show different behaviours depending on whether they decline with radius (NFW, Einasto, and Hubble modified) or if they remain constant (pseudo isothermal sphere).

Nevertheless, the velocity fields and their residuals are strikingly similar no matter the difference in the inner parts. Similarly, the DM masses computed within $R_{22}$ are nearly identical with $M_{DM} \approx 1.49, 1.47, 1.47$, and $1.53 \times 10^{11}$ M$_\odot$ for the NFW, Einasto, pseudo isothermal sphere, and Hubble modified models, respectively. The main explanation for the fact that these profiles yield similar velocity fields and DM masses is the impact of instrumental effects on the forward velocity field model that smooths out the steep velocity gradient in the inner parts. This effect, also known as beam smearing, will be discussed more attentively in Sect. 5.3.

It is important to note that even though all these models provide similar fits and dynamical masses at $R_{22}$ for this galaxy, some might be favoured with respect to the others using additional independent constraints. One of such constraints that could be used is the stellar-to-halo mass relation (SHMR) that links the stellar mass of the galaxies to their typical halo mass. This relation can be constrained either from simulations (e.g. abundance matching or halo occupation distribution) or from empirical models of galaxy formation (for a review on the topic see Wechsler & Tinker 2018). In this specific case, extrapolating the best-fit dark matter profiles to their virial radius provides virial masses ($M_{vir}$) of the order of $10^{12}$ M$_\odot$, except for the pseudo isothermal sphere that reaches rather $7 \times 10^{12}$ M$_\odot$. Furthermore, this gives a stellar-to-halo mass ratio $M_*/M_{vir} \approx 6\%$ for all the profiles (except the pseudo isothermal sphere) which is slightly above the value of 1-2% that is expected for a halo of mass $M_{vir} \approx 10^{12}$ M$_\odot$ given the various SHMRs at $z \approx 1$ found in the literature (e.g. Guo et al. 2010; Behroozi et al. 2013; Kravtsov et al. 2018).

<footnote>$^4$Approximately $1.5 \times 10^{11}$ M$_\odot$ at $R_{22}$ based on a fit performed in Mercier et al. (2022) with a NFW profile for the DM component</footnote>
Only the pseudo isothermal sphere model does recover a mass ratio around 1% consistent with the SHMR. Hence, if the objective is to probe the shape of DM halo profiles at intermediate and high redshift, the SHMR could be used in combination with such forward kinematic models to favour or disfavour some profiles (e.g. Bouché et al. 2022, where it was used to show evidence for core profiles at $z \approx 1$). Alternatively, authors such as Li et al. (2020) have also tried to directly implement the SHMR as a prior in their Bayesian framework when modelling galaxy rotation curves. Because the constraint is directly implemented in the modelling when marginalising over the fitted parameters, it allows to recover best-fit DM halo profiles that are consistent with the SHMR. In this Thesis, the goal was not to constrain the shape of the DM halo profiles but rather (i) to recover the intrinsic circular velocity of the ionised gas and (ii) to derive the DM fraction, both within the optical radius of the galaxies. Because our MUSE data are restricted to the optical size of the galaxies and since we did not want to extrapolate our kinematics measurements at large distances where we would need additional data (e.g. HI), we did not use the SHMR in the modelling.
<figure><img src="image_29.png" /></figure>

Figure 5.2: Example of various dynamical fits performed on the same galaxy 104_CGR79. The top row shows the different best-fit DM rotation curves using the following profiles: NFW (solid line), Einasto (dashed line), pseudo isothermal sphere (dotted line), and Hubble modified (dashed-dotted line). The bulge and disk components used in the fit are also shown as solid blue and red lines, respectively. The vertical dashed blue line represents the position where the DM mass for each model was computed and the red solid line represents the extent of the MUSE velocity field, that is beyond which the rotation curves are extrapolated. The second row represents the best-fit velocity fields for each model and the last row shows the residuals of each fit. As a comparison, we also show the MUSE velocity field as an inset in the first row and the radius where the DM mass is computed as an ellipse in the second row.
## 5.3 Ionised gas kinematics with IFU data

In what follows, I will describe how I modelled the galaxies’ ionised gas kinematics for the [OII] emitters in the MAGIC survey. To begin with, I will discuss in Sect. 5.3.1 how I extracted the galaxies’ kinematics maps from their MUSE cubes using the line fitting tool CAMEL<fnref n="5" /> and then I will quickly go through the theory behind the modelling of the gas velocity field in Sect. 5.3.2. Finally, I will conclude in Sect. 5.3.3 by discussing how this methodology was used in Mercier et al. (2022) and in Chapter 6 in combination with mass models using the kinematics fitting tool MOCKING<fnref n="6" />.

### 5.3.1 Extraction of kinematics maps

The first part of the kinematics modelling is to extract the kinematics maps from the MUSE cubes. This step was done for the MAGIC survey during the first year of my PhD using the line fitting tool CAMEL. CAMEL is a python code developed by Epinat, B. (Epinat et al. 2009, 2012) that is dedicated to fitting emission lines from data cubes in order to derive kinematics maps such as velocity fields, velocity dispersion maps, line fluxes, and so on. The whole process works by writing one configuration file per set of emission lines that has to be fitted simultaneously. The configuration file contains all the information CAMEL requires to perform the fit. Several lines can be fitted simultaneously using the same constraints for the velocity and eventually the velocity dispersion, as well as additional constraints on the flux between a set of lines (e.g. Balmer lines, doublet flux ratio, etc.). Emission lines that are integrated by default in CAMEL are $H\alpha$, $H\beta$, $H\gamma$, $H\delta$, $[OII]$, $[OIII]$, $[NII]$, and $[SII]\lambda6716,6731$, but it is also possible to add any additional line using the EXTRAL keyword.

#### 5.3.1.1 Sub-cube extraction

To begin with, a sub-data cube centred on the galaxy’s location is extracted from the MUSE cube using the following information given in the configuration file: (i) the redshift of the galaxy, (ii) the emission lines to fit, (iii) the spatial dimensions of the sub-data cube, and (iv) the spectral dimension of the sub-data cube. The spatial dimensions of the sub-cube can be defined through the use of the four following keywords: XMIN, XMAX, YMIN, and YMAX. Similarly, the spectral range can be defined in the configuration file with the keywords ZMIN and ZMAX.
To improve the fitting process, it is possible to clip as well as spatially and/or spectrally smooth the cube. For the clipping, three values must be provided: (i) the $\sigma$ value of the $\sigma$-clipping, (ii) the number of passes, and (iii) the size of the box within which the median value is computed. The clipping is performed channel per channel using the usual $\sigma$-clipping method of replacing by the median value in the given channel all the spaxels that are $\sigma$ times above the median. The latter can either be computed in the entire channel or in a box of a given size. The process is repeated until the clipping condition is not met any more or if the maximum number of passes is reached. By default, a 10$\sigma$-clipping is performed with three passes and a box to compute the median value of 3 × 3 spaxels. On the other hand, for both spatial and spectral smoothing the user only has to provide the FWHM of the Gaussian window used by CAMEL. By default, it does not apply any spectral smoothing but does apply a sub-PSF (at least for a MUSE cube) spatial smoothing of 2 spaxels.

<footnote n="5">⁵https://gitlab.lam.fr/bepinat/CAMEL</footnote>

<footnote n="6">⁶https://gitlab.lam.fr/bepinat/MockinG</footnote>
<figure><img src="image_30.png" /><figcaption>Figure 5.3: Illustration of the extraction process of the kinematics maps from a MUSE cube using CAMEL for galaxy 104_CGR79. On the left-hand side is shown the sub-data and variance cubes extracted by CAMEL with three different spaxels that are highlighted. Their spectra (solid line) around the [O II] doublet with their single (orange dashed line) and total (green dotted line) best-fit Gaussian lines are represented in the middle of the figure. The variance spectrum used during the fit is also shown with a black dotted line. On the right-hand side is shown a few important kinematics maps extracted from the cube by CAMEL for all the spaxels: the velocity field, velocity dispersion map, and the two flux maps. The corresponding error maps are also represented on their right with the same scale.</figcaption></figure>

5.3.1.2 Emission line fitting

Once sub-data cubes are extracted, CAMEL will fit spaxel by spaxel each emission line using Gaussian profiles combined with a continuum modelled as a polynomial function. The degree of the polynomial can be provided by the user (zero by default) and CAMEL will automatically set the zeroth degree to the median value of the spaxel's spectrum. All the lines given in the configuration file will then be fitted at the same time by adjusting the redshift of the spaxel as well as the lines' amplitude and width. Furthermore, CAMEL uses by default an option that forces the lines to all have the same velocity dispersion. The user also has the possibility to provide additional keywords to further improve the fit such as giving lower and upper bounds on the dispersion of the line (WMIN, WMAX), the redshift range used to fit the lines (ZRMIN, ZRMAX) or the redshift range used to search for the best redshift (REDMIN, REDMAX).

After fitting each spaxel separately, taking into account the variance cube to weight the fits,
CAMEL will produce various kinematics maps, the most important ones being: (i) the integrated
lines' fluxes, (ii) their velocity offset with respect to the given systemic redshift, (iii) their velocity dispersion, (iv) their S/N, and (v) the uncertainty associated to any of these values, except for the S/N. Other kinematics maps that are less important for the kinematics modelling but that can still be useful are also produced such as the best-fit central wavelength of the line, its width in spaxels rather than in velocity, a model cube that only contains the best-fit Gaussian profiles, a continuum cube that only contains the best-fit continuum polynomial fit or a white light image. The process of fitting emission lines from the sub-data and variance cubes is illustrated in Fig. 5.3. Examples of spectral fits performed by CAMEL on the [OII] doublet is shown in the middle column for three spaxels with, from top to bottom, decreasing S/N. Each spectrum is represented as a solid black line, its uncertainty as a black dotted line, the two single line fits as orange dashed lines and the total fit of the doublet as a green dotted line. Doublet, such as [OII], are treated by CAMEL differently than singlets in the sense that the width of each line in the doublet is always forced to be equal to the that of the other one. The right-most two columns of Fig. 5.3 represent the raw kinematics maps (parameter and then uncertainty) that came out of the fitting process with on the first row the velocity field, on the second the velocity dispersion map, and on the third and last rows the flux of each line in the doublet.

### 5.3.1.3 Maps cleaning

<figure><img src="image_31.png" /><figcaption>Figure 5.4: Illustration of the cleaning process for the kinematics maps using the velocity field of 114_CGR79 as an example.</figcaption></figure>

A last but important step after the maps were extracted is to clean them. Indeed, a velocity field as presented in the right-most part of Fig. 5.3 will not be correctly fitted during the kinematics modelling process because of the contribution of all the pixels around that have large velocities. Most of these pixels tend to have a velocity which is equal to the lower or upper bound of the velocity offset allowed during the fit of the line. Their origin can be twofold: (i) they are noise dominated or (ii) there are sky residuals that were produced during the data reduction of the cube, in which case it was not the emission line that was fitted and therefore the pixel's velocity is not meaningful. This is more likely to happen when the fitted line falls in the reddest part of the spectrum beyond roughly 8000 Å (where sky lines are the most numerous and the brightest), that is at $z \gtrsim 1$ for [OII] emitters.
In Abril-Melgarejo et al. (2021) and Mercier et al. (2022) (see also Chapters 6 and 7), this cleaning process was done in two steps. First, we performed an automatic cleaning by keeping spaxels with a S/N larger than five and with a line width larger than 80% of the LSF FWHM computed using Eq. 3.1. The first criterion ensures that we are not keeping spaxels with a large noise contribution to the flux and the second criterion removes those that contain sky residuals. We used a conservative threshold of 80% for the latter criterion to ensure that we might not be removing spaxels with signal but with a low velocity dispersion given that there is an avoidable uncertainty on both the LSF and on our estimate of the line width that can both affect the dispersion value. Second, we visually inspected the velocity fields and the spectra for each remaining spaxel and removed (i) the few that were isolated from the bulk of the velocity field (typically 3-4 spaxels at most) and (ii) those that were connected to the velocity field but that had large velocity discontinuities with respect to their neighbours. These latter spaxels typically correspond to sky residuals that were not removed because of the conservative 80% threshold on the line width used during the automatic cleaning step that was discussed above. This step was done with PYQUBEVis<fnref n="7" />, a light-weight visualisation tool that allows to open at the same time a kinematics map and its corresponding cube to easily inspect the spectrum associated to each spaxel. As an illustration, we show the cleaning process for galaxy 104_CGR79 applied to its velocity field in Fig. 5.4.

## 5.3.2 Theory behind the kinematics modelling

Below we describe quickly the formalism that is used in MOCKING to model the velocity field extracted from the MUSE cubes. Because these derivations were not done during this Thesis and because the details are already given in full account in Epinat et al. (2008, 2010), we only provide as minimum information to describe the theory behind the kinematics modelling and we refer the interested reader to the appendices of these two papers for more information.

### 5.3.2.1 Line-of-sight velocity

The baseline for modelling a velocity field is to correctly model the LOS velocity in each spaxel. In this chapter, we have already discussed how one can go from a given mass distribution to derive a rotation velocity under the assumption of dynamical equilibrium. However, this rotation or circular velocity is not enough to model a velocity field. Indeed, the only component of the velocity that is accessible for intermediate to high redshift galaxies is the component projected along the line-of-sight: the LOS velocity. In general terms, for a kinematics tracer located in a plane (e.g. gaseous disk) and assuming the same cylindrical coordinates as in Sect. 5.1.1, the LOS velocity writes

$$V_{\text{LOS}} = V_s + V_{\theta} \cos \theta \sin i + V_R \sin \theta \sin i + V_z \cos i, \qquad (5.35)$$

where $V_s$ is the systemic velocity of the galaxy, $i$ is the inclination of the disk ($i = 0^{\circ}$ means face-on and $i = 90^{\circ}$ means edge-on) and $\theta$ is the azimuthal angle in the plane of the disk. The conversion from intrinsic galactic coordinates ($R, \theta$) to those on the plane of the sky ($R', \psi$) is given by (Epinat et al. 2008)

<footnote n="7">⁷https://gitlab.lam.fr/bepinat/PyQubeVis</footnote>
$$R \cos \theta = R' \cos \psi, \qquad (5.36)$$

$$R \sin \theta = R' \sin \psi / \cos i, \qquad (5.37)$$

$$R = R' \sqrt{\cos^2 \psi + \sin^2 \psi / \cos^2 i}, \qquad (5.38)$$

where $R' = \sqrt{x'^2 + y'^2}$ is the distance on the plane of the sky to the centre of the galaxy, with $x'$ and $y'$ the spaxel's position with respect to the galaxy centre ($y'$ positive to the north and $x'$ positive to the west), and $\psi$ is the counter-clockwise azimuthal angle with respect to the galaxy's major axis that measures the position of the spaxel on the plane of the sky and that writes

$$R' \cos \psi = y' \cos \text{PA} - x' \sin \text{PA}, \qquad (5.39)$$

$$R' \sin \psi = -x' \cos \text{PA} - y' \sin \text{PA}, \qquad (5.40)$$

with PA the position angle of the disk plane with respect to the North<fnref n="8" />. If there is only rotation, as is usually assumed, then Eq. 5.35 is simplified with $V_R = V_z = 0$ and we can use the above set of equations to compute the LOS velocity for any pixel given a rotation curve model. However, this is not sufficient in itself as there are strong instrument effects that impact the velocity fields and velocity dispersion maps of intermediate and high redshift galaxies.

### 5.3.2.2 Instrumental impact on the velocity field

The complete description of how the instrumental impact is taken into account when modelling the velocity field can be found in Appendix A of Epinat et al. (2010). However, since it is an important aspect of the dynamical modelling, I provide a full account below. To begin with, we note $S(x', y', \lambda)$ the total flux density at position $(x', y')$ on the plane of the sky and at an observed wavelength $\lambda$. It contains both the continuum and the line flux density noted $F_\lambda(x', y', \lambda)$. The flux of the line $F(x', y')$ is simply given by integrating $F_\lambda$ along the spectral dimension and the high-resolution velocity field corresponds to the first order moment of the line, that is

$$\bar{V}(x', y') = \frac{\int_{\lambda \in \mathbb{R}_+} d\lambda F_\lambda(x', y', \lambda)v(\lambda)}{F(x', y')}, \qquad (5.41)$$

where $v(\lambda)$ is the velocity offset that corresponds to the line being located at an observed wavelength $\lambda$. Equation 5.41 is therefore just saying that we measure the LOS velocity as the average velocity offset weighted by the flux of the line. An example of such a model is shown on the first row and leftmost column of Fig. 5.5 for galaxy 114_CGR79.

First, we consider the effects of the spectral resolution and binning. For the former, its effect will be to convolve the line flux density $F_\lambda$ with the LSF, while for the latter it will be equivalent to integrating the convolved flux density within a spectral window whose dimensions corresponds to those of the spexels. Hence, we have

$$S_2(x', y', \lambda_i) = F_{\lambda,2}(x', y', \lambda_i) + C_2(x', y', \lambda_i), \qquad (5.42)$$

<footnote n="8">⁸The notation has been slightly altered compared to Epinat et al. (2008) to be more consistent with previous chapters.</footnote>
<figure><img src="image_32.png" /><figcaption>Figure 5.5: Illustration of the instrumental impact on the kinematics maps for galaxy 114_CGR79. The first row shows the velocity field, the second row the velocity dispersion, and the last column the [O II] flux map at the original resolution and after interpolating with an oversampling factor of four. To illustrate the effect of the PSF on the shape of the velocity field, we represent four iso-velocity curves at -50, -100, 0, 50, and 100 km s⁻¹ on the high-resolution model (i.e. without PSF nor spatial binning, leftmost column), low-resolution model (after PSF and binning, second column), and on the velocity field extracted from the cube (third column).</figcaption></figure>

where $F_{\lambda,2}$ is the convolved plus spectrally binned line flux density and $C_2$ is the same for the continuum<fnref n="9" />. If we derive the velocity field by integrating over a small spectral window around the line, we might consider the continuum to be constant in which case it will still remain constant after the convolution and the binning. Since there is no flux loss during either steps, the line flux can be computed by summing up the contribution of $F_{\lambda,2}$ in each spectral channel $\lambda_i$. Thus, we can approximate the expression for the velocity field given in Eq. 5.41 by changing the integral into a sum as

$$ \bar{V}(x', y') \approx \sum_i \frac{F_{\lambda,2}(x', y', \lambda_i)v(\lambda_i)}{F(x', y')} \quad (5.43) $$

Now, we turn to the first of the two spatial impacts of the instruments: the PSF. Similarly to the LSF, the PSF will also convolve the flux density but spatially this time. We note this new spatially convolved line flux density $F_{\lambda,3}(x', y', \lambda_i)$. The integrated flux, computed by summing up the contribution of $F_{\lambda,3}$ in each channel $\lambda_i$ will also be impacted by the PSF and thus we note it $F_0(x', y')$. We can derive the expression for the velocity field when using $F_{\lambda,3}$ this time. To do so, we start again from Eq. 5.43 and we replace $F_{\lambda,2}$ by $F_{\lambda,3}$. Because we do not measure $F$ any more

<footnote n="9">⁹Following Epinat et al. (2010), $S_1$ corresponds to the spectrum only affected by the spectral convolution.</footnote>
but $F_0$, we also need to replace the normalisation by the new value. This writes

$$ 
\begin{aligned}
\overline{V_0}(x', y') &= \sum_i \frac{[F_{\lambda,2} \otimes \text{PSF}](x', y', \lambda_i)v(\lambda_i)}{F_0(x', y')} \\
&= \left[ \left( F(x', y') \sum_i \frac{F_{\lambda,2}(x', y', \lambda_i)v(\lambda_i)}{F(x', y')} \right) \otimes \text{PSF} \right] (x', y') / F_0(x', y') \\
&= \frac{[(\overline{V}F) \otimes \text{PSF}](x', y')}{F_0(x', y')},
\end{aligned}
\quad (5.44) $$

where $\otimes$ denotes the spatial convolution operation and $\overline{V}$ is defined in Eq. 5.41. The first line in Eq. 5.44 is just its expression when using $F_{\lambda,3}$ and $F_0$ instead of $F_{\lambda,2}$ and $F$. It is possible to rewrite it into the second line because $v(\lambda_i)$ is independent on both $x'$ and $y'$ and because $F$ is independent of $\lambda_i$. The last line is then just the second line rewritten using Eq. 5.43. Finally, there is one last instrument effect to take into account: spatial binning. It will be equivalent to integrating $F_{\lambda,3}$ within a spatial window corresponding to the dimensions of a spaxel. We note this new line flux density $F_{\lambda,4}(x'_j, y'_j, \lambda_i)$, where $(x'_j, y'_j)$ are now the discrete spaxels' coordinates. The integrated flux will also be impacted and we note it $F_1(x'_j, y'_j)$. Exactly in the same way as for the derivation of Eq. 5.44, we can compute the expression for this new velocity field when replacing $F_{\lambda,2}$ by $F_{\lambda,4}$ in the numerator of Eq. 5.43 and $F$ by $F_1$ in its denominator. The only difference with Eq. 5.44 is that there will be two additional integrals since the flux in each spaxel will be integrated over its surface because of the binning. The final expression writes

$$ \overline{V_1}(x'_j, y'_j) = \frac{\int_{x'_j - \Delta x}^{x'_j + \Delta x} dx' \int_{y'_j - \Delta y}^{y'_j + \Delta y} dy' [(\overline{V}F) \otimes \text{PSF}](x', y')}{F_1(x'_j, y'_j)}, \quad (5.45) $$

where the size of a spaxel along the $x'$ and $y'$ directions is equal to $2\Delta x$ and $2\Delta y$, respectively. Such a model is represented for 114_CGR79 on the first row and second column from the left in Fig. 5.5. It also corresponds to the best-fit model that was found by MOCKING. Thus we also show as a comparison the velocity field extracted from the MUSE cube with CAMEL on the third column.

### 5.3.2.3 Instrumental impact on the velocity dispersion map

The same procedure can be applied to derive the instrumental impact on the velocity dispersion map. Here, the point is not to model the velocity dispersion as for the velocity field but rather to correct it from instrumental effects. The velocity dispersion $\sigma(x', y')$ corresponds to the width of the line and therefore is defined as the second order moment of the flux density:

$$ 
\begin{aligned}
\sigma^2(x', y') &= \frac{\int_{\lambda \in \mathbb{R}_+} d\lambda F_{\lambda}(x', y', \lambda) [v(\lambda) - \overline{V}(x', y')]^2}{F(x', y', \lambda)} \\
&= \overline{V}^2(x', y') - \overline{V}^2(x', y'),
\end{aligned}
\quad (5.46) $$

We already know how $\overline{V}$ is impacted by the LSF and PSF convolutions and the spectral and spatial binning, so we only need to derive their impact on $\overline{V}^2$. As before, we can approximate its
value by considering the LSF convolved and spectrally binned spectrum and summing along the spectral channels instead of integrating:

$$ \overline{V^2}(x', y') \approx \frac{\sum_i F_{\lambda,2}(x', y', \lambda_i) v^2(\lambda_i)}{F(x', y')} \qquad (5.47) $$

The effect of the PSF will also write similarly by replacing $\overline{V}$ by $\overline{V^2}$ in Eq. 5.44 and so will the spatial binning. Thus, once all the instrumental effects are taken into account and using Eq. 5.46, the square of the velocity dispersion writes

$$ \sigma_1^2(x'_j, y'_j) = \frac{\int_{x'_j-\Delta x}^{x'_j+\Delta x} dx' \int_{y'_j-\Delta y}^{y'_j+\Delta y} dy' \left[ \left( \overline{V^2 F} \right) \otimes \text{PSF} \right] (x', y')}{F_1(x'_j, y'_j)} - \left[ \frac{\int_{x'_j-\Delta x}^{x'_j+\Delta x} dx' \int_{y'_j-\Delta y}^{y'_j+\Delta y} dy' \left[ \left( \overline{V F} \right) \otimes \text{PSF} \right] (x', y')}{F_1(x'_j, y'_j)} \right]^2 \qquad (5.48) $$

which reduces to

$$ \sigma_1^2(x'_j, y'_j) = \frac{\iint_{x'_j, y'_j} dx'dy' \left[ \sigma^2 F \otimes \text{PSF} \right] (x', y')}{F_1(x'_j, y'_j)} + \frac{\iint_{x'_j, y'_j} dx'dy' \left[ \overline{V^2} F \otimes \text{PSF} \right] (x', y')}{F_1(x'_j, y'_j)} - \left[ \frac{\iint_{x'_j, y'_j} dx'dy' \left[ \overline{V F} \otimes \text{PSF} \right] (x', y')}{F_1(x'_j, y'_j)} \right]^2 \qquad (5.49) $$

where the integrals are still to be understood within a spaxel. The first term in Eq. 5.49 corresponds to the velocity dispersion truly impacted by instrumental effects (LSF, PSF, and binning) whereas the other terms are purely instrumental effects that are referred to as beam smearing. It is caused by the velocity shear and the PSF that mixes signal from neighbouring spaxels with different velocities and therefore spuriously increases the velocity dispersion measured from the line width.

### 5.3.2.4 Brief comments on modelling the kinematics

We now have everything to model a velocity field taking into account the instrumental effects using Eq. 5.45. We also have an expression given by Eq. 5.49 that tells us how we can correct the velocity dispersion map extracted from a cube to recover the beam-smearing corrected map. Beside requiring a precise knowledge of the LSF and the PSF, each equation requires two important information: (i) a high-resolution velocity field $\overline{V}$ and (ii) the high-resolution flux map $F$. When modelling the velocity field with MOCKING the former is easy to get given that it corresponds to the model being generated. However, this is not true for the flux map that must be known a priori. If there is no high-resolution flux map available for the kinematics tracer being used, then MOCKING will compute an approximation by interpolating a low-resolution map. The interpolated flux map that was used during the modelling of 114_CGR79 is shown on the second
row and rightmost column of Fig. 5.5, with the original flux map extracted from the cube shown above. This interpolation was performed with an oversampling factor of four.

Once the best-fit high-resolution velocity field model has been found by *MOCKING*, it will compute the beam smearing component appearing in Eq. 5.49 and it will subtract it quadratically from the velocity dispersion map to recover the “true” velocity dispersion. The beam smearing correction is represented on the second row of Fig. 5.5 with the original velocity dispersion map on the leftmost column, the beam smearing model that was derived using the high resolution velocity field and the interpolated high-resolution flux map on the second column, and the beam smearing and LSF corrected dispersion map on the third column. We see that the beam smearing is higher where the flux and the velocity shear are large (by a factor of roughly 1.5 with respect to the outer parts of the map), which is near the centre. In addition, about half the spaxels have a null velocity dispersion after the correction was applied. The reason is that the correction is actually larger than the measured value from the line fit and therefore *MOCKING* clips it to zero. It is unlikely that the intrinsic gas dispersion is exactly zero in these pixels, though it is probably not very high (below roughly 30 km s⁻¹, e.g. see Appendix A of Boselli et al. 2021). What is more likely is that because we do not have access to the real high-resolution flux distribution, we are probably overestimating the flux in these pixels when using the interpolated map. Another possibility is due to the fact that we use the high-resolution velocity field model to derive the beam smearing but, as can be seen from the iso-density contours in Fig. 5.5, even the best-fit model does not perfectly fit the velocity field. Thus, the velocity shear is also an approximation of the underlying shear which will therefore also impact the derivation of the beam smearing. A last explanation for this apparent overcorrection could be that we are overestimating the value of the LSF since we have removed it quadratically from the velocity dispersion map.

Before presenting the application of the dynamical modelling to the MAGIC survey, I must first quickly discuss the efficiency of our method compared to other methods found in the literature. Indeed, our modelling relies on the moment map formalism (Epinat et al. 2008, 2010) which requires to extract kinematics maps from a data cube first and then fit them using the formalism. There exists another type of modelling that does not require the extraction of kinematics maps since it directly fits the data cube: 3D fitting codes such as Galaxy parameters and kinematics (GALPAK³D) described in Bouché et al. (2015) or 3D-Based Analysis of Rotating Object via Line Observation (³DBAROLO) described in Di Teodoro & Fraternali (2015). There are actually pros and cons for both methods and it is important to cite them. As discussed in Bouché et al. (2015) and as shown in subsequent analyses that used GALPAK³D (e.g. Bouché et al. 2021), 3D fitting codes are more efficient than our method to model the kinematics of low S/N and thus very low-mass galaxies. The main reason is that the moment map formalism requires kinematics maps which must be extracted by fitting spaxel-per-saxel one or more emission lines in the cube. Thus, for low S/N spaxels this will lead to uncertainties which will propagate to the kinematics model. Nevertheless, there are also some limitations to full 3D fitting codes. The main one is that such codes must produce a 3D galaxy model, hence a mock data cube, to fit the observed data. Therefore, this requires to produce (i) a kinematics model with underlying assumptions (e.g. rotation only, model of rotation curves, etc.), as is the case in our modelling, (ii) a model of the surface brightness distribution of the kinematics tracer that is fitted in the cube, and (iii) a model of the intrinsic velocity dispersion. For tilted ring models, which is what is used in ³DBAROLO, there no need to model the surface brightness distribution directly since it is a free parameter that is fitted in each ring. However, it does assume that the surface brightness is the same along the ring. On the contrary, our method does not need to make assumptions on the spatial distribution
of the kinematics tracer that is used nor on the velocity dispersion. The only assumption that we
make concerning the distribution of the ionised gas ([OII] doublet in our case) used as kinematics
tracer is that it is located in a disk, which is also an assumption made in other codes such as
GALPAK$^{3D}$. Furthermore, for both types of codes the net result is to recover the intrinsic
kinematics by taking into account the impact of beam smearing and to produce model kinematics
maps. Few analyses have tried to compare the results between 2D and 3D kinematics fitting codes
but those that have done so (e.g Contini et al. 2016) did find consistent results between the two
techniques. The only cases where such codes might yield different results are typically nearly
face-on and barely resolved galaxies where the large uncertainty on the inclination renders such
modellings difficult. In particular, our method has the advantage of producing a beam smearing
and LSF corrected velocity dispersion map without making assumptions on its spatial distribution.

5.3.3 Gas kinematics and dynamical modelling of the MAGIC sample

5.3.3.1 Contribution to MocKinG

At the beginning of the Thesis, the kinematics fitting code MOCKING was still in development. Instead, two prior codes written in Interactive Data Language (IDL) between 2008 and 2009 (Epinat et al. 2008, 2010) and which incorporated most of the features in the current code were available. The first focussed on local galaxies whereas the second was designed to model high-redshift galaxies by including the effect of beam smearing. These codes, especially the one specialised in modelling high-redshift galaxies, had been extensively tested and used in previous kinematics studies (e.g. Epinat et al. 2008, 2010, 2012; Vergani et al. 2012; Contini et al. 2016) and more recently in our previous analysis of the impact of the environment on a sub-sample of galaxies in MAGIC presented in Abril-Melgarejo et al. (2021)<fnref n="10" />. Meanwhile, there was an increasing need to rewrite the code in Python in a more modern and flexible manner that would be adapted to handle both local and high-redshift galaxies. A first version of this new code had been written by J. Dumoulin during his Master project at Laboratoire d'Astrophysique de Marseille (LAM) in 2016 and, in the following years, a few additional features were added by B. Epinat such as the possibility to have multiple rotation curve components for a single object. Therefore, one of the first aspects during this Thesis was to perform tests using galaxies from the MAGIC survey to compare the new results with the previous ones from the IDL code. I used the same rotation curve model (flat model, described in Sect. 5.3.3.3) as well as the same MUSE PSF and LSF FWHM values. MOCKING was available since the start with two different fitting techniques (see Sect. 5.3.3.3 for more details): a gradient descent-like algorithm and a Bayesian inference technique. On the contrary, the previous IDL code only worked with the former technique. Hence, the first test that was carried out was to compare the outputs of the two codes using the gradient descent-like algorithm. Then, the second test consisted in comparing the outputs of MOCKING between the gradient descent-like and Bayesian approaches for the same sample of galaxies in MAGIC. The result of the comparison is shown in Appendix A of Mercier et al. (2022) (see also Chapter 6). We found nearly identical results between MOCKING (in gradient descent mode) and the previous IDL code, as was expected given that the two codes fit the data in the same way. We also found consistent results (of the order of 10%) between the Bayesian and gradient descent-like approaches with MOCKING (as shown in Fig. A.1 of Mercier

<footnote n="10">¹⁰It is also this code that I used to perform a first kinematics analysis of the field galaxies in MAGIC during my Master project</footnote>
et al. 2022), with the Bayesian approach that seemed less prone to converge towards a local minimum solution in cases where the velocity field is perturbed. These tests allowed the code to be used for the first time in Boselli et al. (2021) and Adamczyk (2021).

In parallel, I started implementing mass models in `MOCKING`. In particular, the focus was put on properly integrating a rotation curve model for a double exponential disk as discussed in Sect. 5.2.4.3. After further tests, this new model was used for the first time in the analyses of the galaxy scaling relations discussed in Mercier et al. (2022) and in Chapter 6 and of the stellar angular momentum in MAGIC presented in Chapter 7.

### 5.3.3.2 Kinematics maps extraction

In Mercier et al. (2022) and in Chapter 6 we used the kinematics fitting tool `MOCKING` that implements the moment maps formalism presented in Sect. 5.3.2 to model the dynamics of the galaxies in the MAGIC survey. This required first to extract the kinematics maps using the [OII] doublet with `CAMEL` (see Sect. 5.3.1) for the 890 galaxies that are part of the morphological sample, that is [OII] emitters in the MAGIC sample for which it was possible to model the morphology with `GALFIT`. This was done in an automatic fashion with `CAMEL` by extracting sub-data cubes of dimensions 30 × 30 MUSE spaxels around the observed wavelength of the [OII] doublet, which was then spatially smoothed with a Gaussian profile with a sub-resolution FWHM of 2 spaxels. The continuum was fitted with a constant value and the doublet with two Gaussian profiles with the same velocity offset and line width. An example of configuration file used by `CAMEL` to extract the maps for the galaxy 114_CGR79 is shown below:

<table><tr><td>1</td><td>FITSFILE= CGr79_d.fits</td><td>/ Path of the cube FITS file</td></tr><tr><td>2</td><td>OUTPUT = CGr79_d_114_o2</td><td>/ Output directory</td></tr><tr><td>3</td><td>SKYFILE = CGr79_d.fits</td><td>/ Sky spectrum cube FITS file</td></tr><tr><td>4</td><td>HALPHA = False</td><td>/ Fit Halpha ? (default:TRUE)</td></tr><tr><td>5</td><td>OII = True</td><td>/ Fit OII ? (default:FALSE)</td></tr><tr><td>6</td><td>COMMW = True</td><td>/ Use a common width ? (default:FALSE)</td></tr><tr><td>7</td><td>REDSHIFT= 0.6684</td><td>/ Initial redshift</td></tr><tr><td>8</td><td>REDMIN = 0.665617408317857</td><td>/ Minimum redshift for line fit</td></tr><tr><td>9</td><td>REDMAX = 0.671182591682143</td><td>/ Maximum redshift for line fit</td></tr><tr><td>10</td><td>ZRMIN = 0.6572696332714281</td><td>/ Minimum redshift for spectral range around lines</td></tr><tr><td>11</td><td>ZRMAX = 0.6795303667285719</td><td>/ Maximum redshift for spectral range around lines</td></tr><tr><td>12</td><td>INITW = 50.0</td><td>/ Initial line width in km/s</td></tr><tr><td>13</td><td>WMIN = 30.0</td><td>/ Minimum line width in km/s</td></tr><tr><td>14</td><td>WMAX = 250.0</td><td>/ Maximum line width in km/s</td></tr><tr><td>15</td><td>DFIT = 100.0</td><td>/ Bin to fit the line center in km/s</td></tr><tr><td>16</td><td>DGCTNUM = 0</td><td>/ Polynomial degree for the continuum fit</td></tr><tr><td>17</td><td>SCLIP = 10</td><td>/ Sigma-clipping threshold to clean the cube</td></tr><tr><td>18</td><td>XYCLIP = 3</td><td>/ Box width for the median for the sigma-clipping</td></tr><tr><td>19</td><td>NCLIP = 3</td><td>/ Number of passes for clipping (default:3)</td></tr><tr><td>20</td><td>WSMOOTH = 0</td><td>/ Spectral smoothing</td></tr><tr><td>21</td><td>SSMOOTH = 2.0</td><td>/ Spatial smoothing FWHM in pixel (default:0)</td></tr><tr><td>22</td><td>XMIN = 197.0</td><td>/ X minimum value for spatial cut of the input cube</td></tr><tr><td>23</td><td>XMAX = 227.0</td><td>/ X maximum value for spatial cut of the input cube</td></tr><tr><td>24</td><td>YMIN = 251.0</td><td>/ Y minimum value for spatial cut of the input cube</td></tr><tr><td>25</td><td>YMAX = 281.0</td><td>/ Y maximum value for spatial cut of the input cube</td></tr><tr><td>26</td><td>ZMIN = 1180.0</td><td>/ Z minimum value for spectral cut of the input cube</td></tr><tr><td>27</td><td>ZMAX = 1250.0</td><td>/ Z maximum value for spectral cut of the input cube</td></tr></table>

Code 5.1: Example of a `CAMEL` configuration file for galaxy 114_CGR79.
After the maps were extracted, they were cleaned as described in Sect. 5.3.1.3. Since not every galaxy in the [OII] emitters redshift range necessarily has enough [OII] signal to be detected, the cleaning phase led to the removal of 297 galaxies from the morphological sample, thus producing a kinematics sample of 593 galaxies. Because the cleaning of spurious spaxels is mostly driven by S/N considerations, as discussed in Mercier et al. (2022), it is roughly similar to removing galaxies whose integrated [OII] flux is equal to $2 \times 10^{-18} \text{ erg s}^{-1} \text{ cm}^{-2}$ at the median redshift of the groups in MAGIC ($z \approx 0.7$).

### 5.3.3.3 Dynamical modelling

Finally, the last step was to perform the kinematics modelling with MOCKING. First, in Abril-Melgarejo et al. (2021), we used a simple flat model for the rotation curve that writes

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

where $V_t$ is the plateau velocity and $R_t$ is the transition radius between the linearly rising part and the outer plateau. Even though this model is not physically motivated, it is robust to model and to recover the galaxies’ rotation velocity at intermediate redshift (e.g. Wright et al. 2007; Epinat et al. 2009, 2010, 2012; Contini et al. 2016). The ionised gas was assumed to be located in a rotating razor-thin disk. To remove any degeneracies between the velocity field centre’s location and the systemic velocity of the galaxy, as well as between $V_t$ and the gas disk’s inclination, we fixed both the centre’s position and the inclination to the values derived with GALFIT during the morphological modelling. Furthermore, we also fixed the kinematics parameters of the stellar disk and bulge components when performing a mass modelling by assuming that they are robustly constrained from the morphological modelling. All the other parameters (rotation curve parameters of the DM halo component, kinematics position angle (PA), and systemic velocity) were let free during the fitting process. MOCKING implements two methods to perform the fit: (i) CAT_MPFIT<fnref n="11" />, a python implementation of MPFIT (Markwardt 2009) that is a gradient descent kind of technique based on the Levenberg-Marquardt algorithm and (ii) a python implementation of MULTINEST (Feroz & Hobson 2008; Buchner et al. 2014) that is a Bayesian inference tool based on a multi-nodal nested sampling algorithm. As discussed in Sect. 5.3.3.1, we originally used the former algorithm but, after comparing the results between the two methods in Mercier et al. (2022), we relied on the latter in subsequent modellings as it is less prone to finding a local rather than a global minimum.

A first modelling was performed in Mercier et al. (2022) using the same assumptions as in Abril-Melgarejo et al. (2021) but implementing a mass modelling approach (see Sect. 5.2.1). We used a double exponential disk model for the stellar disk component as described in Sect. 5.2.4.3, a Hernquist’s model (see Sect. 5.2.4.1) for the stellar bulge, and we also included an unconstrained NFW DM halo (see Sect. 5.2.4.4). By definition the only constraint on the DM halo is the kinematics itself so that we let free its maximum rotation velocity $V_{h,max}$ and its scale radius $R_s$ during the kinematics modelling. In practice, we should have also included the contribution of the cold gas components (in particular the molecular gas in the inner regions) that have a non-negligible impact, especially at intermediate and high redshift (e.g. Tacconi et al. 2020). To do so, we would have needed to have an idea of its spatial distribution and total mass either from

<footnote n="11">¹¹https://www-astro.physics.ox.ac.uk/~cappellari/software/</footnote>
additional observations (e.g. ALMA) or from scaling laws (e.g. Kennicutt–Schmidt relation). Afterwards, we could have included it as an additional component to the kinematics model (e.g. razor-thin disk). Another option would have been to assume a spatial distribution for the cold gas (e.g. the same exponential distribution as that of the stellar disk or a flat distribution as in Bouché et al. 2022) and then to let the amplitude of this component free to vary (i.e. add the M/L of the gas component as an additional free parameter). Because we had no constraints on the gas mass and distribution in the galaxies of the MAGIC survey, we have decided in Mercier et al. (2022) to not take explicitly into account the contribution of the cold gas to the ionised gas kinematics. In practice, this means that it will be implicitly included in the DM halo component during the fitting process and that it must be accounted for during the analysis when estimating the mass of the NFW profile if one wants to estimate the mass of the DM component alone (e.g. see Sect. 7.4 of Mercier et al. 2022).

Examples of such models can be found in Sect. 5 of Mercier et al. (2022) and in Chapter 6¹². Finally, the same modelling was performed once more for the analysis of the angular momentum (see Chapter 7). Indeed, after the derivation of new stellar mass estimates with CIGALE, as well as more precise MUSE PSF profiles using Moffat profiles (both for the survey paper, see Epinat et al., in prep.), it was decided to perform again the kinematics modelling using these latest values in order to have the best constraints on the galaxies’ kinematics parameters. Besides, the morphology was updated for 17 galaxies in the kinematics sample for the analysis of the angular momentum, mainly to derive better constraints on their bulge component. Since the morphological models directly affect the mass modelling, these galaxies also required to update their kinematics models. An example of a MOCKING configuration file for galaxy 114_CGR79 is shown below. This file corresponds to the latest modelling, including the mass models and the Moffat profile for the PSF.

<footnote n="12">¹²The full morpho-dynamical models can also be found as supplementary material of Mercier et al. (2022) at the following address: https://www.aanda.org/articles/aa/olm/2022/09/aa43110-22/aa43110-22.html</footnote>
1 config fit:
2 method: multinest # multinest or mpfit
3 multinest:
4 evidence tolerance: 0.5
5 max_iter: 0
6 max modes: 100
7 mode tolerance: -1.0e+60
8 n_iter before update: 100
9 n live points: 50
10 null log evidence: -1.0e+90
11 pltstats: false
12 sampling efficiency: 0.8
13 verbose: false
14 files:
15 disp: CGr79_d_114_o2_ssmooth disp common mclean5.0. fits # Dispersion map
16 errvel: CGr79_d_114_o2_ssmooth evel common mclean5.0. fits # Error on Vfield
17 psf: null # PSF image (null=Gaussian/Moffat)
18 vel: CGr79_d_114_o2_ssmooth vel common mclean5.0. fits # Vfield
19 name: CGr79_d_114 # ID
20 objects:
21 obj1: # Only one object implemented for now
22 components: # Free number of components
23 comp1:
24 model: bovy corr # Name of the model in MocKinG - here Bovy curve
25 params:
26 q0: # Disk thickness required for a corrected Bovy curve
27 desc: q0 # Parameter name in the output file
28 fixed: 1 # 1 for fixed, 0 for free
29 limits:
30 - 0.17890397942393682 # Lower bound
31 - 0.17890397942393682 # Upper bound
32 value: 0.17890397942393682 # Initial value
33 rd: # Disk scale length
34 desc: rd
35 fixed: 1
36 limits:
37 - 2.189654476613038
38 - 2.189654476613038
39 value: 2.189654476613038
40 vm_RT: # Razor-thin maximum velocity
41 desc: vm_RT
42 fixed: 1
43 limits:
44 - 165.89090254236532
45 - 165.89090254236532
46 value: 165.89090254236532
47 vm corr: # Maximum velocity of the correction
48 desc: vm corr
49 fixed: 1
50 limits:
51 - 68.40760325005385
52 - 68.40760325005385
53 value: 68.40760325005385
54 comp2:
55 model: hernquist # Hernquist bulge model
56 params:
57 rt: # Radius where vm is reached
58 desc: a
59 fixed: 1
60 limits:
61 - 0.1699294100661157
62 - 0.1699294100661157
63 value: 0.1699294100661157
64 vm: # Maximum rotation velocity
65 desc: vm hern
66 fixed: 1
67 limits:
68 - 121.70174860901213
69 - 121.70174860901213