The stellar remnants of high redshift nuclear starburst discs: a potential origin for nuclear star clusters?

TL;DR

This study models the stellar remnants of nuclear starburst discs at low redshift, suggesting they could be the origins of nuclear star clusters observed in galaxies today, especially those with large supermassive black holes.

Contribution

It provides a detailed comparison of simulated stellar remnants of NSDs with observed nuclear star clusters, supporting the hypothesis that NSDs at high redshift are progenitors of present-day NSCs.

Findings

Good match between simulated and observed NSC properties.

Approximately 20% of remnants have sizes consistent with observed NSCs.

Remnants follow similar size-luminosity relations as observed in nearby NSCs.

Abstract

Nuclear starburst discs (NSDs) are very compact star-forming regions in the centers of galaxies that have been studied as a possible origin for the absorbing gas around a central active galactic nucleus. NSDs may be most relevant at $z\sim 1$ when obscured accretion onto supermassive black holes (SMBHs) is common. This paper describes the characteristics of the stellar remnants of NSDs at $z=0.01$, taking into account the evolution from $z=1$. Using a stellar synthesis model, the colours, masses, and luminosities of the stellar remnants are computed for a suite of 192 two-dimensional NSD models. These properties are compared to observations of local nuclear star clusters (NSCs), and a good match is found between the predicted and observed properties. Dynamical effects will likely cause the final remnant to be a rotating, nearly spherical distribution. In addition, $\approx 20$% of the NSD remnants have half-light radii <~ 10 pc, consistent with NSCs hosted in both late-type and early-type galaxies, and all the remnants follow similar size-luminosity relationships as observed in nearby NSCs. NSDs require the presence of a central SMBH and the most massive and compact stellar remnants are associated with the most massive SMBHs, although stellar clusters with a variety of sizes can be produced by all considered SMBH masses. Overall, NSDs at $z\sim 1$ appear to be a promising origin for the $\gg 1$ Gyr NSC population in early- and late-type galaxies with large SMBHs.

Tables

Table 1. Table 1: All the bins of input parameters for NSDs modeling are listed: M ​ bh 𝑀 bh M\textsubscript{\text{bh}} , R ​ out 𝑅 out R\textsubscript{\text{out}} , f ​ g,out 𝑓 g,out f\textsubscript{\text{g,out}} , and m 𝑚 m . A total of 192 models are computed with all the possible combination of these bins (Gohil & Ballantyne, 2017 ) .

log($M\text{<sub><span>bh</span></sub>}/$M⊙)	$R\text{<sub><span>out</span></sub>}$(pc)	$f\text{<sub><span>g,out</span></sub>}$	$m$
6.5	240	0.2	0.1
7.0	180	0.4	0.3
7.5	120	0.6	0.5
8.0	60	0.8	–

Table 2. Table 2: Median properties of the 192 NSDs and their remnant star clusters. The 192 models are separated into 4 samples: (i) all models, (ii) those without pc-scale starbursts, (iii) those with pc-scale bursts and R ​ 1/2 > 10 𝑅 1/2 10 R\textsubscript{\text{1/2}}>10  pc, and (iv) those with pc-scale bursts and R ​ 1/2 < 10 𝑅 1/2 10 R\textsubscript{\text{1/2}}<10  pc. The columns are the number of models in the sample, the dynamical time (Eq.  1 ), the half-light radius (Eq.  2 ), the remnant mass at z = 0.01 𝑧 0.01 z=0.01 and its V-band absolute magnitude.

Sample	Number	$\tau \text{<sub><span>dyn</span></sub>}$(Myr)	$R\text{<sub><span>1/2</span></sub>}$(pc)	$\log (M_{\ast }/M_{\odot }$)	$M_V$(mag)
All models	192	0.92	50	6.3	$-9.6$
No pc-scale burst	88	0.98	87	6.1	$-9.1$
pc-scale burst ($R\text{<sub><span>1/2</span></sub>}>10$ pc)	69	0.98	43	6.6	$-10.3$
pc-scale burst ( pc)	35	0.4	1.2	6.3	$-9.7$

Equations

$$\tau \text{textsubscript}\text{dyn}\equiv \frac{1}{\Omega \text{textsubscript}\text{out}},$$

$${\int }_{R\text{textsubscript}\text{in}}^{R\text{textsubscript}\text{1/2}}{\int }_{{\lambda }_{\text{min}}}^{{\lambda }_{\text{max}}}f(\lambda ,r)drd\lambda =\frac{1}{2}{\int }_{R\text{textsubscript}\text{in}}^{R\text{textsubscript}\text{out}}{\int }_{{\lambda }_{\text{min}}}^{{\lambda }_{\text{max}}}f(\lambda ,r)drd\lambda ,$$

$$f(r)=\frac{4\pi^{2}hc^{2}r}{\lambda^{5}\exp{\Big{[}\frac{hc}{\lambda k_{B}T\textsubscript{\text{eff}}(r)}}-1\Big{]}}$$

$$⟨{\dot{M}}_{\ast }⟩=\frac{{\int }_{R\text{textsubscript}\text{in}}^{R\text{textsubscript}\text{out}}{\dot{M}}_{\ast }(r)2\pi rdr}{{\int }_{R\text{textsubscript}\text{in}}^{R\text{textsubscript}\text{out}}2\pi rdr}.$$

$$\dot{M}_{*}(t)=A\exp\Bigg{[}-\Bigg{(}\frac{(t-\mu)}{\sqrt{2}\sigma}\Bigg{)}^{2}\Bigg{]},$$

Full text

††pubyear: 2019††pagerange: The stellar remnants of high redshift nuclear starburst discs: a potential origin for nuclear star clusters?–The stellar remnants of high redshift nuclear starburst discs: a potential origin for nuclear star clusters?

1 Introduction

Nuclear starburst discs (NSDs) — extremely compact, Eddington-limited starbursts in the cores of galaxies — have been proposed as a potential method for fueling and obscuring active galactic nuclei (AGNs; Ballantyne 2008; Gohil & Ballantyne 2017). A key prediction of both 1-dimensional (1D) and 2-dimensional (2D) NSD models is the presence of a pc-scale burst of star-formation caused by the loss of dust opacity which inflates the disc atmosphere and results in significant obscuration along most lines-of-sight to the central AGN (Thompson, Quataert & Murray, 2005; Ballantyne, 2008; Gohil & Ballantyne, 2017). This pc-scale starburst is more likely to occur when galaxies have high gas fractions, and therefore NSDs are a possible candidate for obscuring AGNs at $z\sim 1$ and influencing the shape of the X-ray Background (Gohil & Ballantyne, 2018).

If NSDs are obscuring AGNs at high redshifts, then the remnants of these discs may be detectable in the nuclei of local galaxies. Interestingly, over the last several years Hubble Space Telescope (HST) observations have revealed that nuclear star clusters (NSCs) are found at the center of most ($>70$%) nearby early and late-type galaxies (e.g., Carollo et al., 1997; Matthews et al., 1999; Böker et al., 2002; Côté et al., 2006; Seth et al., 2008a; Georgiev & Böker, 2014; Carson et al., 2015; Georgiev et al., 2016). The effective radii (i.e., the half-light radii) of many NSCs is smaller than a few tens of pc (e.g., Böker et al., 2004; Côté et al., 2006; Turner et al., 2012; Georgiev et al., 2016), and they show evidence for a complex star-formation history with the presence of multiple generations of stars (e.g., Walcher et al., 2005; Rossa et al., 2006; Seth et al., 2006, 2010; Lyubenova et al., 2013; Kacharov et al., 2018), including an older population with ages $\gg 1$ Gyr. In addition, NSCs are also observed to be rotating (Seth et al., 2008b), and many of them co-exist with a central SMBH (e.g., Seth et al., 2008a; Graham & Spitler, 2009; Georgiev et al., 2016) in their host galaxies. Two broad class of mechanisms have been suggested for the formation of NSCs. The first focuses on the migration of dense clusters via dynamical friction (e.g., Andersen et al., 2008; Antonini, 2013; Arca-Sedda & Capuzzo-Dolcetta, 2014), while the second one considers in situ nuclear star-formation from gas inflows (e.g., Milosavljević, 2004; Seth et al., 2006). Given the wide range of NSC properties observed in both early- and late-type galaxies, both of these mechanisms are likely relevant and operating with various efficiencies in different galaxies (e.g., Antonini, Barausse & Silk, 2015; Guillard, Emsellem & Renaud, 2016).

A very compact ($\sim 10$ pc), rotating, multi-Gyr-old stellar population associated with a SMBH is very suggestive as a potential remnant from a NSD that was obscuring an AGN at $z\sim 1$. However, to test this idea the stellar content of a NSD must be evolved to the present day to compare with observations of nearby NSCs. This paper presents the result of such an experiment. The next section provides a brief overview of the NSD models, and then Section 3 describes the procedure to compute various properties of the stellar remnants of NSDs. The predicted properties of the resulting star clusters are presented in Section 4, where a comparison to observed NSCs is also performed. Finally, Section 5 contains a short discussion and a summary of the results. A standard $\Lambda$CDM model cosmology is assumed with $H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{\Lambda}=0.72$, and $\Omega_{m}=0.28$ (Hinshaw et al., 2013).

2 Overview of Nuclear Starburst Disc Models

The 1D theory of NSDs was first presented by Thompson et al. (2005), and was extended to 2D by Gohil & Ballantyne (2017). Those papers present a complete description of the model, but a brief overview of the fundamental properties and assumptions of the 2D model is presented here.

The 2D NSD models used in this paper are time-independent, equilibrium calculations of a Keplerian gas disc orbiting in the potential of a SMBH and galactic bulge. The disc midplane is always assumed to be critically unstable to star-formation (i.e., the Toomre $Q$ parameter is very close to one), and the radiation pressure of starlight on dust is fundamental to the vertical support of the disc. The calculation does not compute the formation of individual stars, but rather determines the radial and vertical structure of the NSD by finding the radial distribution of the star-formation rate that leads to hydrostatic, thermal and energy balance.

An interesting aspect of the NSD theory is that gas is assumed to slowly accrete through the disc with a radial velocity $mc_{s}$, where $c_{s}$ is the local sound speed and $m$ is a constant Mach number. This results in a constant competition between star-formation and accretion at every radius. If gas reaches a radius $\sim$pc from the SMBH, then the temperature of the gas may be high enough to sublimate dust. This results in a large burst of star-formation as well as a vertical expansion of an atmosphere. The vertically inflated atmosphere is supported by the radiation pressure from the bursty star-formation (Thompson et al., 2005; Gohil & Ballantyne, 2017). Gohil & Ballantyne (2017) computed the two-dimensional structure of 192 NSDs across its input parameter space and showed that 52% of these models possess the starburst phenomenon at the parsec-scale.

The structure of an individual NSD depends on four input parameters: the disc size ($R\textsubscript{\text{out}}$), the gas fraction at the outer radius ($f\textsubscript{\text{g,out}}$), the Mach number ($m$), and the SMBH mass ($M\textsubscript{\text{bh}}$). The distribution of parameters used in the Gohil & Ballantyne (2017) model set is shown in Table 1, and were not chosen nor optimised for describing local NSCs. Therefore, the NSD remnants derived from this parameter set will provide an unbiased sample to compare against the observed NSC properties.

Each of these 192 models calculated from the different combination of parameters yields a radial distribution of star-formation rates (SFR) needed to vertically support the NSD. The distribution of SFRs can vary widely among the models depending on the input parameters. Two examples of the derived equilibrium SFRs are shown in Figure 1.

The profile with the black stars exhibits a pc-scale burst of star-formation resulting from the loss of dust opacity as described above, while the SFR distribution shown with the cyan circles has most of the star-formation at the outer edge of the disc. In general, all of the 192 NSD models fall into one of the two profiles shown in Fig. 1, but with significant variations in the magnitude of the SFRs (e.g., Ballantyne, 2008; Ballantyne, Armour & Indergaard, 2013). These 192 time-independent SFR profiles derived from the NSD models are then used to derive the predicted properties of the relic star cluster long after the end of the starburst.

3 Calculation of Properties of Remnant Star Cluster

3.1 Timescales and Half-light Radius

As the SFR profiles computed by Gohil & Ballantyne (2017) are time-independent, a characteristic timescale must be identified in order to predict the final star cluster properties. Here, we use the dynamical time of the NSD, $\tau\textsubscript{\text{dyn}}$, which appears to explain observed star-formation efficiencies in both low-$z$ and high-$z$ environments (e.g., Silk, 1997; Elmegreen, 1997; Kennicutt, 1998; Genzel et al., 2010). The dynamical time is computed for all 192 NSD models by

[TABLE]

where $\Omega\textsubscript{\text{out}}$ is the Keplerian angular frequency at $R\textsubscript{\text{out}}$.

If the NSD structure does not change significantly over $\tau\textsubscript{\text{dyn}}$, then the final distribution of light in the remnant cluster will follow the SFR profile predicted by the NSD theory (e.g., Fig. 1). Therefore, the half-light radius of the final remnant, $R\textsubscript{\text{1/2}}$ (also known as the effective radius), can be predicted directly from the NSD model. As described by Thompson et al. (2005) and Gohil & Ballantyne (2017), the total radiative flux due to star-formation is known at each radius in a NSD, so the bolometric spectrum can be computed by assuming each annulus emits as a blackbody with an effective temperature $T\textsubscript{\text{eff}}(r)$. Then $R\textsubscript{\text{1/2}}$ for each NSD is given by the following equation:

[TABLE]

where

[TABLE]

and $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$ are computed based on the maximum and minimum $T\textsubscript{\text{eff}}(r)$. This bolometric $R\textsubscript{\text{1/2}}$ can then be compared to measurements of NSCs in the local Universe since the observed values depend weakly on wavelength (e.g., Georgiev & Böker, 2014).

Figure 2 plots $R\textsubscript{\text{1/2}}$ and $\tau\textsubscript{\text{dyn}}$ for all 192 NSD models considered here, including those with (black stars) and without (cyan circles) the pc-scale burst of star-formation.

The figures show that the NSDs without pc-scale starbursts all have $R\textsubscript{\text{1/2}}\ga 20$ pc, but that the ones with the small scale bursts are clustered in two groups. There is one group of black stars with intermediate sizes ($\sim$10s of pc), and a second group with $R\textsubscript{\text{1/2}}\sim 1$ pc. The first group are models where there is roughly equal amounts of star-formation in the outer disc and in the pc-scale burst, leading to a $R\textsubscript{\text{1/2}}$ that is roughly in the middle of the disc. The NSDs with the small $R\textsubscript{\text{1/2}}$ are ones where the SFR in the pc-scale burst overwhelms all other star-formation in the disc. The dynamical times are $<3$ Myrs for all models and are spread in vertical stripes because the Keplerian frequency of a disc at the outer radius depends only on $R\textsubscript{\text{out}}$ and $M\textsubscript{\text{bh}}$. For a fixed pair of $R\textsubscript{\text{out}}$ and $M\textsubscript{\text{bh}}$, there are 12 models exhibiting the same $\tau\textsubscript{\text{dyn}}$ due to degeneracy into other parameters. Interestingly, the plot shows that the NSDs with pc-scale bursts are found at all $\tau\textsubscript{\text{dyn}}$, with a preference for the small $R\textsubscript{\text{1/2}}$ NSDs to be found at the smaller $\tau\textsubscript{\text{dyn}}$, corresponding to lower $M\textsubscript{\text{bh}}$ and/or smaller disc sizes.

To help determine if a smaller $M\textsubscript{\text{bh}}$ or $R\textsubscript{\text{out}}$ is more important in producing the pc-scale bursts, Figure 3 plots $R\textsubscript{\text{1/2}}$ versus the $M\textsubscript{\text{bh}}$ at the center of the NSD.

Compact NSDs with $R\textsubscript{\text{1/2}}<10$ pc are seen at all $M\textsubscript{\text{bh}}$, but are actually more common at larger $M\textsubscript{\text{bh}}$ with 15 at $M\textsubscript{\text{bh}}=10^{8}$ $M_{\odot}$, as compared to only 4 at $M\textsubscript{\text{bh}}=10^{6.5}$ $M_{\odot}$. This implies that a larger $M\textsubscript{\text{bh}}$ is a stronger factor than $R\textsubscript{\text{out}}$ in producing the most powerful pc-scale bursts, in agreement with the results of the 1D calculations (Ballantyne, 2008). More massive black holes cause hotter and denser NSDs which more readily lead to dust sublimation and a small scale burst of star-formation. However, Fig. 3 also shows that the smallest $R\textsubscript{\text{1/2}}$ are found in NSDs around the lowest mass SMBHs, so the size of the stellar remnant will be weakly correlated with $M\textsubscript{\text{bh}}$.

3.2 Mass, Luminosity and Colours

The predicted mass and spectral properties of the remnant cluster at low redshift can only be calculated by assuming a star-formation history for the NSD models and evolving the stellar population in time. As before, the equilibrium radial structure computed for each NSD is assumed to be fixed during its lifetime. The radially-dependent SFR (e.g., Fig. 1) is then averaged to yield $\langle\dot{M}_{*}\rangle$ for each model, i.e.,

[TABLE]

The star-formation history of each NSD model is then assumed to be a narrow Gaussian in time with a lifetime $\tau\textsubscript{\text{dyn}}$:

[TABLE]

where $A=\langle\dot{M}_{*}\rangle\tau\textsubscript{\text{dyn}}/\sigma\sqrt{2}\pi$, $\mu=t_{0}+\tau\textsubscript{\text{dyn}}/2$, and $\sigma=\tau\textsubscript{\text{dyn}}/5$. The normalization constant ($A$) is chosen so that the stellar mass formed in the NSDs is equal to $\langle\dot{M}_{*}\rangle\tau\textsubscript{\text{dyn}}$. The time $t_{0}$ corresponds to the age of the Universe at $z=1$ (i.e., $5.9$ Gyrs), when such NSDs may be common due to larger galaxy gas fractions (Gohil & Ballantyne, 2018).

The star cluster formed by this star-formation history is then evolved forward in time using Python-FSPS, the Python wrapper of the Flexible Stellar Population Synthesis code (FSPS v3.0) developed by Conroy, Gunn & White (2009) and Conroy & Gunn (2010). All FSPS calculations assume solar metallicity and a Kroupa (2001) initial mass function. The Kriek & Conroy (2013) model is used for the attenuation by dust around old stars. The characteristics of the remnant star clusters are computed assuming an observed redshift $z\textsubscript{\text{obs}}=0.01$, giving an age of the stellar population of $7.7$ Gyr111The properties of the remnant cluster show negligible differences for other starting redshifts (e.g., $z=0.8$ or $1.5$) that yield final cluster ages $\gg 1$ Gyr.. Each FSPS calculation yields the survived stellar mass at $z=0.01$, $M_{\ast}$, and absolute magnitudes in a number of HST and Johnson filters.

3.3 Dynamical Evolution of $R\textsubscript{\text{1/2}}$

The calculation of $R\textsubscript{\text{1/2}}$ described above (Eq. 2) is based on the radial distribution of stars in each model NSD. It is possible that dynamical effects could significantly alter this radial distribution in the $\approx 8$ Gyrs between the NSD event and $z=0.01$. Therefore, it is important to consider how dynamics within the remnant stellar discs may impact the values of $R\textsubscript{\text{1/2}}$ before comparing them to observed NSCs at low redshifts.

The stellar cluster produced by a NSD will be a rotating disc spanning from $\sim 0.1$ pc to $\sim 100$–$200$ pc from the galactic center. The dynamical effects relevant to stars in the disc are nicely summarized by Kocsis & Tremaine (2011) in their Figure 1. It shows that a number of dynamical processes have timescales below $\sim 8$ Gyrs for discs with these sizes, including standard two-body relaxation, disc eccentricity relaxation, Newtonian precession, relativistic precession, vector resonant relaxation (VRR), and warps by a massive perturber (e.g., a molecular torus). With larger SMBH masses ($\sim 10^{8}M_{\odot}$), Lense-Thirring precession timescale can also be shorter than $\sim 8$ Gyrs in the inner edge of the disc.

However, most of these dynamical processes do not change the radial distribution of the stars significantly. Specifically, Newtonian precession and relativistic precession only change the pericenter orientation of the stellar orbits in their orbital planes. Similarly, VRR, Lense-Thirring precession and the warp of the disc redistribute orbital plane orientations of stars in the disc, and have only a minimal effect on the semi-major axes and the value of $R\textsubscript{\text{1/2}}$. In particular, a recent numerical study of VRR in stellar discs by Szöglén & Kocsis (2018) showed an interesting mass dependence in the final stellar distribution. Their simulations found that massive young stars were in a warped disc, but the less massive older stars were in a rotating, mostly spherical distribution. If VRR is important for the NSD remnants, then these results imply that at $z=0.01$ the stellar cluster will more closely resemble a rotating spheroidal system.

The effects of non-resonant relaxation (standard two-body relaxation and disc eccentricity relaxation) on stellar discs have been investigated numerically by Šubr & Haas (2014). These authors find that the relaxation processes scatter stars both inwards and outwards, which causes peaks in the semi-major axis distribution to be smoothed out to both higher and lower values. This numerical result can be rescaled for different SMBH masses and disc sizes, as noted by Šubr & Haas (2014). Therefore, even in cases where there is a peak in stellar density, such as in NSDs with $R\textsubscript{\text{1/2}}\la 10$ pc, there will be no systematic shift in stellar density over time.

The stellar discs produced by NSDs may also be subject to spiral-arm instabilities (e.g., Lee & Goodman, 1999) at different points in its evolution. These instabilities cause mixing of stellar orbits as they pass through the gas, but, as with two-body relaxation, roughly equal numbers of stars are scattered inwards as outwards and there is no significant change in the surface density profile of the disc (Sellwood & Binney, 2002).

From the above discussion we conclude that the values of $R\textsubscript{\text{1/2}}$ computed by eq. 2 will be an adequate estimate of the effective radii of the stellar remnants of NSDs at $z=0.01$ and can be compared to the measurements of nearby NSCs. In addition, the effects of VRR will likely produce a rotating, quasi-spherical remnant at low redshift.

4 Comparison to Low Redshift Nuclear Star Clusters

In this section the computed properties of the NSD remnants at $z=0.01$ are compared to observations of NSCs observed by HST in the local Universe. For reference, Table 2 provides the median $M_{\ast}$, $R\textsubscript{\text{1/2}}$, $\tau\textsubscript{\text{dyn}}$, and absolute V-band magnitude, $M_{V}$, for four different samples of the NSD models.

We begin in Figure 4 by comparing $R\textsubscript{\text{1/2}}$ and $M_{\ast}$ computed from the NSD models to a large number of NSC observations collected from the literature (Carollo et al., 1997; Carollo, Stiavelli & Mack, 1998; Böker et al., 2004; Seth et al., 2006, 2010; Leigh, Böker & Knigge, 2012; Georgiev et al., 2016).

The majority of the NSC masses are computed photometrically using mass-to-light ratios (e.g., Georgiev et al., 2016), which could lead to an underestimate of the mass if a young stellar population was dominating the observed flux. However, dynamical NSC mass measurements (e.g., Seth et al., 2008a; Nguyen et al., 2018) agree with the photometric ones to within a factor of about $2$, indicating that there is no large systematic error with the photometric masses. Therefore, the NSC masses plotted in Fig. 4 are likely good estimates of the mass of the older stellar population in NSCs, and the comparison to the NSD remnant masses will be valid.

Figure 4 shows that $M_{\ast}$ for all 192 models is in good agreement with the majority of the observed masses of local NSCs. In fact, the median remnant mass of $\log(M_{\ast}/M_{\odot})=6.3$ predicted by the NSD models (Table 2) is in exact agreement with the median mass of the NSCs shown in the figure. However, there are disagreements at the edges of the mass distribution. In particular, the maximum NSD remnant mass is $M_{\ast}\approx 3\times 10^{7}$ M*⊙*, more than an order of magnitude smaller than the maximum observed NSC masses. This supports the view that some NSCs are likely built up by a variety of processes over time, and a high-$z$ NSD may only be one of the possible mechanisms. It is also possible that the mismatch in maximum mass is due to our limited sampling of the potential parameter space and considering only a single NSD event during the lifetime of a galaxy. Overall, the strong agreement between $M_{\ast}$ and the observed NSC masses shows that the NSDs have the SFRs and lifetimes necessary to explain the old population of NSCs.

In contrast to the masses, most of the NSD remnants have $R\textsubscript{\text{1/2}}$ significantly larger than the observed NSCs. The median predicted $R\textsubscript{\text{1/2}}$ for the NSD remnants is $50$ pc (Table 2), an order of magnitude larger than the median of the observational sample ($R\textsubscript{\text{1/2}}=5$ pc). However, Fig. 4 shows that the group of NSD remnants with $R\textsubscript{\text{1/2}}<10$ pc are in very good agreement with the properties of NSCs observed in predominantly late-type galaxies (LTGs; red ‘X’s; Georgiev et al. 2016) with median $R\textsubscript{\text{1/2}}$ and $M_{\ast}$ very similar to the observations (last line of Table 2). These compact star clusters are predicted to arise from NSDs with the strongest pc-scale starbursts (e.g., Fig. 1). The scarcity of models at $3\ \mathrm{pc}\la R\textsubscript{\text{1/2}}\la 10\ \mathrm{pc}$ is likely due to the choice of model parameters. Fig. 3 shows that the smallest $R\textsubscript{\text{1/2}}$ values obtained from the NSD models increase with the central $M\textsubscript{\text{bh}}$ mass. Therefore, NSD remnants with $3\ \mathrm{pc}\la R\textsubscript{\text{1/2}}\la 10\ \mathrm{pc}$ would be produced from models with $M\textsubscript{\text{bh}}=10^{8.5}$ M*⊙* and $10^{9}$ M*⊙* and would closely match the observed NSC sizes in early-type galaxies (ETGs; blue ‘Y’s; Georgiev et al. 2016). We conclude that NSDs with the strongest pc-scale starbursts can produce remnants with $R\textsubscript{\text{1/2}}$ similar to those observed in local galaxies.

Figure 4 shows that the majority of NSD remnants have $R\textsubscript{\text{1/2}}\ga 30$ pc which is consistent with only a few observed NSCs hosted in both ETGs and LTGs. In general, such large $R\textsubscript{\text{1/2}}$ values will be a common outcome from NSD models due to the relatively specific conditions needed to produce a strong pc-scale starburst (Ballantyne, 2008). This result implies that only the fraction of high-$z$ NSDs that produce the strongest smaller-scale bursts are a potential origin for the old population observed in local NSCs.

Another important observational property of local NSCs is a size-luminosity relation (Georgiev & Böker, 2014). Figure 5 compares this observed relation (red line; Georgiev & Böker 2014) with the predictions of the NSD remnants. The shaded uncertainty region encompasses the majority of the observed NSCs from Georgiev & Böker (2014).

As seen in Fig. 4, most of the NSD remnants have $R\textsubscript{\text{1/2}}$ too large for their luminosity, but there is a group of $\approx 40$ remnants that are predicted to lie within or very close to the observed size-luminosity relationship.

It is clear that all three sub-groups of NSD remnants follow a size-luminosity relationship of the form $\log R\textsubscript{\text{1/2}}\propto\alpha M_{v}$, with $\alpha$ measured to be $\approx-0.11$. This value is roughly $2\times$ smaller than the observed slope of $\alpha=-0.25\pm 0.02$ found by Georgiev & Böker (2014). As the predicted slope does not take into account the effect of the local galaxy SMBH function on the weighting of the model parameters, the agreement of the two slopes is considered satisfactory. The values of $M_{V}$ predicted for the NSD remnants are also well matched with the observed values (median $M_{V}=-9.6$; Table 2) which is not surprising since the masses are in good agreement and the luminosities are calculated using standard mass-to-light ratios (Conroy, Gunn & White, 2009).

The NSD models require the presence of a SMBH, and, as seen in Sect. 3, the more massive the central SMBH, the more likely a compact starburst will develop at high-$z$. It is therefore interesting to consider the relationship between the observed properties of the predicted NSD remnant and the SMBH. Figure 6 plots $M\textsubscript{\text{bh}}/M_{\ast}$ as a function of $M_{\ast}$ for the NSD remnants.

The four stripes on the figure correspond to the four SMBH masses used in the calculations (Table 1), with the values decreasing from $M\textsubscript{\text{bh}}=10^{8}$ M*⊙* to $10^{6.5}$ M*⊙* as one goes from top to bottom in the plot. The most massive and therefore the most luminous stellar remnants are produced by NSDs around the most massive SMBHs. However, these SMBHs can also produce much less massive remnants. Indeed, except at the most massive end, the NSD theory admits stellar remnants with similar masses and sizes to local NSCs at all SMBH masses. Measurements of NSC and SMBH masses from the same galaxies are relatively rare, but current data shows that ratios $M\textsubscript{\text{bh}}/M_{\ast}>1$ are frequently found in ETGs, but are observed for only a few LTGs (Georgiev et al., 2016; Nguyen et al., 2018).

As a final sanity check, the WFPC2 colours F606W$-$F814W and F450W$-$F814W were computed for the stellar remnants and were found to be $0.58$ and $1.38$, respectively. These values are consistent with the NSC colours observed by Georgiev & Böker (2014) and the expectations for a $\sim 10$ Gyr-old population.

5 Discussion & Summary

The origin of NSCs and their relationship to SMBHs have long been an interesting question in galaxy formation studies. NSCs are frequently present in many types of galaxies, and possess a significant older stellar population, indicating that a large fraction of the cluster has been in place for several Gyrs. This paper argues that this old stellar population in NSCs may be the remnant of a NSD that was present in the galaxy at $z\sim 1$. NSDs may be a common feature around growing SMBHs during an AGN phase at high-$z$ due to the large gas fractions within galaxies at this time. Previous work has shown that they provide a natural origin for the absorbing gas that obscures most AGNs at $z\sim 1$ (Gohil & Ballantyne, 2017, 2018). Importantly, NSDs require the presence of a SMBH, but the size and mass of the remnant cluster depend only weakly on the SMBH mass (e.g., Fig. 6). Therefore, the NSD theory provides an explanation for the complex observed relationship between NSCs and a central SMBH (e.g., Georgiev et al., 2016; Nguyen et al., 2018). However, as there are (typically, low-mass) galaxies with a NSC and no SMBH (e.g., NGC 205; Nguyen et al. 2018) the NSD mechanism proposed here will not be relevant for all galaxies.

We presented several quantitative comparisons in support of a NSD origin for the old stellar population of NSCs in galaxies with SMBHs. First, the predicted masses, colours and luminosities of the remnant stellar cluster correspond nearly exactly to those observed in NSCs (Figs. 4 and 5). This result only relied on the calculated equilibrium SFRs in each NSD and the dynamical time of the discs. Therefore, these predictions are a natural outcome of the star-forming properties of the disc. If NSDs are common in galaxies at $z\sim 1$ than remnant clusters will likely exist with these properties. Secondly, the predicted $R\textsubscript{\text{1/2}}$ of the remnant clusters follow a similar size-luminosity relations as observed in local NSCs (Fig. 5). Finally, the $M\textsubscript{\text{bh}}/M_{\ast}$ ratio of the remnants is also consistent with observations of local NSCs (Fig. 6).

Only $\approx 20$% of the remnants had $R\textsubscript{\text{1/2}}$ values consistent with the majority of local NSCs (Fig. 4), with the remainder being too large by about an order of magnitude, although this fraction depends on the sampling of the NSD parameter space that was explored (for example, models with larger $M\textsubscript{\text{bh}}$ would result in more remnants with $R\textsubscript{\text{1/2}}<10$ pc). This result suggests that NSDs are a viable origin for NSCs only when there is a very powerful pc-scale starburst in the NSD at $z\sim 1$. This scenario is predicted to be most common in galaxies with a large SMBH mass (Fig. 3). Therefore, NSDs may be a viable origin for NSCs in early-type galaxies (as supported by the large $M\textsubscript{\text{bh}}/M_{\ast}$ ratios) and in late-type galaxies with large $M\textsubscript{\text{bh}}$ masses (e.g., M31). While dynamical evolution of the stellar distribution in NSD remnants will likely not significantly impact $R\textsubscript{\text{1/2}}$ (Sect. 3.3), VRR will cause the final remnant to be both rotating and roughly spherical, consistent with the observed shape of the old populations of NSCs in some galaxies (Seth et al., 2008b; Carson et al., 2015).

In summary, it appears that the NSDs with strong pc-scale starbursts at $z\sim 1$ could be the origin of the old populations of local NSCs around massive SMBHs. Indeed, the stellar remnants predicted by this type of NSDs have sizes, colours and luminosities similar to many nearby NSCs. The stellar remnants are also expected to be rotating and have evolved into a near spherical distribution. Testing the predictions presented here will require both new observational and theoretical work. A clear prediction of our NSD model is that the most massive NSCs will be associated with the most massive SMBHs. Based on this result, ETGs may provide the best test cases for the NSD origin hypothesis, as these galaxies appear to consistently have the large $M\textsubscript{\text{bh}}/M_{\ast}$ ratios predicted by the theory. Furthermore, the lower gas content in the nuclei of ETGs will allow for easier charcterization of the old NSC population. Increasing the number of accurate NSC and SMBH mass measurements in both early- and late-type galaxies will also be important in testing the theory. However, the most compelling test of the model will be to search for evidence of NSDs in galaxies at $z\sim 1$, either through a combination of multi-wavelength measurements (Ballantyne, 2008), or by direct imaging using instruments carried by the James Webb Space Telescope. Once the incidence of NSDs is known at high-$z$, then the origin of many nearby NSCs may become much clearer. Lastly, it remains likely that most nearby NSCs will have undergone several episodes of growth, possibly driven by different mechanisms. Therefore, careful observational decomposition of the star clusters (e.g., Kacharov et al., 2018) will continue to be an important tool in understanding the origin of NSCs.

Acknowledgments

The authors thank A. Seth and D. Nguyen for providing the NSC observational data used in this paper, and the anonymous referee for helpful comments that improved the manuscript.