# Could kilomasers pinpoint supermassive stars?

Katarzyna Nowak, Martin Krause, Daniel Schaerer

► To cite this version:

Katarzyna Nowak, Martin Krause, Daniel Schaerer. Could kilomasers pinpoint supermassive stars?. Monthly Notices of the Royal Astronomical Society, 2022, 516 (4), pp.5507-5520. 10.1093/mn-ras/stac2547. hal-04894683

HAL Id: hal-04894683

https://hal.science/hal-04894683v1

Submitted on 19 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Could kilomasers pinpoint supermassive stars?

Katarzyna Nowak<sup>1,1*</sup> Martin G. H. Krause<sup>1,1</sup> and Daniel Schaerer<sup>2,3</sup>

Accepted 2022 September 2. Received 2022 August 29; in original form 2022 July 19

## ABSTRACT

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³ M☉), 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 M☉ 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 M☉ might fit even better, with 10 000 M☉ clearly giving too large velocities. Our investigations thus support the hypothesis that kilomasers could pinpoint supermassive stars.

**Key words:** accretion, accretion discs – hydrodynamics – masers – stars: abundances – globular clusters: general – galaxies: star clusters: general.

## 1 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 high-resolution 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, Baan & Takano (2006)). The association of this kilomaser and the ones found in NGC 4038 and NGC 4039 (Brogan, Johnson & Darling 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³ 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 second parameter problem (Sandage & Wildey 1967) arose in the 1960s when a sample of globular cluster colour–magnitude diagrams revealed that GCs with the same

metallicity have different horizontal branch morphologies. GCs with peculiar chemical compositions, including large variations in certain elements, were detected already in the 1970s (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 2000s 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. 2007, 2015; Milone et al. 2012, 2013; Martocchia et al. 2018), considered to be a result of a spread in helium abundance (ΔY; Norris 2004; D’Antona et al. 2005; Charbonnel 2016; Chantereau, Charbonnel & Meynet 2016) and proving that globular clusters host multiple stellar populations (Bastian & Lardo 2018; Milone et al. 2018). Most globular clusters show no spread in iron abundance but display a similar maximum sodium enhancement. Helium abundance spreads vary (ΔY ~ 0.1) from cluster to cluster, but are generally low (Bastian, Cabrera-Ziri & Salaris 2015; Bastian & Lardo 2018; Lardo et al. 2018; Milone et al. 2018). GCs display large variations in light elements: Na-O, C-N, and Mg-Al anticorrelations (Denissenkov & Denisenkova 1990; Langer, Hoffman & Sneden 1993; Ventura et al. 2001; Prantzos, Charbonnel & Iliadis 2007; Gratton, Carretta & Bragaglia 2012; Charbonnel 2016; Prantzos, Charbonnel & Iliadis 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).

<footnote n="1"><sup>1</sup>Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK</footnote>

<footnote n="2"><sup>2</sup>Observatoire de Genève, Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="3"><sup>3</sup>CNRS, IRAP, 14 Avenue E. Belin, F-31400 Toulouse, France</footnote>

<footnote>*E-mail: k.nowak@herts.ac.uk</footnote>
Most models, seeking to explain the anomalies in globular clusters, 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’ (Bastian & Lardo 2018; 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 (Bastian et al. 2015; 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, Marino & Milone 2022); furthermore, it releases He-burning products that are not widely detected in GCs (Karakas et al. 2006; Decressin et al. 2009; Yong, Grundahl & Norris 2014). FRMSs on the other hand could be able to produce Mg-Al anticorrelations<fnref n="1" /> but simultaneously would show a strong He enrichment (Decressin et al. 2007; Martins, Chantereau & Charbonnel 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 on to 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 10000 K to ≈ 40000 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

<figure><img src="image_2.png" /><figcaption>**Figure 1.** W1 kilomaser spectrum from Gorski et al. (2019). The dotted line indicates the systemic velocity of 116 km s⁻<fnref n="1" />. The low-velocity feature is highlighted by the shaded box. The black arrows on both sides of the systemic velocity point to the ‘high-velocity’ features at 31.6 and 199 km s⁻<fnref n="1" />.</figcaption></figure>

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 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. 2017, 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 summarize the observational arguments, why the W1 kilomaser is a candidate for an SMS accretion disc (Section 2), we then outline our hydrodynamics simulations (Section 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 (Section 4). We discuss our findings in Section 5, arguing that model details are consistent with an SMS mass of a few 1000 M☉ for the W1 system and summarize our conclusions in Section 6.

## 2 THE W1 KILOMASERS AS A CANDIDATE FOR AN ACCRETION DISC AROUND A SUPERMASSIVE STAR

The W1 kilomaser is an H₂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 \times 10^5$ M☉ (compare their table 2). Other kilomasers have likewise been associated with intense star formation, e.g. in the Antennae galaxies (Darling, Brogan & Johnson 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⁻<fnref n="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

<footnote n="1">¹Assuming a significant increase of the ²⁴Mg(*p*, *γ*) reaction with respect to current predicted estimates (Decressin et al. 2007).</footnote>
<figure><img src="image_3.png" /><figcaption>**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 MYSOs, purple circles show kilomasers, the blue star is the W49N and the black star is the W1 kilomaser. See Table 1 for references.</figcaption></figure>

clean disc maser (Pesce et al. 2015), as typically observed for AGN megamasers. The spectrum of W1 looks like a disc maser spectrum, 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, Moran & Reid 1992). This region produces a large number of highly variable 22 GHz H$_2$O maser spots with the total luminosity of ~1 L$_\odot$ (Zhang et al. 2013; Shakhvorostova et al. 2019; Volvach, Volvach & Larionov 2020). The spectrum consists of 316 individual narrow lines between velocities -352.1 and 375.5 km s$^{-1}$ (McGrath, Goss & De Pree 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.

Fig. 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 super-linear, 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$_\odot$ and AGNs ≳ 10$^6$ M$_\odot$.

## 3 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.

### 3.1 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 $\Delta\nu = 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

$$R_{\text{out}} = \frac{GM_{\star}}{(\Delta\nu)^2}, \qquad (1)$$

where G is gravitational constant and M$_\star$ is mass of the central star. For a 1000 M$_\odot$ star, we would have $R_{\text{out}} \sim 129$ au (~6.2 × 10$^{-4}$ pc). 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$_\odot$ with a radius of 1600 au around a massive (7 M$_\odot$) 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 $R_{\text{in}} = 10$ au and $R_{\text{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 M$_\star$.

Following Cuello et al. (2018), we assumed a power law surface density at the start of the calculation

$$\Sigma = \Sigma_0 \left(\frac{R}{R_{\text{in}}}\right)^{-1} \qquad (2)$$

$$\Sigma_0 = \frac{M_{\text{disc}}}{2\pi R_{\text{in}}(R_{\text{out}} - R_{\text{in}})} \qquad (3)$$

giving a height varying with radius as

$$H = 0.05 R_{\text{in}} \left(\frac{R}{R_{\text{in}}}\right)^{1.25} \qquad (4)$$

The density of the disc is then given by

$$\rho = \frac{\Sigma}{2H} = \frac{\Sigma_0 R_{\text{in}}}{2RH}, \qquad (5)$$

<footnote>$^2$Although the authors argue that those high values of radius would suggest it to be a massive rotating structure called ‘toroid’, rather than an accretion disc around a single star. ‘Toroids’ are believed to host a stellar cluster and might not be in equilibrium, whilst accretion discs are in Keplerian orbits around their central stars (Cesaroni et al. 2007).</footnote>
<table><caption>**Table 1.** H₂O masers in different astronomical sources. Kilomasers have their host galaxy indicated in the brackets in the first column. All the stellar masers are from MYSOs, whilst megamasers are from AGN sources.</caption><thead><tr><th>Name of the source</th><th>Luminosity (L<sub>☉</sub>)</th><th>Max vel. spread (km s<sup>-1</sup>)</th><th>Classification</th><th>References</th></tr></thead><tbody><tr><td>NGC 4258</td><td>80</td><td>890</td><td>Megamaser</td><td>(Greenhill et al. 1996)</td></tr><tr><td>NGC 3079</td><td>500</td><td>126</td><td>Megamaser</td><td>(Yamauchi et al. 2004)</td></tr><tr><td>Circinus</td><td>20</td><td>260</td><td>Megamaser</td><td>(Greenhill et al. 2003)</td></tr><tr><td>UGC 3789</td><td>370</td><td>725</td><td>Megamaser</td><td>(Reid et al. 2009)</td></tr><tr><td>NGC 2960 (Mrk 1419)</td><td>400</td><td>465</td><td>Megamaser</td><td>(Henkel et al. 2002; Kuo et al. 2011)</td></tr><tr><td>IC 2560</td><td>122</td><td>323</td><td>Megamaser</td><td>(Ishihara et al. 2001)</td></tr><tr><td>NGC 3393</td><td>320</td><td>601</td><td>Megamaser</td><td>(Kondratko, Greenhill & Moran 2008)</td></tr><tr><td>J0437 + 2456</td><td>178</td><td>338.85</td><td>Megamaser</td><td>(Gao et al. 2017)</td></tr><tr><td>NGC 6323</td><td>500</td><td>600</td><td>Megamaser</td><td>(Kuo et al. 2011, 2015)</td></tr><tr><td>NGC 6926</td><td>340</td><td>200</td><td>Megamaser</td><td>(Sato et al. 2005)</td></tr><tr><td>UGC 6093</td><td>770</td><td>801.28</td><td>Megamaser</td><td>(Zhao et al. 2018)</td></tr><tr><td>W1 (NGC 253)</td><td>1.02</td><td>84.4</td><td>Kilomaser</td><td>(Gorski et al. 2019)</td></tr><tr><td>W49N</td><td>1</td><td>363.8</td><td>Kilomaser</td><td>(McGrath et al. 2004)</td></tr><tr><td>H<sub>2</sub>O-East (NGC 4038/NGC 4039)</td><td>1.3</td><td>13</td><td>Kilomaser</td><td>(Brogan et al. 2010)</td></tr><tr><td>H<sub>2</sub>O-SE (NGC 4038/NGC 4039)</td><td>4.1</td><td>50</td><td>Kilomaser</td><td>(Brogan et al. 2010)</td></tr><tr><td>H<sub>2</sub>O-West (NGC 4038/NGC 4039)</td><td>7.7</td><td>25</td><td>Kilomaser</td><td>(Brogan et al. 2010)</td></tr><tr><td>He 2-10</td><td>0.68</td><td>70</td><td>Kilomaser</td><td>(Darling et al. 2008)</td></tr><tr><td>G229.5711 + 00.1525</td><td>8.5 × 10<sup>-5</sup></td><td>10.25</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G220.4587-00.6081</td><td>2.29 × 10<sup>-6</sup></td><td>5</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G083.7071 + 03.2817</td><td>4.79 × 10<sup>-7</sup></td><td>7.75</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G080.8624 + 00.3827</td><td>2.88 × 10<sup>-6</sup></td><td>6.9</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G081.8652 + 00.7800</td><td>4.27 × 10<sup>-4</sup></td><td>37.5</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G081.7131 + 00.5792</td><td>3.55 × 10<sup>-6</sup></td><td>11.5</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G081.7624 + 00.5916</td><td>6.03 × 10<sup>-7</sup></td><td>23.85</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G084.1940 + 01.4388</td><td>5.13 × 10<sup>-7</sup></td><td>20.6</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G084.9505-00.6910</td><td>1.55 × 10<sup>-6</sup></td><td>2.45</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G095.0531 + 03.9724</td><td>1.62 × 10<sup>-4</sup></td><td>4.7</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G094.2615-00.4116</td><td>9.55 × 10<sup>-6</sup></td><td>20.6</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G094.4637-00.8043</td><td>1.26 × 10<sup>-5</sup></td><td>2.55</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr><tr><td>G094.6028-01.7966</td><td>2.29 × 10<sup>-4</sup></td><td>10.2</td><td>Stellar maser</td><td>(Urquhart et al. 2011)</td></tr></tbody></table>

where we have used equation (2) for the second equality. The density outside the disc is set to $10^{-20}$ g cm⁻³.

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 power-law 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, $T_{\text{eff}}$, of the star of 9 kK and a radius of 100 R<sub>☉</sub>. We hence set the temperature as a function of radius to a scaled up law from Cuello et al. (2018),

$$T = 1920 \text{ K} \left( \frac{R}{R_{\text{in}}} \right)^{-1/2}, \qquad (6)$$

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 equation (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

$$\frac{1}{\rho} \frac{\partial p}{\partial R} = \Omega^2 R - \frac{GM_*}{R^2}, \qquad (7)$$

where pressure is $p = \frac{\rho k_B T}{\mu m_p}$, with Boltzmann constant $k_B$, mean molecular weight $\mu = 2.35$ (Kimura, Kunitomo & Takahashi 2016), $m_p$ as a mass of a proton and T, temperature expression mentioned above. Therefore, the angular velocity for the SMS disc is

$$\Omega^2 = \frac{GM_*}{R^3} - 0.5 \frac{1920Kk_B}{\mu m_p R} \left(\frac{R}{R_{\text{in}}}\right)^{-0.5}. \qquad (8)$$

According to Clarke & Pringle (1993) the most destructive encounter, where the disc loses around 50 per cent of its mass (if the ratio of the outer edge of the disc and the periastron distance, $R_{\text{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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Time evolution of the reaction of an accretion disc around an SMS with $M_{\text{SMS}} = 1000$ M$_{\odot}$, $R_{\text{out}} = 500$ au, and $M_{\text{disc}}/M_{\text{SMS}} = 1$ per cent to a flyby rate of one stellar perturber per year. The left column shows density plots with corresponding maser spectra on the right for selected time steps. The model spectrum is shown in blue, whilst the W1 kilomaser from Gorski et al. (2019) for comparison is plotted in orange. See caption for Fig. 1 for the definition of different features in the maser spectrum.</figcaption></figure>

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 Fig. 3 in Gieles et al. (2018). We took the ‘metal-poor’ case with the number of stars $N = 10^7$. At the time when $m_{\text{sms}}$ is equal $10^3$ M$_{\odot}$, the total mass of all the stars

in the cluster is $10^{6.2}$ M$_{\odot}$. Thus the average mass of the pertuber is $M_{\text{pert}} \sim 0.2$ M$_{\odot}$.

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 $M_0 = 10^6$ M$_{\odot}$, the cluster half-mass radius of $R_{h0} \simeq 2.3$ pc and the initial mass of each star, $m_0 \simeq 0.1$ M$_{\odot}$. The accretion factor,
a, for the mass in the cluster i.e. current mass of the stars (when $M_{SMS} = 1000 \text{ M}_\odot$) divided by the initial mass of the stars, is ~ 2. The value higher than that would correspond to two-body relaxation becoming more important (Gieles et al. 2018). Therefore, the total mass of the stars in the cluster is

$$M = aM_0, \quad (9)$$

which gives the value of $2 \times 10^6 \text{ M}_\odot$. The radius, in which the mass of the stars is contained is $r_0 = a^{-3} R_{h0} = 0.085 \text{ pc}$ (17 500 au). We are only interested in the region of $r_i = 1000 \text{ au}$, as this is the highest value of the accretion disc's radius we are simulating. The density of stars,

$$n = \frac{3M}{4\pi r_0^3 a m_0}, \quad (10)$$

gives $4.45 \times 10^{-7}$ stars per (au)$^3$. Hence the number of stars, N, in the region of interest, $r_i$:

$$N = \frac{4\pi r_i^3 n}{3}, \quad (11)$$

is 1864. The crossing times is defined as $t_{cross} = 2r_i/V_{rms}$, where $V_{rms}$, is the one-dimensional velocity dispersion of the cluster

$$V_{rms} \approx \sqrt{\frac{GM}{6R_{ho}}} = 159 \text{ km s}^{-1}. \quad (12)$$

Therefore, the crossing time, $t_{cross}$, is 60 yr. Applying the values calculated above, the flyby rate $= N/t_{cross} \approx 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 yr. 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}$ to 1 per year.

The initial position of each perturber is 10 times the value of the SMS disc's outer radius, $R_{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 y and x direction, and leaves the grid towards negative y-direction but positive x-direction, allowing for the prograde encounter. The initial y and x values of the perturber are calculated as follows:

$$y_i = -10R_{out} \quad (13)$$

$$x_i = -2R_{\text{peri}} \sqrt{1 - \frac{y_i}{R_{\text{peri}}}}. \quad (14)$$

In order for the perturber to follow a parabolic orbit we use Barker's equation,

$$\Delta t = \sqrt{\frac{2R_{\text{peri}}^3}{GM_t}} \left( D_f + \frac{1}{3} D_f^3 - D_i - \frac{1}{3} D_i^3 \right), \quad (15)$$

where $M_t$ is the sum of the mass of the central star, $M_*$ and the mass of the perturber, $M_{\text{pert}}$. Subscripts 'i' and 'f' indicate the initial and final positions of the perturber,

$$D_{i,f} = \tan \left( \frac{v_{i,f}}{2} \right) \quad (16)$$

and true anomaly is

$$v_{i,f} = \arctan \left( \frac{x_{i,f}}{y_{i,f}} \right). \quad (17)$$

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:

$$A = \frac{3}{2} \sqrt{\frac{2GM_*}{2R_{\text{peri}}^3}} (t - t_0) \quad (18)$$

$$B = \left( A + \sqrt{A^2 + 1} \right)^{\frac{1}{3}}. \quad (19)$$

Hence the true anomaly of the perturber's orbit is

$$v = 2 \arctan \left( B - \frac{1}{B} \right) \quad (20)$$

and its radius is determined as follows:

$$R = \frac{2R_{\text{peri}}}{1 + \cos \nu} \quad (21)$$

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 yr 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 $R_{\text{peri}} = 1.1 R_{\text{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 $y_i$ and $x_i$, we rotated the initial perturber positions using the rotation matrix,

$$R = \begin{bmatrix} x_p \cos \theta - y_p \sin \theta \\ x_p \sin \theta + y_p \cos \theta \end{bmatrix}, \quad (22)$$

where $\theta$ is the orbital angle and its values are generated using random numbers between 0 and $2\pi$.

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, $R_{\text{peri}} \le 10$ au. In this instance the position of perturber, $x_p$ and $y_p$, as well as its radius, $R$, are set to 0. The mass of the perturber, $M_{\text{pert}}$ is then added to the mass of the central star, $M_*$. Direct collisions are purely determined by the values of $R_{\text{peri}}$ for each individual stellar perturber. Therefore, their number is random. For the case of flyby rate of one stellar perturber per 100 yr there is no direct collision ($R_{\text{peri}} \le 10$ au).

Similarly to the work of Vorobyov et al. (2017, 2020) the perturber potentials are smoothed (Klahr & Kley 2006), such that the total potential is given by

$$\Phi = \Phi_* + \Phi_p = -\frac{GM_*}{r} + \Phi_p, \quad (23)$$

where the cubic perturber's potential $\Phi_p$ is defined as

$$\Phi_p = \begin{cases} -\frac{m_p G}{d} \left[ \left( \frac{d}{r_{\text{sm}}} \right)^4 - 2 \left( \frac{d}{r_{\text{sm}}} \right)^3 + 2 \frac{d}{r_{\text{sm}}} \right] & \text{for } d \le r_{\text{sm}} \\ -\frac{m_p G}{d} & \text{for } d > r_{\text{sm}}, \end{cases} \quad (24)$$

with the distance from the perturber $d$ and the smoothing radius $r_{\text{sm}} = 0.13r$, where $r$ 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
<table><caption>**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.</caption><thead><tr><th>Label</th><th><i>M</i><sub>SMS</sub><br>(M<sub>☉</sub>)</th><th><i>R</i><sub>out</sub><br>(au)</th><th><i>M</i><sub>disc</sub><br>(<i>M</i><sub>SMS</sub>)</th><th>Flyby rate<br>(per year)</th></tr></thead><tbody><tr><td>SMS-mrd-hf</td><td>1000</td><td>500</td><td>0.01</td><td>1</td></tr><tr><td>SMS-mrd-mf</td><td>1000</td><td>500</td><td>0.01</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-mrd-lf</td><td>1000</td><td>500</td><td>0.01</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-mrD-hf</td><td>1000</td><td>500</td><td>0.1</td><td>1</td></tr><tr><td>SMS-mrD-mf</td><td>1000</td><td>500</td><td>0.1</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-mrD-lf</td><td>1000</td><td>500</td><td>0.1</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-mRd-hf</td><td>1000</td><td>1000</td><td>0.01</td><td>1</td></tr><tr><td>SMS-mRd-mf</td><td>1000</td><td>1000</td><td>0.01</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-mRd-lf</td><td>1000</td><td>1000</td><td>0.01</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-mRD-hf</td><td>1000</td><td>1000</td><td>0.1</td><td>1</td></tr><tr><td>SMS-mRD-mf</td><td>1000</td><td>1000</td><td>0.1</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-mRD-lf</td><td>1000</td><td>1000</td><td>0.1</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-Mrd-hf</td><td>10 000</td><td>500</td><td>0.01</td><td>1</td></tr><tr><td>SMS-Mrd-mf</td><td>10 000</td><td>500</td><td>0.01</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-Mrd-lf</td><td>10 000</td><td>500</td><td>0.01</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-MrD-hf</td><td>10 000</td><td>500</td><td>0.1</td><td>1</td></tr><tr><td>SMS-MrD-mf</td><td>10 000</td><td>500</td><td>0.1</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-MrD-lf</td><td>10 000</td><td>500</td><td>0.1</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-MRd-hf</td><td>10 000</td><td>1000</td><td>0.01</td><td>1</td></tr><tr><td>SMS-MRd-mf</td><td>10 000</td><td>1000</td><td>0.01</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-MRd-lf</td><td>10 000</td><td>1000</td><td>0.01</td><td>10<sup>-2</sup></td></tr><tr><td>SMS-MRD-hf</td><td>10 000</td><td>1000</td><td>0.1</td><td>1</td></tr><tr><td>SMS-MRD-mf</td><td>10 000</td><td>1000</td><td>0.1</td><td>10<sup>-1</sup></td></tr><tr><td>SMS-MRD-lf</td><td>10 000</td><td>1000</td><td>0.1</td><td>10<sup>-2</sup></td></tr></tbody></table>

potential inside $r_{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, minimizing the errors.

## 3.2 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

$$ \frac{\partial U}{\partial t} = - \nabla \cdot T(U) + S(U), \qquad (25) $$

where the state vector $U$ represents a set of conservative variables, $T(U)$ describes fluxes of each component of the state vector, and $S(U)$ represents the source terms.

For the case of ordinary hydrodynamics, equation (25) reduces to the following Euler equations:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0 \qquad (26) $$

$$ \frac{\partial \rho \mathbf{v}}{\partial t} + \nabla \cdot (\rho \mathbf{v} \mathbf{v} + p \mathbf{I}) = -\rho \nabla \Phi, \qquad (27) $$

where $\rho$ is the mass density, $\mathbf{v}$ is the velocity, $p$ is the thermal pressure, and the gravitational potential is $\Phi$.

The equation of state provides closure,

$$ p = \rho c_s^2, \qquad (28) $$

where the sound speed is a function of radius $c_s(r)$, as given by equation (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_{\odot}$ 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 (HLLC) discontinuity Riemann solver is used to solve the numerical fluxes. The radial grid is defined from 10 to 5000 au for $R_{\text{out}} = 500$ au and from 10 to 10 000 au for $R_{\text{out}} = 1000$ au using a logarithmic scale. The azimuthal, uniform grid is set from 0 to $2\pi$ with periodic boundary. We used a resolution of 50 × 100 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 yr, respectively. We run the simulations without any perturber for 5000 yr, 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 yr for simulations with $M_{\text{SMS}} = 10 000 M_{\odot}$, allowing the accretion disc to fully rotate at least 13 times.

## 3.3 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<sup>-1</sup>. 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 dy. 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<sup>-1</sup> (Kartje, Konigl & Elitzur 1999) and for the Galactic kilomaser W49N is 0.5–1.9 km s<sup>-1</sup> (Liljeström & Gwinn 2000). As an example we take 1 km s<sup>-1</sup>, 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<sup>-1</sup>, dy is multiplied by the number of these regions, $n$, 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),

$$ F = 4.7 \times 10^{17} a_{10} \left(\frac{l}{D}\right)^2 \text{Jy}, \qquad (29) $$
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Time evolution of the maser spectra for SMS with parameters same as Fig. 3 (left) and Fig. 5 (right). The unperturbed maser spectrum, i.e. at time 0 yr, is plotted green and at the end of the simulation run is plotted pink.</figcaption></figure>

where $a_{10} = a/10$, and $a$ defines the effective aspect ratio,

$$a = \frac{l}{\sqrt{A/\pi}}. \qquad (30)$$

For $a$ larger than a few the maser is saturated. We assumed $a_{10} = 1$. The half-length of the maser emission region that lies in a disc along the line of sight is defined as $l$, $D$ is the distance to the source with the observed area $A$. That results in equation (29), reducing to

$$F = 4.7 \times 10^{17} \left( \frac{ndy}{2D} \right)^2 Jy, \qquad (31)$$

where $D = 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 equation (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 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, McKee & Hollenbach (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.

# 4 RESULTS

## 4.1 Disc stability for a 1000 M$_{\odot}$ star

Fig. 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 5000 yr. 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 2500 yr (top row). Using the equation (1) with $\Delta v \sim 43$ km s$^{-1}$, gives a radius for the maser spots of $\sim 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, Fig. 3) between 60 and 75 km s$^{-1}$, giving the radius of the maser spots at around 380. It can also be noticed that the ‘high-velocity’ feature at around 66 km s$^{-1}$ at the time step of 10 000 yr 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 $\sim 350$ au, whilst the maser spots corresponding to the lower-flux feature at 153 km s$^{-1}$ are located at the radius of $\sim 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 Fig. 4 where we compare the unperturbed maser spectrum, i.e. at the very beginning of the simulation run, to the spectrum produced at 5000 yr.

## 4.2 Disc stability for a 10 000 M$_{\odot}$ star

Comparing the results for the case of a 10 000 M$_{\odot}$ SMS (simulation run SMS-Mrd-hf) in Fig. 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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same as Fig. 3 for run SMS-Mrd-hf, which has the same parameters as the simulation in Fig. 2 except the SMS mass is now 10 000 Msun. See caption for Fig. 1 for the definition of different features in the maser spectrum.</figcaption></figure>

in the first snapshot (top row in Fig. 5), symmetrically located at
around −13 and 240 km s⁻¹. Those values gives Δv ~ 126 km s⁻¹,
with equation (1) we get the radius of maser spot of ~ 550 au,
slightly inwards of the T = 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 Fig. 5) and the ‘high-velocity’
features in the corresponding maser spectrum move slightly outwards

from −13 to −26 km s⁻¹ and from 240 to 263 km s⁻¹, leading to
the maser spots being now located at the radius of 424 au therefore
slightly closer to the central star. At timestep of 3000 yr (bottom
row) the ‘high-velocity’ 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
Fig. 4.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison of density plots for simulations with different parameters. Each column represents different flyby rates, left: 1 per 100 yr (lf), middle: 1 per 10 yr (mf), right: 1 per 1 yr (hf). Simulation run in top row: SMS-mRD, second row: SMS-mRd, third row: SMS-MrD, and bottom row: SMS-MRD.</figcaption></figure>

## 4.3 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.

Fig. 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 (second 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 (third 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 (third 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).

## 5 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 on to 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 (Fig. 3) with the single, equal-mass, stellar flyby seen in Cuello et al. (2018, their Fig. 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, $c_s$, for our simulation model is about 1.0 km s⁻¹, whilst the corresponding Keplerian velocity, $v_k$, ranges between ~200 km s⁻¹ and ~42.0 km s⁻¹. Since $c_s \ll v_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 ‘high-velocity’ 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, $Q_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 equation (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 follow 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 summarizing the above findings, the model spectra from simulations for the $M_{SMS} = 1000 M☉$ are a better match with the W1 kilomaser spectrum from Gorski et al. (2019).

The H₂O molecule is dissociated close to the star where $T > 1500$ K. Those conditions do not allow for the water to stay molecular, hence the 22 GHz H₂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 centre, as the temperature increases closer to the centre and decreases with the distance approaching the outer radius. We used temperature profile of $T(r) = T_{in}(R/R_{in})^{-q}$, where $T_{in}$ indicates the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Examples of the maser spectrum evolution for a time interval of 4 yr for simulation run SMS-mrd-hf. Green lines represent the model maser spectrum at the beginning of the chosen time interval and the pink dashed lines show the spectrum at the end of the 4 yr time interval. Both time intervals have been chosen to present the biggest visible change in the spectrum during the simulation time. This is to show what changes could potentially be expected from repeated observations of the W1 kilomaser, presented in Gorski et al. (2019).</figcaption></figure>

effective temperature at the inner radius, $R_{\text{in}}$ and $q$ denotes the power-law index; in our simulations $q = 0.5$. Changing the value of $q$ affects the radius of the maser spots, e.g. $q = 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. Fig. 7 displays two samples of maser spectrum evolution within chosen periods of 4 yr 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$_{\odot}$. 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 $M_{\text{SMS}} = 1000$ M$_{\odot}$ show a smaller velocity spread between the peaks whilst for $M_{\text{SMS}} = 10\,000$ M$_{\odot}$ 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 $M_{\text{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, $R_{\text{out}} = 465$ au. The mean value of velocity spread between systemic and peak velocity of W1 kilomaser spectrum is $\Delta v = 83$ km s$^{-1}$. Using those values and the equation (1), we estimate that $M_{\text{SMS}} \approx 4000$ M$_{\odot}$ would best reproduce features of the W1 kilomaser.

# 6 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 Section 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 $M_{\text{SMS}} = 1000$ M$_{\odot}$ 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$_{\odot}$ 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$_{\odot}$ 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.

Summarizing 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☉.

## ACKNOWLEDGEMENTS

We thank the anonymous reviewer for useful comments and Mark Gieles, Corinne Charbonnel, and DJB Smith for helpful comments on an earlier version of the manuscript. This work was supported by Science and Technology Facilities Council (ST/V506709/1).

## DATA AVAILABILITY

The simulation data presented in this article will be shared on reasonable request to the corresponding author.

## REFERENCES

* Akiyama E., Momose M., Kitamura Y., Tsukagoshi T., Shimada S., Koyamatsu S., Hayashi M., 2013, PASJ, 65, 123

* Anderson J., 2002, in van Leeuwen F., Hughes J. D., Piotto G., eds, ASP Conf. Ser. Vol. 265, Omega Centauri, A Unique Window into Astrophysics. Astron. Soc. Pac., San Francisco. p. 87

* Bastian N., Cabrera-Ziri I., Salaris M., 2015, MNRAS, 449, 3333

* Bastian N., Lardo C., 2018, ARA&A, 56, 83

* Bedin L. R., Piotto G., Anderson J., Cassisi S., King I. R., Momany Y., Carraro G., 2004, ApJ, 605, L125

* Begelman M. C., 2010, MNRAS, 402, 673

* Billington S. J. et al., 2019, MNRAS, 490, 2779

* Brogan C., Johnson K., Darling J., 2010, ApJ, 716, L51

* Brunthaler A., Castangia P., Tarchi A., Henkel C., Reid M. J., Falcke H., Menten K. M., 2009, A&A, 497, 103

* Cesaroni R., 2005, Ap&SS, 295, 5

* Cesaroni R., Galli D., Lodato G., Walmsley C. M., Zhang Q., 2007, in Reipurth B., Jewitt D., Keil K., eds, Protostars and Planets V. The University of Arizona Press, Tucson. p. 197

* Chantereau W., Charbonnel C., Meynet G., 2016, A&A, 592, A111

* Charbonnel C., 2016, EAS Publ. Ser., 80, 177

* Clarke C. J., Pringle J. E., 1993, MNRAS, 261, 190

* Cuello N. et al., 2018, MNRAS, 483, 4114

* D'Antona F., Bellazzini M., Caloi V., Pecci F. F., Galleti S., Rood R. T., 2005, ApJ, 631, 868

* Darling J., Brogan C., Johnson K., 2008, ApJ, 685, L39

* Decressin T., Charbonnel C., Siess L., Palacios A., Meynet G., Georgy C., 2009, A&A, 505, 727

* Decressin T., Meynet G., Charbonnel C., Prantzos N., Ekström S., 2007, A&A, 464, 1029

* Denissenkov P. A., Denisenkova S. N., 1990, Sov. Astron. Lett., 16, 275

* Denissenkov P. A., Hartwick F. D. A., 2014, MNRAS, 437, L21

* Denissenkov P. A., Herwig F., 2003, ApJ, 590, L99

* Denissenkov P. A., VandenBerg D. A., Hartwick F. D. A., Herwig F., Weiss A., Paxton B., 2015, MNRAS, 448, 3314

* Doherty C. L., Gil-Pons P., Lau H. H. B., Lattanzio J. C., Siess L., 2014, MNRAS, 437, 195

* Dullemond C., Monnier J., 2010, ARA&A, 48, 205

* Elitzur M., McKee C. F., Hollenbach D. J., 1991, ApJ, 367, 333

* Ellingsen S. P., Voronkov M. A., Breen S. L., Caswell J. L., Sobolev A. M., 2018, MNRAS, 480, 4851

* Forestini M., Charbonnel C., 1997, A&AS, 123, 241

* Gao F. et al., 2017, ApJ, 834, 52

* Gieles M. et al., 2018, MNRAS, 478, 2461

* Gieles M., Charbonnel C., 2019, Proc. Int. Astron. Union, 14, 297

* Gorski M. D., Ott J., Rand R., Meier D. S., Momjian E., Schinnerer E., Ellingsen S. P., 2019, MNRAS, 483, 5434

* Gratton R. G., Carretta E., Bragaglia A., 2012, A&A Rev., 20, 50

* Greenhill L. J. et al., 2003, ApJ, 590, 162

* Greenhill L. J., Gwinn C. R., Antonucci R., Barvainis R., 1996, ApJ, 472, L21

* Gwinn C. R., Moran J. M., Reid M. J., 1992, ApJ, 393, 149

* Haemmerlé L., Meynet G., Mayer L., Klessen R. S., Woods T. E., Heger A., 2019, A&A, 632, L2

* Henkel C., Braatz J. A., Greenhill L. J., Wilson A. S., 2002, A&A, 394, L23

* Hofner P., Baan S. A., Takano S., 2006, AJ, 131, 2074

* Hollyhead K., Bastian N., Adamo A., Silva-Villa E., Dale J., Ryon J. E., Gazak Z., 2015, MNRAS, 449, 1106

* Ilee J. D., Cyganowski C. J., Brogan C. L., Hunter T. R., Forgan D. H., Haworth T. J., Clarke C. J., Harries T. J., 2018, ApJ, 869, L24

* Ishihara Y., Nakai N., Iyomoto N., Makishima K., Diamond P., Hall P., 2001, PASJ, 53, 215

* Karakas A. I., Fenner Y., Sills A., Campbell S. W., Lattanzio J. C., 2006, ApJ, 652, 1240

* Karakas A., Lattanzio J. C., 2007, PASA, 24, 103

* Kartje J. F., Konigl A., Elitzur M., 1999, ApJ, 513, 180

* Kimura S. S., Kunitomo M., Takahashi S. Z., 2016, MNRAS, 461, 2257

* Klahr H., Kley W., 2006, A&A, 445, 747

* Kondratko P. T., Greenhill L. J., Moran J. M., 2008, ApJ, 678, 87

* Krause M. et al., 2020, Space Sci. Rev., 216, 64

* Krause M., Charbonnel C., Decressin T., Meynet G., Prantzos N., 2013, A&A, 552, A121

* Kuo C. Y. et al., 2011, ApJ, 727, 20

* Kuo C. Y. et al., 2015, ApJ, 800, 26

* Langer G. E., Hoffman R., Sneden C., 1993, PASP, 105, 301

* Lardo C., Salaris M., Bastian N., Mucciarelli A., Dalessandro E., Cabrera-Ziri I., 2018, A&A, 616, A168

* Leroy A. K. et al., 2018, ApJ, 869, 126

* Lesniak M. V., Desch S. J., 2011, ApJ, 740, 118

* Levy R. C. et al., 2022, ApJ, 935, 19

* Liljestöm T., Gwinn C. R., 2000, ApJ, 534, 781

* Martins F., Chantereau W., Charbonnel C., 2021, A&A, 650, A162

* Martins F., Schaerer D., Haemmerlé L., Charbonnel C., 2020, A&A, 633, A9

* Martocchia S. et al., 2018, MNRAS, 473, 2688

* McGrath E. J., Goss W. M., De Pree C. G., 2004, ApJS, 155, 577

* Meire R., 1985, J. Br. Astron. Assoc., 95, 113

* Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228

* Milone A. P. et al., 2013, ApJ, 767, 120

* Milone A. P. et al., 2018, MNRAS, 481, 5098

* Milone A. P., Marino A. F., Piotto G., Bedin L. R., Anderson J., Aparicio A., Cassisi S., Rich R. M., 2012, ApJ, 745, 27

* Norris J. E., 2004, ApJ, 612, L25

* Osborn W., 1971, Observatory, 91, 223

* Pesce D. W., Braatz J. A., Condon J. J., Gao F., Henkel C., Litzinger E., Lo K. Y., Reid M. J., 2015, ApJ, 810, 65

* Piotto G. et al., 2007, ApJ, 661, L53

* Piotto G. et al., 2015, AJ, 149, 91

* Prantzos N., Charbonnel C., Iliadis C., 2007, A&A, 470, 179

* Prantzos N., Charbonnel C., Iliadis C., 2017, A&A, 608, A28

* Qiao L., Haworth T. J., Sellek A. D., Ali A. A., 2022, MNRAS, 512, 3788

* Reid M. J., Braatz J. A., Condon J. J., Greenhill L. J., Henkel C., Lo K. Y., 2009, ApJ, 695, 287

* Rekola R., Richer M. G., McCall M. L., Valtonen M. J., Kotilainen J. K., Flynn C., 2005, MNRAS, 361, 330
Renzini A., Marino A. F., Milone A. P., 2022, MNRAS, 513, 2111

Sandage A., Wildey R., 1967, ApJ, 150, 469

Sanna A. et al., 2021, A&A, 655, A72

Sato N., Yamauchi A., Ishihara Y., Sorai K., Kuno N., Nakai N., Balasubramanyam R., Hall P., 2005, PASJ, 57, 587

Schleicher D. R. G., Palla F., Ferrara A., Galli D., Latif M., 2013, A&A, 558, A59

Shakhvorostova N., Sobolev A., Moran J., Alakoz A., Imai H., Avdeev V., 2019, Adv. Space Res., 65, 772

Siess L., 2010, A&A, 512, A10

Urquhart J. S. et al., 2011, MNRAS, 418, 1689

Ventura P., D'Antona F., Mazzitelli I., Gratton R., 2001, ApJ, 550, L65

Ventura P., Di Criscienzo M., Carini R., D'Antona F., 2013, MNRAS, 431, 3642

Volvach A. E., Volvach L. N., Larionov M. G., 2020, MNRAS, 496, L147

Vorobyov E. I., Skliarevskii A. M., Elbakyan V. G., Takami M., Liu H. B., Liu S.-Y., Akiyama E., 2020, A&A, 635, A196

Vorobyov E. I., Steinrueck M. E., Elbakyan V., Guedel M., 2017, A&A, 608, A107

Vural J. et al., 2014, A&A, 564, A118

Williams G. M. et al., 2022, MNRAS, 509, 748

Yamauchi A., Nakai N., Sato N., Diamond P., 2004, PASJ, 56, 605

Yong D., Grundahl F., Norris J. E., 2014, MNRAS, 446, 3319

Zhang B., Reid M. J., Menten K. M., Zheng X. W., Brunthaler A., Dame T. M., Xu Y., 2013, ApJ, 775, 79

Zhao W. et al., 2018, ApJ, 854, 124

This paper has been typeset from a *TEX/LATEX* file prepared by the author.