© The Authors 2026

1. Introduction

Significant advancements in the observational study of Active Galactic Nuclei (AGNs) have provided growing evidence that massive black holes (MBHs) with masses above 10⁵ M_⊙ form naturally in galaxies and power AGN activity through gas accretion at their centres (Schmidt 1963; Merloni & Heinz 2008; Ueda et al. 2014; Hopkins et al. 2007; Aird et al. 2015). Moreover, AGN demographic studies and dynamical analyses of stars and gas in the central regions of nearby galaxies have shown that a large fraction of massive galaxies in the local Universe host MBHs in their nuclei (Genzel & Townes 1987; Kormendy 1988; Dressler & Richstone 1988; Kormendy & Richstone 1992; Genzel et al. 1994; Peterson et al. 2004; Vestergaard & Peterson 2006).

Over the past few decades, numerous studies have advanced our understanding of MBH formation and growth, revealing a series of tight correlations between MBH mass and key properties of their host galaxies. One of the first relations was reported by Kormendy & Richstone (1995), who found a connection between MBH mass and bulge luminosity. Shortly thereafter, Magorrian et al. (1998) established that MBH mass correlates with bulge stellar mass, an association later refined by Marconi & Hunt (2003) and Häring & Rix (2004). An even tighter correlation was identified by Ferrarese & Merritt (2000), who showed that MBH mass scales strongly with the stellar velocity dispersion (see also Gebhardt et al. 2000). Additional relations have been proposed, such as those linking MBH mass to the Sérsic index of the galaxy (Graham et al. 2001) and to the total stellar mass of the host (Reines & Volonteri 2015; Erwin & Gadotti 2012; Capuzzo-Dolcetta & Tosta e Melo 2017; Harikane et al. 2023). All of these scaling relations carry significant physical implications, as they point to a close co-evolution between galaxy assembly and MBH growth. On the one hand, they suggest that star formation (SF), bulge build-up, and gas accretion onto the MBH are interconnected processes (see, e.g. Silk & Rees 1998; Croton 2006; Wild et al. 2010; Ishibashi & Fabian 2012; DeGraf et al. 2015; Carraro et al. 2020; Habouzit et al. 2019, 2021). On the other hand, they imply that feedback from actively accreting MBHs can influence the baryonic evolution of the host galaxy, regulating SF and shaping its stellar mass build-up (see, e.g. Di Matteo et al. 2005; Croton 2006; Croton et al. 2006; DeGraf et al. 2017; Martín-Navarro et al. 2018).

The estimation of MBH masses used in the scaling relations has been commonly carried out using dynamical modelling of the motions of gas and stars in the galactic centre (Kormendy & Ho 2013), reverberation mapping techniques (Peterson et al. 2004; Lu et al. 2025), or X-ray variability methods (e.g. McHardy et al. 2006). In contrast, the galaxy properties such as stellar mass are typically inferred from mass-to-light ratios or from spectral energy distribution (SED) fitting (Conroy 2013). Although each of these approaches carries its own systematic and statistical uncertainties (see, for instance, Vestergaard & Peterson 2006; Shen et al. 2013), the resulting measurements consistently show that the established scaling relations exhibit an intrinsic scatter, with a population of overmassive MBHs (lying above the median relation) and undermassive MBHs (lying below it). Regarding the former population, overmassive MBHs have been observed across different cosmological epochs. Recent results at z > 5 from the James Webb Space Telescope (JWST) have unveiled a population of MBHs that are tentatively up to two orders of magnitude more massive than their local counterparts (see e.g. Maiolino et al. 2024; Harikane et al. 2023; Kocevski et al. 2023; Pacucci et al. 2023; Matthee et al. 2024; Jones et al. 2025). These discoveries have challenged the standard paradigm of MBH formation and evolution, suggesting that MBHs either originate from heavy seeds (∼10⁴ M_⊙) or grow from light seeds (∼100 M_⊙) through episodic super-Eddington accretion (Inayoshi 2025; Bonoli et al. 2025). Similar overmassive MBHs are also observed at lower redshifts (z ≤ 3), implying that such systems are not necessarily outliers in the context of the galaxy–MBH correlation but a natural outcome (Ferré-Mateu et al. 2021; Mezcua et al. 2023, 2024). Concerning the population of undermassive MBHs, several studies have reported their presence in the local Universe (Hu 2008; Kormendy & Ho 2013; Savorgnan et al. 2016), often associated with pseudobulge-dominated galaxies. Similar cases are observed in the AGN population (Graham & Li 2009), suggesting that active, undermassive MBHs may be in a phase of rapid growth, potentially evolving towards the median scaling relation.

All of the recent observational discoveries motivate the identification and theoretical characterisation of MBH outliers in the scaling relation, as they offer a powerful probe of the physical processes governing MBH growth and galaxy evolution. For instance, overmassive MBHs that lie above the median M_BH − M_* relation, might indicate systems where black hole growth has proceeded more rapidly than the assembly of the host galaxy stellar mass (see e.g. Weller et al. 2023). This accelerated growth may result from prolonged or highly efficient accretion episodes. Conversely, the population of undermassive MBHs, found below the median relation, may correspond to galaxies in which stellar mass assembly has overpassed the MBH growth rate (either through efficient secular or external processes), systems where MBH growth is still in an early phase, or cases where the MBH assembly mechanism is hindered by different processes that limit gas inflow (e.g. Volonteri & Natarajan 2009; Anglés-Alcázar et al. 2013). Alternatively, external mechanisms such as ejections and displacements from the center of the galaxy after gravitational recoils have also been proposed as a mechanism to produce undermassive MBHs (Blecha et al. 2011; Izquierdo-Villalba et al. 2020; Dong-Páez et al. 2025). In this regard, several observational studies have recently reported candidates of recoiled MBHs (e.g. van Dokkum et al. 2023; Barrows et al. 2025; Islam et al. 2026), including what may represent the first confirmed case (van Dokkum et al. 2025). However, the interpretation of some systems remains debated (Sánchez Almeida 2023), and the overall frequency and impact of recoil events are still uncertain.

Within the framework of hierarchical structure formation, all the processes described above are expected to operate concurrently, collectively shaping the full distribution of MBHs around the median M_BH − M_* relation and naturally producing both overmassive and undermassive MBH systems. Understanding the emergence of these two populations requires, therefore, tracing how the relative importance of MBH accretion efficiency, galaxy assembly, feedback, and dynamical processes evolves with redshift. In this work, we address this by exploring the emergence of overmassive and undermassive MBHs in the M_BH − M_* relation at different galaxy masses and cosmological epochs. To this end, we employ the state-of-the-art semi-analytical model (SAM) L-GalaxiesBH, a flexible and well-tested galaxy formation framework that can be run on top of the merger trees from the Millennium suite of simulations (Bonoli et al. 2025). We stress that our goal is to focus on the physical conditions that give rise to the over- and under-massive systems, irrespective of their observability. We are aware that many low-mass MBHs may fall below current detection thresholds, limiting their identification and preventing a full characterisation of the true scaling relation. As a result, it would be challenging to determine the over- or under-massive nature of MBHs. The present work examines exclusively the physical origin of the overmassive and undermassive MBH regimes, while a dedicated assessment of their detectability and observational biases is deferred to future studies. The paper is organised as follows: In Section 2 we present the merger trees from the Millennium simulation suite, the main features of the L-GalaxiesBH semi-analytical model, and the new prescription implemented for gradual baryonic stripping. In Section 3 we present the definition of overmassive and undermassive MBHs and characterise their formation scenarios. Finally, in Section 4 we summarise the key findings. A Λ Cold Dark Matter (ΛCDM) cosmology with parameters Ω_m = 0.315, Ω_Λ = 0.685, Ω_b = 0.045, σ₈ = 0.9 and H₀ = 67.3 km s⁻¹ Mpc⁻¹ is adopted throughout the paper (Planck Collaboration XVI 2014).

2. The galaxy formation model: L-GalaxiesBH SAM

In this section, we summarise the main physics implemented in the L-GalaxiesBH (Bonoli et al. 2025). L-GalaxiesBH builds on the public L-Galaxies SAM (Henriques et al. 2015), which follows the cosmological assembly of galaxies by solving analytical equations along the dark matter (DM) halo merger trees from N-body simulations. Updates to L-Galaxies (Izquierdo-Villalba et al. 2020, 2022; Spinoso et al. 2023) have led to L-GalaxiesBH, which allows a detailed treatment of the formation and evolution of single and binary MBHs.

2.1. Dark matter merger trees

L-GalaxiesBH is designed to work on the merger trees derived from DM only N-body simulations. Here, we use the merger trees extracted from the Millennium (MS, Springel 2005) and Millennium-II (MSII, Boylan-Kolchin et al. 2009) simulations. The MS follows the evolution of 2160³ DM particles of 8.6 × 10⁸ M_⊙/h in a 500 Mpc/h box from z = 127 to the present. MSII tracks the same number of particles but in a 100 Mpc/h box with 125 times higher mass resolution (6.885 × 10⁶ M_⊙/h). Snapshots were stored at 63 (MS) and 68 (MSII) epochs, with typical time intervals of ∼300 Myr (L-GalaxiesBH will make an internal internal interpolation of ∼5 − 20 Myr). Halos and substructures were identified using the friend-of-friend and SUBFIND algorithms and arranged into merger trees with L-HALOTREE (Springel et al. 2001). Both simulations were rescaled to Planck Collaboration XVI (2014) cosmology following Angulo & White (2010) methodology.

Taking into account the above, the MS enables the exploration of large cosmological volumes but at the expense of a limited mass resolution, which limits the study of low-mass galaxies. In contrast, the merger trees of MSII provide sufficient resolution to trace the formation and evolution of small galaxies. However, its smaller volume hinders its capacity to statistically sample rare high-redshift and massive low-redshift systems. To overcome these limitations and exploit the complementary strengths of both simulations, we do not present independent results from MS and MSII. Instead, we combine their merger trees into the ones called MS+Grafting, following the grafting procedure described in Bonoli et al. (2025). This approach proceeds in two steps. First, the merger tree branches of MS are extended below the resolution limit by incorporating higher-resolution information from MSII, thereby recovering merger events that would otherwise be unresolved. Second, when the MS haloes (above and below resolution) are first resolved, they are initialised with representative galaxies drawn from the population of objects generated by L-GalaxiesBH with MSII trees. This ensures that the systems are modelled as already evolved rather than as pristine baryonic reservoirs. Therefore, the MS+Grafting enables us to construct merger trees with the same high resolution as the MSII but within the larger volume (and statistics) of the MS. Therefore, unless otherwise stated, all the results presented here will correspond to the ones obtained with the MS+Grafting merger trees.

2.2. Galaxy and MBH formation and evolution

This section describes the galaxy and MBH formation framework of L-GalaxiesBH (Henriques et al. 2015) and introduces a new prescription for gradual galaxy stripping.

2.2.1. Formation and assembly of galaxies: Gradual galaxy stripping

The L-GalaxiesBH model builds on the seminal works of White & Rees (1978) and White & Frenk (1991). In this framework, a collapsing DM halo captures a fraction of cosmic baryons, forming a hot gaseous halo. This one can cool down via radiative processes and settle at the halo centre, forming a rotationally supported disc. Continuous cold gas accretion fuels SF, forming a stellar disc component and triggering supernova (SN) feedback, which regulates the galaxy assembly by heating and expelling cold gas. In parallel, gas accretion of the central MBH from the hot atmosphere moderates further gas cooling via the injection of energy. Galactic bulges can form via secular processes, such as disc instabilities (DI) triggered by intense and sustained episodes of SF, or through external events caused by galaxy mergers. Specifically, mergers follow the coalescence of parent DM halos on timescales estimated by Binney & Tremaine (2008) and are divided between major and minor: while major mergers (m_R > 0.2) produce spheroidal remnants, minor mergers (m_R < 0.2) preserve the primary disc while growing the bulge via satellite accretion.

L-GalaxiesBH also accounts for environmental effects, including ram-pressure stripping and tidal disruption, as key mechanisms shaping the galaxy population (Henriques et al. 2015, 2020). Regarding tidal disruption, in the fiducial L-GalaxiesBH implementation, satellite galaxies are assumed to be instantaneously disrupted once tidal forces from the host DM halo become strong enough to unbind them. In this work, we modify this treatment by allowing for the gradual stripping of the satellite stellar and cold-gas components before any complete disruption. Specifically, we assume that the material of the satellite can remain self-bound only out to a radius R_t, defined as (King 1962; Taylor & Babul 2001):

$$\begin{array}{c}\hfill R_t=\left(\frac{G{M}_{\mathrm{gal}}}{{\omega }^2-{\mathrm{\Phi }}_{\mathrm{DM}}^{″}(r)}\right),\end{array}$$(1)

where r is the instantaneous radial position of the satellite within its host DM halo, ω is its orbital angular velocity, and $\varphi {″}_{\mathrm{DM}}(r)$ is the second radial derivative of the gravitational potential generated by the DM halo, which is assumed to follow a Navarro–Frenk–White (NFW) density profile (Navarro et al. 1996). Here, M_gal denotes the total baryonic mass of the satellite galaxy (stars plus cold gas). Assuming this R_t definition, any material located beyond this radius is assumed to become unbound from the satellite and subsequently removed from the galaxy. As in the full disruption case, the stripped cold gas is transferred to the hot-gas reservoir of the host halo, while the stripped stellar component is deposited into the diffuse stellar halo. For simplicity, after the stripping process, the remaining cold gas and stellar mass of the galaxy are determined just as the mass enclosed within R_t, assuming an exponential profile for discs and a Hernquist profile for the bulge (Freeman 1970; Hernquist 1990):

$$\begin{array}{cc}& M_i^{\mathrm{disc},\mathrm{N}}=M_i^{\mathrm{disc}}[1-(1+\frac{R_t}{R_i})e^{-R_t/R_i}],\hfill \end{array}$$(2)

$$\begin{array}{cc}& M_{\ast }^{\mathrm{bulge},\mathrm{N}}=M_{\ast }^{\mathrm{bulge}}\left[\frac{R_t^2}{{(R_t+\frac{R_b}{1+\sqrt{2}})}^2}\right],\hfill \end{array}$$(3)

where $M_i^{\mathrm{disc}}=\{M_{\mathrm{g}}^{\mathrm{disc}},M_{\ast }^{\mathrm{disc}}\}$, and $M_{\ast }^{\mathrm{bulge}}$ are the initial masses of the cold gas disc, stellar disc, and stellar bulge, while $M_i^{\mathrm{disc},\mathrm{N}}=M_{\mathrm{g}}^{\mathrm{disc},\mathrm{N}},M_{\ast }^{\mathrm{disc},\mathrm{N}}$, and $M_{\ast }^{\mathrm{bulge},\mathrm{N}}$ are the corresponding new masses after stripping. R_i = {R_g, R_*} and R_b are the effective radii of the cold gas, stellar disc, and stellar bulge, respectively. The new effective radii of these components after the stripping (R_g^N, R_*^N and R_b^N) are derived under the assumption that their central densities are preserved:

$$\begin{array}{cc}& R_i^N={\left(\frac{M_i^{\mathrm{disc},\mathrm{N}}}{M_i^{\mathrm{disc}}}\right)}^{1/2}{R}_i,i=\{\mathrm{g},\ast \},\hfill \end{array}$$(4)

$$\begin{array}{cc}& R_b^N={\left(\frac{M_{\ast }^{\mathrm{bulge},\mathrm{N}}}{M_{\ast }^{\mathrm{bulge}}}\right)}^{1/3}\left(\frac{R_b}{1+\sqrt{2}}\right).\hfill \end{array}$$(5)

2.2.2. The formation of the first MBHs

The MBH seeding framework in L-GalaxiesBH is described in Spinoso et al. (in prep.) and Bonoli et al. (2025). In summary, it represents an extension of the model proposed by Spinoso et al. (2023), which introduces multiple formation channels for MBH seeds. The updated scheme allows for the coexisting emergence of four seed types: Population III remnants (PopIII, ∼100 M_⊙), runaway stellar mergers (RSM, ∼10^3 − 4 M_⊙), direct-collapse black holes (DCBHs, ∼10⁵ M_⊙), and merger-induced DCBHs (miDCBH, ∼8 × 10⁴ M_⊙). To provide a physical basis for DCBH and RSM seeds, Spinoso et al. (2023) connected their formation with the intergalactic metal enrichment and radiative feedback in the form of Lyman-Werner radiation, both driven by star formation and supernova activity. Given the coarse resolution of the Millennium simulations, Spinoso et al. (in prep.) and Bonoli et al. (2025) implemented a sub-grid probabilistic model for PopIII remnant formation, assigning seeds based on halo mass and redshift, with parameters calibrated to reproduce the MBH occupation statistics of Spinoso et al. (2023). The initial seed mass of PopIII remnants is inferred from the stellar mass formed during the halo unresolved evolution, assuming a Larson initial mass function (IMF) and a single, short-lived starburst. Between the time of seed formation and the resolution of its host halo, the seed is allowed to grow at the Eddington rate, modelling the unresolved accretion phase. Finally, similar to RSM and DCBH seeds, the PopIII seeding is terminated once the mean IGM metallicity exceeds a critical value.

2.2.3. Gas accretion onto MBHs: The emergence of the M_BH − M_* relation

In L-GalaxiesBH, newly formed MBHs are initially assigned a random spin, which is then evolved self-consistently through gas accretion and mergers (Izquierdo-Villalba et al. 2020). The growth of MBHs is primarily driven by the accretion of cold gas funneled towards the galactic centre after galaxy mergers or DIs. This inflowing gas accumulates in a central reservoir (M_res) around the MBH and is accreted in two stages (Izquierdo-Villalba et al. 2020, 2024). In the first one, accretion can be Eddington-limited or super-Eddington, depending on the local conditions. Super-Eddington episodes occur when the gas reservoir is large (ℛ = M_res/M_BH > 2 × 10⁴) and the inflow rate is high (Ṁ_inflow > 10 M_⊙/yr). During such events, the Eddington ratio (f_Edd = L_bol/L_Edd) follows the form:

$$\begin{array}{c}\hfill f_{\mathrm{Edd}}=B(\chi )[\frac{0.985}{{\dot{M}}_{\mathrm{Edd}}/\dot{M}+C(\chi )}+\frac{0.015}{{\dot{M}}_{\mathrm{Edd}}/\dot{M}+D(\chi )}],\end{array}$$(6)

where Ṁ and Ṁ_Edd are the MBH and Eddington accretion rates, and L_bol and L_Edd the corresponding luminosities. The functions B(χ), C(χ), and D(χ) are adopted from Madau et al. (2014). When these criteria are not met, accretion proceeds at the standard Eddington limit (f_Edd = 1). The second stage begins once the MBH consumes part of the reservoir (Bonoli et al. 2025), with accretion declining as:

$$\begin{array}{c}\hfill f_{\mathrm{Edd}}={[1+{((t-t_0)/t_Q)}^{1/2}]}^{-2/\beta },\end{array}$$(7)

where t_Q = t_d ξ^β/(βln10) sets the characteristic decline timescale, with t_d = 1.26 × 10⁸ yr, β = 0.4, and ξ = 0.3, consistent with the self-regulated MBH growth models of Hopkins & Hernquist (2009). All free parameters are calibrated to reproduce local M_BH − M_* scaling relations and the stochastic gravitational wave background measured by PTAs (EPTA Collaboration et al. 2023; Agazie et al. 2023; Reardon et al. 2023). In this framework, the M_BH − M_* relation naturally emerges from the coupled evolution of galaxies and their central MBHs, as accretion is intrinsically linked to processes regulating stellar mass (mergers and DI), producing correlated growth across cosmic time.

2.2.4. Formation and evolution of massive black hole binaries

The L-GalaxiesBH model does not assume that, after a galaxy merger, the satellite MBH instantaneously reaches the centre of the remnant galaxy, forms a bound binary with the central MBH, and coalesces. Instead, it follows the dynamical evolution of the MBH from kiloparsec separations down to parsec scales, where massive black hole binaries (MBHBs) form (Izquierdo-Villalba et al. 2022). In brief, after a merger, the satellite MBH is typically deposited several kiloparsecs from the nucleus, where it undergoes a dynamical friction phase (Binney & Tremaine 2008) that gradually brings it to the center (hereafter referred to as pairing MBHs in the pairing phase), eventually forming a gravitationally bound MBHB with the central MBH. Once bound, the MBHB orbital separation and eccentricity evolve via coupled equations that account for the surrounding environment (Izquierdo-Villalba et al. 2022). In gas-rich systems, the inspiralling and eventual coalescence are driven by circumbinary gas interactions and gravitational wave (GW) emission (Dotti et al. 2015; Peters & Mathews 1963), whereas in stellar-dominated systems, three-body interactions with stars (following a Sérsic profile) dominate alongside GWs (Quinlan & Hernquist 1997; Sesana & Khan 2015; Peters & Mathews 1963). In cases of repeated mergers, a third MBH may reach the nucleus before the pre-existing binary coalesces. In such scenarios, L-GalaxiesBH models the resulting triplet dynamics following Bonetti et al. (2018). The model also accounts for MBH mass growth during both the pairing and hardening phases. Specifically, the satellite MBH can both accumulate gas for accretion during the merger and consume it as it migrates toward the galactic nucleus (see Izquierdo-Villalba et al. 2022, 2023a,b). Once a bound MBHB forms, accretion proceeds according to the preferential accretion approach supported by hydrodynamical simulations (e.g., D’Orazio et al. 2013; Farris et al. 2014; D’Orazio & Duffell 2021).

2.2.5. Gravitational wave recoils: Ejections from the galactic centre

When an MBHB merges, L-GalaxiesBH computes the remnant spin following Barausse & Rezzolla (2009), which depends on the mass ratio, spin magnitudes, and relative orientations of the MBH spins and the orbital angular momentum. Since L-GalaxiesBH tracks only spin magnitude, orientations remain unconstrained. To overcome this limitation, the SAM distinguishes two merger types (see also Barausse 2012; Volonteri et al. 2013): wet and dry. In wet mergers (where the binary mass is smaller than the available gas reservoir), the MBH spins align within approximately 10° (Dotti et al. 2010), while in dry mergers they are assumed to be randomly oriented. The magnitude and orientation of the progenitor spins is also important to determine the gravitational recoil received by the remnant MBH due to the anisotropic emission of GWs. L-GalaxiesBH computes this kick velocity using the prescription of Lousto et al. (2012). If the recoil velocity exceeds the escape speed of the host galaxy, the MBH is expelled and deposited within the DM halo, where L-GalaxiesBH subsequently tracks its orbit as presented in Izquierdo-Villalba et al. (2020).

3. Results

This section investigates how different physical processes shape the population of overmassive and undermassive MBHs in the M_BH − M_* plane within the L-GalaxiesBH model. We stress that our focus is on the physical conditions that give rise to these systems, irrespective of their observability. We are aware that many low-mass MBHs may fall below current detection thresholds, limiting their identification and preventing a full characterisation of the true scaling relation. As a result, it would be challenging to determine the real over- or under-massive nature of MBHs.

To investigate the populations of overmassive and undermassive MBHs at all redshifts, we define five galaxy samples, each corresponding to a specific percentile range of the global M_BH − M_* distribution:

Median sample:

• (i)

  1σ sample (hereafter, 1σ): Galaxies lying within the central 68% of the global M_BH − M_* distribution (16^th − 84^th percentiles). Overmassive samples:

• (ii)

  +2σ sample (hereafter, +2σ): Galaxies in the upper intermediate tail (84^th − 97^th percentiles).

• (iii)

  +3σ sample (hereafter, +3σ): Galaxies in the extreme upper tail (97^th − 100^th percentiles). Undermassive samples:

• (iv)

  −2σ sample (hereafter, −2σ): Galaxies in the lower intermediate tail (3^rd − 16^th percentiles).

• (v)

  −3σ sample (hereafter, −3σ): Galaxies in the extreme lower tail (0^th − 3^rd percentiles).

The distinction described above enables a direct comparison between galaxies hosting over- and under- massive MBHs, providing a standardized framework to investigate the physical processes driving deviations from the average M_BH − M_* relation¹. The main results are presented in the following sections, with a summary of their impact on the formation of over- and under- massive populations given in Table 1. In this context, Fig. 1 illustrates the distribution of the five samples in the M_BH − M_* plane at z = 0, 1, 3, and 5². Overall, the model shows good agreement with the z ∼ 0 observations of Capuzzo-Dolcetta & Tosta e Melo (2017), Erwin & Gadotti (2012) and Reines & Volonteri (2015) as well as with the measurements of ultra-compact dwarf galaxies presented by Ferré-Mateu et al. (2021) and the MBHs population unveiled via TDE events by Ramsden et al. (2026). At z = 1 − 3, it also remains consistent with the sample of Suh et al. (2020) and is capable of generating a population of very massive MBHs in low-mass galaxies, as detected by Mezcua et al. (2023, 2024). At z ∼ 5, the model also shows reasonable agreement with the massive galaxy population unveiled by JWST observations. It is noteworthy that, compared to the model version presented in Bonoli et al. (2025), the inclusion of gradual stripping leads to a depletion of z ∼ 5 MBHs with M_BH > 10⁶ M_⊙ in galaxies of M_* ∼ 10⁸ M_⊙. This hinders the production of a large population of overmassive MBH in low-mass galaxies, as recently observed by JWST. We have verified that this difficulty arises because gas stripping systematically reduces the masses of satellite galaxies at the time of merger, relative to a model without stripping. Consequently, the mass ratios in galaxy mergers are smaller, leading to less efficient gas inflows. Therefore, fewer low-mass galaxies experience super-Eddington accretion events. Indeed, Bonoli et al. (2025) showed that such events in L-GalaxiesBH are responsible for generating the JWST-detected MBHs in low-mass galaxies. Our results therefore, indicate that, if current JWST observations are confirmed, the implementation of gradual stripping in our model must be less efficient to allow the formation of overmassive MBHs in high-z low-mass galaxies.

Fig. 1. MBH − M* plane at z = 0, 1, 3, 5. Different colours represent different populations. The results have been compared with the z ∼ 0 sample of Erwin & Gadotti (2012) (squares) Reines & Volonteri (2015) (diamonds) and Capuzzo-Dolcetta & Tosta e Melo (2017) (circles), Ferré-Mateu et al. (2021) (pentagons) and Ramsden et al. (2026) (crosses). The z = 1 and z = 3 results presents the observational sample of Suh et al. (2020) (triangles) and the overmassive MBHs of Mezcua et al. (2023, 2024), respectively. The z = 5 results are compared with Maiolino et al. (2024) (circles), Harikane et al. (2023) (squares), Ding et al. (2023) (stars) and Lupi et al. (2024) (pentagons). At this redshift, the black (grey) points correspond to z ≥ 5 (z ≤ 5). The horizontal line highlights ≤104 M⊙, a mass range where electromagnetic information about the central MBH of the galaxy is expected to be undetectable by current observational facilities.

To further characterise the behaviour of our samples, Fig. 2 presents the redshift evolution of the M_BH/M_* ratio restricted to galaxies with M_* > 10⁸ M_⊙. The 1σ sample displays a pronounced rise in M_BH/M_* from z ∼ 6 (M_BH/M_* ∼ 5 × 10⁻⁶) to z ∼ 2 (M_BH/M_* ∼ 10⁻³), followed by a decline toward the present epoch, reaching M_BH/M_* ∼ 2 × 10⁻⁴ at z = 0. The +2σ and +3σ samples exhibit similar evolutionary trends, although their M_BH/M_* ratios are systematically higher by factors of ∼5 − 7 relative to the 1σ population. In contrast, the −2σ and −3σ samples show little to no evolution: the −3σ population maintains a nearly constant ratio of ∼7 × 10⁻⁶, while the −2σ sample remains flat at ∼10⁻⁶ up to z ∼ 2, followed by an increase of about one order of magnitude at lower redshifts. To assess whether these trends persist across different stellar masses, Fig. 2 also shows the results for galaxies with M_* > 10¹⁰ M_⊙. The general behaviour of the 1σ, −2σ, and −3σ samples remains broadly consistent with that seen in the lower-mass sample, although their M_BH/M_* ratios at z > 3 are systematically lower. However, the +2σ and +3σ populations show a markedly stronger evolution: while their ratios at z < 2 are comparable to those of less massive galaxies, they increase dramatically at higher redshifts (z > 3), by up to four orders of magnitude. For instance, at z ∼ 5 the +3σ (+2σ) sample reaches M_BH/M_* ∼ 0.1 (∼0.01), consistent with the overmassive MBH population recently uncovered by JWST observations (see e.g. Matthee et al. 2024; Harikane et al. 2023).

Fig. 2. Redshift evolution of the median MBH/M* ratio for the five different samples: +3σ (red), +2σ (coral), 1σ (orange), −2σ (blue), and −3σ (purple). Solid lines with circles represent the results for galaxies with M* > 108 M⊙ while dashed lines with squares represent the same but for galaxies with M* > 1010 M⊙. The shaded areas correspond to the 16th − 84th percentiles.

The results shown in Fig. 2 indicate that the scaling relation predicted by L-GalaxiesBH exhibits a clear evolution with redshift. At high redshift, the relation also displays a pronounced dependence on galaxy stellar mass. These trends contrast with several previous theoretical studies, which report little to no evolution in the M_BH − M_* relation across cosmic time (see Volonteri et al. 2016; Huang et al. 2018; Habouzit et al. 2019; Marshall et al. 2020). Such discrepancies are likely driven by differences in numerical resolution and in the modelling of key physical processes, including MBH growth, MBH seeding, AGN/SN feedback and SF prescriptions (see, e.g. Habouzit et al. 2021). A brief comparison of the M_BH/M_* ratio between L-GalaxiesBH and five hydrodynamical simulations is provided in Appendix A. However, a detailed investigation of the origin of differences is beyond the scope of this work.

In the next sections, we investigate the physical processes that drive deviations from the median M_BH − M_* relation, focusing on how overmassive and undermassive MBH populations emerge within the L-GalaxiesBH framework.

3.1. The role of the environment

As described in Section 2, following a halo–halo merger, the galaxy hosted by the smaller halo sinks toward the center of the new host dark matter halo. During this process, it undergoes gradual stripping of stars and cold gas, which can displace the system from its original location in the M_BH − M_* relation. This effect is illustrated in the first column of Fig. 3, which shows the evolutionary track of a galaxy and its central MBH extracted from L-GalaxiesBH. As shown, once the galaxy becomes satellite and experiences gradual stripping (indicated by the vertical dashed line), the MBH–galaxy system shifts along the M_BH − M_* relation: the MBH mass remains unchanged, while the stellar mass decreases by roughly a factor of eight.

Fig. 3. Examples of evolutionary pathways for overmassive and undermassive MBHs. The upper panels show the galaxy merger trees, the middle panels illustrate the assembly of stellar mass and MBH mass, and the lower panels present the evolution of the ratio MBH/M*. The first column corresponds to an overmassive MBH whose growth is driven by the gradual stripping of the host galaxy. The second column shows an undermassive MBH produced by a GW recoil event. The third column depicts an undermassive MBH resulting from a quiescent merger history, while the fourth column presents an overmassive MBH formed through an active merger history combined with a super-Eddington accretion episode.

To evaluate the broader role of gradual stripping in generating over- or under-massive MBHs, Fig. 4 presents the fraction of galaxies in each sample affected by this process. For the −3σ, −2σ, 1σ, and +2σ populations, less than 10% of galaxies experience stellar stripping, independently of redshift or stellar mass. However, the +3σ sample exhibits a markedly different behaviour, independently of mass: less than 5% of galaxies are affected at z > 3, but this fraction rises steadily toward lower redshifts, from ∼10% at z ∼ 3 to ∼50% at z ∼ 0. A mild dependence on stellar mass is also evident, with smaller galaxies being slightly more likely to undergo stripping³. Beyond the occurrence of this environmental process, the magnitude of stellar mass loss further differentiates the impact of stripping across samples. This is shown in the lower panel of Fig. 4, which reports the median fraction of stellar mass lost by each galaxy since the halo–halo merger. In the −3σ, −2σ, 1σ, and +2σ samples, stripping has a negligible effect, with mass losses rarely exceeding 10%, independent of stellar mass or redshift. Conversely, the +3σ sample experiences significantly larger losses, which increase with decreasing redshift and decreasing stellar mass. For galaxies with 10⁸ < M_* < 10⁹ M_⊙, the stripped fraction grows from ∼20% at z ∼ 3 to ∼90% by z = 0, with a comparable trend for galaxies in the range 10⁹ < M_* < 10¹⁰ M_⊙. Galaxies with M_* > 10¹⁰ M_⊙ are less affected, with mass losses reaching at most ∼50% by z = 0. We stress that these fractions might remain sensitive to the underlying modelling assumptions about stripping, particularly the predicted galaxy sizes in L-GalaxiesBH and the adopted galaxy mass profiles, which together determine how much mass lies beyond R_t (see Eq. (1)). Regarding the former, in L-GalaxiesBH galaxy scale lengths are computed self-consistently from angular momentum and energy conservation during galaxy assembly, and are in good agreement with both high-redshift JWST constraints (Herrero-Carrión et al. 2026) and low-redshift observations (Izquierdo-Villalba et al. 2019), supporting the adopted structural properties. By contrast, the assumed mass profiles (see Section 2.2.1) can more directly regulate the stripped fraction in +3σ systems: more extended stellar or gas distributions increase the mass beyond R_t, enhancing stripping, while more compact profiles have the opposite effect. Despite these dependencies, the qualitative result about gradual stripping leading to the build-up of extreme systems at low redshift is expected to remain robust.

Fig. 4. Upper panel: Redshift evolution of the fraction of galaxies that are undergoing a gradual stripping (fSatellite). Lower panel: Redshift evolution of the median fraction of stellar mass lost by galaxies due to gradual stripping (fMass Stripped). In all panels, there are 5 different samples: +3σ (red), +2σ (coral), 1σ (orange), −2σ (blue) and −3σ (purple). The left column corresponds to galaxies with 108 < M* < 109 M⊙, the middle column to those with 109 < M* < 1010 M⊙, and the right column to galaxies with M* > 1010 M⊙.

Taken together, the results presented in this section indicate that stellar stripping contributes to the scatter in the M_BH − M_* relation by producing uniquely (extreme) overmassive MBHs via the reduction of their host galaxies stellar mass (> 20%), while leaving the MBH mass unaffected. These results are consistent with the findings of Weller et al. (2023), who showed that in the TNG50 hydrodynamical simulation more than 95% of satellite galaxies in overdense environments (where tidal stripping is likely to occur) lie above the global M_BH − M_* relation. On the observational side, Ferré-Mateu et al. (2021) found that ultracompact dwarfs and compact ellipticals (known to be low-mass systems because of their in situ formation or the stripping of massive galaxies) systematically host overmassive MBHs. Despite these results, our work indicates that environmental process are significant only for ∼40% of the +3σ population at z < 3 (∼10% of the whole overmassive sample), and it cannot account for the overmassive MBHs observed in the +2σ population at the same redshift, nor in the +2σ and +3σ samples at z > 3. These findings demonstrate that gradual stellar stripping alone is insufficient to explain the full population of overmassive MBHs, implying that additional physical mechanisms must contribute throughout cosmic time.

3.2. The role of the ejections after gravitational wave recoils

Several theoretical studies have demonstrated that gravitational recoil following MBH coalescences can eject MBHs from their host galaxies, leaving a measurable imprint on galaxy scaling relations (Blecha et al. 2011; Gerosa & Sesana 2015; Izquierdo-Villalba et al. 2020; Dong-Páez et al. 2025). Building on these studies, and to investigate the impact of gravitational recoils and subsequent MBH ejections on the M_BH − M_* relation across cosmic time, the upper panel of Fig. 5 shows the fraction of galaxies that have experienced at least one central MBH ejection over their lifetime. The results reveal a clear dependence on both stellar mass and sample type. For galaxies with 10⁸ < M_* < 10⁹ M_⊙, fewer than 20% of systems in the 1σ, +2σ, and +3σ samples undergo an ejection event. The −2σ and −3σ samples display moderately higher ejection fractions, increasing from ∼20% at z ≥ 3 to ∼40% by z ∼ 0. The relatively low ejection fractions in these M_* < 10⁹ M_⊙ galaxies are primarily driven by their reduced MBH occupation fraction, which L-GalaxiesBH estimates to be in the range 50%−80% depending on redshift (see Figure 7 of Bonoli et al. 2025). This lowers the probability that galaxy mergers result in MBH coalescence and, therefore, the frequency of producing recoil events. In the intermediate stellar mass range, 10⁹ < M_* < 10¹⁰ M_⊙, the ejection fractions in the 1σ, +2σ, and +3σ samples show a moderate raise from 10% at z > 2 to 30% at lower redshifts. In contrast, the −2σ and −3σ samples exhibit stronger redshift evolution. In the −2σ sample, the fraction of galaxies experiencing at least one MBH ejection rises from ∼40% at z ∼ 5 to ∼80% at z < 1. For the −3σ sample, the increase is more rapid, with the ejection fraction already saturating at ∼80% by z ∼ 4. Finally, for galaxies with M_* > 10¹⁰ M_⊙, the overall behaviour is qualitatively similar to that observed at lower stellar masses. However, in the −2σ and −3σ samples, the ejection fractions show little to no redshift dependence, having already reached saturation levels of 90%−95% by z ∼ 6.

Fig. 5. First panel: Redshift evolution of the fraction of galaxies that underwent an ejection of their central MBH due to GW recoils (fejection). Second panel: Median time spent in the dynamical-friction phase (tDF) by MBHs that re-filled an empty galactic nucleus due to a GW ejection. Third panel: Redshift evolution of the median redshift at which the MBH was ejected from the center of the galaxy (zejection). The colour coding and mass dependence are the same as in Fig. 4.

The results above indicate that gravitational recoil and the resulting MBH ejection are common processes primarily affecting undermassive MBHs. This has important implications for the evolution of galaxy–MBH systems in the M_BH − M_* plane, as the absence of a central MBH after an ejection creates an opportunity for any pairing MBH orbiting within the galaxy (coming from a past galaxy merger) to migrate to the galaxy nucleus and establish itself as the new central MBH⁴. This process can reposition the galaxy-MBH system either upward or downward in the M_BH − M_* plane, depending on the mass of the incoming MBH. An illustrative example is shown in the second panel of Fig. 3, which shows a galaxy-MBH system in which, following an MBH merger at z ∼ 0.15, the remnant MBH is ejected. After ∼1 Gyr of the central MBH ejection, the galaxy nucleus is filled by another MBH that was deposited in the galaxy at z ∼ 1 (i.e ∼6 Gyr ago) and is much lighter (∼10⁴ M_⊙) than expected for the host galaxy stellar mass (∼10¹⁰ M_⊙). This drastically reduces the galaxy M_BH/M_* ratio from ∼10⁻⁴ at z ∼ 0.25 to ∼10⁻⁷ at z ∼ 0.

The undermassive outcome presented in Fig. 3 is not unique to the illustrative example, but instead reflects a general behaviour of pairing MBHs in our model. Although L-GalaxiesBH allows pairing MBHs acquired through galaxy mergers to grow, their off-center location causes them to rapidly exhaust their available gas reservoirs, stalls their growth and prevents any sustained co-evolution with their new host galaxies. As a result, when these MBHs eventually settle in a devoid galaxy nucleus, they have masses close to those they had when they were initially acquired during the galaxy merger, typically lower than what would be expected for the new host galaxy. Therefore, the duration of the pairing phase plays a key role in setting the degree of desynchronized growth between the host galaxy and the MBH that eventually repopulates its empty nucleus, and thus determines the final location of the galaxy–MBH system in the M_BH − M_* plane (see Appendix B for the measurement of the galaxy-MBH growth desynchronization via duty cycle analysis). Motivated by this, it is important to examine how externally acquired MBHs in different samples migrate toward the nucleus in galaxies that temporarily lack a central MBH due to gravitational ejections. This is presented in the second row of Fig. 5, which depicts the time these MBHs spend in the pairing phase within the galaxy before settling at the center. The model predicts that the pairing phases of the −2σ and −3σ samples varies from 0.3 Gyr at z ∼ 5 to 3 Gyr at z ∼ 1, independently on the galaxy mass. Conversely, the pairing phase in the 1σ, +2σ, and +3σ samples is typically much shorter, with characteristic durations of ∼0.2 − 0.4 Gyr. Interestingly, the duration of the pairing phase tends to be slightly longer in more massive systems, as these galaxies are more extended and MHBs can be deposited at larger galactocentric distances than in lower-mass galaxies. In addition to the duration of the pairing phase, the timing of the gravitational ejection itself constitutes a second critical factor in determining whether a galaxy–MBH system becomes undermassive, overmassive, or converges toward the median relation: even when the pairing phase is short, an ejection occurring close to the galaxy observed redshift may leave insufficient time for the newly established central MBH to grow in synchrony with its host and converge toward the median M_BH − M_* relation. To quantify this effect, the third row of Fig. 5 shows the redshift of the last MBH ejection for the different samples. In the 1σ, +2σ, and +3σ samples, ejections predominantly occur well before the galaxies reach their observed redshift, with little dependence on stellar mass. Combined with the short pairing phases, this explains why galaxies in the 1σ, +2σ, and +3σ samples can experience gravitational recoils yet remain consistent with the average or overmassive populations: the newly acquired MBHs sink efficiently to the nucleus and have enough time to co-evolve with their hosts before observation. Conversely, Fig. 5 indicates that in the −2σ and −3σ samples, ejections typically occur much closer to the galaxy observed redshift. In these systems, the long pairing phases together with the limited time available for co-evolution inhibit significant convergence toward the median M_BH − M_* relation, leaving the systems as undermassive.

To link the dynamical processes discussed above with galaxy population statistics, Fig. 6 compares the predicted M_BH/M_* ratios in models with and without MBH ejections. For galaxies with M_* > 10⁸ M_⊙, the absence of gravitational recoils systematically increases the M_BH/M_* ratios in the 1σ, −2σ, and −3σ samples, reflecting a reduced fraction of strongly undermassive systems. This effect is particularly pronounced at z < 3, where the no-recoil model produces M_BH/M_* ratios in the −2σ and −3σ samples up to 1 dex higher than in the model that includes gravitational ejections. In more massive galaxies (M_* > 10¹⁰ M_⊙), the discrepancy between the two models grows even larger, highlighting that recoil-driven ejections have the strongest impact at the high-mass end of the M_BH − M_* relation.

Fig. 6. Redshift evolution of the median MBH/M* ratio for a model in which no GW kicks and ejections are included. The upper (lower) panel corresponds to galaxies with M* > 108 M⊙ (M* > 1010 M⊙) and solid lines correspond to the model without (with) ejections after gravitational recoils. The colour coding and mass dependence are the same as in the upper panel of Fig. 2.

The analysis presented in this section enables us to conclude that gravitational recoils play a key role in producing exclusively undermassive MBHs in galaxies with M_* > 10¹⁰ M_⊙ at all redshifts, as well as in ∼40% of systems with M_* < 10⁹ M_⊙ at z < 3. The physical mechanism behind the formation of this undermassive population involves ejection events of MBHs after gravitational recoils and the subsequent repopulation of the galactic nucleus by pairing MBHs left over from early galaxy mergers. These newly settled MBHs are usually light as they spend extended periods outside the central regions and only settle into the nucleus close to the observed epoch, leaving them insufficient time to grow in step with their host galaxies. Nonetheless, our results reveal that recoils alone cannot fully explain all undermassive MBHs: even in a no-recoil model, a population of undermassive MBHs persists, particularly among M_* < 10⁹ M_⊙ galaxies. In the next section, we analyse how different evolutionary paths contribute to the development of undermassive and the emergence of overmassive MBHs over cosmic time.

3.3. From a deficit to an excess: Gas accretion process after mergers and secular evolution

Previous sections have demonstrated that stellar mass stripping can account for the emergence of a subset of overmassive MBHs at z < 3, while gravitational recoil events explain the presence of undermassive MBHs in galaxies with M_* > 10⁹ M_⊙ across all redshifts. However, these processes alone cannot fully reproduce the overmassive MBH population observed at high redshift, nor the undermassive population in galaxies with M_* < 10⁹ M_⊙. To investigate the origin of these remaining populations, we now turn to the growth histories of MBHs, which are shaped by both DIs (internal secular processes) and galaxy mergers (external events)⁵.

– Galaxy mergers: Focusing on the role of external events, Fig. 7 presents the fraction of MBHs in each sample that have experienced at least one growth episode triggered by major or minor mergers. As shown, major mergers are subdominant across all redshifts and stellar mass bins, with fewer than 5% of MBHs in any sample undergoing such events. In contrast, minor mergers play a much more significant role. In galaxies with M_* > 10⁹ M_⊙, over 90% of MBHs in the +3σ, +2σ, 1σ, and −2σ samples, and over 60% in the −3σ sample, have experienced growth episodes triggered by minor mergers, independent of redshift. Despite such a high occurrence rate, the impact of minor mergers on MBHs in these galaxies varies across samples. This is illustrated in Fig. 8 which shows that MBHs in the 1σ, +2σ, and +3σ samples typically undergo a median of 10 − 50 minor merger driven growth episodes, a rate 5 − 10 times higher than that of the −2σ and −3σ populations. This higher frequency of minor merger driven growth is accompanied by systematically larger baryonic mass ratios in the +3σ, +2σ, and 1σ samples, which increase with redshift and saturate at 0.02 − 0.03 at z = 0, roughly a factor of three larger than in the −3σ and −2σ samples. In the case of MBHs hosted in M_* < 10⁹ M_⊙ galaxies, we observe trends broadly consistent with those found in MBHs hosted in more massive systems, but only for the +3σ, +2σ, and 1σ samples. Conversely, the −2σ and −3σ populations exhibit significantly fewer minor merger driven growth episodes, with only ∼60% and ∼20% of cases. Concerning the frequency of MBH growth via minor mergers, the +3σ and +2σ samples undergo an average of 5 − 10 events, whereas the −2σ and −3σ samples experience typically only 1 − 3 events. To illustrate the impact of such limited merger activity on MBH growth, the third row of Fig. 3 shows the merger history of a galaxy (M_* ∼ 10⁹ M_⊙) hosting an undermassive MBH at z = 5. As seen, the galaxy undergoes few mergers, which restricts the growth of its central MBH, leaving it close to its initial seed mass (∼100 M_⊙). Taking together these results imply that minor mergers contribute proportionally more to the growth of the central MBH than to the build-up of the host galaxy stellar mass, thereby having a relevant role in the production of over- and under-massive MBHs.

Fig. 7. Upper panel: Redshift evolution of the fraction of MBHs in a given stellar mass bin that underwent a growth episode via minor merger (fmin merg). Middle panel: Redshift evolution of the fraction of MBHs in a given stellar mass bin that experienced a growth phase after a major merger (fmaj,  merg). Lower panel: Redshift evolution of the fraction of MBHs in a given stellar mass bin that had a growth phase after a disc instability (fDI). The colour coding and mass dependence are the same as in Fig. 4. A major merger is defined as a galaxy interaction whose baryonic (gas+stars) mass ratio is larger than 0.2. Otherwise, the interaction is tagged as a minor merger.

Fig. 8. Upper panel: Redshift evolution of the median number of minor mergers that triggered the growth of MBHs hosted in different galaxies (Nminor,  merg). Lower panel: Redshift evolution of the median maximum baryonic mass ratio of minor mergers ($m_{R,\mathrm{minor}}^{\max }$) that triggered the MBH growth in different galaxies. The colour coding and mass dependence are the same as in Fig. 4.

– Disc instabilities: Beyond the influence of galaxy interactions, secular evolution within galactic discs can trigger DI events, directing gas toward the galactic center and promoting central MBH growth. To quantify the impact of these events in our population of MBHs, the lower row of Fig. 7 shows the fraction of MBHs in different samples that underwent a growth phase driven by DIs. As shown, MBHs residing in M_* < 10⁹ M_⊙ galaxies receive only a minor contribution from DIs, with fewer than 20% of MBHs in this mass range growing via such events across all five samples, independently of redshift. For MBHs in galaxies with M_* > 10⁹ M_⊙, the 1σ, −2σ, and −3σ samples continue to exhibit the same low contribution from DIs. However, DI events play a substantially larger role in the +2σ and +3σ samples and show a clear dependence on galaxy mass: the fraction of MBHs growing via DIs rises from 30%−40% in galaxies with 10⁹ < M_* < 10¹⁰ M_⊙ to nearly 90% in galaxies with M_* > 10¹⁰ M_⊙, independent of redshift. These results emphasise that, alongside minor mergers, DIs constitute a major channel for MBH growth in the +2σ and +3σ samples for galaxies with M_* > 10⁹ M_⊙ across all redshifts.

The trends highlighted above reveal that the specific growth histories of MBHs largely determine the emergence of under- and over-massive systems, accounting for populations that gravitational recoils and stellar stripping alone cannot explain. Undermassive MBHs, primarily comprising the −2σ and −3σ samples, have experienced exceptionally quiescent evolutionary pathways, characterised by negligible contributions from DIs and only a few galaxy merger events. Such quiet growth histories are most prevalent in low-mass galaxies (M_* < 10⁹ M_⊙), where the presence of undermassive MBHs cannot be fully explained by gravitational recoil events alone. In contrast, overmassive MBHs, represented by the +2σ and +3σ samples, are associated with highly active growth histories, driven by both minor mergers and DI events. This vigorous evolutionary activity becomes increasingly common at higher redshifts, particularly for DI-driven episodes, where stellar stripping alone cannot account for the formation of overmassive MBHs. In fact, in these rapidly evolving environments, efficient gas inflows can drive accretion rates approaching or exceeding the Eddington limit, with super-Eddington phases likely playing a critical role in enabling rapid MBH growth at early cosmic times (see, e.g., spectroscopically identified MBHs in JWST observations, Maiolino et al. 2024; Harikane et al. 2023). This behaviour is illustrated in the fourth column of Fig. 3, which shows the merger history of a galaxy hosting an overmassive MBH at z = 4. As shown, thanks to the rich merger history, the central MBH undergoes a super-Eddington accretion episode, pushing the M_BH/M_* ratio from 10⁻⁶ to ∼2 × 10⁻³ in just a few Myr.

To quantify the role of super-Eddington accretion in shaping overmassive MBHs, Fig. 9 presents the fraction of MBHs experiencing at least one super-Eddington phase. As shown, overmassive MBHs in galaxies with M_* < 10⁹ M_⊙ rarely undergo such episodes, indicating that the emergence of these systems in low-mass hosts is driven by sustained and frequent growth episodes rather than by short-lived phases of extremely rapid accretion. By contrast, overmassive MBHs in more massive galaxies (M_* > 10⁹ M_⊙) exhibit a significantly higher incidence of super-Eddington episodes. Specifically, these events are concentrated at z > 3, with more than 50% of MBHs in the +2σ and +3σ samples experiencing at least one super-Eddington phase. In comparison, the overall MBH population (1σ sample) shows systematically lower fractions, rarely exceeding 30−40%. Examining the redshift distribution of super-Eddington episodes (second and third panels of Fig. 9), we can see that MBHs in the +2σ and +3σ samples begin their rapid growth earlier than those in the 1σ sample. Specifically, for overmassive MBHs observed at z > 3, both the first and last super-Eddington episodes occur at z > 5, whereas for overmassive MBHs observed at z < 3, these episodes typically occur around z ∼ 4. Furthermore, the lower panels of Fig. 9 show that the maximum nuclear inflow rates associated with mergers or DI events triggering super-Eddington phases are consistently higher in the +2σ and +3σ samples than in the 1σ sample. This underscores the importance of super-Eddington accretion in driving the growth of overmassive MBHs. Independent of galaxy mass, high-redshift episodes are more intense, with inflow rates reaching up to twice those observed at lower redshift. Finally, it is worth noticing that the super-Eddington model of L-GalaxiesBH is controlled by different parameters calibrated to reproduce the high-z quasar bolometric luminosity functions (Izquierdo-Villalba et al. 2024). Variations that make the criteria for super-Eddington accretion more stringent would reduce the number density of overmassive black holes (see Appendix 1 of Izquierdo-Villalba et al. 2024). However, these changes will not alter the qualitative conclusion that high-z overmassive MBHs assemble primarily through active growth episodes in dense, gas-rich environments, where large inflows of fresh gas fuel rapid MBH accretion.

Fig. 9. First panel: Redshift evolution of the fraction of galaxies whose central MBH underwent at least one (solid line) or five (pale lines) super-Eddington phases (fSupEdd). Second panel: Redshift evolution of the epoch at which MBHs within a given stellar mass bin experienced their most recent super-Eddington accretion phase ($z_{\mathrm{SupEdd}}^{\mathrm{Last}}$). Third panel: Redshift evolution of the epoch at which MBHs within a given stellar mass bin experienced their first super-Eddington accretion phase ($z_{\mathrm{SupEdd}}^{\mathrm{First}}$). Fourth panel: Median value of the largest mass inflow (Ṁinflow) that triggered the super-Eddington phase in the MBH. The colour coding and mass dependence are the same as in Fig. 4.

4. Conclusions

In this paper, we explore the physical mechanisms driving the formation of overmassive and undermassive MBHs in the M_BH − M_* scaling relation across cosmic time and stellar mass. To this end, we have used the state-of-the-art L-GalaxiesBH SAM, an extension of L-Galaxies code that self-consistently follows the formation and evolution of single and binary MBHs on top of cosmological merger trees of the Millennium suite of simulations. The model incorporates improved prescriptions for MBH seeding, gas accretion (including super-Eddington phases), stellar and gas stripping, and the full dynamical evolution of MBHs and MBH binaries. The main findings of this work can be summarised as follows:

• The model predicts a redshift evolution in the M_BH/M_* ratio. For the global population (M_* > 10⁸ M_⊙), the value rises from ∼10⁻⁵ at z > 4 to ∼10⁻³ at z ∼ 2, and then declines to ∼10⁻⁴ by z = 0. The undermassive population shows little evolution, with ratios remaining between 10⁻⁶ − 10⁻⁷. The overmassive population follows the same trend as the global population, but with ratios up to 1 dex higher. Notably, for this population, there is a mass dependence: high-z galaxies with M_* > 10¹⁰ M_⊙ exhibit larger ratios (∼0.1) than their z < 2 counterparts.

• Gradual stripping of the galaxy stellar mass contributes significantly to the emergence of overmassive MBHs at z < 3, independently of the galaxy stellar mass. However, this mechanism can explain only ∼10% of the low-redshift overmassive population and fails to account for the presence of overmassive MBHs at high redshift.

• A rich history of mergers and secular evolution drives frequent and sustained episodes of MBH growth, naturally leading to the formation of overmassive MBHs across all redshifts and stellar masses. In M_* < 10⁹ M_⊙ galaxies, overmassive MBHs are primarily associated with highly active merger histories, experiencing on average ∼10 interactions (∼5 times more than the average population). In M_* > 10⁹ M_⊙ galaxies, overmassive MBHs arise from the combination of both vigorous merger activity and frequent disc instabilities driven by secular evolution. Unlike their low-mass counterparts, these massive systems often create the conditions necessary for super-Eddington accretion events, allowing MBHs to grow rapidly and reach overmassive status.

• MBH ejections driven by gravitational recoils play a crucial role in the formation of undermassive MBHs in galaxies with M_* > 10⁹ M_⊙. In particular, the absence of a central MBH after a gravitational ejection creates an opportunity for any pairing MBH orbiting within the galaxy (coming from a past galaxy merger) to migrate to the galaxy nucleus and establish itself as the new central MBH. The long sinking timescales of these MBHs (often exceeding 0.3 Gyr at z ∼ 5 and 3 Gyr at z ∼ 1), combined with the fact that gravitational ejections tend to occur shortly before the galaxy is observed, prevent the newly established MBH from growing in step with its host. As a result, these systems remain persistently undermassive relative to the expected M_BH − M_* relation.

• The rise of undermassive MBHs in galaxies with M_* < 10⁹ M_⊙ cannot be attributed to gravitational-driven MBH ejections, but instead reflects the intrinsically quiescent merger and secular evolutionary histories of their host galaxies. These systems rarely experience disc instability events, and fewer than 60% undergo galaxy interactions capable of feeding the central MBH. Moreover, both the frequency and the intensity of such mergers are significantly reduced, as these galaxies typically experience only a single interaction, compared to an average of three mergers in the broader population of MBHs hosted in M_* > 10⁹ M_⊙ galaxies.

Our results indicate that the over- and under-massive systems in the M_BH − M_* relation arises from the combined, and often competing, effects of: (i) stellar stripping, which shifts systems upward in the M_BH − M_* plane (ii) gravitational recoils and MBH ejections, which shift systems downward and (iii) the diversity of MBH fueling histories, where super-Eddington bursts and quiescent growth move MBHs upward and downward, respectively. The interplay of these processes is essential to reproduce both the overmassive MBHs recently observed by JWST and the undermassive MBHs placed in high-mass galaxies. These findings support a picture in which the co-evolution of galaxies and MBHs is fundamentally non-linear and merger-driven, with environmental and dynamical processes playing a key role in shaping the full diversity of MBH populations across cosmic time.

Acknowledgments

DIV thanks Monica Colpi, Alberto Sessana, Alessandro Lupi, Mar Mezcua, Daniele Spinoso and Silvia Bonoli for their valuable feedback. The project that gave rise to these results received the support of a fellowship from “la Caixa” Foundation (ID 100010434). The fellowship code is LCF/BQ/PI25/12100024. DIV also acknowledges the financial support provided under the European Union’s H2020 ERC Consolidator Grant “Binary Massive Black Hole Astrophysics” (B Massive, Grant Agreement: 818691) and the European Union Advanced Grant “PINGU” (Grant Agreement: 101142079).

References

Appendix A: Comparison of the M_BH/M_* ratio with Hydrodynamical Simulations

In this appendix we compare the evolution of the M_BH/M_* ratio (for M_* > 10¹⁰ M_⊙ galaxies) predicted by L-GalaxiesBH with that obtained from five cosmological hydrodynamical simulations: TNG-50, TNG-100, TNG-300 (Nelson et al. 2019b,a; Pillepich et al. 2019), SIMBA (Flagship, Davé et al. 2019), and Astrid (Bird et al. 2022). These simulations differ in their galaxy formation prescriptions, numerical resolution, and adopted initial black hole seed masses.

As shown in Fig. A.1, the largest differences between L-GalaxiesBH and the hydrodynamical simulations appear at z > 3, where the latter predict typical M_BH/M_* ratios between a factor of ∼ 1  −  3 dex larger than those obtained with L-GalaxiesBH. This discrepancy is primarily driven by the assumed seed masses. In L-GalaxiesBH, MBHs are predominantly seeded through a light-seed channel associated with PopIII remnants, with typical initial masses of 10  −  100 M_⊙. In contrast, most hydrodynamical simulations adopt seed masses of 10⁴  −  10⁵ M_⊙, which naturally facilitates the formation of systems with higher M_BH/M_* ratios at early times. To enable a more direct comparison and reduce the impact of the seed-mass assumption, we recomputed the M_BH/M_* evolution in L-GalaxiesBH considering only systems hosting MBHs with masses >  10⁶ M_⊙, comparable to the effective minimum black hole mass present in hydrodynamical simulations after seeding. Under this selection, the predicted M_BH/M_* ratios in L-GalaxiesBH shift toward higher values and become more consistent with those obtained from the hydrodynamical simulations. The remaining differences are likely driven by other aspects of the galaxy and MBH assembly models, such as the treatment of supernova and MBH feedback or star formation. At z < 3, L-GalaxiesBH shows good agreement with the predictions of SIMBA and Astrid, while the TNG simulations still yield M_BH/M_* ratios larger by approximately a factor of seven. Nevertheless, regardless of the specific hydrodynamical simulation considered, all models exhibit a qualitatively similar trend to that found in L-GalaxiesBH: the M_BH/M_* ratio increases toward lower redshift, although the amplitude of this evolution varies among simulations.

Fig. A.1. Redshift evolution of the median MBH/M* ratio for galaxies with M∗ > 1010 M⊙ in TNG-50, TNG-100, TNG-300, SIMBA (flagship), and Astrid hydrodynamical simulations. The results are compared with L-GalaxiesBH accounting for all the galaxies hosting a MBH (solid line) and galaxies hosting only MBHs with MBH > 106 M⊙. The shaded areas correspond to the 16th  −  84th percentiles.

Appendix B: The duty cycle of the average, overmassive and unedermassive population

In this appendix, we quantify the degree of desynchronization between the growth of galaxies and their central MBHs in the undermassive population. To this end, Fig. B.1 presents the redshift evolution of the MBH duty cycle for the different samples. We define the duty cycle as the fraction of time since MBH formation during which the MBH is actively accreting (including super-Eddington and Eddington phases). The median population (1σ) exhibits a clear redshift dependence: the duty cycle is relatively high (∼ 30  −  40%) at z > 3, while it declines to ∼ 10  −  20% at lower redshifts. A similar evolutionary trend is found for the +2σ and +3σ samples, although with systematically higher duty cycles at all redshifts, reaching values of ∼ 60  −  70% at z > 3. In contrast, the undermassive samples (−2σ and −3σ) show little to no redshift evolution, maintaining low duty cycles of ∼ 5  −  20% across all epochs. This difference is particularly pronounced at stellar masses M_* < 10¹⁰ M_⊙, where, at fixed stellar mass, the duty cycles of the undermassive systems are lower by a factor of ∼ 4  −  5 compared to the median population.

Fig. B.1. Duty cycle as a function of stellar mass for the different samples. Each row corresponds to a given redshift bin (from top to bottom: z = 0, 1, 3, 5). Points indicate median values, while error bars show the 16th  −  84th percentiles. For clarity, the different samples are slightly offset along the M* axis (by 0.05 dex) to avoid overlap.

Taken together, these trends provide a quantitative measure of the degree of desynchronization between galaxy growth and MBH accretion in the undermassive population. At fixed stellar mass, MBHs in such a sample spend substantially less time in active accretion phases than their counterparts in the median population. This indicates that, although their host galaxies continue to build up stellar mass, MBH growth is significantly delayed relative to the general population. Consequently, these MBHs remain underdeveloped and follow evolutionary pathways that are largely decoupled from the canonical co-evolutionary track.

All Tables

All Figures

Fig. 1. MBH − M* plane at z = 0, 1, 3, 5. Different colours represent different populations. The results have been compared with the z ∼ 0 sample of Erwin & Gadotti (2012) (squares) Reines & Volonteri (2015) (diamonds) and Capuzzo-Dolcetta & Tosta e Melo (2017) (circles), Ferré-Mateu et al. (2021) (pentagons) and Ramsden et al. (2026) (crosses). The z = 1 and z = 3 results presents the observational sample of Suh et al. (2020) (triangles) and the overmassive MBHs of Mezcua et al. (2023, 2024), respectively. The z = 5 results are compared with Maiolino et al. (2024) (circles), Harikane et al. (2023) (squares), Ding et al. (2023) (stars) and Lupi et al. (2024) (pentagons). At this redshift, the black (grey) points correspond to z ≥ 5 (z ≤ 5). The horizontal line highlights ≤104 M⊙, a mass range where electromagnetic information about the central MBH of the galaxy is expected to be undetectable by current observational facilities.

In the text

	Fig. 2. Redshift evolution of the median MBH/M* ratio for the five different samples: +3σ (red), +2σ (coral), 1σ (orange), −2σ (blue), and −3σ (purple). Solid lines with circles represent the results for galaxies with M* > 108 M⊙ while dashed lines with squares represent the same but for galaxies with M* > 1010 M⊙. The shaded areas correspond to the 16th − 84th percentiles.
In the text

Fig. 3. Examples of evolutionary pathways for overmassive and undermassive MBHs. The upper panels show the galaxy merger trees, the middle panels illustrate the assembly of stellar mass and MBH mass, and the lower panels present the evolution of the ratio MBH/M*. The first column corresponds to an overmassive MBH whose growth is driven by the gradual stripping of the host galaxy. The second column shows an undermassive MBH produced by a GW recoil event. The third column depicts an undermassive MBH resulting from a quiescent merger history, while the fourth column presents an overmassive MBH formed through an active merger history combined with a super-Eddington accretion episode.

In the text

Fig. 4. Upper panel: Redshift evolution of the fraction of galaxies that are undergoing a gradual stripping (fSatellite). Lower panel: Redshift evolution of the median fraction of stellar mass lost by galaxies due to gradual stripping (fMass Stripped). In all panels, there are 5 different samples: +3σ (red), +2σ (coral), 1σ (orange), −2σ (blue) and −3σ (purple). The left column corresponds to galaxies with 108 < M* < 109 M⊙, the middle column to those with 109 < M* < 1010 M⊙, and the right column to galaxies with M* > 1010 M⊙.

In the text

Fig. 5. First panel: Redshift evolution of the fraction of galaxies that underwent an ejection of their central MBH due to GW recoils (fejection). Second panel: Median time spent in the dynamical-friction phase (tDF) by MBHs that re-filled an empty galactic nucleus due to a GW ejection. Third panel: Redshift evolution of the median redshift at which the MBH was ejected from the center of the galaxy (zejection). The colour coding and mass dependence are the same as in Fig. 4.

In the text

	Fig. 6. Redshift evolution of the median MBH/M* ratio for a model in which no GW kicks and ejections are included. The upper (lower) panel corresponds to galaxies with M* > 108 M⊙ (M* > 1010 M⊙) and solid lines correspond to the model without (with) ejections after gravitational recoils. The colour coding and mass dependence are the same as in the upper panel of Fig. 2.
In the text

Fig. 7. Upper panel: Redshift evolution of the fraction of MBHs in a given stellar mass bin that underwent a growth episode via minor merger (fmin merg). Middle panel: Redshift evolution of the fraction of MBHs in a given stellar mass bin that experienced a growth phase after a major merger (fmaj,  merg). Lower panel: Redshift evolution of the fraction of MBHs in a given stellar mass bin that had a growth phase after a disc instability (fDI). The colour coding and mass dependence are the same as in Fig. 4. A major merger is defined as a galaxy interaction whose baryonic (gas+stars) mass ratio is larger than 0.2. Otherwise, the interaction is tagged as a minor merger.

In the text

	Fig. 8. Upper panel: Redshift evolution of the median number of minor mergers that triggered the growth of MBHs hosted in different galaxies (Nminor,  merg). Lower panel: Redshift evolution of the median maximum baryonic mass ratio of minor mergers ($m_{R,\mathrm{minor}}^{\max }$) that triggered the MBH growth in different galaxies. The colour coding and mass dependence are the same as in Fig. 4.
In the text

Fig. 9. First panel: Redshift evolution of the fraction of galaxies whose central MBH underwent at least one (solid line) or five (pale lines) super-Eddington phases (fSupEdd). Second panel: Redshift evolution of the epoch at which MBHs within a given stellar mass bin experienced their most recent super-Eddington accretion phase ($z_{\mathrm{SupEdd}}^{\mathrm{Last}}$). Third panel: Redshift evolution of the epoch at which MBHs within a given stellar mass bin experienced their first super-Eddington accretion phase ($z_{\mathrm{SupEdd}}^{\mathrm{First}}$). Fourth panel: Median value of the largest mass inflow (Ṁinflow) that triggered the super-Eddington phase in the MBH. The colour coding and mass dependence are the same as in Fig. 4.

In the text

	Fig. A.1. Redshift evolution of the median MBH/M* ratio for galaxies with M∗ > 1010 M⊙ in TNG-50, TNG-100, TNG-300, SIMBA (flagship), and Astrid hydrodynamical simulations. The results are compared with L-GalaxiesBH accounting for all the galaxies hosting a MBH (solid line) and galaxies hosting only MBHs with MBH > 106 M⊙. The shaded areas correspond to the 16th  −  84th percentiles.
In the text

	Fig. B.1. Duty cycle as a function of stellar mass for the different samples. Each row corresponds to a given redshift bin (from top to bottom: z = 0, 1, 3, 5). Points indicate median values, while error bars show the 16th  −  84th percentiles. For clarity, the different samples are slightly offset along the M* axis (by 0.05 dex) to avoid overlap.
In the text