Investigating black hole accretion disks as potential polluter sources for the formation of enriched stars in globular clusters

Accretion disks surrounding stellar mass black holes (BHs) have been suggested as potential locations for the nucleosynthesis of light elements, which are our primary observational discriminant of multiple stellar populations within globular clusters. The population of enriched stars in globular clusters are enhanced in N14, Na23, and sometimes in Al27 and/or in K39. In this study, our aim is to investigate the feasibility of initiating nucleosynthesis for these four elements in BH accretion disks, considering various internal parameters such as the temperature of the gas and timescale of the accretion. To achieve this, we employed a 132-species reaction network. We used the slim disk model, suitable for the Super-Eddington mass accretion rate and for geometrically and optically thick disks. We explored the conditions related to the mass, mass accretion rate, viscosity, and radius of the BH-accretion disk system that would allow for the creation of N14, Na23, Al27, and K39 before the gas is accreted onto the central object. Our findings reveal that there is no region in the parameter space where the formation of Na23 can occur and only a very limited region where the formation of N14, Al27, and K39 is plausible. Specifically, this occurs for BHs with masses lower than 10 solar masses, with a preference toward even lower mass values and extremely low viscosity parameters ($\alpha<10^{-3}$). Such values are highly unlikely based on current observations of stellar mass BHs. However, such low mass BHs could actually exist in the early universe, as so-called primordial BHs. In conclusion, our study suggests that the nucleosynthesis within BH accretion disks of four elements of interest for the multiple stellar populations is improbable, but not impossible, using the slim disk model.


Introduction
The origin of multiple stellar populations (MSPs) in globular clusters (GCs) is still an unresolved puzzle.These populations appear as separate sequences in color-magnitude diagrams and they also show very specific patterns in the distribution of stellar abundances, the three most studied ones being the anti-correlation between carbon (C) and nitrogen (N), oxygen (O) and sodium (Na), and magnesium (Mg) and aluminium (Al).Two main stellar populations have been identified: so-called "pristine stars," with pristine chemical abundances, and "enriched stars," showing enhancement in N, Na, and Al, and depletion in C, O, and Mg.
At first, Norris & Pilachowski (1985) found that giant stars with a high nitrogen abundance are also enhanced in sodium.Denisenkov & Denisenkova (1989, 1990) mentioned evolutionary mixing to explain these chemical peculiarities.Sodium could be synthesized through the 22 Ne(p,γ) 23 Na reaction in red giants, at the same location inside the star, where N is created and O destroyed, thus resulting in an overabundance of both N and Na and a depletion of O. Later on, Langer & Hoffman (1995) suggested that Al could be considerably enhanced through very deep mixing and the activation of the MgAl reaction chain within bright giants.However, observations (Gratton et al. 2004) have provided evidence against this type of evolutionary model as main ⋆ e-mail: laurane.freour@univie.ac.at sequence stars were found to also exhibit abundance variations (Cannon et al. 1998).
Many other theories have been suggested to shed light on the origin of these abundance patterns (see Bastian & Lardo 2018 for an extensive review).Most of them involve a mixing between "pristine" and "enriched" material coming from different polluters.Already 40 years ago AGB stars were found to show abundance anomalies in some chemical elements (Cottrell & Da Costa 1981).Some 20 years later, Ventura et al. (2001) proposed that low-mass stars might have been polluted by the envelope of AGB stars, leading to the observed abundances; p-capture elements can indeed be synthesized during the hot-bottom burning phase in AGBs.However, Denissenkov et al. (2015) highlighted that the temperatures reached during this phase are too high (T ≳ 100 MK) to reproduce the observed isotopic ratios of Mg, requiring a narrow range of temperature around 80 MK, as shown by Prantzos et al. (2007).Denissenkov et al. (2015) suggested that supermassive stars (Denissenkov & Hartwick 2014) would be suitable polluter candidates as these stars can reach the needed central temperature, contrary to main sequence stars with a mass of M ≲ 10 3 M ⊙ (such as fast-rotating massive stars, Decressin et al. 2007).Still, such supermassive stars currently remain hypothetical and have not been observed yet.De Mink et al. (2009) suggested interacting binaries as a polluter source, where the primary star could eject most of its envelope in the surrounding medium due to non-conservative mass transfer.Using interacting binaries and fast-rotating massive stars, Bastian et al. (2013) proposed the early disk accretion scenario.Lowmass protostars would sweep up the enriched material ejected by these two polluters through their protoplanetary disk.However, hydrodynamical simulations carried out by Wijnen et al. (2016) revealed an incompatibility between the different timescales involved.Protoplanetary disks would evaporate much faster than the timescale required for the creation of a sufficient amount of enriched material.None of the polluter sources proposed so far have been able to reproduce all the observational constraints on MSPs at the same time.
In particular, most of them are facing the so-called mass budget problem.Observations have revealed a high percentage of enriched stars, ranging from 40% in low-mass clusters to 90% in high-mass clusters (Milone et al. 2016).The amount of enriched material required to form the enriched stars largely exceeds the amount of material that can be created by the polluters.To solve this issue, a modified IMF has been invoked (Prantzos & Charbonnel 2006;Charbonnel et al. 2014) as well as a heavy loss of polluters from the GCs (Renzini 2013), but both solutions are controversial (Bastian 2015).Breen (2018) suggested that light elements could be created in accretion disks around black holes if the right conditions in terms of temperature and density of the gas were reached.These chemical elements could then be expelled through outflows from the accretion disk, mix with the surrounding pristine gas, and form the enriched stars in globular clusters.A sufficient amount of enriched material could be created in only ∼ 3 Myr, depending on the initial number of black holes, assuming a super-Eddington accretion flow, overcoming the mass-budget problem.The observed abundance trends in MSPs are consistent with H-burning through the activation of the CNO cycle as well as the NeNa and MgAl chains (Langer et al. 1993).The temperature involved in such reactions is a key parameter to obtain the expected yields.
To reproduce the light element abundance patterns, the polluter should have a narrow range of temperature between 70 MK and 80 MK, (Prantzos et al. 2007).Heavier elements such as K require higher temperatures, namely, T > 180 MK, to be created, which can often be incompatible with the proposed stellar polluter sources (Iliadis et al. 2016;Prantzos et al. 2017).Black hole accretion disks can have more versatile properties than stars, depending on the optical depth, the viscosity, and the state of the black hole, resulting in broad ranges of densities, temperature, and radial velocities (see Abramowicz & Fragile 2013 for a review of black hole accretion disks).A temperature profile covering a broad range could be an advantage for the creation of the enriched material.If the accretion disk has a radial temperature profile between 200 MK and 50 MK, it would be possible to synthesize not only Na, O, C, N, Mg, and Al, but also heavier elements like K. The diverse range of physical properties exhibited by accretion disks may account for the variety of properties among stellar populations in different globular clusters.
Several studies have set their primary focus on investigating heavy element nucleosynthesis within black hole accretion disks.Chakrabarti et al. (1987) explored the conditions under which nucleosynthesis could occur inside thick accretion disks around black holes.At the center of the thick disk, products of the PP-chain, the CNO cycle, and rp-process (rapid proton process) could be created.They found that the amount of the elements produced was strongly dependent on the physical parameters of the disk, such as the mass of the black hole and the viscosity parameter, with smaller masses boosting the nucleosynthesis of such elements.This is promising as such reactions produce elements lighter than Z=40.Using the formalism presented in Chakrabarti & Mukhopadhyay (1999), Mukhopadhyay & Chakrabarti (2000) investigated three different types of accretion flow: hot, moderately hot and cooler flows.Depending on the inputs of their accretion disk model (seven variables in total), they find that the nucleosynthesis of some chemical elements is possible.Hu & Peng (2008) investigated the nucleosynthesis in advection-dominated accretion flow (ADAF) and outflow regions around a 10M ⊙ black hole.They found that for a heating factor f ≈ 0.001, representing the fraction of the viscously dissipated energy, a significant amount of 26 Al, 51 Cr, 53 Mn, and 55 Fe could be created in the accretion disk.These results suggest that nucleosynthesis could occur in accretion disks around black holes and be an effective way to increase the abundance of certain chemical elements in the interstellar medium.
This work focuses on the issue of multiple stellar populations in globular clusters.As a starting point, our attention is only directed towards specific (light) elements whose abundance pattern is observed in most of the GCs: C, N, O, Na, Al, and Mg.We also considered the case of K, challenging current theories as it requires very high temperatures to be enhanced, which are often unreachable by the classical polluter sources.Several points should be considered when reckoning whether black hole accretion disks could be the pollution source for the enriched material in globular clusters, for instance, whether the conditions in the temperature can be reached in an accretion disk; what the different timescales involved in the formation and accretion of the material are; and whether they are compatible with observations.In this paper, we address the two first points.In Section 2, we introduce the reaction network that we employed and analyze the effects of temperature and time on the yields of C, N, O, Na, Al, Mg, and K in nucleosynthesis.Then, in Section 3, we discuss the accretion disk model used and examine alternative models.Our approach for determining the feasibility of enriched material formation in an accretion disk surrounding a stellar mass black hole and the corresponding results are presented and discussed in Section 4. We summarize our findings in Section 5.

Nucleosynthesis of the enriched material
The stars in most GCs show chemical anomalies in light elements (He, C, N, O, Na, Al, and sometimes Mg).In other GCs, the MSPs also show a spread in heavier elements such as silicon (Si) or K.In this study, we restricted the computation to C, N, O, Na, Al, Mg, Si, and K.The temperatures required to initiate the depletion of O and formation of N from the CNO cycle are ≥ 40 MK, matching the temperature activating the NeNa chain.The MgAl chain requires higher temperatures ≥ 70 MK to lead to an enhancement in Al and a depletion in Mg.Si and K are created for T > 100 MK and T > 180 MK, respectively.Thus, we varied the temperature between 50 MK and 250 MK.
The abundance of a given species over time is described by a differential equation that depends on its interactions with other chemical elements.If we consider two particles j and k, then the mean lifetime of particle j against destruction by particle k is (Hix & Meyer 2006): where ⟨σv⟩ j,k is the reaction rate between j and k, depending on the temperature, and n k is the number density of particle k.Thus, a higher density leads to a lower mean lifetime for the particle.In the reaction network simulations, this translates into unchanged final abundance when increasing (decreasing) the density by a certain factor and simultaneously decreasing (increasing) the simulation time by the same factor, as highlighted by Prantzos et al. (2007).Running two simulations with densities ρ 1 and ρ 2 = 100 × ρ 1 for a time t 1 and t 2 = 0.01t 1 , respectively, leads to similar results.Considering the coupling between density and time, we maintain a constant value for ρ and investigate the evolution of the mass ratio of chemical elements with respect to time.The main drivers in our simulations are thus the temperature and the timescale, 50 < T < 250 MK and 10 −6 < t < 10 14 s ≈ 3 Myr.The upper limit for the timescale corresponds to the approximate time needed for a cluster to be gasfree according to observations of young massive clusters (Bastian et al. 2014).
To run the nuclear reaction network simulations, we used the python package Pynucastro (Smith Clark et al. 2022), an opensource project providing a user-friendly interface to solve nuclear reaction network equations.For a two-body reaction involving species A and B, the differential equation describing the evolution of the molar fraction, Y A , of species A is given by: where the term (1+δ AB ) in the numerator accounts for double destruction in case species A and B are identical and N A ⟨σv⟩ is the reaction rate.The reaction network is characterized by a set of differential equations similar to Eq. 2. Due to the stiffness of the system, specific numerical methods should be employed to solve such differential equations.We used the backward differentiation formula (BDF) integration method (Byrne & Hindmarsh 1975) from the SciPy module.Most reaction rates are taken from the REACLIB rate library (Cyburt et al. 2010).In our reaction network, we included all the species producing or consuming the chemical elements involved in the CNO cycle, the MgAl, NeNa, and K reaction chains.In such way we are sure that no important reaction is omitted from the network.A total of 132 isotopes are included.
In Sect.2.1, we introduce the observed chemical abundances necessary to initialize our simulations and to compare our results with observations.As we are specifically interested in enhanced species (N, Na, Al, K), we define the enrichment timescale for the enhancement of these chemical elements in Sect.2.2.Such an enrichment timescale is particularly relevant when studying the formation of enriched material in accretion disks or in gaseous environments that are not stars, because the timescales associated with stellar objects are determined by stellar evolution models, and the nucleosynthesis in stars is different.Finally, we investigate the conditions in temperature and timescale required to create the abundance patterns observed in GCs in Sect.2.3.The work presented in this section is independent of the polluter source (AGBs, FRMS, black hole's accretion disks, ...) and of the formation scenario for multiple stellar populations.

Observed pristine abundance pattern
The first step to set up our simulations is to define an initial set of abundances (given as mass ratios).We used the pristine composition of NGC 6752 (see Table 1) because this cluster has been widely studied over the past few years (Carretta 2013;Lapenna et al. 2016;Mucciarelli et al. 2017;Pancino et al. 2017;Lee 2018;Martins et al. 2021;Mucciarelli et al. 2019).It is one of the few globular clusters for which observations of C, N, Na, O, Al, Mg (including its isotopes), Si, and K are available (Yong et al. 2003;Carretta et al. 2005Carretta et al. , 2007Carretta et al. , 2012)).For these chemical elements, we estimated an average value for the pristine mass ratio.For helium, we use a mass fraction of 0.246 (Milone et al.Table 1.Pristine mass ratio used as input in the nuclear reaction code.The chemical elements are listed in the first column and their corresponding mass ratio in the second column.

Element
M  (Gratton et al. 2005).For the other elements (Ar, Ca, Sc, and Ni), we proceeded as follows.First, we took the solar mass fraction of chemical element X from Table 4 in Lodders (2021).Then, as done by Prantzos et al. (2007), following Goswami & Prantzos (2000), we recovered the abundance of element X at the metallicity of NGC 6752.These abundances are given as [X/Fe] and we apply the following formula to convert from abundance to mass ratio: To derive the abundances of isotopes for a chemical element, we further multiplied the right-hand side of Eq. 3 by the isotopic ratio.The adopted pristine mass ratios are presented in Table 1.

Enrichment timescale
In nuclear physics, a widely used timescale is the burning time, τ burn , characterizing the timescale on which an individual chemical element is destroyed (Hix & Meyer 2006): where Y(Z) is the nuclear abundance of specie Z and Ẏ(Z) the time variation of this abundance.The burning time is proportional to the inverse of the density ρ and the reaction rate ⟨σv⟩: (5) However, for our purpose we are particularly interested in the timescale required for a chemical element to be enhanced with respect to the initial (or pristine) value.Thus, we define the enrichment timescale, τ e .We use the parameter ξ to set a threshold for the definition of the timescale such that: The enrichment timescale of a chemical element, X, is determined by such factors as ξ, the initial mass ratio, M X 0 , and the temperature of the simulation, the latest playing a crucial role in obtaining the mass ratio, M X , in the nuclear reaction network.Figure 1 illustrates the definition of τ e .We plot the evolution of the mass ratio of 27 Al for two temperatures.Then, we set the enrichment threshold ξ = 2. Finally, we determine the enrichment timescale of 27 Al for each simulation by identifying the point in time when the 27 Al mass ratio surpasses twice the pristine mass ratio.If ξ = 2, as in the case of Fig. 1, the enrichment timescale is thus the time necessary to obtain a mass ratio for the chemical element X greater than two times the initial mass ratio, M X 0 .Similarly, we define a depletion timescale in Eq. 7, using the same parameter ξ.This timescale is used to inspect the conditions in temperature and time required to obtain the observed anti-correlations between chemical elements in MSPs: In order to provide some context for the values of ξ, we first need to introduce the notion of dilution.A widely accepted explanation for the emergence of MSPs in globular clusters involves the intervention of a polluter, which generates and expels enriched material containing a specific chemical element in a mass ratio denoted as M e .This enriched material subsequently mixes with the surrounding pristine gas, characterized by a mass ratio of M p , leading to the formation of an enriched star.The final composition of the resulting star is: where ω is the fraction of enriched material.In Fig. 2 we plot the observed abundances of seven chemical elements taken from Carretta et al. (2005) for C and N, from Carretta et al. (2007) for O and Na, from Gruyters et al. (2014) for Al and Mg, and from Mucciarelli et al. (2017) for K, and represented as black stars.Our aim is to illustrate the influence of ξ on the achieved enrichment after dilution to compare it with observations.Here, we focus on the impact of ξ and ω combined.The effect of ω alone has already been extensively discussed in previous works (e.g., Prantzos et al. 2007;Di Criscienzo et al. 2018).We used four distinct values for ω, going from ω = 0.1, where the final star is mainly made of pristine material to ω = 0.7, where the star is mainly made of enriched material.For each value of the mixing weight ω we vary the enrichment parameter ξ from ξ = 1 (the enriched material has the same composition as the pristine) to ξ = 17 (enrichment of the material seventeen times higher than the pristine composition).The final abundances obtained are plotted as colored circles.We used the same value of ξ for both chemical elements in the same panel.For example, in the upper left panel of Fig. 2, the orange point on the dilution trail ω = 0.7 corresponds to an enrichment threshold of ξ = 13 for Na and a depletion threshold of ξ = 13 for O.An important observation from Fig. 2 is the discrepancy in abundance ranges among different chemical elements, with ∆[O/Fe] being approximately 1, while ∆[Mg/Fe] is around 0.3.Consequently, while a value of ξ = 2 (represented by blue dots, indicating that the expelled material has a mass ratio twice its pristine value) can reproduce the extreme magnesium and potassium compositions with a mixing weight of ω =0.5, replicating the extreme oxygen abundance necessitates a significantly higher value of ξ even with a substantial mixing weight.However, it is important to keep in mind that the output of our simulation, which is the mass ratio, is strongly influenced by the temperature.Consequently, the enrichment timescale is also affected by the temperature.Figure 3 illustrates how the enrichment timescale changes with the parameter ξ for different temperatures, represented by the colored lines and ranging from 50 MK (blue line) to 250 MK (red line).When ξ = 8, the minimum value for the enrichment timescale is τ e = 10s, usually quite large compared to a typical viscous timescale (see Sect. 3).Such value of ξ only makes it possible to recover the extreme composition in Mg and Al, assuming a mixing weight w > 0.5 (see Fig. 2, pink dots).Overall, ξ should be smaller than 3 to give an enrichment timescale τ e < 1s, and the temperature should also have specific values (T > 150 MK).As shown on Fig. 2, only the extreme Mg abundance can be recovered with ξ = 3 (blue dots).The main conclusion drawn from the figure is that as the desired enrichment level increases, the corresponding timescale also increases for all temperatures.However, it is important to note that for certain temperatures (e.g., T > 200 MK), achieving a high enrichment level (e.g., ξ = 4) is not possible within the duration of the simulation ∆t = 10 14 s.

Parameter space for the anti-correlations
Using the initial chemical abundances given in Table 1, varying the temperature between 50 MK< T < 250 MK and the time 10 −6 < t < 10 14 s, we obtain for each temperature and time the final mass ratios for the 132 chemical elements.The density is kept constant to 1000 g /cm 3 .This value has been chosen arbitrarily and has only a small impact on the results, because of the coupling between time and density mentioned at the beginning of Sect. 2. If we would decrease the density, for instance to get a better agreement with accretion disks around compact objects in  The lines are colored according to the temperature (in MK).The higher is ξ, the more enriched is the final composition in the corresponding chemical element, here 27 Al.The enrichment timescale is directly dependent on the temperature.a binary system, the primary impact would be a rightward shift of the results shown in Figs. 1 and 4, toward larger timescales.
Using the definition of the enrichment and depletion timescales from Sect.2.2 with ξ = 2, we plot in Fig. 4 the region of parameter space in temperature and time required to obtain the specified anti-correlations, key characteristic of MSPs.The empty parameter space in the lower right panel stands out.This means that in the entire temperature range explored and for 10 14 order of magnitude in time, Potassium never be produced together with the observed light chemical elements anti-correlations.For these light elements, only a very narrow range of temperature allows the simultaneous production of N, Na, Al, and Si and destruction of C, O, and Mg, in good agreement with Prantzos et al. (2007).These results are valid regardless of the initial value for the density.As mentioned in Sect.2, changing the density will only change the timescale in which chemical elements can be formed, thus shifting the parameter space above to the left or to the right.These results, independent from polluter type and scenarios for the formation of MSPs, suggest a different origin for the heavy element patterns, pointing toward either multiple polluter types or a single polluter able to reach lower temperatures close to 70 MK as well as T > 150 MK in a short amount of time.Motivated by these results we explore the concept of nucleosynthesis occurring in accretion disks surrounding stellar mass black holes as potential sources of pollution.We delve into this topic in the subsequent two sections.

Accretion disk model
The Shakura-Sunyaev (S-S) disk model (Shakura & Sunyaev 1973) is one of the most commonly used analytical model to describe accretion disk.Using a few assumptions, (Shakura & Sunyaev 1973) showed that it is possible to solve analytically the system of thin disk equations.The S-S model is stable for sub-Eddington mass accretion rates ( ṁ << ṁEdd ).The disk is geometrically thin (and is therefore referred to as a "thin disk" model), and optically thick.The accretion disk is divided in three domains, depending on different physical conditions.In the inner part, the radiation pressure and the electron scattering are dominant.In the middle part, electron scattering is still the main source of the opacity, but gas pressure dominates.Finally, in the outer part of the disk, the temperature and density are low.The gas pressure still dominates but the main contributor to the opacity is free-free absorption.
The Novikov-Thorne (N-T) model (Novikov & Thorne 1973) is the relativistic generalization of the S-S model.The black hole is considered to be a rotating Kerr black hole.Relativistic correction terms are added to the Newtonian equations.These corrections are particularly important close to the inner edge of the rotating Kerr black hole where the gravitational pull is big.Far away from the center, N-T equations are equivalent to the S-S equations.In our case, we are mainly interested in the material close to the inner edge of the black hole, where outflows arise.
Breen (2018) emphasized the necessity of a super-Eddington accretion flow to generate a sufficient amount of enriched material within a limited timeframe.This requirement stems from the need for the material to be produced in a matter of a few million years, aligning with the observational constraints observed in gas-free young massive clusters.The slim disk model, origi-nally introduced by Abramowicz et al. (1988) serves as a valuable expansion of the S-S model in scenarios involving a super-Eddington accretion flow, as visually represented in Fig. 1 from Abramowicz & Fragile (2013).For accretion rates ṁ ≈ 0.1 ṁEdd , the S-S disk may become unstable due to the radiation pressure contribution to the viscous torque (Lightman & Eardley 1974;Shakura & Sunyaev 1976;Janiuk et al. 2002).For mass accretion rates close to or larger than the Eddington limit, the disk becomes thick and advection comes in the balance as an extra cooling mechanism.The slim disk can be used as an alternative disk model.
In an optically thick accretion flow and for high accretion rates (i.e., in the framework of the slim disk model), a portion of the photons produced within the disk becomes unable to exit its surface because they are confined by the inward flow of matter.This critical radius, at which this confinement happens, is commonly referred to as the "trapping radius," r trap .Such a radius is only dependent on the mass accretion rate.It delimits the accretion disk into two regions: the inner region (in the slim disk regime) and the outer region (where cooling is radiation-dominated), and the accretion disk can be described by the standard model.The existence of such radius is one of the main difference with the advection-dominated accretion flow (ADAF) model in which the disk is optically thin and photons can therefore escape the disk immediately.As highlighted by Narayan & Yi (1995), who introduced self-similar solutions for ADAF, the temperature profile is then almost virial.In such an accretion flow, the ion and electron temperature is decoupled at small radii.The temperature profile plotted in their Fig. 3 matches the temperature criteria for nucleosynthesis at large radius R > 5000 R sch .We suggest that ADAF cannot be the process generating the enriched material in GCs for two main reasons.Firstly, the density is significantly low, on the order of magnitude much smaller than that of S-S disks, as illustrated, for instance, in Hu & Peng (2008).Given the interdependence of time and density in two-body reactions within a chemical reaction network (see Sect. 2.3), the generation of enriched material would necessitate an extensive duration to take place in an environment characterized by such a low density.Secondly, the reaction chains involved in the formation of such material are activated for temperatures that are reached at very distant radius.As mentioned previously, outflows occur very close to the inner edge of the accretion disk.Thus, even if nucleosynthesis occurs in ADAF as suggested by Hu & Peng (2008) and Zhang et al. (2017), the synthesis of the chemical elements of interest for MSPs would happen at a large radius and could not be expelled into the surrounding medium to mix with the pristine gas.
Therefore, we chose to use the slim disk model, covering accretion rates ranging from 10 −4 ṀEdd to 10 5 ṀEdd .Although it has been acknowledged that the N-T model loses accuracy under super-Eddington accretion rates, we nonetheless included it in our simulations for verification purposes.We describe the slim accretion disk by using the detailed expressions presented by Watarai (2006).Aspects that are of particular importance to our discussion include: the temperature, density, and radial velocity of the flow, v r .The relevant equations used to describe the three quantities are characterized by two groups of equations, denoted by the subscript in and out, depending on the position within the accretion disk, with respect to a threshold radius.Such radius threshold is set by the trapping radius, defined in Watarai (2006) as the radius where the accretion time is shorter that the photon diffusion time, preventing photons from escaping due to the radial flow of matter.This radius marks the transition point between the standard disk, dominated by radiative cooling, and the slim disk, incorporating advection effects.The solutions within the trapping radius are as follows: where α is the viscosity of the gas, r g = 2GM/c 2 is the Schwarzschild radius, and m = M/M ⊙ and ṁ = Ṁ/(L E /c 2 ) are the normalized mass and mass accretion rate, expressed in unit of solar mass and Eddington accretion rate, respectively.These equations are slightly different from the canonical self-similar solutions (see, e.g., Wang & Zhou 1999), because of the treatment of energy transport.In Watarai (2006), the ratio of the advective cooling rate to the viscous heating rate, f, is not constant and given by: where D is a numerical coefficient taken to be equal 2.18, and x =(r/r g )/ ṁ.The quantities outside the trapping radius are similar to the solutions from Shakura & Sunyaev (1973): T out ≈ 4.79 × 10 7 m 10 It is important to highlight that the density is provided as (projected) surface density.The conventional conversion formula from 2D surface density to a 3D density profile is as follows: Here, H represents the disk's height.However, for the sake of simplicity and due to the coupling between time and density in nuclear reaction network simulations, we opted for the simplified vertical one-zone approximation commonly used in the standard disk model: Σ = 2Hρ.This choice is further discussed in Sect. 5. Finally, the viscous timescale can be defined as the time needed for the material at a given radius, R, to be accreted onto the compact object:

Nucleosynthesis in black hole accretion disks
In Sect.2, we showed the challenge to simultaneously obtain the anti-correlations between C and N, O and Na, Mg and Al, and Mg and K from the same temperature range.In this section, we aim to investigate whether the thermodynamic conditions of accretion disks are favorable for the production of N, Na, Al, and K.

Method
Table 2 summarizes the method followed in this work to obtain the quantity of interest represented by the viscous and enrichment timescales from the inputs of the accretion disk model.Firstly, we coded the accretion disks' equations.The outputs of the disk model are the density, temperature, and radial velocity of the accretion flow.The temperature is used as an input for the reaction network, as all the reaction rates are temperature dependent.By varying the parameters of the black hole and its accretion disk (black hole mass, m, mass accretion rate, ṁ, radial distance to the black hole, R, viscosity of the gas, α), we obtained different densities, temperatures, and radial velocities for the gas.Using these temperatures and densities as input of the reaction network, we can compute the enrichment timescale and compare it to the viscous timescale in order to determine whether accretion disks can host the required nucleosynthesis of different (light) chemical elements.The range of mass, mass accretion rate, radius, and viscosity investigated is given in the second column of Table 2.The range of black hole masses and accretion rates is very broad, covering 11 and 8 orders of magnitude, respectively.To ensure the completeness of our investigation, we sampled our parameter space using a two-step method designed to optimally cover it.First, we divided the space into a regular grid of range of [l × 10 i , l × 10 i+1 ] for i varying between 0 and log 10 (u/l) − 1 , where l and u are the lower and upper bound of the interval.Then, within each grid cell, we use latin hypercube sampling (McKay et al. 1979) to take 100 points randomly sampled.A Latin hypercube is a sampling method that ensures each variable in a set is sampled uniformly across its range, and each sample is independent of the others.We have a total of 11 × 8 × 2 × 4 different intervals within our four-dimensional (4D) parameter space.Within each interval, we randomly sampled 100 values for each parameter.We thus obtained a total of 704 × 100 × 4 = 281, 000 Table 2.This table provides an overview of the methodology employed to obtain the quantities of interest presented in the last column.The first three columns provide the input parameters, the ranges explored, and the outputs of the accretion disk model.The three next columns present the inputs, their range, and the outputs of the reaction network used.Two steps lead to the computation of the important timescales: the viscous timescale only depends on the accretion disk model, while the burning/enriched timescale depends on the evolution of the mass ratio of chemical element X (dependent on the temperature, thus, on the disk model and on the reaction network).
values of (m, ṁ,r,α).We used these values as inputs for the analytical model of the accretion disk and obtained 281,000 tuples of (T ,ρ,v r ).For the nuclear reaction network, we only used the tuple with T > 2 MK as no nucleosynthesis of N, Na, Al, and K can occur for smaller temperatures.We are left with a total of 25,730 values.For each value we run the nuclear reaction code, solving a stiff system of equations (Eq. 3 in Smith Clark et al. 2022).The results are presented in Sect.4.2.

Results
In the following, we study the ratio τ nuc = τ visc τ e .If τ visc < τ e , thus τ nuc < 1, the material in the accretion disk is accreted before nucleosynthesis can occur.If τ nuc > 1, the enrichment timescale is larger than the viscous timescale, and nucleosynthesis could occur in the accretion disk before the gas is accreted.The results are plotted in Figs. 5 and 6.
Figure 5 shows the entire range of parameter values leading to a temperature high enough for nucleosynthesis to occur (T > 2 MK), in the logarithmic plane, m versus ṁ, and colorscaled according to α.Some combinations of parameters, such as a high black hole mass and low accretion rate, result in temperatures that are too low to allow for nucleosynthesis.We ran the simulation for different enrichment thresholds (values of ξ).
In the left panel of Fig. 6, we show the values of (m, ṁ, r, α) resulting in an enrichment of 14 N, 27 Al, and 39 K greater than two times the pristine value (ξ = 2).Only a very light black holes with m < 1M ⊙ and very small viscosity parameter α < 0.003 makes the enrichment of 27 Al and 39 K possible.Regarding 14 N, a few models around m ≈ 1 M ⊙ and moderate mass accretion rate 1 Ṁ Edd < ṁ < 100 ṀEdd makes the enrichment possible.No combination of (m, ṁ, r, α) create the right conditions in temperature and density to enrich 23 Na with a threshold of ξ = 4.
We further investigated the potential formation of 14 N and 27 Al focusing on an important aspect that has not been addressed so far: the duration of the enrichment.In Fig. 1, the green-blue curve demonstrates a defined duration for the enrichment, which we estimate to be approximately ∆τ e ≈ 100 seconds.Conversely, the blue curve indicates an indefinite duration for the enrichment.
If ∆τ e is small, the chemical element is very quickly destroyed.Thus, even if it would be created, it would be very unlikely that it can be expelled from the accretion disk, pollute the surrounding medium, and explain the abundances in enriched stars in globular clusters.We define the duration of the enrichment: Here, ∆τ e can be visualized in Fig. 7, showing the evolution of 14 N and 27 Al mass ratios as a function of time for two different successful runs, as the time length when the different colored curves are above the enrichment thresholds represented by the green and red lines.We compute the duration of the enrichment ∆τ e for each of the points shown in Fig. 6.For 27 Al, ∆τ e barely exceeds 100 seconds and never exceeds 500 seconds in the 113 successful runs.Regarding 14 N, once the mass ratio is exceeding the threshold, it stays at a constant value until the end of the simulation.

Discussion
The results presented in Section 4.2 indicate the potential for nucleosynthesis of 14 N, 27 Al, and, interestingly, 39 K within an accretion disk surrounding sub-solar mass black holes using the slim disk model.These possibilities arise when there is a reasonable accretion rate, close to the Eddington value, and an extremely low viscosity parameter.In the following section, we delve into a more detailed discussion regarding the impact of enrichment duration on the ejection of enriched material and examine the feasibility of such values for m and α.
Given the small duration of the enrichment in 27 Al, less than a few hundred seconds, even if a tiny black hole with mass M = 0.01M ⊙ and an accretion disk with an extremely low viscous parameter (α < 10 −3 ) exists, the ejection of the newly formed 27 Al would have to occur in an incredibly fine-tuned way in order to pollute the surrounding medium.We find the same results after examining the duration of the enrichment for 39 K, even when setting a smaller value for the enrichment (ξ = 2, justified by the range of [K/Fe] abundances being much smaller than the range of [Al/Fe], as shown in Fig. 2).
The low value for the viscous parameter α is a very challenging constraint for the nucleosynthesis.While numerical simulations suggest a value close to α ≈ 0.02 (Hawley et al. 2011), observations are better represented with α ≈ 0.1 (King et al. 2007).In any case, the values required to enrich the gas in 14 N, 27 Al, and 39 K is almost two orders of magnitude smaller than what is suggested by numerical simulations.
The last parameter putting strong constraints on the nucleosynthesis is the black hole mass, which also tends to be very small, although the mass required for the nucleosynthesis of 14 N is considerably more plausible compared to the cases of 27 Al and 39 K.When considering a modeling approach based on X-ray observations, a black hole a mass of 1.3M ⊙ still falls within  (Carr & Kuhnel 2021) and evidence for their existence has recently been claimed using LIGO/Virgo gravitational-wave data (Franciolini et al. 2022).If such primordial black holes exist, they could be promising for globular clusters as they are among the oldest structures in the Universe.However, there are still many uncertainties around these objects such as: how long they survive across cosmic time; whether they can have an accretion disk and, if yes, what would be a suitable model to describe them; and whether they could have a much lower viscosity parameter, α, that is compatible with the value found in this work.Undoubtedly, the progress of research on PBHs will be significantly propelled by upcoming gravitational wave telescopes such as Einstein and LISA (Franciolini et al. 2023).

Conclusion
We have explored the possibility of accretion disks around stellar mass black holes as the polluter creating the enriched material at the origin of the phenomenon of multiple stellar population in globular clusters, following up on the study by Breen (2018).Using the Python package pynucastro (Smith Clark et al. 2022), we performed nuclear reaction simulations to investigate the condition in temperature and timescale required to create some of the anti-correlations between chemical elements observed in GCs.
We find that a very small range of temperature and time (equivalent to density) could lead to the simultaneous 12 C -14 N, 16 O -23 Na, and 24 Mg -27 Al anti-correlation patterns.However, the 24 Mg -39 K anti-correlations observed in some GCs requires a different region of the parameter space.These results are in very good agreement with previous works (Prantzos et al. 2007;Iliadis et al. 2016) and suggest that the enriched material is either coming from multiple polluters at different temperatures or from a single polluter having a large temperature range, for example black hole accretion disks.
In this study, we used the slim disk model with the framework outlined in the work by Watarai (2006).This optically and geometrically thick model proves to be particularly suitable for scenarios involving super-Eddington mass accretion rates.Here, advection emerges as a significant cooling mechanism and can no longer be neglected.Unlike the S-S model (Shakura & Sunyaev 1973), the slim disk model introduces less restrictive assumptions, notably by not imposing a Keplerian gas rotation profile.However, it is essential to acknowledge a potential limitation in our approach.To derive the volume density of the gas, an input for the nuclear reaction network, we employed a straightforward technique of dividing surface density by scale height.This method is typically employed under the one-zone approximation, which remains valid for geometrically thin accretion disks where the vertical variations of primary variables are averaged due to a significantly smaller scale height compared to the radial distance (H << r).However, this assumption does not hold true for the slim disk model.Nevertheless, given that the density is intricately linked with time in the nuclear reaction network and because the enrichment timescale generally exceeds the viscous timescale by several orders of magnitude, we anticipate that our results remain robust despite this assumption.Another limitation of our work lies in the absence of time dependence in the black hole accretion disk system.While using a constant value for the mass and viscosity is acceptable, the mass accretion rate can vary with time (see, e.g., van der Klis 2004).Depending on the timescale of the variation, this may or may not impact nucleosynthesis in the accretion disk.Conducting further hydrodynamical simulations in conjunction with a nuclear reaction net-work within accretion disks could offer additional insights and constraints regarding the potential for nucleosynthesis.
Using the slim disk model, we probed a very wide range of parameters for the black hole and accretion disk system, varying the black hole mass, m, the mass accretion rate, ṁ, the viscosity parameter, α, and the radius, r, over many orders of magnitude.Each combination of (m, ṁ, α, r) leads to a temperature, T , density, ρ, and radial velocity, v R , of the accretion flow.We then use T and ρ as input to the nuclear reaction network and obtain mass ratios for 132 chemical elements.We are particularly interested in chemical elements in over-abundance in the enriched stars of GCs, namely 14 N, 23 Na, 27 Al, and 39 K.
We defined the enrichment timescale as the time needed for a chemical element to be enriched with respect to the pristine value, using a free parameter ξ quantifying the level of enrichment.Subsequently, we compared the enrichment timescale to the viscous timescale, characterizing the lifetime of the gas before accretion onto the black hole, depending on the radius and the radial velocity of the flow.We generated 281,000 sets of values for the parameters (m, ṁ, α, r) and we extracted 25,730 (T , ρ) pairs where the temperature is high enough for the nucleosynthesis to happen (T > 2 × 10 6 K).Among theses pairs, only 82 lead to an enrichment timescale for 27 Al smaller than the viscous timescale, using an enrichment threshold of ξ = 2, meaning that 27 Al can be twice more enriched than the pristine mass ratio before the gas is accreted by the black hole in only ≈ 0.3% of the simulated points.The number of points reduces to 20 for 14 N, to 6 for 39 K, and to 0 for 23 Na.The few quadruplets making an enrichment in 27 Al or 39 K possible all correspond to very light black holes m < 0.01M ⊙ and very low viscosity parameter α < 10 −3 .Even under the assumption that stars consist entirely of enriched material (ω = 1), it is important to note that such a low enrichment threshold (ξ = 2) will still be unable to recover extreme [Al/Fe] values.The enrichment in 14 N requires a larger black hole mass, around 1 M ⊙ , but (as with 27 Al or 39 K) a very low viscosity parameter.Such a low value (α < 10 −3 ) is disfavored by current observations and hydrodynamical simulations.The question of whether this kind of black hole could really have existed in the past is yet to be answered.

Fig. 1 .
Fig. 1.Illustration of the concept of enrichment timescale.The green and blue lines represent the 27 Al mass ratio obtained from nuclear reaction simulations at temperatures of T = 80 MK and T = 150 MK, respectively.The dash-dotted line represents the threshold set to determine enrichment, with ξ = 2.In this case, 27 Al is considered enriched if M27 Al > 2 × M27 Al 0 .The two orange dashed lines correspond to the enrichment timescales, represented by their respective x-values.

Fig. 2 .
Fig. 2. Chemical abundance of O versus Na (upper-left), Mg versus Al (upper-right), O versus K (lower-left), and C versus N (lower-right).Black stars represent the observations taken from Carretta et al. (2005) for C and N, from Carretta et al. (2007) for O and Na, from Gruyters et al. (2014) for Al and Mg, and from Mucciarelli et al. (2017) for K.The green rectangle indicates the location of the most extreme abundances.The colored circles representing the abundances obtained using the reaction network (T = 80 MK and ρ = 1000 g/cm 3 ) and the dilution model introduced in Eq. 8 with ω = 0.1, 0.3, 0.5, and 0.7 are colored according to the enrichment level parameterized by the parameter ξ.While the extreme potassium abundances require only a small ξ, other chemical elements like oxygen require a large mixing weight and enrichment level to reproduce the extreme population.

Fig. 3 .
Fig.3.Illustration of the evolution of the enrichment timescale, τ e , as a function of the parameter ξ.The lines are colored according to the temperature (in MK).The higher is ξ, the more enriched is the final composition in the corresponding chemical element, here 27 Al.The enrichment timescale is directly dependent on the temperature.

Fig. 4 .
Fig. 4. Temperature and time generating the observed anti-correlations 12 C-14 N (upper-left), 16 O-23 Na (upper-right), 24 Mg-27 Al (middle-left), and 24 Mg-39 K (middle-right).The lower middle and right panels represent the intersection of the parameter space reproducing the observed abundance patterns in light elements (middle) and in all the elements studied (right).Only a small region in temperature and time makes it possible to reproduce the light element patterns.At such temperature and time, 39 K can not be enhanced, so a second set of temperature and time values are needed.

Fig. 5 .
Fig. 5. Parameter space exploration for the black hole and accretion disk system.A total of 25,730 points are plotted out of the initial 281,000 parameter, representing approximately 9% of the samples (m, ṁ, r, α) that yield temperatures exceeding the threshold of 2 MK required as input into the nuclear reaction network.The points are displayed in a logarithmic scale on the m versus ṁ plane and color-scaled based on α values.

Fig. 6 .
Fig.6.Identification of (m, ṁ, r, α) values resulting in an enrichment of 14 N, 27 Al, and 39 K larger than two (ξ = 2, left panel), and greater than four (ξ = 4, right panel).The size of the points corresponds to the distance to the center of the black hole, r.Only 20 (upper-left), 82 (middle-left), and 6 (lower-left) out of the 25,730 points enable an enrichment by a factor of 2, translating into 0.07%, 0.3%, and 0.02% cases of enrichment.When ξ = 4, the number of points becomes 19 (upper-left), 6 (middle-left), and 4 (lower-left), translating into 0.07%, 0.02%, and 0.015% case of enrichment.The largest distance to the black hole making an enrichment possible occurs for a value of r = 61r sch , visible in the upper and middle left panels.

Fig. 7 .
Fig.7.Variation of27 Al and 14 N mass ratio over time for two runs where the enrichment timescale exceeds the viscous timescale.The label details the parameters of the black hole and accretion disk system.The mass, M, is expressed in M ⊙ , the mass accretion rate in terms of the Eddington rate, Ṁ Edd , and the radius in unit of the Schwarzschild radius.The two thresholds for considering different enrichment level are represented by the dash-dot and the dotted lines.While the 27 Al mass ratio only remain above the lines for a brief duration, the enrichment in 14 N is rather stable and remains at a high level until the end of the simulation.