# Constraining the stellar populations of ultra-diffuse galaxies in the MATLAS survey using spectral energy distribution fitting

Maria Luisa Buzzo, Duncan A. Forbes, Thomas H. Jarrett, Francine R.

Marleau, Pierre-Alain Duc, Jean P. Brodie, Aaron J. Romanowsky, Jonah S.

Gannon, Steven R. Janssens, Joel Pfeffer, et al.

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

► To cite this version:

Maria Luisa Buzzo, Duncan A. Forbes, Thomas H. Jarrett, Francine R. Marleau, Pierre-Alain Duc, et al. Constraining the stellar populations of ultra-diffuse galaxies in the MATLAS survey using spectral energy distribution fitting. Monthly Notices of the Royal Astronomical Society, 2024, 529, pp.3210-3234. 10.1093/mnras/stae564. insu-04538495

HAL Id: insu-04538495

https://insu.hal.science/insu-04538495

Submitted on 11 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining the stellar populations of ultra-diffuse galaxies in the MATLAS survey using spectral energy distribution fitting

Maria Luisa Buzzo, Duncan A. Forbes, Thomas H. Jarrett, Francine R.

Marleau, Pierre-Alain Duc, Jean P. Brodie, Aaron J. Romanowsky, Jonah S.

Gannon, Steven R. Janssens, Joel Pfeffer, et al.

► To cite this version:

Maria Luisa Buzzo, Duncan A. Forbes, Thomas H. Jarrett, Francine R. Marleau, Pierre-Alain Duc, et al.. Constraining the stellar populations of ultra-diffuse galaxies in the MATLAS survey using spectral energy distribution fitting. Monthly Notices of the Royal Astronomical Society, 2024, 529, pp.3210-3234. 10.1093/mnras/stae564. insu-04538495

HAL Id: insu-04538495

https://insu.hal.science/insu-04538495v1

Submitted on 11 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Constraining the stellar populations of ultra-diffuse galaxies in the MATLAS survey using spectral energy distribution fitting

Maria Luisa Buzzo<sup>D, 1,2*</sup> Duncan A. Forbes,<sup>1,2</sup> Thomas H. Jarrett<sup>D, 3</sup> Francine R. Marleau,<sup>4</sup>

Pierre-Alain Duc,<sup>5</sup> Jean P. Brodie,<sup>1,2,6</sup> Aaron J. Romanowsky<sup>D, 7,8</sup> Jonah S. Gannon,<sup>1,2</sup>

Steven R. Janssens<sup>D, 1,2</sup> Joel Pfeffer<sup>D, 1,2</sup> Anna Ferré-Mateu<sup>D, 1,9,10</sup> Lydia Haacke,<sup>1,2</sup>

Warrick J. Couch,<sup>1</sup> Sungsoon Lim<sup>11</sup> and Rubén Sánchez-Janssen<sup>12</sup>

Accepted 2024 February 16. Received 2024 January 19; in original form 2023 October 10

## ABSTRACT

We use spectral energy distribution fitting to place constraints on the stellar populations of 59 ultra-diffuse galaxies (UDGs) in the low-to-moderate density fields of the MATLAS survey. We use the routine PROSPECTOR, coupled with archival data in the optical from the Dark Energy Camera Legacy Survey, and near- and mid-infrared imaging from the Wide-field Infrared Survey Explorer, to recover the stellar masses, ages, metallicities, and star formation time-scales of the UDGs. We find that a subsample of the UDGs lies within the scatter of the mass-metallicity relation (MZR) for local classical dwarfs. However, another subsample is more metal-poor, being consistent with the evolving MZR at high redshift. We investigate UDG positioning trends in the mass-metallicity plane as a function of surface brightness, effective radius, axis ratio, local volume density, mass-weighted age, star formation time-scale, globular cluster (GC) counts, and GC specific frequency. We find that our sample of UDGs can be separated into two main classes: Class A: comprised of UDGs with lower stellar masses, prolonged star formation histories (SFHs), more elongated, inhabiting less dense environments, hosting fewer GCs, younger, consistent with the classical dwarf MZR, and fainter. Class B: UDGs with higher stellar masses, rapid SFHs, rounder, inhabiting the densest of our probed environments, hosting on average the most numerous GC systems, older, consistent with the high-redshift MZR (i.e. consistent with early-quenching), and brighter. The combination of these properties suggests that UDGs of Class A are consistent with a ‘puffed-up dwarf’ formation scenario, while UDGs of Class B seem to be better explained by ‘failed galaxy’ scenarios.

**Key words:** galaxies: formation – galaxies: fundamental parameters – galaxies: star clusters: general – galaxies: stellar content.

## 1 INTRODUCTION

Ultra-diffuse galaxies (UDGs), despite having been heavily studied for almost a decade now (van Dokkum et al. 2015), are still the topic of much debate regarding their formation mechanisms, including dark matter content (e.g. Toloba et al. 2018, 2023; Mancera Piña et al. 2019a, b, 2022; Forbes et al. 2020; Gannon et al. 2020, 2021, 2022, 2023a; Iodice et al. 2023; Zaritsky et al. 2023), stellar populations (e.g. Román & Trujillo 2017; Ferré-Mateu et al. 2018, 2023; Ruiz-

Lara et al. 2018; Chilingarian et al. 2019; Barbosa et al. 2020; Buzzo et al. 2022; Villaume et al. 2022; Heesters et al. 2023; Iodice et al. 2023), and their globular cluster (GC) systems (e.g. van Dokkum et al. 2018; Trujillo et al. 2019; Forbes et al. 2020; Müller et al. 2020; Gannon et al. 2021, 2022, 2023a; Danieli et al. 2022).

While resembling classical dwarf galaxies in terms of their stellar masses ($M_{\star} < 10^{9} M_{\odot}$) and surface brightnesses ($\mu_{g,0} > 24$ mag arcsec$^{-2}$), they by definition exhibit much larger effective radii ($R_{e} > 1.5$ kpc). Some of them were found to have populous GC systems and evidence for more massive dark matter haloes ($M_{\mathrm{halo}} > 10^{11} M_{\odot}$) than their classical dwarf counterparts (Beasley et al. 2016; van Dokkum et al. 2019; Forbes et al. 2020; Gannon et al. 2020, 2022,

<footnote n="1"><sup>1</sup>Centre for Astrophysics and Supercomputing, Swinburne University, John Street, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="2"><sup>2</sup>ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote n="3"><sup>3</sup>Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Road, Cape Town 7535, South Africa</footnote>

<footnote n="4"><sup>4</sup>Institut für Astro- und Teilchenphysik, Universität Innsbruck, Technikerstraße 25/8, A-6020 Innsbruck, Austria</footnote>

<footnote n="5"><sup>5</sup>CNRS, Observatoire astronomique de Strasbourg Université de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="6"><sup>6</sup>University of California Observatories, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, San José State University, One Washington Square, San Jose, CA 95192, USA</footnote>

<footnote n="8"><sup>8</sup>Department of Astronomy & Astrophysics, University of California Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="9"><sup>9</sup>Instituto Astrofísica de Canarias, Av. Via Lactea s/n, E-38205 La Laguna, Spain</footnote>

<footnote n="10"><sup>10</sup>Departamento de Astrofísica, Universidad de La Laguna, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="11"><sup>11</sup>Department of Astronomy, Yonsei University, 50 Yonsei-ro Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="12"><sup>12</sup>UK Astronomy Technology Centre, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote>*E-mail: luisa.buzzo@gmail.com</footnote>
# Constraining the stellar populations of ultra-diffuse galaxies in the MATLAS survey using spectral energy distribution fitting

Maria Luisa Buzzo<sup>D, 1,2*</sup> Duncan A. Forbes,<sup>1,2</sup> Thomas H. Jarrett<sup>D, 3</sup> Francine R. Marleau,<sup>4</sup>

Pierre-Alain Duc,<sup>5</sup> Jean P. Brodie,<sup>1,2,6</sup> Aaron J. Romanowsky<sup>D, 7,8</sup> Jonah S. Gannon,<sup>1,2</sup>

Steven R. Janssens<sup>D, 1,2</sup> Joel Pfeffer<sup>D, 1,2</sup> Anna Ferré-Mateu<sup>D, 1,9,10</sup> Lydia Haacke,<sup>1,2</sup>

Warrick J. Couch,<sup>1</sup> Sungsoon Lim<sup>11</sup> and Rubén Sánchez-Janssen<sup>12</sup>

Accepted 2024 February 16. Received 2024 January 19; in original form 2023 October 10

## ABSTRACT

We use spectral energy distribution fitting to place constraints on the stellar populations of 59 ultra-diffuse galaxies (UDGs) in the low-to-moderate density fields of the MATLAS survey. We use the routine PROSPECTOR, coupled with archival data in the optical from the Dark Energy Camera Legacy Survey, and near- and mid-infrared imaging from the Wide-field Infrared Survey Explorer, to recover the stellar masses, ages, metallicities, and star formation time-scales of the UDGs. We find that a subsample of the UDGs lies within the scatter of the mass-metallicity relation (MZR) for local classical dwarfs. However, another subsample is more metal-poor, being consistent with the evolving MZR at high redshift. We investigate UDG positioning trends in the mass-metallicity plane as a function of surface brightness, effective radius, axis ratio, local volume density, mass-weighted age, star formation time-scale, globular cluster (GC) counts, and GC specific frequency. We find that our sample of UDGs can be separated into two main classes: Class A: comprised of UDGs with lower stellar masses, prolonged star formation histories (SFHs), more elongated, inhabiting less dense environments, hosting fewer GCs, younger, consistent with the classical dwarf MZR, and fainter. Class B: UDGs with higher stellar masses, rapid SFHs, rounder, inhabiting the densest of our probed environments, hosting on average the most numerous GC systems, older, consistent with the high-redshift MZR (i.e. consistent with early-quenching), and brighter. The combination of these properties suggests that UDGs of Class A are consistent with a ‘puffed-up dwarf’ formation scenario, while UDGs of Class B seem to be better explained by ‘failed galaxy’ scenarios.

**Key words:** galaxies: formation – galaxies: fundamental parameters – galaxies: star clusters: general – galaxies: stellar content.

## 1 INTRODUCTION

Ultra-diffuse galaxies (UDGs), despite having been heavily studied for almost a decade now (van Dokkum et al. 2015), are still the topic of much debate regarding their formation mechanisms, including dark matter content (e.g. Toloba et al. 2018, 2023; Mancera Piña et al. 2019a, b, 2022; Forbes et al. 2020; Gannon et al. 2020, 2021, 2022, 2023a; Iodice et al. 2023; Zaritsky et al. 2023), stellar populations (e.g. Román & Trujillo 2017; Ferré-Mateu et al. 2018, 2023; Ruiz-

Lara et al. 2018; Chilingarian et al. 2019; Barbosa et al. 2020; Buzzo et al. 2022; Villaume et al. 2022; Heesters et al. 2023; Iodice et al. 2023), and their globular cluster (GC) systems (e.g. van Dokkum et al. 2018; Trujillo et al. 2019; Forbes et al. 2020; Müller et al. 2020; Gannon et al. 2021, 2022, 2023a; Danieli et al. 2022).

While resembling classical dwarf galaxies in terms of their stellar masses ($M_{\star} < 10^{9} M_{\odot}$) and surface brightnesses ($\mu_{g,0} > 24$ mag arcsec$^{-2}$), they by definition exhibit much larger effective radii ($R_{e} > 1.5$ kpc). Some of them were found to have populous GC systems and evidence for more massive dark matter haloes ($M_{\mathrm{halo}} > 10^{11} M_{\odot}$) than their classical dwarf counterparts (Beasley et al. 2016; van Dokkum et al. 2019; Forbes et al. 2020; Gannon et al. 2020, 2022,

<footnote n="1"><sup>1</sup>Centre for Astrophysics and Supercomputing, Swinburne University, John Street, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="2"><sup>2</sup>ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Australia</footnote>

<footnote n="3"><sup>3</sup>Department of Physics and Astronomy, University of the Western Cape, Robert Sobukwe Road, Cape Town 7535, South Africa</footnote>

<footnote n="4"><sup>4</sup>Institut für Astro- und Teilchenphysik, Universität Innsbruck, Technikerstraße 25/8, A-6020 Innsbruck, Austria</footnote>

<footnote n="5"><sup>5</sup>CNRS, Observatoire astronomique de Strasbourg Université de Strasbourg, UMR 7550, F-67000 Strasbourg, France</footnote>

<footnote n="6"><sup>6</sup>University of California Observatories, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="7"><sup>7</sup>Department of Physics and Astronomy, San José State University, One Washington Square, San Jose, CA 95192, USA</footnote>

<footnote n="8"><sup>8</sup>Department of Astronomy & Astrophysics, University of California Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA</footnote>

<footnote n="9"><sup>9</sup>Instituto Astrofísica de Canarias, Av. Vía Lactea s/n, E-38205 La Laguna, Spain</footnote>

<footnote n="10"><sup>10</sup>Departamento de Astrofísica, Universidad de La Laguna, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="11"><sup>11</sup>Department of Astronomy, Yonsei University, 50 Yonsei-ro Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="12"><sup>12</sup>UK Astronomy Technology Centre, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote>*E-mail: luisa.buzzo@gmail.com</footnote>
2023a; Toloba et al. 2023; Zaritsky et al. 2023). This mixture of dwarf- and massive-like galaxy properties has led to many different proposed formation scenarios for these galaxies.

A scenario relying on internal processes suggests that UDGs start their lives as classical dwarf galaxies, and through an extended sequence (i.e. repeated episodes over a long period of time) of supernovae feedback (Di Cintio et al. 2017; Chan et al. 2018), they get enlarged to the sizes we observe today. Another scenario suggests that dwarfs with high-spin haloes can evolve into UDGs (Amorisco & Loeb 2016; Benavides et al. 2021). We refer to such scenarios as ‘puffed-up dwarf’ formation scenarios. Other scenarios rely on external processes to form UDGs via e.g. tidal stripping and heating (e.g. Carleton et al. 2019) and tidal interaction (e.g. Tremmel et al. 2020), which account for the apparent radial size excess. Some UDGs in the past have been suggested to be consistent with a tidal dwarf galaxy origin as well (see e.g. Duc et al. 2014; Buzzo et al. 2023; Gannon et al. 2023b). We refer to this formation scenario throughout the paper as ‘tidal UDGs’.

In addition to these, one scenario suggests that UDGs started their lifetimes destined to be large and massive but had their star formation truncated early on. As a consequence, their stellar masses do not increase at the same rate as their haloes, resulting in galaxies with dwarf-like stellar masses enclosed by halo masses similar to those of more massive M33-like galaxies (van Dokkum et al. 2015; Peng & Lim 2016; Danieli et al. 2022; Gannon et al. 2023a). For the remainder of the manuscript, we refer to this formation scenario as ‘failed galaxy’ scenario.

One crucial discriminant between these formation scenarios is the stellar populations of UDGs. For example, if they were formed by a puffed-up dwarf scenario, they would overall be expected to have similar stellar populations to those of classical dwarfs. One expectation, for example, is that they may follow the same scaling relations as classical dwarfs, e.g. the mass-metallicity relation (MZR). It is important to mention, however, that different puffed-up dwarf scenarios may lead to very different star formation histories (SFHs), i.e. supernovae feedback would lead to a bursty SFH, while high-spin haloes or tidal interactions may be more similar to an exponentially decaying SFH, but in both cases the UDGs are expected to preserve the stellar populations of their progenitors. Tidal UDGs are expected to be young, metal-rich, and gas-dominated for their stellar mass (Duc et al. 2014; Haslbauer et al. 2019). On the other hand, if UDGs are failed galaxies, the predictions for their stellar populations are less clear. They are not expected to have stellar populations similar to classical dwarfs (as classical dwarfs were not their progenitors), nor to more massive M33-like galaxies (as they did not evolve to become such). The only stellar population prediction for galaxies formed by this scenario is that they would have old or even ancient stellar populations (≥10 Gyr), as they by definition have suffered from early-quenching. For the same reason, they may be expected to have lower metallicities and possibly be alpha enhanced as they have had shorter time-scales to form stellar masses comparable to those of classical dwarfs that had prolonged SFHs (see e.g. Forbes et al. 2020), and also because, as suggested by Danieli et al. (2022), these galaxies could be mainly formed of disrupted GCs, thus being more metal-poor. As a key point in the comparison, the MZR is known to evolve with redshift (Ma et al. 2016; Chartab et al. 2024). Galaxies that follow the MZR at high redshift have stellar populations consistent with the chemical enrichment up to that point, thus, either galaxies at high redshift themselves or galaxies that have not had much chemical evolution since then (i.e. early quenched). We thus may expect failed galaxies to follow the MZRs of high-redshift galaxies.

Some of these predictions were recently confirmed by Ferré-Mateu et al. (2023) (hereafter, FM23), who used the largest spectroscopic study of UDGs to date (25) to show that there is a correlation between α-enrichment and the SFH of UDGs, further emphasized by the positioning of UDGs in the cluster phase-space diagram. They found that early-infall UDGs are the ones with the fastest star formation episodes and have higher [α/Fe], consistent with early-quenching scenarios, such as failed galaxy ones. Conversely, UDGs with prolonged SFHs are the late infallers, having lower [α/Fe] ratios, being better explained by puffed-up dwarf scenarios.

Additionally, Buzzo et al. (2022) (hereafter, B22) used spectral energy distribution (SED) fitting to show that there is a correlation between the age, metallicity, and environment that UDGs reside in. The older ones display the lowest metallicities, are consistent with the MZR at high redshift (z ~ 2), and are found in the densest environments. The younger ones (≤8 Gyr) follow the classical dwarf MZR, being more metal-rich and are found in less dense environments, such as the field and groups. These findings are in general agreement with the findings of Barbosa et al. (2020), who used SED fitting to study 100 UDGs in the field and to show that they are mostly young (~6 Gyr) and follow the classical dwarf MZR. In their sample, some UDGs were found to be more metal-rich than the classical dwarf MZR ([Fe/H] ≥ −0.5 dex) and extremely young (≤1 Gyr), similar to what has been found for some cluster UDGs (Ferré-Mateu et al. 2018; Ruiz-Lara et al. 2018). Their ages and metallicities are consistent with them being tidal UDGs (Collins & Read 2022).

B22 have also shown that the stellar populations of UDGs seem to correlate with their GC-richness. Old and metal-poor UDGs have on average the highest number of GCs, while younger and more metal-rich ones (consistent with the classical dwarf MZR) have the lowest GC numbers. This trend, however, was not found by FM23 using spectroscopy. Selection effects can be the cause of the differences, since both samples are not fully representative of the population of UDGs. Thus, to understand if this trend holds and whether it can or cannot be extrapolated to other UDGs, it is necessary to test it on a larger sample of UDGs. These would preferably be spread across different environments and span a variety of GC numbers so that selection effects can be diminished and conclusive interpretations made.

In this study, we extend the work of B22 to a sample of 59 UDGs in the MATLAS low- and moderate-density environments (Marleau et al. 2021), to help balancing the sample that was biased towards higher density environments in B22, and to thus start building up a representative sample of UDGs. 38 out of these 59 MATLAS UDGs have GC counts from *Hubble Space Telescope (HST)* imaging. Similar to B22, we employ optical to mid-infrared SED fitting to explore the stellar populations of these MATLAS UDGs. Buzzo et al. (in preparation) will combine both samples (i.e. B22 and this study) with a control sample of classical dwarf galaxies to perform an statistically meaningful and comprehensive study of UDGs spread across environments.

The paper is structured as follows: in Section 2 we present a summary of the UDG sample studied in this work and the data available for each UDG. In Section 3 we describe our methods to obtain the photometry and morphology of the galaxies and our SED fitting methodology. In Section 4 we provide our results. In Section 5 we discuss the implications of our results within the theoretical predictions for UDGs and as compared to the literature. In Section 6 we present the summary and the conclusions of the paper.

Throughout this paper, when converting distances to redshifts (or vice versa), we assume the cosmological parameters from the
2023a; Toloba et al. 2023; Zaritsky et al. 2023). This mixture of dwarf- and massive-like galaxy properties has led to many different proposed formation scenarios for these galaxies.

A scenario relying on internal processes suggests that UDGs start their lives as classical dwarf galaxies, and through an extended sequence (i.e. repeated episodes over a long period of time) of supernovae feedback (Di Cintio et al. 2017; Chan et al. 2018), they get enlarged to the sizes we observe today. Another scenario suggests that dwarfs with high-spin haloes can evolve into UDGs (Amorisco & Loeb 2016; Benavides et al. 2021). We refer to such scenarios as ‘puffed-up dwarf’ formation scenarios. Other scenarios rely on external processes to form UDGs via e.g. tidal stripping and heating (e.g. Carleton et al. 2019) and tidal interaction (e.g. Tremmel et al. 2020), which account for the apparent radial size excess. Some UDGs in the past have been suggested to be consistent with a tidal dwarf galaxy origin as well (see e.g. Duc et al. 2014; Buzzo et al. 2023; Gannon et al. 2023b). We refer to this formation scenario throughout the paper as ‘tidal UDGs’.

In addition to these, one scenario suggests that UDGs started their lifetimes destined to be large and massive but had their star formation truncated early on. As a consequence, their stellar masses do not increase at the same rate as their haloes, resulting in galaxies with dwarf-like stellar masses enclosed by halo masses similar to those of more massive M33-like galaxies (van Dokkum et al. 2015; Peng & Lim 2016; Danieli et al. 2022; Gannon et al. 2023a). For the remainder of the manuscript, we refer to this formation scenario as ‘failed galaxy’ scenario.

One crucial discriminant between these formation scenarios is the stellar populations of UDGs. For example, if they were formed by a puffed-up dwarf scenario, they would overall be expected to have similar stellar populations to those of classical dwarfs. One expectation, for example, is that they may follow the same scaling relations as classical dwarfs, e.g. the mass-metallicity relation (MZR). It is important to mention, however, that different puffed-up dwarf scenarios may lead to very different star formation histories (SFHs), i.e. supernovae feedback would lead to a bursty SFH, while high-spin haloes or tidal interactions may be more similar to an exponentially decaying SFH, but in both cases the UDGs are expected to preserve the stellar populations of their progenitors. Tidal UDGs are expected to be young, metal-rich, and gas-dominated for their stellar mass (Duc et al. 2014; Haslbauer et al. 2019). On the other hand, if UDGs are failed galaxies, the predictions for their stellar populations are less clear. They are not expected to have stellar populations similar to classical dwarfs (as classical dwarfs were not their progenitors), nor to more massive M33-like galaxies (as they did not evolve to become such). The only stellar population prediction for galaxies formed by this scenario is that they would have old or even ancient stellar populations (≥10 Gyr), as they by definition have suffered from early-quenching. For the same reason, they may be expected to have lower metallicities and possibly be alpha enhanced as they have had shorter time-scales to form stellar masses comparable to those of classical dwarfs that had prolonged SFHs (see e.g. Forbes et al. 2020), and also because, as suggested by Danieli et al. (2022), these galaxies could be mainly formed of disrupted GCs, thus being more metal-poor. As a key point in the comparison, the MZR is known to evolve with redshift (Ma et al. 2016; Chartab et al. 2024). Galaxies that follow the MZR at high redshift have stellar populations consistent with the chemical enrichment up to that point, thus, either galaxies at high redshift themselves or galaxies that have not had much chemical evolution since then (i.e. early quenched). We thus may expect failed galaxies to follow the MZRs of high-redshift galaxies.

Some of these predictions were recently confirmed by Ferré-Mateu et al. (2023) (hereafter, FM23), who used the largest spectroscopic study of UDGs to date (25) to show that there is a correlation between α-enrichment and the SFH of UDGs, further emphasized by the positioning of UDGs in the cluster phase-space diagram. They found that early-infall UDGs are the ones with the fastest star formation episodes and have higher [α/Fe], consistent with early-quenching scenarios, such as failed galaxy ones. Conversely, UDGs with prolonged SFHs are the late infallers, having lower [α/Fe] ratios, being better explained by puffed-up dwarf scenarios.

Additionally, Buzzo et al. (2022) (hereafter, B22) used spectral energy distribution (SED) fitting to show that there is a correlation between the age, metallicity, and environment that UDGs reside in. The older ones display the lowest metallicities, are consistent with the MZR at high redshift (z ~ 2), and are found in the densest environments. The younger ones (≤8 Gyr) follow the classical dwarf MZR, being more metal-rich and are found in less dense environments, such as the field and groups. These findings are in general agreement with the findings of Barbosa et al. (2020), who used SED fitting to study 100 UDGs in the field and to show that they are mostly young (~6 Gyr) and follow the classical dwarf MZR. In their sample, some UDGs were found to be more metal-rich than the classical dwarf MZR ([Fe/H] ≥ −0.5 dex) and extremely young (≤1 Gyr), similar to what has been found for some cluster UDGs (Ferré-Mateu et al. 2018; Ruiz-Lara et al. 2018). Their ages and metallicities are consistent with them being tidal UDGs (Collins & Read 2022).

B22 have also shown that the stellar populations of UDGs seem to correlate with their GC-richness. Old and metal-poor UDGs have on average the highest number of GCs, while younger and more metal-rich ones (consistent with the classical dwarf MZR) have the lowest GC numbers. This trend, however, was not found by FM23 using spectroscopy. Selection effects can be the cause of the differences, since both samples are not fully representative of the population of UDGs. Thus, to understand if this trend holds and whether it can or cannot be extrapolated to other UDGs, it is necessary to test it on a larger sample of UDGs. These would preferably be spread across different environments and span a variety of GC numbers so that selection effects can be diminished and conclusive interpretations made.

In this study, we extend the work of B22 to a sample of 59 UDGs in the MATLAS low- and moderate-density environments (Marleau et al. 2021), to help balancing the sample that was biased towards higher density environments in B22, and to thus start building up a representative sample of UDGs. 38 out of these 59 MATLAS UDGs have GC counts from *Hubble Space Telescope (HST)* imaging. Similar to B22, we employ optical to mid-infrared SED fitting to explore the stellar populations of these MATLAS UDGs. Buzzo et al. (in preparation) will combine both samples (i.e. B22 and this study) with a control sample of classical dwarf galaxies to perform an statistically meaningful and comprehensive study of UDGs spread across environments.

The paper is structured as follows: in Section 2 we present a summary of the UDG sample studied in this work and the data available for each UDG. In Section 3 we describe our methods to obtain the photometry and morphology of the galaxies and our SED fitting methodology. In Section 4 we provide our results. In Section 5 we discuss the implications of our results within the theoretical predictions for UDGs and as compared to the literature. In Section 6 we present the summary and the conclusions of the paper.

Throughout this paper, when converting distances to redshifts (or vice versa), we assume the cosmological parameters from the
Planck 2020 collaboration ($H_0 = 67.4 \pm 0.5$ km s$^{-1}$ Mpc$^{-1}$; $\Omega_m = 0.315 \pm 0.007$, Planck Collaboration VI 2020).

## 2 DATA SAMPLE

The MATLAS UDGs were identified by the deep optical imaging of the large observing program Mass Assembly of early-Type gaLAXies with their fine Structures (MATLAS) survey (Duc et al. 2014, 2015; Błek et al. 2020; Duc 2020; Marleau et al. 2021; Poulain et al. 2021). The sample of 59 UDGs used in this work was selected by Marleau et al. (2021), using the same UDG criteria as described above (i.e. $\mu_{g,0} > 24$ mag arcsec$^{-2}$, $R_e > 1.5$ kpc). These UDGs are all in low-to-moderate density environments. There are no UDGs in high-density environments such as clusters of galaxies in our sample. The environment that the galaxies reside in was calculated using a K-nearest neighbour algorithm. Using the 10 nearest galaxies to the UDGs, the local surface density ($\Sigma_{10}$) and local volume density ($\rho_{10}$) are calculated as a proxy for the environment. It therefore does not separate central from satellite UDGs. More details about the environment determination can be found in Duc et al. (2014) and Marleau et al. (2021).

For simplicity, throughout this paper we follow the MATLAS convention and refer to the closest massive galaxy to the UDGs as hosts, but that does not necessarily imply that the UDGs are satellites of the massive nearby galaxy.

No selection against star-forming UDGs was made, meaning that our sample has both star-forming and quiescent UDGs, differently from the sample of B22 and FM23, who focused only on quiescent UDGs.

As previously mentioned, the identification of the low surface brightness sources was carried out on the imaging fields around massive early-type galaxies (ETGs) in the nearby Universe. The UDGs are, thus, classified as being part of the group/in the vicinity of the massive ETG. However, since the UDG definition is based on physical size, it is imperative to understand the real distance to the galaxy to confirm whether or not it is a bonafide UDG. Few galaxies to date have been followed-up spectroscopically to test if they are indeed associated with the group or suffering from projection effects. The largest follow-up was carried out by Heesters et al. (2023) who studied 56 MATLAS dwarfs, out of which three were UDGs (MATLAS-585, MATLAS-2019, and MATLAS-2103). They found that 75 per cent of the studied galaxies were at the same distances as their hosts, including the three studied UDGs. In this study, in addition to the SED fitting, we present spectroscopic follow-up of three MATLAS UDG candidates with Keck/DEIMOS (MATLAS-342, MATLAS-368, and MATLAS-1059) to recover their redshifts and probe their association with the group/massive neighbour. The spectroscopic observations and results are described in Appendix A. We find that all three UDGs are at the same redshifts as their hosts, in agreement with the findings of Heesters et al. (2023). We therefore consider it reasonable to assume that the rest of the sample is also at the same redshift as their hosts. We further discuss the assumed redshifts for the MATLAS UDGs in Section 4.1.

The description of the imaging data used to perform SED fitting is in Section 2.1. The data used to perform the GC number counts are described in Section 2.2.

### 2.1 Photometric data

In this work, we use data from the optical to mid-infrared to study the stellar population properties of 59 UDGs in the MATLAS survey. Next, we present the data used for each galaxy, along with how the

photometry was measured in each band. We note that we tried to recover the photometry in the ultraviolet for the MATLAS UDGs using *Galaxy Evolution Explorer (GALEX)* imaging. However, all of our measurements turned out to be fainter than the $3\sigma$ limit of the survey, independently if the imaging was observed with the Deep (DIS, $3\sigma_{DIS}(NUV) = 24.9$ mag), Medium (MIS, $3\sigma_{MIS}(NUV) = 23.2$ mag), or All-Sky (AIS, $3\sigma_{AIS}(NUV) = 21.3$ mag) survey (Martin et al. 2005; Morrissey et al. 2007).

Optical, near- and mid-IR magnitude measurements are in AB magnitudes and were corrected for Galactic extinction using the two-dimensional dust maps of Schlegel, Finkbeiner & Davis 1998 (recalibrated by Schlafly & Finkbeiner 2011) and the extinction law of Calzetti et al. (2000).

#### 2.1.1 DECaLS optical imaging

Although the MATLAS UDGs were identified using Canada-France-Hawaii Telescope (CFHT, surface brightness limit of 28.5–29 mag arcsec$^{-2}$ in the g band, Duc et al. 2014; Poulain et al. 2021) data, in this study, we use data from the Dark Energy Camera Legacy Survey (DECaLS, surface brightness limit of 28.5 mag arcsec$^{-2}$ in the r band, Li et al. 2022) to perform the SED fitting. This is because DECaLS has imaging available in the g, r, i, and z bands for most galaxies in our sample, while CFHT has only the g and r bands (with a few rare cases where the i band is also available). We caution, none the less, that the DECaLS coadded imaging is noisy in the outskirts of galaxies (Li et al. 2022), which may lead to slight differences in the recovered photometry and morphological parameters when compared to CFHT, although consistent within uncertainties.

Differently from B22, in this study we do not use aperture photometry of our galaxies. Interested in obtaining total magnitudes in addition to morphological properties of the galaxies, we use the multiwavelength galaxy fitting (GALFITM, Häußler et al. 2013; Vika et al. 2013) routine to study the galaxies in the optical. Detailed explanation of the input parameters and configuration of GALFITM are given in Section 3.1.

For all of the galaxies, archival optical coadded data were obtained from the DECaLS (Dey et al. 2019). 36 out of our 59 UDGs were observed as part of the DECaLS Data Release 10, which includes the g, r, i, and z bands. The remaining 23 galaxies were observed as part of the DECaLS Data Release 9 and have imaging available in the g, r, and z bands. The reduction and calibration of the DECaLS data are described in Dey et al. (2019).

As discussed in B22, the coadded DECaLS data have shallower depths and more uncertain sky subtractions than other optical surveys focused on low surface brightness galaxies. None the less, tests carried out by B22 have shown that the photometry obtained with the DECaLS coadded data are consistent within $1\sigma$ with the photometry obtained by Lim et al. (2020) and Pandya et al. (2018) using deeper data reduced with a low surface brightness-appropriate pipeline. The photometry recovered from DECaLS is, however, on average 0.1 mag fainter than the ones in the literature. This systematic difference was incorporated into our final magnitudes as described in B22. The final photometry in the optical for all MATLAS UDGs was compared to that obtained by Poulain et al. (2021) using CFHT data and was found to provide similar results.

#### 2.1.2 WISE near-IR and mid-IR imaging

The *Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010)* is a space telescope that has imaged the entire sky in four filters
Planck 2020 collaboration ($H_0 = 67.4 \pm 0.5$ km s$^{-1}$ Mpc$^{-1}$; $\Omega_m = 0.315 \pm 0.007$, Planck Collaboration VI 2020).

## 2 DATA SAMPLE

The MATLAS UDGs were identified by the deep optical imaging of the large observing program Mass Assembly of early-Type gaLAXies with their fine Structures (MATLAS) survey (Duc et al. 2014, 2015; Błek et al. 2020; Duc 2020; Marleau et al. 2021; Poulain et al. 2021). The sample of 59 UDGs used in this work was selected by Marleau et al. (2021), using the same UDG criteria as described above (i.e. $\mu_{g,0} > 24$ mag arcsec$^{-2}$, $R_e > 1.5$ kpc). These UDGs are all in low-to-moderate density environments. There are no UDGs in high-density environments such as clusters of galaxies in our sample. The environment that the galaxies reside in was calculated using a K-nearest neighbour algorithm. Using the 10 nearest galaxies to the UDGs, the local surface density ($\Sigma_{10}$) and local volume density ($\rho_{10}$) are calculated as a proxy for the environment. It therefore does not separate central from satellite UDGs. More details about the environment determination can be found in Duc et al. (2014) and Marleau et al. (2021).

For simplicity, throughout this paper we follow the MATLAS convention and refer to the closest massive galaxy to the UDGs as hosts, but that does not necessarily imply that the UDGs are satellites of the massive nearby galaxy.

No selection against star-forming UDGs was made, meaning that our sample has both star-forming and quiescent UDGs, differently from the sample of B22 and FM23, who focused only on quiescent UDGs.

As previously mentioned, the identification of the low surface brightness sources was carried out on the imaging fields around massive early-type galaxies (ETGs) in the nearby Universe. The UDGs are, thus, classified as being part of the group/in the vicinity of the massive ETG. However, since the UDG definition is based on physical size, it is imperative to understand the real distance to the galaxy to confirm whether or not it is a bonafide UDG. Few galaxies to date have been followed-up spectroscopically to test if they are indeed associated with the group or suffering from projection effects. The largest follow-up was carried out by Heesters et al. (2023) who studied 56 MATLAS dwarfs, out of which three were UDGs (MATLAS-585, MATLAS-2019, and MATLAS-2103). They found that 75 per cent of the studied galaxies were at the same distances as their hosts, including the three studied UDGs. In this study, in addition to the SED fitting, we present spectroscopic follow-up of three MATLAS UDG candidates with Keck/DEIMOS (MATLAS-342, MATLAS-368, and MATLAS-1059) to recover their redshifts and probe their association with the group/massive neighbour. The spectroscopic observations and results are described in Appendix A. We find that all three UDGs are at the same redshifts as their hosts, in agreement with the findings of Heesters et al. (2023). We therefore consider it reasonable to assume that the rest of the sample is also at the same redshift as their hosts. We further discuss the assumed redshifts for the MATLAS UDGs in Section 4.1.

The description of the imaging data used to perform SED fitting is in Section 2.1. The data used to perform the GC number counts are described in Section 2.2.

### 2.1 Photometric data

In this work, we use data from the optical to mid-infrared to study the stellar population properties of 59 UDGs in the MATLAS survey. Next, we present the data used for each galaxy, along with how the

photometry was measured in each band. We note that we tried to recover the photometry in the ultraviolet for the MATLAS UDGs using *Galaxy Evolution Explorer (GALEX)* imaging. However, all of our measurements turned out to be fainter than the $3\sigma$ limit of the survey, independently if the imaging was observed with the Deep (DIS, $3\sigma_{DIS}(NUV) = 24.9$ mag), Medium (MIS, $3\sigma_{MIS}(NUV) = 23.2$ mag), or All-Sky (AIS, $3\sigma_{AIS}(NUV) = 21.3$ mag) survey (Martin et al. 2005; Morrissey et al. 2007).

Optical, near- and mid-IR magnitude measurements are in AB magnitudes and were corrected for Galactic extinction using the two-dimensional dust maps of Schlegel, Finkbeiner & Davis 1998 (recalibrated by Schlafly & Finkbeiner 2011) and the extinction law of Calzetti et al. (2000).

#### 2.1.1 DECaLS optical imaging

Although the MATLAS UDGs were identified using Canada-France-Hawaii Telescope (CFHT, surface brightness limit of 28.5–29 mag arcsec$^{-2}$ in the g band, Duc et al. 2014; Poulain et al. 2021) data, in this study, we use data from the Dark Energy Camera Legacy Survey (DECaLS, surface brightness limit of 28.5 mag arcsec$^{-2}$ in the r band, Li et al. 2022) to perform the SED fitting. This is because DECaLS has imaging available in the g, r, i, and z bands for most galaxies in our sample, while CFHT has only the g and r bands (with a few rare cases where the i band is also available). We caution, none the less, that the DECaLS coadded imaging is noisy in the outskirts of galaxies (Li et al. 2022), which may lead to slight differences in the recovered photometry and morphological parameters when compared to CFHT, although consistent within uncertainties.

Differently from B22, in this study we do not use aperture photometry of our galaxies. Interested in obtaining total magnitudes in addition to morphological properties of the galaxies, we use the multiwavelength galaxy fitting (GALFITM, Häußler et al. 2013; Vika et al. 2013) routine to study the galaxies in the optical. Detailed explanation of the input parameters and configuration of GALFITM are given in Section 3.1.

For all of the galaxies, archival optical coadded data were obtained from the DECaLS (Dey et al. 2019). 36 out of our 59 UDGs were observed as part of the DECaLS Data Release 10, which includes the g, r, i, and z bands. The remaining 23 galaxies were observed as part of the DECaLS Data Release 9 and have imaging available in the g, r, and z bands. The reduction and calibration of the DECaLS data are described in Dey et al. (2019).

As discussed in B22, the coadded DECaLS data have shallower depths and more uncertain sky subtractions than other optical surveys focused on low surface brightness galaxies. None the less, tests carried out by B22 have shown that the photometry obtained with the DECaLS coadded data are consistent within $1\sigma$ with the photometry obtained by Lim et al. (2020) and Pandya et al. (2018) using deeper data reduced with a low surface brightness-appropriate pipeline. The photometry recovered from DECaLS is, however, on average 0.1 mag fainter than the ones in the literature. This systematic difference was incorporated into our final magnitudes as described in B22. The final photometry in the optical for all MATLAS UDGs was compared to that obtained by Poulain et al. (2021) using CFHT data and was found to provide similar results.

#### 2.1.2 WISE near-IR and mid-IR imaging

The *Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010)* is a space telescope that has imaged the entire sky in four filters
with effective wavelengths of 3.4, 4.6, 12, and 22 µm (near to mid-infrared). For this study, we gathered *WISE* data for all the galaxies in our sample. These data are a mix of archival ALLWISE data and bespoke data construction and analysis, including custom mosaic construction from *WISE* single frames. The reduction, calibration, and photometric measurement processes are thoroughly described in B22. We use *WISE* in all of our SED fits as it has been shown by B22 to significantly improve the stellar population results and to help breaking the age-metallicity degeneracy. All photometric measurements can be found in Table B1.

## 2.2 GC numbers

Total GC numbers have been obtained for 38 galaxies in our sample using *HST/ACS* imaging. The reduction, source detection, GC candidate selection, and final GC counts are thoroughly discussed in Marleau et al. (submitted). Here we use the total (background and completeness corrected) GC numbers of the MATLAS UDGs to test the trend between GC-richness and metallicity found by B22 on our more complete sample. In B22 and FM23, UDGs considered to be GC-rich were those with more than 20 GCs, as this value roughly corresponds to a halo mass of $10^{11} M_{\odot}$ (Burkert & Forbes 2020). UDGs with less than 20 GCs were considered GC-poor. In this study, we do not impose this hard separation at $N_{GC} = 20$. Instead, we analyse GC numbers in a continuous manner so that evolutionary trends can be identified.

# 3 ANALYSIS

## 3.1 GALFIT

In this study, we use GALFIT (Häußler et al. 2013; Vika et al. 2013) to perform simultaneous galaxy fitting of our optical bands. This technique allows the combination of lower signal-to-noise (S/N) images, such as in the z band, with images with high S/N, and balances them out, resulting in consistent and well-constrained parameters.

We note that GALFIT was only run on the optical images, not on the *WISE* ones. This is because our methodology (described in Jarrett et al. 2012, 2013, 2019; Buzzo et al. 2022) was already shown to provide reliable total magnitudes for the *WISE* imaging.

For all galaxies, a single Sérsic model was fitted in the optical bands, using as initial guesses the morphological parameters obtained by Poulain et al. (2021). In our GALFIT set-up, the magnitudes were allowed to vary freely amongst bands. Alternatively, the effective radius ($R_e$), Sérsic index ($n$), axis ratio ($b/a$), and position angle (PA) were allowed to vary, with the condition that they are constant in all bands. A discussion on appropriate degrees of freedom to use for each free parameter in GALFIT is provided by Buzzo et al. (2021).

The sky value for each wavelength band was obtained from the image header. To run GALFIT, we created synthetic PSFs of each image using the *PSFex* routine (Bertin 2002). The parameters used to create such PSFs were obtained using the *Source Extractor* routine (*SEXTRACTOR*, Bertin 2002, version 2.19.5). We masked all interlopers, including GCs and nuclear star clusters, using the segmentation image output of *SEXTRACTOR*. The routine was run in single-mode in each band. We defined a source detection minimum area of 3 pixels and a threshold of $3\sigma$ above the background to detect sources. *SEXTRACTOR* was run with an automatic background fitting, with a background cell size of 64 pixels.

For the masking process, we visually inspected every mask to ensure no star-forming region was being masked as this could bias the final results to older stellar populations. To do this, we carefully looked into *GALEX* data on every one of the galaxies to see if any star formation was detected. As previously mentioned, for none of the galaxies UV emission was detected, indicating that the process was not masking any star-forming region.

In Fig. 1 we show a GALFIT model for the four optical bands of the galaxy MATLAS-2019 (NGC 5846, UDG1). This example shows the highest mask fraction applied to our sample, as MATLAS-2019 is the galaxy with the most GCs within the sample of MATLAS UDGs. All other galaxies have had none to very few sources masked in front of them, therefore not compromising the final recovered photometry of the galaxies. GALFIT models for the other galaxies are of similar quality to the one shown in Fig. 1.

Results from these fits are given in Table B2. We find that the MATLAS UDGs have a median Sérsic index of $n = 0.95 \pm 0.14$, and a median axis ratio of $b/a = 0.62 \pm 0.13$. Poulain et al. (2021) also obtained morphological parameters for all of the MATLAS dwarfs. Their results for the MATLAS UDGs show a median Sérsic index of $n = 0.79 \pm 0.10$ and a median axis ratio of $b/a = 0.67 \pm 0.14$. Thus, their morphological findings for the MATLAS UDGs are consistent with ours. GALFIT, similarly to GALFIT, provides unrealistically small uncertainties. In Appendix C, we discuss some of the implications of using these underestimated uncertainties to the results obtained with SED fitting.

## 3.2 PROSPECTOR

For the SED fitting, we run the fully Bayesian Markov Chain Monte Carlo (MCMC) based inference code PROSPECTOR (version 1.2.0¹ Leja et al. 2017; Johnson et al. 2021b), complemented by the FLEXIBLE STELLAR POPULATION SYNTHESIS package (FSPS; Conroy, Gunn & White 2009; Conroy & Gunn 2010; Conroy, White & Gunn 2010, version 3.2). To sample the posteriors, we used the dynamic nestled sampling (Skilling 2004; Higson et al. 2019) algorithm *dynesty* (Speagle 2020). *dynesty* was configured using 100 samples, 100 live points, and a 0.05 tolerance when finishing the baseline run.

A full and thorough description of the main configuration and models used in PROSPECTOR is available in B22. In the latter, PROSPECTOR models including dust attenuation provide results closer to those found with spectroscopy for UDGs. Thus, in this study, we only use models that include dust attenuation, assuming the Gordon et al. (2003) attenuation curve. Whether this dust is artificially added to improve the models or the galaxies actually have dust can only be probed using far infrared or ultraviolet deep data, which will be pursued in the future. We note, none the less, that tests with PROSPECTOR found a dust attenuation of $A_V \sim 0.1-0.2$ mag for Milky Way GCs when there should be none, suggesting an artificial addition of dust in PROSPECTOR models (Johnson et al. 2021b, fig. 6).

For our fits of all galaxies, we include upper limit fluxes coming from the 12 and 22 µm bands from *WISE*. For cases of extreme low S/N in the z band of DECaLS, this band was also incorporated in the fitting as an upper limit. For further understanding of how PROSPECTOR deals with upper limits, see appendix A of Sawicki (2012). The inclusion of the upper limits was shown in Buzzo et al. (2022) to be helpful in constraining the amount of dust attenuation

<footnote n="1">¹We note that tests carried out using PROSPECTOR version 1.0 have delivered systematically younger ages and longer star formation time-scales.</footnote>
with effective wavelengths of 3.4, 4.6, 12, and 22 µm (near to mid-infrared). For this study, we gathered *WISE* data for all the galaxies in our sample. These data are a mix of archival ALLWISE data and bespoke data construction and analysis, including custom mosaic construction from *WISE* single frames. The reduction, calibration, and photometric measurement processes are thoroughly described in B22. We use *WISE* in all of our SED fits as it has been shown by B22 to significantly improve the stellar population results and to help breaking the age-metallicity degeneracy. All photometric measurements can be found in Table B1.

## 2.2 GC numbers

Total GC numbers have been obtained for 38 galaxies in our sample using *HST/ACS* imaging. The reduction, source detection, GC candidate selection, and final GC counts are thoroughly discussed in Marleau et al. (submitted). Here we use the total (background and completeness corrected) GC numbers of the MATLAS UDGs to test the trend between GC-richness and metallicity found by B22 on our more complete sample. In B22 and FM23, UDGs considered to be GC-rich were those with more than 20 GCs, as this value roughly corresponds to a halo mass of $10^{11} M_{\odot}$ (Burkert & Forbes 2020). UDGs with less than 20 GCs were considered GC-poor. In this study, we do not impose this hard separation at $N_{GC} = 20$. Instead, we analyse GC numbers in a continuous manner so that evolutionary trends can be identified.

# 3 ANALYSIS

## 3.1 GALFIT

In this study, we use GALFIT (Häußler et al. 2013; Vika et al. 2013) to perform simultaneous galaxy fitting of our optical bands. This technique allows the combination of lower signal-to-noise (S/N) images, such as in the z band, with images with high S/N, and balances them out, resulting in consistent and well-constrained parameters.

We note that GALFIT was only run on the optical images, not on the *WISE* ones. This is because our methodology (described in Jarrett et al. 2012, 2013, 2019; Buzzo et al. 2022) was already shown to provide reliable total magnitudes for the *WISE* imaging.

For all galaxies, a single Sérsic model was fitted in the optical bands, using as initial guesses the morphological parameters obtained by Poulain et al. (2021). In our GALFIT set-up, the magnitudes were allowed to vary freely amongst bands. Alternatively, the effective radius ($R_e$), Sérsic index ($n$), axis ratio ($b/a$), and position angle (PA) were allowed to vary, with the condition that they are constant in all bands. A discussion on appropriate degrees of freedom to use for each free parameter in GALFIT is provided by Buzzo et al. (2021).

The sky value for each wavelength band was obtained from the image header. To run GALFIT, we created synthetic PSFs of each image using the *PSFex* routine (Bertin 2002). The parameters used to create such PSFs were obtained using the *Source Extractor* routine (*SEXTRACTOR*, Bertin 2002, version 2.19.5). We masked all interlopers, including GCs and nuclear star clusters, using the segmentation image output of *SEXTRACTOR*. The routine was run in single-mode in each band. We defined a source detection minimum area of 3 pixels and a threshold of $3\sigma$ above the background to detect sources. *SEXTRACTOR* was run with an automatic background fitting, with a background cell size of 64 pixels.

For the masking process, we visually inspected every mask to ensure no star-forming region was being masked as this could bias the final results to older stellar populations. To do this, we carefully looked into *GALEX* data on every one of the galaxies to see if any star formation was detected. As previously mentioned, for none of the galaxies UV emission was detected, indicating that the process was not masking any star-forming region.

In Fig. 1 we show a GALFIT model for the four optical bands of the galaxy MATLAS-2019 (NGC 5846, UDG1). This example shows the highest mask fraction applied to our sample, as MATLAS-2019 is the galaxy with the most GCs within the sample of MATLAS UDGs. All other galaxies have had none to very few sources masked in front of them, therefore not compromising the final recovered photometry of the galaxies. GALFIT models for the other galaxies are of similar quality to the one shown in Fig. 1.

Results from these fits are given in Table B2. We find that the MATLAS UDGs have a median Sérsic index of $n = 0.95 \pm 0.14$, and a median axis ratio of $b/a = 0.62 \pm 0.13$. Poulain et al. (2021) also obtained morphological parameters for all of the MATLAS dwarfs. Their results for the MATLAS UDGs show a median Sérsic index of $n = 0.79 \pm 0.10$ and a median axis ratio of $b/a = 0.67 \pm 0.14$. Thus, their morphological findings for the MATLAS UDGs are consistent with ours. GALFIT, similarly to GALFIT, provides unrealistically small uncertainties. In Appendix C, we discuss some of the implications of using these underestimated uncertainties to the results obtained with SED fitting.

## 3.2 PROSPECTOR

For the SED fitting, we run the fully Bayesian Markov Chain Monte Carlo (MCMC) based inference code PROSPECTOR (version 1.2.0¹ Leja et al. 2017; Johnson et al. 2021b), complemented by the FLEXIBLE STELLAR POPULATION SYNTHESIS package (FSPS; Conroy, Gunn & White 2009; Conroy & Gunn 2010; Conroy, White & Gunn 2010, version 3.2). To sample the posteriors, we used the dynamic nestled sampling (Skilling 2004; Higson et al. 2019) algorithm *dynesty* (Speagle 2020). *dynesty* was configured using 100 samples, 100 live points, and a 0.05 tolerance when finishing the baseline run.

A full and thorough description of the main configuration and models used in PROSPECTOR is available in B22. In the latter, PROSPECTOR models including dust attenuation provide results closer to those found with spectroscopy for UDGs. Thus, in this study, we only use models that include dust attenuation, assuming the Gordon et al. (2003) attenuation curve. Whether this dust is artificially added to improve the models or the galaxies actually have dust can only be probed using far infrared or ultraviolet deep data, which will be pursued in the future. We note, none the less, that tests with PROSPECTOR found a dust attenuation of $A_V \sim 0.1-0.2$ mag for Milky Way GCs when there should be none, suggesting an artificial addition of dust in PROSPECTOR models (Johnson et al. 2021b, fig. 6).

For our fits of all galaxies, we include upper limit fluxes coming from the 12 and 22 µm bands from *WISE*. For cases of extreme low S/N in the z band of DECaLS, this band was also incorporated in the fitting as an upper limit. For further understanding of how PROSPECTOR deals with upper limits, see appendix A of Sawicki (2012). The inclusion of the upper limits was shown in Buzzo et al. (2022) to be helpful in constraining the amount of dust attenuation

<footnote n="1">¹We note that tests carried out using PROSPECTOR version 1.0 have delivered systematically younger ages and longer star formation time-scales.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. GALFITM model of MATLAS-2019 (NGC 5846.UDG1). The first column shows the input image. The second column is the mask used in the model. The third column is the model. The fourth column is the residual (input – model). The rows show the modelling in the g, r, i, and z bands, respectively.</figcaption></figure>

found for these galaxies. In the cases where no dust is found, the upper limits do not play an important role in the fit.

We assume a delayed exponentially declining SFH, as this is consistent with the recovered shape of the SFHs of UDGs by FM23. The delayed model is an extension of the regular exponentially declining SFH. The latter usually assumes that star formation jumps from zero to its maximum value at the time $t_{age}$ and then declines exponentially within a time-scale $\tau$. The delayed model, alternatively, multiplies the exponentially declining SFH by the time since $t_{age}$. This is capable of removing the discontinuity in the SFH at $t_{age}$ and the condition that star formation can only decline after that point. This

results in a more flexible, robust, and physical SFH (see Carnall et al. 2018; Leja et al. 2019, for thorough discussions and comparisons between parametric SFHs). Thus, the form of the chosen SFH is

$$SFR(t) \propto \begin{cases} (t - t_{age}) \exp(-\frac{t-t_{age}}{\tau}), & \text{if } t > t_{age} \\ 0, & \text{if } t < t_{age} \end{cases} \quad (1)$$

According to this definition, $t_{age}$ measures the onset of star formation within a galaxy and $\tau$ is an approximation of the star formation time-scale, i.e. how long does a galaxy take to quench after reaching peak star formation. In this study, we are interested in
<figure><img src="image_2.png" /><figcaption>Figure 1. GALFITM model of MATLAS-2019 (NGC 5846.UDG1). The first column shows the input image. The second column is the mask used in the model. The third column is the model. The fourth column is the residual (input – model). The rows show the modelling in the g, r, i, and z bands, respectively.</figcaption></figure>

found for these galaxies. In the cases where no dust is found, the upper limits do not play an important role in the fit.

We assume a delayed exponentially declining SFH, as this is consistent with the recovered shape of the SFHs of UDGs by FM23. The delayed model is an extension of the regular exponentially declining SFH. The latter usually assumes that star formation jumps from zero to its maximum value at the time $t_{age}$ and then declines exponentially within a time-scale $\tau$. The delayed model, alternatively, multiplies the exponentially declining SFH by the time since $t_{age}$. This is capable of removing the discontinuity in the SFH at $t_{age}$ and the condition that star formation can only decline after that point. This

results in a more flexible, robust, and physical SFH (see Carnall et al. 2018; Leja et al. 2019, for thorough discussions and comparisons between parametric SFHs). Thus, the form of the chosen SFH is

$$SFR(t) \propto \begin{cases} (t - t_{age}) \exp(-\frac{t - t_{age}}{\tau}), & \text{if } t > t_{age} \\ 0, & \text{if } t < t_{age} \end{cases} \quad (1)$$

According to this definition, $t_{age}$ measures the onset of star formation within a galaxy and $\tau$ is an approximation of the star formation time-scale, i.e. how long does a galaxy take to quench after reaching peak star formation. In this study, we are interested in
the mass-weighted ages ($t_M$) for the galaxies, which are calculated analytically within PROSPECTOR using as input parameters the age since the onset of star formation ($t_{age}$) and the star formation time-scale ($\tau$). We note that B22 has incorrectly reported $t_{age}$ as the mass-weighted age, not the true $t_M$ calculated analytically with PROSPECTOR. $t_M$ is only strongly different from $t_{age}$ in the cases where $\tau$ is long, which is usually not the case for UDGs. Thus, the bulk of the comparisons between $t_{age}$ from B22 and spectroscopic mass-weighted ages from the literature hold. However, some individual cases may have larger differences and be thus less consistent with spectroscopy than reported. Buzzo et al. (in preparation) is working on combining the data set in this study with the one from B22, where the photometry for all galaxies will be obtained consistently with GALFITM and the true $t_M$ will be reported for every galaxy so that we have a homogeneous data set, where fair comparisons can be made and conclusions can be assessed.

For this study, we use two different configurations in PROSPECTOR.

(i) $A_V \neq 0$; $z = z_{host}$: five free parameters. Stellar mass (log $M_*/M_\odot$), total metallicity ([M/H]$^2$), age ($t_{age}$), star formation time-scale ($\tau$), and diffuse interstellar dust ($A_V$). Redshifts ($z$) are fixed to the redshift of the group/massive galaxy where they were identified. The redshifts used come from Marleau et al. (2021).

(ii) $A_V \neq 0$; $z \neq z_{host}$: six free parameters (log($M_*/M_\odot$), [M/H], $t_{age}$, $\tau$, $A_V$, and $z$). In this case, we leave the redshifts of the UDGs free.

For these two scenarios, we placed linearly uniform priors on our free parameters. The range of the priors in our models are determined by the coverage of the Padova isochrones used in FSPS (Marigo & Girardi 2007; Marigo et al. 2008). These are log($M_*/M_\odot$) = 6–10, [M/H] = -2.0 to 0.2 dex, $\tau$ = 0.1–10 Gyr, $t_{age}$ = 0.1–14 Gyr, $A_V$ = 0–4.344 mag, and redshift $z$ = 0–0.045. This redshift range translates to a luminosity distance range of 0 Mpc < $D_L$ < 200 Mpc. Since all of the galaxies were identified in groups out to 45 Mpc, we assume that this redshift range is representative of the distance of all UDGs. We note that the default stellar masses obtained from PROSPECTOR are the total masses formed in the galaxy throughout its life. These were converted to surviving masses (most commonly quoted in the literature) in this work using a built-in function within PROSPECTOR.

As mentioned in B22, we emphasize that different prior assumptions do not significantly alter our results or conclusions, with the exception of the prior assumption on the shape of the SFH, as discussed by Webb et al. (2022).

# 4 RESULTS

Here we present our results on the stellar population properties of the 59 MATLAS UDGs studied in this work. As explained in Section 3.2, we carried out two model configurations in PROSPECTOR, one where the redshift is a free parameter and one with the redshift fixed. In this section we show our results for both configurations, discuss the differences between them and the implications of using one or the other.

The resulting SEDs from both configurations of every galaxy were inspected to ensure that every photometric point was well modelled and that the shape was as expected. We found a median $\chi^2_{reduced}$ of $1.41 \pm 2.65$. The best fits are for MATLAS-2103 and MATLAS-585, both with $\chi^2_{reduced} = 1.01$ and the worst is for MATLAS-1589, with

$\chi^2_{reduced} = 7.34$. Fig. 2 shows a typical fit of an UDG in our sample, with $\chi^2_{reduced} = 1.32$. Upper limits were not included in the calculation of the $\chi^2$. Figs 3 and 4 show the distribution of the parameters for all of the galaxies in our sample in both configurations.

It is important to notice that choosing a parametric SFH for the galaxies comes with many limitations. One of such is that the SFHs of the galaxies may have multiple components, i.e. an old quiescent and another with a recent burst of star formation. In such cases, a model with a parametric SFH that smooths over two distinct star formation episodes will be the most statistically correct, but will not reflect the real stellar populations present in the galaxies. Fitting the galaxies with a non-parametric SFH may be a solution for this problem, but it adds another one: too many free parameters for a limited amount of data, likely resulting in overfitting. To try and avoid this type of problem, we have gathered GALEX data for all of the MATLAS UDGs to be able to correctly fit the star-forming regions if the galaxies had any. None of the galaxies in our sample were detected in the FUV or NUV bands of GALEX, and provided only very faint upper limits. We have tested including these UV upper limits on a subsample of our galaxies to understand the effect they would have on the recovered stellar populations. For all galaxies, the results are consistent with the ones obtained without the usage of GALEX, indicating that PROSPECTOR is correctly tracing the dominant stellar populations within the galaxies only by using the optical to mid-IR data.

## 4.1 Comparison between models with free and fixed redshifts

In Fig. 2, we show a comparison between the output best-fitting SED and covariance matrix for the two configurations of PROSPECTOR applied to one of the galaxies in our sample, MATLAS-2019. It can be seen that the model where the redshift is free has a much larger spread in the determination of the stellar mass, while the results with the redshift fixed seem to be better constrained. Apart from the stellar mass, the other parameters do not seem to be strongly affected by the freedom of the redshift.

The same behaviour can be observed if we look into the whole population of galaxies, instead of an individual case. Figs 3 and 4 provide a comparison of all of the stellar population properties of the MATLAS UDGs obtained with both configurations. As can be seen in these figures, the metallicity, star formation time-scale, age, and dust attenuation are not strongly affected by the freedom of fitting or lack thereof to fit redshift in the models. Galaxies are found to be on average younger and more metal-rich in the models where the redshift is fixed, but consistent within the uncertainties with the stellar populations found in the models with the redshift free. However, the stellar mass of the galaxies shows a larger difference when comparing models. Overall, galaxies have higher stellar masses in the models where the redshift is a free parameter. That is because in integrated SEDs at non-cosmological distances, the redshift and stellar mass are highly degenerate parameters.

As explained in Section 2, we expect that most UDGs in our sample are at the same redshift as their hosts in the models with the redshift free, if our redshift estimates are reliable. We explore the difference between the redshift of the host and our recovered photometric redshift ($\delta_z$) in Fig. 5. As can be seen, while there is a peak at $\delta_z = 0$, a significant portion (>70 per cent of the UDGs have $|\delta_z|>0.0067$, equivalent to a difference in velocity of $|\delta_v|>200$ km s$^{-1}$). Thus, they do not have photometric redshifts consistent with those of their hosts. Most of them are in fact found in their background (i.e. $\delta_z < 0$). Since photometric redshifts at

<footnote n="2">²We use [M/H] rather than [Z/H] (as used in B22), because Z denotes a mass fraction, and therefore [M/H] is less confusing.</footnote>
the mass-weighted ages ($t_M$) for the galaxies, which are calculated analytically within PROSPECTOR using as input parameters the age since the onset of star formation ($t_{age}$) and the star formation time-scale ($\tau$). We note that B22 has incorrectly reported $t_{age}$ as the mass-weighted age, not the true $t_M$ calculated analytically with PROSPECTOR. $t_M$ is only strongly different from $t_{age}$ in the cases where $\tau$ is long, which is usually not the case for UDGs. Thus, the bulk of the comparisons between $t_{age}$ from B22 and spectroscopic mass-weighted ages from the literature hold. However, some individual cases may have larger differences and be thus less consistent with spectroscopy than reported. Buzzo et al. (in preparation) is working on combining the data set in this study with the one from B22, where the photometry for all galaxies will be obtained consistently with GALFITM and the true $t_M$ will be reported for every galaxy so that we have a homogeneous data set, where fair comparisons can be made and conclusions can be assessed.

For this study, we use two different configurations in PROSPECTOR.

(i) $A_V \neq 0$; $z = z_{host}$: five free parameters. Stellar mass (log $M_*/M_\odot$), total metallicity ([M/H]$^2$), age ($t_{age}$), star formation time-scale ($\tau$), and diffuse interstellar dust ($A_V$). Redshifts ($z$) are fixed to the redshift of the group/massive galaxy where they were identified. The redshifts used come from Marleau et al. (2021).

(ii) $A_V \neq 0$; $z \neq z_{host}$: six free parameters (log($M_*/M_\odot$), [M/H], $t_{age}$, $\tau$, $A_V$, and $z$). In this case, we leave the redshifts of the UDGs free.

For these two scenarios, we placed linearly uniform priors on our free parameters. The range of the priors in our models are determined by the coverage of the Padova isochrones used in FSPS (Marigo & Girardi 2007; Marigo et al. 2008). These are log($M_*/M_\odot$) = 6–10, [M/H] = -2.0 to 0.2 dex, $\tau$ = 0.1–10 Gyr, $t_{age}$ = 0.1–14 Gyr, $A_V$ = 0–4.344 mag, and redshift $z$ = 0–0.045. This redshift range translates to a luminosity distance range of 0 Mpc < $D_L$ < 200 Mpc. Since all of the galaxies were identified in groups out to 45 Mpc, we assume that this redshift range is representative of the distance of all UDGs. We note that the default stellar masses obtained from PROSPECTOR are the total masses formed in the galaxy throughout its life. These were converted to surviving masses (most commonly quoted in the literature) in this work using a built-in function within PROSPECTOR.

As mentioned in B22, we emphasize that different prior assumptions do not significantly alter our results or conclusions, with the exception of the prior assumption on the shape of the SFH, as discussed by Webb et al. (2022).

# 4 RESULTS

Here we present our results on the stellar population properties of the 59 MATLAS UDGs studied in this work. As explained in Section 3.2, we carried out two model configurations in PROSPECTOR, one where the redshift is a free parameter and one with the redshift fixed. In this section we show our results for both configurations, discuss the differences between them and the implications of using one or the other.

The resulting SEDs from both configurations of every galaxy were inspected to ensure that every photometric point was well modelled and that the shape was as expected. We found a median $\chi^2_{reduced}$ of $1.41 \pm 2.65$. The best fits are for MATLAS-2103 and MATLAS-585, both with $\chi^2_{reduced} = 1.01$ and the worst is for MATLAS-1589, with

$\chi^2_{reduced} = 7.34$. Fig. 2 shows a typical fit of an UDG in our sample, with $\chi^2_{reduced} = 1.32$. Upper limits were not included in the calculation of the $\chi^2$. Figs 3 and 4 show the distribution of the parameters for all of the galaxies in our sample in both configurations.

It is important to notice that choosing a parametric SFH for the galaxies comes with many limitations. One of such is that the SFHs of the galaxies may have multiple components, i.e. an old quiescent and another with a recent burst of star formation. In such cases, a model with a parametric SFH that smooths over two distinct star formation episodes will be the most statistically correct, but will not reflect the real stellar populations present in the galaxies. Fitting the galaxies with a non-parametric SFH may be a solution for this problem, but it adds another one: too many free parameters for a limited amount of data, likely resulting in overfitting. To try and avoid this type of problem, we have gathered GALEX data for all of the MATLAS UDGs to be able to correctly fit the star-forming regions if the galaxies had any. None of the galaxies in our sample were detected in the FUV or NUV bands of GALEX, and provided only very faint upper limits. We have tested including these UV upper limits on a subsample of our galaxies to understand the effect they would have on the recovered stellar populations. For all galaxies, the results are consistent with the ones obtained without the usage of GALEX, indicating that PROSPECTOR is correctly tracing the dominant stellar populations within the galaxies only by using the optical to mid-IR data.

## 4.1 Comparison between models with free and fixed redshifts

In Fig. 2, we show a comparison between the output best-fitting SED and covariance matrix for the two configurations of PROSPECTOR applied to one of the galaxies in our sample, MATLAS-2019. It can be seen that the model where the redshift is free has a much larger spread in the determination of the stellar mass, while the results with the redshift fixed seem to be better constrained. Apart from the stellar mass, the other parameters do not seem to be strongly affected by the freedom of the redshift.

The same behaviour can be observed if we look into the whole population of galaxies, instead of an individual case. Figs 3 and 4 provide a comparison of all of the stellar population properties of the MATLAS UDGs obtained with both configurations. As can be seen in these figures, the metallicity, star formation time-scale, age, and dust attenuation are not strongly affected by the freedom of fitting or lack thereof to fit redshift in the models. Galaxies are found to be on average younger and more metal-rich in the models where the redshift is fixed, but consistent within the uncertainties with the stellar populations found in the models with the redshift free. However, the stellar mass of the galaxies shows a larger difference when comparing models. Overall, galaxies have higher stellar masses in the models where the redshift is a free parameter. That is because in integrated SEDs at non-cosmological distances, the redshift and stellar mass are highly degenerate parameters.

As explained in Section 2, we expect that most UDGs in our sample are at the same redshift as their hosts in the models with the redshift free, if our redshift estimates are reliable. We explore the difference between the redshift of the host and our recovered photometric redshift ($\delta_z$) in Fig. 5. As can be seen, while there is a peak at $\delta_z = 0$, a significant portion (>70 per cent of the UDGs have $|\delta_z|>0.0067$, equivalent to a difference in velocity of $|\delta_v|>200$ km s$^{-1}$). Thus, they do not have photometric redshifts consistent with those of their hosts. Most of them are in fact found in their background (i.e. $\delta_z < 0$). Since photometric redshifts at

<footnote n="2">²We use [M/H] rather than [Z/H] (as used in B22), because Z denotes a mass fraction, and therefore [M/H] is less confusing.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** PROSPECTOR best-fitting SED and covariance matrix for MATLAS-2019. In all panels, purple colours stand for the models with the redshift free ($z \neq z_{\text{host}}$) and green colours are the models with the redshift fixed to the host ($z = z_{\text{host}}$). Top right: Best-fitting spectra of MATLAS-2019 in both configurations. The red points are the observed magnitudes of the galaxy in the *g*, *r*, *i*, *z*, W1, W2, W3, and W4 bands, respectively. The yellow curves show the DECaLS transmission curves. The blue curves are the WISE transmission curves. The remaining panels show the MCMC corner plot comparing the posterior distribution of the two configurations. The first panel in each column shows the 1D posterior distribution of the fitted parameter, while the remaining panels show the correlation between parameters. This image can be read and interpreted as a covariance matrix. Columns stand for redshift, stellar mass, metallicity, star formation time-scale, mass-weighted age, and interstellar diffuse dust extinction. Metallicity, star formation time-scale, age, and dust attenuation are similar in both configurations. However, the stellar mass is strongly affected by the freedom of the redshift, as these two parameters are highly degenerate.</figcaption></figure>

non-cosmological distances are very uncertain, we consider $\delta_z$ being different than zero as a confirmation that redshifts with this method cannot be trusted rather than an indication that MATLAS UDGs are not at the redshifts of their hosts.

Going forward, we treat the models with fixed redshift as the best representation of our galaxies. Thus, from now on, we will only

discuss the results of this configuration. We reinforce, none the less, that even if the galaxies are not at this assumed distance, changes in the stellar populations (except stellar mass) between models are small and well within their uncertainties (see Figs 2, 3, and 4). We carefully checked that the changes in stellar mass do not strongly affect any of our conclusions that follow. Specifically, we have checked that
<figure><img src="image_3.png" /><figcaption>**Figure 2.** PROSPECTOR best-fitting SED and covariance matrix for MATLAS-2019. In all panels, purple colours stand for the models with the redshift free ($z \neq z_{\text{host}}$) and green colours are the models with the redshift fixed to the host ($z = z_{\text{host}}$). Top right: Best-fitting spectra of MATLAS-2019 in both configurations. The red points are the observed magnitudes of the galaxy in the g, r, i, z, W1, W2, W3, and W4 bands, respectively. The yellow curves show the DECaLS transmission curves. The blue curves are the WISE transmission curves. The remaining panels show the MCMC corner plot comparing the posterior distribution of the two configurations. The first panel in each column shows the 1D posterior distribution of the fitted parameter, while the remaining panels show the correlation between parameters. This image can be read and interpreted as a covariance matrix. Columns stand for redshift, stellar mass, metallicity, star formation time-scale, mass-weighted age, and interstellar diffuse dust extinction. Metallicity, star formation time-scale, age, and dust attenuation are similar in both configurations. However, the stellar mass is strongly affected by the freedom of the redshift, as these two parameters are highly degenerate.</figcaption></figure>

non-cosmological distances are very uncertain, we consider $\delta_z$ being different than zero as a confirmation that redshifts with this method cannot be trusted rather than an indication that MATLAS UDGs are not at the redshifts of their hosts.

Going forward, we treat the models with fixed redshift as the best representation of our galaxies. Thus, from now on, we will only

discuss the results of this configuration. We reinforce, none the less, that even if the galaxies are not at this assumed distance, changes in the stellar populations (except stellar mass) between models are small and well within their uncertainties (see Figs 2, 3, and 4). We carefully checked that the changes in stellar mass do not strongly affect any of our conclusions that follow. Specifically, we have checked that
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Age and metallicity of the 59 MATLAS UDGs recovered with PROSPECTOR using the two different configurations discussed in Section 3.2. In all panels, purple colours stand for the models with the redshift free ($z \neq z_{\text{host}}$) and green colours are the models with the redshift fixed to the host ($z = z_{\text{host}}$). The bottom-right panel shows the age–metallicity distribution of the 59 MATLAS UDGs in both configurations. Histograms at the top and on the left show the marginal distributions of the data displayed on the age–metallicity plane. We conclude that models with the redshift free (purple) recover on average younger ages and more metal-rich populations for the UDGs.</figcaption></figure>

the mass–metallicity bimodality found for the models with fixed redshift and thoroughly discussed in Section 5.1 is also present in the models where the redshift is a free parameter. We have also checked that the change in stellar mass is small and not enough to change the classification of the UDGs according to their positioning on the mass–metallicity plane (i.e. which MZR they follow).

## 4.2 Median stellar populations

The results of the models with the redshift fixed at the redshift of the host galaxy are fully presented in Table B3. The median and median

absolute deviations of the stellar populations of the whole sample of MATLAS UDGs are presented below. We find that the MATLAS UDGs have intermediate-to-old ages, with a median mass-weighted age of $t_M = 7.1 \pm 1.8$ Gyr. The galaxies display an average metal-poor population with [M/H] = -1.2 ± 0.2 dex. These stellar populations are equivalent to those found for UDGs with spectroscopy by FM23. We find an average $\tau$ of 1.6 ± 0.7 Gyr. Finally, the median interstellar diffuse dust attenuation coming from the SED fitting of the galaxies in our sample is $A_V = 0.12 \pm 0.07$ mag, consistent with the expectation that these galaxies should have little-to-no dust component and with the caveat that the non-zero dust content could be a small artificial
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Age and metallicity of the 59 MATLAS UDGs recovered with PROSPECTOR using the two different configurations discussed in Section 3.2. In all panels, purple colours stand for the models with the redshift free ($z \neq z_{\text{host}}$) and green colours are the models with the redshift fixed to the host ($z = z_{\text{host}}$). The bottom-right panel shows the age–metallicity distribution of the 59 MATLAS UDGs in both configurations. Histograms at the top and on the left show the marginal distributions of the data displayed on the age–metallicity plane. We conclude that models with the redshift free (purple) recover on average younger ages and more metal-rich populations for the UDGs.</figcaption></figure>

the mass–metallicity bimodality found for the models with fixed redshift and thoroughly discussed in Section 5.1 is also present in the models where the redshift is a free parameter. We have also checked that the change in stellar mass is small and not enough to change the classification of the UDGs according to their positioning on the mass–metallicity plane (i.e. which MZR they follow).

## 4.2 Median stellar populations

The results of the models with the redshift fixed at the redshift of the host galaxy are fully presented in Table B3. The median and median

absolute deviations of the stellar populations of the whole sample of MATLAS UDGs are presented below. We find that the MATLAS UDGs have intermediate-to-old ages, with a median mass-weighted age of $t_M = 7.1 \pm 1.8$ Gyr. The galaxies display an average metal-poor population with [M/H] = -1.2 ± 0.2 dex. These stellar populations are equivalent to those found for UDGs with spectroscopy by FM23. We find an average $\tau$ of 1.6 ± 0.7 Gyr. Finally, the median interstellar diffuse dust attenuation coming from the SED fitting of the galaxies in our sample is $A_V = 0.12 \pm 0.07$ mag, consistent with the expectation that these galaxies should have little-to-no dust component and with the caveat that the non-zero dust content could be a small artificial
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Stellar mass, star formation time-scale, and dust attenuation of the 59 MATLAS UDGs recovered with PROSPECTOR using the two different configurations discussed in Section 3.2. In all panels, purple colours stand for the models with the redshift free and green colours are the models with the redshift fixed to the host. We conclude that models with the redshift free (purple) result in slightly more massive UDGs and a much narrower distribution centred at $log(M_*/M_\odot) \sim 7.8$ than models with the redshift fixed. No significant difference is observed in the star formation time-scales and dust attenuation recovered in both configurations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Accuracy in recovering photometric redshifts as a function of redshift. $\delta_z$ measures the difference between the spectroscopic redshifts of the candidate hosts and the recovered photometric redshifts of the UDGs. Marginal histogram shows the distribution of $\delta_z$. Positive values mean that UDGs were found to be in the foreground of their hosts. Negative values mean that UDGs were found in the background of their hosts. From previous findings (Heesters et al. 2023), the expectation is that at least 75 per cent of the dwarf galaxies in the MATLAS survey are at the same redshifts as their hosts. Our results show a much larger portion of the UDGs at different redshifts (>70 per cent of UDGs with $|\delta_z|>0.0067$, equivalent to a difference in velocity of $|\delta_v|>200$ km s$^{-1}$), indicating that photometric redshifts are not reliable at non-cosmological distances.</figcaption></figure>

addition that compensates for inaccuracies in the SED models. It is interesting to notice that although our dust priors extend out to 4.3 mag, the highest $A_V$ value found is 1.17 mag, highlighting the importance of the inclusion of the WISE upper limits from the 12 and 22 $\mu$m $\mu$m bands to constrain the amount of dust in the galaxies.

This power of the near- and mid-infrared bands in constraining the dust amount found had been discussed previously both by Pandya et al. (2018) and B22. None the less, as expected, UDGs with higher amounts of dust show higher degrees of degeneracy between the dust extinction and the metallicity and age.

Trends between the recovered stellar populations and the galaxy morphologies, environments, scaling relations, and GC-richnesses are investigated in Section 5.1.

## 4.3 Comparison with the literature

In this section, we compare our SED fitting results with stellar populations recovered with spectroscopy in the literature for some of the MATLAS UDGs.

### (i) MATLAS-2019

MATLAS-2019, also known as NGC 5846_UDG1 (Forbes et al. 2021), has been heavily studied in the past years due to its unusually large GC population (Müller et al., 2021; Danieli et al. 2022, Haacke et al., in preparation). In terms of its stellar populations, MATLAS-2019 has been spectroscopically studied by Müller et al. (2020), Heesters et al. (2023), and FM23.

Müller et al. (2020) found, using VLT/MUSE data, that MATLAS-2019 has a metallicity of $[M/H] = [Fe/H] = -1.33_{-0.01}^{+0.19}$ dex, an age of $11.2_{-0.8}^{+1.8}$ Gyr, and a stellar mass of $3.6 \times 10^8 M_\odot$ ($\log(M_*/M_\odot) = 8.56$). Morphological parameters for MATLAS-2019 were also obtained by Müller et al. (2020), who found a Sérsic index of $n = 0.73 \pm 0.01$, an axis ratio of $b/a = 0.90 \pm 0.01$, and an effective radius of $R_e = 17.2$ arcsec ($\sim 2.2$ kpc at the assumed distance of 26.3 Mpc). Heesters et al. (2023) studied MATLAS-2019 using the exact same VLT/MUSE data as Müller et al. (2020). They found, none the less, much older and more metal-poor populations, i.e. $[M/H] = [Fe/H] = -1.88_{-0.06}^{+0.13}$ dex and an age of $13.5_{-0.2}^{+0.5}$ Gyr. A reason for this discrepancy is not discussed by Heesters et al. (2023). Finally, FM23 has studied MATLAS-2019 using Keck/KCWI data. They found a metallicity of $[M/H] = -1.15 \pm 0.25$ dex, an age of $8.2 \pm 3.05$ Gyr and a stellar mass of $\log(M_*/M_\odot) = 8.1$, assuming a distance of 24.89 Mpc to the UDG (Forbes et al. 2021). The cause
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Stellar mass, star formation time-scale, and dust attenuation of the 59 MATLAS UDGs recovered with PROSPECTOR using the two different configurations discussed in Section 3.2. In all panels, purple colours stand for the models with the redshift free and green colours are the models with the redshift fixed to the host. We conclude that models with the redshift free (purple) result in slightly more massive UDGs and a much narrower distribution centred at $log(M_*/M_\odot) \sim 7.8$ than models with the redshift fixed. No significant difference is observed in the star formation time-scales and dust attenuation recovered in both configurations.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Accuracy in recovering photometric redshifts as a function of redshift. $\delta_z$ measures the difference between the spectroscopic redshifts of the candidate hosts and the recovered photometric redshifts of the UDGs. Marginal histogram shows the distribution of $\delta_z$. Positive values mean that UDGs were found to be in the foreground of their hosts. Negative values mean that UDGs were found in the background of their hosts. From previous findings (Heesters et al. 2023), the expectation is that at least 75 per cent of the dwarf galaxies in the MATLAS survey are at the same redshifts as their hosts. Our results show a much larger portion of the UDGs at different redshifts ($>70$ per cent of UDGs with $|\delta_z|>0.0067$, equivalent to a difference in velocity of $|\delta_V|>200$ km s$^{-1}$), indicating that photometric redshifts are not reliable at non-cosmological distances.</figcaption></figure>

addition that compensates for inaccuracies in the SED models. It is interesting to notice that although our dust priors extend out to 4.3 mag, the highest $A_V$ value found is 1.17 mag, highlighting the importance of the inclusion of the WISE upper limits from the 12 and 22 $\mu$m $\mu$m bands to constrain the amount of dust in the galaxies.

This power of the near- and mid-infrared bands in constraining the dust amount found had been discussed previously both by Pandya et al. (2018) and B22. None the less, as expected, UDGs with higher amounts of dust show higher degrees of degeneracy between the dust extinction and the metallicity and age.

Trends between the recovered stellar populations and the galaxy morphologies, environments, scaling relations, and GC-richnesses are investigated in Section 5.1.

## 4.3 Comparison with the literature

In this section, we compare our SED fitting results with stellar populations recovered with spectroscopy in the literature for some of the MATLAS UDGs.

### (i) MATLAS-2019

MATLAS-2019, also known as NGC 5846 UDG1 (Forbes et al. 2021), has been heavily studied in the past years due to its unusually large GC population (Müller et al., 2021; Danieli et al. 2022, Haacke et al., in preparation). In terms of its stellar populations, MATLAS-2019 has been spectroscopically studied by Müller et al. (2020), Heesters et al. (2023), and FM23.

Müller et al. (2020) found, using VLT/MUSE data, that MATLAS-2019 has a metallicity of [M/H] = [Fe/H] = $-1.33_{-0.01}^{+0.19}$ dex, an age of $11.2_{-0.8}^{+1.8}$ Gyr, and a stellar mass of $3.6 \times 10^8 M_\odot$ (log($M_*/M_\odot$) = 8.56). Morphological parameters for MATLAS-2019 were also obtained by Müller et al. (2020), who found a Sérsic index of $n = 0.73 \pm 0.01$, an axis ratio of $b/a = 0.90 \pm 0.01$, and an effective radius of $R_e = 17.2$ arcsec ($\sim 2.2$ kpc at the assumed distance of 26.3 Mpc). Heesters et al. (2023) studied MATLAS-2019 using the exact same VLT/MUSE data as Müller et al. (2020). They found, none the less, much older and more metal-poor populations, i.e. [M/H] = [Fe/H] = $-1.88_{-0.06}^{+0.13}$ dex and an age of $13.5_{-0.2}^{+0.5}$ Gyr. A reason for this discrepancy is not discussed by Heesters et al. (2023). Finally, FM23 has studied MATLAS-2019 using Keck/KCWI data. They found a metallicity of [M/H] = $-1.15 \pm 0.25$ dex, an age of $8.2 \pm 3.05$ Gyr and a stellar mass of log($M_*/M_\odot$) = 8.1, assuming a distance of 24.89 Mpc to the UDG (Forbes et al. 2021). The cause
for this difference is yet to be understood; however, we note that the S/N of the spectra used by FM23 is higher (S/N ~ 15–20 pix⁻¹) than that used by both Müller et al. (2020) and Heesters et al. (2023) (S/N = 12.4 pix⁻¹).

With PROSPECT, we recover a metallicity of [M/H] = −1.40<sup>+0.10</sup><sub>−0.12</sub> dex, a mass-weighted age of 11.22<sup>+1.79</sup><sub>−3.23</sub> Gyr and a stellar mass of log(M<sub>⋆</sub>/M<sub>☉</sub>) = 8.01<sup>+0.03</sup><sub>−0.05</sub> for MATLAS-2019, for a distance of 25.2 Mpc (i.e. the distance to the host NGC 5846). Our study is consistent within the uncertainties with both Müller et al. (2020) and FM23. The results, none the less, in terms of age and metallicity, are closer to the ones found by Müller et al. (2020), indicating a prevalence of a metal-poor and old population. We found MATLAS-2019 to be younger and more metal-rich than the results obtained by Heesters et al. (2023). Additionally, the morphological parameters obtained with GALFITM for MATLAS-2019 were n = 0.6 ± 0.1, b/a = 0.97 ± 0.03, and R<sub>e</sub> = 15.7 ± 0.4 arcsec. These are all consistent within the uncertainties with the findings of Müller et al. (2020).

(ii) MATLAS-585 and MATLAS-2103

These two MATLAS UDGs were recently studied using VLT/MUSE data by Heesters et al. (2023). MATLAS-585 was found to have a metallicity of [M/H] = [Fe/H] = −1.88<sup>+0.17</sup><sub>−0.13</sub> dex and an age of 9.0<sup>+2.9</sup><sub>−3.1</sub> Gyr. With PROSPECT, we find that MATLAS-585 has a metallicity of [M/H] = −1.36<sup>+0.27</sup><sub>−0.16</sub> dex and an age of 8.91<sup>+3.44</sup><sub>−4.31</sub> Gyr. Thus, our results in terms of age are consistent with those of Heesters et al. (2023), but they find a much more metal-poor population for MATLAS-585 than we do with PROSPECT.

MATLAS-2103 was found by Heesters et al. (2023) to have a metallicity of [M/H] = [Fe/H] = −1.69<sup>+0.01</sup><sub>−0.22</sub> dex and an age of 11.3<sup>+2.2</sup><sub>−2.5</sub> Gyr. In this study, we find that MATLAS-2103 has a metallicity of [M/H] = −1.40<sup>+0.14</sup><sub>−0.13</sub> dex and an age of 9.79<sup>+2.81</sup><sub>−3.84</sub> Gyr. Similar to the case of MATLAS-585, we find ages consistent with those of Heesters et al. (2023), but their stellar populations are more metal-poor than ours. Notably, Heesters et al. (2023) tend to find systematically more metal-poor stellar populations for galaxies if we take this study, FM23 and Müller et al. (2020) as baselines.

Overall, we show that our results are consistent with those found using spectroscopy. Our SED fitting technique had already been shown to provide consistent results with spectroscopy in the studies of B22 and FM23. In both cases our results were shown to have a metallicity offset of ~0.2 dex, which is well within our uncertainties. We reinforce that while SED fitting results have much larger uncertainties, stellar populations obtained from broad-band photometry are far less restricted than spectroscopic observations in depth and spatial coverage. Additionally, the ability of this technique to gather information from a wide wavelength range makes up for the loss of detailed features provided by spectroscopic data. SED fitting also allows for complete sample studies such as this one, which are prohibitively expensive to do spectroscopically on big telescopes. Apart from the larger uncertainties, another downside of SED fitting is the difficulty in obtaining [α/Fe] and dynamical information from the galaxies. Thus, spectral and photometric information are complementary and should be jointly used to better understand UDGs.

5 DISCUSSION

As mentioned in the Introduction, the stellar populations of UDGs
can help understanding their formation histories. This can be made
by means of comparing these populations to those of different classes
of galaxies to try and identify similarities/differences between them.
In Section 5.1, we use scaling relations (i.e. the MZR) to compare

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

**Figure 6.** Stellar mass–metallicity distribution of UDGs. **Top panel:** MATLAS UDGs are shown in red. Field UDGs from Barbosa et al. (2020) are shown in blue. UDGs in different environments from the spectroscopic sample of FM23 are shown in green. UDGs in a variety of environments from B22 are shown in magenta. The Simon (2019) MZR for classical Local Group dwarf galaxies is shown with the black dashed line. The dash-dotted line is the evolving MZR at redshift z = 2.2 from Ma et al. (2016). **Middle panel:** Difference in metallicity between the MATLAS UDGs and the Simon (2019) classical dwarf MZR. UDGs at δ<sub>MZR</sub> = 0 follow perfectly the MZR, if δ<sub>MZR</sub> > 0, then the UDG lies above the MZR, if δ<sub>MZR</sub> < 0, then the UDG lies below the MZR. **Bottom panel:** Difference in metallicity between the MATLAS UDGs and the Ma et al. (2016) MZR at z = 2.2. MATLAS UDGs show a bimodality in the mass–metallicity distribution, with the first mode being consistent with the classical dwarf MZR and the second being more consistent with the high-redshift MZR.

UDGs to classical dwarfs (observations) and high-redshift galaxies (simulations) to try and elucidate their origins. In Section 5.2 we discuss if there is evidence for two classes of UDG based on their stellar populations and morphological properties and associate these classes with their most likely formation scenarios.

5.1 MZR: clues to the origins of UDGs

We explore in this section the positioning of MATLAS UDGs on the stellar mass–metallicity plane, especially in comparison with known MZRs. We investigate trends with different morphological and stellar population parameters, as well as GC-richness.

For reference, in Figs 6 and 7, we plot two MZRs: Simon (2019) for classical dwarfs, and the simulated MZR from Ma et al. (2016) for evolving sources at a redshift of z = 2.2. To be consistent, we plot total metallicities from PROSPECT on top of the converted to [M/H] Simon (2019) relation, from its original in [Fe/H] (Vazdekis et al. 2015). A full discussion on how to make this conversion and
for this difference is yet to be understood; however, we note that the S/N of the spectra used by FM23 is higher (S/N ~ 15–20 pix⁻¹) than that used by both Müller et al. (2020) and Heesters et al. (2023) (S/N = 12.4 pix⁻¹).

With PROSPECT, we recover a metallicity of [M/H] = −1.40<sup>+0.10</sup><sub>−0.12</sub> dex, a mass-weighted age of 11.22<sup>+1.79</sup><sub>−3.23</sub> Gyr and a stellar mass of log(M<sub>⋆</sub>/M<sub>☉</sub>) = 8.01<sup>+0.03</sup><sub>−0.05</sub> for MATLAS-2019, for a distance of 25.2 Mpc (i.e. the distance to the host NGC 5846). Our study is consistent within the uncertainties with both Müller et al. (2020) and FM23. The results, none the less, in terms of age and metallicity, are closer to the ones found by Müller et al. (2020), indicating a prevalence of a metal-poor and old population. We found MATLAS-2019 to be younger and more metal-rich than the results obtained by Heesters et al. (2023). Additionally, the morphological parameters obtained with GALFITM for MATLAS-2019 were n = 0.6 ± 0.1, b/a = 0.97 ± 0.03, and R<sub>e</sub> = 15.7 ± 0.4 arcsec. These are all consistent within the uncertainties with the findings of Müller et al. (2020).

(ii) MATLAS-585 and MATLAS-2103

These two MATLAS UDGs were recently studied using VLT/MUSE data by Heesters et al. (2023). MATLAS-585 was found to have a metallicity of [M/H] = [Fe/H] = −1.88<sup>+0.17</sup><sub>−0.13</sub> dex and an age of 9.0<sup>+2.9</sup><sub>−3.1</sub> Gyr. With PROSPECT, we find that MATLAS-585 has a metallicity of [M/H] = −1.36<sup>+0.27</sup><sub>−0.16</sub> dex and an age of 8.91<sup>+3.44</sup><sub>−4.31</sub> Gyr. Thus, our results in terms of age are consistent with those of Heesters et al. (2023), but they find a much more metal-poor population for MATLAS-585 than we do with PROSPECT.

MATLAS-2103 was found by Heesters et al. (2023) to have a metallicity of [M/H] = [Fe/H] = −1.69<sup>+0.01</sup><sub>−0.22</sub> dex and an age of 11.3<sup>+2.2</sup><sub>−2.5</sub> Gyr. In this study, we find that MATLAS-2103 has a metallicity of [M/H] = −1.40<sup>+0.14</sup><sub>−0.13</sub> dex and an age of 9.79<sup>+2.81</sup><sub>−3.84</sub> Gyr. Similar to the case of MATLAS-585, we find ages consistent with those of Heesters et al. (2023), but their stellar populations are more metal-poor than ours. Notably, Heesters et al. (2023) tend to find systematically more metal-poor stellar populations for galaxies if we take this study, FM23 and Müller et al. (2020) as baselines.

Overall, we show that our results are consistent with those found using spectroscopy. Our SED fitting technique had already been shown to provide consistent results with spectroscopy in the studies of B22 and FM23. In both cases our results were shown to have a metallicity offset of ~0.2 dex, which is well within our uncertainties. We reinforce that while SED fitting results have much larger uncertainties, stellar populations obtained from broad-band photometry are far less restricted than spectroscopic observations in depth and spatial coverage. Additionally, the ability of this technique to gather information from a wide wavelength range makes up for the loss of detailed features provided by spectroscopic data. SED fitting also allows for complete sample studies such as this one, which are prohibitively expensive to do spectroscopically on big telescopes. Apart from the larger uncertainties, another downside of SED fitting is the difficulty in obtaining [α/Fe] and dynamical information from the galaxies. Thus, spectral and photometric information are complementary and should be jointly used to better understand UDGs.

5 DISCUSSION

As mentioned in the Introduction, the stellar populations of UDGs
can help understanding their formation histories. This can be made
by means of comparing these populations to those of different classes
of galaxies to try and identify similarities/differences between them.
In Section 5.1, we use scaling relations (i.e. the MZR) to compare

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

**Figure 6.** Stellar mass–metallicity distribution of UDGs. **Top panel:** MATLAS UDGs are shown in red. Field UDGs from Barbosa et al. (2020) are shown in blue. UDGs in different environments from the spectroscopic sample of FM23 are shown in green. UDGs in a variety of environments from B22 are shown in magenta. The Simon (2019) MZR for classical Local Group dwarf galaxies is shown with the black dashed line. The dash-dotted line is the evolving MZR at redshift z = 2.2 from Ma et al. (2016). **Middle panel:** Difference in metallicity between the MATLAS UDGs and the Simon (2019) classical dwarf MZR. UDGs at δ<sub>MZR</sub> = 0 follow perfectly the MZR, if δ<sub>MZR</sub> > 0, then the UDG lies above the MZR, if δ<sub>MZR</sub> < 0, then the UDG lies below the MZR. **Bottom panel:** Difference in metallicity between the MATLAS UDGs and the Ma et al. (2016) MZR at z = 2.2. MATLAS UDGs show a bimodality in the mass–metallicity distribution, with the first mode being consistent with the classical dwarf MZR and the second being more consistent with the high-redshift MZR.

UDGs to classical dwarfs (observations) and high-redshift galaxies (simulations) to try and elucidate their origins. In Section 5.2 we discuss if there is evidence for two classes of UDG based on their stellar populations and morphological properties and associate these classes with their most likely formation scenarios.

5.1 MZR: clues to the origins of UDGs

We explore in this section the positioning of MATLAS UDGs on the stellar mass–metallicity plane, especially in comparison with known MZRs. We investigate trends with different morphological and stellar population parameters, as well as GC-richness.

For reference, in Figs 6 and 7, we plot two MZRs: Simon (2019) for classical dwarfs, and the simulated MZR from Ma et al. (2016) for evolving sources at a redshift of z = 2.2. To be consistent, we plot total metallicities from PROSPECT on top of the converted to [M/H] Simon (2019) relation, from its original in [Fe/H] (Vazdekis et al. 2015). A full discussion on how to make this conversion and
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Distribution of the physical properties of the MATLAS UDGs in the stellar mass–metallicity plane. The plus signs show the results obtained from SED fitting with PROSPECTOR for the 59 UDGs in our sample, with the two confirmed tidal UDGs in our sample being highlighted with star symbols (MATLAS-1824 and MATLAS-478). Top row: UDGs are colour-coded by surface brightness, effective radii, axis ratio, and local environment (Marleau et al. 2021), respectively. Bottom row: UDGs are colour-coded by their mass-weighted ages, star formation time-scales, number of GCs (Marleau et al., submitted), and GC specific frequencies (per unit stellar mass), respectively. The Simon (2019) MZR for classical Local Group dwarf galaxies is shown with the dashed black line. The dash-dotted black line is the evolving MZR at redshift $z = 2.2$ from Ma et al. (2016). An alternative view of this figure is shown in Appendix D1. These results indicate a possible bimodality in the MATLAS UDGs population: Population (i) UDGs are fainter, smaller, elongated, live in less dense environments, are younger, have prolonged SFHs, have a smaller number of GCs and on average smaller specific frequencies. Population (ii) UDGs in this population are brighter, bigger, rounder, live in denser environments, are older, have shorter SFHs, host on average more GCs and have on average larger specific frequencies.</figcaption></figure>

the caveats that arise in doing so are given in B22. In Fig. 6 we show the distribution of the MATLAS UDGs in the mass–metallicity plane in comparison to three other works in the literature. In blue we show the study of 100 field UDGs using SED fitting from Barbosa et al. (2020). The study of 29 UDGs in a variety of environments also using SED fitting from B22 is shown in magenta. Lastly, the largest spectroscopic study of UDGs to date (25) from FM23 is shown in green. The UDGs in FM23 are in multiple environments, but biased towards clusters.

The MATLAS UDGs, at first glance, seem to form a new MZR, with a metallicity that is anticorrelated with the stellar mass. However, the population shows a bimodality, indicating that instead of an anticorrelation, the MATLAS UDGs split into two samples. One of them is well within the classical dwarf MZR and another has lower metallicities and lies closer to the simulated high-$z$ MZR. To make this bimodality clearer, in the bottom panels of Fig. 6, we show the residual MZR (i.e. the difference between the metallicities of the UDG population and the value expected from the MZR). It can be seen that the low-mass metal-rich UDGs are well within the scatter of the classical dwarf MZR, while the high-mass metal-poor UDGs are further down and in fact are more consistent with the MZR at high redshift. This bimodality may be indicative of different formation scenarios for the UDGs populating each mode.

When comparing to the works in the literature, we see that while all have a population of UDGs that follows well the classical dwarf

MZR, they all seem to probe different regions of the parameter space. Barbosa et al. (2020) and the MATLAS UDGs are biased towards lower mass UDGs in low-density environments, while B22 and FM23 are biased towards more massive cluster UDGs. FM23, in addition, probes a population of more massive, metal-rich, cluster UDGs, although overall their distribution scatters around the MZR relations at $log(M_*/M_\odot) \sim 8$.

Barbosa et al. (2020) find many UDGs to be more metal-rich than the classical dwarf MZR. These were found to be the youngest in their sample. FM23 found a similar subpopulation of UDGs more metal-rich than the classical dwarf MZR, as well as another subpopulation that is more metal-poor. We recover the population of metal-poor UDGs found by them, but no UDGs in the same mass range that are metal-rich. The explanation for this is that the sample of FM23 is biased towards the cluster environment where tidal effects are more common (Sales, Wetzel & Fattahi 2022) and result in the higher metallicities found.

Finally, the sample covered by B22, although more representative than the current sample in terms of environments, is more constrained in terms of stellar masses, showing only a small range around $8.0 \le \log(M_*/M_\odot) \le 8.5$. Similar to the findings of FM23 and the ones of this study, they found that a subsample of the UDGs follows the classical dwarf MZR, while another population is much more metal-poor and is more consistent with the high-$z$ MZR.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Distribution of the physical properties of the MATLAS UDGs in the stellar mass–metallicity plane. The plus signs show the results obtained from SED fitting with PROSPECTOR for the 59 UDGs in our sample, with the two confirmed tidal UDGs in our sample being highlighted with star symbols (MATLAS-1824 and MATLAS-478). Top row: UDGs are colour-coded by surface brightness, effective radii, axis ratio, and local environment (Marleau et al. 2021), respectively. Bottom row: UDGs are colour-coded by their mass-weighted ages, star formation time-scales, number of GCs (Marleau et al., submitted), and GC specific frequencies (per unit stellar mass), respectively. The Simon (2019) MZR for classical Local Group dwarf galaxies is shown with the dashed black line. The dash-dotted black line is the evolving MZR at redshift $z = 2.2$ from Ma et al. (2016). An alternative view of this figure is shown in Appendix D1. These results indicate a possible bimodality in the MATLAS UDGs population: Population (i) UDGs are fainter, smaller, elongated, live in less dense environments, are younger, have prolonged SFHs, have a smaller number of GCs and on average smaller specific frequencies. Population (ii) UDGs in this population are brighter, bigger, rounder, live in denser environments, are older, have shorter SFHs, host on average more GCs and have on average larger specific frequencies.</figcaption></figure>

the caveats that arise in doing so are given in B22. In Fig. 6 we show the distribution of the MATLAS UDGs in the mass–metallicity plane in comparison to three other works in the literature. In blue we show the study of 100 field UDGs using SED fitting from Barbosa et al. (2020). The study of 29 UDGs in a variety of environments also using SED fitting from B22 is shown in magenta. Lastly, the largest spectroscopic study of UDGs to date (25) from FM23 is shown in green. The UDGs in FM23 are in multiple environments, but biased towards clusters.

The MATLAS UDGs, at first glance, seem to form a new MZR, with a metallicity that is anticorrelated with the stellar mass. However, the population shows a bimodality, indicating that instead of an anticorrelation, the MATLAS UDGs split into two samples. One of them is well within the classical dwarf MZR and another has lower metallicities and lies closer to the simulated high-$z$ MZR. To make this bimodality clearer, in the bottom panels of Fig. 6, we show the residual MZR (i.e. the difference between the metallicities of the UDG population and the value expected from the MZR). It can be seen that the low-mass metal-rich UDGs are well within the scatter of the classical dwarf MZR, while the high-mass metal-poor UDGs are further down and in fact are more consistent with the MZR at high redshift. This bimodality may be indicative of different formation scenarios for the UDGs populating each mode.

When comparing to the works in the literature, we see that while all have a population of UDGs that follows well the classical dwarf

MZR, they all seem to probe different regions of the parameter space. Barbosa et al. (2020) and the MATLAS UDGs are biased towards lower mass UDGs in low-density environments, while B22 and FM23 are biased towards more massive cluster UDGs. FM23, in addition, probes a population of more massive, metal-rich, cluster UDGs, although overall their distribution scatters around the MZR relations at $log(M_*/M_\odot) \sim 8$.

Barbosa et al. (2020) find many UDGs to be more metal-rich than the classical dwarf MZR. These were found to be the youngest in their sample. FM23 found a similar subpopulation of UDGs more metal-rich than the classical dwarf MZR, as well as another subpopulation that is more metal-poor. We recover the population of metal-poor UDGs found by them, but no UDGs in the same mass range that are metal-rich. The explanation for this is that the sample of FM23 is biased towards the cluster environment where tidal effects are more common (Sales, Wetzel & Fattahi 2022) and result in the higher metallicities found.

Finally, the sample covered by B22, although more representative than the current sample in terms of environments, is more constrained in terms of stellar masses, showing only a small range around $8.0 \le \log(M_*/M_\odot) \le 8.5$. Similar to the findings of FM23 and the ones of this study, they found that a subsample of the UDGs follows the classical dwarf MZR, while another population is much more metal-poor and is more consistent with the high-$z$ MZR.
### 5.1.1 Dependence with morphology and local environment

The fact that the morphology of galaxies is correlated with the environment they reside in has been known since the 1970's, giving origin to the so-called 'Morphology-Density Relation' (Oemler 1974; Davis & Geller 1976; Dressler 1980; Dressler et al. 1997), first proposed by Edwin Hubble in the 1930s (Hubble & Humason 1931). In a very simplistic explanation, ETGs are found to occupy the densest environments such as the centres of clusters of galaxies, while late-type galaxies are observed in their outskirts and to be prevalent in less dense environments, such as the field. The morphology-density relation was also found to hold for classical dwarf galaxies (Phillipps et al. 1998; Penny et al. 2011), showing that they behave in a similar way to their massive counterparts in terms of the environments they inhabit. Noteworthy, Poulain et al. (2021) also found a morphology-density relation for the dwarfs on the MATLAS fields comparable to that seen in dense clusters, and Li et al. (2023) found a similar behaviour for UDGs surrounding Milky Way analogues.

Various works in the literature have looked into the morphological parameters of UDGs (e.g. Kado-Fong et al. 2020, 2021; Rong et al. 2020; Van Nest et al. 2022; Li et al. 2023), such as their ellipticities, which can hold much information about the physical processes involved in forming UDGs. Rong et al. (2020), for example, found that the morphologies of cluster UDGs depend on their luminosities, environments, and redshifts. In terms of axis ratios, they showed that there are at least two populations of UDGs, the elongated (b/a ~ 0.4) and the round ones (b/a ~ 0.9). In their study, UDGs at higher redshifts (z ~ 0.3) are more elongated than their low-redshift counterparts, indicating a possible evolutionary path of UDGs becoming rounder with time. They also show that bright ($M_g$ > -15.2 mag) UDGs are rounder than the fainter ones ($M_g$ < -15.2 mag). Finally, they showed that UDGs closer to the centre of clusters are rounder than the ones in the outskirts or outside the cluster virial radius. Their findings apply to both nucleated and non-nucleated UDGs, although it is worth noting that the number of nucleated UDGs increases towards the centre of clusters and these are the roundest sources. A bimodality in the axis ratio, mass, and luminosity of UDGs has also been reported in the Coma and Virgo clusters by Lim et al. (2018) and Lim et al. (2020), respectively. They find that UDGs that are dark-matter-dominated (M/L > 1000) have relatively rounder shapes (higher b/a), while UDGs with M/L ~ 500 are more elongated (low b/a). They also found that the round UDGs have higher GC specific frequencies than the more elongated ones.

In the top row of Fig. 7, we show the MATLAS UDGs colour-coded by their surface brightnesses, effective radii, axis ratio, and local volume density. Similar to the findings of Rong et al. (2020), we also find a correlation between the surface brightness and the mass of the galaxies. Less massive UDGs are, on average, the faintest ones, while more massive UDGs are the brightest. This is expected given the narrow range in surface brightness and size, and the fact that the Sérsic index of the galaxies is close to 1. In addition, a similar correlation exists for the effective radii of the galaxies. Less massive UDGs are smaller, while the more massive ones are bigger. This behaviour is not unique to UDGs, it is in fact expected from the mass-size relation, but it is reassuring to confirm this expectation in our sample of galaxies.

In terms of their axis ratios, we can see a weak correlation with both the mass and metallicity, where more elongated UDGs seem to be positioned in a less massive/ more metal-rich region, while the roundest UDGs occupy a position corresponding to more

massive, more metal-poor populations. There is a non-negligible portion of UDGs also with round morphologies located in the upper parts of the plane. One explanation for this is that the ellipticity of galaxies is highly dependent on their inclinations, which could drive oblate/flattened galaxies to look round if they are face-on. It is also worth noticing that more massive edge-on UDGs would likely be pushed over the surface brightness threshold definition for UDGs due to their higher apparent brightness, which can explain why we see such a correlation between the surface brightness, axis ratio and mass of the galaxies. Additionally, less massive galaxies tend to be rounder (e.g. Sánchez-Janssen et al. 2019), however, these faint galaxies are below our mass range, also explaining why we do not have a population of round galaxies at lower masses.

We note the study of Cardona-Barrero et al. (2020), who used the NIHAO simulations to suggest a correlation between the morphology and kinematics of UDGs. They suggested that more elongated UDGs (i.e. smaller axis ratios) are rotationally supported, while the rounder ones (higher axis ratios) are pressure-dominated. These findings are further supported for classical dwarfs by Pfeffer et al. (submitted).

If we now compare these findings with the panel where the UDGs are coloured according to the density of the environment they reside in (i.e. smaller values of $\log \rho_{10}$ are for galaxies that are fairly isolated, while higher values signify that galaxies are in more dense environments), we can start to see some forms of the morphology-density relation appearing. That is, UDGs that show the most elongated morphologies are the ones located in the less dense environments. Alternatively, UDGs that are rounder tend to be in the most dense environments. If we consider that UDGs in denser environments are prone to more mergers and encounters than their isolated counterparts, then having rounder morphologies is a natural consequence of that (Moore et al. 1996, 1999; Moore, Lake & Katz 1998). We caveat that this interpretation is done for the full sample of galaxies, and individual cases may differ from the ensemble properties.

If we analyse the environment panel on its own, there appears to be a weak trend between the environment and the positioning of the UDGs, where the majority of the UDGs in the most dense of our probed environments are below the classical dwarf MZR. We remind the reader that all of the MATLAS UDGs are, by definition, in low-to-moderate density environments, and thus, the range of densities probed is limited and does not include the most dense environments, such as clusters. Selection effects can be playing a role in our findings and, thus, comparisons with UDGs located in higher density environments are necessary and will be pursued in a follow-up study (Buzzo et al., in preparation). It is interesting, none the less, to see that even with such a small range of densities, a correlation between the environment and the stellar populations can still be weakly identified. Although a correlation was identified with the local volume density, we found no correlation between any properties and the projected groupcentric distance of the UDGs. We emphasize that these are the interpretations for the ensemble of MATLAS UDGs, but as can be seen, individual exceptions may exist (e.g. round UDGs at low-density environments or elongated ones at the highest probed densities). We use colour-coding to show the changes in the properties of the MATLAS UDGs across the mass-metallicity plane in Fig. 7, but for an alternative view, we provide in Fig. D1 direct plots of each property against their positioning on the plane to emphasize the discussed trends.
### 5.1.1 Dependence with morphology and local environment

The fact that the morphology of galaxies is correlated with the environment they reside in has been known since the 1970's, giving origin to the so-called 'Morphology-Density Relation' (Oemler 1974; Davis & Geller 1976; Dressler 1980; Dressler et al. 1997), first proposed by Edwin Hubble in the 1930s (Hubble & Humason 1931). In a very simplistic explanation, ETGs are found to occupy the densest environments such as the centres of clusters of galaxies, while late-type galaxies are observed in their outskirts and to be prevalent in less dense environments, such as the field. The morphology-density relation was also found to hold for classical dwarf galaxies (Phillipps et al. 1998; Penny et al. 2011), showing that they behave in a similar way to their massive counterparts in terms of the environments they inhabit. Noteworthy, Poulain et al. (2021) also found a morphology-density relation for the dwarfs on the MATLAS fields comparable to that seen in dense clusters, and Li et al. (2023) found a similar behaviour for UDGs surrounding Milky Way analogues.

Various works in the literature have looked into the morphological parameters of UDGs (e.g. Kado-Fong et al. 2020, 2021; Rong et al. 2020; Van Nest et al. 2022; Li et al. 2023), such as their ellipticities, which can hold much information about the physical processes involved in forming UDGs. Rong et al. (2020), for example, found that the morphologies of cluster UDGs depend on their luminosities, environments, and redshifts. In terms of axis ratios, they showed that there are at least two populations of UDGs, the elongated ($b/a \sim 0.4$) and the round ones ($b/a \sim 0.9$). In their study, UDGs at higher redshifts ($z \sim 0.3$) are more elongated than their low-redshift counterparts, indicating a possible evolutionary path of UDGs becoming rounder with time. They also show that bright ($M_g > -15.2$ mag) UDGs are rounder than the fainter ones ($M_g < -15.2$ mag). Finally, they showed that UDGs closer to the centre of clusters are rounder than the ones in the outskirts or outside the cluster virial radius. Their findings apply to both nucleated and non-nucleated UDGs, although it is worth noting that the number of nucleated UDGs increases towards the centre of clusters and these are the roundest sources. A bimodality in the axis ratio, mass, and luminosity of UDGs has also been reported in the Coma and Virgo clusters by Lim et al. (2018) and Lim et al. (2020), respectively. They find that UDGs that are dark-matter-dominated (M/L > 1000) have relatively rounder shapes (higher $b/a$), while UDGs with M/L ~ 500 are more elongated (low $b/a$). They also found that the round UDGs have higher GC specific frequencies than the more elongated ones.

In the top row of Fig. 7, we show the MATLAS UDGs colour-coded by their surface brightnesses, effective radii, axis ratio, and local volume density. Similar to the findings of Rong et al. (2020), we also find a correlation between the surface brightness and the mass of the galaxies. Less massive UDGs are, on average, the faintest ones, while more massive UDGs are the brightest. This is expected given the narrow range in surface brightness and size, and the fact that the Sérsic index of the galaxies is close to 1. In addition, a similar correlation exists for the effective radii of the galaxies. Less massive UDGs are smaller, while the more massive ones are bigger. This behaviour is not unique to UDGs, it is in fact expected from the mass-size relation, but it is reassuring to confirm this expectation in our sample of galaxies.

In terms of their axis ratios, we can see a weak correlation with both the mass and metallicity, where more elongated UDGs seem to be positioned in a less massive/ more metal-rich region, while the roundest UDGs occupy a position corresponding to more

massive, more metal-poor populations. There is a non-negligible portion of UDGs also with round morphologies located in the upper parts of the plane. One explanation for this is that the ellipticity of galaxies is highly dependent on their inclinations, which could drive oblate/flattened galaxies to look round if they are face-on. It is also worth noticing that more massive edge-on UDGs would likely be pushed over the surface brightness threshold definition for UDGs due to their higher apparent brightness, which can explain why we see such a correlation between the surface brightness, axis ratio and mass of the galaxies. Additionally, less massive galaxies tend to be rounder (e.g. Sánchez-Janssen et al. 2019), however, these faint galaxies are below our mass range, also explaining why we do not have a population of round galaxies at lower masses.

We note the study of Cardona-Barrero et al. (2020), who used the NIHAO simulations to suggest a correlation between the morphology and kinematics of UDGs. They suggested that more elongated UDGs (i.e. smaller axis ratios) are rotationally supported, while the rounder ones (higher axis ratios) are pressure-dominated. These findings are further supported for classical dwarfs by Pfeffer et al. (submitted).

If we now compare these findings with the panel where the UDGs are coloured according to the density of the environment they reside in (i.e. smaller values of $\log \rho_{10}$ are for galaxies that are fairly isolated, while higher values signify that galaxies are in more dense environments), we can start to see some forms of the morphology-density relation appearing. That is, UDGs that show the most elongated morphologies are the ones located in the less dense environments. Alternatively, UDGs that are rounder tend to be in the most dense environments. If we consider that UDGs in denser environments are prone to more mergers and encounters than their isolated counterparts, then having rounder morphologies is a natural consequence of that (Moore et al. 1996, 1999; Moore, Lake & Katz 1998). We caveat that this interpretation is done for the full sample of galaxies, and individual cases may differ from the ensemble properties.

If we analyse the environment panel on its own, there appears to be a weak trend between the environment and the positioning of the UDGs, where the majority of the UDGs in the most dense of our probed environments are below the classical dwarf MZR. We remind the reader that all of the MATLAS UDGs are, by definition, in low-to-moderate density environments, and thus, the range of densities probed is limited and does not include the most dense environments, such as clusters. Selection effects can be playing a role in our findings and, thus, comparisons with UDGs located in higher density environments are necessary and will be pursued in a follow-up study (Buzzo et al., in preparation). It is interesting, none the less, to see that even with such a small range of densities, a correlation between the environment and the stellar populations can still be weakly identified. Although a correlation was identified with the local volume density, we found no correlation between any properties and the projected groupcentric distance of the UDGs. We emphasize that these are the interpretations for the ensemble of MATLAS UDGs, but as can be seen, individual exceptions may exist (e.g. round UDGs at low-density environments or elongated ones at the highest probed densities). We use colour-coding to show the changes in the properties of the MATLAS UDGs across the mass-metallicity plane in Fig. 7, but for an alternative view, we provide in Fig. D1 direct plots of each property against their positioning on the plane to emphasize the discussed trends.
### 5.1.2 Dependence with age and star formation time-scale

We now turn to understanding how the positioning of the UDGs in the stellar mass–metallicity plane depends on their SFH properties, e.g. their ages and star formation time-scales.

Ferré-Mateu et al. (2018) have shown that there exists an age dependence with projected clustercentric distance, i.e. UDGs are younger at larger projected clustercentric radii. Alabi et al. (2018) and Kadowaki et al. (2021) also found a colour dependence with the environment, with bluer UDGs residing in lower density environments than the redder ones. Finally, B22 found, using SED fitting, that UDGs in the field are on average younger and more metal-rich than cluster UDGs. Group UDGs in their work were found to be in a transitional space in the age–metallicity plane. However, these were based on very small number statistics, hence the need for this study focused mainly on group UDGs. Similar to these recent findings, the ages of the MATLAS UDGs also seem to distill into two main types: younger UDGs following the dwarf MZR and older ones following the high-$z$ MZR. We however caution the reader that ages recovered with PROSPECTOR are very sensitive to the version of the code used (see Section 3.2), as well as to the uncertainties in the optical photometry (see Appendix C).

One interesting recent finding on the stellar populations of UDGs by FM23 is a correlation between the $\alpha$ enrichment ([α/Fe]) and the SFHs of UDGs. Using a phase-space diagram, which may be particularly related to the environment of the UDGs, they show that the UDGs with the highest [α/Fe] are in the region corresponding to early infallers, while UDGs with the lowest [α/Fe] tend to be in the outskirts of clusters or in less dense environments. The $\alpha$-enrichment can be associated with the speed at which star formation starts and quenches within galaxies. Thus, this parameter can very well separate early-quenched galaxies (high [α/Fe]) from galaxies with prolonged SFHs (low [α/Fe]).

With SED fitting we cannot recover [α/Fe]. We can, however, look into the star formation time-scale of galaxies, i.e. the rate of the exponential decay of star formation, to understand how fast SF was quenched within the galaxy. In the bottom row of Fig. 7, we colour-code our UDGs by their mass-weighted ages ($t_M$) and star formation time-scales ($\tau$). We can see that UDGs that have the fastest star formation episodes (i.e. smallest $\tau$) are in the lowest regions of the plane, coinciding with the positioning of the older galaxies in the age panel. They are consistent with the evolving MZR at high redshift, indicative of early-quenching, which is in agreement with the results from FM23, and suggestive of a failed galaxy-like origin. On the other hand, UDGs that show prolonged SFHs (i.e. large $\tau$) are the ones that are younger and more metal-rich, aligned with a puffed-up dwarf formation scenario, as also found by FM23.

### 5.1.3 Dependence with GC-richness and specific frequency

One of the properties of UDGs that draws a lot of attention is the large number of GCs that some of them host in spite of their low stellar masses. Some extreme examples include DF44 (van Dokkum et al. 2017, but see also Saifollahi et al. 2022 and Forbes & Gannon 2024) and NGC 5846. UDG1/MATLAS-2019 (Müller et al. 2020; Forbes et al. 2021; Danieli et al. 2022). Simulations and analytical models have tried to explain the existence of these many GCs in such low-mass galaxies (e.g. Carleton et al. 2021; Trujillo-Gomez et al. 2021; Danieli et al. 2022; van Dokkum et al. 2022; Doppel et al. 2024). In parallel, observational work has focused on understanding what are the main differences between GC-poor and GC-rich UDGs,

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Difference in metallicity ($\delta_{MZR}$) between known MZRs and the MATLAS UDGs colour-coded by their number of GCs. Description is the same as Fig. 6. GC-poor MATLAS UDGs follow the classical dwarf MZR from Simon (2019). On another hand, all UDGs that follow the Ma et al. (2016) MZR at high redshift are GC-rich.</figcaption></figure>

in an attempt to correlate the higher fraction of GCs with other galaxy properties (Forbes et al. 2020; La Marca et al. 2022).

In terms of stellar populations, B22 used SED fitting techniques to find evidence for an inherent difference between GC-rich and GC-poor UDGs. They found that GC-poor UDGs are consistent with the classical dwarf MZR from Simon (2019), having consistently higher metallicities than their GC-rich counterparts. Interestingly, all of the UDGs that are consistent with the evolving MZR at $z = 2.2$ proposed by Ma et al. (2016) are GC-rich, suggesting a formation scenario relying on early-quenching for these galaxies. Following up on this finding, FM23 studied the stellar populations of 25 UDGs, the largest spectroscopic sample to date, using data from Keck/KCWI. They found that GC-poor UDGs consistently follow the classical dwarf MZR. However, they found that the GC-rich UDGs can be found both above and below the Simon (2019) MZR. They suggested that the ones above it are likely tidal UDGs, while the ones below the MZR are consistent with early-quenching (failed galaxy scenarios), similar to what was suggested by B22.

In the bottom row of Fig. 8 we show the MATLAS UDGs colour-coded by their number of GCs. Differently from B22 and FM23, we do not impose the hard separation that GC-rich UDGs are the ones with $N_{GC} \ge 20$ and GC-poor UDGs the ones with fewer. Instead, the colour-scheme has been discretized so that evolutionary trends can be identified. In this figure, we see that GC-poor UDGs are overall more metal-rich galaxies, consistent with the classical dwarf MZR, and consistent with the recent findings of Jones et al. (2023) that gas-rich, field UDGs host few GCs. Alternatively, we see that all UDGs that follow the evolving MZR at high redshift are GC-rich. We do not find the population of GC-rich metal-rich cluster UDGs found by FM23, which may be expected as all of our UDGs are in groups/field, and thus less prone to tidal effects (see Fig. 7).

These results are in agreement with the recent findings of Pfeffer et al. (submitted), who used the EAGLE hydrodynamical simulations to find a correlation between the axis ratio and GC-richness of dwarf galaxies. In their study, round (higher $b/a$) dwarfs host more GCs than their elongated (low $b/a$) counterparts. As a a caveat, the works of Smith et al. (2013) and Smith et al. (2015) have shown that it is not possible to make a galaxy rounder through tidal stripping while also retaining a large number of GCs. If the harassment and stripping are strong enough to significantly alter their shapes, then the bulk of the GC system would be lost as well. In this study, however, we refer to galaxies that are round and GC-rich from their early stages.
### 5.1.2 Dependence with age and star formation time-scale

We now turn to understanding how the positioning of the UDGs in the stellar mass–metallicity plane depends on their SFH properties, e.g. their ages and star formation time-scales.

Ferré-Mateu et al. (2018) have shown that there exists an age dependence with projected clustercentric distance, i.e. UDGs are younger at larger projected clustercentric radii. Alabi et al. (2018) and Kadowaki et al. (2021) also found a colour dependence with the environment, with bluer UDGs residing in lower density environments than the redder ones. Finally, B22 found, using SED fitting, that UDGs in the field are on average younger and more metal-rich than cluster UDGs. Group UDGs in their work were found to be in a transitional space in the age–metallicity plane. However, these were based on very small number statistics, hence the need for this study focused mainly on group UDGs. Similar to these recent findings, the ages of the MATLAS UDGs also seem to distill into two main types: younger UDGs following the dwarf MZR and older ones following the high-$z$ MZR. We however caution the reader that ages recovered with PROSPECTOR are very sensitive to the version of the code used (see Section 3.2), as well as to the uncertainties in the optical photometry (see Appendix C).

One interesting recent finding on the stellar populations of UDGs by FM23 is a correlation between the $\alpha$ enrichment ([α/Fe]) and the SFHs of UDGs. Using a phase-space diagram, which may be particularly related to the environment of the UDGs, they show that the UDGs with the highest [α/Fe] are in the region corresponding to early infallers, while UDGs with the lowest [α/Fe] tend to be in the outskirts of clusters or in less dense environments. The $\alpha$-enrichment can be associated with the speed at which star formation starts and quenches within galaxies. Thus, this parameter can very well separate early-quenched galaxies (high [α/Fe]) from galaxies with prolonged SFHs (low [α/Fe]).

With SED fitting we cannot recover [α/Fe]. We can, however, look into the star formation time-scale of galaxies, i.e. the rate of the exponential decay of star formation, to understand how fast SF was quenched within the galaxy. In the bottom row of Fig. 7, we colour-code our UDGs by their mass-weighted ages ($t_M$) and star formation time-scales ($\tau$). We can see that UDGs that have the fastest star formation episodes (i.e. smallest $\tau$) are in the lowest regions of the plane, coinciding with the positioning of the older galaxies in the age panel. They are consistent with the evolving MZR at high redshift, indicative of early-quenching, which is in agreement with the results from FM23, and suggestive of a failed galaxy-like origin. On the other hand, UDGs that show prolonged SFHs (i.e. large $\tau$) are the ones that are younger and more metal-rich, aligned with a puffed-up dwarf formation scenario, as also found by FM23.

### 5.1.3 Dependence with GC-richness and specific frequency

One of the properties of UDGs that draws a lot of attention is the large number of GCs that some of them host in spite of their low stellar masses. Some extreme examples include DF44 (van Dokkum et al. 2017, but see also Saifollahi et al. 2022 and Forbes & Gannon 2024) and NGC 5846. UDG1/MATLAS-2019 (Müller et al. 2020; Forbes et al. 2021; Danieli et al. 2022). Simulations and analytical models have tried to explain the existence of these many GCs in such low-mass galaxies (e.g. Carleton et al. 2021; Trujillo-Gomez et al. 2021; Danieli et al. 2022; van Dokkum et al. 2022; Doppel et al. 2024). In parallel, observational work has focused on understanding what are the main differences between GC-poor and GC-rich UDGs,

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Difference in metallicity ($\delta_{MZR}$) between known MZRs and the MATLAS UDGs colour-coded by their number of GCs. Description is the same as Fig. 6. GC-poor MATLAS UDGs follow the classical dwarf MZR from Simon (2019). On another hand, all UDGs that follow the Ma et al. (2016) MZR at high redshift are GC-rich.</figcaption></figure>

in an attempt to correlate the higher fraction of GCs with other galaxy properties (Forbes et al. 2020; La Marca et al. 2022).

In terms of stellar populations, B22 used SED fitting techniques to find evidence for an inherent difference between GC-rich and GC-poor UDGs. They found that GC-poor UDGs are consistent with the classical dwarf MZR from Simon (2019), having consistently higher metallicities than their GC-rich counterparts. Interestingly, all of the UDGs that are consistent with the evolving MZR at $z = 2.2$ proposed by Ma et al. (2016) are GC-rich, suggesting a formation scenario relying on early-quenching for these galaxies. Following up on this finding, FM23 studied the stellar populations of 25 UDGs, the largest spectroscopic sample to date, using data from Keck/KCWI. They found that GC-poor UDGs consistently follow the classical dwarf MZR. However, they found that the GC-rich UDGs can be found both above and below the Simon (2019) MZR. They suggested that the ones above it are likely tidal UDGs, while the ones below the MZR are consistent with early-quenching (failed galaxy scenarios), similar to what was suggested by B22.

In the bottom row of Fig. 8 we show the MATLAS UDGs colour-coded by their number of GCs. Differently from B22 and FM23, we do not impose the hard separation that GC-rich UDGs are the ones with $N_{GC} \ge 20$ and GC-poor UDGs the ones with fewer. Instead, the colour-scheme has been discretized so that evolutionary trends can be identified. In this figure, we see that GC-poor UDGs are overall more metal-rich galaxies, consistent with the classical dwarf MZR, and consistent with the recent findings of Jones et al. (2023) that gas-rich, field UDGs host few GCs. Alternatively, we see that all UDGs that follow the evolving MZR at high redshift are GC-rich. We do not find the population of GC-rich metal-rich cluster UDGs found by FM23, which may be expected as all of our UDGs are in groups/field, and thus less prone to tidal effects (see Fig. 7).

These results are in agreement with the recent findings of Pfeffer et al. (submitted), who used the EAGLE hydrodynamical simulations to find a correlation between the axis ratio and GC-richness of dwarf galaxies. In their study, round (higher $b/a$) dwarfs host more GCs than their elongated (low $b/a$) counterparts. As a a caveat, the works of Smith et al. (2013) and Smith et al. (2015) have shown that it is not possible to make a galaxy rounder through tidal stripping while also retaining a large number of GCs. If the harassment and stripping are strong enough to significantly alter their shapes, then the bulk of the GC system would be lost as well. In this study, however, we refer to galaxies that are round and GC-rich from their early stages.
As a sidenote, the effective radius of galaxies, not only UDGs, has been shown to correlate with the number of GCs by Harris, Harris & Alessi (2013), also in agreement with our findings.

To understand if this trend between the metallicity and GC number is caused solely by the GC-richness or if it is being influenced by the mass of the galaxy, we also calculated the specific frequency per unit stellar mass of the galaxies, $T_N$ (Zepf & Whitmore 1993). This specific frequency is defined as:

$$T_N = \frac{N_{GC}}{\left(\frac{M_*}{10^9 M_\odot}\right)}. \qquad (2)$$

The MATLAS UDGs are colour-coded by their specific frequencies in the rightmost panel in the bottom row of Fig. 7. We can see now that the trend between GC-richness and metallicity is not so clear anymore. Although on average metal-poor UDGs have higher specific frequencies, there is a non-negligible population of more metal-rich UDGs that also have high $T_N$ values. This indicates that the trend found between the metallicity and the GC number may be incorporating correlations with the stellar mass. It is important to note, none the less, that the UDGs that follow the high-redshift MZR all show high values of $T_N$. This suggests that at least part of the reason that we do not find GC-rich, metal-rich UDGs is our lower stellar mass (on average) than FM23. When controlling for the stellar mass of the galaxies using $T_N$, this population becomes apparent. These GC-rich metal-rich galaxies are likely tidal UDGs, as suggested by FM23. This finding is reinforced by the fact that one confirmed tidal UDG in our sample, MATLAS-1824, has a similar metallicity and $T_N$ as these other GC-rich metal-rich UDGs.

Further supporting these findings, Lim et al. (2018) have studied the specific frequencies ($S_N$) of UDGs in the Coma cluster and also found two different populations of UDGs in terms of their GC content. In their work, UDGs with high $S_N$ were found to be more dark matter dominated, with mass-to-light (M/L) ratios > 1000, while UDGs with low $S_N$ have M/L ~500, consistent with the predictions of failed galaxy and puffed-up dwarf scenarios, respectively.

Interestingly, a transitional population is also identified that has on average ~10 GCs ($10 < T_N < 200$); and these do not follow either the classical dwarf MZR nor the high-z MZR (see light blue points in Fig. 8). Seeing a continuous trend in the GC numbers immediately leads to questioning whether GC-poor UDGs could evolve into GC-rich ones (or vice versa). However, these galaxies, if we put together all of the properties studied so far, seem to be inherently different, not only in terms of their GC systems. GC-poor UDGs are primarily found in the field or low-density groups, they are young, are elongated, and show prolonged SFHs. GC-rich UDGs, on the other hand, are found in clusters or high-density groups, are primarily round, have old to ancient ages, and have quenched long ago, having formed the bulk of their stars in early and quick star formation episodes. Given their different properties, an evolutionary path between the two populations seems unlikely and, thus, two (or more) formation scenarios seem to be necessary to explain these different types of UDGs (Prole et al. 2019; Jones et al. 2023).

## 5.2 Multiple classes of UDGs?

The results in the previous sections seemed to indicate that there exists at least two populations of UDGs. Driven by the idea that UDGs with different properties were formed by different processes, in what follows we investigate and compare the distribution and possible correlation between their global properties. Our goal is to identify

classes of objects with similar properties and possibly associate them with one of the prominent formation scenarios for UDGs.

Assuming the existence of distinct populations of UDGs (Forbes et al. 2020), and in order to better associate each galaxy to one of them, we applied the centroid-based clustering algorithm KMeans (MacQueen et al. 1967) to the data. KMeans has been applied to Astronomy in various contexts and was proven to be a valuable methodology to classify large data sets (see e.g. Sánchez Almeida et al. 2010; García-Dias et al. 2018; Coccato et al. 2022). The KMeans algorithm is an unsupervised learning method and is part of the Scikit-learn library (Pedregosa et al. 2011), responsible for implementing machine learning algorithms in Python.

In our implementation, KMeans groups together galaxies that have similar properties in the multiparameter space defined in our data set and flags them accordingly into classes. The initial parameters considered and compared here are stellar mass (log ($M_*/M_\odot$)), star formation time-scale ($\tau$), mass-weighted age ($t_M$), local volume density (log $\rho_{10}$), number of GCs ($N_{GC}$), GC specific frequency ($T_N$), central surface brightness ($\mu_{g,0}$), and the residual metallicity between the UDGs and the classical dwarf MZR ($\delta_{dwarf MZR}$). The latter was included instead of metallicity in order to remove the degeneracy between metallicity and stellar mass so both properties could be analysed separately. Positive values of $\delta_{dwarf MZR}$ mean that UDGs are above the dwarf MZR (i.e. more metal-rich), 0 is when they follow exactly the MZR, and negative values are for the cases where UDGs lie below the MZR (i.e. more metal-poor). The Sérsic index ($n$) and effective radius ($R_e$), were not included in the analysis as they are degenerate with the surface brightness. Amongst our total sample of 59 galaxies, only 38 UDGs have all of these physical quantities available (i.e. 21 UDGs do not have GC counts). Therefore, these 38 are the ones used in our analysis. This is a small sample of UDGs and we caution that results obtained from it need to be further tested in a larger sample of galaxies, as these may not be representative of a diverse populations of UDGs and should not be extrapolated to such. All properties were linearly scaled down to a range between 0 and 1 so that good measures of distance between classes could be assessed for every parameter analysed. The KMeans clustering algorithm is isotropically applied to all parameters in the multiparameter space and, therefore, tends to produce round (rather than elongated) clusters. In this case, if unscaled data is provided (i.e. with unequal variances), then more weight is given to variables with smaller variance. Thus, rescaling the data is a necessary step to allow all properties to be equally examined by the algorithm.

We applied KMeans to the data set allowing it to freely find the number of classes that best represents our data. We used the silhouette score quantity to measure the goodness of the clustering technique, where a high score means that the classes are well apart from each other and can be clearly distinguished. The score ranges from -1 to 1, where -1 means that sources are being associated with the wrong class, 0 means that classes are indistinguishable, and 1 means that the classes are perfectly separated. Two classes were found as the best number to describe our sample of UDGs, delivering the highest silhouette score (0.4) when compared to other numbers of classes. KMeans allocates every single data point to the closest class centre in the multiparameter space, which means that all UDGs analysed in this study are allocated to one of the classes. This does not mean that every UDG belonging to a particular class will have all of the properties compatible with that class, but that overall, most of their properties would be consistent. These classes were further assessed using principal component analysis, which showed that two clusters are able to recover 60 per cent of the variance of the data set. DBScan, another machine learning clustering algorithm, was also applied to
As a sidenote, the effective radius of galaxies, not only UDGs, has been shown to correlate with the number of GCs by Harris, Harris & Alessi (2013), also in agreement with our findings.

To understand if this trend between the metallicity and GC number is caused solely by the GC-richness or if it is being influenced by the mass of the galaxy, we also calculated the specific frequency per unit stellar mass of the galaxies, $T_N$ (Zepf & Whitmore 1993). This specific frequency is defined as:

$$T_N = \frac{N_{GC}}{\left(\frac{M_*}{10^9 M_{\odot}}\right)}. \qquad (2)$$

The MATLAS UDGs are colour-coded by their specific frequencies in the rightmost panel in the bottom row of Fig. 7. We can see now that the trend between GC-richness and metallicity is not so clear anymore. Although on average metal-poor UDGs have higher specific frequencies, there is a non-negligible population of more metal-rich UDGs that also have high $T_N$ values. This indicates that the trend found between the metallicity and the GC number may be incorporating correlations with the stellar mass. It is important to note, none the less, that the UDGs that follow the high-redshift MZR all show high values of $T_N$. This suggests that at least part of the reason that we do not find GC-rich, metal-rich UDGs is our lower stellar mass (on average) than FM23. When controlling for the stellar mass of the galaxies using $T_N$, this population becomes apparent. These GC-rich metal-rich galaxies are likely tidal UDGs, as suggested by FM23. This finding is reinforced by the fact that one confirmed tidal UDG in our sample, MATLAS-1824, has a similar metallicity and $T_N$ as these other GC-rich metal-rich UDGs.

Further supporting these findings, Lim et al. (2018) have studied the specific frequencies ($S_N$) of UDGs in the Coma cluster and also found two different populations of UDGs in terms of their GC content. In their work, UDGs with high $S_N$ were found to be more dark matter dominated, with mass-to-light (M/L) ratios > 1000, while UDGs with low $S_N$ have M/L ~500, consistent with the predictions of failed galaxy and puffed-up dwarf scenarios, respectively.

Interestingly, a transitional population is also identified that has on average ~10 GCs ($10 < T_N < 200$); and these do not follow either the classical dwarf MZR nor the high-z MZR (see light blue points in Fig. 8). Seeing a continuous trend in the GC numbers immediately leads to questioning whether GC-poor UDGs could evolve into GC-rich ones (or vice versa). However, these galaxies, if we put together all of the properties studied so far, seem to be inherently different, not only in terms of their GC systems. GC-poor UDGs are primarily found in the field or low-density groups, they are young, are elongated, and show prolonged SFHs. GC-rich UDGs, on the other hand, are found in clusters or high-density groups, are primarily round, have old to ancient ages, and have quenched long ago, having formed the bulk of their stars in early and quick star formation episodes. Given their different properties, an evolutionary path between the two populations seems unlikely and, thus, two (or more) formation scenarios seem to be necessary to explain these different types of UDGs (Prole et al. 2019; Jones et al. 2023).

## 5.2 Multiple classes of UDGs?

The results in the previous sections seemed to indicate that there exists at least two populations of UDGs. Driven by the idea that UDGs with different properties were formed by different processes, in what follows we investigate and compare the distribution and possible correlation between their global properties. Our goal is to identify

classes of objects with similar properties and possibly associate them with one of the prominent formation scenarios for UDGs.

Assuming the existence of distinct populations of UDGs (Forbes et al. 2020), and in order to better associate each galaxy to one of them, we applied the centroid-based clustering algorithm KMeans (MacQueen et al. 1967) to the data. KMeans has been applied to Astronomy in various contexts and was proven to be a valuable methodology to classify large data sets (see e.g. Sánchez Almeida et al. 2010; García-Dias et al. 2018; Coccato et al. 2022). The KMeans algorithm is an unsupervised learning method and is part of the Scikit-learn library (Pedregosa et al. 2011), responsible for implementing machine learning algorithms in Python.

In our implementation, KMeans groups together galaxies that have similar properties in the multiparameter space defined in our data set and flags them accordingly into classes. The initial parameters considered and compared here are stellar mass (log ($M_*/M_{\odot}$)), star formation time-scale ($\tau$), mass-weighted age ($t_M$), local volume density (log $\rho_{10}$), number of GCs ($N_{GC}$), GC specific frequency ($T_N$), central surface brightness ($\mu_{g,0}$), and the residual metallicity between the UDGs and the classical dwarf MZR ($\delta_{dwarf MZR}$). The latter was included instead of metallicity in order to remove the degeneracy between metallicity and stellar mass so both properties could be analysed separately. Positive values of $\delta_{dwarf MZR}$ mean that UDGs are above the dwarf MZR (i.e. more metal-rich), 0 is when they follow exactly the MZR, and negative values are for the cases where UDGs lie below the MZR (i.e. more metal-poor). The Sérsic index ($n$) and effective radius ($R_e$), were not included in the analysis as they are degenerate with the surface brightness. Amongst our total sample of 59 galaxies, only 38 UDGs have all of these physical quantities available (i.e. 21 UDGs do not have GC counts). Therefore, these 38 are the ones used in our analysis. This is a small sample of UDGs and we caution that results obtained from it need to be further tested in a larger sample of galaxies, as these may not be representative of a diverse populations of UDGs and should not be extrapolated to such. All properties were linearly scaled down to a range between 0 and 1 so that good measures of distance between classes could be assessed for every parameter analysed. The KMeans clustering algorithm is isotropically applied to all parameters in the multiparameter space and, therefore, tends to produce round (rather than elongated) clusters. In this case, if unscaled data is provided (i.e. with unequal variances), then more weight is given to variables with smaller variance. Thus, rescaling the data is a necessary step to allow all properties to be equally examined by the algorithm.

We applied KMeans to the data set allowing it to freely find the number of classes that best represents our data. We used the silhouette score quantity to measure the goodness of the clustering technique, where a high score means that the classes are well apart from each other and can be clearly distinguished. The score ranges from -1 to 1, where -1 means that sources are being associated with the wrong class, 0 means that classes are indistinguishable, and 1 means that the classes are perfectly separated. Two classes were found as the best number to describe our sample of UDGs, delivering the highest silhouette score (0.4) when compared to other numbers of classes. KMeans allocates every single data point to the closest class centre in the multiparameter space, which means that all UDGs analysed in this study are allocated to one of the classes. This does not mean that every UDG belonging to a particular class will have all of the properties compatible with that class, but that overall, most of their properties would be consistent. These classes were further assessed using principal component analysis, which showed that two clusters are able to recover 60 per cent of the variance of the data set. DBScan, another machine learning clustering algorithm, was also applied to
<figure><img src="image_11.png" /><figcaption>**Figure 9.** The two groups of UDGs identified using the KMeans clustering algorithm. In this figure, the radial axis shows the central value of a given property in each class, while the angular axis shows each property being considered by the clustering algorithm. This polar plot is to be read as the relative difference in the properties of UDGs belonging to the Classes A and B. Class A is shown in green and is comprised of UDGs that have lower stellar masses, prolonged SFHs (i.e. higher star formation time-scales), are more elongated (lower *b/a*), live in less dense environments (lower log *ρ*₁₀), host fewer GCs, are younger (smaller *t*ₘ), have higher δdwarf MZR (i.e. closer to zero and consistent with the classical dwarf MZR), and are fainter (higher μg, 0). As these properties agree with the predictions of puffed-up dwarf-like formation scenarios, we associate Class A with UDGs having this origin. Class B is shown in pink and includes UDGs with higher stellar masses, rapid SFHs (i.e. smaller τ), rounder (high *b/a*), living in the densest of our probed environments (higher log *ρ*₁₀), hosting on average the biggest GC systems, older (higher *t*ₘ), with smaller δdwarf MZR (i.e. negative values, being thus more metal-poor than the classical dwarf MZR), and brighter (smaller μg, 0). As these properties are aligned with the predictions of failed galaxy-like formation scenarios, we associate Class B with UDGs having formed through this pathway.</figcaption></figure>

the data and similarly returned that the best number of classes to represent this data set is two, with a silhouette score of 0.3.

The two classes of UDGs are shown in Fig. 9. This plot shows the relative difference between the properties of the UDGs belonging to each class and provides an intuitive way of understanding which properties provide the biggest separation between classes, and properties that contribute less. From Fig. 9, we can see that the best discriminators of UDGs belonging to different classes are log(M*/M☉), τ, tₘ, δdwarf MZR, and μg, 0. Properties that provide weaker separations are b/a, log ρ₁₀, and Nₐ₉. Several iterations of this plot were analysed to remove properties that did not contribute to the separations, or that were degenerate with others properties already included and were, therefore, not adding anything new to the separation. Properties excluded were Rₑ, n and PA for showing no significant difference between groups. At this point, it is important to remind the reader that in machine learning algorithms, cleaning the data set from any properties that are not contributing to the learning process is one of the most important tasks to recover reliable results. One might think that leaving parameters that return a null separation will not make a difference, but these add noise to the clustering process and could pollute and compromise the final recovered

<table><caption>**Table 1.** Mean values of UDG classes obtained with the KMeans clustering results. The first block is the classification done using Nₐ₉ (only 38 galaxies). The second block is the classification without Nₐ₉ (all 59 MATLAS UDGs).</caption><thead><tr><th>Parameter</th><th>Class A<br>(Puffed-up dwarf)</th><th>Class B<br>(Failed galaxy)</th></tr></thead><tbody><tr><td>Number</td><td>23</td><td>15</td></tr><tr><td>log (M<sub>*</sub>/M<sub>☉</sub>)</td><td>7.31 ± 0.17</td><td>7.87 ± 0.16</td></tr><tr><td>τ (Gyr)</td><td>2.27 ± 0.75</td><td>1.58 ± 0.96</td></tr><tr><td>b/a</td><td>0.59 ± 0.19</td><td>0.69 ± 0.17</td></tr><tr><td>log ρ<sub>10</sub> (Mpc<sup>-3</sup>)</td><td>-1.53 ± 0.56</td><td>-1.36 ± 0.59</td></tr><tr><td>N<sub>GC</sub></td><td>2.27 ± 2.30</td><td>13.09 ± 11.01</td></tr><tr><td>t<sub>M</sub> (Gyr)</td><td>5.74 ± 2.02</td><td>7.34 ± 2.36</td></tr><tr><td>δ<sub>dwarf MZR</sub> (dex)</td><td>-0.05 ± 0.23</td><td>-0.33 ± 0.23</td></tr><tr><td>μ<sub>g, 0</sub> (mag arcsec<sup>-2</sup>)</td><td>25.91 ± 0.94</td><td>25.19 ± 0.62</td></tr><tr><td>Number</td><td>33</td><td>26</td></tr><tr><td>log (M<sub>*</sub>/M<sub>☉</sub>)</td><td>7.28 ± 0.19</td><td>7.71 ± 0.27</td></tr><tr><td>τ (Gyr)</td><td>2.68 ± 0.59</td><td>1.54 ± 0.81</td></tr><tr><td>b/a</td><td>0.56 ± 0.21</td><td>0.67 ± 0.16</td></tr><tr><td>log ρ<sub>10</sub> (Mpc<sup>-3</sup>)</td><td>-1.57 ± 0.54</td><td>-1.41 ± 0.53</td></tr><tr><td>t<sub>M</sub> (Gyr)</td><td>5.07 ± 1.79</td><td>7.25 ± 2.20</td></tr><tr><td>δ<sub>dwarf MZR</sub> (dex)</td><td>-0.02 ± 0.23</td><td>-0.31 ± 0.21</td></tr><tr><td>μ<sub>g, 0</sub> (mag arcsec<sup>-2</sup>)</td><td>26.04 ± 1.06</td><td>25.33 ± 0.64</td></tr></tbody></table>

results. We emphasize, none the less, that the recovered classes are consistent independently of the addition of these parameters, however, the classes are found to have much larger uncertainties as the null parameters contaminate the clear separation provided by the significant ones.

On the other hand, it is important to analyse carefully the properties that are providing the biggest separation to ensure these are physically meaningful and not just a result of selection effects or the properties being degenerate/correlated. For example, the surface brightness, effective radius, and number of GCs are known to correlate with the stellar mass. We thus carried out tests excluding every single one of these properties from the clustering algorithm to check whether the classes would change or not. We found that the recovered classes are similar regardless if any of these properties is excluded or included in the process, they however result in larger uncertainties on the class centres. Overall, in all iterations of applying KMeans to our data set (i.e. with the inclusion/omission of properties) the classes maintained, except for one iteration. The only property capable of changing the identified classes was Tₙ, this is because even though all of the GC-rich galaxies have on average high specific frequencies, GC-poor UDGs show a mix of specific frequencies, capable of confusing the clustering algorithm and strongly changing the identified classes. Whether the separation with Tₙ is more correct than the one without it needs to be further assessed in a larger sample of galaxies. In this study, we choose to use the results that are stable in the absence of Tₙ. As a last comment, we note that if we exclude the GC number from the analysis, we can apply KMeans to the full sample of 59 UDGs (as the GC number is the only property that is not available for all galaxies). Once this is done, the exact same classes are identified, and the variance on the properties is smaller as the classes are comprised of a larger sample of galaxies. This gives us confidence that our clustering methodology is capturing the main properties of our data set even when constrained to a smaller subsample.

The properties of the two classes of UDGs identified with KMeans as summarized in Table 1 are as follows:

(i) Class A. UDGs in this class have lower stellar masses, prolonged SFHs (i.e. higher star formation time-scales), are more
<figure><img src="image_10.png" /><figcaption>**Figure 9.** The two groups of UDGs identified using the KMeans clustering algorithm. In this figure, the radial axis shows the central value of a given property in each class, while the angular axis shows each property being considered by the clustering algorithm. This polar plot is to be read as the relative difference in the properties of UDGs belonging to the Classes A and B. Class A is shown in green and is comprised of UDGs that have lower stellar masses, prolonged SFHs (i.e. higher star formation time-scales), are more elongated (lower *b/a*), live in less dense environments (lower log *ρ*₁₀), host fewer GCs, are younger (smaller *t*ₘ), have higher δdwarf MZR (i.e. closer to zero and consistent with the classical dwarf MZR), and are fainter (higher μg, 0). As these properties agree with the predictions of puffed-up dwarf-like formation scenarios, we associate Class A with UDGs having this origin. Class B is shown in pink and includes UDGs with higher stellar masses, rapid SFHs (i.e. smaller τ), rounder (high *b/a*), living in the densest of our probed environments (higher log *ρ*₁₀), hosting on average the biggest GC systems, older (higher *t*ₘ), with smaller δdwarf MZR (i.e. negative values, being thus more metal-poor than the classical dwarf MZR), and brighter (smaller μg, 0). As these properties are aligned with the predictions of failed galaxy-like formation scenarios, we associate Class B with UDGs having formed through this pathway.</figcaption></figure>

the data and similarly returned that the best number of classes to represent this data set is two, with a silhouette score of 0.3.

The two classes of UDGs are shown in Fig. 9. This plot shows the relative difference between the properties of the UDGs belonging to each class and provides an intuitive way of understanding which properties provide the biggest separation between classes, and properties that contribute less. From Fig. 9, we can see that the best discriminators of UDGs belonging to different classes are log(M*/M☉), τ, tₘ, δdwarf MZR, and μg, 0. Properties that provide weaker separations are b/a, log ρ₁₀, and Nₐ₉. Several iterations of this plot were analysed to remove properties that did not contribute to the separations, or that were degenerate with others properties already included and were, therefore, not adding anything new to the separation. Properties excluded were Rₑ, n and PA for showing no significant difference between groups. At this point, it is important to remind the reader that in machine learning algorithms, cleaning the data set from any properties that are not contributing to the learning process is one of the most important tasks to recover reliable results. One might think that leaving parameters that return a null separation will not make a difference, but these add noise to the clustering process and could pollute and compromise the final recovered

<table><caption>**Table 1.** Mean values of UDG classes obtained with the KMeans clustering results. The first block is the classification done using Nₐ₉ (only 38 galaxies). The second block is the classification without Nₐ₉ (all 59 MATLAS UDGs).</caption><thead><tr><th>Parameter</th><th>Class A<br>(Puffed-up dwarf)</th><th>Class B<br>(Failed galaxy)</th></tr></thead><tbody><tr><td>Number</td><td>23</td><td>15</td></tr><tr><td>log (M<sub>*</sub>/M<sub>☉</sub>)</td><td>7.31 ± 0.17</td><td>7.87 ± 0.16</td></tr><tr><td>τ (Gyr)</td><td>2.27 ± 0.75</td><td>1.58 ± 0.96</td></tr><tr><td>b/a</td><td>0.59 ± 0.19</td><td>0.69 ± 0.17</td></tr><tr><td>log ρ<sub>10</sub> (Mpc<sup>-3</sup>)</td><td>-1.53 ± 0.56</td><td>-1.36 ± 0.59</td></tr><tr><td>N<sub>GC</sub></td><td>2.27 ± 2.30</td><td>13.09 ± 11.01</td></tr><tr><td>t<sub>M</sub> (Gyr)</td><td>5.74 ± 2.02</td><td>7.34 ± 2.36</td></tr><tr><td>δ<sub>dwarf MZR</sub> (dex)</td><td>-0.05 ± 0.23</td><td>-0.33 ± 0.23</td></tr><tr><td>μ<sub>g, 0</sub> (mag arcsec<sup>-2</sup>)</td><td>25.91 ± 0.94</td><td>25.19 ± 0.62</td></tr><tr><td>Number</td><td>33</td><td>26</td></tr><tr><td>log (M<sub>*</sub>/M<sub>☉</sub>)</td><td>7.28 ± 0.19</td><td>7.71 ± 0.27</td></tr><tr><td>τ (Gyr)</td><td>2.68 ± 0.59</td><td>1.54 ± 0.81</td></tr><tr><td>b/a</td><td>0.56 ± 0.21</td><td>0.67 ± 0.16</td></tr><tr><td>log ρ<sub>10</sub> (Mpc<sup>-3</sup>)</td><td>-1.57 ± 0.54</td><td>-1.41 ± 0.53</td></tr><tr><td>t<sub>M</sub> (Gyr)</td><td>5.07 ± 1.79</td><td>7.25 ± 2.20</td></tr><tr><td>δ<sub>dwarf MZR</sub> (dex)</td><td>-0.02 ± 0.23</td><td>-0.31 ± 0.21</td></tr><tr><td>μ<sub>g, 0</sub> (mag arcsec<sup>-2</sup>)</td><td>26.04 ± 1.06</td><td>25.33 ± 0.64</td></tr></tbody></table>

results. We emphasize, none the less, that the recovered classes are consistent independently of the addition of these parameters, however, the classes are found to have much larger uncertainties as the null parameters contaminate the clear separation provided by the significant ones.

On the other hand, it is important to analyse carefully the properties that are providing the biggest separation to ensure these are physically meaningful and not just a result of selection effects or the properties being degenerate/correlated. For example, the surface brightness, effective radius, and number of GCs are known to correlate with the stellar mass. We thus carried out tests excluding every single one of these properties from the clustering algorithm to check whether the classes would change or not. We found that the recovered classes are similar regardless if any of these properties is excluded or included in the process, they however result in larger uncertainties on the class centres. Overall, in all iterations of applying KMeans to our data set (i.e. with the inclusion/omission of properties) the classes maintained, except for one iteration. The only property capable of changing the identified classes was Tₙ, this is because even though all of the GC-rich galaxies have on average high specific frequencies, GC-poor UDGs show a mix of specific frequencies, capable of confusing the clustering algorithm and strongly changing the identified classes. Whether the separation with Tₙ is more correct than the one without it needs to be further assessed in a larger sample of galaxies. In this study, we choose to use the results that are stable in the absence of Tₙ. As a last comment, we note that if we exclude the GC number from the analysis, we can apply KMeans to the full sample of 59 UDGs (as the GC number is the only property that is not available for all galaxies). Once this is done, the exact same classes are identified, and the variance on the properties is smaller as the classes are comprised of a larger sample of galaxies. This gives us confidence that our clustering methodology is capturing the main properties of our data set even when constrained to a smaller subsample.

The properties of the two classes of UDGs identified with KMeans as summarized in Table 1 are as follows:

(i) Class A. UDGs in this class have lower stellar masses, prolonged SFHs (i.e. higher star formation time-scales), are more
elongated (lower *b/a*), live in less dense environments (lower log ρ₁₀), host fewer GCs, are younger (smaller *t*<sub>M</sub>), have higher δ<sub>dwarf MZR</sub> (i.e. closer to zero and consistent with the classical dwarf MZR), and are fainter (higher μ<sub>g, 0</sub>).

(ii) Class B. UDGs in this class have higher stellar masses, rapid SFHs (i.e. smaller τ), are rounder (high *b/a*), live in the densest of our probed environments (higher log ρ₁₀), some host the biggest GC systems, are older (higher *t*<sub>M</sub>), have smaller δ<sub>dwarf MZR</sub> (i.e. negative values, being thus more metal-poor than the classical dwarf MZR), and are brighter (smaller μ<sub>g, 0</sub>).

The centres of the two identified classes and the number of UDGs associated with each class are shown in the upper block of Table 1 for the classification including N<sub>GC</sub>. The lower block, on the other hand, shows the centres and number of UDGs following the classification applied to all 59 galaxies (i.e. removing N<sub>GC</sub> from the analysis). The quoted uncertainties are the standard deviation of each analysed property considering the galaxies that belong to that class, i.e. an approximation of how far apart objects belonging to the same class can be. The cluster centres and uncertainties can also be seen in the histogram distribution of the properties in Fig. 10. Properties that provide the best separation between classes will have more clear bimodal distributions, while the ones that are not strong discriminators will have similar Gaussian distributions. The classes each UDG in our sample are associated in the analyses with and without N<sub>GC</sub> are given in Table B2. It is reinforcing to see that 37 out of the 38 UDGs included in both analyses (i.e. with and without N<sub>GC</sub>) have preserved their classes. Only one UDG (MATLAS-984) changed classes when the number of GCs was excluded. We note that changing classes is expected if the classification is being driven by the property that was removed (i.e. N<sub>GC</sub>) and the other properties are in overlapping cluster regions. This is the case for MATLAS-984, which has properties that could fit into any class when N<sub>GC</sub> is not considered (see Tables B2 and B3).

We emphasize that these classes were obtained using a small sample of UDGs, all in low-density environments. Thus, the centres found may not be representative of a vast population of UDGs with strongly divergent properties. These cluster centres should perhaps be interpreted as a measurement of which properties are important discriminators of UDG formation scenarios and which ones are not. Although there is not a net separation or threshold in the properties of UDGs formed via the two scenarios, and possibly not a one-to-one correspondence between a given class and a formation scenario, we recognize that UDGs in Class A have properties (see Table 1) that are more aligned with being formed through puffed-up dwarf scenarios, i.e. properties similar to those of classical dwarf galaxies. On the other hand, UDGs in Class B share similar properties to those of failed galaxy scenarios (see Table 1), i.e. properties consistent with early-quenching. We note that tidal UDGs, due to low number statistics, was not found as a third class in our data set, although these types of UDGs have overall different properties than puffed-up dwarfs or failed galaxies. This class is expected to become clearer in larger samples of UDGs, including more massive and cluster ones, where tidal effects are more present.

This clustering algorithm was applied to a small set of galaxies and needs to be further tested, both with a larger number of UDGs and in comparison with classical dwarfs. This way, a more comprehensive comparison between subpopulations can be made. We note, none the less, that this separation (in terms of GC number, axis ratio, and mass) has also been reported in the Coma cluster (Lim et al. 2018), so these two classes of UDGs may be a general separation regardless of their environments. A study combining and comparing

the MATLAS UDGs, the UDGs studied in B22 and a control sample of classical MATLAS dwarfs is under preparation (Buzzo et al., in preparation) and will focus on deepening these clustering findings. Other properties can also be added to the clustering in the future to enrich the classification, namely [α/Fe] (e.g. Ferré-Mateu et al. 2018; Villaume et al. 2022, FM23), dynamical masses and velocity dispersions (e.g. Gannon et al. 2020, 2021, 2022, 2023a), positioning in the phase-space diagram (e.g. Gannon et al. 2022; Forbes et al. 2023, FM23), etc. As a final remark, once the classes are further tested and if they are proven to be well founded, they could be used to predict unknown properties of UDGs. For example, by knowing the stellar populations and morphology of a given UDG, using these clustering results, one would be able to predict if it is GC-poor or GC-rich, in cases where knowing the GC-richness of the galaxy would otherwise be impossible (e.g. at distance higher > 100 Mpc), further showing the importance of such findings.

# 6 CONCLUSIONS

In this study, we have used the fully Bayesian Monte Carlo Markov Chain inference code PROSPECTOR to perform SED fitting on the 59 MATLAS UDGs using data from the optical to the mid-IR. We use PROSPECTOR to recover stellar populations using two different configurations, one with the redshift as a free parameter, and one with the redshift fixed to the redshift of the group where the UDG was identified.

We find that the models with the redshift fixed to that of the host group are the most appropriate to study these galaxies. This assumption is based on the results of Heesters et al. (2023), who found that 75 per cent of the MATLAS dwarfs studied by them are at the same distance as their massive hosts. This is supported by our study with Keck/DEIMOS of three MATLAS UDGs, also found to be at the same redshift as their hosts.

Our results indicate that the MATLAS UDGs have intermediate-to-old ages, with an average mass-weighted age of 7.1 ± 1.8 Gyr. They are also metal-poor with an average [M/H] of −1.2 ± 0.2 dex. They have an average star formation time-scale τ of 1.6 ± 0.7 Gyr and are consistent with no dust attenuation, displaying an average of $A_V = 0.12 \pm 0.07$ mag.

When studying the MATLAS UDGs in the mass-metallicity plane, we find that their distribution is highly bimodal, with the first mode being well explained by the classical dwarf MZR, while the second one is more consistent with the evolving MZR at a high redshift.

To further investigate the cause of this bimodality, we study the positioning of the UDGs in the mass-metallicity plane according to their surface brightness, effective radii, axis ratios, local volume densities, mass-weighted ages, star formation time-scales, number of GCs and GC specific frequencies. We find that UDGs split into two main classes:

**Class A:** They have lower stellar masses, prolonged SFHs, are more elongated, live in less dense environments, host fewer GCs, are younger, are consistent with the classical dwarf MZR, and are fainter.

**Class B:** UDGs in this class have higher stellar masses, rapid SFHs, are rounder, live in the densest of our probed environments, host on average the biggest GC systems, are older, lie below the classical dwarf MZR, and are brighter. These galaxies are overall better explained by the evolving MZR at high redshift, i.e. consistent with early-quenching.

This overall picture seems to indicate that UDGs belonging to Class A, combining all of the aforementioned properties, are better explained by puffed-up dwarf formation scenarios. Alternatively,
elongated (lower *b/a*), live in less dense environments (lower log ρ₁₀), host fewer GCs, are younger (smaller *t*<sub>M</sub>), have higher δ<sub>dwarf MZR</sub> (i.e. closer to zero and consistent with the classical dwarf MZR), and are fainter (higher μ<sub>g, 0</sub>).

(ii) Class B. UDGs in this class have higher stellar masses, rapid SFHs (i.e. smaller τ), are rounder (high *b/a*), live in the densest of our probed environments (higher log ρ₁₀), some host the biggest GC systems, are older (higher *t*<sub>M</sub>), have smaller δ<sub>dwarf MZR</sub> (i.e. negative values, being thus more metal-poor than the classical dwarf MZR), and are brighter (smaller μ<sub>g, 0</sub>).

The centres of the two identified classes and the number of UDGs associated with each class are shown in the upper block of Table 1 for the classification including N<sub>GC</sub>. The lower block, on the other hand, shows the centres and number of UDGs following the classification applied to all 59 galaxies (i.e. removing N<sub>GC</sub> from the analysis). The quoted uncertainties are the standard deviation of each analysed property considering the galaxies that belong to that class, i.e. an approximation of how far apart objects belonging to the same class can be. The cluster centres and uncertainties can also be seen in the histogram distribution of the properties in Fig. 10. Properties that provide the best separation between classes will have more clear bimodal distributions, while the ones that are not strong discriminators will have similar Gaussian distributions. The classes each UDG in our sample are associated in the analyses with and without N<sub>GC</sub> are given in Table B2. It is reinforcing to see that 37 out of the 38 UDGs included in both analyses (i.e. with and without N<sub>GC</sub>) have preserved their classes. Only one UDG (MATLAS-984) changed classes when the number of GCs was excluded. We note that changing classes is expected if the classification is being driven by the property that was removed (i.e. N<sub>GC</sub>) and the other properties are in overlapping cluster regions. This is the case for MATLAS-984, which has properties that could fit into any class when N<sub>GC</sub> is not considered (see Tables B2 and B3).

We emphasize that these classes were obtained using a small sample of UDGs, all in low-density environments. Thus, the centres found may not be representative of a vast population of UDGs with strongly divergent properties. These cluster centres should perhaps be interpreted as a measurement of which properties are important discriminators of UDG formation scenarios and which ones are not. Although there is not a net separation or threshold in the properties of UDGs formed via the two scenarios, and possibly not a one-to-one correspondence between a given class and a formation scenario, we recognize that UDGs in Class A have properties (see Table 1) that are more aligned with being formed through puffed-up dwarf scenarios, i.e. properties similar to those of classical dwarf galaxies. On the other hand, UDGs in Class B share similar properties to those of failed galaxy scenarios (see Table 1), i.e. properties consistent with early-quenching. We note that tidal UDGs, due to low number statistics, was not found as a third class in our data set, although these types of UDGs have overall different properties than puffed-up dwarfs or failed galaxies. This class is expected to become clearer in larger samples of UDGs, including more massive and cluster ones, where tidal effects are more present.

This clustering algorithm was applied to a small set of galaxies and needs to be further tested, both with a larger number of UDGs and in comparison with classical dwarfs. This way, a more comprehensive comparison between subpopulations can be made. We note, none the less, that this separation (in terms of GC number, axis ratio, and mass) has also been reported in the Coma cluster (Lim et al. 2018), so these two classes of UDGs may be a general separation regardless of their environments. A study combining and comparing

the MATLAS UDGs, the UDGs studied in B22 and a control sample of classical MATLAS dwarfs is under preparation (Buzzo et al., in preparation) and will focus on deepening these clustering findings. Other properties can also be added to the clustering in the future to enrich the classification, namely [α/Fe] (e.g. Ferré-Mateu et al. 2018; Villaume et al. 2022, FM23), dynamical masses and velocity dispersions (e.g. Gannon et al. 2020, 2021, 2022, 2023a), positioning in the phase-space diagram (e.g. Gannon et al. 2022; Forbes et al. 2023, FM23), etc. As a final remark, once the classes are further tested and if they are proven to be well founded, they could be used to predict unknown properties of UDGs. For example, by knowing the stellar populations and morphology of a given UDG, using these clustering results, one would be able to predict if it is GC-poor or GC-rich, in cases where knowing the GC-richness of the galaxy would otherwise be impossible (e.g. at distance higher > 100 Mpc), further showing the importance of such findings.

# 6 CONCLUSIONS

In this study, we have used the fully Bayesian Monte Carlo Markov Chain inference code PROSPECTOR to perform SED fitting on the 59 MATLAS UDGs using data from the optical to the mid-IR. We use PROSPECTOR to recover stellar populations using two different configurations, one with the redshift as a free parameter, and one with the redshift fixed to the redshift of the group where the UDG was identified.

We find that the models with the redshift fixed to that of the host group are the most appropriate to study these galaxies. This assumption is based on the results of Heesters et al. (2023), who found that 75 per cent of the MATLAS dwarfs studied by them are at the same distance as their massive hosts. This is supported by our study with Keck/DEIMOS of three MATLAS UDGs, also found to be at the same redshift as their hosts.

Our results indicate that the MATLAS UDGs have intermediate-to-old ages, with an average mass-weighted age of 7.1 ± 1.8 Gyr. They are also metal-poor with an average [M/H] of -1.2 ± 0.2 dex. They have an average star formation time-scale τ of 1.6 ± 0.7 Gyr and are consistent with no dust attenuation, displaying an average of A<sub>V</sub> = 0.12 ± 0.07 mag.

When studying the MATLAS UDGs in the mass-metallicity plane, we find that their distribution is highly bimodal, with the first mode being well explained by the classical dwarf MZR, while the second one is more consistent with the evolving MZR at a high redshift.

To further investigate the cause of this bimodality, we study the positioning of the UDGs in the mass-metallicity plane according to their surface brightness, effective radii, axis ratios, local volume densities, mass-weighted ages, star formation time-scales, number of GCs and GC specific frequencies. We find that UDGs split into two main classes:

**Class A:** They have lower stellar masses, prolonged SFHs, are more elongated, live in less dense environments, host fewer GCs, are younger, are consistent with the classical dwarf MZR, and are fainter.

**Class B:** UDGs in this class have higher stellar masses, rapid SFHs, are rounder, live in the densest of our probed environments, host on average the biggest GC systems, are older, lie below the classical dwarf MZR, and are brighter. These galaxies are overall better explained by the evolving MZR at high redshift, i.e. consistent with early-quenching.

This overall picture seems to indicate that UDGs belonging to Class A, combining all of the aforementioned properties, are better explained by puffed-up dwarf formation scenarios. Alternatively,
<figure><img src="image_12.png" /><figcaption>**Figure 10.** The two groups of UDGs identified using the KMeans clustering algorithm. Distribution of UDGs in the parameter space defined by stellar mass (log (*M*<sub>⋆</sub>/*M*<sub>☉</sub>)), star formation time-scale (τ), axis ratio (*b*/*a*), local volume density (log ρ<sub>10</sub>), number of GCs (*N*<sub>GC</sub>), mass-weighted age (*t*<sub>M</sub>), δ<sub>dwarf MZR</sub> (i.e. level of consistency between the UDGs and the classical dwarf MZR), and central surface brightness (μ<sub>g, 0</sub>). The green points indicate galaxies that belong to class A (associated with the 'puffed-up dwarf' scenario), pink points indicate galaxies that belong to class B (associated with the 'failed galaxy' scenario). The histograms are smoothed using a Kernel Density Estimation (KDE) technique and indicate the distribution of a given parameter; green and pink histograms indicate those of class A and B, respectively.</figcaption></figure>

UDGs of Class B, taking into consideration all of the morphological
and stellar population properties found, seem to be more aligned with
failed galaxy formation scenarios, in agreement with the spectro-
scopic studies. Further testing of this clustering technique, including

a larger sample of UDGs and classical dwarfs for comparison, is
under preparation (Buzzo et al., in preparation).

This paper provides a continuation of the photometric study of the
stellar populations of UDGs across the sky. We demonstrate that SED
<figure><img src="image_11.png" /><figcaption>**Figure 10.** The two groups of UDGs identified using the KMeans clustering algorithm. Distribution of UDGs in the parameter space defined by stellar mass (log (*M*<sub>☉</sub>/*M*<sub>☉</sub>)), star formation time-scale (τ), axis ratio (*b*/*a*), local volume density (log ρ<sub>10</sub>), number of GCs (*N*<sub>GC</sub>), mass-weighted age (*t*<sub>M</sub>), δ<sub>dwarf MZR</sub> (i.e. level of consistency between the UDGs and the classical dwarf MZR), and central surface brightness (μ<sub>g, 0</sub>). The green points indicate galaxies that belong to class A (associated with the 'puffed-up dwarf' scenario), pink points indicate galaxies that belong to class B (associated with the 'failed galaxy' scenario). The histograms are smoothed using a Kernel Density Estimation (KDE) technique and indicate the distribution of a given parameter; green and pink histograms indicate those of class A and B, respectively.</figcaption></figure>

UDGs of Class B, taking into consideration all of the morphological
and stellar population properties found, seem to be more aligned with
failed galaxy formation scenarios, in agreement with the spectro-
scopic studies. Further testing of this clustering technique, including

a larger sample of UDGs and classical dwarfs for comparison, is
under preparation (Buzzo et al., in preparation).

This paper provides a continuation of the photometric study of the
stellar populations of UDGs across the sky. We demonstrate that SED
fitting techniques, coupled with a broad wavelength coverage, are an
important approach to statistically understand the origins of UDGs
and their impact on the field of galaxy formation and evolution.

Software: ASTROPY (Astropy Collaboration 2013, 2018), FSPS (Conroy & Gunn 2010; Conroy et al. 2010), PYTHON-FSPS (Johnson et al. 2021a), PROSPECTOR (Leja et al. 2017; Johnson et al. 2021b).

ACKNOWLEDGEMENTS

We thank the anonymous referee for suggestions that greatly im-
proved the manuscript.

We are thankful to Patrick Durrell for his contributions to the MATLAS survey. We thank Pavan Uttarkar for all the help with KMeans and other clustering algorithms tested and used in this work. We also thank Amelia Fraser-McKelvie for her colloquium at Swinburne University which gave us the idea of applying clustering algorithms to our data set. We thank Yimeng Tang for important discussions about PROSPECTOR'S configurations and outputs.

This research was supported by the Australian Research Coun-
cil Centre of Excellence for All Sky Astrophysics in 3 Di-
mensions (ASTRO 3D), through project number CE170100013.
DAF, JPB, JP, and WJC thank the ARC for support via
DP220101863. AJR was supported by National Science Foundation
grant AST-2308390. AFM has received support from CEX2019-
00920-S, RYC2021-031099-I, and PID2021-123313NA-100
of MICIN/AEI/10.13039/501100011033/FEDER,UE, NextGenera-
tionEU/PRT. SL acknowledges the support from the Sejong Science
Fellowship Program through the National Research Foundation of
Korea (NRF-2021R1C1C2006790).

This paper is based in part on observations from the Legacy Survey, which consists of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing–Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). This publication uses data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration.

DATA AVAILABILITY

DECaLS data are available via the Legacy survey portal. WISE data
are available via the WISE archive. MATLAS data are available from
the CDS archives and through the MATLAS portal.

REFERENCES

Alabi A. et al., 2018, MNRAS, 479, 3308

Amorisco N. C., Loeb A., 2016, MNRAS, 459, L51

Astropy Collaboration, 2013, A&A, 558, A33

Astropy Collaboration, 2018, AJ, 156, 123

Barbosa C. E. et al., 2020, ApJS, 247, 46

Beasley M. A., Romanowsky A. J., Pota V., Navarro I. M., Martinez Delgado
D., Neyer F., Deich A. L., 2016, ApJ, 819, L20

Benavides J. A. et al., 2021, Nat. Astron., 5, 1255

Bertin E., Mellier Y., Radovich M., Missonnier G., Didelon P., Morin B., 2002, in Bohlender D. A., Durand D., Handley T. H., eds, ASP Conf. Ser. Vol. 281, Astronomical Data Analysis Software and Systems XI. Astron. Soc. Pac., San Francisco, p. 228

Bílek M. et al., 2020, MNRAS, 498, 2138

Burkert A., Forbes D. A., 2020, AJ, 159, 56

Buzzo M. L. et al., 2021, MNRAS, 504, 2146

Buzzo M. L. et al., 2022, MNRAS, 517, 2231(B22)

Buzzo M. L., Forbes D. A., Brodie J. P., Janssens S. R., Couch W. J., Romanowsky A. J., Gannon J. S., 2023, MNRAS, 522, 595

Calzetti D., Armus L., Bohlin R. C., Kinney A. L., Koornneef J., Storchi-Bergmann T., 2000, ApJ, 533, 682

Cappellari M., 2017, MNRAS, 466, 798

Cardona-Barrero S., Di Cintio A., Brook C. B. A., Ruiz-Lara T., Beasley M. A., Falcón-Barroso J., Macciò A. V., 2020, MNRAS, 497, 4282

Carleton T., Errani R., Cooper M., Kaplinghat M., Peñarrubia J., Guo Y., 2019, MNRAS, 485, 382

Carleton T., Guo Y., Munshi F., Tremmel M., Wright A., 2021, MNRAS, 502, 398

Carnall A. C., McLure R. J., Dunlop J. S., Davé R., 2018, MNRAS, 480, 4379

Chan T. K., Kereš D., Wetzel A., Hopkins P. F., Faucher-Giguère C. A., El-Badry K., Garrison-Kimmel S., Boylan-Kolchin M., 2018, MNRAS, 478, 906

Chartab N., Newman A. B., Rudie G. C., Blanc G. A., Kelson D. D., 2024, ApJ, 960, 73

Chilingarian I. V., Afanasiev A. V., Grishin K. A., Fabricant D., Moran S., 2019, ApJ, 884, 79

Coccato L., Fraser-McKelvie A., Jaffé Y. L., Johnston E. J., Cortesi A., Pallero D., 2022, MNRAS, 515, 201

Collins M. L. M., Read J. I., 2022, Nat. Astron., 6, 647

Conroy C., Gunn J. E., 2010, ApJ, 712, 833

Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486

Conroy C., White M., Gunn J. E., 2010, ApJ, 708, 58

Danieli S. et al., 2022, ApJ, 927, L28

Davis M., Geller M. J., 1976, ApJ, 208, 13

Dey A. et al., 2019, AJ, 157, 168

Di Cintio A., Brook C. B., Dutton A. A., Macciò A. V., Obreja A., Dekel A., 2017, MNRAS, 466, L1

Doppel J. E., Sales L. V., Benavides J. A., Toloba E., Peng E. W., Nelson D., Navarro J. F., 2024, MNRAS, 529, 1827

Dressler A., 1980, ApJ, 236, 351

Dressler A. et al., 1997, ApJ, 490, 577

Duc P.-A., 2020, preprint (arXiv:2007.13874)

Duc P.-A., Paudel S., McDermid R. M., Cuillandre J.-C., Serra P., Bournaud F., Cappellari M., Emsellem E., 2014, MNRAS, 440, 1458

Duc P.-A. et al., 2015, MNRAS, 446, 120

Faber S. M. et al., 2003, in Iye M., Moorwood A. F. M., eds, Proc. SPIE Conf. Ser. Vol. 4841, Instrument Design and Performance for Optical/Infrared Ground-based Telescopes. SPIE, Bellingham, p. 1657

Ferré-Mateu A. et al., 2018, MNRAS, 479, 4891

Ferré-Mateu A., Gannon J. S., Forbes D. A., Buzzo M. L., Romanowsky A. J., Brodie J. P., 2023, MNRAS, 526, 4735(FM23)

Forbes D. A., Gannon J., 2024, MNRAS, 528, 608

Forbes D. A., Alabi A., Romanowsky A. J., Brodie J. P., Arimoto N., 2020, MNRAS, 492, 4874

Forbes D. A., Gannon J. S., Romanowsky A. J., Alabi A., Brodie J. P., Couch W. J., Ferré-Mateu A., 2021, MNRAS, 500, 1279

Forbes D. A. et al., 2023, MNRAS, 525, L93

Gannon J. S., Forbes D. A., Romanowsky A. J., Ferré-Mateu A., Couch W. J., Brodie J. P., 2020, MNRAS, 495, 2582

Gannon J. S. et al., 2021, MNRAS, 502, 3144

Gannon J. S. et al., 2022, MNRAS, 510, 946

Gannon J. S., Forbes D. A., Brodie J. P., Romanowsky A. J., Couch W. J., Ferré-Mateu A., 2023a, MNRAS, 518, 3653

Gannon J. S., Buzzo M. L., Ferré-Mateu A., Forbes D. A., Brodie J. P., Romanowsky A. J., 2023b, MNRAS, 524, 2624

Garcia-Dias R., Allende Prieto C., Sánchez Almeida J., Ordovás-Pascual I., 2018, A&A, 612, A98

Gordon K. D., Clayton G. C., Misselt K. A., Landolt A. U., Wolff M. J., 2003, ApJ, 594, 279

Harris W. E., Harris G. L. H., Alessi M., 2013, ApJ, 772, 82

Haslbauer M., Dabringhausen J., Kroupa P., Javanmardi B., Banik I., 2019, A&A, 626, A47