Issue |
A&A
Volume 693, January 2025
|
|
---|---|---|
Article Number | A180 | |
Number of page(s) | 19 | |
Section | Interstellar and circumstellar matter | |
DOI | https://doi.org/10.1051/0004-6361/202452800 | |
Published online | 15 January 2025 |
The formation and stability of a cold disc made out of stellar winds in the Galactic centre
1
Hamburger Sternwarte, Universität Hamburg,
Gojenbergsweg 112,
21029
Hamburg,
Germany
2
Max-Planck-Institut für Astrophysik,
Karl-Schwarzschild-Straße 1,
85748
Garching,
Germany
3
Departamento de Ciencias, Facultad de Artes Liberales, Universidad Adolfo Ibáñez,
Av. Padre Hurtado 750,
Viña del Mar,
Chile
4
Millennium Nucleus on Transversal Research and Technology to Explore Supermassive Black Holes (TITANS),
Chile
5
Department of Physics and Astronomy, Bartol Research Institute, University of Delaware,
Newark,
DE
19716,
USA.
6
Universitäts-Sternwarte, Ludwig-Maximilians-Universität München,
Scheinerstr. 1,
81679
Munich,
Germany
7
Max-Planck Institute for Extraterrestrial Physics,
Giessenbacherstr. 1,
85748
Garching,
Germany
8
Excellence Cluster ORIGINS,
Boltzmannstrasse 2,
85748
Garching,
Germany
9
The Oskar Klein Centre, Department of Astronomy, AlbaNova, Stockholm University,
106 91
Stockholm,
Sweden
10
Department of Astronomy, University of Michigan,
1085 S. University,
Ann Arbor,
MI
48109,
USA
★★ Corresponding author; calderon@mpa-garching.mpg.de
Received:
29
October
2024
Accepted:
5
December
2024
Context. The reported discovery of a cold (~104 K) disc-like structure within the central 5 × 10−3 pc around the super-massive black hole at the centre of the Milk Way, Sagittarius A* (Sgr A*), has challenged our understanding of the gas dynamics and thermodynamic state of the plasma in its immediate vicinity. State-of-the-art simulations do not agree on whether or not such a disc can indeed be a product of the multiple stellar wind interactions of the mass-losing stars in the region.
Aims. The aims of this study are to constrain the conditions for the formation of a cold disc as a natural outcome of the system of the mass-losing stars orbiting around Sgr A*, to investigate whether the disc is a transient or long-lasting structure, and to assess the validity of the model through direct comparisons with observations.
Methods. We performed a set of hydrodynamic simulations of the observed Wolf-Rayet (WR) stars feeding Sgr A* using the finite- volume adaptive mesh refinement code Ramses. We focus, for the first time, on the impact of the chemical composition of the plasma emanating from the WR stars.
Results. The simulations show that the chemical composition of the plasma affects the radiative cooling to a sufficient degree to impact the properties of the medium, such as density and temperature, and, as a consequence, the rate at which the material inflows onto Sgr A*. We demonstrate that the formation of a cold disc from the stellar winds is possible for certain chemical compositions that are consistent with the current observational constraints. However, even in such cases, it is not possible to reproduce the reported properties of the observed disc-like structure, namely its inclination and the fluxes of its hydrogen recombination lines.
Conclusions. We conclude that the stellar winds alone are not sufficient to form the cold disc around Sgr A* inferred from observations. Either relevant ingredients are still missing in the model, or the interpretation of the observed data needs to be revised.
Key words: accretion, accretion disks / hydrodynamics / stars: winds, outflows / stars: Wolf-Rayet / Galaxy: center
© The Authors 2025
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model.
Open Access funding provided by Max Planck Society.
1 Introduction
The Galactic centre (GC) hosts the closest super-massive black hole to us, Sagittarius A* (Sgr A*; see Genzel et al. 2010, for a review). Unlike the black holes present in active galactic nuclei, Sgr A* is highly underluminous. It does not seem to be currently accreting a significant amount of material or to have a standard accretion disc around it (Yuan & Narayan 2014). Nevertheless, Murchikova et al. (2019) detected a disc-like structure around Sgr A*. Using the 1.3 millimetre recombination line H30α, these authors observed a double-peaked emission line with a full velocity width of ~2200 km s−1. The centre of the emission coincides with Sgr A* and extends up to 0.11″ (or ~104 Schwarzschild radii) to both the redshifted and blueshifted sides. Murchikova et al. (2019) interpreted this feature as a rotating disc with a mass of 10−5–10−4 M⊙. Yusef-Zadeh et al. (2020) also reported the presence of broad hydrogen recombination lines – including the double-peaked H30α line – at the position of Sgr A* using the Very Large Array. However, these authors interpreted the signatures as a jet emanating from the central region. So far, there has been no detection in the near-infrared, despite many observations. Ciurlo et al. (2021) report an upper limit for Brγ that is two orders of magnitude below the extrapolation from the reported H30α flux. At present, there is no consensus on how such observations can be reconciled or from where the observed cold material might come.
The gaseous environment around Sgr A* is dominated by the outflows from around 30 Wolf-Rayet (WR) stars, all located within a fraction of a parsec of the black hole (Paumard et al. 2006; Martins et al. 2007). At such close distances, the winds interact strongly in shocks that thermalise them and create a hot and diffuse X-ray-emitting plasma (Quataert 2004). However, given the relatively low velocity of some of the outflows (450–600 km s−1; Martins et al. 2007), the resulting plasma can be prone to radiative cooling and form dense, cold clumps (Cuadra et al. 2005; Calderón et al. 2016), which end up being embedded in the diffuse, hot plasma (~107 K; Baganoff et al. 2003). The complex interplay between the stellar winds around Sgr A* has been the subject of hydrodynamic simulations using a variety of numerical techniques. With smoothed particle hydrodynamics (SPH) simulations, Cuadra et al. (2006) showed that, if most of the slow-wind stars are located in a relatively compact stellar disc, many cold clumps form and quickly coalesce in a cold gaseous disc with a radius of ~1″. However, using a stellar distribution closer to the one observed, Cuadra et al. (2008, 2015) showed that no conspicuous disc appears over the 1800 yr time-span of their models. Calderón et al. (2020b) performed finite-volume hydrodynamic simulations on a Cartesian grid of the same system – which are better suited to modelling shocks, the multi-phase medium, and subsequent cooling – and found that the formation of a cold disc is indeed possible. According to their model, the stellar wind of the star IRS 33E, which is relatively dense (~10−5 M⊙ yr−1) and slow (450 km s−1; Martins et al. 2007), interacts with the medium, creating a shell that is dense enough to quickly radiate its thermal energy and break into denser and smaller structures. These pieces manage to fall inwards, forming a cold (~ 104 K) disc with a total mass of ~5× 10−3 M⊙ within a simulation time of 3500 yr. Based on this result, Ciurlo et al. (2021) found that the properties of the modelled disc are consistent with the non-detection in Brγ. Nevertheless, this scenario has not been confirmed by analogous grid-based simulations. Ressler et al. (2020) revisited the same system using the same numerical approach, although with a different code: ATHENA++ (Stone et al. 2020), and found no disc formation, even when extending their simulation time up to 9000 yr. Building on this model, Solanki et al. (2023) also studied this system through hydrodynamic modelling, but encompassing a larger region and evolving it for much longer timescales. These authors included the presence of the cir- cumnuclear disc (CND), an observed gaseous structure located between 1.5 and 3.0 pc (Becklin et al. 1982; Genzel 1989) and with a total mass of 3–4 × 104 M⊙ (Dinh et al. 2021; James et al. 2021). As a result, they showed that the formation of a cold disc is possible on long timescales (≳300 000 yr) due to the interaction of the innermost boundary of the CND and the WR stellar winds, which results in the transport of CND material to smaller scales. However, the authors warned that the lack of certain physical mechanisms in the model on such timescales (e.g. magnetic fields, supernovae, and thermal conduction) might affect the robustness of this result. Thus, the exact formation process of the observed cold disc remains unknown.
A key factor driving the formation of clumps and a disc in the model of Calderón et al. (2020b) is that radiative cooling allows the gas to get rid of its energy efficiently. The strength of this process depends on the chemical composition of the gas, which is highly unusual in the GC given its origin as winds from evolved massive stars. This work explores, for the first time, the impact of radiative cooling through studying specific chemical compositions: solar with 1 Z⊙, 3 Z⊙, and 5 Z⊙. More importantly, we developed a more realistic chemical setup that considers the atmospheric abundances for the different WR subtypes present in the region in order to follow the thermodynamic evolution of the gas in a more appropriate manner. Our results show that the formation of the disc is indeed determined by the composition of the gas. However, the properties of this disc do not agree with the observations when considering realistic wind compositions compatible with current observational constraints. We also improved the numerical setup in order to better assess the disc orientation and more directly compare our results with those of Ressler et al. (2018, 2020) as well as with SPH models (Cuadra et al. 2008, 2015; Russell et al. 2017).
This article is organised as follows: in Section 2, we present the numerical approach, the setup, and the models investigated. Section 3 presents and describes the results of the numerical simulations. In Section 4, we present the synthetic observables obtained through post-processing the models and contrast them with observations. We compare our grid-based models with SPH models in Section 5. In Section 6, we present an analytic analysis of the stability of the observed cold disc. In Section 7, we discuss our findings and uncertainties in the parameters of our models. Finally, in Section 8, we present conclusions and final remarks. Throughout this paper we use the mass and distance to Sgr A* of 4.3 × 106 M⊙ and 8.33 kpc, respectively (Gillessen et al. 2017; GRAVITY Collaboration 2019, 2021), so that 1″ corresponds to a length of ~0.04 pc ≈105 RSch, where RSch is the Schwarzschild radius.
2 Numerical simulations
2.1 Equations
The simulations were performed using the adaptive-mesh refinement (AMR) hydrodynamic code Ramses (Teyssier 2002). This code uses a second-order Godunov method with a shockcapturing scheme to solve the Euler equations in their conservative form, that is,
(1)
(2)
(3)
(4)
where (ρ, u, P) are the primitive hydrodynamic variables: density, velocity, and pressure, respectively. The set of quantities si correspond to tracer scalar fields that are advected with the fluid whose usage we introduce in Section 2.3. The sink terms on the right-hand side correspond to the effects of the time-independent gravitational potential ϕ = ϕ(x) (in our case solely of Sgr A*; see Section 2.2), and the total radiative losses due to optically thin radiative cooling , with x the position vector, T the fluid temperature, and Xi the chemical composition of the fluid. The total specific energy density e is given by
(5)
where γ is the adiabatic index that is set to 5/3. Additionally, we consider the fluid can be described as an ideal gas so that the temperature can be calculated through P = (ρ/µ)kB T, being µ the mean molecular weight and kB the Boltzmann constant.
2.2 Numerical setup
The model considers the system of WR stars blowing stellar winds while they move on their observed Keplerian orbits around Sgr A*. The setup is analogous to our previous work (Calderón et al. 2020b). The central black hole gravitational field is modelled as a point mass of 4.3 × 106 M⊙ (GRAVITY Collaboration 2019, 2021). The stars are simulated as test particles that only feel the gravitational pull of Sgr A*. Their initial position and velocity vectors are set so that they move on the orbits constrained by observations (Paumard et al. 2006; Cuadra et al. 2008; Gillessen et al. 2019; von Fellenberg et al. 2022). The stellar winds in the simulation are modelled following the procedure by Lemaster et al. (2007), which has been validated and used in our previous work (Calderón et al. 2020a,b). This consists in defining a “masked region” around each star that is a spherical volume where the hydrodynamic variables are reset in each timestep, in order to reproduce the free wind expansion solution of a spherical wind with certain mass-loss rate Ṁw, terminal velocity Vw, and temperature Tw. For these quantities we used the values constrained through modelling of the infrared spectra (Martins et al. 2007; Cuadra et al. 2008). The wind temperature was set to the lowest temperature allowed in the simulation, Tw = 104 K which is determined by the strong ultraviolet radiation produced by the hundreds of massive stars in the region.
The simulations were run in a Cartesian grid in a cubic domain of side 40″ (~1.6 pc) with outflow boundary conditions (zero gradients). We used the exact Riemann solver (e.g. Toro 2009) combined with the MinMod slope limiter, as these choices allow the modelling of hydrodynamic instabilities, which otherwise may be quenched due to numerical diffusion. We investigated other choices such as the Harten-Lax-van LeerContact (HLLC; Toro et al. 1994) and/or the MonCen slope limiter. Our tests showed that the use of the HLLC solver amplified the grid artifacts while the use of the MonCen limiter resulted in unwanted numerical artifacts such as spurious oscillations. For a careful analysis of these choices we refer the reader to the Appendix A. Regarding the numerical resolution of our models, the coarse resolution was 643 cells, allowing adaptive refinement of four levels, that is, effectively Δx ≈ 0.0391″ ≈ 1.56 × 10−3 pc. However, the regions around the stars allowed two extra levels of refinement (Δx ≈ 9.77 mas ≈ 3.91 × 10 pc) in order to guarantee at least eight cells as the radius of the masked regions where the winds are initialised. These regions have a fixed radius of 80 mas ≈ 3.2 × 10−3 pc. Additionally, the vicinity around the central black hole has a fixed nested grid with eight refinement levels above the coarse resolution (Δx ≈ 2.44 mas ≈ 9.77 × 10−5 pc). At the location of Sgr A*, we defined a spherical region where the hydrodynamic variables are reset to low values of density and pressure at rest, in order to avoid artificial accumulation of material. The refinement strategy is based on density gradients on top of the geometric criteria previously defined. In order to explore the role of the AMR potentially affecting the results we tested different values of the smoothing parameter that controls the refinement in transitions between refinement levels. We tested quadrupling the smoothing parameter and the results remained unchanged1.
2.3 Models
The main parameter explored in this work is the optically thin radiative cooling function. This choice is determined by specifying the chemical composition of the fluid. Unfortunately, the metallicity of the young stars in the GC is still poorly known, yet it has been argued that it should be higher than Solar (Z = 2-3 Z⊙; Genzel et al. 2010). Past numerical works have considered that the composition of the gas correspond to Solar abundances but with metallicity three times the Solar value, Z = 3 Z⊙ (Cuadra et al. 2008; Calderón et al. 2016; Ressler et al. 2018, 2020; Solanki et al. 2023). However, more than the ‘bulk’ metallicity of the stars, the most appropriate composition choice would be the abundances of the WR stellar atmospheres, which normally differ significantly from Solar (Herald et al. 2001; Onifer et al. 2008). Although Ressler et al. (2018, 2020) used chemical compositions lacking hydrogen, the rest of the element abundances remained unchanged relative to Solar with 3 Z⊙. As a result, the cooling function varied at low temperature (<105 K) but at higher temperatures it remains unchanged. Moreover, the composition choice not only determines the cooling function but also the ion mean molecular weight as well as the electron to proton number densities ne/np. These quantities are taken into account in the sink term in the energy equation (see Eq. (3)), and following Schure et al. (2009) can be expressed as follows
(6)
where the subscript i stands for a given chemical abundance. Bear in mind that different chemical compositions also correspond to different temperature Ti , as this is obtained through the ideal gas expression that makes uses of the mean molecular weight of the fluid. In the general case, if we consider a fluid element composed of N mixtures of abundances, the total radiative losses will be given as a summation, that is,
(7)
where we have expressed the total radative cooling as a linear combination of the mixtures. Notice that we have introduced the passive scalar fields si , as we used them to quantify the fraction of a given chemical composition. In this work, we introduced three types of compositions that represent the WR subtypes based on their atmosphere abundances that we assume to be constant over the few thousand year duration of the simulations. In the region where the winds are generated, the corresponding passive scalar is set to 1 while the rest are set to 0. By doing so, it is possible to identify how much material of a given cell is supplied by which subgroup of WR stars. The mass fractions of each mixture are shown in Table 1.
In this work, we investigated five models with different compositions: solar composition with metallicities Z⊙, 3 Z⊙, and 5 Z⊙, a mixture of three compositions based on the spectroscopic constraints of the subtypes of WR stars, and a variation of the latter one motivated by the uncertainty in the WR stellar atmospheres and, specifically on the H fraction in them that is set to XH = 40% (see Section 7.2, for a discussion). Figure 1 shows the cooling function Λ = Λ(T) as a function of temperature for different chemical mixtures. Each curve was calculated as a linear combination of the abundance of a given element and its contribution to the total radiative cooling. The values of the composition per element are shown in Table 1. The contribution of each element to the total cooling was taken from the plasma models developed by Štofanová et al. (2021). The cooling functions corresponding to WR subtypes compositions: WC8-9, WN5-7, and WN8-9/Ofpe are shown with dashed blue, dotted orange, and dotted-dashed green lines, respectively. The Solar and Solar with 3 Z⊙ are represented with thin and thick solid black lines, respectively. Furthermore, we added the cooling function used in previous works by Ressler et al. (2018, 2020) as a dotted-dashed red line. Here it can be observed that a Solar composition with 3 Z⊙ (the typical used value) is between the range of the different WR compositions but the subtypes WC89 and WN89/Ofpe are about a factor two higher.
The list of the runs investigated in this work are shown in Table 2. The initial setup was identical to our previous work (Calderón et al. 2020b), i.e. the medium across the whole domain was set to constant and low-enough values of density ρ = 10−24 g cm−3 and pressure with cs,f = 10 km s−1, so that the stellar winds do not encounter impediments and fill quickly the domain. All models were run for a total of 3500 yr but setting the starting state of the system in the past, so that the final state of the simulations corresponds to the current state of the system. This is achieved by integrating the current position and velocity vectors of the stars back in time, assuming that they have followed purely Keplerian orbits (Paumard et al. 2006; von Fellenberg et al. 2022) within this timescale due to the gravitational field of Sgr A* (Cuadra et al. 2008; Ressler et al. 2018; Calderón et al. 2020b).
Atmospheric mass fractions (in percentage) and mean molecular weights.
![]() |
Fig. 1 Comparison of cooling functions Λ(T) for different chemical abundances. The thin and thick lines show the radiative cooling for Solar and three times Solar abundances, respectively. The dashed blue, dotted orange, and dot-dashed green lines represent the cooling functions for atmosphere abundances corresponding to WR subtypes: WC8-9, WN5-7, WN8-9 and Ofpe/WN9 based on the compositions compiled shown in Table 1. The contributions of each element to the radiative cooling were taken from the plasma models by Štofanová et al. (2021). Additionally, the cooling function used in Ressler et al. (2018) is shown as a thick black line. |
Simulation runs and parameters.
3 Results
We proceed to describe the evolution and final state of the simulations at t = 0 that corresponds to the present time. The runs with Solar composition, but varying metallicities, were used for understanding the impact of increasing and decreasing the effect of radiative cooling in general. However, since they do not represent realistic compositions we opted not to discuss them in detail here. Instead, we chose to focus on the more physically motivated models, i.e. WR_f07 and WR_f1 which we also refer as fiducial and enhanced, respectively. For completeness, we give a brief description of the models A1, A2, and A3 in Appendix B.
3.1 Hydrodynamics
The simulation evolution of all models is analogous, especially in the initial phases. At t = −3500 yr, the stars begin to orbit around Sgr A* that is located at the centre of the domain, while their winds quickly fill the whole computational domain. After ~500 yr (t = −3000 yr), the system reaches a quasi-steady state. At this point, the domain is full of diffuse (∼10−22 –10−21 g cm−3) and hot plasma (∼107–108 K) due to the shocked stellar winds and their interactions. This is illustrated in Figure 2 that shows density maps across different simulation times. The maps show projected density fields along the line-of-sight direction (z axis) weighted by density, i.e. ∫ ρ2dz/ ∫ρdz, in order to highlight the densest regions with the highest resolution. The top panels show two models at t = −700 yr, respectively. In them, it is possible to observe the complex structure developed due to the stellar wind interactions such as bow shocks, instabilities, and dense clumps as a result of condensation through radiative cooling. Left- and right-hand side panels display models with different chemical compositions WR_f07 and WR_f1, respectively. On the left-hand side, it can be seen how mildy efficient cooling affects mostly one of the stellar winds. This star corresponds to IRS 33E and its wind is the slowest (∼450 km s−1). This fact causes its shocked temperature to be in an efficient region of the radiative cooling functions that allows its material to cool down and become denser. Also, some of its material manages to reach the vicinity of Sgr A*, as the map shows an elongated clump being accreted and stretched. The right-hand side also portrays this picture but since the radiative cooling is enhanced more dense clumps can be observed, especially close to IRS 33E and Sgr A*. Finally, the bottom panels of Fig. 2 show the models at time t = 0, i.e. the present time. One can clearly see that the amount of dense material that has accumulated around Sgr A* is different in both maps. Although in the model WR_f07 (see bottom left-hand side panel of Fig. 2) some dense material is spiraling towards the black hole overall there are is no clear structure around it. In the case of WR_f1 (see bottom right-hand side panel of Fig. 2), the dense material has settled at the centre in a a disc-like structure. This difference is due to the efficiency of the radiative cooling. Since model WR_f1 has enhanced cooling, more dense clumps and filaments form, and some of them manage to fall onto Sgr A*. Overall, the hydrodynamic evolution of the two models is analogous: the winds fill the domain during the first ~500 yr, then the systems reach a quasi-steady state, and in the last 500 yr they diverge as the model WR_f1 creates much more dense, cool material that settles around Sgr A*. A quantitative analysis of the evolution of the properties of system is presented as follows.
To quantify the accretion rate at different spatial scales we calculated the mass flux as a function of both radial distance and time Ṁ(r, t) = 4πr2ρ(r, t)vr(r, t) averaged over a spherical shell. Here, positive and negatives signs in the radial velocity refer to outflow and inflow mass fluxes, respectively. First, we analyse the net mass flux at the innermost radius we can resolve properly, which is equivalent to five times the inner boundary radius, i.e. 5rin = 5 × 10−4 pc ≈ 1.25 × 10−2″. Figure 3 shows |Ṁ(r = 5rin, t)| as a function of time for models WR_f07 and WR_f1 that are displayed as solid blue and orange lines, respectively. In both cases, at t > −3300 yr the net mass flux reaches levels of ~10−6 M⊙ yr−1 with short variability episodes that increase the rate by factors of two to four. During this quasi-steady phase, both simulations display similar behaviour qualitatively, likely determined by the identical stellar wind configuration. This stage lasts until t ≈ −500 yr where the mass-flow rates start to deviate from each other. The fiducial model continues to display a variability amplitude of the same order of magnitude. However, the WR_f1 model shows a transition to a strongly gas inflow dominated phase, with inflow rates that are enhanced by a factor four to eight. This stage of the evolution is what we refer to as the disc formation phase, analogously to our previous models reported in Calderón et al. (2020b). The spatial distribution of the inflowing material is not simple. The fiducial case shows the same behaviours described in previous work (e.g. Ressler et al. 2018) where there are two components: one that most of the time follows the same orientation of the stars in the clockwise disc, and other that comes from the south pole of such a disc. In the enhanced case, the behaviour is analogous but the inflow component along the disc is more relevant. In both cases, the inflowing material exhibits a significant range in angular momentum orientation as discussed further in this section. Thus, none of them shows a simple thin disc accretion scenario but rather a thick disc with an extra component along the south pole.
Next, we proceed to analyse the time-averaged behaviour of the simulations over the last 500 yr, as this can give us an idea of the general state of the system minimising the effects of the stochastic variability. First, we analyse the (volume-weighted) density and (mass-weighted) temperature radial profiles of the simulations, which are shown on the top and bottom panels in Fig. 4, respectively. The dashed blue and solid orange lines represent the simulations WR_f07 and WR_f1, respectively. The orange lines in both panels clearly highlight the presence of the cold disc. The density profile decays with r2 in both cases at large scales (≳1″). This is the result of most of the material flowing outwards at these scales following a roughly isotropic spherical wind as seen in previous models (Ressler et al. 2018; Calderón et al. 2020b). However, the profiles differ at smaller scales (<1″) where the fiducial case transitions to ρ ∝ r−1, while the model WR_f1 shows a density enhancement due to the presence of the disc. This increase in density is more than one order of magnitude. Regarding the temperature profiles, both models match at larger scales (≳1″) with a constant temperature of the order of 107 K, set by the stellar wind collisions. Again, the profiles differ at smaller scales where the fiducial case follows T ∝ r− 1, and the enhanced cooling case displays a temperature profile about two orders of magnitudes lower on average. The minimum in temperature corresponds to the region where the disc contributes with most of the mass for a given spherical shell.
To quantify and characterise the mass inflow and outflow regimes in the simulation domain we calculated them as radial profiles. Figure 5 displays the absolute value of the mass flow rates as a function of distance from Sgr A*. The solid and dashed lines represent the mass inflow and outflow rates, respectively; while the colours follow the same convention as before. Thus, if the solid line is above the dashed line the net mass flow is inwards and vice-versa. Analogous to the model by Ressler et al. (2018), the fiducial model encompasses different regimes due to the dominant direction of the mass flow at given spatial scales. At r > 3″, the outflow component dominates and the inflow is negligible. We find a net mass outflow rate of ~5 × 10−4 M⊙ yr−1. Notice that the enhanced cooling model exhibits exactly the same behaviour at these scales. At smaller scales (1″ – 3″), there is a transition region where |Ṁin| ~ |Ṁout| that corresponds to the location of most mass-losing stars and wind interactions. For model WR_f1, the change is more abrupt due to the presence of the cold disc. In both models, at r < 1″ the inflow component is larger than the outflow, so the material inflows across these scales and down to the innermost boundary. The net effect makes the mass inflow rate about five times larger which generates the cold disc.
The origin of the infalling material can be traced in two ways: analysing the angular momentum of the material close to the inner boundary, and using the scalar tracer fields that we introduced to label the chemical abundances of the different WR subtypes. Figure 6 shows Hammer projections of the angular momentum direction of the gas enclosed in a sphere of radius 0.25″ (~0.01 pc) for the fiducial and enhanced cooling runs in the top and bottom panels, respectively. Each point represents a different simulation time over the last 1000 yr, and they are connected with dashed lines. For reference, we added the angular momentum direction of the orbits of the mass-losing stars as black star symbols, and the location of the clockwise disc based on the orbits we used as a grey shaded circle. This analysis shows that the angular momentum direction of the infalling material varies stochastically but overall tends to align with the orientation of the orbits in the clockwise disc. However, the degree of variability depends on whether or not the model results in the formation of a disc. The fiducial model shows more variability while the enhanced cooling run displays that the angular momentum aligns more consistently with the orientation of the clockwise disc. The fact that the infalling material at this spatial scales comes from these stars agrees with previous numerical models (Cuadra et al. 2008; Ressler et al. 2018; Calderón et al. 2020b). Furthermore, this is not surprising since these stars orbit closer to Sgr A*, and have winds that are relatively slow (~600 km s), which results into shocks more prone to radiative cooling and with smaller angular momentum. It is also relevant to mention that the average angular momentum direction of the accreted material has a significant dispersion, despite being similar to the direction of the angular momentum of the stars in the clockwise disc. Although not shown here, the angular momentum direction distribution shows a polar thickness of ~40° for the bulk of the inflow material, and in some cases ~100°. This means that the inflow takes place preferentially along the orientation of the clockwise disc (or the cold disc) but also along higher inclinations close to the disc orientation. Such a behaviour can be seen in both the fiducial and enhanced models, and are consistent with previous work results (e.g. Ressler et al. 2018).
Figure 7 shows radial profiles of the scalar tracers. The dashed, solid, and solid-dashed represent the scalars that correspond to the WR subtypes WC89, WN57, and WN89/Ofpe, respectively. The blue and orange lines stand for the results for the models WR_f07 and WR_f1. In this analysis, the scalar tracer values represent the fraction of the mass density that comes from a given WR subtype wind. Thus, in general most of the material supplied into the domain comes from the WN89/Ofpe stars. At smaller scales (r < 1″), the dominance is even clearer as the material from these stars makes about 80–90% of the mass budget. This is also the result of these stars having slow winds and colder shocked material. The rest of the stars contribute mostly to the outflow of material, as their contribution is more relevant at larger scales. Notice that both models, fiducial and enhanced, display roughly the same behaviour across the entire domain.
Overall, the fiducial model displays properties of the medium and gas dynamics consistent with the previous models that do not form a disc (Cuadra et al. 2008; Ressler et al. 2018, 2020; Calderón et al. 2020b, run for 1100 yr). The model with enhanced cooling is consistent with the simulations that show the formation of a cold disc as an outcome of this system (Calderón et al. 2020b). Before discussing if any of the models can be favoured given the current observational constraints we proceed to characterise the cold disc arising in the WR_f1, as its properties will aid us to do so.
![]() |
Fig. 2 Comparison of the simulations with different chemical abundances at two different simulation times. The panels show projected density maps weighed by density along the z axis, i.e. ∫ ρ2dz/ ∫ ρdz, which is parallel to the line of sight. Top and bottom panels show the systems at t = −0.7 kyr and t = 0 (present), respectively. Left- and right-hand side panels show the runs WR_f07 and WR_f1, respectively. All maps display the full computational domain. |
![]() |
Fig. 3 Net mass-flow rate across a sphere of radius 5rin = 5 × 10−4 pc (1.25 × 10−2″). At this radius, the direction of the mass flow is inwards throughout the entire simulation. The models WR_f07 and WR_f1 are shown in solid blue and orange lines, respectively. |
![]() |
Fig. 4 Radial profiles of time-averaged volume-weighted density (top) and mass-weighted temperature (bottom) over the last 500 yr of simulation time. The models WR_f07 and WR_f1 are shown in solid orange and dashed blue lines, respectively. |
![]() |
Fig. 5 Radial profiles of the mass inflow Ṁin (solid lines) and outflow Ṁout (dashed lines) rates averaged over the last 500 yr of the simulation. The models WR_f07 and WR_f1 are shown in solid orange and blue lines, respectively. |
![]() |
Fig. 6 Orientation of the angular momentum of the gas enclosed in a sphere of radius 0.01 pc (~0.25″). The vertical dimension represents inclination i, while the horizontal dimensions stand for the longitude of the ascending node Ω. Thus, a face-on star orbiting clockwise on the sky would be at the north pole of the graph. Top and bottom panels show the runs WR_f07 and WR_f1, respectively. Each dot corresponds to the analysis of a single snapshot. As a reference, the orientation of the orbital angular momentum of the mass-losing stars is shown as black star symbols, and the grey shaded region corresponds to the average direction of the clockwise disc at (104°, 126°) with a dispersion of 16° (Yelda et al. 2014). |
![]() |
Fig. 7 Radial profiles of the scalar tracer fields that represent the WR subtypes: WC89 (dashed lines), WN57 (solid lines), and WN89/Ofpe (dot-dashed) averaged over the last 500 yr of the simulation. The models WR_f07 and WR_f1 are shown in solid orange and blue lines, respectively. |
3.2 Properties of the cold disc
At the present time (t = 0), the model WR_f1 shows the presence of a cold disc. Figure 8 shows zoomed maps of the central region of 2.5″× 2.5″. The left- and right-hand side panels show projected maps of density and line-of-sight velocity (both weighted by density), respectively. Here it is possible to observe that the projected diameter of the disc is roughly 1″. We estimated the mass of the disc by isolating the cold material (T < 105 K) and integrating the cell mass content within a sphere of radius 1.0″. We found that the total mass of the disc is 0.005 M⊙ which is consistent with the value in our previous work (Calderón et al. 2020b). However, the mass accumulated in the disc does not converge and keeps increasing at least for the next hundreds of years in the simulation.
On the right-hand side panel of Fig. 8, we can see that the line-of-sight velocity peaks at ~2000 km s−1 along both directions (towards and outwards the observed). This coincides with the maximum extension of the observed H30α line (Murchikova et al. 2019). Additionally, the disc is observed to be ~45° tilted in projection. This is not in agreement with the observations since, at least in projection the difference is about ~90° (see Figure 1 of Murchikova et al. 2019).
4 Post-processing of observables
In order to assess the validity of our models for the Galactic centre we proceed to synthesise observational quantities that we could compare easily with observations. We focus on the recombination lines H30α and Brγ. Afterwards, we calculate the X-ray spectrum from our models.
4.1 Recombination lines: H30α and Brγ
The cold gas at the location of Sgr A* has been observed through the H30α recombination line (Murchikova et al. 2019; Yusef-Zadeh et al. 2020), and only with an upper limit in the Brγ recombination line (Ciurlo et al. 2021). In order to explore if the models are consistent with these observations we have calculated the expected flux from these emission lines. To compute the emission of the H30α line we used the coefficients given below fitting a piecewise linear function to them as a function of density and considering a decay ∝ T−1 (see Murchikova et al. 2019, supplementary information),
(8)
(9)
(10)
(11)
In the case of the Brγ line, its emissivity can be calculated following Schartmann et al. (2015) as
(12)
Then, in each cell of the domain we computed the respective emissivity values, multiplied them by their npne (assuming full ionisation), and then integrated the data cube along the z coordinate which is perpendicular to the line of sight. It is important to remark that in this calculation it is necessary to take into account the hydrogen nucleus np and electron ne densities appropriately. Since our models consider different amounts of hydrogen in the plasma we include only the material that contains some, i.e. the material coming from the WN89/Ofpe stars. Additionally, we only used the hydrogen fraction of that material which is XH = 11.5% and XH = 40% for the models WR_f07 and WR_f1.
Figure 9 shows the maps over a region of 2″× 2″ centred in Sgr A* that resulted from this procedure. Although the maps were also calculated for the model WR_f07 we chose not to show them here since the emission is negligible (four orders of magnitude smaller) due to the lack of cold gas that contributes to the recombination line flux. From these maps we estimated the flux as seen from Earth, that is, at a distance of 8.33 kpc, so that they could be compared directly with the observed values. These were calculated by integrating within a circular area of a projected radius p = 0.23″ and p = 1″ centred at location of Sgr A*. These radii were motivated by the observations that used an aperture of 0.23″ for extracting the flux observed . However, our models show that the disc actually extends beyond this size, so to get an idea of the whole flux of the model we chose a size of p = 1″ that encloses entirely the disc. The flux values obtained are shown in Table 3 where we also added the upper limit for Brγ from Ciurlo et al. (2021), and the flux of H30α measured from integrating the spectrum across the velocity channels given by Murchikova et al. (2019). In the case of the Brγ line, we can see that only the model WR_f07 has a flux consistent with the upper limit reported for an aperture of p = 0.23″. The model WR_f1 has a flux that is ~100 times higher than the upper limit. Also, we see that the disc in our models extends beyond this projected radius, yet its emission is dominated by the inner p = 0.23″ (see Table 3). The situation is more complex for the H30α line since the observed flux is in both cases higher than the synthetic emission. Although the model that forms a disc is closer it still remains 30% times fainter than the reported value. The model WR_f07 has negligible emission as it is ten order of magnitudes fainter compared to the observation. Thus, it is still not clear how to reconcile the models with both line emission simultaneously as none of them is capable to reproduce the observations. We discuss further the interpretation of these results in Section 7.1.
![]() |
Fig. 8 Density and line-of-sight velocity maps of the central 2.5″× 2.5″ of the model WR_fl at t = 0. The left- and right-hand side panels show projected density weighed by density and line-of-sight velocity weighed by mass, respectively both integrated along the z-axis, which is parallel to the line of sight. |
![]() |
Fig. 9 Line-of-sight integrated maps of the emission of the H30α and Brγ recombination lines shown on the left- and right-hand side panels. The maps correspond to the inner 2″× 2″ of the model WR_f1. |
4.2 Spectral X-ray emission
The X-ray emission from the central parsec is also an observable that has allowed validation of previous numerical models (see Russell et al. 2017; Ressler et al. 2018; Wang et al. 2020). Here we also post-processed our models aiming to obtain the Iron K- alpha spectrum of this region at ~6.7 keV as this is the strongest X-ray feature from Sgr A* (e.g. Russell et al. 2017; Corrales et al. 2020). Our approach is based on the procedure outlined in Balakrishnan et al. (2024b) but adapted to our finite-volume grid-based hydrodynamic model, which is briefly described as follows. First, we assigned energy dependent X-ray emissivities that were taken from the VVAPEC model (Smith et al. 2001), which was obtained with XSPEC (Arnaud 1996). The densities and temperatures were taken from the last output of our simulation, i.e. at t = 0. Additionally, the tracer fields were used to identify the composition fraction of the parcels of gas according to the WR subtypes. Since the optical depth through the simulation domain is optically thin in the X-ray, the radiative transfer equation to solve is simply given by the integration along the line of sight of the emissivities of the respective cells. Given that we can trace the fraction of the material that comes from a given WR subtype this actually corresponds to a linear combination of the tracer field si and the abundance dependent emissivity for each cell, i.e.
(13)
(14)
where ϕk represents a Gaussian function that it is used to obtain the emission across different velocity channels, and to take into account the thermal broadening represented with a velocity of vth. This function is given as follows
(15)
We performed calculations through velocity channels spanning line-of-sight velocities vlos from −3000 km s−1 to 3000 km s−1. These were later converted into energy space via Doppler broadening. Thus, for each energy bin from the emissivity tables we obtained 150 line profiles using a fine resolution of 6400 per dex. In order to sum them appropriately we interpolate them into a unique energy range so that we can obtain a single spectrum that includes both the continuum and line contributions. The result of this procedure is shown in Fig. 10. This contains the X-ray spectrum emitted from a sky-projected annulus of 0.5″ < p < 1.5″ for models WR_f07 (solid blue line) and WR_f1 (dashed orange line). As a reference, we included the estimation from Balakrishnan et al. (2024b) computed from an analogous model that used the SPH approach, which we discuss in the following section. All spectra show qualitative agreement highlighting the different features of the Fe complex at 6.7 keV and 6.9 keV. This is also reflected in the broadening of the lines due to the velocity across the line of sight, which is expected since the models are based on the same stellar orbits and wind properties.
It is relevant to highlight the level of agreement between the finite-volume and SPH models. The qualitative agreement is reasonable and the level of the emission is of the same order of magnitude. The exact quantitative difference among the models could be mainly attributed to the exact cooling table (composition) used in each simulation, which is the result of the different density and temperature distributions that can be observed in Fig. 4. This is why the WR_f07 shows the highest flux level as it has the least efficient cooling and, as a result it has hotter material. Analogously, the run WR_f1 with an enhanced cooling specified in the composition of its plasma results in a reduction of ~40% of its flux. In case of the SPH model, the quantitative differences could be affected due to the intrinsic ability of the generic SPH to capture shocks which impacts directly the density, velocity structure, and the temperature of the medium.
Notice that the flux level of the spectrum of the SPH model by Balakrishnan et al. (2024b) is within the cases explored in this work. Their cooling model is not the same as the ones explored here but this analysis show that it is equivalent to a cooling efficiency between both of our models.
Although not shown here we also compared the X-ray spectra from simulation WR_f1 before and after the disc is formed, and found negligible differences. Finally, we also analysed the slope of the continuum emission within this energy range and found no significant differences among the models. Thus, in this work the presence of the disc does not imprint a feature in the X-ray spectrum that is resolvable by any current or forthcoming X-ray telescope, in agreement with Balakrishnan et al. (2024b).
Radiative flux of the recombination lines from central region.
![]() |
Fig. 10 Synthetic X-ray spectra computed from the numerical models in an annulus of 0.5″ < p < 1.5″. The models WR_f07 and WR_f1 are shown in solid orange and dashed blue lines, respectively. Additionally, the thin solid black line corresponds to the results from the SPH model (Balakrishnan et al. 2024b). |
![]() |
Fig. 11 Density and line-of-sight velocity maps of the central 2.5″× 2.5″ of the SPH model at t = 0. The left- and right-hand side panels show projected density and line-of-sight velocity, respectively. These panels are analogous to the ones shown in Fig. 8. |
5 Comparison with Lagrangian models
Hydrodynamic simulations of the feeding of Sgr A* by the WR stellar winds have been conducted with different numerical tools. In this work, we have presented the results of using a finite- volume approach that solves the hydrodynamic equations in the Eulerian form. However, there has been extensive work on this problem using codes that solve the equations in the Lagrangian form. Specifically, the use of SPH has been the main choice for such a task (Cuadra et al. 2005, 2006, 2008, 2015; Russell et al. 2017; Wang et al. 2020). Despite the limitations of the generic SPH technique (e.g. capturing shocks) the models managed to reproduce the observed accretion rate at the Bondi radius (Cuadra et al. 2008) as well as the X-ray emission (Russell et al. 2017). In this context, we have conducted a comparison of the Eulerian models calculated with Ramses, and the Lagrangian models computed with the SPH code Gadget (Springel 2005). The SPH simulation setup is as similar as it can be to our setup. This corresponds to an updated version from the models in Russell et al. (2017). For more details on the setup of the SPH models, we refer the reader to Balakrishnan et al. (2024a,b). We focus the comparison on the analysis on the final state of the system, specifically in the cold discs formed after 3000 yr.
Figure 11 shows density and line-of-sight velocity maps projected along the line of sight of the SPH model at t = 0 on the left- and right-hand side panels, respectively. Notice that these panels are analogous to the ones shown in Fig. 8 for the Eulerian models presented in this work. First, let us compare the disc density maps between the two approaches. Simply through visual inspection it is possible to see that the structure of the discs is not exactly the same. The Eulerian disc has a more continuous structure towards smaller scales while the Lagrangian disc displays a ring-like structure. This feature can be attributed to the spatial resolution and the inner boundary radius. The Eulerian models indeed manage to resolve smaller scales towards the centre, and the radius of the inner boundary is 5 mas while the Lagrangian models resolve ~ 1 mas and have an inner boundary of 0.1″. As a result, the Lagrangian model obtained a disc 50 times lighter (~10−4 M⊙) than in the Eulerian simulation. This difference could be attributed to the radiative cooling employed in each simulation: WR subtype abundance and Solar with 3 Z⊙ in the Eulerian and the Lagrangian model, respectively. Regarding the projected orientation of the disc both models seem to be in agreement. This can be analysed more carefully on the line-of-sight velocity maps (see right-hand side panels of Figs. 8 and 11). Both maps display that the discs indeed match their projected orientation as well as their most blue- and redshifted velocities are roughly −2000 km s−1 and 2000 km s−1.
Overall, the agreement on the outcome a cold disc around Sgr A* whose properties align between two different approaches indicates that this is a result independent of the numerical technique. Thus, this shows that the determinant factors to obtain this results lies in both the long-enough simulation time (≳3000 yr) as well as efficient radiative cooling.
6 Disc stability analysis
In this section, we present an analysis of the stability of the observed cold disc in the Galactic Centre. Based on its reported properties, the disc should be depleted by accretion after its viscous timescale tν ≈ 43 kyr, assuming an α = 0.1 thin disc (Shakura & Sunyaev 1973), a process much slower than its generation as modelled in our simulations, although faster than its formation out of CND material (Solanki et al. 2023). It is important to bear in mind that our hydrodynamic simulations do not take into account physical viscosity therefore only numerical viscosity is present. Nonetheless, there are several external mechanisms that could be capable of destroying the disc faster than through accretion. Among them are the gravitational and hydrodynamic interactions of the disc with the stars, the stellar-wind-disc collisions due to the presence of the S-star cluster in its location, and/or the evaporation due to heat flowing between the hot medium and the cold disc. Although some of these physical scenarios have been studied analytically and numerically, there has not been any direct application to the cold disc discovered in our Galactic Centre.
6.1 Stellar cluster perturbing a thin disc
The stability of a cold disc perturbed by stellar passages has been investigated by Ostriker (1983). Their work presented an analytical approach to calculate how the disc loses angular momentum due to the passages, considering the integrated effect of a stellar cluster in stationary state. Besides this work, most efforts have been put on following how the stellar dynamics are affected by the disc presence (e.g. Fabj et al. 2020, and references therein) or on the effect of embedded accreting stars and/or black holes on the disc (e.g. Tagawa et al. 2021) rather than on the consequences to the disc properties.
In this section, we summarise the approach of Ostriker (1983) and apply it to the Galactic centre case. Let us start considering a cold, thin accretion disc around a compact object. Within the thin disc we assume that its height H(r) is small compared to the radial coordinate r, i.e. H(r)/r << 1. Also, we consider that the sound speed of the material in the disc is small compared to the speed of the stars at the same radius. In this scenario, a star passing through the disc will produce a torque that removes angular momentum from the disc. We consider that the star interacts with the disc gravitationally and hydrodynamically The gravitational interaction of a star crossing the disc can be described using the Bondi-Hoyle approach (Bondi & Hoyle 1944), while the hydrodynamic interaction is just a physical collision of a solid sphere with a gaseous disc.
Let us consider a system of reference where the x axis is in the direction of rotation and z is in the direction of the rotation pole of the disc. Then, according to Ostriker (1983) the total change of linear momentum along the x direction Δpx for a single stellar passage is
(16)
where R* is the radius of the star, Σd = ʃdiscρdH(r) is the surface density of the disc, ρd is the volume density of the disc, q is the hydrodynamic drag parameter that is set to q = 2 assuming a large Mach number collision, v0 is the local stellar velocity, ΛD ≈ is the Coulomb logarithm, v* is the escape velocity from the surface of the star, η = v*/v0 is the hardness parameter, vd is the disc velocity, vrel is the relative velocity vector, and (θ, ϕ) are the spherical coordinate angles along the polar and azimuthal directions.
In order to consider the effect of the complete stellar cluster interacting with the disc we need to integrate over the cluster phase-space volume. First, we integrate over the local stellar velocity distribution. Let us define dN as the number of stars in the velocity range v → v + dv, i.e
(17)
where f(v) is the velocity distribution function, n* is the stellar number density. Then, the number of crossings per unit time and area in the (θ, ϕ) direction is
(18)
Now, integrating over the phase space, the total transfer of momentum to the disc per unit area at r is
(19)
where µ = cos ϕ.
As the local momentum per unit area in the disc is px,d(r) = Σ(r)/τd(r) the inverse of the drag timescale to remove the angular momentum of the disc is τd = −px,d/ṗx,tot can be written as
(20)
where the quantities (I0, I1) depend solely on the stellar distribution function. For instance, in the case of a Maxwellian distribution the values of these parameters can be obtained numerically and are I0 = 12.553 and I1 = 0.474.
In order to apply this model to the Galactic Centre we need to calculate the stellar density within the size of the disc. Following Schödel et al. (2007) the stellar mass density at r < 6 arcsec is described by
(21)
Assuming that the stellar mass density is composed mainly by stars whose radii are R* = 1 R⊙ and move on their orbits typically at v0 = 1000 km s−1 we can use Eq. (20) to estimate the inverse of the drag timescale,
(22)
Based on this calculation the angular momentum of the cold disc should be removed after ~185 Myr of interactions of the stars. This is much longer than the formation time-scale found in our simulations. Furthermore, it is even much longer than the age of the Wolf-Rayet stars, which is about 6 Myr (Martins et al. 2007). Then, the gravitational and hydrodynamic interactions of the stars onto the disc should not affect the stability of the disc.
6.2 Stellar winds perturbing a disc
In the case of stellar irradiation affecting the state of the accretion flow there have been many studies mainly motivated by the recent pericentre passages of the star S2 around Sgr A*. For instance, Cuadra et al. (2003) could rule out the existence of an optically thick disc based on the lack of its thermally reprocessed emission as it gets illuminated by S2. Nayakshin (2004) estimated that a star as luminous as S2 could heat up and ionise an inner disc, which could enhance the accretion rate. Giannios & Sironi (2013); Christie et al. (2016) calculated the bremsstrahlung emission as a result of the stellar wind shocking a Radiatively Inefficient Accretion Flow (RIAF) around Sgr A*, finding an X-ray luminosity comparable to quiescent emission of Sgr A*. In fact, no noticeable increase was detected during 2018 (see Table 1 of Andrés et al. 2022). Hosseini et al. (2020) constrained the observed L′-band variability of S2 to be about 2–3%, based on a model of bow shock of its stellar wind, implying an ambient density <105 cm−3, which rules out the presence of a standard accretion disc at S2’s pericentre, in agreement with the previous studies. It is important to remark however that such a limit marginally allows the existence of both the disc reported by Murchikova et al. (2019), and also the one in our simulation.
The stars in the S cluster also have relatively powerful stellar winds. As the stars inhabit the black hole in the vicinity of the cold disc, it is expected that the winds interact with it. By either opening a bubble due to the high kinetic energy carried in the winds or depositing thermal energy they affect the disc properties. We proceed to perform analytical estimates to quantify the impact of these processes.
6.2.1 Bubbles in the disc
Let us consider a star blowing an isotropic stellar wind with a mass-loss rate Ṁw at a terminal speed vw. Then, its density ρw at a distance r from the star is given by
(23)
If such a star is immersed in a medium whose number density is nm at rest at a temperature Tm, the radius of the bubble rb opened by the wind can be calculated equating the ram pressure of the wind and the thermal energy of the medium, i.e.
(24)
where we have used as fiducial values the disc parameters measured by Murchikova et al. (2019) and typical stellar wind properties for B stars (Oskinova et al. 2011; Krtička 2014). Then, the radius of the bubble is about 0.012 arcsec, which is approximately one tenth of the disc radius.
In order to open such a bubble a star would need to be within the disc for at least the wind crossing time, which is the minimum time it takes for a bubble of that size, i.e. rb/vw ≈ 0.5 yr. Let us estimate the duration of a stellar passage through the disc. Along the perpendicular direction a star crossing the disc would take at most 0.1 Rd/v0, being H(r)/r = 0.1. One S star moves typically at v0 ≈ 1000 km s−1, then 0.1 Rd/v0 ≈ 1 yr. As both timescales are comparable, we conclude that the bubble can be opened.
For such a perturbation to last we need to check if either shear or sound waves could eliminate it, i.e. close the bubble. In the case of shear, this is given by the orbital speed of the disc and the size of the bubble:
(25)
where Δvkep corresponds to the difference in Keplerian velocity at the centre and edge of the bubble. As this expression increases with r, the largest value will be at the outer edge of the disc, then tclose < 5 yr. This value should be interpreted as the longest duration of a perturbation of the wind onto the disc. As there are 12 S-stars on orbits whose pericentre distances are shorter than the radius of the disc and their orbital periods are of the order of ~10 yr (Gillessen et al. 2017) we expect a bubble to be opened every year on average. However, as shear would close such bubbles in less than five years there would be at most four bubbles opened on average in a stationary state. As the size of the bubbles is very small compared to the size of the disc (two orders of magnitude in area), four of them would not have an impact on its structure. The sound crossing timescale rb/cs ≈ 30 yr is longer than the shear timescale so we do not expect it to have an impact on erasing perturbations.
6.2.2 Winds injecting thermal energy
The S stars are of spectral type B with masses of 10 M⊙ (e.g. Habibi et al. 2017), and thus their mass-loss rate must be ~10−8 M⊙ yr−1 with terminal velocity of 1000 km s−1. The supersonic nature of the winds should be enough to compress the shocked material at high temperature. Potentially, this thermal energy could be deposited onto the disc during each stellar passage. Ifwe consider that most of the kinetic energy of the wind is transformed into thermal energy about ~6.4 × 1033 erg would be deposited in the disc during the ~1 yr long passage (see above). In order to check if the wind energy can significantly impact the thermodynamic state of the disc let us estimate its total internal energy Uint.
(26)
where kB is the Boltzmann constant and Vd is the disc volume. Replacing the observed disc properties we obtain that Uint = 1.63 × 1042 erg. From this calculation it is clear that even if all the energy from the wind in a stellar passage is deposited the contribution is still negligible to modify the state of the disc. Even if we consider the integrated effect of all the S stars passages on a timescale comparable to the viscous timescale of the disc (~43 kyr) the contribution is still <0.1% of the total thermal energy of the disc. Thus, the energy supplied by the stellar winds into the disc is not enough to affect its state.
6.3 Thermal conduction
The vicinity of Sgr A* at the Bondi radius is made out of a diffuse (10 cm−3), hot medium (107 K; Baganoff et al. 2003) due to the shocked stellar winds from the Wolf-Rayet stars in the region. At the disc vicinity (~103 RSch), the conditions are more extreme as the plasma has a density of about ~5 × 103 cm−3 (Gillessen et al. 2019) and temperature of ~108 K. As the disc is significantly colder it is expected that heat flows from the medium to the disc. The large temperature gradient should enable thermal conduction to take place and likely evaporate the disc.
The problem of evaporating a cold structure sitting in a hot environment due to thermal conduction was first studied by Cowie & McKee (1977). In that work, the authors derived an analytical expression for the mass-loss rate and evaporation timescale for a spherically symmetric gas cloud. Unfortunately, these expressions are not valid for more complex geometries. Meyer & Meyer-Hofmeister (1994) studied specifically the evaporation of a cold thin accretion disc in a hot corona. Although the model considered only one radial zone it was able to calculate the full perpendicular structure of the accretion disc and the corona by means of numerical simulations.
Liu et al. (2004) used the model by Meyer & Meyer-Hofmeister (1994) to analyse which solutions were consistent with the observational data available at the time in case there were a cold disc in the center of the Galaxy. They found that any transient disc would evaporate quickly. With their obtained evaporation rate of ~10−4M⊙ yr−1, the life-time of a disc with properties such as found by Murchikova et al. (2019) or formed in our models would be ~1 yr. Similar analysis have argued that small clumps (~M⊕) in the Galactic Centre should also evaporate quickly (~1 yr) due to thermal conduction (Burkert et al. 2012; Calderón et al. 2018). According to this, cold gas should not be able to last long if formed, yet we observe many cold gas structures in the central parsec: the many dusty cold clumps detected in the vicinity of IRS 13E (Fritz et al. 2010; Peißker et al. 2023), and a larger gaseous cloud X7 (Ciurlo et al. 2023). Thus, this regime of classical thermal conduction might not be applicable to the Galactic centre environment. Another possibility has been studied by Nayakshin (2004) who focused on the same problem and found solutions in a different regime, where the electron mean free path is long enough that an accretion disc does not evaporate but rather increases its mass by condensation. A more thorough analysis of the role of thermal conduction is beyond the scope of this work and is deferred to future study.
Our findings shows that gravitational and hydrodynamic effects require timescales that are too long to have an impact on the disc stability. On the other hand, we find that thermal conduction in the classical case should evaporate the observed disc in about one year. Given the time-scale for disc formation found in our model, the Nayakshin (2004) condensation regime remains to our knowledge as the only viable physical scenario that would allow its continued existence and potential identification with Murchikova et al. (2019)’s reported disc.
7 Discussion
In this section, we interpret our findings exploring the uncertainties in the abundances of WR stars in general and of the ones in the central parsec. Additionally, we contrast the properties of the simulated and observed cold discs as well as discuss further implications.
7.1 Is the simulated disc the observed disc?
In order to compare the simulated and the observed discs, first we analyse their physical properties. The simulated disc has a total mass of ~5 × 10−3 M⊙, while the structure observed in H30α was reported to have a mass of 10−5 − 10−4 M⊙. However, comparing these quantities might not be appropriate since the mass of the observed structure is based on three assumptions: a thin disc with uniform density and temperature, and a masing factor of ~ 100. All of them minimise the value of the mass inferred from the observations. For instance, dropping the maser assumption the inferred mass increases in a factor ten (see Eqs. (30) and (31) of Supplementary Material in Murchikova et al. 2019). Bearing this in mind, we could argue that the mass of the simulated and observed disc might be of the same order provided there was no maser. Although removing that assumption would create tension reconciling both the Brγ and H30α observed fluxes.
Regarding its extension, the simulated disc has a projected radius of ~0.5″, while the observed disc has a radius of 0.23″. However, since most of the recombination line emission of the disc is originated from its central region (p < 0.23″, see Table 3) the sizes of the simulated and observed discs are in agreement.
A more direct comparison can be done through the analysis of the mock observational quantities computed from the simulations. In Section 4.1, we computed the H30α and Brγ recombination line emission and estimated their fluxes. Here, we have found that both synthetic fluxes are in tension with the observations. First, the simulated Brγ flux is 100 times higher than the upper limit. It is relevant to remark that this limit is already extinction corrected. Second, the simulated H30α flux is 30% lower than the observed one. These results differ from the estimates by Ciurlo et al. (2021) from our previous models due to, again, the uniform disc assumptions. The simulated disc is not thin or possess uniform density or temperature. Although the mean density and temperature of the simulated disc is of the order of ~10−5 cm−3 and ~104 K its geometry is more complex with a larger variance in its properties. As a result, both the synthetic Brγ and H30α fluxes are much higher than the ones inferred by Ciurlo et al. (2021).
Based on the line-of-sight velocity maps shown in Murchikova et al. (2019), the orientation of the observed disc displays the red- and blueshifted sides to be on the East (left-hand) and West (right-hand) sides, respectively. Although our model WR_f1 also shows this configuration the exact inclination of the discs does not match perfectly. Specifically, the discs are inclined in ~90° among each other, being the simulated disc tilted in the clockwise direction on the sky. In principle, this argues against the WR stars as the main responsible for the formation of the observed disc. But we cannot rule out completely their role in the process, since a combination of the WR stars with other gaseous structures, such as the minispiral or the CND, acting simultaneously might result in a net effect that manages to reproduce the observed angular momentum. Solanki et al. (2023) showed that material from the circumnuclear disc could fall close to Sgr A* though on much longer timescales. Certainly, this scenario would be much more difficult and challenging to simulate with the appropriate chemical abundances as well as with high-enough spatial resolution.
It is also relevant to notice that the observed accretion flow orientation at ≲10 RSch, where RSch corresponds to the Schwarzschild radius of Sgr A*, is consistent with the orientation of the clockwise stellar disc (GRAVITY Collaboration 2020; Event Horizon Telescope Collaboration 2022). Thus, our models reproduce that connection between the dynamics from large (>105 RSch) to small (≲10 RSch) scales. However, the observed cold disc does not seem to follow such a connection.
Finally, an intriguing aspect of our findings is that the WR atmospheres with more H in them would favour the formation of the disc. At the same time, the disc indeed was detected in the H30α recombination line, which obviously indicates the presence of H. Thus, it is sensible to ask how much H can be in the winds of WR stars but, specifically in the ones of the subtype WN89/Ofpe. This point is discussed in the following section.
7.2 Uncertainties in abundances
The results of the hydrodynamic models demonstrate that the chemical abundances of the material in the WR winds are the key factor to determining whether or not a cold disc can be formed due to their action. Unfortunately, there are two large uncertainties in this context. First, the metallicity of the young population of the nuclear star cluster is not well constrained (e.g. Genzel et al. 2010). Although there is agreement that it is above Solar its exact value has not been established so far, with most studies quoting the range Z = 2–3 Z⊙. However, even if the metallicity is determined more precisely, the problem would not be entirely solved, since the abundances of the plasma in the region are dominated by the material supplied by the WR star winds whose compositions are much more uncertain. Theoretically, the evolution of these stars is a topic of active research, and improved prescriptions for mass loss are changing the modelled properties of observed WR stars (see e.g. Gormaz-Matamala et al. 2023, and references therein). Simply based on observations, up to now, it has been reported that some of these stars might possess from XH = 0 to XH = 50% for Galactic WR stars of the WN subtype (Hamann et al. 2006, 2019), and even as high as XH = 60% for extragalactic ones (Todt et al. 2015). In the vicinity of Sgr A*, Martins et al. (2007) characterised the wind properties of the sample of WR stars through the analysis of their infrared spectra. They reported abundance H/He ratios in the range 2–5 for these stars which does not allow us to rule out the scenario regarding this aspect. A variation on the hydrogen abundance in the chemical mixture in a factor five might change the net radiative cooling rate by a factor 1.1–1.4, especially at lower temperatures. In our models, WR_f07 considered a fiducial value of hydrogen mass fraction of 11.5% (Russell et al. 2017) but the model WR_f1 assumed 40%. Such a change was enough to modify the cooling function by about 30%, which ended up affecting significantly the final result and the state of the plasma. Thus, within our current understanding of the abundances in the WR stellar atmospheres both scenario are plausible. Only better quality spectra and more sophisticated models of them could allow us to constrain the wind properties and their composition more precisely.
8 Conclusions
We present a numerical study of the system of mass-losing stars feeding Sgr A*, focusing on the impact of the chemical composition of the plasma. Through studying simulations with different abundances, we find that the system can either form or not form a cold disc around the black hole, provided that the model is run for sufficient timescales (≳3000 yr) to create a dense enough medium in the inner parsec of the Galaxy. As in our previous work, we confirm that, if formed, the cold disc can significantly impact the hydrodynamic and thermodynamic state of the plasma within the central parsec. As a result, the inferred mass-inflow rate at 5 × 10−4 pc can reach up to 8 × 10−6 M⊙ yr−1, which is between four and eight times higher than the value without the presence of such a disc. Our models still point to the mass-losing stars in the clockwise disc being the main source of the accreted material due to their dense and slow winds, as previous models have shown (e.g. Cuadra et al. 2008; Ressler et al. 2018; Calderón et al. 2020b). In this work, we have followed the origin of the inflowing material more carefully through the use of passive scalar fields for different WR subtypes. This allowed us to confirm the origin of the accreted material but more importantly led us to conclude that this is mostly provided by the winds of the stars of the WR subtype WN89/Ofpe. This provides evidence that the infalling material must have a non-negligible fraction of hydrogen, as these WR stars have not entirely lost their atmospheric hydrogen via winds.
Our models also show that the formation of a cold disc depends strictly on the radiative cooling employed, which is determined by the chemical abundances of the plasma. Previous models only considered solar abundances with 3 Z⊙ but do not agree on whether or not a cold disc can be formed. From our results in Appendix B, we speculate that at Z ≈ 3 Z⊙ there is a transition between the regimes of persistent quasi-steady accretion and disc formation, which would imply that small differences in the cooling function or numerical implementation can result in models reproducing either regime. This may explain the differing results obtained by Calderón et al. (2020b) and Ressler et al. (2020) with similar numerical and physical setups.
In the case where a cold disc is formed around Sgr A*, the structure may resemble the extension and the line-of-sight velocity structure of the observed disc. However, we are not able to reproduce direct observational quantities. Specifically, the exact sky-projected inclination of the simulated and observed discs differs by ~90°. The upper limit on the Brγ emission line is 100 lower than the synthetic flux. The observed H30α emission line flux is 30% higher than in our model.
Furthermore, we contrasted the Eulerian (finite-volume grid-based) models of this work with analogous (SPH) Lagrangian models in order to address the potential influence of the chosen approach on the outcome of the simulations. Balakrishnan et al. (2024b) showed that a cold disc is also formed around Sgr A* if the simulation is run for long timescales. Both cold structures agree on their sizes, their velocity structure, and their sky-projected inclination. Although there are differences in the total mass and specific inner structure of the discs, these could be attributed to the resolution employed and the exact assumption to consider the innermost boundary. Despite these differences, the agreement is reasonable among the numerical approaches, especially when analysing the synthetic X-ray emission (see Fig. 10). The spectral features agree qualitatively across the models presented in this work and the Lagrangian model. The levels of the X-ray emission are all within the same order of magnitude but the exact base level depends on the radiative cooling chosen, which is determined by the chemical abundances in the WR atmospheres.
We have also explored the stability and long-term evolution of the putative disc under the effect of perturbing agents: the nuclear star cluster, the wind of all the mass-losing stars (including the S stars) impacting the disc through opening bubbles and/or depositing thermal energy via shocks, and the potential effect of thermal conduction. Our analysis shows that gravitational and hydrodynamic effects are likely to impact the disc on longer timescales than the viscous timescale. However, within the classical thermal conduction framework, a cold disc in such a region should not exist. Based on this and the disc-formation timescale from our simulations, we speculate that the condensation model of Nayakshin (2004) remains a plausible mechanism consistent with the existence of a disc.
In conclusion, through the present work, we have been able to reconcile the discrepancy among the numerical simulations of the system of mass-losing stars feeding Sgr A*. The formation of a cold disc on a ~3000 yr timescale is indeed possible for certain chemical abundances that are consistent with the current observational constraints. Nonetheless, it is not possible to favour or disfavour this scenario due to the related uncertainties. Thus, high-resolution spectra of the WR stars and more sophisticated modelling of them could allow us to discern the validity of such scenarios. We must warn that even in the case where the cold disc is a result of the action of the WR stars, it is not possible to reproduce all of its observed properties, specifically its inclination and recombination line fluxes. More sophisticated models that include infalling material from other sources might be necessary to produce the exact inclination of the structure, which could also impact the emission of the structure.
Acknowledgements
We thank Prof. Sergei Nayakshin and Prof. Stanley Owocki for very helpful discussions about this project. We also thank Dr. Sean M. Ressler for sharing the radiative cooling function shown in Figure 1. Additionally, we are grateful to Dr. Álex Gormaz-Matamala for helping with references and discussions on the WR wind abundances. DC and JC acknowledge the support of the Kavli Foundation through its summer program at the Max Planck Institute for Astrophysics where fruitful discussions took place. In addition, DC thanks the warm hospitality of the Universidad Adolfo Ibañez in Chile and University of Delaware in the USA where part of this project was developed. We acknowledge the support from ANID in Chile through FONDECYT Regular grant 1211429 and Millennium Science Initiative Program NCN2023_002. The research of DC and SR was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2121 – “Quantum Universe” – 390833306. Since 01.10.2024, DC has been funded by the Alexander von Humboldt Foundation immediately. C.M.P.R. acknowledges support from NASA Chandra Theory grant number TM3-24001X. Resources supporting this work’s SPH simulations were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercom-puting (NAS) Division at Ames Research Center. SR is also supported by the European Research Council (ERC) Advanced Grant INSPIRATION under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 101053985), the Swedish Research Council (VR) under grant number 2020-05044, and the research environment grant “Gravitational Radiation and Electromagnetic Astrophysical Transients” (GREAT) funded by the Swedish Research Council (VR) under Dnr 2016-06012, by the Knut and Alice Wallenberg Foundation under grant Dnr. KAW 2019.0112. The numerical simulations of this work were run on the high-performance computing system COBRA of the Max Planck Computing and Data Facility. Most of the analysis of the numerical data was carried out using the PYTHON package YT (Turk et al. 2011). Furthermore, this work made use of PYTHON libraries NUMPY (Harris et al. 2020) and MATPLOTLIB (Hunter 2007), as well as of the NASA’s Astrophysics Data System.
Appendix A Numerical setup choices
Calderón et al. (2020b) simulated the system of WR stars feeding Sgr A* while orbiting around it following an analogous numerical setup compared to the models presented in this work. There, we had found that the resulting cold disc tended to align with the grid after its formation. In order to remedy this numerical artifact we investigated the use of a MinMod slope limiter, instead of MonCen. However, we warned that it was not straightforward to maintain the spherical symmetry of the winds using the Min-Mod slope limiter. In this work, we investigated more carefully these numerical choices in order to be confident that our results are robust.
First, we investigated the choice of the MonCen slope limiter. To do this, we made use of the same setup presented in this work but paying attention to the values of the passive scalar fields. In principle, the values of these fields should be between zero and one. However, this experiment showed that the MonCen slope limiter produced non-physical oscillations that in some cases resulted in negative values of the passive scalar fields. This behaviour was not observed when using the MinMod slope limiter. Thus, we selected this choice throughout this work. To ensure the spherical symmetry of the stellar winds we had to impose a less restrictive refinement criterion based on density gradients as well as a smoother transition between refinement levels through the parameter NEXPAND.
Second, we tested the choice of the Riemann solver employed but using the MinMod slope limiter. Currently, the HLLC Rie-mann solver (Toro et al. 1994) is widely used in grid-based hydrodynamic simulations. In this work, we have found that this choice results in a cold disc that tends to align with the Cartesian grid as seen in Calderón et al. (2020b). This could indicate problems in conserving the angular momentum which is a known issue with simulations in Cartesian grids when the relevant regions are not well resolved. With the usage of the exact Riemann solver (Toro 2009), though more computationally expensive, it was possible to avoid this artificial alignment. Thus, the only sensible combination was the use of the exact Riemann solver with a MinMod slope limiter but we had to ensure to count with enough resolution and smoothness in the transition of the refinement levels. As a reference of these experiments, Figure A.1 shows density projection (weighted by density) maps of the final state of the system for two runs with Solar composition with 3Z⊙ but with different Riemann solvers. The left- and right-hand side panels display the models using the HLLC and the exact Riemann solvers, respectively. Although overall the density structure looks similar there are subtle differences. The most relevant is related to the disc orientation and structure. On the left-hand side, the disc tries to align with one of the Cartesian axes, while on the right-hand side the disc remains consistently with the shown orientation.
![]() |
Fig. A.1 Comparison of the final state (present time) of the simulations with different Riemann solvers. Left-and right-hand side panels correspond to runs with Harten-Lax-van Leer-Contact (HLLC) and exact Riemann solvers, respectively. Both models were run using the MinMod slope limiter. |
Appendix B Solar composition varying metallicity
Before investigating the models with the new implementation of the differential cooling we ran some experiments simply modifying the contribution of the metals to the total radiative cooling. We simulated three models varying only the composition used: A1 used Solar composition, A3 used Solar composition with 3Z⊙, and A5 used Solar composition with 5Z⊙. The density projected (weighted by density) maps of the final state of the simulations are shown in Figure B.1. The left-hand side, central, and right-hand side panels display the runs A1, A3, and A5, respectively. A higher metallicity in the plasma increases the radiative cooling directly. As a result, more and denser filaments, clumps, and a denser central disc are obtained with increasing the metal-licity. On the contrary, if the metallicity is decreased to a Solar value we observe almost no overdensities, and definitively no disc formation. We speculate that at Z ≈ 3Z⊙ there is a transition between the regimes of persistent quasi-steady accretion and disc formation, which would imply that small differences in the cooling function or numerical implementation can result in models reproducing either regime. This may explain the differing results obtained by Calderón et al. (2020b) and Ressler et al. (2020) with similar numerical and physical setups.
In order to perform a more quantitative comparison of these results Figure B.2 presents radial profiles of density (weighted by volume) and temperature (weighted by mass) on the top and bottom panels, respectively. The dashed blue, solid orange, and blue dotted-dashed lines represent the models A1, A3, and A5, respectively. The density profiles clearly shows the impact of cooling in the structure of the medium. More efficient cooling, like in cases A3 and A5, results in an order of magnitude denser medium within 1″ compared to the case with inefficient cooling, i.e. A1. Notice that even higher metallicity in A5 extends the impact the denser region to slightly larger spatial scales. The effect is similar in the temperature profile. In A1, the temperature profile decays with r within the central 1″. But when the metallicity is high enough to make cooling efficient the temperature profile can decrease up to two orders of magnitude depending on the exact value.
In summary, these experiments allowed us to quantify the impact of radiative cooling by simply decreasing or increasing the radiative cooling influence. However, since the chemical composition in reality is much more complicated and differs from Solar composition significantly these models should be interpreted only as general guides to assess the role of radiative cooling. In the main text of this work, we devoted to explore more physically motivated mixtures that, although depend on more uncertain parameters are more robust when attempting to build a physical model that can be constrasted with observational data.
![]() |
Fig. B.1 Comparison of the simulations with Solar composition varying the metallicity at two different simulation times. The panels show projected density maps weighed by density along the ɀ axis, i.e. ∫ ρ2dɀ/ ∫ ρdɀ, which is parallel to the line of sight. Left-hand side, central, and right-hand side panels show the runs A1, A3, and A5, respectively. All maps display the full computational domain. |
![]() |
Fig. B.2 Radial profiles of time-averaged volume-weighted density (top) and mass-weighted temperature (bottom) over the last 500 yr of simulation time. The models A1, A3, and A5 are shown in solid orange, dashed blue, and dashed-dotted green lines, respectively. |
References
- Andrés, A., van den Eijnden, J., Degenaar, N., et al. 2022, MNRAS, 510, 2851 [CrossRef] [Google Scholar]
- Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, 101, Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, 17 [NASA ADS] [Google Scholar]
- Baganoff, F. K., Maeda, Y., Morris, M., et al. 2003, ApJ, 591, 891 [Google Scholar]
- Balakrishnan, M., Corrales, L., Markoff, S., et al. 2024a, ApJ, 974, 98 [NASA ADS] [CrossRef] [Google Scholar]
- Balakrishnan, M., Russell, C. M. P., Corrales, L., et al. 2024b, ApJ, 974, 99 [NASA ADS] [CrossRef] [Google Scholar]
- Becklin, E. E., Gatley, I., & Werner, M. W. 1982, ApJ, 258, 135 [Google Scholar]
- Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [Google Scholar]
- Burkert, A., Schartmann, M., Alig, C., et al. 2012, ApJ, 750, 58 [CrossRef] [Google Scholar]
- Calderón, D., Ballone, A., Cuadra, J., et al. 2016, MNRAS, 455, 4388 [Google Scholar]
- Calderón, D., Cuadra, J., Schartmann, M., et al. 2018, MNRAS, 478, 3494 [CrossRef] [Google Scholar]
- Calderón, D., Cuadra, J., Schartmann, M., et al. 2020a, MNRAS, 493, 447 [Google Scholar]
- Calderón, D., Cuadra, J., Schartmann, M., Burkert, A., & Russell, C. M. P. 2020b, ApJ, 888, L2 [Google Scholar]
- Christie, I. M., Petropoulou, M., Mimica, P., & Giannios, D. 2016, MNRAS, 459, 2420 [NASA ADS] [CrossRef] [Google Scholar]
- Ciurlo, A., Morris, M. R., Campbell, R. D., et al. 2021, ApJ, 910, 143 [NASA ADS] [CrossRef] [Google Scholar]
- Ciurlo, A., Campbell, R. D., Morris, M. R., et al. 2023, ApJ, 944, 136 [NASA ADS] [CrossRef] [Google Scholar]
- Corrales, L., Baganoff, F. K., Wang, Q. D., et al. 2020, ApJ, 891, 71 [NASA ADS] [CrossRef] [Google Scholar]
- Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135 [NASA ADS] [CrossRef] [Google Scholar]
- Crowther, P. A. 2007, ARA&A, 45, 177 [Google Scholar]
- Cuadra, J., Nayakshin, S., & Sunyaev, R. 2003, A&A, 411, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Cuadra, J., Nayakshin, S., Springel, V., & Di Matteo, T. 2005, MNRAS, 360, L55 [CrossRef] [Google Scholar]
- Cuadra, J., Nayakshin, S., Springel, V., & Di Matteo, T. 2006, MNRAS, 366, 358 [NASA ADS] [CrossRef] [Google Scholar]
- Cuadra, J., Nayakshin, S., & Martins, F. 2008, MNRAS, 383, 458 [NASA ADS] [Google Scholar]
- Cuadra, J., Nayakshin, S., & Wang, Q. D. 2015, MNRAS, 450, 277 [CrossRef] [Google Scholar]
- Dinh, C. K., Salas, J. M., Morris, M. R., & Naoz, S. 2021, ApJ, 920, 79 [NASA ADS] [CrossRef] [Google Scholar]
- Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2022, ApJ, 930, L12 [NASA ADS] [CrossRef] [Google Scholar]
- Fabj, G., Nasim, S. S., Caban, F., et al. 2020, MNRAS, 499, 2608 [NASA ADS] [CrossRef] [Google Scholar]
- Fritz, T. K., Gillessen, S., Dodds-Eden, K., et al. 2010, ApJ, 721, 395 [Google Scholar]
- Genzel, R. 1989, in The Center of the Galaxy, 136, ed. M. Morris, 393 [Google Scholar]
- Genzel, R., Eisenhauer, F., & Gillessen, S. 2010, Rev. Mod. Phys., 82, 3121 [Google Scholar]
- Giannios, D., & Sironi, L. 2013, MNRAS, 433, L25 [NASA ADS] [CrossRef] [Google Scholar]
- Gillessen, S., Plewa, P. M., Eisenhauer, F., et al. 2017, ApJ, 837, 30 [Google Scholar]
- Gillessen, S., Plewa, P. M., Widmann, F., et al. 2019, ApJ, 871, 126 [NASA ADS] [CrossRef] [Google Scholar]
- Gormaz-Matamala, A. C., Cuadra, J., Meynet, G., & Curé, M. 2023, A&A, 673, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- GRAVITY Collaboration (Abuter, R., et al.) 2019, A&A, 625, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- GRAVITY Collaboration (Bauböck, M., et al.) 2020, A&A, 635, A143 [CrossRef] [EDP Sciences] [Google Scholar]
- GRAVITY Collaboration (Abuter, R., et al.) 2021, A&A, 647, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Habibi, M., Gillessen, S., Martins, F., et al. 2017, ApJ, 847, 120 [Google Scholar]
- Hamann, W. R., Gräfener, G., & Liermann, A. 2006, A&A, 457, 1015 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Hamann, W. R., Gräfener, G., Liermann, A., et al. 2019, A&A, 625, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
- Herald, J. E., Hillier, D. J., & Schulte-Ladbeck, R. E. 2001, ApJ, 548, 932 [Google Scholar]
- Hosseini, S. E., Zajac?ek, M., Eckart, A., Sabha, N. B., & Labadie, L. 2020, A&A, 644, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
- James, T. A., Viti, S., Yusef-Zadeh, F., Royster, M., & Wardle, M. 2021, ApJ, 916, 69 [NASA ADS] [CrossRef] [Google Scholar]
- Krticka, J. 2014, A&A, 564, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Lemaster, M. N., Stone, J. M., & Gardiner, T. A. 2007, ApJ, 662, 582 [CrossRef] [Google Scholar]
- Liu, B. F., Meyer, F., & Meyer-Hofmeister, E. 2004, A&A, 421, 659 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Lodders, K. 2003, ApJ, 591, 1220 [Google Scholar]
- Martins, F., Genzel, R., Hillier, D. J., et al. 2007, A&A, 468, 233 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Meyer, F., & Meyer-Hofmeister, E. 1994, A&A, 288, 175 [NASA ADS] [Google Scholar]
- Murchikova, E. M., Phinney, E. S., Pancoast, A., & Blandford, R. D. 2019, Nature, 570, 83 [NASA ADS] [CrossRef] [Google Scholar]
- Nayakshin, S. 2004, [arXiv preprint astro-ph/0402469] [Google Scholar]
- Nayakshin, S., Cuadra, J., & Sunyaev, R. 2004, A&A, 413, 173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Onifer, A., Heger, A., & Abdallah, J. 2008, in Astronomical Society of the Pacific Conference Series, 391, Hydrogen-Deficient Stars, eds. A. Werner, & T. Rauch, 305 [NASA ADS] [Google Scholar]
- Oskinova, L. M., Todt, H., Ignace, R., et al. 2011, MNRAS, 416, 1456 [NASA ADS] [CrossRef] [Google Scholar]
- Ostriker, J. P. 1983, ApJ, 273, 99 [NASA ADS] [CrossRef] [Google Scholar]
- Paumard, T., Genzel, R., Martins, F., et al. 2006, ApJ, 643, 1011 [NASA ADS] [CrossRef] [Google Scholar]
- Peißker, F., Zajac?ek, M., Thomkins, L., et al. 2023, ApJ, 956, 70 [CrossRef] [Google Scholar]
- Quataert, E. 2004, ApJ, 613, 322 [NASA ADS] [CrossRef] [Google Scholar]
- Ressler, S. M., Quataert, E., & Stone, J. M. 2018, MNRAS, 478, 3544 [NASA ADS] [CrossRef] [Google Scholar]
- Ressler, S. M., Quataert, E., & Stone, J. M. 2020, MNRAS, 492, 3272 [NASA ADS] [CrossRef] [Google Scholar]
- Russell, C. M. P., Wang, Q. D., & Cuadra, J. 2017, MNRAS, 464, 4958 [CrossRef] [Google Scholar]
- Schartmann, M., Ballone, A., Burkert, A., et al. 2015, ApJ, 811, 155 [Google Scholar]
- Schödel, R., Eckart, A., Alexander, T., et al. 2007, A&A, 469, 125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
- Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91 [Google Scholar]
- Solanki, S., Ressler, S. M., Murchikova, L., Stone, J. M., & Morris, M. R. 2023, ApJ, 953, 22 [NASA ADS] [CrossRef] [Google Scholar]
- Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
- Štofanová, L., Kaastra, J., Mehdipour, M., & de Plaa, J. 2021, A&A, 655, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Stone, J. M., Tomida, K., White, C. J., & Felker, K. G. 2020, ApJS, 249, 4 [NASA ADS] [CrossRef] [Google Scholar]
- Tagawa, H., Kimura, S. S., Haiman, Z., et al. 2021, arXiv e-prints [arXiv:2112.01544] [Google Scholar]
- Teyssier, R. 2002, A&A, 385, 337 [CrossRef] [EDP Sciences] [Google Scholar]
- Todt, H., Sander, A., Hainich, R., et al. 2015, A&A, 579, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Toro, E. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction (Berlin: Springer) [CrossRef] [Google Scholar]
- Toro, E. F., Spruce, M., & Speares, W. 1994, Shock Waves, 4, 25 [Google Scholar]
- Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9 [NASA ADS] [CrossRef] [Google Scholar]
- von Fellenberg, S. D., Gillessen, S., Stadler, J., et al. 2022, ApJ, 932, L6 [NASA ADS] [CrossRef] [Google Scholar]
- Wang, Q. D., Li, J., Russell, C. M. P., & Cuadra, J. 2020, MNRAS, 492, 2481 [NASA ADS] [CrossRef] [Google Scholar]
- Yelda, S., Ghez, A. M., Lu, J. R., et al. 2014, ApJ, 783, 131 [Google Scholar]
- Yuan, F., & Narayan, R. 2014, ARA&A, 52, 529 [NASA ADS] [CrossRef] [Google Scholar]
- Yusef-Zadeh, F., Royster, M., Wardle, M., et al. 2020, MNRAS, 499, 3909 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
![]() |
Fig. 1 Comparison of cooling functions Λ(T) for different chemical abundances. The thin and thick lines show the radiative cooling for Solar and three times Solar abundances, respectively. The dashed blue, dotted orange, and dot-dashed green lines represent the cooling functions for atmosphere abundances corresponding to WR subtypes: WC8-9, WN5-7, WN8-9 and Ofpe/WN9 based on the compositions compiled shown in Table 1. The contributions of each element to the radiative cooling were taken from the plasma models by Štofanová et al. (2021). Additionally, the cooling function used in Ressler et al. (2018) is shown as a thick black line. |
In the text |
![]() |
Fig. 2 Comparison of the simulations with different chemical abundances at two different simulation times. The panels show projected density maps weighed by density along the z axis, i.e. ∫ ρ2dz/ ∫ ρdz, which is parallel to the line of sight. Top and bottom panels show the systems at t = −0.7 kyr and t = 0 (present), respectively. Left- and right-hand side panels show the runs WR_f07 and WR_f1, respectively. All maps display the full computational domain. |
In the text |
![]() |
Fig. 3 Net mass-flow rate across a sphere of radius 5rin = 5 × 10−4 pc (1.25 × 10−2″). At this radius, the direction of the mass flow is inwards throughout the entire simulation. The models WR_f07 and WR_f1 are shown in solid blue and orange lines, respectively. |
In the text |
![]() |
Fig. 4 Radial profiles of time-averaged volume-weighted density (top) and mass-weighted temperature (bottom) over the last 500 yr of simulation time. The models WR_f07 and WR_f1 are shown in solid orange and dashed blue lines, respectively. |
In the text |
![]() |
Fig. 5 Radial profiles of the mass inflow Ṁin (solid lines) and outflow Ṁout (dashed lines) rates averaged over the last 500 yr of the simulation. The models WR_f07 and WR_f1 are shown in solid orange and blue lines, respectively. |
In the text |
![]() |
Fig. 6 Orientation of the angular momentum of the gas enclosed in a sphere of radius 0.01 pc (~0.25″). The vertical dimension represents inclination i, while the horizontal dimensions stand for the longitude of the ascending node Ω. Thus, a face-on star orbiting clockwise on the sky would be at the north pole of the graph. Top and bottom panels show the runs WR_f07 and WR_f1, respectively. Each dot corresponds to the analysis of a single snapshot. As a reference, the orientation of the orbital angular momentum of the mass-losing stars is shown as black star symbols, and the grey shaded region corresponds to the average direction of the clockwise disc at (104°, 126°) with a dispersion of 16° (Yelda et al. 2014). |
In the text |
![]() |
Fig. 7 Radial profiles of the scalar tracer fields that represent the WR subtypes: WC89 (dashed lines), WN57 (solid lines), and WN89/Ofpe (dot-dashed) averaged over the last 500 yr of the simulation. The models WR_f07 and WR_f1 are shown in solid orange and blue lines, respectively. |
In the text |
![]() |
Fig. 8 Density and line-of-sight velocity maps of the central 2.5″× 2.5″ of the model WR_fl at t = 0. The left- and right-hand side panels show projected density weighed by density and line-of-sight velocity weighed by mass, respectively both integrated along the z-axis, which is parallel to the line of sight. |
In the text |
![]() |
Fig. 9 Line-of-sight integrated maps of the emission of the H30α and Brγ recombination lines shown on the left- and right-hand side panels. The maps correspond to the inner 2″× 2″ of the model WR_f1. |
In the text |
![]() |
Fig. 10 Synthetic X-ray spectra computed from the numerical models in an annulus of 0.5″ < p < 1.5″. The models WR_f07 and WR_f1 are shown in solid orange and dashed blue lines, respectively. Additionally, the thin solid black line corresponds to the results from the SPH model (Balakrishnan et al. 2024b). |
In the text |
![]() |
Fig. 11 Density and line-of-sight velocity maps of the central 2.5″× 2.5″ of the SPH model at t = 0. The left- and right-hand side panels show projected density and line-of-sight velocity, respectively. These panels are analogous to the ones shown in Fig. 8. |
In the text |
![]() |
Fig. A.1 Comparison of the final state (present time) of the simulations with different Riemann solvers. Left-and right-hand side panels correspond to runs with Harten-Lax-van Leer-Contact (HLLC) and exact Riemann solvers, respectively. Both models were run using the MinMod slope limiter. |
In the text |
![]() |
Fig. B.1 Comparison of the simulations with Solar composition varying the metallicity at two different simulation times. The panels show projected density maps weighed by density along the ɀ axis, i.e. ∫ ρ2dɀ/ ∫ ρdɀ, which is parallel to the line of sight. Left-hand side, central, and right-hand side panels show the runs A1, A3, and A5, respectively. All maps display the full computational domain. |
In the text |
![]() |
Fig. B.2 Radial profiles of time-averaged volume-weighted density (top) and mass-weighted temperature (bottom) over the last 500 yr of simulation time. The models A1, A3, and A5 are shown in solid orange, dashed blue, and dashed-dotted green lines, respectively. |
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.