© The Authors 2026

1. Introduction

The stellar mass–size relation of galaxies encodes key information about the processes that regulate galaxy growth across cosmic time. Observational studies have shown that this relation evolves markedly with redshift, such that galaxies of a given stellar mass appear systematically more compact at earlier epochs and grow in size towards the present day (Shen et al. 2003; Trujillo et al. 2006; Somerville et al. 2008; van der Wel et al. 2014; Buitrago & Trujillo 2024). In particular, Buitrago & Trujillo (2024) report that Milky Way-like disc galaxies (M_* ∼ 5 × 10¹⁰ M_⊙) have increased their characteristic sizes by a factor of ∼2 over the last 8 Gyr, accompanied by a significant decrease in their stellar mass surface density at the disc edge. This evidence highlights the dramatic inside-out growth of discs, with an average radial expansion rate of ∼1.5 kpc Gyr⁻¹. This evolution is widely interpreted as the result of an interplay between in situ star formation, merger-driven mass assembly, and the influence of the dark matter (DM) halos in which galaxies reside.

Traditionally, galaxy sizes have been quantified using parameters such as the effective radius (R_e), defined as the radius enclosing half of the galaxy total light (de Vaucouleurs 1948; Graham et al. 2005). While widely adopted, these metrics can be sensitive to the depth of the observations, assumed light profile, and effects of inclination, making them vulnerable to significant scatter, especially for galaxies with irregular or extended morphologies (Capaccioli & de Vaucouleurs 1983; Trujillo et al. 2006; Mosleh et al. 2013; Costantin et al. 2023). More recently, alternative size definitions such as the R₁ radius–defined as the galactocentric distance where the stellar mass surface density reaches 1 M_⊙ pc⁻²–have been introduced as empirical proxies for the physical edge of in situ star formation (Trujillo et al. 2020). This threshold is motivated by the critical gas density required for star formation (Σ_c ∼ 3–10 M_⊙ pc⁻²; Schaye 2004) combined with typical star formation efficiencies, which yield stellar mass surface densities of order 1–3 M_⊙ pc⁻². For galaxies with stellar masses comparable to that of the Milky Way, this ‘star formation edge’ is indeed close to Σ_* ∼ 1 M_⊙ pc⁻² (Martínez-Lombilla et al. 2019). While the exact stellar mass surface density at the disc edge is expected to vary with both galaxy mass and morphology (Chamba et al. 2020, 2022; Buitrago & Trujillo 2024), adopting Σ_* = 1 M_⊙ pc⁻² provides a useful proxy for defining the size of the galaxy. This approach enables a practical and consistent comparison of galaxy sizes across different morphological types, particularly in the low surface brightness regime. Nevertheless, R₁ should not be regarded as a strict physical boundary, as its interpretation depends on the underlying morphology and is less meaningful for spheroid-dominated systems, whose growth is largely merger-driven.

Several recent studies have sought to investigate the physical origin of the mass–size relation and its scatter using cosmological simulations. For example, Du et al. (2024) and Ma et al. (2024) have highlighted the role of dark matter halo spin, concentration, and assembly history in regulating galaxy size. However, their analyses rely on conventional size metrics such as R_e, which may obscure or dilute underlying correlations due to intrinsic measurement uncertainties. More recently, Arjona-Gálvez et al. (2025) tested the R₁ size definition in state-of-the-art cosmological simulations, showing that rotation-dominated systems, such as bulgeless galaxies, tend to occupy the upper envelope of the mass–size relation.

Theoretically, in the framework of hierarchical structure formation, the formation of galactic discs is tightly coupled to the properties of their host dark matter halos. Following the classical picture of Fall & Efstathiou (1980), and later developments by Mo et al. (1998) and Somerville et al. (2008), the size of a disc is expected to scale with the virial radius and angular momentum of its halo, such that r_disc ∝ λr_vir, where λ is the halo spin parameter. This provides a simple physical basis for connecting the dark and luminous components of disc galaxies. However, while this classical picture offers a compelling first-order framework, it does not fully capture the complex role of mergers in shaping both the luminous and dark matter components of galaxies (Governato et al. 2009, 2010; Brook et al. 2011; Naab & Ostriker 2017). In particular, the existence of massive, bulgeless galaxies poses a challenge to the Λ cold dark matter (ΛCDM) paradigm. In a hierarchical Universe, galaxies grow through frequent mergers and accretion, the survival-or even the formation-of dynamically cold, pure stellar discs without central bulges is not straightforward (Hopkins et al. 2009). Mergers, especially those with high mass ratios, are expected to trigger central mass concentration through gas inflows and violent relaxation, promoting bulge growth (Naab et al. 2014). Yet, observational surveys find large, rotationally supported galaxies with no significant bulge, indicating that some systems may evolve through more quiescent or finely tuned assembly histories (Kautsch 2009; Barazza et al. 2009; Kormendy et al. 2010; Costantin et al. 2020).

Addressing this tension is one of the central goals of the Bulgeless Evolution And the Rise of Discs (BEARD, Méndez-Abreu et al., in prep.) survey. This project is designed to systematically characterise the structure, kinematics, stellar content, and formation pathways of a statistically significant sample of nearby bulgeless galaxies. BEARD is a multi-facility survey that combines deep broadband photometry, narrowband Hα imaging, and long-slit and integral field unit (IFU) spectroscopy to constrain the merger history and stellar mass growth of massive bulgeless galaxies. Combined with the study of accurately matched samples from numerical simulations, the survey aims to uncover the physical conditions that allow bulgeless galaxies to survive in a hierarchical Universe.

In this work, we present a pioneering analysis of a sample of nearby bulgeless galaxies from the deep images obtained as part of the BEARD survey. Deep g and r band observations are used to derive surface brightness, colour, and stellar mass density radial profiles, enabling the measurement of precise structural parameters such as R₁ and central stellar mass density, Σ_1, kpc, defined as the stellar mass surface density enclosed within the central kiloparsec. To place our observational results in a theoretical context, we also make use of simulated galaxies from the IllustrisTNG50 cosmological simulation (Nelson et al. 2019; Pillepich et al. 2019), enabling a direct comparison between observed and simulated structural properties. By applying the same analysis techniques to both observed and simulated datasets, we are able to perform a direct, one-to-one comparison of the structural parameters. We focus on the connection between the stellar structure and properties of the dark matter halo and merger history, exploring whether halo spin and accretion-driven processes can explain the diversity in disc extents at fixed stellar mass. This approach provides a robust framework to test whether current ΛCDM-based models can reproduce the existence and characteristics of bulgeless galaxies in the local Universe.

The paper is organised as follows. In Sect. 2, we describe the sample selection and observational data used in this study. Sect. 3 outlines the data reduction procedures and the methodology employed to derive structural parameters, including the R₁ radius and central stellar mass density. In Sect. 4, we present the stellar mass–size relation for bulgeless galaxies and investigate its connection to halo properties and merger histories using both observations and simulations. Finally, our main conclusions are summarised in Sect. 5.

Throughout this study, we assume a standard ΛCDM cosmology with Ω_m = 0.3, Ω_Λ = 0.7, and a Hubble constant of H₀ = 70 km s⁻¹ Mpc⁻¹. For the simulated galaxies from IllustrisTNG50, whose internal cosmology slightly differs from these values, all distances have been corrected to ensure consistency with the adopted cosmology. All magnitudes are given in the AB photometric system.

2. Sample and observations

The galaxies presented in this work were observed as part of the BEARD survey (Méndez-Abreu et al., in prep.), where the volume-limited BEARD parent sample consists of 75 spiral galaxies selected from the 13th Data Release of the Sloan Digital Sky Survey (SDSS-DR13) spectroscopic catalogue (Albareti et al. 2017) and HyperLEDA database (Makarov et al. 2014). The sample includes galaxies with stellar mass log(M_*/M_⊙) > 10 located within a distance of 40 Mpc. To select systems without prominent bulges, we applied a concentration cut of C = R₉₀/R₅₀ < 2.5 (Strateva et al. 2001; Graham & Driver 2005) to galaxies with low inclination (i < 60°), where R₉₀ and R₅₀ are the radii enclosing 90% and 50% of the total Petrosian flux, respectively. A further refinement of the sample, based on multi-component photometric decompositions (Zarattini et al., in prep.), was used to define the fiducial sample of 54 bulgeless Milky Way-like galaxies of the BEARD parent sample. In these decompositions – which include separate components for bulge, disc, bar, and unresolved central sources – we imposed B/T < 0.1 (bulge-to-total light ratio) and B/D < 0.08 (bulge-to-disc ratio), ensuring a robust selection of galaxies with negligible bulges.

To robustly estimate the location of the R₁ radius (Trujillo et al. 2020), information on the colour of the stellar population is necessary. For each galaxy, exposures were taken in both the g and r bands of the Sloan Digital Sky Survey (SDSS; York et al. 2000) system. The observations were obtained using the Wide Field Camera (WFC) mounted on the 2.5 m Isaac Newton Telescope (INT) at the Roque de los Muchachos Observatory, spanning three observing proposals between 2019 and 2021 (see Table A.1 in Appendix A). In total, 33 nights were allocated to the survey. The WFC consists of four CCDs arranged in an L-shaped mosaic, providing a total field of view of approximately 34 × 34 arcmins, with a pixel scale of 0.33 arcsecs/pixel. The data reach limiting depths in surface brightness ranging from μ_g = 27.7–30.6 mag arcsec⁻² and μ_r = 28.3–30.2 mag arcsec⁻². These measurements corresponded to 3σ fluctuations relative to the image background, using 10 × 10 arcsec² boxes, following surface brightness limit definition by Román et al. (2020).

From the volume-limited BEARD parent sample, a total of 40 galaxies were observed with this deep imaging strategy. However, we discarded those whose analysis was significantly affected by contamination from bright foreground stars in their outskirts and/or by the presence of galactic cirri, as both effects compromise the reliability of the surface brightness profile measurements. After this filtering process, we retained a final working sample of 22 galaxies, which are listed in Table A.2. Of these, 17 belong to the fiducial sample of 54 galaxies defined in Zarattini et al. (in prep.), while the remaining 5 do not.

3. Methodology

3.1. The IllustrisTNG50 simulated sample

For the simulated comparison, we adopted the same galaxy sample described in Cardona-Barrero et al. (2025), extracted from the IllustrisTNG50-1 run (Pillepich et al. 2019; Nelson et al. 2019). This simulation follows the evolution of dark matter, gas, and stars in a (51.7 Mpc)³ volume, with a dark matter particle mass of 4.5 × 10⁵ M_⊙ and baryonic mass resolution of 8.5 × 10⁴ M_⊙. The assumed cosmology is consistent with Planck Collaboration XXVII (2016), but distances were corrected a posteriori to match the cosmology adopted in this work.

The high spatial resolution of IllustrisTNG50, with a gravitational softening length of 575 comoving pc until z = 1 and fixed to 288 physical pc at lower redshifts, enables a detailed study of galaxy structure. We selected central galaxies with stellar masses 10¹⁰ < M_* < 5 × 10¹¹ M_⊙, applying an isolation criterion that excludes companions of comparable mass within 0.5 Mpc. This resulted in a sample of 537 simulated galaxies. For further details on the simulation setup and sample selection, we refer the reader to Cardona-Barrero et al. (2025).

3.2. Data reduction

Image processing was carried out with the purpose of providing the highest reliability in the low surface brightness regime. For data reduction, we separated between image sets. Each set of images was formed by the blocks of nights indicated in Table A.1, which are considered to have a uniform and similar flat-field between them, due to the temporal proximity of the observations.

Each set of images was bias-subtracted by standard procedures. For flat-fielding, we heavily masked the bias-subtracted images using SExtractor (Bertin & Arnouts 1996) and NoiseChisel (Akhlaghi & Ichikawa 2015). The coaddition of these masked images after normalisation provided a high-quality night-sky flat. Once this night sky flat was obtained, the images were flat-fielded and reduced. This procedure was performed independently for the g and r bands.

The process of obtaining the coadds – final stacked images created by combining multiple individual exposures to increase depth and signal-to-noise (S/N) – for each galaxy and each band was performed as follows. First, we computed astrometric solutions by using Astrometry.net (Lang et al. 2010) to obtain an initial solution, which was then refined with SCAMP (Bertin 2006) for greater accuracy. The next step was to produce a starting seed coadd, which served as the first step in an iterative process. This first coadd was made by stacking all the images of the field, performing a sky subtraction with a constant value. This preserved all the information in the low surface brightness regime, however, strong gradients were present in it. To provide a more efficient sky subtraction and to remove residual gradients while preserving the low surface brightness information we performed an iterative process. This process consisted of obtaining a mask from this coadd, and applying this mask to all individual exposures. Once the individual exposures were masked, we calculated the sky background using Zernike polynomials (Zernike 1934). This sky background was subtracted from all individual exposures, which after photometric calibration referenced to SDSS, were combined to produce a new coadd. This process was repeated several times until a satisfactory coadd was obtained in which the residual gradients were minimised, while preserving the actual flux of the sources at low surface brightness. Typically, we used orders of Zernike polynomials between 2 and 4, depending on the intensity of the gradients present in the images. These gradients are variable between epochs and are especially intense in moderate moonshine conditions.

Similar coaddition processes have been carried out in the past (Sánchez-Alarcón et al. 2023; Román et al. 2023; Pérez et al. 2024; Jiménez-Teja et al. 2025), with highly efficient results. We thus made sure to avoid oversubtraction problems around sources, being particularly problematic the halos around galaxies of large apparent size, as the case of BEARD galaxies.

3.3. Obtaining the extended PSF of the complete sample

It is well established that the PSF plays a critical role in the analysis of low surface brightness structures (Capaccioli & de Vaucouleurs 1983; Sandin 2014; Spavone et al. 2020; Infante-Sainz et al. 2020; Sedighi et al. 2025). In particular, the outskirts of galaxies may be significantly affected by scattered light originating from nearby stars as well as from the central regions of the galaxies themselves. Accurate PSF characterisation is therefore essential to ensure the reliability of the derived results. To this end, we followed a methodology consistent with previous works (Infante-Sainz et al. 2020; Sedighi et al. 2025; Marrero-de la Rosa et al., in prep.). Rather than adopting a single, static PSF model, we constructed a tailored PSF for each observed field, accounting for variations in atmospheric seeing and instrumental response across the observing campaign. The PSF was constructed by stacking isolated stars from each field, grouped into concentric radial zones to accurately trace both the core and the extended wings. Unsaturated stars were used to define the inner PSF profile, while progressively brighter (and even saturated) stars were employed at larger radii, where their high flux allows us to recover the faint outer wings with improved S/N despite saturation in the central pixels. This multi-region strategy allowed us to recover the PSF structure out to large radii (> 1000 arcsec), minimising background contamination and maximising radial coverage. As a result, a total of 22 distinct PSF profiles were derived in each filter (g and r), one for each galaxy in the sample (Fig. 1).

Fig. 1. Normalised radial profile of the PSFs for all fields observed in the g and r bands for the sample galaxies. The PSF total flux is scaled to unity. Vertical shaded regions indicate structural zones defined in the galaxy profile analysis (Inner, Subintermediate, Intermediate, Outer and Outer Ext) (see Sect. 3.3 for details). The y-axis represents normalised intensity in arbitrary units.

To optimise both spatial resolution and S/N at all radii, the radial characterisation of the PSF was divided into five structural regions, following a procedure analogous to the one that will be presented in Marrero-de la Rosa et al. (in prep.). Each of these radial zones dominates the PSF profile in S/N within its respective range and boundaries between regions were chosen to ensure that adjacent profiles overlap at comparable S/N levels, guaranteeing smooth transitions across the full radial extent (Infante-Sainz et al. 2020). While stellar magnitudes are expressed in the GaiaG-band (Gaia Collaboration 2016, 2023), the PSFs were built from stars selected directly within each observed galaxy field, ensuring consistency with the specific observing conditions. The five regions, highlighted in Fig. 1 with distinct colours, are defined as follows:

• Inner (0–3.3 arcsec, highlighted in green): constructed using stars with magnitudes between 19 < G < 20 mag.

• Subintermediate (3.3–5.8 arcsec, in yellow): based on stars between 17 < G < 18 mag.

• Intermediate (5.8–9.9 arcsec, in salmon): based on stars between 15 < G < 16 mag.

• Outer (9.9–115.5 arcsec, in orange): modelled using the brightest field stars (G < 15 mag), subdivided into five bins of magnitude to produce multiple independent stacks and improve robustness.

• Outer External (> 115.5 arcsec, in grey): constructed using a single bright star (HD 8648) of G = 7.24 mag located in the field of NGC 521. Such a bright source is essential to reliably trace the faintest outskirts of the PSF profile, where the signal from fainter stars becomes indistinguishable from the background noise.

Each radial region overlaps with the next within transition zones defined to ensure that both profiles share comparable S/N levels at the overlapping radius. This guaranteed smooth and continuous PSF profiles across the full dynamic range. The final PSFs resulted from the sequential stitching of these regions, calibrated to match flux at the overlap boundaries, and validated to avoid background oversubtraction or residual artefacts. Details of the corrections applied for scattered light from foreground stars, for self-scattered light from the galaxies themselves, and for masking contaminant sources and subtracting the local background can be found in Appendices B, C and D, respectively.

3.4. Radial surface-brightness, colour, and stellar mass profiles

Once all necessary image corrections have been applied, reliable surface brightness profiles can be derived. To ensure the robustness of these profiles, only the regions of the elliptical isophotes that intersect the wedge profiles–previously used for local background estimation–were considered (see Appendix D for details). This procedure was carried out independently for both the g- and r-band images.

To place all galaxies on a common footing and mitigate projection effects, the surface brightness profiles were corrected for inclination following the prescription from Trujillo et al. (2020). This correction assumes that face-on profiles are systematically fainter than inclined ones due to geometric projection, and is implemented via a polynomial adjustment of the form

$$\begin{array}{c}\hfill \mathrm{\Delta }\mu ={\displaystyle \sum _{j=0}^4}{\alpha }_j{(b/a)}^j,\end{array}$$(1)

where Δμ is the inclination correction in mag arcsec⁻², b/a is the observed axial ratio of the galaxy, and the coefficients α_j are determined empirically for a given z₀/h ratio, with z₀ and h being the vertical and radial scale lengths of the stellar disc, respectively. For this work, we adopted the coefficients corresponding to z₀/h = 0.12, as listed in Table 1 of Trujillo et al. (2020), which are representative of the late-type galaxies in our sample.

Additionally, all profiles were corrected for foreground Galactic extinction using the A_γ values provided by the NASA NED Foreground Galactic Extinctions service (Schlafly & Finkbeiner 2011). These extinction corrections were independently applied for each photometric band and galaxy, ensuring a consistent comparison of intrinsic galaxy properties across the full sample. After applying both the inclination and extinction corrections, the final corrected surface brightness profiles used throughout this work are given by the following:

$$\begin{array}{c}\hfill {\mu }_{\mathrm{corr},\gamma }={\mu }_{\mathrm{obs},\gamma }-A_{\gamma }+\mathrm{\Delta }\mu .\end{array}$$(2)

This correction allows all galaxies to be consistently analysed as if observed face-on. From the inclination-corrected profiles, both the (g − r)_o colour profiles and the stellar mass surface density profiles (Σ_*) were computed (Fig. 2).

Fig. 2. Left panel: colour image of NGC 3486. A 1 arcmin scale bar is shown, together with the area used to extract the profiles, highlighted as a wedge. The R1 radius is shown as a cyan dashed-lined ellipse. Middle panel: surface brightness profiles in the g (dashed blue) and r (solid red) bands, along with the corresponding stellar mass surface density profile (black points). The vertical and horizontal dashed grey lines marks the R1 radius and 1 M⊙ pc−2, respectively. Right panel: Colour profile (g − r)o corrected from extinction.

To estimate the stellar mass surface density, we followed the approach described in Bakos et al. (2008) and Roediger & Courteau (2015), using the relation

$$\begin{array}{c}\hfill {\log }_{10}{\mathrm{\Sigma }}_{\ast }={\log }_{10}{(M/L)}_{\gamma }-0.4({\mu }_{\mathrm{corr},\gamma }-m_{\odot ,\gamma }^{\mathrm{abs}})+8.629,\end{array}$$(3)

where μ_corr, γ is the corrected surface brightness in a given γ band, $m_{\odot ,\gamma }^{\mathrm{abs}}$ is the absolute magnitude of the Sun in that band, and (M/L)_γ is the mass-to-light ratio calculated as

$$\begin{array}{c}\hfill {\log }_{10}{(M/L)}_{\gamma }=a_{\gamma }+b_{\gamma }\times \mathrm{colour}.\end{array}$$(4)

In analogy with the methodology of Chamba et al. (2022), we adopted γ = g as it provides the deepest photometric data. The colour term used was (g − r)_o, along with coefficients a_g = −0.984 and b_g = 2.029 taken from Roediger & Courteau (2015). The absolute magnitude of the Sun is $m_{\odot ,g}^{\mathrm{abs}}=5.11$ from Willmer (2018).

The R₁ radius was computed along the elliptical isophote where the stellar mass surface density reaches a value of 1 M_⊙ pc⁻², with the galactocentric distance measured along the semi-major axis, i.e. not projected onto a wedge profile. This metric serves as an empirically motivated proxy for the extent of the stellar component in low surface brightness regimes, and it is particularly useful for comparing galaxies across a wide range of morphologies in both observations and simulations (Trujillo et al. 2020; Arjona-Gálvez et al. 2025). An example of the measurement of R₁, highlighted as a vertical grey line, can be seen in Fig. 2.

For the simulated galaxies from the IllustrisTNG50 project (Pillepich et al. 2019; Nelson et al. 2019), the R₁ radius was derived directly from the stellar mass density profiles. These were computed after centring each galaxy using the shrinking-sphere method of Power et al. (2003), reorienting the system so that the stellar disc is face-on with respect to its angular momentum vector, and then measuring the stellar mass density in logarithmically spaced radial bins. The R₁ radius was identified as the radius at which this profile reaches 1 M_⊙ pc⁻², ensuring consistency with the observational definition while accounting for the three-dimensional nature of the simulation.

The total stellar mass of each galaxy, M_*, was obtained by integrating the stellar mass surface-density profile out to the radius where the g-band surface brightness profile reaches a limiting value of 29 mag arcsec⁻². This threshold was adopted to ensure consistency with the methodology employed in Trujillo et al. (2020) and it represents a conservative limit, beyond which the photometric uncertainties become significant. In those particular cases where the surface brightness limit of 29 mag arcsec⁻² was not reached, the stellar mass was integrated up to the R₁ radius, which still encompasses the bulk of the galaxy stellar mass. Once the total mass had been computed using these criteria, the stellar mass–size relation could be derived.

3.5. Environmental characterisation

The potential impact of the surrounding environment on the scatter of the stellar mass–size relation was investigated by quantifying the local galaxy density for both the BEARD and IllustrisTNG50 samples. As a first step, we adopted the projected surface density to the fifth nearest neighbour (Σ₅), defined as ${\mathrm{\Sigma }}_5=5/(\pi r_5^2)$, where r₅ represents the projected distance to the fifth closest galaxy. This metric follows standard practice and has been employed in previous works (e.g. Aguerri et al. 2009). Its specific application to the BEARD sample will be described in detail in the forthcoming survey paper Méndez-Abreu et al. (in prep.).

For the BEARD sample, Σ₅ was calculated using two complementary estimates: (i) a photometric version, based on SDSS imaging, and (ii) a spectroscopic version employing redshift information from SDSS-DR13 (Albareti et al. 2017) to minimise contamination from background or foreground sources. In both cases, neighbouring galaxies were required to have magnitudes within ±2 mag in the r band relative to the host galaxy, ensuring that only systems of comparable luminosity were considered. The same magnitude selection was applied to the simulated galaxies when identifying neighbours within the IllustrisTNG50 volume.

The simulated galaxies were analysed using the same definition of Σ₅, applied to all central systems in the IllustrisTNG50 run. It is worth noting that these simulated galaxies were selected according to an isolation criterion that excludes companions of similar stellar mass within a projected radius of 0.5 Mpc. This selection ensures that the simulated sample resembles the typically isolated nature of the majority of BEARD galaxies, although it may bias the environmental characterisation towards lower-density regimes.

In addition to Σ₅, we also computed the instantaneous tidal radius (r_{tidal, inst}) for the first five nearest neighbours of each central galaxy, defined as

$$\begin{array}{c}\hfill r_{\mathrm{tidal},\mathrm{inst}}=D\times {\left[\frac{m_{\mathrm{sat}}}{3M_{\mathrm{host}}}\right]}^{1/3},\end{array}$$(5)

where D is the instantaneous separation between galaxies and a circular orbit (e_orb = 0) is assumed. This formulation follows the classical tidal approximation (e.g. King 1962; Johnston et al. 2002), in which the ratio between the satellite and host enclosed masses determines the strength of the tidal field. Similar approaches have been applied in recent studies of diffuse tidal features and low-surface-brightness systems (e.g. Montes et al. 2021). In this work, r_{tidal, inst} was computed using the stellar masses and separations of galaxies in the IllustrisTNG50 simulation, considering up to the five closest neighbours of each BEARD analogue.

4. Results and discussion

In this section, we explore the stellar mass–size relation using R₁ as a physically motivated size indicator. Fig. 3 shows the stellar mass–size relation for the 22 galaxies from our BEARD sample (light-green squares). To place our results in a broader context and improve completeness, we include a supplementary dataset of 1005 galaxies compiled by Trujillo et al. (2020), shown as grey points. This ancillary sample comprises 279 ellipticals, 464 spirals, and 262 dwarfs, spanning a stellar mass range of 10⁷ < M_* < 10¹² M_⊙ and redshifts between 0.01 < z < 0.1. A detailed description of the sample selection and analysis can be found in the original publication. Although the BEARD and Trujillo et al. (2020) samples differ in certain aspects of the data reduction and analysis, both followed broadly similar procedures to those outlined in Sect. 3.2. For the BEARD sample, the reductions were carried out using an updated version of the pipeline, and an additional wavelet-based deconvolution was applied to correct for self-scattered light from the galaxies themselves. In the overlapping mass range, the two datasets show good agreement.

Fig. 3. Stellar mass–size relation defined by the R1 radius. The main panel shows the global distribution of galaxies from different datasets: light-green squares show the BEARD sample of observed galaxies, grey points correspond to the observational sample from Trujillo et al. 2020 (TCK20), and cyan circles represent simulated galaxies from the IllustrisTNG50 run of the IllustrisTNG project (Cardona-Barrero et al. 2025). The right panels separate the IllustrisTNG50 sample into three categories based on their dynamically decomposed bulge-to-disc ratio: BL galaxies in green, BI galaxies in blue, and BD galaxies in orange. To compute the coloured lines for each morphological class, the stellar mass range was divided into 11 equally spaced bins; in each bin, the average R1 value was calculated, and the shaded regions indicate the corresponding standard deviation. For comparison, the TCK20 (grey dots) and BEARD (light-green squares) samples are also shown in each panel.

Despite the fact that using R₁ as a size proxy yields a typical scatter in the mass–size relation that is approximately 2.5 times lower than that obtained with conventional size estimators based on R_e (Trujillo et al. 2020), there is a remaining scatter in the mass–size relation. In order to investigate it, we complement our analysis with a subsample of galaxies from the IllustrisTNG50 cosmological simulation (Cardona-Barrero et al. 2025). The aim is to provide a theoretical framework capable of interpreting the location and dispersion of the observed bulgeless galaxies within the mass–size plane (cyan circles, see Sect. 3.1).

Simulated galaxies are classified into three categories based on their bulge-to-disc (B/D) ratios: bulgeless (BL; also BEARD-Like, 134 galaxies), bulge-intermediate (BI, 269 galaxies), and bulge-dominated (BD, 134 galaxies). This classification is established by dividing the B/D distribution into quartiles, with thresholds at B/D = 0.057 and B/D = 0.383, where the bulge and disc components are dynamically identified through their orbital properties (Cardona-Barrero et al. 2025). This approach allows us to directly compare the structural diversity of galaxies across a continuous range of morphologies and to link the observed scatter with underlying physical parameters, such as halo properties and merger histories.

The BEARD sample is composed exclusively of bulgeless galaxies with B/D ≲ 0.08. This allows for a direct comparison with the BL population in simulations and provides a stringent observational benchmark for the formation and structure of pure disc galaxies.

4.1. The mass–size relation of Milky Way analogues

Within the stellar mass range covered by our sample of IllustrisTNG50 simulated spiral galaxies, a clear structural trend emerges in the mass–size relation (Fig. 3). Galaxies classified as BL predominantly occupy the upper envelope of the relation, while BD systems tend to cluster on the lower side. Although this overall separation is evident, the dispersion bands shown in Fig. 3 partly overlap. This indicates that, at fixed stellar mass, BL galaxies are generally more extended than their BD counterparts, introducing an intrinsic scatter in the relation. Such a trend aligns well with the results from numerical simulations presented by Arjona-Gálvez et al. (2025), where the fastest rotating galaxies are preferentially located in the upper part of the mass–size diagram.

To quantify these trends, we performed linear fits to the mass–size relation for BEARD, BL, BI, and BD galaxies (Table 1; see Fig. E.1 in Appendix E). For the IllustrisTNG50 sample, there are clear differences in both slope and intercept among the three morphological categories: BL galaxies define the upper envelope with shallower slopes and larger radii at fixed stellar mass, BD galaxies occupy the most compact regime with the steepest slope, and BI systems trace an intermediate sequence. When comparing BL galaxies in IllustrisTNG50 and BEARD, the fitted slopes are fully consistent, but a small offset in intercept is evident, with BEARD galaxies lying slightly below the simulated BL sequence. This offset is reinforced by the presence of a number of IllustrisTNG50 BL galaxies that deviate from the observed mass–size relation, lying outside the region populated by both the BEARD and Trujillo et al. (2020) samples.

To further test the consistency between BEARD and IllustrisTNG50, we restricted both samples to the stellar mass range 10.0 < log(M_*/M_⊙) < 10.5 and generated 10⁴ random realisations of the BL subsample, each time drawing the same number of galaxies as in BEARD within this range. The resulting distribution of scatters is shown in Fig. E.2 in Appendix E. The observed BEARD scatter (σ_R1 = 0.11 ± 0.02) is fully consistent with the mean scatter of BL galaxies in IllustrisTNG50 (0.11 ± 0.02). Kolmogorov–Smirnov (KS) tests (Stephens 1974; Fasano & Franceschini 1987), however, indicate that the distributions are statistically distinct: for the 1D comparison in log R₁, we find p = 5.2 × 10⁻³, and for the 2D (log M_*, log R₁) plane, p = 1.1 × 10⁻³. These values suggest that, while the level of scatter is comparable, the detailed distributions of the observed and simulated galaxies remain different.

Our BEARD galaxies do not preferentially lie on the upper envelope of the mass–size relation; rather, they span a broad range of R₁ values at fixed stellar mass, covering the full vertical extent of the relation. The error bars shown in Fig. 3 are computed by combining two main sources of uncertainty: first, the sensitivity of R₁ to the local background estimation in the wedge profiles – quantified by adding or subtracting the measured background uncertainty to the surface brightness profiles – and second, the uncertainty in stellar mass estimation. The latter follows the prescription of Trujillo et al. (2020), who compared their stellar mass estimates (derived using g − r colour profiles and the M/L relation of Roediger & Courteau 2015) with independent values from the Portsmouth SED-fitting catalogue (Maraston et al. 2013). For massive spiral galaxies analogous to those in our sample, they report an rms scatter of 0.24 dex, which we adopt here as a representative uncertainty in stellar mass. This choice is justified by the similarity of the adopted M/L prescription and by the fact that both studies compute stellar masses by integrating the stellar mass surface density up to the limit of 29 mag arcsec⁻² in the g-band. Despite the intrinsic tightness of the mass–size relation when using R₁, these combined error bars do not fully account for the observed vertical scatter in most cases.

Therefore, in the following sections, we investigate potential physical drivers of this dispersion. First, we examine the connection between galaxy size and inner structure by measuring the stellar mass surface density within the central kiloparsec (Σ_1, kpc) for both the BEARD observations and the IllustrisTNG50 simulations. We then turn to the simulated sample to explore how global halo properties – including total mass, concentration and spin parameter – may influence galaxy size. Finally, we analyse the assembly histories of the simulated galaxies, focusing on their merger rates and gas accretion patterns, to evaluate their impact on the structural diversity of present-day discs.

4.2. Connecting the central and outer parts of the galaxy

A commonly used parameter to trace the central structure of galaxies is Σ_1, kpc. This parameter has been shown to be a robust tracer of the internal structure of galaxies and their evolutionary stage (Cheung et al. 2012; van Dokkum et al. 2014; Tacchella et al. 2015; Lee et al. 2018; Luo et al. 2020). One of the key advantages of using Σ_1, kpc is that, unlike parameters that rely on bulge–disc decompositions or high-resolution imaging, it can be consistently measured across both observations and simulations using straightforward aperture photometry.

However, to avoid the built-in correlation between total stellar mass and central mass density (Luo et al. 2020), we normalise Σ_1, kpc by the total stellar mass of each galaxy, defining a specific central stellar mass density ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}\equiv {\mathrm{\Sigma }}_{1,\mathrm{kpc}}/M_{\ast }$. This allows us to better isolate structural differences between galaxies independently of their stellar mass. Unlike our approach, which removes the stellar-mass dependence through this normalisation, Luo et al. (2020) describe the correlation between Σ_1, kpc and M_* using a quadratic polynomial fit in log–log space.

When connecting the inner structure of galaxies, traced by ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$, with their global properties, we observe clear differences as a function of bulge prominence (Fig. 4). BD galaxies follow a well-defined trend in both panels: in the left panel, ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ decreases with stellar mass with a slope of −0.39 and a Pearson correlation coefficient of r = 0.73, while in the right panel the correlation with R₁ is similarly strong, with a slope of −0.68 and r = 0.77. This highlights the tight link between central density and overall structure for bulge-dominated systems. By contrast, BL galaxies do not show a dominant correlation, with their central densities contributing mainly to the scatter of the mass–size relation. Interestingly, our BEARD galaxies exhibit a stronger correlation between ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and R₁ (slope −1.43, r = 0.82) than with stellar mass (slope −0.50, r = 0.62), suggesting that for bulgeless discs the extent of the system is more closely connected to its central density. This result, however, may be affected by completeness biases in our observational sample. A more detailed exploration of these scenarios is presented in Sect. 4.5.

Fig. 4. Relation between the specific central mass density ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and the global properties of galaxies. Left panel: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ as a function of stellar mass. Right panel: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ versus R1. In both panels, simulated galaxies from IllustrisTNG50 are colour-coded according to their bulge prominence, with BD galaxies (orange diamonds) occupying the upper regions of the relations, while BL systems (green circles) lie systematically below. BI galaxies (blue crosses) span a wide range in both ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and size, bridging the two extremes. Observational data from the BEARD sample are shown as light-green squares, consistently occupying the region where BL galaxies are located.

4.3. The impact of the dark matter halo

To explore how the properties of dark matter halos influence the internal structure of spiral galaxies, we examine the relations between the halo mass ($M_{200}^{\mathrm{DM}}$), spin parameter (λ; defined as λ = (J |E|^1/2)/(G M^5/2), where J, E, and M are the total angular momentum, energy, and mass of the system, and G is the Newtonian constant of gravitation, following Bullock et al. 2001), concentration ($c_{200,\mathrm{c}}^{\mathrm{DM}}$), and key stellar structural properties such as ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and R₁, using the simulated galaxies from IllustrisTNG50. Fig. 5 reveals multiple bimodal trends across this parameter space when these galaxies are separated by morphology.

Fig. 5. Relations between stellar and dark matter properties for the simulated galaxies from IllustrisTNG50. The figure is organised as a matrix, where each row corresponds to a stellar quantity (from top to bottom: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$, M* and R1) and each column corresponds to a dark matter halo property (from left to right: λ, $M_{200}^{\mathrm{DM}}$, and $c_{200,\mathrm{c}}^{\mathrm{DM}}$). BL, BI, and BD galaxies are shown as green circles, blue isodensity contours, and orange diamonds, respectively. On the top and right panels, histograms of the corresponding quantities are shown for BL and BD galaxies, where vertical lines represent the mean of each distribution.

A strong correlation is found between ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and $M_{200}^{\mathrm{DM}}$ for BD galaxies, which follow a tight sequence (panel b; Fig. 5). BL and BI galaxies tend to deviate from this trend, populating the lower-density regime at fixed halo mass. This behaviour mirrors the correlation between ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and M_*. When comparing ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ with λ (panel a; Fig. 5), the negative correlation for BD galaxies becomes less clear due to larger scatter. However, a clear separation persists: at fixed spin, BD galaxies exhibit significantly higher central densities than BL systems. A similar pattern is observed at fixed $c_{200,\mathrm{c}}^{\mathrm{DM}}$, where BD galaxies again present systematically higher ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ values, suggesting that spin, concentration and mass jointly modulate the ability of halos to support dense stellar cores (panel c; Fig. 5).

Regarding stellar mass versus halo mass, a correlation with $M_{200}^{\mathrm{DM}}$ is evident (panel e; Fig. 5), with a small offset, that would result in a slight difference in the stellar-to-halo mass, similar to what we found in Rosas-Guevara et al. (2026). There is no robust trend for stellar mass versus concentration. But, analogously to the findings of Proctor et al. (2025), we also observe that halos with ${log}_{10}(M_{200}^{\mathrm{DM}}/M_{\odot })\gtrsim 12.5$ are predominantly populated by BD galaxies. Interestingly, the stellar mass versus spin relation shows that between 10.0 < log₁₀(M_*/M_⊙) < 10.5, BL galaxies tend to have higher spin values at fixed stellar mass.

We find a strong positive correlation between the R₁ radius and $M_{200}^{\mathrm{DM}}$, indicating that galaxy size increases with halo mass across the sample. Notably, at fixed $M_{200}^{\mathrm{DM}}$, BL galaxies tend to exhibit systematically larger R₁ values compared to their BD counterparts (panel h; Fig. 5). In contrast, neither the R₁–$c_{200,\mathrm{c}}^{\mathrm{DM}}$ nor the M_*–$c_{200,\mathrm{c}}^{\mathrm{DM}}$ relation display strong global trends, though morphological differences remain (panels i and f; Fig. 5).

These findings are consistent with both theoretical and observational works. Simulations by Rodriguez et al. (2024) demonstrate that, at fixed stellar mass, disc-dominated galaxies tend to reside in lower-mass halos than spheroid-dominated ones. In addition, Ma et al. (2024) report a clear link between the dynamical support of simulated disc galaxies and angular momentum of their host halos, highlighting the role of halo spin in sustaining rotationally supported structures. Observational studies support these trends through satellite counts, gravitational lensing, and dynamical tracers such as HI and globular clusters (e.g. Wang & White 2012; Mandelbaum et al. 2016; Posti et al. 2019). Altogether, our analysis reinforces the view that both the mass and spin of dark matter halos-shaped in part by their assembly history-are fundamental in explaining the diversity of internal galaxy structures, especially the dichotomy between systems with and without prominent bulges.

4.4. The role of merger history

To assess the impact of merger-driven processes on the internal structure of spiral galaxies, we analyse the relations between ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$, M_* and R₁ with respect to merger-related quantities: the mean ‘cold’ (i.e. star-forming) gas fraction of all the objects that have merged with the galaxy, weighted by their maximum stellar masses, total number of mergers, both major (stellar mass ratio 1:4) and minor mergers since z = 5, and accreted stellar mass fraction since z = 5. All the information was taken from the publicly available IllustrisTNG50 catalogues (Rodriguez-Gomez et al. 2017; Eisert et al. 2023). Fig. 6 reveals several notable trends, many of which exhibit clear morphological bimodalities.

Fig. 6. Relations between stellar and merger properties for the simulated galaxies from IllustrisTNG50. The figure is organised as a matrix, where each row corresponds to a stellar quantity (from top to bottom: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$, M* and R1) and each column corresponds to a merger property (from left to right: mean gas fraction, total mergers since z = 5 and accreted mass fraction). BL, BI, and BD galaxies are shown as green circles, blue isodensity contours, and orange diamonds, respectively. On the top and right panels, histograms of the corresponding quantities are shown for BL and BD galaxies, where vertical lines represent the mean of each distribution.

When comparing ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ with the mean gas fraction (panel a; Fig. 6), BD galaxies tend to exhibit lower gas fractions and higher central densities, while BL galaxies consistently occupy the lower-density, higher-gas-fraction regime. In particular, we do not find BL galaxies associated with mergers that had little cold gas, which is consistent with the idea that gas-poor (dry) mergers tend to produce spheroidal systems (Naab et al. 2009). At higher gas fractions, both BD and BL morphologies coexist; however, BL systems remain systematically less centrally concentrated than their BD counterparts. This suggests that central compaction likely plays a key role in building dense inner regions, particularly in gas-poor merger events. The relation between Σ_1, kpc^spec and the total number of mergers (panel b; Fig. 6) echoes the behaviour seen with stellar and halo mass: BD galaxies follow a tight sequence where central density decreases with merger count, whereas BL systems show no clear trend and cluster towards lower Σ_1, kpc^spec values. A similar dichotomy is seen when comparing Σ_1, kpc^spec with the accreted stellar mass fraction (panel c; Fig. 6): BD galaxies show a slight negative trend, while BL galaxies exhibit a positive correlation. This may indicate different evolutionary paths, with mergers playing contrasting roles depending on morphology.

In the M_*-mean gas fraction plane (panel d; Fig. 6), no strong correlations are apparent, but a sharp morphological transition is evident: galaxies with log₁₀(Mean Gas Fraction)≲ − 0.6 are exclusively BD, highlighting the role of gas depletion in bulge formation or preservation. When examining stellar mass versus total merger count (panel e; Fig. 6), we find a clear positive correlation, with BD galaxies being more relevant at log₁₀(M_*/M_⊙)≳11 and total mergers ≳100. However, no significant trends emerge when stellar mass is plotted against the accreted mass fraction.

In terms of size, R₁ shows similar behaviour to M_* when plotted against the mean gas fraction (panel g; Fig. 6). However, the R₁-total mergers plane (panel h; Fig. 6) reveals an important distinction: at fixed merger count, BL galaxies tend to be more extended than their BD counterparts. This supports the idea that merger-driven evolution contributes not only to central concentration but also to the scatter in the mass–size relation. A weak correlation is observed between R₁ and the accreted mass fraction (panel i; Fig. 6), although no clear separation by morphology is found.

These results highlight the complex interplay between merger activity and galaxy structure. Our findings are broadly consistent with prior simulation-based studies that emphasise the importance of mergers in shaping both the morphology and compactness of galaxies. For instance, Genel et al. (2015) and Rodriguez et al. (2024) find that galaxies with richer merger histories develop more prominent bulges and denser cores. Proctor et al. (2025) show that discs subjected to higher merger mass ratios become thicker and more extended, whereas more quiescent histories promote the survival of thin discs. Similarly, Di Cintio et al. (2019) demonstrate that the orbital alignment of the last major merger can result in either high- or low-surface-brightness systems, producing two distinct morphologies at fixed stellar mass. Finally, Jiang et al. (2025) report that giant BL galaxies at high redshift experience quieter merger histories compared to bulge-bearing systems. Furthermore, Martig et al. (2009) suggest that mergers can also induce compaction in gas-rich systems, leading to high central densities. Our results support this view while extending it by showing how the role of mergers varies across morphology and may be a key driver of the observed scatter in structural scaling relations. The relation between merger history and dark matter halo properties has been further explored and is discussed in Appendix F.

4.5. Intrinsic scatter within the bulgeless population

While the stellar mass–size relation exhibits a bimodality between BD and BL galaxies, a level of intrinsic scatter is also present within the BL population itself. As shown in Fig. 3, most of the simulated BL galaxies occupy the upper envelope of the relation, exhibiting large R₁ values at fixed M_*. This population is in good agreement with predictions from Arjona-Gálvez et al. (2025), where rotation-dominated galaxies are expected to exhibit more extended stellar discs due to their higher angular momentum content (Du et al. 2024).

However, Fig. 3 also shows a significant number of BEARD BL galaxies falling below the simulated BL population. Besides any possible systematics in the comparison between observations and simulations, we also find a non-negligible fraction of BL systems in the IllustrisTNG50 sample that lie below the mean mass–size trend defined by BL galaxies (green line in Fig. 3). This raises the question of what physical processes could account for such diversity within a morphologically homogeneous class. To investigate this, we selected a subpopulation of simulated BL galaxies from IllustrisTNG50 that systematically fall below the average R₁ value in bins of stellar mass by more than 1σ. A visual representation of this selected subsample can be found in Appendix G, where its location in the mass–size plane is shown more clearly.

The analysis of this subset reveals several key trends (Fig. 7). Firstly, these ‘below-trend’ BL galaxies exhibit, by definition, systematically smaller R₁ radii, as expected from the mass–size relation. Secondly, they show higher ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ values, similar to those found in BD systems, confirming their more compact nature. This configuration – characterised by smaller R₁ and higher central densities – mirrors the correlation between ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and R₁ already identified in the BEARD galaxies (Sect. 4.2), reinforcing the view that central density is also a key factor behind their location in the mass–size plane. Finally, their mean gas fractions are slightly higher than those of typical BL galaxies, suggesting that these systems may be in an earlier or more gas-rich phase of disc evolution. This latter result is consistent with previous simulations where gas-rich mergers tend to create more compact and concentrated galaxies (Naab et al. 2009).

Fig. 7. Normalised distributions of three structural parameters for the simulated BL galaxies from IllustrisTNG50 (grey) and simulated BL galaxies that lie significantly below the average mass–size trend (hatched cyan). Left: R1 radius,. Middle: central stellar mass density within 1 kpc, ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$. Right: mean gas fraction since z = 5.

These findings point towards a complex and heterogeneous evolutionary pathway for BL galaxies. The presence of compact, gas-rich discs with high central densities in the BL population may reflect different angular momentum acquisition histories, residual effects of minor mergers, or ongoing secular evolution. Disentangling these scenarios will require a more detailed dynamical characterisation of their stellar and gaseous components.

4.6. Environmental trends

Fig. 8 shows the distribution of Σ₅ for the BEARD and IllustrisTNG50 samples. The simulated populations (BL, BI, and BD) display no clear morphological segregation in Σ₅, suggesting that large-scale environment does not play a dominant role in setting the structural differences among simulated galaxies. For the BEARD sample, the photometric estimate of Σ₅ tends to overpredict local densities, most likely due to projection effects, whereas the spectroscopic version yields values more consistent with those derived from the simulations. A small secondary peak is observed in the spectroscopic BEARD distribution, corresponding to galaxies located within the Virgo Cluster region, naturally reflecting their denser environments.

Fig. 8. Distribution of the projected environmental density parameter, Σ5, for the simulated BL, BI, and BD galaxies compared with the BEARD-Phot and BEARD-Spec samples.

The distributions of the instantaneous tidal radius, log₁₀(r_{tidal, inst}), are presented in Fig. 9. No significant differences are found among the three morphological categories (BL, BI, and BD) at any neighbour rank. The median values and overall dispersion of the distributions largely overlap, indicating that the tidal field strength experienced by these systems is similar across morphologies. Only for the first neighbour (i.e. the closest companion) is there a marginal trend for BD galaxies to reach slightly smaller r_{tidal, inst} values, which may indicate that their outer structures are more susceptible to the influence of nearby companions. However, this trend is weak and no significant differences are observed between the subsamples.

Fig. 9. Violin diagrams showing the distribution of the instantaneous tidal radius, rtidal, inst, for the first five nearest neighbours of each galaxy. The three morphological classes (BL, BI, BD) are displayed in different colours.

Overall, these results suggest that neither the large-scale environment, as traced by Σ₅, nor the local tidal field, as traced by r_{tidal, inst}, exerts a dominant influence on the structural differences observed between BL, BI, and BD galaxies. The morphological segregation observed in the stellar mass–size relation therefore appears to be primarily governed by intrinsic factors such as the assembly history.

5. Conclusions

In this work, we have carried out a comprehensive structural analysis of a representative sample of nearby massive BL galaxies from the BEARD survey. Using very deep g- and r-bands imaging, we extracted surface brightness, colour, and stellar mass radial profiles for each system, applying rigorous corrections for PSF effects, scattered light contamination, and inclination. These corrections enabled a homogeneous and accurate measurement of key parameters, most notably the R₁ radius and the central stellar mass density Σ_1, kpc. The inclusion of a matched galaxy sample from the IllustrisTNG50 simulation provides an important theoretical counterpart. We find that the simulations successfully reproduce the observed trends.

Simulated galaxies reveal that BL galaxies systematically occupy the upper envelope of the stellar mass–size relation, exhibiting significantly larger R₁ values at fixed M_* compared to BD systems. This offset is accompanied by a markedly lower central stellar mass density, suggesting a fundamentally different evolutionary path for these galaxies. This trend is in agreement with the findings of Di Cintio et al. (2019), who showed that galaxies with lower central mass densities (i.e. lower surface brightness) develop more extended stellar distributions, resulting in larger effective radii. Although their simulations probed slightly lower stellar masses, at the low end of Milky Way–like systems, the qualitative consistency strengthens the interpretation of our results. Interestingly, while BD galaxies follow a tight sequence in the Σ_1, kpc–M_* plane, BL systems show much greater diversity in central density and disc extent. This trend is only marginally followed by our observed BEARD BL galaxies. Some BEARD galaxies are located in the upper part of the mass–size relation (Fig. 3), as predicted by simulations, while others contribute to the scatter. Despite our limited number of observed galaxies, these differences might be attributed to variations in the central mass surface density of the galaxies.

Further analysis reveals that although we do not find strong global correlations between halo spin and merger-related parameters, BL galaxies tend to reside in halos with slightly higher spin. This may suggest that angular momentum still plays a role in enabling disc growth and preventing central classical bulge formation when merger configurations are favourable. However, the scatter in the mass–size relation does not appear to be primarily driven by the number of mergers experienced by each system, as both BL and BD galaxies show comparable merger counts since z = 5. Instead, the key factor may lie in the specific configuration of those mergers – such as the relative orientation, halo spin, and orbital parameters – which could determine whether the disc structure is preserved or disrupted (Rosas-Guevara et al. 2026). Moreover, our environmental analysis reveals no significant differences in local density or tidal field strength among the BL, BI, and BD populations, suggesting that the environment is unlikely to be the main driver of their structural diversity. However, this result may be partially influenced by the sample selection, particularly the isolation criteria applied to both the observed and simulated galaxies.

The BEARD survey provides a robust observational testbed for evaluating such scenarios and highlights the need to incorporate both dark matter halo properties and cosmological assembly paths when interpreting the structural diversity of disc galaxies in the local Universe.

Acknowledgments

The authors thank the anonymous referee for their thoughtful and constructive comments, which helped to improve the quality and clarity of this manuscript. The authors also sincerely thank Dr. Ignacio Trujillo for his valuable insights and constructive feedback throughout the development of this work. This project is possible thanks to financial support from the Spanish Ministry of Science and Innovation (MICINN) to the CoBEARD project (PID2021-128131NB-I00). JMA acknowledges the support of the Viera y Clavijo Senior programme funded by ACIISI and ULL and the support of the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033) under grant nos. PID2021-128131NB-I00 and CNS2022-135482 and the European Regional Development Fund (ERDF) ‘A way of making Europe’ and the ‘NextGenerationEU/PRTR’. AdLC acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) through RYC2022-035838-I and PID2021-128131NB-I00 (CoBEARD project). Data for this paper has been obtained under the International Time Programme of the CCI (International Scientific Committee of the Observatorios de Canarias of the IAC) with the telescope operated on the island of La Palma by the operators in the Observatorio del Roque de los Muchachos. EAG acknowledges support from the Agencia Espacial de Investigación del Ministerio de Ciencia e Innovación (AEI-MICIN) and the European Social Fund (ESF+) through a FPI grant PRE2020-096361. ADC kindly thanks the Spanish Ministerio de Ciencia, Innovación y Universidades through grant CNS2023-144669, programa Consolidación Investigadora. MCC acknowledges the support of AC3, a project funded by the European Union’s Horizon Europe Research and Innovation programme under grant agreement No 101093129. MCC acknowledges financial support from the Spanish Ministerio de Ciencia, Innovación y Universidades (MCIU) under the grant PID2021-123417OB-I00 and PID2022-138621NB-I00. EMC and AP acknowledge the support by the Italian Ministry for Education University and Research (MUR) grant PRIN 2022 2022383WFT “SUNRISE” (CUP C53D23000850006) and Padua University grants DOR 2022–2024. SZ acknowledges the financial support provided by the Governments of Spain and Aragón through their general budgets and the Fondo de Inversiones de Teruel, the Aragonese Government through the Research Group E16_23R, and the Spanish Ministry of Science and Innovation and the European Union – NextGenerationEU through the Recovery and Resilience Facility project ICTS-MRR-2021-03-CEFCA. The project that gave rise to these results received the support of a fellowship from the “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/PR24/12050015. LC acknowledges support from grants PID2022-139567NB-I00 and PIB2021-127718NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by “ERDF” A way of making Europe. FP acknowledges support from the Horizon Europe research and innovation programme under the Maria Skłodowska-Curie grant “TraNSLate” No 101108180, and from the Agencia Estatal de Investigación del Ministerio de Ciencia e Innovación (MCIN/AEI/10.13039/501100011033) under grant (PID2021-128131NB-I00) and the European Regional Development Fund (ERDF) “A way of making Europe”. VC acknowledges ANID (FONDECYT grant no. 11250723).

References

Appendix A: Observing Log and galaxy sample summary

This appendix provides additional information on the observational campaign and the final galaxy sample used in this work. Table A.1 summarises the observations carried out with the Isaac Newton Telescope (INT) as part of three proposals between 2020 and 2021, including the number of nights allocated and those effectively observed. Table A.2 lists the 22 galaxies that compose the final BEARD deep imaging sample. For each object, we report the J2000 coordinates, limiting surface brightness depths in the g and r bands, total stellar mass, and the R₁ radius.

Appendix B: Subtraction of scattered light from foreground stars

Following the characterisation of the PSF in both g and r bands, presented in Sect. 3.3, we modelled and subtracted the scattered light field produced by foreground stars in each field of the survey. This process is essential, as previous studies showed that the extended halos of bright stars can act as a diffuse illumination screen with effective surface brightness levels of μ_V∼ 28–29 mag arcsec⁻², contaminating the detection and measurement of low surface brightness structures in galaxy outskirts (Slater et al. 2009; Trujillo & Fliri 2016).

Fig. B.1. Example of stellar contamination and subtraction using the extended PSF. Left: colour composite image of NGC2543, created using the SDSS g and r bands and the astscript-color-faint-gray algorithm described in Infante-Sainz & Akhlaghi (2024). Middle: Same field after subtraction of bright stars. Right: Scattered-light map, showing the extended halos of field stars.

To ensure accurate PSF modelling, the first step involved constructing segmentation maps and masks for all detected sources in each field. These segmentation maps are pixel-based classifications that distinguish astronomical sources from the background, identifying which pixels belong to stars, galaxies, or artefacts. This process allowed us to mask contaminating sources and reliably isolate the stars used for building the PSF. This was achieved using NoiseChisel and Segment (Akhlaghi & Ichikawa 2015; Akhlaghi 2019), two tools from the GNU Astronomy Utilities (Gnuastro) suite. These algorithms allow for robust detection of both compact and extended sources and assign unique labels to each detection, enabling their exclusion from the PSF analysis. The resulting masks are essential to prevent contamination from background galaxies or blended sources, particularly in the diffuse light regime.

From the segmentation maps, a preliminary catalogue of detected sources was generated. This was then cross-matched with the Gaia DR3 catalogue (Gaia Collaboration 2016, 2023) using a 5 arcsec matching radius to account for potential positional discrepancies in saturated stars. Following the cross-match, an additional filter was applied to retain only those sources with either parallax or proper motion measurements exceeding three times their associated uncertainties, ensuring a clean stellar sample for PSF construction.

Stars were selected based on their Gaia G-band magnitudes, retaining only those brighter than 16 mag to ensure the removal of all bright stars whose extended scattered light may significantly contaminate the field. Each selected star was modelled individually using the extended PSF and subtracted from the image. The subtraction followed a hierarchical approach: brighter stars were addressed first, as their halos dominate the scattered light distribution. This ordering prevents contamination from bright sources during the modelling of nearby, fainter ones. Further technical details on the star selection and subtraction procedure can be found in Marrero-de la Rosa et al. (in prep.).

Each star was modelled using a two-step process: first, the stellar centroid was fixed using Gaia coordinates; second, the PSF model was scaled in flux to match the observed radial surface brightness profile of the star. To define a reliable normalisation region, we estimated the saturation radius R_sat as a function of stellar magnitude. This radius corresponds to the point beyond which the detector saturation no longer affects the observed profile. It was determined by analysing the first derivative of the radial surface brightness profile and identifying the radius at which the inner saturated region begins to deviate from the expected smooth decline of the PSF. Once R_sat was defined, the flux normalisation was performed over an annular region spanning 1.5 R_sat to 4 R_sat, ensuring that the selected portion of the profile is both unsaturated and has a sufficient S/N for accurate scaling.

A scaling factor was computed as the 2σ-clipped mean ratio between the observed and model PSF within this annulus. Once scaled, the PSF model was subtracted from the image. The final scattered light map was constructed by summing the contributions of all modelled stars. The surface brightness levels of the resulting maps were found to be consistent with typical values reported in the literature for diffuse stellar halos caused by scattered light, with effective levels of μ_g, r ∼ 28–29 mag arcsec⁻² (e.g. Slater et al. 2009; Trujillo & Fliri 2016). An illustrative example of this procedure, showing the original image, the star-subtracted field, and the corresponding scattered-light map, is presented in Fig. B.1. This procedure was applied uniformly across the full BEARD survey. To implement these procedures, we have developed dedicated Python tools, described in Marrero-de la Rosa et al. (in prep.), which will be publicly released through GitHub repositories. The first one, LISAN (Layered Intensity Spread and Analysis for Night-sky structures), is dedicated to the construction and calibration of extended PSFs over wide-field imaging datasets. The second, MAHDI (Mitigation Algorithm for Halo and Diffuse Illumination), performs scattered light subtraction by modelling and removing the halo of each star using the previously characterised PSF. Both codes will be publicly available via the following GitHubhttps://github.com/CarlosMDLR/LISAN

https://github.com/CarlosMDLR/MAHDI repositories.

Appendix C: Correction of self-scattered light via wavelet-based deconvolution

In addition to the contamination from foreground stars, the extended wings of the PSF can scatter light from the galaxy itself into its outskirts. This effect can artificially flatten the surface brightness profiles and bias measurements of outer disc structures and stellar mass distributions.

To correct for this internal scattering, we applied a deconvolution algorithmhttps://github.com/aasensio/Wavelet_deconvolution based on the isotropic undecimated wavelet transform (IUWT) (Starck & Murtagh 1994; Starck et al. 2010; Carrillo et al. 2012), as introduced in Golini et al. (2025) and that will be explained in detail in Marrero-de la Rosa et al. (in prep.). This method leverages the multiscale decomposition of the image to separate smooth, extended structures (associated with PSF scattering) from compact and coherent galaxy features. The algorithm minimises a loss function that includes terms for the fidelity of the observed convolved image, smoothness of the reconstruction, and wavelet-based regularisation to suppress noise amplification.

The output of the deconvolution yields a clean, PSF-corrected image in which light redistributed by the PSF has been reconcentrated towards its physical origin, effectively restoring the intrinsic structure of the galaxy. Residual analysis shows a consistent pattern: central surface brightness increases while the outer profiles steepen, matching theoretical expectations for PSF broadening effects. This correction is essential for an accurate estimation of parameters such as the central mass surface density within 1 kpc (Σ_1, kpc) and the R₁ radius, and was systematically applied to all galaxies in the sample.

To evaluate the impact of the point spread function (PSF) on our surface brightness measurements, we applied a wavelet-based deconvolution algorithm to the g- and r-band images of a representative galaxy in our sample, NGC 7606. The resulting profiles, shown in Fig. C.1, illustrate the expected redistribution of light due to PSF broadening: after deconvolution, the central regions of the galaxy become noticeably brighter, while the outer parts exhibit a steeper decline in surface brightness. This effect highlights how the PSF tends to spread central light towards the outskirts, artificially flattening observed profiles if left uncorrected. The test validates the importance of PSF treatment, particularly for the analysis of low surface brightness features.

Fig. C.1. Effect of wavelet-based PSF deconvolution on the surface brightness profiles of NGC 7606 in the g- and r-bands. The solid curve shows the profile prior to deconvolution while the dashed curve shows the profile after deconvolution.

Appendix D: Masking the data and subtracting the local background

In the analysis of low surface brightness images, the construction of accurate masks is of paramount importance, as even faint, unmasked sources can introduce significant biases in photometric measurements. In this work, a conservative masking strategy was adopted by combining the outputs of multiple algorithms: NoiseChisel and Segment (Akhlaghi & Ichikawa 2015; Akhlaghi 2019) along with Max-Tree Objects (MTO; Teeninga et al. 2016; Hashem Faezi et al. 2024). The final masks were obtained by merging the individual masks produced by these methods, ensuring a robust identification of contaminant sources.

With regard to background estimation and following the methodology described in Golini et al. (2025), we accounted for possible local background variations that may not be fully captured by the global sky subtraction applied during the initial data reduction. Prior to the application of the wavelet-based deconvolution, elliptical isophotes centred on each galaxy were constructed, sharing the same position angle (PA) and axis ratio. These structural parameters were derived using the photutils package (Jedrzejewski 1987; Bradley et al. 2019), where elliptical fits were performed iteratively until a satisfactory solution was reached. The resulting isophotes were visually inspected and manually adjusted in cases where the automated procedure failed to converge or yielded poor fits.

To robustly estimate the local background, we made use of wedges, which we defined as narrow, sector-shaped regions extending radially outwards from the galaxy centre. These wedges trace the galaxy light along specific angular directions and are particularly useful to study the outskirts while minimising contamination. Wedge-shaped profiles were extracted in directions carefully selected to avoid bright stars or background galaxies, allowing us to analyse local background fluctuations in well-defined sectors of the image. In each case, the sky level was estimated by computing the median pixel value within a region deemed free of contaminating flux.

Appendix E: Additional analysis of the mass–size relation

In this appendix we present complementary figures supporting the analysis of the scatter in the mass–size relation. Fig. E.1 shows the linear fits to the log R₁–log M_* relation for BEARD, BL, BI, and BD galaxies. The fitted slopes and intercepts for each population are listed in Table 1, highlighting the systematic differences between the morphological categories.

In addition, Fig. E.2 illustrates the statistical comparison between BEARD and IllustrisTNG50 BL galaxies. In this test, 10⁴ random realisations of BL galaxies were drawn in the stellar mass range 10.0 < log(M_*/M_⊙) < 10.5, each time matching the number of BEARD galaxies in this interval. The resulting distribution of scatters demonstrates that the observed scatter in BEARD is fully compatible with the expectations from simulations, although Kolmogorov–Smirnov tests indicate that the two samples are statistically distinguishable in detail. 0

Fig. E.1. Linear fits to the log R1–log M* relation for BEARD, BL, BI, and BD galaxies. Dashed lines show the best fits for each population; see Table 1 for fitting parameters.

Fig. E.2. Distribution of scatter in log R1 for 104 random realisations of IllustrisTNG50 BL galaxies in the stellar mass range 10.0 < log(M*/M⊙) < 10.5, each time matching the number of BEARD galaxies in this interval. The vertical solid line shows the BEARD scatter, while the dashed line marks the BL mean.

Appendix F: Merger parameters and their relation to halo spin and mass

To further explore the interplay between merger history and dark matter halo properties, we analyse the correlations between three key merger parameters–total number of mergers since z = 5, mean gas fraction, and accreted mass fraction–and three fundamental halo quantities: the spin parameter (λ), the concentration ($c_{200,\mathrm{c}}^{\mathrm{DM}}$) and the halo mass ($M_{200}^{\mathrm{DM}}$). The results are summarised in Fig. F.1.

Fig. F.1. Relations between merger and dark matter properties for the simulated galaxies in our study. The figure is organised as a matrix of scatter plots, where each row corresponds to a merger property (from left to right: mean gas fraction, total mergers since z = 5 and accreted mass fraction), and each column corresponds to a dark matter quantity (from top to bottom: spin parameter λ, halo mass $M_{200}^{\mathrm{DM}}$, and halo concentration $c_{200,\mathrm{c}}^{\mathrm{DM}}$). BL, BI and BD are shown as green circles, blue isodensity contours, and orange diamonds, respectively. On the top and right panels, histograms of the corresponding quantities are shown for BL and BD galaxies, where vertical lines represent the mean of each distribution.

Focusing first on the spin parameter, we find no strong global correlation with any of the merger-related variables. Nevertheless, the previously observed morphological separation persists: bulgeless galaxies (BL) tend to exhibit slightly higher spin values than their bulge-dominated (BD) counterparts, regardless of the merger parameter considered. This subtle separation reinforces the idea that high-spin halos are more likely to host galaxies with reduced central mass concentrations, potentially due to their merger histories being less conducive to bulge growth.

When examining the merger parameters as a function of halo mass, a noticeable positive correlation emerges between the total number of mergers and $M_{200}^{\mathrm{DM}}$. More massive halos tend to have undergone a larger number of mergers since z = 5, consistent with expectations from hierarchical structure formation. In contrast, the remaining panels–mean gas fraction and accreted mass fraction versus halo mass, and also the ones about the concentration–do not show significant trends or morphological distinctions.

Appendix G: Subsample of bulgeless galaxies below the mass–size trend

In this appendix, we provide additional figures related to the analysis of a specific subsample of bulgeless galaxies from the IllustrisTNG50 simulation. As discussed in the main text, a notable fraction of bulgeless systems appears to deviate from the typical location of this morphological class in the mass–size plane, falling below the average 1 σ trend traced by simulated bulgeless galaxies. To investigate their properties in more detail, we isolate a subsample of such galaxies–referred to as below-trend–which lie systematically below the mean R₁ values of bulgeless galaxies in bins of stellar mass. Fig. G.1 shows the location of this below-trend population in the mass–size diagram, highlighting their offset relative to the main bulgeless sequence.

Fig. G.1. Stellar mass–size relation using R1. The BEARD observational sample is shown in light-green squares and the TCK20 reference sample in grey dots. The population of simulated bulgeless galaxies from the IllustrisTNG50 sample is shown in green (BL), while the subsample lying below the mean mass–size trend is highlighted in black squares (BL-below-trend).

All Tables

All Figures

	Fig. 1. Normalised radial profile of the PSFs for all fields observed in the g and r bands for the sample galaxies. The PSF total flux is scaled to unity. Vertical shaded regions indicate structural zones defined in the galaxy profile analysis (Inner, Subintermediate, Intermediate, Outer and Outer Ext) (see Sect. 3.3 for details). The y-axis represents normalised intensity in arbitrary units.
In the text

Fig. 2. Left panel: colour image of NGC 3486. A 1 arcmin scale bar is shown, together with the area used to extract the profiles, highlighted as a wedge. The R1 radius is shown as a cyan dashed-lined ellipse. Middle panel: surface brightness profiles in the g (dashed blue) and r (solid red) bands, along with the corresponding stellar mass surface density profile (black points). The vertical and horizontal dashed grey lines marks the R1 radius and 1 M⊙ pc−2, respectively. Right panel: Colour profile (g − r)o corrected from extinction.

In the text

Fig. 3. Stellar mass–size relation defined by the R1 radius. The main panel shows the global distribution of galaxies from different datasets: light-green squares show the BEARD sample of observed galaxies, grey points correspond to the observational sample from Trujillo et al. 2020 (TCK20), and cyan circles represent simulated galaxies from the IllustrisTNG50 run of the IllustrisTNG project (Cardona-Barrero et al. 2025). The right panels separate the IllustrisTNG50 sample into three categories based on their dynamically decomposed bulge-to-disc ratio: BL galaxies in green, BI galaxies in blue, and BD galaxies in orange. To compute the coloured lines for each morphological class, the stellar mass range was divided into 11 equally spaced bins; in each bin, the average R1 value was calculated, and the shaded regions indicate the corresponding standard deviation. For comparison, the TCK20 (grey dots) and BEARD (light-green squares) samples are also shown in each panel.

In the text

Fig. 4. Relation between the specific central mass density ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and the global properties of galaxies. Left panel: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ as a function of stellar mass. Right panel: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ versus R1. In both panels, simulated galaxies from IllustrisTNG50 are colour-coded according to their bulge prominence, with BD galaxies (orange diamonds) occupying the upper regions of the relations, while BL systems (green circles) lie systematically below. BI galaxies (blue crosses) span a wide range in both ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$ and size, bridging the two extremes. Observational data from the BEARD sample are shown as light-green squares, consistently occupying the region where BL galaxies are located.

In the text

Fig. 5. Relations between stellar and dark matter properties for the simulated galaxies from IllustrisTNG50. The figure is organised as a matrix, where each row corresponds to a stellar quantity (from top to bottom: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$, M* and R1) and each column corresponds to a dark matter halo property (from left to right: λ, $M_{200}^{\mathrm{DM}}$, and $c_{200,\mathrm{c}}^{\mathrm{DM}}$). BL, BI, and BD galaxies are shown as green circles, blue isodensity contours, and orange diamonds, respectively. On the top and right panels, histograms of the corresponding quantities are shown for BL and BD galaxies, where vertical lines represent the mean of each distribution.

In the text

Fig. 6. Relations between stellar and merger properties for the simulated galaxies from IllustrisTNG50. The figure is organised as a matrix, where each row corresponds to a stellar quantity (from top to bottom: ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$, M* and R1) and each column corresponds to a merger property (from left to right: mean gas fraction, total mergers since z = 5 and accreted mass fraction). BL, BI, and BD galaxies are shown as green circles, blue isodensity contours, and orange diamonds, respectively. On the top and right panels, histograms of the corresponding quantities are shown for BL and BD galaxies, where vertical lines represent the mean of each distribution.

In the text

	Fig. 7. Normalised distributions of three structural parameters for the simulated BL galaxies from IllustrisTNG50 (grey) and simulated BL galaxies that lie significantly below the average mass–size trend (hatched cyan). Left: R1 radius,. Middle: central stellar mass density within 1 kpc, ${\mathrm{\Sigma }}_{1,\mathrm{kpc}}^{\mathrm{spec}}$. Right: mean gas fraction since z = 5.
In the text

	Fig. 8. Distribution of the projected environmental density parameter, Σ5, for the simulated BL, BI, and BD galaxies compared with the BEARD-Phot and BEARD-Spec samples.
In the text

	Fig. 9. Violin diagrams showing the distribution of the instantaneous tidal radius, rtidal, inst, for the first five nearest neighbours of each galaxy. The three morphological classes (BL, BI, BD) are displayed in different colours.
In the text

	Fig. B.1. Example of stellar contamination and subtraction using the extended PSF. Left: colour composite image of NGC2543, created using the SDSS g and r bands and the astscript-color-faint-gray algorithm described in Infante-Sainz & Akhlaghi (2024). Middle: Same field after subtraction of bright stars. Right: Scattered-light map, showing the extended halos of field stars.
In the text

	Fig. C.1. Effect of wavelet-based PSF deconvolution on the surface brightness profiles of NGC 7606 in the g- and r-bands. The solid curve shows the profile prior to deconvolution while the dashed curve shows the profile after deconvolution.
In the text

	Fig. E.1. Linear fits to the log R1–log M* relation for BEARD, BL, BI, and BD galaxies. Dashed lines show the best fits for each population; see Table 1 for fitting parameters.
In the text

	Fig. E.2. Distribution of scatter in log R1 for 104 random realisations of IllustrisTNG50 BL galaxies in the stellar mass range 10.0 < log(M*/M⊙) < 10.5, each time matching the number of BEARD galaxies in this interval. The vertical solid line shows the BEARD scatter, while the dashed line marks the BL mean.
In the text

Fig. F.1. Relations between merger and dark matter properties for the simulated galaxies in our study. The figure is organised as a matrix of scatter plots, where each row corresponds to a merger property (from left to right: mean gas fraction, total mergers since z = 5 and accreted mass fraction), and each column corresponds to a dark matter quantity (from top to bottom: spin parameter λ, halo mass $M_{200}^{\mathrm{DM}}$, and halo concentration $c_{200,\mathrm{c}}^{\mathrm{DM}}$). BL, BI and BD are shown as green circles, blue isodensity contours, and orange diamonds, respectively. On the top and right panels, histograms of the corresponding quantities are shown for BL and BD galaxies, where vertical lines represent the mean of each distribution.

In the text

	Fig. G.1. Stellar mass–size relation using R1. The BEARD observational sample is shown in light-green squares and the TCK20 reference sample in grey dots. The population of simulated bulgeless galaxies from the IllustrisTNG50 sample is shown in green (BL), while the subsample lying below the mean mass–size trend is highlighted in black squares (BL-below-trend).
In the text