High-soft to Low-hard State Transition in Black Hole X-Ray Binaries with GRMHD Simulations

To understand the decaying phase of outbursts in the black hole (BH) X-ray binaries (BH-XRBs), we performed very long general relativistic magnetohydrodynamic simulations of a geometrically thin accretion disk around a Kerr BH with slowly rotating matter injected from outside. We thoroughly studied the flow properties, dynamical behavior of the accretion rate, magnetic flux rate, and jet properties during the temporal evolution. Due to the interaction between the thin disk and injected matter, the accretion flow near the BH goes through different phases. The sequence of phases is: soft state → soft-intermediate state → hard-intermediate state → hard state → quiescent state. For the accretion rate (and hence the luminosity) to decrease (as observed) in our model, the mass injection should not decay slower than the angular momentum injection. We also observed quasiperiodic oscillations (QPOs) in the accretion flow. Throughout the evolution, we observed low-frequency QPOs (∼10 Hz) and high-frequency QPOs (∼200 Hz). Our simple unified accretion flow model for state transitions is able to describe outbursts in BH-XRBs.


INTRODUCTION
The black hole (BH) X-ray binaries (BH-XRBs) are primarily dormant, with occasional outbursts triggered by a significant increase in accretion activity.In BH-XRBs, a hardness-intensity diagram, also known as the 'Q' diagram, is typically used to characterize an outburst as well as the features of its various spectral states (hard state, hard/soft intermediate state, and soft state) (see for detail, Fender et al. 2004;Remillard & McClintock 2006;Done et al. 2007;Dunn et al. 2010;Belloni 2010;Belloni & Motta 2016).
In the soft state, the thermal component dominates the spectrum, which is explained by a geometrically thin but optically thick accretion disk in the vicinity of the BH (Shakura & Sunyaev 1973;Novikov & Thorne 1973).A power-law component dominates the spectrum in the hard state, and on the contrary, the contributions from Corresponding author: Indu K. Dihingia ikd4638@sjtu.edu.cn,ikd4638@gmail.com both thermal and power-law components can be seen in the intermediate state (e.g., Homan & Belloni 2005;McClintock et al. 2006;Belloni 2010).The power-law component is caused by the Compton upscattering of soft X-ray photons due to the hot electrons present in the corona close to the BH (Thorne & Price 1975;Sunyaev & Truemper 1979;Chakrabarti & Titarchuk 1995).
In the last few decades, the consensus has been that the combination of hot accretion flow and standard thin-disk in different degrees can explain different states of the outburst quite nicely (e.g., Narayan & Yi 1994, 1995;Esin et al. 1997;Narayan & McClintock 2012;Chakrabarti & Titarchuk 1995;Chakrabarti 2018).However, it is still unknown how the different states are connected and how or why they transition from one to another.Realistic simulations of accretion flow may help understand the state transitions.Nonetheless, simulation for a realistic timescale (hours to tens of days) of an outburst is challenging, especially using general relativistic magneto-hydrodynamics (GRMHD).The framework of GRMHD is indispensable for studying physics around BHs (see, Mościbrodzka 2019; Davis & Tchekhovskoy 2020;Mizuno 2022;Dexter et al. 2021;Dihingia et al. 2022, references therein).
In this work, we report a unified time evolution model for the decaying phase of the outburst in BH-XRBs by GRMHD simulations.Our model can show continuous transitions of spectral states from one to another in the decaying part of the outburst.In the decaying phase, we can safely neglect the contribution of radiative cooling to the structure of the accretion flow.Radiative cooling is more important in triggering the formation of the soft state (SSD) in the rising phase (see Meyer et al. 2007;Das & Sharma 2013;Wu et al. 2016;Liu & Qiao 2022).
In section 2, we briefly discuss the numerical setup, and in the following sections (3,4), we discuss the results obtained from our simulations.Finally, in section 5, we summarise our results and discuss future plans.

NUMERICAL SETUP
This work performs a set of axisymmetric (2D) ideal GRMHD simulations using the GRMHD code BHAC (Porth et al. 2017;Olivares et al. 2019) in Modified Kerr-Schild coordinates.We apply a spherical polar grid (r, θ, ϕ) where the grid spacing is logarithmic in the radial direction and linear in the polar direction.The simulations use a generalized unit system with G = M BH = c = 1, where G, M BH , and c are the universal gravitational constant, the mass of the BH, and the speed of light, respectively.The distance and time are in units of r g = GM BH /c 2 and t g = GM BH /c 3 , respectively.Here simulations are performed considering M BH = 10M ⊙ , for that, the simulation time unit t g corresponds to 4.91 × 10 −5 sec.
We set up an unmagnetized geometrically thin disk in rotational equilibrium around a Kerr BH (spin parameter a = 0.94) initially (Dihingia et al. 2021;Dihingia & Vaidya 2022).It is based on the standard thin disk model by Novikov & Thorne (1973).The density distribution of the thin disk in Boyer-Lindquist coordinates (r, θ) is given by, where z = r cos(θ).The density distribution is kept in a thin disk configuration by choosing α = 2, and an appropriate disk height profile (H) following Riffert & Herold (1995) and Peitz & Appl (1997).Here ρ e (r) is the density profile of the standard thin disk in the equatorial plane (Eq.7 of Dihingia et al. (2021)).As seen in Eq. 1, the initial thin disk is exponentially decaying to an external medium at r out (= 30r g , fixed).The external medium is set by following floor treatment for GRMHD (e.g., Porth et al. (2017)).For the initial velocity profiles of the thin disk, we set u r = 0, u θ = 0 and provide initial azimuthal velocity u ϕ Kep following Dihingia et al. (2021) (see Eq. 12 of the reference).
We consider that the accretion flow is subjected to radiative cooling following Dihingia et al. (2023a,b).However, for simplicity, we do not consider electron thermodynamics in this study.Accordingly, we apply the single temperature approximation (T = m p c 2 Θ/k B , Θ = p/ρ) to calculate the radiative cooling rates.We consider Bremsstrahlung (Q br ) and synchrotron (Q cs ) radiative cooling processes to be dominated in the optically thin regime.With this, the total radiative cooling rate is calculated as: where the Compton enhancement factor (η), which is calculated following (Narayan & Yi 1995).We do not provide the explicit expression of the cooling rates to avoid repetition.The detailed expression can be found in Esin et al. (1996); Dihingia et al. (2023a,b).
In a thin disk, the optically thick components of the radiative losses can not be neglected due to the high optical depth.Therefore, we consider generalized cooling formula capable of choosing radiative cooling processes depending on the optical depth as suggested by Narayan & Yi (1995) and Esin et al. (1996), where τ = τ es + τ abs , τ es = 2σ T n e H, and τ abs = (HQ − thin /σ T T 4 ).Here, H is the scale height of the accretion flow, which is obtained as H = T 4 /|∇T 4 | (Fragile & Meier 2009).Finally, n e is the number density of the electrons, and σ T is the Thomson cross-section of the electron.
Along with the thin disk around the BH, we inject matter continuously from the region away from the BH with a radius r = r inject .The properties of the injected material are considered as ) where, t scale,ρ and t scale,u ϕ define the timescales of the injection for density and angular velocity, and f 0 (= 0.4, fixed for current work to facilitate sub-Keplerian injection of matter) denotes the fraction of the ratio of the azimuthal velocity of injected matter with respect to the Keplerian velocity at the given radius.For this pilot study, we use t scale,u ϕ = 50, 000 t g = 2.455 M BH /10 M ⊙ sec, and choose three values of t scale,ρ = 50, 000 t g , 75, 000 t g , and 100, 000 t g , i.e., 2.455, 3.68, and 4.91 M BH /10 M ⊙ sec.We label them as MOD1, MOD2, and MOD3, respectively.We perform a simulation of MOD1 with a higher resolution (MOD1HR).Furthermore, we choose, ρ in as, (5) For the current study, we choose α in = 8 and r inject = 400 r g .We consider the injected material to be magnetized, and for simplicity, it contains a poloidal singleloop magnetic field which is given by the following vector potential, The strength of the initial magnetic field is set by a minimum value of plasma-β (= p gas /p mag ) at the injection region as β in = p in gas /p in mag | min = 200.Note that the choice of the parameters for this study is not unique.We expect similar qualitative results for different sets of parameters in similar ranges.We discuss the time evolution of the initial setup without the injected matter and only with injected matter (without the thin disk) in Appendix A, which justifies our choice of initial setup and the parameters.
The simulation domain is r ∈ [0.96 r g , 1000 r g ] and θ ∈ [0, π].We perform two simulations with different maximum effective resolutions of 640 × 288 (base resolution 320×144; MOD1-3) and 1280×640 (base resolution 640×288; MOD1HR).To resolve the equatorial thin disk better, the maximum resolution is concentrated near the equator within θ = π/2.4 to θ = π − π/2.4.Simulations run up to time t = 10 6 t g = 49.1 M BH /10M ⊙ sec.Since we need to evolve our simulations for a long time, even our higher resolution simulation cannot afford to adequately resolve the magnetorotational instability (MRI) in the later time due to the weakening of the magnetic field strength over time.However, in the initial phase, we resolved the fastest-growing MRI wavelength, at least with quality factor Q θ ≳ 6 for t ≲ t scale,u ϕ at the disk surface, where magnetized injected matter and unmagnetized thin-disk interacts.

FLOW STRUCTURE
In this section, we discuss the accretion flow structures in different phases of evolution for MOD1.At time, t = 0.25, 2. 455, 12.27, 24.55, and 49.1 M BH /10 M ⊙ sec, we plot the time-averaged density and temperature (p/ρ) distribution in Fig. 1.The time average is performed within 0.1 sec around the marked time.The white red (upper)/green (lower), and grey solid lines in the figure show the Bernoulli parameter −hu t = 1, magnetization σ = b 2 /ρ = 1, and Lorentz factor γ = 1.1 contours.In addition to these panels, follow the supplementary material, where we show a continuous transition of flow structure up to the end of the simulation for MOD1.
Initially, we started with a cold (low temperature) and geometrically thin accretion disk around the BH, which can be seen in panels (a and f) of Fig. 1.In this phase, the matter of the disk has Keplerian angular momentum.From the start of simulations, the injected matter interacts with the thin disk, creating a hot corona surrounding the thin disk.The matter in the hot corona has lower angular momentum than the high-density thin disk.We show the density and temperature snapshots at this time range in panels (b and g) of Fig. 1.This suggests a Keplerian thin disk sandwiched by sub-Keplerian components on both sides of the equatorial plane.Such a structure of accretion flow is usually identified as a two-component accretion flow (TCAF).In such a flow, most of the radiation comes from the thin disk at the equatorial plane.However, due to the hot corona surrounding the thin disk, some hard X-rays are expected to be observed by the inverse-Compton process.Therefore, such a TCAF may explain the soft-intermediate states of an outburst (Chakrabarti & Titarchuk 1995;Chakrabarti 2018, and references therein).Unlike the pure-thin disk case, in this stage, we see outflows from the disk (see −hu t > 1 and γ > 1.1 region around the polar region).In this phase, a weak jet is launched.This weak jet region becomes stronger during the transition to the next phase of the evolution (see section 4 for detail).
With further evolution, the thin disk loses its angular momentum due to the onset of magnetized wind, and matter in the accretion flow heats up.The wind redistributes the high-angular momentum matter far from the BH and creates a very large torus (see panels (c and h) of Fig. 1).It suggests that the density of the flow in most of the region is very low compared to the earlier two phases (panels (a) and (b) of Fig. 1).On the contrary, the temperature of the flow in most of the region is very high compared to the earlier two phases (panels (f) and (g) of Fig. 1).These facts strongly suggest a flow behavior in the hot accretion flow regime (e.g., advection-dominated accretion flow (ADAF)).In this phase, we see a well-evolved steady jet, which is more compact and collimated than in the earlier stages (see the −hu t > 1 and σ = 1 contours in the panel).In this phase, the whole accretion flow acts as a hot corona around the BH.Therefore, we expect the dominant ra- diation component to be mostly hard X-rays produced from synchrotron self-Comptonization and Comptonization of bremsstrahlung photons.However, the optical depth of the flow close to the BH may be slightly higher due to the large size of the torus.This may contribute to some soft X-ray components in the spectra.Accordingly, this phase may be associated with the hard-intermediate state of an outburst.
After creating a large torus, the matter in the torus keeps reducing due to accretion, and the torus becomes smaller with time.We show the corresponding density and temperature distribution in panels (d and i) of Fig. 1.The maximum density in the torus decreases by two orders of magnitudes, which keeps reducing at the end of simulation time.
The reduction of torus size and density decreases the effective optical depth of the flow close to the BH.Accordingly, we may expect an abundance of hard X-rays in the emission spectra.However, the total luminosity will start to decrease as the total mass of the accretion flow drops with time.In this phase, the jet becomes much narrower than in earlier stages (see the −hu t > 1 and γ = 1.1 contours in the panel).This phase resembles a radiatively inefficient quiescent state (RIAF) (e.g., Narayan et al. 1996Narayan et al. , 1997;;Hameury et al. 1997).Thus, we observe that the accretion flow makes a smooth transition from a cold (SSD) disk to a hot torus structure during the evolution.

TIME SERIES ANALYSIS
Along with the flow structure, it is also important to study the time series from the runs to characterize different states further.Following this, Fig. 2 presents the time evolution of magnetic flux (Φ H ) at the horizon for MOD1 (black), MOD1HR (red), MOD2 (blue), and MOD3 (green), (for detail definition, see Eq. ( 96) of Porth et al. (2017)).
We can classify the evolution of the flow into magnetic-driven high-angular momentum flow (MD-HAF) and infall-driven low-angular momentum flow (IDLAF) phases by looking at these evolutions.In the MDHAF, high angular momentum is transported outwards due to active magnetic winds (e.g., Dihingia et al. (2021)), causing accretion towards the BH.During this phase, the magnetic flux increases with time.For MOD1, we observe that up to t ∼ 10 M BH /10 M ⊙ sec (see the black vertical line), the magnetic flux increases.The maximum magnetic flux observed for MOD1 is Φ H ∼ 0.04.After that, magnetic flux remains roughly constant, which we identify as IDLAF.However, we observe slight drops in the magnetic flux within time t ∼ 25 − 40 M BH /10 M ⊙ sec.
By comparing the magnetic flux evolution for two different numerical resolution cases, MOD1 and MOD1HR, we find that they follow a similar increasing trend initially and remain roughly constant after time t > 14 M BH /10 M ⊙ sec.The transition is shown by the vertical red line in the figure.However, we see distinct quantitative differences between the values of the magnetic flux at both resolutions.The maximum magnetic  flux for MOD1HR is about Φ H ∼ 0.08.Accordingly, we expect a stronger relativistic jet in MOD1HR as compared to MOD1.We discuss this in the next section in detail.Interestingly, unlike the low-resolution model, we do not observe any drops in between the evolution.
Comparing simulation models with different t scale,ρ = 50, 000, 75, 000, 100, 000 t g (MOD1, MOD2, MOD3) in Fig. 2, we find that magnetic flux profiles evolve slightly differently for them.Although they show the transition from MDHAF to IDLAF in exactly a similar manner and at the same location, i.e., t ∼ 10 M BH /10 M ⊙ sec.For a higher value of t scale,ρ , more low-angular momentum matter is supplied from the injection radius.It facilitates faster mixing of injected matter and the high angular momentum matter in the thin disk.Interestingly, we note that the magnetic flux drops significantly with time after t > 25 M BH /10 M ⊙ sec for MOD2 and MOD3.The magnetic flux at the end of the simulation is lower for a higher value of t scale,ρ .The average value of magnetic flux for MOD1, MOD2, and MOD3, at the end of simulation (t = 10 6 t g ) is about Φ H ∼ 0.040, 0.019, and 0.017, respectively.This suggests that t scale,ρ plays a role in deciding the characteristics throughout the simulation as well as the final stage.
Next, we would like to study the accretion rate as a function of time.To do that, we plot the absolute value accretion rate for MOD1 at different radii | Ṁ50 |(r = 50r g ), in panels Fig. 3(a)-(d), respectively.In panels Fig. 3(a) and (c), we also plot the same evolution for the high-resolution run (MOD1HR) for comparison.In Fig. 3(b) and (d), we also plot the accretion rate evolution for different values of t ρ , i.e., for models MOD2 and MOD3 for comparison.For MOD1, we observe that for higher radii (r = 20, 50r g ), the accretion rate always shows a decreasing trend.However, at radius r = 10, the accretion rate at t = 20 − 25 M BH /10 M ⊙ sec shows an increasing trend, but after that, it drops monotonically.On the contrary, at a smaller radius (r = 5r g ), the initial accretion rate is smaller but the trend is mostly the same.
Ideally, the accretion rate in the thin disk is supposed to be very high due to the viscosity present in the thin disk.Such viscosity could be triggered in a thin disk due to magneto-rotational instabilities (MRI) if a weak magnetic field is present in the flow (Balbus & Hawley 1998).We performed some 3D test simulations with a magnetic field and presented the results in Appendix B. In these simulations, we observed three orders of magnitude higher accretion rate in the thin disk with resolved MRI as compared to unresolved MRI.The numerical cost of such simulations is extremely high.On top of that, we require simulations of temporal length t ∼ 10 6 t g , which is impossible to achieve with current numerical resources.Apart from this caveat, axisymmetric consideration is the most efficient framework to perform longer simulations.Accordingly, we neglect any such viscosity of the flow and continue with axisymmetric consideration in this paper.Accordingly, we observe that the accretion rate close to the black hole is smaller at the beginning of the simulation because it is in equilibrium.
Note that, during the transition from MDHAF to ID-LAF, the accretion flow makes a transition from a cold thin disk to a hot thick disk (torus) (see Fig. 1).Accordingly, the accretion rate of the thin disk due to the viscosity should also vanish (not studied here).Although Fig. 3(a-b) shows a slightly different shape of the ac-cretion rate evolution with time, it is expected to be half of the bell-shaped curve after the accretion of the thin disk component, with its maximum at the beginning.Such accretion rate evolution can be clearly seen at radii r = 20, 50r g (see Fig. 3(c-d)).This behavior of the accretion rate evolution is expected for the accretion rate in the decaying phase of an outburst (Esin et al. 1997;Ferreira et al. 2006;Kylafis & Belloni 2015).
By comparing the evolution of model MOD1 (black) and MOD1HR (red) in Fig. 3(a-c), we observe that they behave mostly similarly with some quantitative differences.The primary driver of accretion in our simulation is the interaction between unmagnetized highangular momentum and magnetized low-angular momentum matter.This process creates the magnetized wind, which transports angular momentum and leads to accretion in the low-resolution models.Apart from this process, the high-resolution model (MOD1HR) subjected to resolved MRI at least t ≲ t scale,u ϕ .Due to that, the high-resolution model accumulates higher magnetic flux, which consequently results in stronger winds from the thin disk compared to the low-resolution model.Notably, these differences are most prominent near the black hole (see Fig. 3a).However, farther from the black hole, the differences in the evolution become negligible (see Fig. 3c).Additionally, as the matter in the thin disk depletes over time, the disparity between the accretion evolution of MOD1 and MOD1HR diminishes, even close to the black hole.
Finally, we want to compare the accretion rate evolution for different simulation models with different t scale,ρ .As we mentioned earlier, a higher value of t scale,ρ ensures more injected low angular momentum.In panel Fig. 3(d), we observe that the accretion rate for all the models shows an overall decrease with time.However, for the higher value of t scale,ρ (MOD3), the decrease in accretion rate is slower as compared to the lower value of t scale,ρ (MOD1).On the contrary, in Fig. 3(b), the accretion rate at r = 10r g for MOD2 shows a more or less constant profile and slightly decreases towards the end, unlike MOD1, where the accretion rate shows a roughly decreasing tendency throughout.However, in model MOD3 with the highest value of t scale,ρ , the accretion rate evolution shows an increasing trend with time.It is needless to mention that with more time the accretion rate will certainly decrease with time.We observe a decaying trend of accretion rate in the decaying phase of an outburst (Esin et al. 1997;Ferreira et al. 2006;Kylafis & Belloni 2015).Therefore, it implies that higher values of t scale,ρ than those of t scale,u ϕ , which leads to an increase in accretion rate (and hence luminosity), are not preferable to explain commonly known features of outburst in BH-XRBs.The equal decaying time-scale for angular momentum (t scale,u ϕ ) and mass (t scale,ρ ) is preferred.Accordingly, in the next two sections, we will only focus on MOD1 and MOD1HR for the analysis.

Jet properties
Next, we discuss the properties of the jet in terms of the jet radius and the Lorentz factor.We consider any outflow with a Lorentz factor γ > 1.1 (or jet velocity v j ∼ 0.42 c) as a jet for simplicity.With this definition, Fig. 4 shows jet radius and Lorentz factor (γ) as a function of time for MOD1 and MOD1HR.The black, red, and blue lines in the figure refer to the quantities calculated at radius r = 50, 100, and 200 r g , respectively, for MOD1 model.Subsequently, we also plot the same quantities for the high-resolution (MOD1HR) run at r = 200 r g in green for comparison.We first explain, the jet properties for MOD1, and after that we will compare that with MOD1HR.
We do not see any jets initially because the magnetic flux accumulated around the BH is not high enough to launch a jet.After t > 2 M BH /10 M ⊙ sec, we start to see a weak but very broad jet when magnetic flux starts to accumulate on the BH (see panel (c) of Fig. 3).With time, the jet radius monotonically decreases, suggesting the jet becomes more narrow and compact with time.Although, within time t ∼ 10 − 42 M BH /10 M ⊙ sec, we observe a slight increase in the jet radius, but eventually the jet radius drops significantly after t > 42 M BH /10 M ⊙ sec.However, the Lorentz factor does not show a monotonic behavior with time.It roughly follows the magnetic flux evolution as shown in Fig. 2. Initially, the Lorentz factor increases with the increase of magnetic flux, giving a maximum of γ ∼ 1.5 for MOD1 at r = 200r g .After flow changes its nature from MD-HAF to IDLAF (t ∼ 108 M BH /10 M ⊙ sec), the Lorentz factor starts to decrease with a minimum value γ ∼ 1.3 for MOD1 at all radii.After t > 42 M BH /10 M ⊙ sec, the Lorentz factor starts to increase again but in an episodic fashion.The jet radius keeps decreasing with time.Note that the compact jet is expected in the radiatively inefficient hard state (see for discussion Kylafis et al. 2012).
Comparing the results from high-resolution (MOD1HR) and low-resolution runs (MOD1) in Fig. 4, we see pretty similar behavior in the Lorentz factor.In the higher resolution run, the magnetorotational instabilities (MRI) could be resolved marginally, which helps in stronger disk winds than that of MOD1.This affects the maximum magnetic flux regions, where we observe the strongest jet.Accordingly, the maximum Lorentz factor for MOD1HR is about γ ∼ 1.9, which is much larger than that of the maximum value for model MOD1 (γ ∼ 1.4).This shows that the jet physics depends on the resolution of the simulation, and therefore, higher resolution runs are necessary to capture physics related to the jet.On the contrary, the jet radius initially starts with exactly similar evolution, but after some time, we observe a higher increasing trend in comparison to the low-resolution run (blue line).This is due to the fact that the numerical viscosity reduces with resolution.Accordingly, the size of the torus reduced slowly in comparison with the low-resolution model (MOD1).Also, magnetic flux is much higher for MOD1HR, which can expand the jet with time due to the reduction of the supports from the sheath region.Eventually, we expect the jet radius to decrease with time after a sufficiently long simulation time.

Oscillations in the accretion flow
Due to the instabilities between infalling matter and the outflow, the accretion flow exhibits some oscillations.To analyze the oscillations in more detail, we calculate the power density spectra (PDS) of accretion rate at different radii r = 5r g , 20r g , 40r g and 50r g for MOD1 at simulation time t = 2 − 2.5M BH /10M ⊙ sec and plot them in Fig. 5, where the corresponding times are written on each panel in units of M BH /10M ⊙ .At this time, the PDS shows clear peaks in all the radii, suggesting the existence of quasi-periodic oscillations (QPO) in the accretion flow.However, the frequency of oscillations (ν QPO ) differs at different radii.The value of the ν QPO is about ∼ 200, 100, 10, and 8Hz ×10M ⊙ /M BH at radii r = 5, 20, 40, and 50r g , respectively.Note that the frequencies of oscillation are neither equal to the local Keplerian frequency nor equal to radial epicyclic frequency.Although, at radii r = 40 and 50r g , they are quite close to the Keplerian frequency (∼ radial epicyclic frequency for r ≫ 1).However, the frequency of oscillation changes with time at a given radius.
Next, we calculate the PDS at different simulation times for MOD1 and we plot them in Fig. 6.For this figure, we only calculate the PDS of accretion rate radius r = 50 ( Ṁ50 ).We observe that the shape of the PDS changes over time as the flow changes its characteristics.The frequency of oscillation also changes with time.These QPO oscillations in the accretion flow can modulate the emission from the accretion flow, resulting in the observed QPOs in BH-XRBs.Note that, in this phase of evolution, the observed QPO frequency will depend on the observed energy of the emission and the dominant region (radii) of that emission range of the disk.phase, the ratio of injected angular momentum to the Keplerian angular momentum at the innermost stable circular orbit (ISCO) λ inj /λ Kep (ISCO) ≫ 1, we observe low-frequency QPOs (∼ 10 Hz).However, during the phase λ inj /λ Kep (ISCO) ∼ 1, we do not observe prominent QPO features.Finally, for the phase with λ inj /λ Kep (ISCO) ≪ 1, we observe high-frequency QPOs (∼ 200 Hz) in the accretion flow.Accordingly, after a sufficiently long simulation time, we expect a similar frequency of oscillation at different radii.To show that, we plot PDS calculated in a similar way as Finally, we would like to see the impact of high resolution for consistency.To do that, we plot PDS at time t = 2 − 2.5 M BH /10 M ⊙ sec Fig. 8(a,b) and at time t = 14 − 14.5 M BH /10 M ⊙ sec (c,d) for model MOD1HR and MOD1 for comparison.The PDS shows excellent consistency in both resolutions.This result also suggests that the production of QPOs in our simulation does not depend on unresolved/partially resolved MRI.On the contrary, the observed QPO would be due to the mixing of low-angular momentum and high-angular momentum matter, which we expect to act similarly in both resolutions.

SUMMARY AND DISCUSSION
In summary, as the injected matter interacts with the thin disk around the BH, the thin cold disk transitions to a two-component accretion flow.With time, more matter is injected from the outer edge.As a result, the accretion flow becomes a geometrically thick hot torus.Eventually, the torus size decreases, and finally, it is left with very low-density hot matter around the BH.Although we have not calculated the spectra throughout the evolution, still looking at the density, temperature, magnetization, accretion rate evolution, and jet properties, we can comment on the possible spectral evolution.
Throughout our simulation, flow transitions as follows: soft state → soft-intermediate state → hardintermediate state → hard state → quiescent state, which exactly resembles the decaying stage of an outburst.Quite fascinatingly, we could detect the jet properties as expected by many observations, and they suggest jets only in the intermediate states (see Belloni et al. 2011;Kylafis et al. 2012;Kylafis & Reig 2018).Also, we observed a distinct variation in QPO frequency throughout the evolution and observed both kinds of frequencies, i.e., low-frequency QPOs (∼ 10Hz) and high-frequency QPOs (∼ 200Hz).
Based on our knowledge, there are no such reports in the literature showing continuous transitions between spectral states with a unified model.However, it is always anticipated that the accretion rate will take a bell-shaped nature with time for an outburst (e.g., Esin et al. 1997;Ferreira et al. 2006;Kylafis & Belloni 2015;Kumar & Yuan 2022).For our simulation, we only supply half of the bell-shaped curve for the injected matter at the outer edge, following our interest in simulating only the decaying phase.Earlier, Das & Sharma (2013), have also shown with hydrodynamic simulations that the transition from a thin disk to a RIAF happens because of mass exhaustion due to accretion (with different injection models and without the jets and QPO signatures).This suggests that the decaying injected matter is the only essential ingredient to trigger the decaying phase of an outburst, which is very generic in nature, and its qualitative features are independent of the detailed modeling.
A longer decaying time scale for mass as compared to angular momentum results in an increasing accretion rate close to the black hole, which is not expected for commonly known decaying phase outburst.Accordingly, we suggest that an equal decaying time scale for mass and angular momentum is more appropriate to explain the decaying phase of commonly known outbursts in BH-XRBs.
Additionally, our high-resolution run shows quite similar results for accretion rates and QPO frequencies.However, for the high-resolution case, we observe a higher accumulation of magnetic flux around the event horizon.Thus, we observe stronger jets in the highresolution run than that in the low-resolution run.This essentially suggests that to handle the physics of the jet precisely, the setup that MRI is well-resolved will be better.On top of that, the asymmetric nature of our simulation may also impact the jet physics significantly.It is well known that asymmetric MRI can not sustain dynamo for the long term (e.g.,anti-dynamo theorem, Cowling (1933)).Accordingly, a full 3D simulation would be much better for understanding the jet.However, performing such a very long, highly resolved 3D run is extremely challenging with current computational resources.A mean-field dynamo approach will be the most efficient way to include in asymmetric calculations, where it is possible to accomplish 3D dynamo effects (e.g., Parker 1955;Krause & Raedler 1980;Sadowski et al. 2015).Nonetheless, with 3D or with a mean-field dynamo, we expect that the morphology of the Lorentz factor or the magnetic flux evolution remains unaltered.However, we expect some changes in the maximum value of the Lorentz factor and the jet radius obtained from the current results.
Finally, we want to mention some additional caveats and the future studies we plan as a follow-up to this work.In this work, we only demonstrate the decaying stage of an outburst with GRMHD simulations.It will be in our priorities to simulate the full cycle of an outburst (including the rising phase).Note that a realistic outburst in BH-XRBs could be of a few hours to months long and rarely even years (for e.g., Yan & Yu 2015), which is impossible to perform under the current framework.Therefore, more quantitative and parameter space surveys are required to correlate the realistic decay time with the physical properties of injected matter (β in , f 0 , t scale,ρ ), the initial thin accretion disk (r out , initial thickness), or the properties of the central BH (a, M BH , etc.).Moreover, throughout our study, we consider a single temperature approximation, which overestimates the temperature of electrons.Therefore, a robust two-temperature approach would be worth trying for self-consistent results (e.g., Sádowski et al. 2017;Dihingia et al. 2023a).However, we expect that the key features of our study will not alter with it.This research is supported by the National Natural Science Foundation of China (Grant No. 12273022) and the Shanghai Municipality orientation program of basic research for international scientists (grant no.22JC1410600).PS acknowledges a Swarnajayanti Fellowship from the Department of Science and Technology, India (DST/SJF/PSA-03/2016-17).The simulations were performed on Pi2.0 and Siyuan Mark-I at Shanghai Jiao Tong University.This work has made use of NASA's Astrophysics Data System (ADS).We appreciate the thoroughness and thoughtful comments provided by the reviewers that have improved the manuscript.To assess the robustness of our initial setup, we ran two simulation setups: one without any injected matter and another with only injected matter, excluding the initial thin disk.The average density profiles in between t = 4.5 − 5.0 × 10 4 is shown in panels (a) and (b) of Fig. 9, respectively.This simulation time corresponds to t ∼ 2.455 M BH /10 M ⊙ sec.Note that the result of the combination of the thin disk and the injected matter is presented in Fig. 1(b).Without magnetized injected matter, the thin disk can not transport angular momentum and maintain its thin disk structure throughout the simulation.On the other hand, without the thin disk, the injected matter forms a thick torus around the black hole, as shown in Fig. 9(b).The size and shape of the torus are influenced by the value of f 0 , which is 0.4 in our test.This value determines the angular momentum of the injected matter.A higher f 0 value means the matter takes longer to reach the black hole, and the developed torus position becomes farther away.If we kept injecting similar matter over a long period, radiative cooling might lead to the formation of a thin disk from the thick torus.This process represents the rising phase of an outburst in BH-XRBs.However, due to simulation time limits, we added an initial thin disk near the black hole to represent only the declining phase of the outburst in BH-XRBs.We plan to explore the full outburst cycle in future studies.We used a normalization factor of density (0.1) in equation 5, ensuring that the injected density is much lower than the density of the thin disk.Otherwise, the thin disk will be depleted very fast.With this constant, the density maximum for the torus in Fig. 9(b) is three orders of magnitude lower than that of the density maximum of the thin disk in Fig. 9(a).Our main goal is to show a qualitative representation of the outburst's decaying phase, so we kept this factor constant.

B. ACCRETION RATE WITH VISCOSITY
To study the effects of viscosity in the accretion disk, we consider three 3D runs of the thin-disk setup with different resolutions: 1. 160 × 64 × 64 (LR) 2. 320 × 128 × 128 with two refinement level(MR) and 3. 640 × 256 × 256 with two refinement level (HR).For these runs, we put a weak loop magnetic field inside the thin disk with β min = 10.Also, to reduce the extent of the thin disk by reducing the value of r out = 20, which is r out = 30 for the 2D models.Since we only wanted to study the evolution of the disk and therefore, we simulate within θ = π/3 to θ = 2/3π and maximum resolution is concentrated within θ = π/2.4 to θ = π − π/2.4 (same as 2D models).Due to high numerical cost, we only perform simulation up to t = 3000t g .In Fig. 10 and 11, we compare the accretion rate at the event horizon ( ṀH ) and the density profile for the three models (LR, MR, and HR).For these models, the average value of MRI quality factors (< Q θ >, < Q r >, < Q ϕ >) are of the order of ∼ 1(LR), ∼ 5 − 6 (MR), and ∼ 10 (HR) close to the black hole r < 10, respectively.Note that to resolve MRI, the quality factor must be Q ∼ 6 or greater (see Sano et al. 2004).Accordingly, in the model LR, MR, and HR, the MRI is under-resolved, marginally-resolved, and well-resolved, respectively.As a result, angular momentum transport due to MRI is the most efficient in model HR  and inefficient in model LR.This result leads to two orders of magnitude higher accretion rates at the horizon for model HR as compared to LR, which can be clearly observed in Fig. 10.
For completeness, we also show the time average density distribution for models LR(a), MR(b), and HR(c) within t = 2000 − 3000t g in Fig. 11.By comparing these panels, we see that for model LR accretion disk is detached from the event horizon, which is quite similar to our 2D runs (Fig. 1a).For higher resolution (HR), the accretion matter is attached up to the horizon, which supplies the black hole with a consistent accretion rate.

Figure 1 .
Figure 1.Logarithmic time-average density (upper panels) and temperature (lower panels) distribution at different times marked on the top of the figures.The time averaging is performed within 0.1 sec around the marked time.The white, red (upper)/green (lower), and grey solid lines show −hut = 1, σ = 1, and γ = 1.1 contours.

Figure 2 .
Figure 2. Evolution of the magnetic flux at the event horizon (ΦH ) as a function of time.The black, red, blue, and green lines correspond to models MOD1, MOD1HR, MOD2, and MOD3, respectively The corresponding vertical lines indicate the transition from MDHAF to IDLAF.

Figure 3 .
Figure 3. Evolution of the accretion rate at different radii r = 5rg, 10rg, 20rg, and 50rg from top to bottom.In panels (a) and (c), a comparison is shown between MOD1 and MOD1HR, and in panels (b) and (d), a comparison between MOD1, MOD2, and MOD3 is shown.

Figure 4 .
Figure 4. Evolution of jet radius and the Lorentz factor (γ) as a function of time.The black, red, and blue lines in the figure refer to the quantities calculated at radius r = 50, 100, and 200rg.The green line corresponds to the high-resolution run.

Figure 5 .
Figure 5. Plots of PDS of accretion rates at different radii marked on each panel within simulation time t = 2 − 2.5 MBH/10M⊙ sec.

Figure 7 .
Figure 7. Plots of PDS of accretion rates at different radii marked on each panel within simulation time t = 40−40.5MBH/10M⊙ sec.
Fig. 5 but at later time t = 40 − 40.5 M BH /10 M ⊙ sec in Fig. 7.All the panels of the figure show prominent ν QP O ∼ 200 × 10 M ⊙ /M BH Hz.This essentially also suggests that the structure of the accretion flow in a low-angular momentum would be a torus structure at this time.The density maximum of the torus at simulation time t ∼ 40 M BH /10 M ⊙ sec is around r max ∼ 6r g , and the Keplerian frequency of that radius is around ν Kep ∼ 200 × 10 M ⊙ /M BH Hz, which show good resemblance with the ν QP O of the PDS at all radii in Fig. 7.

Figure 10 .
Figure 10.Evolution of accretion rate at the event horizon for model LR, MR, and HR.

Figure 11 .
Figure 11.Time averaged density snapshot in the poloidal plane for model (a) LR, (b) MR, and (c) HR at simulation time t = 2000 − 3000tg.