Could kilomasers pinpoint supermassive stars?

A strong nuclear kilomaser, W1, has been found in the nearby galaxy NGC 253, associated with a forming super star cluster. Kilomasers could arise from the accretion disc around supermassive stars (>10^3 Msun), hypothetical objects that have been proposed as polluters responsible for the chemical peculiarities in globular clusters. The supermassive stars would form via runaway collisions, simultaneously with the cluster. Their discs are perturbed by stellar flybys, inspiralling and colliding stars. This raises the question if an accretion disc would at all be able to survive in such a dynamic environment and mase water lines. We investigated what the predicted maser spectrum of such a disc would look like using 2D hydrodynamic simulations and compared this to the W1 kilomaser. We derived model maser spectra from the simulations by using a general maser model for appropriate disc temperatures. All our model discs survived. The model maser spectra for the most destructive case for the simulations of M = 1000 Msun are a reasonable match with the W1 kilomaser spectrum in terms of scaling, flux values and some of the signal trends. Details in the spectrum suggest that a star of a few 1000 Msun might fit even better, with 10,000 Msun clearly giving too large velocities. Our investigations thus support the hypothesis that kilomasers could pinpoint supermassive stars.


INTRODUCTION
Kilomasers are much more luminous than masers from normal massive star formation sites in the Milky Way, but much less luminous than megamasers in Active Galactic Nuclei (AGN). Recently, a highresolution spectrum became available of a nuclear kilomaser, W1, found in the nearby galaxy NGC 253 linked to a young massive cluster (Gorski et al. 2019, previously observed by Brunthaler et al. (2009 and Hofner et al. (2006)). The association of this kilomaser and the ones found in NGC 4038 and NGC 4039 (Brogan et al. 2010) with young massive clusters is quite striking. Krause et al. (2020) proposed that kilomasers could arise in accretion discs around supermassive stars (SMSs). Denissenkov & Hartwick (2014) and Gieles et al. (2018) suggested that the SMSs, believed to be of a mass of at least 10 3 M , form via runaway collisions, simultaneously with a young massive cluster, usually in a very dense central region. They also proposed these so far hypothetical stars as polluters responsible for the chemical peculiarities in old globular clusters (GCs).
Globular clusters were traditionally believed to be simple, single population objects born in one coeval formation event with no internal chemical evolution. This view was completely revised in recent decades. For example, the 2 nd parameter problem (Sandage & Wildey 1967) arose in the 1960's when a sample of globular cluster colourmagnitude diagrams revealed that GCs with the same metallicity have different horizontal branch morphologies. GCs' with peculiar chem-★ E-mail: k.nowak@herts.ac.uk ical compositions, including large variations in certain elements, were detected already in the 1970's (Osborn 1971), however their origins were then assigned to internal deep mixing processes as a consequence of the evolution of stars.
Spectroscopic measurements in the 2000's showed that most Galactic and extra-galactic GCs demonstrate multiple sequences in the colour-magnitude diagram (Anderson 2002;Bedin et al. 2004;Piotto et al. 2007Piotto et al. , 2015Milone et al. 2012Milone et al. , 2013Martocchia et al. 2018), considered to be a result of a spread in helium abundance (Δ , Norris 2004;D'Antona et al. 2005;Charbonnel 2016;Chantereau et al. 2016) and proving that globular clusters host multiple stellar populations Milone et al. 2018). Most globular clusters show no spread in iron abundance but display a similar maximum sodium enhancement. Helium abundance spreads vary (Δ 0.1) from cluster to cluster, but are generally low Milone et al. 2018;Lardo et al. 2018). GCs display large variations in light elements: Na-O, C-N and Mg-Al anticorrelations (Denissenkov & Denisenkova 1990;Langer et al. 1993;Ventura et al. 2001;Prantzos et al. 2007;Gratton et al. 2012;Charbonnel 2016;Prantzos et al. 2017). Still the most noticeable feature in most globular cluster is the Na-O anticorrelation. A hot-hydrogen burning environment is needed to vary those abundances with the concurrent p-capturing reactions of the CNO-cycle ( 20 MK), NeNa ( 45 MK) and MgAl ( 70 MK) chains leading to the rise of those anticorrelations (Gratton et al. 2012;Prantzos et al. 2017).
Most models, seeking to explain the anomalies in globular clus-ters, refer to self-enrichment, where certain stars, polluters, within a cluster are capable of enriching other stars within the same cluster. It is also vital to include in the models how the observed amount of sodium can be produced and subsequently be accreted by the low mass proto-stars in the GCs; a requirement that multiple generation models struggle to meet. This issue is generally referred to as the 'mass budget problem' Gieles et al. 2018). In order to explain the above-mentioned anomalies three potential polluters have received a lot of attention in the literature: Asymptotic Giant Branch (AGB) stars (Ventura et al. 2001), fast-rotating massive stars (FRMS, Decressin et al. 2007;Krause et al. 2013) and supermassive stars (Denissenkov & Hartwick 2014). The nucleosynthesis of the first two proposed polluters does not correspond to the ones of GCs Prantzos et al. 2017). AGB stars display O-Na correlation instead of the anticorrelation observed (Forestini & Charbonnel 1997;Denissenkov & Herwig 2003;Karakas & Lattanzio 2007;Siess 2010;Ventura et al. 2013;Doherty et al. 2014;Renzini et al. 2022), furthermore it releases He-burning products, that are not widely detected in GCs (Karakas et al. 2006;Decressin et al. 2009;Yong et al. 2014). FRMSs on the other hand could be able to produce Mg-Al anticorrelations 1 but simultaneously would show a strong He enrichment (Decressin et al. 2007;Martins et al. 2021). The essential central temperature to activate the MgAl chain is reached by a supermassive star at the very beginning of its evolution when the abundance of He is still low (Prantzos et al. 2017). Thereupon, in its early evolutionary phase, the H-burning products show agreement with various anticorrelations observed in the GCs (Denissenkov & Hartwick 2014;Denissenkov et al. 2015). The SMS is assumed to be fully convective and therefore releases the material at the very beginning of the main sequence phase in a radiatively driven wind. The ejecta would then mix with star-forming gas that either accretes onto proto-stars or collapses to form stars independently (Krause et al. 2020). The model of concurrent formation of proto-GCs and SMS proposed by Gieles et al. (2018) provides the correct chemical patterns through the 'conveyor-belt' production of hot-H burning yields that also nicely solves the mass budget problem. The biggest disadvantage of the SMS model is arguably that no such objects are known to date (cf. Renzini et al. 2022).
Supermassive stars have also been proposed as a candidates for progenitors of supermassive black holes, but in this scenario they are expected to form through very rapid accumulation of gas (Begelman 2010;Schleicher et al. 2013;Haemmerlé et al. 2019).
Supermassive stars would be extragalactic objects, surrounded by gas and dust as frequently observed for young, massive clusters (e.g., Hollyhead et al. 2015), resulting in higher extinction. Effective temperature range between less than 10,000 K to ≈ 40,000 K (Gieles & Charbonnel 2019;Martins et al. 2020). For example for the hot case the star would be classified as blue, which would make such a star similar to normal massive stars, hard to resolve in far-away massive enough clusters and likely heavily absorbed in the gas-rich environment. This could be better for the case of a bloated and cooler star.
Proto-stars in general have accretion disc and hence we would expect an SMS to also have one. An alternative method to detect them could therefore be MASER emission from any accretion disc around an SMS. GHz MASERs are often associated with massive star formation (Ellingsen et al. 2018;Billington et al. 2019). It is unclear, however, if an accretion disc would survive in an environment as 1 Assuming a significant increase of the 24 Mg( , ) reaction with respect to current predicted estimates (Decressin et al. 2007). dynamic as expected in the centre of a forming massive star cluster. A high rate of flybys, inspirals and collisions with the SMS might well inhibit the formation of the large molecular column with low velocity shear that are required for the emission of the maser line.
Flybys are known to be able to disrupt accretion discs down to the periastron radius of the flyby (Clarke & Pringle 1993;Cuello et al. 2018;Vorobyov et al. 2017Vorobyov et al. , 2020. For the SMS case, many stars will even collide with the SMS, but they should have much smaller mass than the SMS. Understanding the final outcome requires numerical simulation. Here, we first summarise the observational arguments, why the W1 kilomaser is a candidate for an SMS accretion disc (Sect. 2), we then outline our hydrodynamics simulations (Sect. 3) with which we are able to demonstrate that SMS are expected to maintain accretion discs and compare the results to observations of the W1 kilomaser (Sect. 4). We discuss our findings in Sect. 5, arguing that model details are consistent with an SMS mass of a few 1000 M for the W1 system and summarise our conclusions in Sect. 6.

THE W1 KILOMASERS AS A CANDIDATE FOR AN ACCRETION DISC AROUND A SUPERMASSIVE STAR
The W1 kilomaser is an H 2 O maser observed at 22.2 GHz. It is located in the starburst galaxy NGC 253 at a distance of 3.5 Mpc (Rekola et al. 2005) and associated with source 11 in Leroy et al. (2018). The still embedded super star cluster has an age of 1-2 Myr and a mass of 4×10 5 M (compare their Table 2). Other kilomasers have likewise been associated with intense star formation, e.g., in the Antennae galaxies (Darling et al. 2008, and references therein), and where spatial resolution was sufficient, they have been directly associated with super star clusters (Brogan et al. 2010).
The spectrum of the W1 kilomaser is shown in Fig. 1. It has three distinct line systems: the prominent one in systemic velocity of 116 km s −1 and two 'high-velocity' features on either side of the systemic velocity at substantially lower flux. If a maser shows two or the three corresponding lines (or line systems), it is called a clean disc maser (Pesce et al. 2015), as typically observed for AGN megamasers. The spectrum of W1 looks like a disc maser spectrum, Figure 2. Luminosity versus velocity spread plot for water masers from different sources. Red circles indicate megamasers from AGN sources, green one denotes stellar masers from massive young stellar objects (MYSOs), purple circles show kilomasers, the blue star is the W49N and the black star is the W1 kilomaser. See Table 1 for references. but is much weaker than the typical AGN megamasers, about two orders of magnitude lower in luminosity.
Some extragalactic kilomasers in super star clusters have been compared to the Galactic high mass star forming region W49N, located roughly 11.1 kpc away (Gwinn et al. 1992). This region produces a large number of highly variable 22 GHz H 2 O maser spots with the total luminosity of ∼ 1 L (Zhang et al. 2013;Shakhvorostova et al. 2019;Volvach et al. 2020). The spectrum consists of 316 individual narrow lines between velocities -352.1 and 375.5 km s −1 (McGrath et al. 2004), but in all cases, the extragalactic kilomasers appear to have a more peaked and narrower spectrum. The best available kilomaser spectrum is probably from W1 (Gorski et al. 2019). The spectrum of W1 is clearly different from the one of W49N and instead consistent with the clean disc maser spectrum as argued above. Figure 2 compares velocity spread versus luminosity of the W1 kilomaser with water masers from different types of sources, including AGNs, massive young stellar objects (MYSOs) and other kilomasers, associated with other forming super star clusters and W49N. The spread in velocity space is a few hundred km s −1 for megamasers and W49N, whereas for the W1 kilomaser, it is around 80 km s −1 . We have extracted this data from various references displayed in Table 1. The relationship between luminosity and velocity spread is superlinear, and extragalactic kilomasers including the W1 kilomaser are in the middle between massive YSOs and AGNs. Hence, it would seem reasonable that kilomasers belong to an object class that is in between normal massive stars 10 2 M and AGNs 10 6 M .

AN ACCRETION DISC MODEL FOR COLLISIONALLY MAINTAINED SUPERMASSIVE STARS
In the following, we first derive constraints on the model parameters from relevant observations and then describe the setup of the hydrodynamic simulations used to test the stability of an SMS accretion disc in the dynamical environment in the centre of a massive young cluster. Our aim is to model the collisionally pumped 'high-velocity' features of the W1 kilomaser, which contain information on the mass of the central object. We do not attempt to model radiatively pumped features near systemic velocity.

General model
The similarity of the W1 maser spectrum to AGN megamaser discs suggests that it may be useful to apply the physics learnt from AGN to kilomasers. The 'high-velocity' features can be used to calculate the radius of the disc, as they represent emission closer to the outer edges of the disc in AGNs. Hence assuming this also applies to the disc around a supermassive star and the kilomaser it produces, we take the velocity spread out to both 'high-velocity' features in the W1 kilomaser to be Δ = 83 km s −1 . This measurement was performed at a digital copy of the original data kindly provided by Mark Gorski. We assume that all the high-velocity maser spots are on circular Keplerian orbits on the midline of the disc (Pesce et al. 2015). Keplerian dynamics then gives the radius of the maser spot as: where is gravitational constant and M ★ is mass of the central star.
The high-velocity maser spots in the W1 kilomaser might be located in the closest area to the inner disc that can mase the water line and not necessarily in the vicinity of the outer edges of the disc. We compare this to candidate discs of high-mass protostars by Cesaroni et al. (2007). One of the candidates listed, IRAS 20126+4104, is the best-studied case of a Keplerian accretion disc around high mass star. The disc's mass is about 4 M with a radius of 1600 au around a massive (7 M ) YSO (Cesaroni 2005). The most common value of radius of the disc is 500 au and some values are as high as 20,000 au 2 (see also Ilee et al. 2018;Sanna et al. 2021;Williams et al. 2022). Taking all the above arguments into consideration we simulated accretion discs with the inner and outer radius of in = 10 au and out = 500 au or 1000 au. We stress that the inner radius chosen is not meant to imply that the disc would not continue much further towards the star. The disc in this region would, however, likely to be too hot to allow the formation of maser spots (compare below).
Masses of accretion discs vary substantially between objects and do not scale with the radius of the disc or the mass of the central star (Cesaroni et al. 2007). We decided for 1 per cent (similar to the flyby simulations of Cuello et al. (2018)) and 10 per cent of ★ .
Following Cuello et al. (2018), we assumed a power law surface density at the start of the calculation: giving a height varying with radius as: The density of the disc is then given by: where we have used Eq. 2 for the second equality. The density outside the disc is set to 10 −20 g cm −3 . Temperature in discs around massive stars is known from observations and radiative transfer calculations to be similar to the stellar surface temperature close to the star and then drops off with a powerlaw index between -0.5 and -1, where observations seem to indicate an index closer to -0.5 (Lesniak & Desch 2011;Akiyama et al. 2013;Vural et al. 2014;Qiao et al. 2022). Supermassive stars were modelled to have temperatures between about 10 kK to 40 kK and radii between 10 2 and few times 10 3 (Gieles et al. 2018;Gieles & Charbonnel 2019). Our resulting estimate for the mass of the hypothetical SMS in the W1 system depends linearly on the radius of the maser spots and hence linearly on the temperature and the radius of the SMS. We therefore chose a conservative estimate of the effective temperature, eff , of the star of 9 kK and a radius of 100 R . We hence set the temperature as a function of radius to a scaled up law from Cuello et al. (2018): where the first value is the temperature at the inner edge of the disc in our simulations. Up to around 100 au radius the temperature is higher than 300 K and hence allows for the water lines to be mased.
A locally isothermal equation of state is used with the temperature as defined by Eq. 6. The angular velocity is calculated from the below equation of radial hydrostatic equilibrium equation, which is preserved by the balance of gravitational acceleration with centrifugal acceleration and the pressure gradient: where pressure is = B p , with Boltzmann constant B , mean molecular weight = 2.35 (Kimura et al. 2016), p as a mass of a proton and , temperature expression mentioned above. Therefore the angular velocity for the SMS disc is: According to Clarke & Pringle (1993) the most destructive encounter, where the disc looses around 50 per cent of its mass (if the ratio of the outer edge of the disc and the periastron distance, peri , is 0.8) is prograde and coplanar, whilst the retrograde and coplanar one has almost no impact on the disc. In the scenario, where orbital and disc planes are close to orthogonal, the perturber plunges through the disc without having any previous interaction with. Therefore the disc particles located within periastron remain bound and sustain their alignment, while the rest of the material is transferred out on to orbits inclined to the original orbital plane. Therefore we will only consider the case where the encounter affects the disc in the most destructive manner, namely prograde and coplanar (inclination of the orbit of 0 • angle).
Flyby simulations are typically done with similar masses of the perturbers and the star with the accretion disc. For the case we consider expected perturber mass is much lower than the expected mass of the SMS. However, there are many perturbers expected which can add up to a mass comparable to the mass of the SMS. The mass of the perturbers is estimated from Figure 3 in Gieles et al. (2018). We took the 'metal-poor' case with the number of stars = 10 7 . At the time when sms is equal 10 3 M , the total mass of all the stars in the cluster is 10 6.2 M . Thus the average mass of the pertuber is pert ∼ 0.2 M . Using the same plot from Gieles et al. (2018) we estimated the number of collisions and their frequency. We took the initial mass in stars 0 = 10 6 M , the cluster half-mass radius of h0 2.3 pc and the initial mass of each star, 0 0.1 M . The accretion factor, , for the mass in the cluster i.e. current mass of the stars (when SMS = 1000 M ) divided by the initial mass of the stars, is ∼ 2. The value higher than that would corresponds to two-body relaxation becoming more important (Gieles et al. 2018). Therefore the total mass of the stars in the cluster is: which gives the value of 2×10 6 M . The radius, in which the mass of the stars is contained is 0 = a −3 h0 = 0.085 pc (17500 au). We are only interested in the region of i = 1000 au, as this is the highest value of the accretion disc's radius we are simulating. The density of stars: gives 4.45 × 10 −7 stars per (au) 3 . Hence the number of stars, , in the region of interest, : is 1864. The crossing times is defined as cross = 2 i / rms , where rms , is the one-dimensional velocity dispersion of the cluster: Therefore, the crossing time, cross , is 60 years. Applying the values calculated above, the flyby rate = / cross 30 per year. We are only interested in the collisions that occur in the plane of the disc, 10 per cent of the solid angle, and only in the prograde ones. Therefore the flyby rate is around 1 per year to the order of the magnitude. Mark Gieles (private communication) has kindly provided his estimate of the collision rate of 0.025 per year, and for the ones that occur in the plane of the disc is equal to 0.0025 flybys per year, therefore around 1 per 400 years. This is certainly a lower limit, as it only takes into account the actual hits, not the near misses that can still impact the accretion disc significantly. Guided by these constraints, we simulated flyby rates from 10 −2 yr −1 to 1 yr −1 The initial position of each perturber is 10 times the value of the SMS disc's outer radius, out . As mentioned in Cuello et al. (2018), by setting this initial value it removes any artificial effects caused by a sudden introduction of the perturber close to the disc. The perturber is entered in the negative and direction, and leaves the grid towards negative -direction but positive -direction, allowing for the prograde encounter. The initial and values of the perturber are calculated as follows: In order for the perturber to follow a parabolic orbit we use Barker's equation: where t is the sum of the mass of the central star, ★ and the mass of the perturber, pert . Subscripts 'i' and 'f' indicate the initial and final positions of the perturber, and true anomaly is: We use equations 15, 16 and 17 to calculate the position of the perturber as a function of time. According to Meire (1985) the most practical solution to the Barker's equation follows the below calculations: Hence the true anomaly of the perturber's orbit is: and its radius is determined as follows: for each perturber. Flyby rate as well as the time of the simulation allows for the number of perturbers to be calculated i.e. for 5000 years of simulation time and a flyby rate of 1 per year, we included 5000 perturbers. Periastron distance values are generated using random numbers up to the maximum value of periastron distance set as peri = 1.1 out . The start time of the perturbers is also random. That ensures the perturbers are spread out in time according to the chosen flyby rate. To ensure perturbers are coming from different angles, without changing the starting values of i and i , we rotated the initial perturber positions using the rotation matrix: where is the orbital angle and its values are generated using random numbers between 0 and 2 . We simulated only prograde encounters, hence only half of the encounters calculated above, as they are more destructive than the retrograde ones.
We defined direct stellar collision of the perturber with the central star when the periastron distances, peri ≤ 10 au. In this instance the position of perturber, p and p , as well as its radius, , are set to 0. The mass of the perturber, pert is then added to the mass of the central star, ★ . Direct collisions are purely determined by the values of peri for each individual stellar perturber. Therefore their number is random. For the case of flyby rate of one stellar perturber per 100 years there is no direct collision ( peri ≤ 10 au).
Similarly to the work of Vorobyov et al. (2017Vorobyov et al. ( , 2020 the perturber potentials are smoothed (Klahr & Kley 2006), such that the total potential is given by: where the cubic perturber's potential Φ p is defined as: with the distance from the perturber and the smoothing radius sm = 0.13 , where is the distance of the given cell from the SMS. The top term in the above equation has an additional term in the square brackets, which artificially reduces the strength of the potential inside sm . Since the resolution is limited and the region is unresolved, the value of the acceleration would be wrong. The application of the above term allows for the correct calculation of the accelerations on the scales that can be resolved, minimising the errors.

Computational setup
To simulate the accretion disc of the supermassive star, the finite volume fluid dynamics code PLUTO is used (Mignone et al. 2007). It has been designed to integrate a system of conservation laws given by: where the state vector represents a set of conservative variables, ( ) describes fluxes of each component of the state vector and ( ) represents the source terms. For the case of ordinary hydrodynamics, Eq. 25 reduces to the following Euler equations: where is the mass density, is the velocity, is the thermal pressure and the gravitational potential is Φ. The equation of state provides closure: where the sound speed is a function of radius s ( ), as given by Eq. 6. We run PLUTO on the University of Hertfordshire high performance computing cluster (UHHPC). The code is set up to work with non-dimensional units and hence the dimensionalisation to c.g.s. units is essential for the simulation. We tested the simulation setup for a single parabolic stellar flyby of both perturber and central star of mass 1 M and successfully reproduced the results published by Cuello et al. (2018). The HD equations are solved in two dimensions in spherical coordinates. The spatial integration is performed using linear reconstruction of the fluxes to the cell boundaries and time evolution is computed using a second order Runge-Kutta scheme. The Harten-Lax-van Leer-Contact discontinuity (HLLC) Riemann solver is used to solve the numerical fluxes. The radial grid is defined from 10 au to 5000 au for out = 500 au and from 10 au to 10,000 au for out = 1000 au using a logarithmic scale. The azimuthal, uniform Table 2. Set of simulations of accretion disc around the supermassive star, varying its mass, mass and outer radius of the accretion disc, periastron distance and flyby rate. grid is set from 0 to 2 with periodic boundary. We used a resolution of 50x100 grid elements for radial and azimuthal direction respectively. We performed a set of simulations, varying different parameters outlined in Table 2. For the simulation runs SMS-mrd and SMS-Mrd the orbital period at our typical outer disc radius of 500 au is 355 and 112 years, respectively. We run the simulations without any perturber for 5000 years, which allows the disc to settle from initial high-amplitude oscillations. Smaller oscillations are ignored since the high number of stellar flybys and collisions will have a more substantial impact on the disc than the oscillations itself. After this initial period the actual run-time of the simulations starts, where the shortest run is around 1500 years for simulations with SMS = 10,000 M ★ , allowing the accretion disc to fully rotate at least 13 times.

Model maser spectrum
We derived the model maser spectra and plotted it together with the W1 kilomaser from Gorski et al. (2019), kindly provided in electronic form by Mark Gorski, for direct comparison. The systemic velocity of W1 in NGC 253 is 116 km s −1 . For the integration along the line of sight, we map the simulation output on to a Cartesian grid and assume the source to be observed edge-on. Each of these Cartesian cells has one particular velocity and is regarded as a velocity coherent maser region or a part thereof with length d . We need to adapt values for the linewidth and the maser spot size to determine the spectrum. The values for these can vary: for extragalactic kilomasers those values have not been measured, for AGN, a typical value is 3 km s −1 (Kartje et al. 1999) and for the Galactic kilomaser W49N is 0.5-1.9 km s −1 (Liljeström & Gwinn 2000). As an example we take 1 km s −1 , as with this value we already obtain enough flux to reproduce the observations, with a higher linewidth we get even higher fluxes. For regions along the line of sight that have the same velocity within a fiducial velocity coherence range of 1 km s −1 , d is multiplied by the number of these regions, , to obtain the total length of the velocity coherent region. The resolution of the grid is approximated based on the maser spot size of ∼ 10 au which is also roughly consistent with the marginally resolved maser spots in W49N from Zhang et al. (2013). If we decrease the spot size, the flux increases.
Fluxes are calculated using the flux equation to model AGN maser discs derived by Kartje et al. (1999): where 10 = /10, and defines the effective aspect ratio: For a larger than a few the maser is saturated. We assumed 10 = 1. The half-length of the maser emission region that lies in a disc along the line of sight is defined as , is the distance to the source with the observed area . That results in Eq. 29, reducing to: where = 3.5 Mpc is the distance to the starburst galaxy NGC 253. The 22 GHz H 2 O line requires dense gas of at least 10 7 cm −3 and temperatures larger than 300 K (Gorski et al. 2019), hence we restrict the density for masing regions accordingly. The temperature defined by Eq. 6 is also restricted between the already mentioned minimum of 300 K up to 1500 K, as to allow the gas to remain molecular (Kartje et al. 1999). Those parameters are required for the water to be mased if it is collisionally pumped.
Since we only model collisionally pumped masers, the model spectra do not include the 'low-velocity' features, which are radiatively pumped by the the background infrared radiation, usually by the central object. For the purpose of this study, we only focus on the 'high-velocity' features in maser spectra. The model maser spectrum shows some residual positive flux at systemic velocity for each simulation run. This is the result of flux being calculated for each area in the disc that meet the density and temperature requirements for maser emission in the line of sight, thus computing the emission as a ring in the disc. According to Elitzur et al. (1991) around 51 per cent of the maser emission is lost through the sides, which would actually result in such a positive flux around zero velocity. As the resulting fluxes are low, we make no attempt to model the sideways emission accurately. This approximation needs to be taken into consideration when interpreting the results. Figure 3 (left column) shows density plots for three time steps for run SMS-mrd-hf with small, low-mass accretion disc and high perturber frequency. A slow, steady and smooth dispersion of the disc is well visible. Some gas is being spread out over 2000 au distance within 5,000 years. Spiral features develop, but the main part of the original disc remains intact. In the right column we show the comparison of the modelled H 2 O maser spectra for each time step with the W1 kilomaser (Gorski et al. 2019) in orange. The system velocity was added to the modeled velocities for direct comparison to the observed spectrum. The 'high-velocity' features are clearly visible on each side, with comparable flux value but slightly closer to the systemic velocity than the W1 kilomaser. Those values are approximately 74 and 160 km s −1 at timestep 2,500 years (top row). Using the Eq. 1 with Δ ∼ 43 km s −1 , gives a radius for the maser spots of ∼ 480 au. As expected, this is the largest radius in the model disc where the gas temperature is still high enough for water masers. As the disc evolves with time the 'high-velocity' feature on the left becomes broader with lower flux (middle row, Figure 3) between 60 to 75 km s −1 , giving the radius of the maser spots at around 380. It can also be noticed that the the 'high-velocity' feature at around 66 km s −1 at the time step of 10,000 years has a slightly higher flux to the feature at 153 km s −1 . Looking at the corresponding density plot it can be observed that a spiral arm has formed, moving clockwise and creating a long, velocity-coherent region near ∼ 350 au, whilst the maser spots corresponding to the lower-flux feature at 153 km s −1 are located at the radius of ∼ 650 au. Additionally it can be noticed from the time evolution of the maser spectra in the same figure that the spikes of the 'high-velocity' features are moving slightly inwards or outwards, depending on the location of the spiral arm, visible in the corresponding density plots. That effect is clearly visible on the left plot in Figure 4 where we compare the unperturbed maser spectrum, i.e. at the very beginning of the simulation run, to the spectrum produced at 5,000 years.

Disc stability for a 10,000 M star
Comparing the results for the case of a 10,000 M SMS (simulation run SMS-Mrd-hf) in Figure 5 to the lower mass one, keeping the other parameters the same, we notice that the disc spreads out faster, with more distinct spiral arms and more mass reaching a larger radius. It is hard to establish what maximum distance the mass could reach since we are limited by the computational grid (up to 5000 au). Since the gas rotates much faster than in the previous case, the velocity range in the maser spectrum (right column in Fig. 5) had to be extended to 800 km s −1 in order to see the 'high-velocity' features, which appear in the first snapshot (top row in Fig. 5), symmetrically located at around -13 and 240 km s −1 . Those values gives Δ ∼ 126 km s −1 , with Eq. 1 we get the radius of maser spot of ∼ 550 au, slightly inwards of the = 300 K limit above which we allow maser emission. This, again, reflects the quick build-up of structure in the gas distribution. As the disc evolves with time, the mass and spiral arms spread out more (middle row in Figure 5) and the 'high-velocity' features in the corresponding maser spectrum move slightly outwards from -13 to -26 km s −1 and from 240 to 263 km s −1 , leading to the maser spots being now located at the radius of 424 au, therefore slightly closer to the central star. At timestep of 3,000 years (bottom row) the 'highvelocity' feature on the left of the systemic velocity has broadened its peak whilst the spike on the right has increased its flux due to the movement of the spiral arm and increasing column density along the line of sight, as clearly seen on the right plot in Figure 4.

Effects of changing disc mass, radius and flyby rate
The results, especially with same parameters but varying flyby rates have a moderate impact on how the density plots and model maser spectra look like. Similarly we do not see strong differences in simulations where the mass or radius of the disc is varied. Figure 6 shows comparison of the density maps of the snapshot for selected simulation runs for similar stages of disc evolution. In the simulations with higher disc mass (top row, simulation run SMS-mRD) compared to the lower one (2nd row, simulation run SMS-mRd) it is noticeable that more disc mass is being spread out over larger distances, especially for the low flyby rate (left column). Since the gravitational acceleration is only proportional to the potential of the central star, the amount of mass in the disc does not influence how fast the gas rotates. There is, however, a hydrodynamic interaction of the disc with the background gas, which allows a disc with more inertia to spread somewhat faster. Especially for the 10,000 M , where the disc gas rotates much faster, the spiral arm spreads out more vigorously and becomes more distinct i.e. dense (3rd and bottom rows, simulation runs SMS-MrD and SMS-MRD, respectively). This also shows that the background gas, with its chosen low value, is still of some significance for the spreading of the gas. We do not notice any substantial change in the disc with radius of 500 au (3rd row) compared to the 1000 au one (bottom row), apart from the disc spreading faster for the lower radius case. This is a consequence of the lower value of the maximum periastron distance for the perturbers' orbit. We notice that the higher flyby rates make the disc spread out faster (compare each column in Fig. 6, representing different flyby rates).

DISCUSSION
We have simulated the stability of accretion discs around so far hypothetical supermassive stars in the dynamic environment of a young massive cluster. All our discs preserve a dense central region, and spread out at different rates, depending on the exact parameters we chose.
In the multiple stellar flybys around the supermassive star, a single 0.2 M stellar perturber does not have any influence on the central star and its disc. This agrees with the results from Vorobyov et al. (2020), where with decreasing perturber mass the effect on the disc is gradually lessening. However multiple perturbers of that low mass do have consequences on the disc, for all the chosen simulations parameters (compare different columns in Fig. 6). Taking the case of the most destructive flyby rate of one stellar perturber per year for the 1000 M star, where we have entered 5000 perturbers of 0.2 M mass onto the computational domain, the total mass adds up to 1000 M . The total mass of the perturbers therefore equals the mass of the central star, the SMS. Comparing the case for the 10 M accretion disc (Figure 3) with the single, equal-mass, stellar flyby seen in Cuello et al. (2018, their Figure 2, top row), the overall magnitude of the effects on the accretion discs are very comparable, in particular, the disc spreads out. We also noticed some differences: the density plots of a single perturber parabolic flyby show the disc mass being pulled out as one track, often the shape is reminiscent of spiral arms, whilst multiple flybys show much smoother structure with no preferred direction of the ejected mass from the disc. This is of course an expected consequence of the isotropic bombardment of the disc by the large number of stellar perturbers hitting over a long time allowing the mass distribution to become more homogeneous over larger distances.
Including radiative transfer in our simulation would have some effect on the disc dynamics, changing opacity and therefore heating different parts of the disc at different times. This might contribute to some expansion of the disc. The sound speed, s , for our simulation model is about 1.0 km s −1 , whilst the corresponding Keplerian velocity, k , ranges between ∼ 200 km s −1 and ∼ 42.0 km s −1 . Since s « k we would not expect the disturbances introduced by live radiative transfer to be significant.
Incorporating magnetic fields to the model, would result in the magneto-rotational instability (MRI), that would in turn lead to viscosity, disc diffusion and accretion. As a consequence the disc would spread, decrease its density and hence have a lower flux in the 'highvelocity' features of our maser model. We model a star that grows by collisions and the accretion rate by gas is assumed to be low. Hence it is implicit in the model assumption that the MRI should have small effects during the simulation time. In a realistic scenario spreading by MRI would likely be offset by the inflow of additional gas into the circumstellar environment.
Self-gravity should not play a major role for the evolution of our model discs as the ratio of the mass of the disc to the mass of the central star is small. Toomre's stability parameter, T , (Krause et al. 2013, their Equation 15) evaluated for our disc parameters is in range ∼ 3.4 -9.2. Therefore the disc is stable against gravitational collapse.
The model maser spectra produced here successfully reproduce the 'clean' disc maser case, as they have two of the features mentioned before, namely, both the 'high-velocity' features (compare Pesce et al. 2015). Moreover, as we are only simulating collisionally pumped maser emissions, therefore the systemic feature is not present in our model. In order to model those features, the simulation would have to include background infrared radiation from the central source.
To calculate the flux in the maser model we used Eq. 29 taken from Kartje et al. (1999). This expression is considered to be a general formula, assuming optimal collisional pumping conditions. The authors of the paper derive the equation, by considering a velocity-coherence 'box' in an edge-on disc, to which the maser emission is confined to (Kartje et al. 1999, their Equation 9). The results that we obtained from the simulations are in agreement with the observations, while there is some freedom in choosing spot sizes and linewidth, which both influence the predicted flux. We regard this as a success of the model.
The spectrum of the unperturbed disc in both simulations, seen in Fig 4, is not smooth. This is due to the resolution from the computational grids for simulations and spectral modelling. This leads to substantial flux changes from cell to cell as well as between snapshot times, when the maser beam direction has moved on.
The 1000 M simulations generally reproduce the spectral features of the W1 kilomaser better. Apart from similar fluxes, for both 'high-velocity' features, there is a general trend for the signal to fol- low an increase towards the systemic velocity for the left hand side of the maser for many models. However looking at the right hand side, the model does not follow the trend of the signal of the W1 kilomaser for the models in all the cases of different parameters. The 'high-velocity' features in the W1 kilomaser look different on the respective sides of the systemic velocity. This is similar to our model spectra, which also show similar differences between the two sides.
The model spectra for the higher mass SMS, 10,000 M fit worse to the W1 kilomaser. The velocity range was substantially increased in order for the two 'high-velocity' features to be visible. The model maser for the higher mass of the central star hence starts to resemble spectra for AGN megamasers, in terms of the offset of those features from systemic velocity. Thus summarising the above findings, the model spectra from simulations for the SMS = 1000 M are a better match with the W1 kilomaser spectrum from Gorski et al. (2019).
The H 2 O molecule is dissociated close to the star where > 1500 K. Those conditions do not allow for the water to stay molecular, hence the 22 GHz H 2 O lines cannot be mased. Similarly the required water level population inversion does not occur at a larger distances, where temperature remains below 300 K. The model maser spots move inwards and outwards within that distance range, depending on the movement of the spiral arm, and its density. The higher the disc mass, the more it spreads out and the denser the spiral arm is.
The brightest maser spots are located roughly in the distance between 350-650 au, for all the models and this would easily accommodate bloated SMS models that could have radii of tens of au. The 'high-velocity' features do not appear closer to the inner radius and the center, as the temperature increases closer to the center and decreases with the distance approaching the outer radius. We used temperature profile of ( ) = in ( / in ) − , where in indicates the effective temperature at the inner radius, in and denotes the power-law index; in our simulations = 0.5. Changing the value of affects the radius of the maser spots, e.g., = 0.75 (Dullemond & Monnier 2010;Vural et al. 2014) would move the brightest maser emission towards the inner radius. However, we have taken care to chose our temperature and radius for the star very conservatively. We believe it is unlikely that the mass of the central object would be much lower than our estimate.
We assumed that the disc is edge-on to the observer, but if the disc is inclined then the maser spectrum might start showing slight changes, i.e. decrease in the velocity and flux values, especially of the 'high-velocity' features. As the edge-on disc is 90 • , if the disc would then be observed at some small angle i.e. 85 • , then the flux and velocity values in the maser spectrum will decrease. The inclination is well constrained by the fact that we see maser spots, and the observed velocities are hardly affected by a quite possible small deviation from an 'edge-on' disc.
The presence of the dense spiral and their movement around the computational grid causes 'high-velocity' features to modulate the peak heights and move them inwards or outwards in the model maser spectra. The variability is in particular a consequence of the beamed nature of maser emission. Figure 7 displays two samples of maser spectrum evolution within chosen periods of 4 years for simulation SMS-mrd-hf, where either one of the 'high-velocity' features has visibly increased its flux or has slightly moved outwards from the systemic velocity. Hence it would be beneficial to observe the W1 kilomaser with the Very-Long-Baseline Interferometry (VLBI) to obtain more data for another epoch in order to see if there is a slight displacement for the peaks of the 'high-velocity' features. Those findings could potentially confirm the presence of the dense spiral arms moving around the central massive object in the forming superstar cluster in NGC 253.
We note that Levy et al. (2022) presented new 350 GHz dust emission observations of the W1 host cluster. These observations seem to indicate that the source 11 from Leroy et al. (2018) splits up into four different objects. In their analysis the largest subsystem is 20,000 M . This substructure could also relate to a forming superbubble where the dust is beginning to be displaced from the host cluster.
The 'high-velocity' peaks in the W1 kilomaser are observed at 31.6 and 199 km s −1 . The model spectra for SMS = 1000 M show a smaller velocity spread between the peaks whilst for SMS = 10,000 M it is larger. We can already conclude that in order for the model peaks to appear at the same velocity values as the W1 kilomaser, the SMS value needs to be set between these two extreme values, if the disc temperature declines in a similar way as we assumed in our simulations. We mentioned before that the maser spots for both cases are then located between 380-550 au, with the average value, out = 465 au. The mean value of velocity spread between systemic and peak velocity of W1 kilomaser spectrum is Δ = 83 km s −1 . Using those values and the Eq. 1, we estimate that SMS ≈ 4000 M would best reproduce features of the W1 kilomaser.

CONCLUSIONS
We have shown that the W1 kilomaser is a clean disc maser, according to the definition of Pesce et al. (2015), showing all three expected maser line systems (compare Sect. 2). Flux and velocity spread are intermediate between massive star and AGN masers. The flux expected for a model SMS accretion disc agrees well with the one of the W1 kilomaser.
Using a range of plausible parameters for flybys and collisions of perturbers with SMS, we show that in all cases an accretion disc survives for the entire simulation time and is able to produce prominent water masers with similar strength as observed in W1.
The results showed that maser models produced for the simulations where SMS = 1000 M are a better match for the W1 kilomaser from Gorski et al. (2019). The model exhibits two 'high-velocity' features with similar fluxes and the left-hand side of the model spectrum shows a similar trend of the flux with the velocity. For larger stars with mass of 10,000 M the 'high-velocity' peaks extend their velocity values beyond the observed spectrum, starting to resemble megamasers from AGNs. Hence we can conclude that, within our assumptions, the SMS with 10,000 M is too high a mass for the W1 kilomaser.
We included plausible upper and lower limits for the number of perturbers entering the grid at a given time, the mass of the accretion disc and its outer radius as well as the mass of the central supermassive star. The simulations show the expected variations, e.g., for the rate the disc spreads out, but overall, the model maser spectra did not show significant differences between them.
We focused on the most destructive case for the disc i.e. the flyby rate of one stellar perturber per year, as well as including only prograde and coplanar encounters. We conclude that dynamic flybys would not impact any SMS accretion disc too strongly and hence would not prevent the formation of maser spots.
Summarising the results presented here, we can confirm that the results support the hypothesis for a supermassive star being present in the forming massive star cluster, and potentially being able to produce the best match to the observed maser spectrum of the W1 kilomaser obtained with an SMS around 4000 M .