Top-heavy stellar mass distribution in galactic nuclei inferred from the universally high abundance ratio of [Fe/Mg]

The existence of quasars and active galactic nuclei (AGNs) provides solid evidence of mass accretion onto supermassive black holes (SMBHs) (Lynden-Bell, 1969; Ho, 2008). Accretion disks powering bright AGNs are expected to be unstable by their self-gravity in the outer region at $\gtrsim 0.01$ pc and a good fraction of the gas is consumed by star formation (Paczynski, 1978; Shlosman & Begelman, 1987; Goodman, 2003; Inayoshi & Haiman, 2016). Theoretical studies have suggested that stars embedded in the self-gravitating regimes of AGN disks undergo gas accretion after their formation and thus the stellar mass distribution tends to be top-heavy (Goodman & Tan, 2004; Nayakshin et al., 2007; Mapelli et al., 2012). In fact, observations of the central nuclear region ($\sim 0.5$ pc) of the Milky Way have shown the existence of hundreds of massive, young OB stars and Wolf-Rayet stars (Krabbe et al., 1995; Genzel et al., 2003; Ghez et al., 2003; Levin & Beloborodov, 2003), and the absence of low-mass stars compared to the number expected with a standard Salpeter initial mass function (IMF) (Nayakshin & Sunyaev, 2005; Paumard et al., 2006; Alexander et al., 2007; Bartko et al., 2010).

One possible way to study the stellar mass function in nuclear regions of external galaxies is to measure chemical abundance of interstellar media (ISM) in AGNs since the production efficiency of heavy elements depends on the progenitor mass sensitively. Observations of high-$z$ quasars have extensively measured the line flux ratio of FeII and MgII emitted from the broad-line regions (BLRs) as a proxy of the abundance ratio of iron to $\alpha$ elements, ${\rm[Fe/\alpha]}$ (e.g., Dietrich et al., 2002; Iwamuro et al., 2002; Barth et al., 2003; Freudling et al., 2003; Kurk et al., 2007). In general, ${\rm[Fe/\alpha]}$ is expected to decline toward higher redshifts because $\alpha$ elements are predominantly produced in core-collapse supernovae (CCSNe), whereas iron is mainly released in type Ia SNe with a long delay time of $\gtrsim$ 0.5 Gyr (e.g., Matteucci & Greggio, 1986; Matteucci & Recchi, 2001; Maoz et al., 2010; Homma et al., 2015). However, quasar observations show that the Fe II/Mg II flux ratio does not decrease but appears to be almost constant up to $z\sim 7.5$, by which type Ia SNe hardly make a contribution to iron production (Iwamuro et al., 2002; Dietrich et al., 2003; Maiolino et al., 2003; De Rosa et al., 2011; Mazzucchelli et al., 2017; Shin et al., 2019; Onoue et al., 2020; Schindler et al., 2020; Yang et al., 2021). Moreover, Sameshima et al. (2017); Sameshima et al. (2020) have extensively analyzed the observed spectra of high-z quasars and found that the high Fe II/Mg II flux ratio is caused by overproduction of iron, namely an extremely high iron to magnesium abundance ratio ${\rm[Fe/Mg]}>0.2$. Therefore, those observations require the early production of iron in AGNs through a pathway different from normal CCSNe and type Ia SNe.

We propose a possible explanation to understand the origin of fully-established chemical environments in AGNs at a wide range of redshift, considering a top-heavy mass distribution of stars born in AGNs. This hypothesis leads to the frequent occurrence of massive CCSNe and pair-instability SNe (PISNe) and thus increases the ${\rm[Fe/Mg]}$ abundance ratio in the ISM because energetic explosions caused by heavier progenitors generally provide a larger amount of iron relative to magnesium in their ejecta. Such short-lived stellar progenitors establish iron-enriched BLRs in a short time of $<$ 10 Myr, which would naturally explain the apparent universality of the high flux ratio of Fe II/Mg II. In this paper, we quantify the shape of the stellar mass distribution that reproduces ${\rm[Fe/Mg]}>0.2$ as observed in quasar BLRs. Moreover, if a substantial fraction of those massive stars end up in binary BHs in galactic nuclei, their mergers would be detectable as gravitational-wave sources for both ground- and space-based interferometers such as advanced LIGO/Virgo and LISA/Tianqin.

The rest of the paper is organized as follows. We first introduce the stellar yields we assume in Section 2.1, and present the shape of the stellar IMF that explains the observed abundance ratio of ${\rm[Fe/Mg]}$ in Section 2.2. In Section 3, we model the history of mass growth of stars embedded in AGN disks and argue how the required top-heavy IMFs are established. In Section 4, the cosmic rate of binary BH mergers is calculated based on the top-heavy IMF hypothesis. In Section 5, we briefly discuss alternative explanations of the high-${\rm[Fe/Mg]}$ ratios observed in BLRs. Finally, our conclusion is summarized in Section 6. Throughout this paper, we adopt ${\rm[Fe/Mg]}=0.2$ as a reference value of the observed abundance ratio in quasars.

Luminous AGNs are powered by mass accretion through dense gaseous disks, where newly-born stars grow in mass via continuous gas accretion. In this section, we model the growth history of stars embedded in AGN disks so that their mass distribution explains the high ${\rm[Fe/Mg]}$ abundance ratio observed in AGNs (see Section 2.2).

The time evolution of the stellar mass distribution is calculated by solving the continuity equation of

where the total number of stars is fixed (i.e., the source term in the right hand side is set to zero) and $\dot{M}_{\ast}$ is the growth rate of a star with a mass of $M_{\ast}$. The growth rate of stars embedded in an AGN disk has not been understood properly owing to the uncertainties of stellar feedback (radiation and/or winds) and its impacts on the disk structure (Cantiello et al., 2021; Dittmann et al., 2021; Jermyn et al., 2021). To avoid numerous uncertainties in modeling of those effects, we characterize the mass accretion rate using a simple functional form with three free parameters:

where $\dot{M}_{\rm 0.1}$ is the mass growth rate at $M_{\ast}=0.1~{}{\rm M_{\odot}}$. In this model, low-mass stars with $M_{\ast}\ll M_{\rm c}$ undergo the Bondi-like accretion, $\dot{M}_{\ast}\propto M_{\ast}^{2}$, whereas for high mass stars with $M_{\ast}\gg M_{\rm c}$, the mass growth rate is moderated to $M_{\ast}^{\alpha}~{}(\alpha<2)$. This modification is motivated by the effects suppressing stellar growth owing to stellar radiation/mechanical feedback ($\alpha\simeq 1$; Eddington-limited growth), and the disk rarefication, gas shear motion, and tidal effects ($\alpha\lesssim 1$).

We initiate stars whose mass distribution follows a Salpeter IMF characterized with ${M_{\rm max}}=50~{}{\rm M_{\odot}}$ and $\Gamma=-2.35$. The maximum mass of the IMF is adopted so that the ratio of ${\rm[Fe/Mg]}$ becomes consistent with the upper bound for metal-poor stars in our Galaxy, ${\rm[Fe/Mg]}\simeq-0.2$ (see also the right panel of Fig. 3). We evolve the IMF by solving Eq. (4) for given three parameters and terminate its evolution at $t=5$ Myr, which corresponding to the lifetime of massive stars. We evaluate the difference between the resultant mass distribution ($\phi_{\rm cal}$) and a reference top-heavy IMF ($\phi_{\rm TH}$) that produces ${\rm[Fe/Mg]}\simeq 0.2$, calculating the $\chi^{2}$ value defined by

where the subscript $i$ indicates the $i$-th mass bin ($0.5\leq M_{i}/{\rm M_{\odot}}\leq 500$). We find the best-fitted parameter set of ($\dot{M}_{\rm 0.1}$, $M_{\rm c}$ and $\alpha$) to minimize the $\chi^{2}$ value with the Markov Chain Monte Carlo (MCMC) method.

We perform this fitting procedure for the two reference IMF models shown in Table 1. The top panel of Figure 4 presents the stellar mass distribution obtained for the IMF-A (red) and IMF-B (blue) from the MCMC realizations. For both the cases, the generated mass functions have two peaks at $M_{\ast}\lesssim 1$ and $\gtrsim 100~{}{\rm M_{\odot}}$, but the best-fit solutions (solid curves) reproduce shapes consistent with the reference IMFs. The $M_{\ast}$–$\dot{M}_{\ast}$ relation generated by the fitting method is shown in the bottom panel of Figure 4. The mean values with standard deviations of $\dot{M}_{\rm 0.1}$, $M_{\rm c}$ and $\alpha$, are also listed in Table 1. In the low mass regime ($M_{\ast}<10~{}{\rm M_{\odot}}$), the mass growth rates for the two cases are almost identical. Remarkably, the normalization of the mass growth rate at $M_{\ast}\lesssim 1~{}{\rm M_{\odot}}$ appears to be consistent with the Bondi accretion rate onto a star embedded in star-forming regions of an AGN disk, where the typical density and temperature of the gas is $n\sim 10^{7\mathcal{-}9}~{}{\rm cm^{-3}}$ and gas temperature $T\sim 10^{3\mathcal{-}4}~{}{\rm K}$, respectively (e.g., Cantiello et al., 2021). Above $M_{\ast}\sim 10~{}{\rm M_{\odot}}$, the best-fit models require the mass accretion rate to decrease with stellar mass. The critical mass for the transition is consistent with that predicted by analytical models for stellar growth in an AGN disk, mainly owing to stellar feedback by radiation and/or winds and the tidal effect caused by the central SMBHs (e.g., Fukushima et al., 2018; Dittmann et al., 2021). Moreover, at this mass range, the IMF-B requires a weaker suppression in $\dot{M}_{\ast}$ (i.e., a higher value of $\alpha$) to reproduce a larger fraction of stars with $M_{\ast}>100~{}{\rm M_{\odot}}$. This implies that the IMF slope and maximum mass depend sensitively on the nature of stellar mass growth and self-regulation mechanisms in the nuclear regions. Despite a large uncertainty of the power-law index $\alpha$ based on theoretical studies, our MCMC procedure suggests that a reasonable range of $-1.5\lesssim\alpha\lesssim 0$ successfully reproduces the reference top-heavy IMFs.

In summary, the top-heavy IMF that reproduces the chemical abundance of BLRs can be reasonably explained with the star-forming AGN disk model. We also claim that the $M_{\ast}$–$\dot{M}_{\ast}$ relation obtained in this paper gives a meaningful comparison to future theoretical studies for the stellar mass growth mechanism around AGNs.

As shown in Section 2.2, the high ${\rm[Fe/Mg]}$ abundance ratio seen in quasars implies that massive stars and remnant BHs are born more abundantly than those in normal star-forming environments. Stellar-mass BHs in AGN disks are expected to migrate inward due to gaseous dynamical friction and form close binaries with other BHs, and a fraction of them would merge via GW emission within a Hubble time. In this section, we quantitatively examine the merger rate of binary BHs (BBHs) that form in the galactic nuclei and compare it with the event rate inferred by LIGO/Virgo observations. In the following discussion, we consider the IMF-A as an example that leads to ${\rm[Fe/Mg]}>0.2$ and argue the quantitative difference in the BBH merger rate from the IMF-B case.

The merger rate of BBHs formed in an AGN disk (per unit time) is calculated by

where $\bar{\epsilon}_{\rm m}$ is the BH merger fraction, $\bar{f}_{\rm bh}$ is the number fraction of stars that end up as massive BHs for a given stellar mass function, and ${\rm SFR}_{\rm AGN}$ is the star formation rate (SFR) in the AGN disk. We set $\bar{\epsilon}_{\rm m}=0.5$ based on a recent study modeling binary hardening processes in AGN disks (Tagawa et al., 2020) and assume that those binaries coalesce immediately in a delay time much shorter than a Hubble time. To evaluate the SFR in the AGN, we suppose that the surface density of SFR is approximated as a single power-law distribution of $\dot{\Sigma}_{\rm\ast}=r^{-\gamma}$ and thus the SFR within a radius of $r$ is given by ${\rm SFR}(<r)\propto r^{2-\gamma}$ for $\gamma\neq 2$. Observationally, the SFR within $\sim 100~{}{\rm pc}$ is known to correlate with the BH accretion rate $\dot{M}_{\rm AGN}$ estimated from the radiative luminosity by assuming a 10% radiative efficiency ($\eta_{\rm rad}=0.1$) and the ratio is $\dot{M}_{\rm AGN}/{\rm SFR(<100~{}{\rm pc})}\sim O(0.1)$ (Esquej et al., 2014). Using the empirical relation, we estimate the SFR in the AGN disk as ${\rm SFR}_{\rm AGN}\sim 10\dot{M}_{\rm AGN}(r_{\rm d}/100~{}{\rm pc})^{2-\gamma}$, where $r_{\rm d}$ is the size of the AGN disk. We set the AGN disk size to $r_{\rm d}=10~{}{\rm pc}$, within which BHs can migrate inward and merge with other ones within $10^{8}$ yr, the typical lifetime of AGN disks (e.g, Stone et al., 2017; Dittmann & Miller, 2020). In this calculation, we assume $\gamma=1$, leading to ${\rm SFR}_{\rm AGN}\sim\dot{M}_{\rm AGN}(r_{\rm d}/10~{}{\rm pc})$, and thus consider that the cosmic star formation rate within $r_{\rm d}$ equals to the cosmic mass accretion rates onto the SMBH measured by Delvecchio et al. (2014). This is also motivated by the fact that star-forming galaxies tend to have nuclear stellar clusters as massive as the central SMBHs within a few pc from the galactic centers (Georgiev et al., 2016). Note here that with a larger value of $\gamma$, the SFR is weighted in the inner part of the AGN disk. We discuss how the choice of the index $\gamma$ affects the estimate of BH merger rates below.

Next, we divide the AGN disk into two regions. In the inner region at $r<r_{\rm in}=0.1~{}{\rm pc}$, which corresponds to the typical size of the BLRs of AGNs (e.g, Peterson, 1993; Kaspi et al., 2000), we assume that the stellar mass distribution follows the IMF-A owing to efficient stellar growth in the AGN disk (Thompson et al., 2005; Cantiello et al., 2021). As argued in Section 2.2, the mass function shape is adopted so that the ${\rm[Fe/Mg]}$ abundance ratio becomes as high as $\simeq 0.2$. In the outer region at $r_{\rm in}\leq r\leq r_{\rm d}$, we consider that the stellar population follows a Salpeter IMF characterized by $\Gamma=-2.35$ and ${M_{\rm max}}=50~{}{\rm M_{\odot}}$. Under these conditions, we calculate the BH formation rate in the AGN disk as

The number fraction of BHs is calculated by integrating the IMF over a mass range where BHs are left behind. The stellar yield model we adopt in this paper predicts that only 10% of the initial progenitor mass end up in their compact remnants as a consequence of enormous iron ejection from their collapsing cores (\al@Umeda2008ApJ, Nomoto2013ARA&A; \al@Umeda2008ApJ, Nomoto2013ARA&A). We here assume that the minimum stellar mass leaving a BH is $M_{\ast}=30~{}{\rm M_{\odot}}$ and thus the minimum BH mass is $M_{\rm\bullet,min}=3~{}{\rm M_{\odot}}$ (e.g., Özel et al., 2010; Spera et al., 2015). Massive stars with $M_{\ast}=140$–$260~{}{\rm M_{\odot}}$ are considered to take place PISNe without forming compact remnants. Note that while stars with $M_{\ast}>260~{}{\rm M_{\odot}}$ would leave BHs, the IMF-A and the Salpeter IMF we consider do not cover such an extremely massive population. Thus, over the mass range of $30\leq M_{\ast}/{\rm M_{\odot}}\leq 140$, we obtain $f_{\rm bh,in}=1.2\times 10^{-2}~{}{\rm M_{\odot}}^{-1}$ and $f_{\rm bh,out}=6.6\times 10^{-4}~{}{\rm M_{\odot}}^{-1}$ in the inner and outer region of the AGN disk, respectively. It is also worth noting that the 10% of the converting fraction associated with the HN and SLSN models leads to the mean remnant BH mass of $\bar{M_{\bullet}}=7.4~{}{\rm M_{\odot}}$ $(3.8~{}{\rm M_{\odot}})$ in the inner (outer) region. These values are compatible with the low-mass peak in the mass spectra of merging BHs inferred from LIGO-Virgo observations, but do not explain the rate for heavier BBHs at $M_{\bullet}\sim 30~{}{\rm M_{\odot}}$ (Abbott et al., 2021). Therefore, further mass growth processes of remnant BHs, such as continuous gas accretion and repeated mergers in dense AGN disks, are essential to simultaneously reproduce the chemical abundance ratio in BLRs and the mass distribution of merging BHs (Inayoshi et al., 2016; Safarzadeh & Haiman, 2020; Tagawa et al., 2020; Tagawa et al., 2021; Wang et al., 2021).

Finally, combining Eqs. (7) and (8), the cosmic merger rate of BBHs (in units of ${\rm Gpc^{-3}~{}yr^{-1}}$) is calculated as

where the mass accretion rate onto an SMBH in an AGN disk is replaced with the cosmic BH accretion rate density (BHAD) estimated from the AGN bolometric luminosity function of $dN/dL_{\rm AGN}$ as

We here adopt the function form of $\Psi_{\rm BHAD}(z)$ quantified by Delvecchio et al. (2014), where the AGN bolometric luminosity is calculated by a broad-band spectral energy distribution decomposition for obscured AGNs at $0<z<4$ found by Herschel and the Wide-field Infrared Survey Explorer (WISE).

Figure 5 shows the redshift-dependent BH merger rate at $0<z<3$ originating from the inner (orange) and outer (blue) regions of AGN disks. Since the two rates scale with the cosmic BHAD, the merger rates peak at the cosmic noon at $z\sim 2$ and decline toward lower redshifts. In the high-redshift universe ($z>1$), the BBH merger rate associated with the AGN outer regions is within the range inferred from GW observations in the LIGO/Virgo O3a run (the gray region; The LIGO Scientific Collaboration et al., 2021a). On the other hand, the BBH merger rate in the AGN inner regions, where the BBH formation efficiency is higher owing to a top-heavy mass distribution but the total SFR in the smaller area is lower, is about five times lower than that in the outer regions. This overall trend holds as long as the index is $\gamma\lesssim 1.5$. At $z<0.5$, the contribution from the two populations would explain about 10 % of the local merger rate of $\mathcal{R}_{\rm BBH}\simeq 19.1~{}{\rm Gpc}^{-3}~{}{\rm yr}^{-1}$ (The LIGO Scientific Collaboration et al., 2021a). Although the merger rates we calculate are lower than the LIGO/Virgo prediction in the local universe, this discrepancy could be resolved by taking into account a further delay-time of BH mergers, i.e., allowing a fraction of those BHs to coalesce even after host AGN disks have diminished.

Here, we briefly argue the case with the IMF-B. As shown in Table 1, the number fraction of BHs owing to stars with $M_{\ast}=30$–$140~{}{\rm M_{\odot}}$ for the IMF-B is 4 times lower than that for the IMF-A because the total number of stars per unit mass is lower for a top-heavier mass distribution. Therefore, the BBH merger rate in the inner part of AGN disks scales down accordingly from the rate expected for the IMF-A case in Figure 5. In the case of the IMF-B, we also consider BH formation from massive stars with $M_{\ast}>260~{}{\rm M_{\odot}}$, where ejection of stellar yields is negligible in the model we adopt. Therefore, the mass of BH remnants becomes as high as the progenitor mass in the pre-collapsing phase. As a result, intermediate massive BHs (IMBHs) with $M_{\bullet}\gtrsim 100~{}{\rm M_{\odot}}$ is considered to form from massive stars with $M_{\ast}>260~{}{\rm M_{\odot}}$, and the number fraction of IMBHs is estimated as $f_{\rm imbh}=1.1\times 10^{-3}~{}{\rm M_{\odot}}^{-1}$. This value leads to the cosmic merger rate of $\mathcal{R}_{\rm BBH}\sim 0.04~{}{\rm Gpc}^{-3}~{}{\rm yr}^{-1}$, which is slightly lower than the upper limit of IMBH merger rates of $0.056~{}{\rm Gpc}^{-3}~{}{\rm yr}^{-1}$ constrained by the LIGO/Virgo observations (The LIGO Scientific Collaboration et al., 2021b).

The model with the IMF-B yields another notable difference from the IMF-A case in the event rate of PISNe, whose progenitors are considered to be as massive as $M_{\ast}=140-260~{}{\rm M_{\odot}}$. Therefore, a stellar population that follows the IMF-B leads to a larger number of PISNe at a fraction of $f_{\rm pisn}=2.7\times 10^{-3}~{}{\rm M_{\odot}}^{-1}$ (the rightmost column in Table 1), while the IMF-A with $M_{\rm max}=130~{}{\rm M_{\odot}}$ does not produce a PISN ($f_{\rm pisn}=0$). Replacing $f_{\rm bh}$ with $f_{\rm pisn}$ in Eq. (9) and setting $\bar{\epsilon}_{\rm m}=1$, we estimate the PISN rate in AGN disks as $R_{\rm pisn}\sim 0.1~{}{\rm Gpc^{-3}~{}yr^{-1}}$ at $z=0.2$ and $\sim 3~{}{\rm Gpc^{-3}~{}yr^{-1}}$ at $z=2$. Those expected rates are consistent with the fact that there has not been a clear detection of PISNe. This gives an upper limit of the event rate at $\sim 1~{}{\rm Gpc^{-3}~{}yr^{-1}}$ in the local universe (Pan et al., 2012). It is worth noting that the upcoming observations with Vela C. Rubin Observatory’s LSST will enable us to explore PISN events up to $z\sim 2$ (Kasen et al., 2011), and could detect a few hundreds of PISNe per year if massive stars universally formed in AGNs following extremely top-heavy IMFs. Therefore, future observations of luminous transient events would be useful to constrain the stellar IMF in AGN disks.

In this section, we briefly mention two other possible scenarios that a high abundance ratio of ${\rm[Fe/Mg]}$ is achieved in AGN disks. The first one to increase the iron abundance of the interstellar medium is associated with type Ia SNe produced by thermonuclear explosions of white dwarfs (WDs) that exceed the Chandrasekhar mass limit owing to mass accretion through the AGN disk. If the growth timescale of WDs is fast enough to produce type Ia SNe in a shorter delay time since star formation ($\ll 1$ Gyr), the nuclear regions of AGNs would be quickly enriched by the yields composed of iron, enabling us to explain the high abundance ratio of ${\rm[Fe/Mg]}$ in quasars observed at $z>6$. For instance, the $M_{\ast}$–$\dot{M}_{\ast}$ relation shown in Fig. 5 yields $\dot{M}_{\ast}\sim 10^{-6}~{}{\rm M_{\odot}}~{}{\rm yr}^{-1}$ at $M_{\ast}=1~{}{\rm M_{\odot}}$, implying that WDs in a star forming AGN disk grow to the Chandrasekhar mass within a few Myr. However, the progenitors of those WDs would grow further in mass and become heavy enough to end up as NSs or BHs instead of WDs (e.g., Cantiello et al., 2021; Dittmann et al., 2021). Moreover, an analytical study by Pan & Yang (2021) claimed that accreting WDs spin up and thus prevent mass accretion before reaching the critical mass for the onset of thermonuclear explosions. Therefore, type Ia SNe induced by accreting WDs in dense AGN disks would not make a significant contribution to iron-enrichment of BLRs.

Another possible production of iron is through GRB events caused by massive stellar progenitors of $M_{\ast}>260~{}{\rm M_{\odot}}$. While such massive stars are considered to directly collapse to BHs at the end of their lifetime without producing iron-rich yields, relativistic jets associated with GRB-like explosions, if any, would release a large fraction of newly synthesized metals into the interstellar medium. Ohkubo et al. (2006) has conducted nucleosynthesis calculations for massive stars with $M_{\ast}=500$ and $1000~{}{\rm M_{\odot}}$ and found that massive GRB-like explosions preferentially produce iron-rich ejecta with ${\rm[Fe/Mg]}\sim 0.5$. Therefore, an extremely top-heavy IMF that leads to a large number of massive GRBs would be a possible solution to reproduce the high abundance ratio of ${\rm[Fe/Mg]}>0.2$ seen in high-$z$ quasars. However, the stellar yield model owing to massive GRBs depends sensitively on the explosion energy and the jet opening angles. We leave this issue for future investigation and make a robust conclusion on chemical enrichment in quasars using yield models calibrated by various SN observations.

In this paper, we argue the properties of the stellar mass distribution in galactic nuclear regions to explain the high iron to magnesium abundance ratio of ${\rm[Fe/Mg]}\gtrsim 0.2$, which is observed in broad-line regions (BLRs) of active galactic nuclei (AGNs). Nuclear synthesis models suggest that massive explosive events such as super-luminous supernovae (SLSNe) with progenitor mass of $M_{\ast}\sim 60$–$140~{}{\rm M_{\odot}}$ and pair-instability supernovae (PISNe) with $M_{\ast}\sim 220$–$260~{}{\rm M_{\odot}}$ produce highly iron-enriched ejecta and effectively enhance ${\rm[Fe/Mg]}$ of the interstellar medium. We calculate the mass-integrated stellar yields varying the maximum mass ${M_{\rm max}}$ and power-law index $\Gamma$ of the stellar initial mass function (IMF), and find that the iron-enriched environments with ${\rm[Fe/Mg]}\gtrsim 0.2$ form when the IMF is characterized with ${M_{\rm max}}\sim 100$–$150~{}{\rm M_{\odot}}$ and $\Gamma\gtrsim-1$ or ${M_{\rm max}}\gtrsim 250~{}{\rm M_{\odot}}$ and $\Gamma\gtrsim 0$. This suggests that massive stars preferentially form in galactic nuclei and promote chemical enrichment in BLRs.

Next, we discuss how those top-heavy IMFs are established in the galactic nuclear regions. Following a theoretical model of massive star formation via rapid accretion through dense AGN disks, we calculate the time-evolution of the stellar mass distribution and model the stellar mass growth rate so that the abundance ratio of ${\rm[Fe/Mg]}\gtrsim 0.2$ observed in AGNs is explained by their explosive events (see Eqs. 4 and 5). The model requires Bondi-like stellar growth that is regulated by radiative/mechanical stellar feedback at $M_{\ast}>10~{}{\rm M_{\odot}}$.

Finally, we calculate the cosmic merger rate of binary BHs formed in AGN disks, assuming the top-heavy IMF in the inner BLRs ($r<0.1~{}{\rm pc}$) to reproduce ${\rm[Fe/Mg]}>0.2$ and a Salpeter IMF in the outer part ($0.1~{}{\rm pc}<r<10~{}{\rm pc}$). We find that the merging BH population in the outer region is responsible for the total merger rate, which is generally consistent with the rate inferred by the LIGO and Virgo observations at $z\sim 1$–$3$. The BH merger rate in AGN disks declines toward lower redshifts, as the cosmic AGN activity fades out, but still accounts for about 10 % of the total merger rate observed in the local universe. We also suggest that top-heavy IMFs achieved in the inner region do not necessarily promote formation of massive BHs with $M_{\bullet}\gtrsim 10~{}{\rm M_{\odot}}$ since enormous mass ejection from their HNe and SLSNe is required to reproduce the high iron to magnesium abundance ratio in BLRs. Thus, further mass growth mechanisms of remnant BHs, such as continuous gas accretion and repeated BH mergers in AGN disks, would be essential to establish a high-mass peak around $M_{\bullet}\sim 30~{}{\rm M_{\odot}}$ in the observed mass spectra of merging BHs. In addition, the top-heavy hypothesis would require a high detection rate of PISNe at $z\lesssim 2$ in the upcoming LSST observations. This will enable us to put a meaningful constraint on the stellar mass distribution in the nuclear regions of AGNs.

We thank K. Hotokezaka and M. Onoue for useful discussions. This work was supported in part by JSPS KAKENHI Grant Numbers 17H06363, 17K14249, 20H00179, 20H05855, 21H04499, 21H04997, and 21K20378. K. I. acknowledges support from the National Natural Science Foundation of China (12073003, 12003003, 11721303, 11991052, 11950410493), the National Key R&D Program of China (2016YFA0400702).

The data underlying this article will be shared on reasonable request to the corresponding author.