© The Authors 2025

1 Introduction

Ultraluminous X-ray sources (ULXs) have been a subject of intense study in the field of high-energy astrophysics for several decades. These enigmatic objects, characterized by X-ray luminosities exceeding 10³⁹ erg/s, challenge our understanding of accretion physics and compact object formation (see e.g., Kaaret et al. 2017; King et al. 2023, for a recent review). With over 1800 known ULXs identified to date (Walton et al. 2022), they represent a significant population of extreme X-ray emitters in the Universe, yet their fundamental nature remains a topic of ongoing debate.

The diversity of ULXs has been highlighted in recent classification schemes (e.g., Wiktorowicz et al. 2017), which attempt to categorize these objects based on their observational properties and potential underlying physical mechanisms. While the majority of ULXs were initially thought to harbor intermediate-mass black holes (IMBHs), our understanding has evolved considerably with the detection of pulsar accretors (Bachetti et al. 2014).

The importance of studying ULXs, particularly in the context of globular clusters (GC), cannot be overstated. ULXs serve as laboratories for extreme accretion physics, potentially shedding light on the formation and evolution of double compact objects – the progenitors of gravitational wave sources. Moreover, ULXs continue to be considered as candidates for elusive IMBHs with masses in the range of 10²–10⁵ M_⊙ (Kaaret et al. 2017, and references therein). The existence of IMBHs may provide crucial insights into the formation pathways of supermassive BHs observed at the centers of most galaxies. These intermediate-mass objects could also explain the presence of supermassive BHs observed at high redshift (Greene et al. 2020; Askar et al. 2023).

In this paper, we present a comprehensive study of ULXs in GCs, focusing on their stellar-mass manifestations. The dense stellar environment introduces additional complexities and possibilities for exotic binary formations. By examining their properties, distribution, and potential formation mechanisms, we aim to contribute to the broader understanding of ULX physics and their role in compact object evolution.

Dense, quasi-spherical collections of stars known as GCs orbit the centers of galaxies. These typically old stellar systems, often containing millions of stars within a relatively small volume, are hubs of dynamical interactions. The masses of Milky Way GCs range from a few 10³ to a few 10⁶ M_⊙ , with half-mass radii typically between 1.5 pc and a few tens of pc, and ages spanning from 8 to 13.5 Gyr (Harris 1996; Salaris & Weiss 2002; VandenBerg et al. 2013; Baumgardt & Hilker 2018). Detailed properties of specific Galactic GCs can be found in the catalogs of Harris (1996, updated 2010)¹ and Baumgardt & Hilker (2018)².

The evolution of GCs is governed by the relaxation process connected with continuous distant gravitational interactions between all objects in the cluster. This leads to global phenomena like core collapse, that is contraction of the cluster core and gradual increase in central density, and mass segregation, that is the tendency of heavier stars to occupy the central regions while pushing lighter stars outwards. Close dynamical interactions also influence stars and binaries individually, affecting binary parameters (e.g., separation and eccentricity), or through exchanges (e.g., binary changes one of it’s component into another star), ejections (star or binary becomes unbound with the cluster), and binary formation (single stars became bound and form a binary).

Such dynamical encounters in dense star clusters can lead to the formation of close binary systems, including X-ray binaries (Pooley et al. 2003; Ivanova et al. 2005, 2008; Kremer et al. 2018, 2019), cataclysmic variables (Oh et al. 2024; Belloni et al. 2019), and gravitational wave progenitors (e.g., Benacquista & Downing 2013; Mandel & Broekgaarden 2022, and references therein). Recent studies have emphasized the importance of GC dynamics in shaping exotic stellar populations (Di Carlo et al. 2021; Oh et al. 2024), highlighting the role of these environments in producing and evolving compact object binaries.

While the majority of ULXs are detected in star-forming regions of galaxies, there is increasing evidence for their presence in GCs as well. These GCULXs are found exclusively in extragalactic GCs, all of which are unresolved at distances of >16 Mpc, with no ULXs detected in Milky Way GCs to date. Pioneering observations of GCULXs have been reported in several galaxies, with varied GC environments, spectral behavior, and variability patterns.

Maccarone et al. (2007) discovered the first ULX in a GC (NGC 4472). It exhibits an X-ray luminosity of 4 × 10³⁹ erg/s and rapid variability. Similar detections were reported in NGC 1399 (Shih et al. 2010; Irwin et al. 2010), with more GCULXs in NGC 4472 (Maccarone et al. 2011) and NGC 4649 (Roberts et al. 2012). Some sources identified in GCs exhibited flares just above the ULX-defining limit of 10³⁹ erg/s (Irwin et al. 2016; Sivakoff et al. 2005).

Recent observations have extended the known GCULX population, including seven additional GCULXs associated with M87 GCs (Dage et al. 2020), three in NGC 1316 (Dage et al. 2021), 10 new candidates in GCs of massive (>10^11.5 M_⊙) early- type galaxies (Thygesen et al. 2023), and two more in NGC 4261 (Nair et al. 2023). The current population of approximately 30 known GCULXs, identified through systematic searches of thousands of GC systems in elliptical galaxies, represents a remarkably rare phenomenon compared to the ~1800 ULXs identified in other galactic environments (Walton et al. 2022).

In dense stellar systems, compact objects can be readily ejected, which can considerably impact the formation channels of such objects (e.g., Kulkarni et al. 1993). Furthermore, GCs typically exhibit high ages and, consequently, fewer highermass stars. Given their presence in older stellar populations, GCULXs are likely powered by low-mass X-ray binaries, in contrast to field ULXs which are typically associated with massive donor stars but often lack optical counterparts. Given their presence in older stellar populations, GCULXs are likely powered by low-mass X-ray binaries, in contrast to field ULXs which are expected to be also associated with low-mass as well as massive donor stars, but often lack optical counterparts. These characteristics contribute to the intriguing nature of GCULXs.

2 Methodology

2.1 MOCCA Monte Carlo code

The MOCCA (MOnte Carlo Cluster evolution Code; Hypki & Giersz 2013) is a state-of-the-art numerical tool that combines the Monte Carlo method for stellar dynamics (Hénon 1971; Stodolkiewicz 1986) with detailed stellar and binary evolution algorithms from the Binary Star Evolution code (BSE; Hurley et al. 2000, 2002) with recent updates (level-C from Kamlah et al. 2022, and references therein) and low-N scattering code FEWBODY (Fregeau et al. 2004) to follow close dynamical interactions. This hybrid approach makes MOCCA particularly well-suited for studying the formation and evolution of exotic objects like ULXs in dense stellar environments, offering an optimal balance between computational efficiency and physical accuracy.

Recent improvements to the MOCCA code (Hypki et al. 2022, 2025; Giersz et al. 2024) have significantly enhanced its capabilities for ULX studies. First of all, the mass transfer calculations now incorporate a more sophisticated treatment of super-Eddington accretion, crucial for modeling ULX systems (Shakura & Sunyaev 1973; Wiktorowicz et al. 2015). Additionally, improvements in core radii calculations improved the Roche-lobe overflow calculations. The introduction of detailed evolutionary tracking provide unprecedented insight into the spatial distribution and evolutionary pathways of ULX systems within clusters equivalent to detailed outputs from codes like startrack (Belczynski et al. 2008).

A major advancement relevant to this study is MOCCA’s capability to handle multiple stellar populations, a phenomenon now recognized as ubiquitous in GCs (e.g., Piotto et al. 2015; Renzini et al. 2015; Gratton et al. 2019; Bastian & Lardo 2018; Milone & Marino 2022). While a comprehensive analysis of multiple stellar populations in our models is presented in Hypki et al. (2022, 2025); Giersz et al. (2024), here we specifically investigate their impact on ULX formation and evolution. The code now tracks different stellar populations with distinct chemical compositions, ages, and spatial distributions, allowing us to examine how properties specific to multiple stellar populations affect compact object formation and ULX characteristics. For detailed analysis of simulational results on multiple stellar populations and all recent changes to the MOCCA code see Giersz et al. (2024, and references therein).

The X-ray luminosity (L_X) of accretion disks is calculated using the Shakura-Sunyaev model (Shakura & Sunyaev 1973), accounting for both sub-Eddington and super-Eddington accretion regimes. The Eddington mass accretion rate is defined as ${\stackrel{.}{\text{M}}}_{\text{Edd}}=2.2\times {10}^{-8}{\text{M}}_{\text{acc}}{\eta }_{0.1}^{-1}\left[{\text{M}}_{\odot }{\text{yr}}^{-1}\right]$, where M_acc is the accre- tor mass and η(=10η_0.1) is the accretion efficiency. We assumed η_0.1 = 1. The Eddington luminosity is then L_Edd ≈ 5·66 × 10⁴⁵ Ṁ_Edd [erg s⁻¹]. The X-ray luminosity is calculated as: $$L_X=\{\begin{array}{ll}{\text{L}}_{\text{Edd}}\left(\stackrel{.}{\text{M}}/{\stackrel{.}{\text{M}}}_{\text{Edd}}\right),\hfill & \text{ for }\stackrel{.}{\text{M}}\le {\stackrel{.}{\text{M}}}_{\text{Edd}}\hfill \\ {\text{L}}_{\text{Edd}}\left[1+\ln \left(\stackrel{.}{\text{M}}/{\stackrel{.}{\text{M}}}_{\text{Edd}}\right)\right],\hfill & \text{ for }\stackrel{.}{\text{M}}>{\stackrel{.}{\text{M}}}_{\text{Edd}}\hfill \end{array}$$(1)

where Ṁ is the mass accretion rate. See Lasota & King (2023) for recent discussion of this formulation.

2.2 Simulations

Figure 1 illustrates three pathways for ULX formation in relation to GCs. In the “IN-CLUSTER” scenario, ULX progenitors form and remain bound to the cluster, often through multiple dynamical interactions, including exchanges. “ESCAPER” ULXs originate in the GC but are ejected through interactions, emitting outside the cluster while retaining their GC-imprinted dynamical history. “FIELD” ULXs, by contrast, form and evolve in isolation. The key distinction lies in the dynamical imprint: “ESCAPER” ULXs bear the hallmarks of GC interactions, unlike “FIELD” systems, which provide a baseline for understanding the unique contribution of clusters to ULX formation.

The MOCCA code employs hundreds of parameters to govern the dynamical and stellar evolution of GCs. Based on previous studies of compact object populations – including cataclysmic variables (Belloni et al. 2017), ULXs in the field (Wiktorowicz et al. 2019), and double white dwarf systems (Hellström et al. 2024) – we identified key parameters that most significantly influence ULX formation and evolution. Our parameter selection strategy focused on those that control dynamics and initial stellar population properties. This systematic approach allows us to explore the most relevant parameter space while maintaining computational feasibility. In this section we present discussion of the utilized parameters.

In all simulations we adopted the Kroupa (2001) initial mass function (IMF) for both stellar populations, with mass ranges of 0.08-150 M_⊙ for the first population and 0.08–20 M_⊙ for the second population. For binary systems, we implemented a pairing mechanism, that combines a uniform mass ratio distribution (0.1 < q < 1.0) for massive stars (M > 5 M_⊙) following Kiminki & Kobulnicky (2012); Sana et al. (2012); Kobulnicky et al. (2014), with random pairing for lower-mass stars. All clusters were initialized in virial equilibrium (Q_vir = 0.5). For the underlying density distribution, we used King models with concentration parameters of W₀ = 3.0 and W₀ = 7.0 for the first and second populations, respectively, representing moderately concentrated initial configurations. More details on the MOCCA parameters can be found in Kamlah et al. (2022); Hypki et al. (2022, 2025); Giersz et al. (2024).

For all simulations with dynamics enabled, we tracked the evolution of escaped systems to analyze the properties of their populations and potential ULXs that can form after their progenitors are ejected from the cluster. These escapers evolve as isolated binaries, similar to field binary evolution (e.g., Fragos et al. 2015; Wiktorowicz et al. 2019; Zuo et al. 2021). Systems can escape the cluster through multiple mechanisms: (1) dynamical interactions, particularly strong few-body encounters, (2) natal kicks from supernovae explosions, including both direct and Blaauw kicks, (3) gradual two-body relaxation (Spitzer & Hart 1971), and (4) tidal stripping of outer cluster regions by the host galaxy’s gravitational field (Gnedin & Ostriker 1997).

The galactocentric distance (R_gc) for each simulation was scaled to maintain an initial tidal radius of R_tid ≈ 43 pc across all models. This normalization of the tidal radius facilitates meaningful comparisons between simulations, particularly when analyzing the degree of tidal filling and cluster structural parameters. The tidal radius follows the relation ${\text{R}}_{\text{tid}}\propto {\text{R}}_{\text{gc}}^{2/3}{\text{M}}_{\text{tot}}^{1/3}$ (e.g., Webb et al. 2013), where M_tot is the cluster mass. This scaling ensures that clusters experience comparable relative tidal forces despite different masses and orbital parameters. We note that while the initial tidal radii are identical, the clusters’ subsequent evolution may lead to different filling factors³ depending on their internal dynamics and mass-loss history.

We performed a comprehensive set of simulations exploring different initial conditions and physical parameters. Table A.1 summarizes our simulation grid, where each model is identified by a unique label used throughout this paper. The simulations vary in the number of stellar populations (N_pop), with varied spatial distribution of different generations through the concentration parameter (R_h,2/R_h,1), defined as the ratio between the half-mass radii of subsequent generation relative to the first generation. The initial binary fraction (f_bin) was either 10%, or 95%, which affected also other parameters (see below). We explored both tidally filling (TF) and non-tidally filling (nTF) configurations (Tidal). Some simulations incorporate the new features (NF) described below, particularly relevant for multiple population scenarios. For each simulation, we provide key structural parameters at t = 0 Myr, including the galactocen- tric radius (R_gc), half-mass radius (R_h), total stellar mass (M_tot), core radius (Rc; according to Casertano & Hut 1985), core mass (M_c, i.e., mass inside R_c), and central density (ρ_c; according to Casertano & Hut 1985). Additionally, we computed corresponding field populations with dynamics turned off (denoted by nodyn in the model labels) to serve as control cases, though these are not shown separately in the table.

To better compare simulations with different binary fractions, we made several adjustments to the parameter space. Since a lower binary fraction results in lower total stellar mass (as binary systems typically have higher masses than single stars), we adjusted the initial number of objects to have a similar halfmass relaxation time for both models: simulations with f_bin = 95% used n = 600 000 objects, while those with f_bin = 10% used n = 1 063 635 objects. We also modified the initial distribution of semi-major axes. For simulations with f_bin = 95%, we used the modified version of the Kroupa (1995) period distribution (see Belloni et al. 2017) for stars with M < 5 M_⊙, and the distributions from Sana et al. (2012); Oh et al. (2015) for stars with M > 5 M_⊙. For simulations with f_bin = 10%, we employed a log-uniform distribution for M < 5 M_⊙ and the Sana et al. (2012) period distribution for M > 5 M_⊙. Additionally, for the low binary fraction simulations, we disabled eigenevolution (Belloni et al. 2017). The separation was limited to the maximal value of 100 AU.

In a part of our simulations (marked with NF in label), we introduced NF to investigate the influence of multiple stellar populations (see Giersz et al. 2024, for details). Specifically, we implemented a delay time (t_delay) of 100 Myr for the second population, during which these stars act as gas particles without undergoing stellar evolution, relaxation, or dynamical interactions. The gas accumulation for the second population begins at t_{delay_fraction} = 0.5, meaning halfway through the delay period, coinciding with the asymptotic giant branch (AGB) wind contribution phase. To accommodate this, we activated AGB wind accumulation for the second population formation.

The cluster’s orbital evolution was modified by implementing an artificial orbit change at t = 1000 Myr, where the galactocentric distance increases by a factor of 2 (see Giersz et al. 2024, for details). Finally, we employed single-population scaling, where the scaling parameters are applied exclusively to the first population, while the second population’s velocities are independently normalized. First and second populations are separately in the virial equilibrium.

Except dynamical evolution parameters, our simulations incorporate several key prescriptions for stars and binaries, some of them particularly relevant to ULX formation, specifically: we employed the rapid supernova mechanism as described by Fryer et al. (2012) for both neutron star (NS) and black hole (BH) formation; Neutron star and BH natal kicks were drawn from a Maxwellian distribution with σ = 265 km/s (Hobbs et al. 2005); For BHs, these kicks were scaled by the fallback factor, which depends on the supernova prescription; Neutron stars formed through electron-capture supernovae were assigned substantially lower kicks with σ = 3 km/s; Common envelope evolution was modeled using the parameters from the BSE code: α = 0.5 and λ = 0.0. All simulated models in this study were set to a fixed metallicity of Z = 0.001, which represents the typical metallicity of Galactic globular clusters⁴.

Fig. 1 Schematic representation of the origins and eventual locations of ULX progenitors and ULXs in GCs. Progenitors can either form and remain bound to the cluster, leading to in-cluster ULX emission (“INCLUSTER”), or be dynamically ejected and emit as ULXs outside the cluster (“ESCAPER”). Additionally, ULXs originating and emitting in the field (“FIELD”) represent progenitors formed without direct interaction with cluster dynamics.

3 Results

ULXs in this study are defined as binary systems with neutron star or black hole accretors exhibiting X-ray luminosities exceeding 10³⁹ erg s⁻¹. White dwarf-powered ULXs are excluded from this analysis.

3.1 Number evolution of ULXs

Figure 2 illustrates the temporal evolution of the expected number of ULXs⁵ (E[N_ULX]) across all simulations, accompanied by the count of simulations yielding non-zero E[N_ULX] for each time bin. The highest E[N_ULX] is observed in the initial evolutionary stages, typically within the first 300 Myr (first bin in Figure 2). This trend persists across all simulations, irrespective of parameter variations. These findings align with previous studies on field populations, where dynamical interactions are negligible, which reported that the ULX phase predominantly occurs in the early evolutionary stages of binaries (e.g., Wiktorowicz et al. 2017). Observational evidence also supports this conclusion (e.g., Wolter et al. 2018).

This result suggests either:

• Early ULX evolution is largely independent of cluster dynamics parameters, or

• Insufficient time has elapsed for dynamics to significantly influence ULX formation.

The number of simulations with E[N_ULX] , 0 decreases over time, resembling predictions for populations formed in burstlike star formation episodes (e.g., Wiktorowicz et al. 2017). In GCs, this trend may result from a combination of system age effects on ULX formation (as observed in field populations) and the interplay of positive and negative effects related to dynamics. Generally, our results indicate a higher likelihood of observing ULXs in younger stellar clusters (≲300 Myr) compared to older ones, which is consistent with observational results (e.g., Dage et al. 2025).

Figure 2 reveals substantial variation in E[N_ULX] between simulations, particularly in later evolutionary phases. This variability likely stems from stochastic processes (simulation independent, that is resulting from random number generator seed) and varying dynamical effects (simulation dependent). The high variability in later evolutionary phases in comparison to the first bin (0–300 Myr) where low variation is observed suggest that, E[N_ULX] is affected by dynamics in at least a significant way.

The variations in the otherwise monotonically decreasing number of non-zero predictions can be attributed to randomness. The peak at 5 Gyr and even more pronounced peak at 12 Gyr are not correlated to any significant changes in the GC structure and are not statistically significant, therefore can be treated as anomalies.

Please note that the provided E[N_ULX] values represent predictions per cluster. To derive observational expectations, one must sum these values across all observed clusters and account for observational limitations, such as magnitude limits and resolution constraints. A detailed discussion of observational predictions will be addressed in a separate study.

Fig. 2 Evolution of ULX numbers in GCs. Expected number of ULXs (E[NULX]) as a function of time since GC formation. All simulations with non-zero expected rates are presented. The line indicates the number of simulations with non-zero predictions (i.e., number of dots in this time bin). Each time bin spans 300 Myr and the points are centered on these bins.

Fig. 3 Comprehensive overview of the E[NULX] evolution since the formation of a GC. Each panel represents a set of simulations with identical parameters: number of populations (Npop), concentration (Rh,2/Rh,1), binary fraction (fbin), and the presence of new features (NF). Within each panel, ULX populations from various environments are presented: tidally filling (TF) GC, non-tidally filling (nTF) GC, nodyn (binaries evolving without dynamical interactions), and escapers (progenitors formed in a GC, but experiencing the ULX phase after ejection from the cluster). See Table A.1 for simulations parameters overview. Bars and points represent the E[NULX] in 300 Myr bins.

3.2 Temporal evolution and parameter dependence of E[N_ULX]

Figure 3 presents a comprehensive analysis of E[N_ULX] dependence on model parameters and its temporal evolution. The most significant effects are observed in parameters related to the environment, specifically whether the GC is tidally filling (TF) or non-tidally filling (nTF), and the general impact of dynamics (incluster vs. field formation). These results support the hypothesis that dynamical interactions strongly influence the formation of ULXs in dense stellar environments. In particular, the exchange process may lead to the formation of binaries that are not possible through regular formation channels. The specific trends and patterns are described below.

Non-tidally filling clusters exhibit significantly larger ULX populations compared to tidally filling clusters. This discrepancy can be attributed to the extended time frame density available for dynamical interactions to enhance progenitor formation in non- tidally filling clusters, whereas tidally filling clusters experience more early ejections. More generally, the expected number of ULXs formed in GCs is substantially higher than in corresponding field populations, particularly for non-tidally filling clusters. Typically, ULXs formed in GCs demonstrate higher E[N_ULX] rates in specific time bins, and these rates are more continuous, whereas field populations mostly exhibit isolated occurrences throughout the timeline. These relations do not apply to early evolutionary stages where all simulations show high and similar values of E[N_ULX] (see Section 3.1).

Simulations with two populations (middle column) generally produce more ULXs than those with a single population (left column). This can be attributed to the higher concentration of the second population, which enhances ULX progenitor formation. However, when the concentration is excessively high (right column), this effect is less pronounced, which supports the claim that too frequent and strong dynamical interactions can substantially destroy ULX progenitors.

The initial binary fraction appears to have a limited influence on ULX formation in GCs. For f_bin = 95% wide binaries are quickly destroyed in dynamical interactions and mostly remain binaries similar to ones for f_bin = 10%. The population of ULX progenitors remains sparse and emerges regardless of the initial relative number of binaries. Dynamical effects effectively destroy existing binaries and create new ones, a process primarily responsible for the formation of ULX progenitors at later evolutionary times. Conversely, the binary fraction influences field populations, with a positive correlation between initial binaries and ULX progenitors, as there are no external processes to destroy ULX progenitors or create new ones in isolation.

The incorporation of NF has a noticeable effect on two-population simulations with moderate concentration (R_h,2/R_h,1 = 0.2; Figure 3, middle column). Simulations with NF exhibit higher E[N_ULX] values compared to those without NF. The NF do not affect field populations as they primarily relate to the time delay for the second population.

The escaper populations generally align with the ULX population formed in situ, suggesting ongoing ejection (relaxation and interactions), and thus a continuous supply of ULX progenitors to the escaper population, which later become ULXs while unbound to the cluster. There are generally fewer ULXs among escapers than in in-cluster populations. However, ULXs from escapers can dominate field populations. In general, there are no significant differences in the formation time of the ULXs from escapers to those in the field. An exception occurs in simulations with relatively strong ULX formation (N_pop = 2, R_h,₂/R_h,1 = 0.2), where ULXs from escapers form at evolutionary ages not covered by field populations. For tidally filling clusters, the rate of ULX formation from escapers can exceed that of in-cluster formation, especially in simulations without NF.

3.3 Formation of ULXs

This section focuses on ULX progenitors – systems that will evolve into ULXs – and their evolution up to the initial ULX phase. It is important to note that some systems may experience multiple ULX phases interspersed with periods of quiescence or X-ray binary phases (L_X < 10³⁹ erg/s). These phase transitions typically occur on timescales of at least hundreds of years, driven by dynamical interactions or stellar evolution processes. Additionally, individual stars can participate in multiple ULXs, each with a different companion. This phenomenon is particularly prevalent for IMBHs (see Sect. 3.5). While our models focus on these long-term changes, it’s worth noting that observations have shown both globular cluster and field ULXs can exhibit significant short-term variability on timescales of days to months (e.g., Maccarone et al. 2007; Shih et al. 2010; Brightman et al. 2023).

Table 1 provides a comprehensive overview of ULX progenitors on zero-age main sequence (ZAMS), that is stars that later became components of ULXs. The most common progenitors are those of pristine ULXs, that is binaries that do not undergo exchanges or disruptions prior to the ULX phase, although still they can dynamically interact and change their properties. In these systems, the primary typically becomes the accretor and the secondary the donor (A-D(p)). These progenitors evolve rapidly, initiating the ULX phase after approximately med(t_ULX,start) ≈ 6 Myr, in stark contrast to other progenitors which generally begin after med(t_ULX,start) ≳ 3 Gyr. Pristine ULXs are observed in all simulations.

Interestingly, some pristine ULXs form where the secondary, less massive on ZAMS, becomes the accretor (D-A(p)). These progenitors are considerably rarer (about two orders of magnitude less common than A-D(p)) and typically form ULXs only in later evolutionary phases (∼7 Gyr).

Non-pristine ULXs can originate from both binary and single star progenitors on ZAMS. However, stars born in binaries more frequently become components of ULXs, even in non-pristine cases. These systems invariably undergo significant interactions, such as binary formation, exchanges or disruptions, before forming a ULX. While some of these ULXs can form relatively early (∼100 Myr), the majority emerge in later evolutionary phases (>3 Gyr). In non-pristine ULXs, both the primary and secondary can act as either the accretor or donor with comparable probability. Our simulations did not reveal any non-pristine ULXs where both the primary and secondary are parts of ULXs (i.e., “A-D” or “D-A” type), though this is likely due to limited statistics rather than a physical constraint.

3.4 ULX properties

In this section we discuss the properties of in-cluster ULXs, that is these which reside inside the cluster. For a discussion of ULXs among escapers see Section 3.8.

Table 2 summarizes the key properties of ULXs in our sample, categorized into different groups based on the types of accretors and donors. The majority of ULXs in our sample are powered by BH accreting from various companion types. BH accretors paired with MS donors (G_BH,MS) show median masses of 27.64 M_⊙, while those with HG and CHeB companions (G_BH,HG and G_BH,CHeB, respectively) have slightly lower median masses (22.70 M_⊙ and 16.88 M_⊙, respectively). The maximum X-ray luminosities for these systems (median log₁₀L_X,max ≈ 40.18–40.36 erg/s) are well above the ULX defining limit of 10³⁹ erg/s. The observed sample of GCULXs, with the brightest having L_X ≈ 4 × 10³⁹ erg/s (Dage et al. 2019), appears inconsistent with such high luminosities. However, it is important to note that the simulated ULXs exhibiting these extreme luminosities primarily form in very young clusters (t_phys ≲ 10 Myr), whereas the observed GCULXs are found in older globular clusters. Our simulations predict <0.1 E[N_ULX] in older clusters (8–13 Gyr) similar to Milky Way globular clusters, and comparable limit for intermediate-age clusters (2–8 Gyr). Additionally, high-luminosity ULXs are typically short-lived, further reducing their detection probability.

Within ULXs, NS accretors (G_NS) have median masses around 1.27 M_⊙. These systems typically have lower X-ray luminosities (log₁₀ L_X,max ≈ 39.22) compared to their BH counterparts, but still clearly exceed the Eddington limit for a typical NS.

The orbital characteristics of ULXs vary significantly across different groups. GBH,MS systems have median separations of 40.05 R_⊙, while G_NS systems show much tighter orbits (median 0.45 R_⊙). Most ULXs in our sample have low eccentricities, suggesting that tidal forces or a common envelope phase have circularized their orbits before the onset of the ULX phase.

Duty cycles vary widely among ULX groups, with G_BH,MS systems showing the highest median value (0.41), while other groups have much lower values. This variability in duty cycles may explain the transient nature observed in many ULXs and has implications for their detectability (cf. Wiktorowicz et al. 2017).

A notable subset of our sample consists of IMBH-powered systems (G_IMBH). These ULXs are characterized by extremely high accretor masses (median 9351 M_⊙) and the highest X-ray luminosities in our sample (log₁₀L_X,max ≈ 40.84). Such luminous systems, called extreme-ULXs (Wiktorowicz et al. 2019), are one of the main observational candidates for IMBHs and can provide insights into their formation and properties.

The diversity in donor types across ULX groups highlights the various evolutionary pathways that can lead to ULX formation. Donors of types MS, HG, and RG are all well-represented, suggesting that ULXs can form and persist across a wide range of stellar evolutionary stages.

Huge majority are ULXs with BH accretors (typically >80%) formed early in cluster evolution.

Table A.2 shows a modest anticorrelation between the relative size of the G_IMBH and BH ULXs (Kendall⁶'s τ = −0.49, −0.38, −0.37, for G_BH,MS, G_BH,HG, and G_BH,CHeB, respectively). The presence of IMBH results in a continuous formation of IMBH powered ULX though the cluster history which, being an independent source to other BH ULX groups present mostly in the early cluster evolution phases, lowers their fractional input to the ULX population. IMBH quickly removes BHs and possibly NSs from the system (e.g., Hong et al. 2020).

3.5 IMBH ULXs

In this study, we define IMBHs as those with masses exceeding 500 M_⊙. IMBH ULXs were observed in 6 simulations, all featuring multiple stellar populations (npop2) and non-tidally filling initial conditions (nTF). When present, IMBHs contribute significantly to the ULX population, accounting for 16–44% of ULXs over a Hubble time (Table A.2).

IMBHs form rapidly, within several Myr after the ZAMS, classifying them as “fast IMBHs” (Giersz et al. 2015). In simulations with delayed second populations, IMBH formation occurs 100+ Myr after the initial population. Notably, IMBH ULXs form exclusively in non-tidally filling (nTF) simulations, strenghtening the preference for denser stellar environments.

Fisher’s exact tests revealed no significant correlations between simulation parameters and IMBH formation at the 90% confidence level. However, when considering only simulations with multiple populations (pop=2) and active dynamics, a statistically significant relationship (95% confidence) emerged between the filling parameter and IMBH ULX presence, favoring non-tidally filling clusters.

Table 3 presents key properties of the six IMBHs observed in our simulations. These IMBHs display a substantial mass range, from as low as ~500 M_⊙ to as high as ~13500 M_⊙, demonstrating their significant growth potential within globular cluster environments. IMBH ULXs typically exhibit intermittent ULX phases (duty-cycle between 0.02 and 0.87%) with various donors, interspersed with periods of non-ULX mass accretion or mergers. Donors are mostly MS stars, with white dwarfs (particularly COWD) being the second most common companion type.

IMBH ULXs formed exclusively in simulations with two stellar populations and non-tidally filling initial conditions. The centrally concentrated second population facilitates stellar interactions and mergers, while nTF clusters retain more stars, particularly BHs, in central regions, extending the IMBH formation window.

IMBH ULXs typically exist throughout the entire GC life, with first ULX phases appearing as early as ~8 Myr in some models and continuing until the end of our simulations at ~14 Gyr, indicating that IMBH ULXs can form very late in GC evolution.

3.6 ULX progenitors

The initial properties of binary systems that evolve into ULXs provide crucial insights into their formation channels and evolutionary pathways. Our results indicate that the primary distinguishing factor between ULX progenitors and the general stellar population is their higher initial mass. No other significant constraints on ZAMS properties were identified for ULX progenitors.

Table 4 summarizes the ZAMS properties of ULX progenitors, categorized by ULX group and component role (accretor or donor). Accretor progenitors typically have higher ZAMS masses (M_ZAMS ) than their companions (M_comp), consistent with being predominantly formed from primary stars (f_prim ≥ 0.89).

The majority of ULX progenitors originate in binary systems (f_bin ≥ 0.67), with most maintaining their original pairing to form the ULX (f_pris ≥ 0.66, except for the G_IMBH group). For accretors in the G_IMBH and G_NS groups, f_bin ≈ 0.66, indicating no preference between single and binary star origins.

Initial separations (a) exhibit a highly skewed distribution, with G_NS group progenitors typically having much wider separations than other groups. Values for the G_IMBH group are not particularly informative as these stars do not become components of ULXs without undergoing strong dynamical interactions. All groups show a preference for lower initial eccentricities (med(e) ≲ 0.3).

The G_IMBH group shows distinct characteristics compared to other groups. It has the lowest ZAMS masses for both accretors and donors, a unique dominant formation modes (A-n for accre- tors and D-n for donors), and the lowest f_bin . The G_BH,MS and G_BH,HG groups exhibit similar patterns, with high ZAMS masses and a strong preference for the A-D(p) formation mode. The G_NS group has the highest fraction of ULXs formed from escapers (f_esc = 0.39).

3.7 ULX descendants

The fate of ULX systems after their active phase provides valuable insights into their evolutionary pathways and potential contributions to other astrophysical phenomena.

Table 5 summarizes the various outcomes for ULX systems across different simulation configurations. The table shows, that ULX systems can undergo several fates: disruption, escape from the cluster, exchange interactions, mergers, or other outcomes. The distribution of these fates appears largely stochastic across different simulation parameters. Systems categorized as “other” typically experience no significant events for the remainder of their evolution and are predominantly NS/BH + WD or double compact objects at the end of simulation.

Merger events typically occur ${0.79}_{-0.53}^{+0.74}$ Myr after the ULX phase. Double compact objects form in ∼63% of cases, with ∼34% of them being ejected from the cluster before merger or the end of simulation. These double compact objects that remain in the cluster may contribute to the population of gravitational wave sources with merger rate of ∼11%.

Disruptions may result from dynamical interactions or binary evolution. In exchanges, at least one of the components continues evolution as part of some other binary. Our results suggest, that in initially more centrally concentrated environment (nTF) there are more disruptions and exchanges after the ULX phases, then in other simulations.

3.8 ULXs among escapers

Cluster-formed binaries are frequently ejected, suggesting that some field ULXs may originate from GCs. These ULX progenitors become unbound from their parent clusters and initiate the ULX phase while residing in the galactic field, evolving as isolated binaries.

Figure 3 illustrates the expected numbers and ages of ULXs among escapers. Generally, escaper ULXs coexist in similar numbers and at comparable evolutionary times as their in-cluster counterparts, except the very early phase of cluster evolution (≲ 300 Myr). Table 6 presents a detailed comparison of properties between in-cluster and escaper ULXs.

Escaper ULXs exhibit several distinct characteristics compared to their in-cluster counterparts. They tend to have lower- mass accretors and donors, with median masses of 13.00 and 6.05 M_⊙, respectively, compared to 22.24 and 39.34 M_⊙ for incluster ULXs. Escaper ULXs also have a significantly higher fraction of NS accretors (40% vs. 4%) and tighter orbits (median 57.42 vs. 119.93 R_⊙). This orbital difference largely stems from dynamical hardening processes that occur prior to ejection, whereby close encounters cause binaries to become more tightly bound (Heggie 1975). Additionally, the ejection mechanisms preferentially remove harder binaries with sufficient kinetic energy to escape the cluster potential. The median maximum X-ray luminosity of escaper ULXs (log₁₀L_X,max = 39.89 erg/s) is slightly lower than that of in-cluster ULXs (log₁₀L_X,max = 40.26 erg/s).

Notably, escaper ULXs are less numerous, with a median count of 31 compared to 176 for in-cluster ULXs. They also typically have more evolved companion stars, with 65% having CHeB donors compared to 55% HG donors for in-cluster ULXs. The duty cycles and ULX phase onset times are similar for both populations.

The ejection processes preferentially remove less massive systems, explaining the lower masses observed in escaper ULXs. This mass-dependent retention occurs because more massive systems require stronger dynamical encounters to achieve escape velocity from the cluster potential. The higher fraction of NS accretors among escapers could be due to the higher natal kicks and retention of more massive BHs within the cluster potential.

These findings have important implications for understanding the origin and evolution of field ULXs. A significant fraction of observed field ULXs may have originated in GCs, which could explain some of the observed diversity in ULX populations. While our current analysis captures the properties of escaper ULXs at the moment of escape and during the ULX phase, tracking their complete galactic trajectories would require additional post-processing to integrate their orbits through the galactic potential (e.g., Cabrera & Rodriguez 2023). Furthermore, the differences between escaper and in-cluster ULXs provide valuable insights into how environment and dynamical history can influence the properties of these extreme systems.

In conclusion, our simulations reveal that escaper ULXs form a significant and distinct population with properties that differ from their in-cluster counterparts. These differences provide valuable insights into the formation and evolution of ULXs in various environments and highlight the need for comprehensive observational surveys to fully characterize the ULX population across different galactic settings.

4 Discussion

The results presented in this study are specific to the simulations employed and may not capture general trends applicable to all GC environments. The parameter choices were carefully tailored to conditions characteristic of GCs, informed by insights from analogous studies.

A key limitation in the broader study of ULXs in GCs is the dearth of confirmed observations within the Milky Way and its immediate vicinity. The brightest Milky Way X-ray binaries fall well below the ULX luminosity threshold (Fortin et al. 2023, 2024). The substantial distances involved often complicate the detection and characterization of low-luminosity donor stars, which are believed to comprise a significant fraction of ULX donors (e.g., Wiktorowicz et al. 2017).

Furthermore, a critical challenge for this study lies in the difficulty of unambiguously differentiating GC-associated ULXs from those in the field environment. Observational uncertainties frequently obscure the distinction, as localization on the plane of a GC does not inherently confirm membership within the cluster. Determining the precise distance to the source is often fraught with uncertainties, leaving room for the possibility that the ULX is a foreground or background source, merely coincident with the GC in projection. This ambiguity highlights the importance of complementary approaches, such as proper motion studies or radial velocity measurements, to confirm GC membership conclusively.

4.1 Comparison with field populations

We find a high fraction of NS accretors among escapers (40%) compared to in-cluster ULXs (4%). This is similar to the predictions of Wiktorowicz et al. (2019), who found that NS ULXs outnumber BH ULXs in regions with constant star formation and solar metallicity for ages above ∼ 1 Gyr.

The properties of our escaper NS ULXs show some similarities with typical field NS ULXs described by Wiktorowicz et al. (2017). They found that field NS ULXs typically have ~1.3 M_⊙ NS accretors and ~1.0 M_⊙ red giant donors. This is comparable to our escaper ULXs with NS accretors, which have low-mass ≲2 M_⊙ evolved (HeMS) donors. Our systems are very compact (a ≲15 R_⊙) which may explain the lower fraction of giant donors.

Our escaper NS ULXs have a median maximum X-ray luminosity of log₁₀L_X,max = 39.22 erg/s, with a few BH ULXs exceeding 10⁴⁰ erg/s. This stays in contrast to the field populations, where a significant fraction of more luminous sources is also expected (e.g., Wiktorowicz et al. 2015). In the case of GCs, we find that such extreme ULXs are more likely to be retained within the cluster.

The formation pathways of escaper ULXs may differ from those in the field. For instance, Fragos et al. (2015) studied the formation of NS ULXs like M82 X-2, finding that systems with 8–10 M_⊙ donor stars could produce such ULXs. In contrast, our escaper NS ULXs have lower donor masses (≲2 M_⊙), suggesting a different formation pathway.

The temporal distribution of ULXs in our simulations differs from that observed in field populations. Kuranov et al. (2021) found that the maximum number of ULXs (∼10 for a star formation rate of 10 M_⊙ /yr) is reached ~1 Gyr after the beginning of star formation in field populations. Our escaper ULXs show a different temporal distribution with a more uniform appearance of escaper ULXs through the GC lifespan.

It is worth noting that even field populations are not entirely free from dynamical influences. Klencki et al. (2017) showed that wide binaries in the field are affected by multiple weak interactions (fly-bys) which can impact binary evolution. This suggests that the distinction between cluster-formed and field- formed ULXs may not be as clear-cut as previously thought, and that dynamical effects could play a role in shaping ULX populations across various environments.

4.2 Counterparts

Identifying ULX counterparts is crucial yet challenging (Heida et al. 2016; Wiktorowicz et al. 2021). While field ULXs face issues of faintness and crowded environments (López et al. 2020; Allak et al. 2022; Heida et al. 2019), extragalactic GC ULXs present unique challenges due to their unresolved stellar populations. However, GCs’ well-understood stellar demographics suggest these ULXs are likely LMXBs. Rare cases allow indirect donor identification, such as oxygen emission indicating a WD donor in a NGC 4472 GC ULX (Steele et al. 2014). Accurate astrometry remains essential for all ULX studies (Allak et al. 2022), and most counterparts, in both field and GCs, remain tentative due to observational biases (Heida et al. 2019; Wiktorowicz et al. 2021).

4.3 Expectations for Antennae galaxy

The correlation between X-ray source positions and stellar clusters has been previously noted in starburst galaxies (e.g., Kaaret et al. 2004). Poutanen et al. (2013)’s study of the Antennae galaxies reveals a significant association between bright X-ray sources and stellar clusters, with most X-ray sources located outside clusters. This supports the idea that many ULXs are massive X-ray binaries ejected from their birth clusters, rather than IMBHs.

The displacements of ULXs from cluster centers, ranging up to 300 parsecs, are statistically significant and likely due to ejection of massive binaries. The ejection mechanisms suggested by Poutanen et al. (2013), including supernova kicks and few-body encounters, are consistent with our results showing higher NS accretor fractions and smaller separations among escaped ULXs.

The distribution of ULXs likely evolves with stellar population age. While the Antennae study focus on young populations, older stellar systems are expected to decrease their X-ray binary populations (e.g., Fragos et al. 2008), highlighting the importance of considering population age in ULX distribution studies.

4.4 Future prospects

Future studies of ULXs in GCs should include more simulations to better estimate statistical effects, errors, and the significance of conclusions. A more consistent survey of the parameter space is needed, as well as the incorporation of evolutionary parameters (both stellar and binary evolution) beyond dynamical ones.

In the future work we plan to include the effect of beaming in ULXs (e.g., Lasota & King 2023). Beaming increases the apparent luminosity of ULXs, but decreases the probability of observation due to possible misalignment (Wiktorowicz et al. 2019; Khan et al. 2022). Population synthesis simulations, Wiktorowicz et al. (e.g., 2019) indicate that the majority of NS ULXs are beamed. The beaming factor is dependent on the mass transfer rate, with higher mass transfer rates resulting in stronger beaming (King 2009).

Additionally, the importance of wind Roche-lobe overflow (see e.g., Wiktorowicz et al. 2021; Zuo et al. 2021, and references therein) for NS ULXs with (super)giant donors should be explored. Wind Roche-lobe overflow can lead to stable mass transfer even with large mass ratios, and can boost mass transfer rates to reach ULX luminosity levels.

5 Conclusions

This study presents the first numerical investigation into the formation and evolution of ULXs within GCs, utilizing a subset of simulations to explore parameter dependencies.

Our simulations reveal that dynamical interactions may play a critical role in ULX formation, particularly in environments with high stellar density. Even if the initial binary fraction is low, the interactions can produce many ULXs. On average, we find that approximately 96% of ULXs in our simulations have BH companions and the number can be even higher for very young clusters (≲300 Myr). Among escaper, ULXs have a much higher fraction of NS accretors (∼40%).

Our simulations show that the ratio of escaper ULXs to incluster ULXs is approximately 1:7, but nearly 2:1 for ULX with NS accretors. This ratio is the highest in tidally filling clusters, even with high initial binary fractions, where rapid ejections can hinder in-cluster ULX formation by ejecting the progenitors.

Our findings suggest that the relative scarcity of ULXs observed in GCs may be attributed to the advanced age of their stellar populations. In contrast, field populations usually have continuous star formation.

The apparent absence of ULXs in Milky Way GCs aligns more closely with our models of initially tidally filling clusters rather than non-tidally filling ones. This distinction is crucial: tidally filling clusters experience significant early mass loss through tidal stripping by their host galaxy’s gravitational field, leading to the ejection of potential ULX progenitors before they can evolve into ULXs. In contrast, non-tidally filling clusters maintain higher central densities for most of their life, allowing for enhanced ULX formation through dynamical interactions. The presence of IMBH ULXs exclusively in non-tidally filling simulations further emphasizes the importance of maintaining dense stellar environments for certain ULX formation channels.

Furthermore, our results suggest that field populations may be significantly “polluted” by ULXs ejected from GCs (escapers), which could contribute to the relative underrepresentation of ULXs in GCs compared to field environments. Escapers can have properties very similar to the field populations or form distinct configurations.

Future observational efforts should prioritize the comparison of ULX properties across diverse environments, including GCs and the galactic field. Such studies will be critical for validating our theoretical predictions and advancing our understanding of the diverse formation channels of ULXs.

Acknowledgements

We are grateful to anonymous referee for their insightful comments which helped to improve the paper. G.W., M.G., A.H., L.H. were supported by the Polish National Science Center (NCN) through the grant 2021/41/B/ST9/01191 A.A. acknowledges support for this paper from project No. 2021/43/P/ST9/03167 co-funded by the Polish National Science Center (NCN) and the European Union Framework Programme for Research and Innovation Horizon 2020 under the Marie Skłodowska-Curie grant agreement No. 945339. For the purpose of Open Access, the authors have applied for a CC- BY public copyright license to any Author Accepted Manuscript (AAM) version arising from this submission.

Data availability

The data underlying this article are available in Zenodo at https://doi.org/10.5281/zenodo.14953837. This includes simulation data on ULX systems, globular cluster properties, and initial binary conditions, along with corresponding header files containing column definitions.

Appendix A Additional tables

References

All Tables

All Figures

Fig. 1 Schematic representation of the origins and eventual locations of ULX progenitors and ULXs in GCs. Progenitors can either form and remain bound to the cluster, leading to in-cluster ULX emission (“INCLUSTER”), or be dynamically ejected and emit as ULXs outside the cluster (“ESCAPER”). Additionally, ULXs originating and emitting in the field (“FIELD”) represent progenitors formed without direct interaction with cluster dynamics.

In the text

	Fig. 2 Evolution of ULX numbers in GCs. Expected number of ULXs (E[NULX]) as a function of time since GC formation. All simulations with non-zero expected rates are presented. The line indicates the number of simulations with non-zero predictions (i.e., number of dots in this time bin). Each time bin spans 300 Myr and the points are centered on these bins.
In the text

Fig. 3 Comprehensive overview of the E[NULX] evolution since the formation of a GC. Each panel represents a set of simulations with identical parameters: number of populations (Npop), concentration (Rh,2/Rh,1), binary fraction (fbin), and the presence of new features (NF). Within each panel, ULX populations from various environments are presented: tidally filling (TF) GC, non-tidally filling (nTF) GC, nodyn (binaries evolving without dynamical interactions), and escapers (progenitors formed in a GC, but experiencing the ULX phase after ejection from the cluster). See Table A.1 for simulations parameters overview. Bars and points represent the E[NULX] in 300 Myr bins.

In the text