Bridging relativistic jets from black hole scales to long-term electromagnetic radiation distances: A moving-mesh general relativistic hydrodynamics code with the HLLC Riemann solver

Relativistic jets accompany the collapse of massive stars, the merger of compact objects, or the accretion of gas in active galactic nuclei. They carry information about the central engine and generate electromagnetic radiation. No self-consistent simulations have been able to follow these jets from their birth at the black hole scale to the Newtonian dissipation phase, making the inference of central engine property through astronomical observations undetermined. We present the general relativistic moving-mesh framework to achieve the continuity of jet simulations throughout space and time. We implement the general relativistic extension for the moving-mesh relativistic hydrodynamic code, JET, and develop a tetrad formulation to utilize the Harten-Lax-van Leer Contact (HLLC) Riemann solver in the general relativistic moving-mesh code. The new framework is able to trace the radial movement of relativistic jets from central regions where strong gravity holds all the way to distances of jet dissipation.


I. INTRODUCTION
Relativistic collimated outflows, known as jets, are associated with many astrophysical systems of vastly different scales, from stellar to galactic and even to extragalactic levels.Phenomena like microquasars, young stellar objects, gamma-ray bursts (GRBs), active galactic nuclei (AGN), and quasars demonstrate the prevalence of relativistic jets and highlight the ubiquity of the underlying physical processes that give rise to these phenomena.
A central aspect shared by these varied astrophysical systems is the phenomenon of accretion, in which matter is attracted and pulled into a dense celestial body, like a black hole or neutron star.As matter falls onto these objects, gravitational and magnetic forces play crucial roles in launching and collimating the relativistic jets.Studying relativistic jets across different scales provides astronomers with a unique opportunity to probe fundamental astrophysical processes and test our understanding of high-energy physics in extreme environments.
Commencing with the Penrose process [1,2], numerous theoretical investigations have been undertaken to explore jets and mass outflows near black holes.The Penrose process initially elucidates energy extraction from infalling matter into a rotating black hole.Subsequently, the seminal work by Blandford and Znajek (BZ) demonstrated that jet energy could be extracted from the rotational energy of large-scale magnetic fields surrounding spinning black holes.Later, Blandford and Payne (BP) highlighted that matter could also depart from the surface of the accretion disk due to magneto-centrifugal acceleration.
One of the fundamental questions in accretion disk physics is how the angular momentum transfers within the disk.Initially, Shakura and Sunyaev introduced the 'α-disc' model in a groundbreaking paper.However, the source of the ad hoc viscosity in this model remains questionable.In contrast, recent years have seen widespread acceptance of magneto-rotational instability (MRI; Balbus and Hawley) as the primary mechanism for angular momentum transport in accretion flows.
Another fundamental question in accretion disk physics is the generation of the large poloidal magnetic field as it is pretty natural to assume a toroidal field configuration for accretion flows.To begin with, the orbital differential shear would predominantly amplify the toroidal magnetic field by the shearing of seed poloidal magnetic field, the so-called Ω effect.It took simulators many years to achieve the necessary resolutions and finally report the self generation of the large-scale poloidal magnetic field in black hole accretion disk due to the αeffect (which relies on the buoyancy and Coriolis forces to convert toroidal into poloidal magnetic flux) [7,8].The general mean-field dynamo theory (see, e.g., [9][10][11][12][13]) has been widely used to investigate the generation of largescale magnetic fields from small-scale turbulence.
Recent long-term general-relativistic (GR) neutrinoradiation magnetohydrodynamics (MHD) simulations of the merger of the binary neutron star and black hole neutron star have shown that effective viscous processes, α − Ω magnetic dynamo can lead to the generation of large-scale magnetic field, and post-merger mass ejection [14][15][16].The analysis of the binary neutron star (BNS) merger remnant and post-merger ejecta has been investigated in detail (see, e.g., [17][18][19]).Still, the process of successfully launching a relativistic jet is undoubtedly complex.For a comprehensive understanding of the launching mechanism, general relativistic magnetohydrodynamic (GRMHD) simulations that integrate intricate microphysical processes are imperative.On the other hand, relativistic outflows play a pivotal role in a multitude of astronomical phenomena.For example, it has been speculated that the BNS merger remnants and relativistic ejecta are the central engines of gamma-ray bursts [20][21][22][23] and kilo-nova [24][25][26][27][28][29].Relativistic outflows or jets are instrumental in shaping the emission profiles and contributing significantly to the high-energy radiation observed.Understanding these electromagnetic observations requires tracking the propagation of relativistic jets and their interaction with the ambient medium for a long period of time.However, simulating the complete journey of relativistic jets and the related emission process is numerically challenging.Studies in literature split focus on various parts of the whole process.Many studies conduct MHD/GRMHD simulations to investigate the jet launching process and early propagation (see, e.g., [30][31][32][33][34][35][36][37][38][39][40][41]).Some other studies use special relativistic MHD/HD simulations to investigate the jet's interaction with the ambient medium, away from the central compact region (see, e.g., [42][43][44][45][46][47][48][49][50][51][52][53]).In this study, we propose a formulation to achieve the continuum of jet simulations throughout space and time and potentially bridge these two research domains.The formulation is built upon the development of the moving-mesh technique [54][55][56][57][58][59][60][61][62], which has demonstrated its efficiency in simulating ultrarelativistic jets (see, e.g., [51,63,64]).The extension of the moving-mesh technique to the general relativistic hydrodynamics only appears in recent years.We have seen several moving-mesh codes been extended to include GR effects [60,62,65].Most of these movingmesh codes use Harten-Lax-van Leer (HLL) or Harten-Lax-van Leer-Einfeldt (HLLE) Riemann solver [66,67].However, the HLLC approximate Riemann solver [68] resolves not only the extremal waves but also the contact discontinuity in the Riemann fan and is useful for maintaining contact discontinuities with high precision.Its implementation in fixed-mesh GR employs a local frame transformation [69,70].In this study, we provide the mathematical formulation of incorporating the HLLC Riemann solver into a general relativistic moving-mesh code and demonstrate its robustness in simulating fluid flows under strong gravity.In Sec.II, we implement the general relativistic extension to the special relativistic moving-mesh hydrodynamic code JET [57] using the reference metric formulation [71][72][73][74].In Sec.III, we illustrate the tetrad formulation for solving the HLLC Riemann problem in general relativity and the procedures to incorporate it into the moving-mesh framework.Section IV presents several code implementation techniques.In Sec.V, we conduct several simulations with fixed mesh to test the robustness of the GR extension in the code.In Sec.VI, we conduct numerical tests with the movingmesh grid demonstrating the code's capability to track and resolve the relativistic outflow.For the first time in literature, we successfully launch a relativistic jet from the black hole-torus system and simulate its complete propagation to the dissipation distance.Such simulation provides additional evidence supporting the feasibility of full-time-domain jet simulations, as discussed in our earlier research [64].Conclusions and future work are discussed in Sec.VII.
Throughout this paper, we use the Greek indices (α, β, µ, ν, . . . ) running from 0 to 3 to denote the spacetime components, and the Latin indices (i, j, k, . . . ) running from 1 to 3 to denote the space components.We adopt the geometric units G = c = M ⊙ = 1 throughout this paper.All the length scales and timescales are expressed in units of the gravitational radius r g = GM ⊙ /c 2 and t g = r g /c, respectively, unless stated otherwise.

II. GENERAL RELATIVISTIC HYDRODYNAMICS IN A REFERENCE METRIC FORMULATION
The 2D special relativistic moving-mesh hydrodynamic code JET adopts spherical coordinates assuming axisymmetry.The cell interfaces orthogonal to the radial direction are allowed to move radially.The code is essentially Lagrangian in the radial direction, coupled laterally by transverse flux.This setup is particularly suitable for modeling relativistic radial outflows [55].To minimize the modifications for the code, we derive the general relativistic hydrodynamic equations in a way that resembles the special relativistic counterparts.In the following, we lay out the implementation steps for clarity.Despite of the axisymmetry property of the JET code, throughout this paper we will show all the derivations without imposing any symmetry for completeness.
In the standard 3+1 decomposition (see, e.g., [75][76][77]), the spacetime is foliated by a family of spatial hypersurface Σ t with future-pointing timelike unit normal vector denoted by n µ , which decomposes the line element as where α is the lapse function, β i is the shift vector, and γ ij is the spatial metric induced on Σ t .In terms of the lapse and shift, the normal vector n µ can be expressed as We adopt a conformal decomposition of the spatial metric γ ij where ψ := (γ/γ) 1/12 is the conformal factor, γij is the conformal spatial metric, and γ and γ are the determinants of γ ij and γij respectively.Following the referencemetric formulation (as shown in [78]), we define the residual metric ϵ ij as where γij is a time-independent background reference metric.For our purpose, we specialize γij to be a flat metric in spherical coordinate (r, θ, ϕ) as γij := diag(1, r 2 , r 2 sin 2 θ).To make the conformal scaling unique, we set γ = γ := det(γ ij ) (see, e.g., [79]).We denote ∇ µ , D i , Di and Di as the covariant derivatives of spacetime metric g µν , γ ij , γij and γij respectively.
The equations of relativistic hydrodynamics are based on conservation of rest mass and conservation of energy-momentum where ρ is the rest-mass density and u µ is the fluid fourvelocity and T µν is the stress-energy tensor.Here we assume perfect fluid for T µν in the form where P is the pressure, ε is the specific internal energy and h := 1 + ε + P/ρ is the specific enthalpy.In 3+1 decomposition, T µν can be decomposed as where W := αu t is the Lorentz factor and v i := u i /W + β i /α is the fluid velocity measured by the normal observer.
We adopt the Valencia formulation in reference metric formulation following [72,80] to rewrite the hydrodynamics equations in conservative form as with state vectors q being the conserved variables where (D, S i , τ ) := (ρW, ρhW u i , ρhW 2 − P − D) are the density, momentum density and energy density variables in Valencia form respectively. f j and s represent the flux and source terms respectively written as The detailed derivation is shown in Appendix A for the readers' interests.
One key ingredient of the reference metric method is to evolve tensorial quantities in an orthonormal basis with respect to the background metric.In this way, all tensor components are explicitly free of coordinate singularities.We will follow the notation of [80] to distinguish between coordinate-basis and orthonormal-basis components.The plain Latin indices represent the tensor components in the standard coordinate basis, while the Latin indices surrounded with curly braces denote the components in the background orthonormal basis.We also introduce a set of basis vector ê{k} i that are orthonormal with respect to the background metric γij , For the flat background metric in spherical coordinates, this leads to So any tensor A i j defined in the standard coordinate basis can be decomposed into its orthonormal basis counterpart A {i} {j} as As an example, the residual metric ϵ ij can be expressed in terms of the components in the orthonormal basis ê{k} i as while for the conserved momentum we have (q Sr , q S θ , q S ϕ ) = ψ 6 γ/γ S {r} , rS {θ} , r sin θS {ϕ} .
The complete set of general relativistic hydrodynamic equations in 3D spherical coordinates under reference metric formalism (9) where K ij is the extrinsic curvature, M := ψ 6 γ/γ and v{i} := v {i} − β {i} /α, alongside with the special relativistic Eq. (A25) in Appendix A 4 for comparison (see also [81]).Noted that in practice we compute ∂ k ϵ {i}{j} numerically in source terms instead of ∂ k ϵ ij as while the second term ∂ k ê{l} i ê{m} j is evaluated analytically.

III. TETRAD FORMATION AND THE HLLC RIEMANN SOLVER
To evaluate the numerical flux through cell interfaces, HLL-type (HLLE/HLLC) Riemann solvers have been designed for relativistic hydrodynamics in Minkowski spacetime [68,82].Most of the GRHD/GRMHD codes in the literature use HLLE Riemann solver in curved spacetime (see, e.g.[83][84][85][86]).The HLLC Riemann solver that captures the contact discontinuity in the wave fan has recently been added for GR codes [69,70,87].We follow previous works for the implementation of the HLLC Riemann solver in general relativity [69,70,88].The basic idea is based on the equivalence principle: physical laws in a local inertial frame of a curved spacetime have the same form as in special relativity.When we define such inertial frame, we can then use the solution of Riemann problems in a local Minkowskian frame to construct the corresponding solution in curved spacetime.The previous section derives the general relativistic hydrodynamic equations in a reference metric formulation.For the benefit of the coming discussion, we will revert to the original formulation [84] in this section with g = det(g µν ) satisfying √ −g = α √ γ.The state vector F 0 and the flux vector F i are given by and the source term in this formulation is denoted by S.
Since the source is irrelevant to the tetrad formulation in following discussions, we here omit the explicit form of S.
Let us consider a single computational cell of our discrete spacetime Ω, bounded by a closed threedimensional surface ∂Ω.We take the 3-surface ∂Ω as the standard-oriented geometric object made up of two spacelike surfaces {Σ x 0 , Σ x 0 +∆x 0 } plus timelike surfaces {Σ x i − , Σ x i + } that join the two temporal slices together, where x i ± are the cell boundaries of Ω in ±x i directions.The integral form of the system (19) is where is the volume element of cell Ω.From now we will drop the wedge symbol ∧ for simplicity.The integral form (21) can be rewritten in the following conservation form where F0 x 0 is the volume integral of F 0 at x 0 given by F0 x 0 := and F i ± is the integrated spatial flux across the cell interfaces x i ± given by A. Tetrad formulation Instead of attempting a direct resolution of the Riemann problem within the curved spacetime, our approach entails deliberately converting the left and right states at a given interface into a local Minkowskian frame of reference.This methodology enables the utilization of developments in the realm of special relativistic Riemann problems, as proposed by [88,89].
To begin with, we define a new tetrad basis e (μ) that satisfies a list of properties as shown in [69]: 1. e (μ) must be orthogonal to e (ν) for all µ ̸ = ν.
2. Each e (μ) must be normalized to have an inner product of ±1 with itself, with e ( 0) being timelike and e ( î) being spacelike.
3. e ( 0) must be orthogonal to surfaces of constant x 0 .
4. The projection of e ( î) onto to the surfaces of constant x 0 is orthogonal to the surface of constant x i within that submanifold.

Without loss of generality, let us only consider the conversion of the volume integral F0
x 0 in Eq. ( 24) and the first spatial flux integral F 1 + in Eq. ( 23).We define the following tetrad basis in the spherical coordinates with x µ = (t, r, θ, ϕ) (the detailed derivation can be found in Appendix of [69,70]) as where the coefficients are given by The covariant components of the tetrad basis are given by e (μ)µ = g µν e (μ) ν .Specifically The transformation of vector and tensor between the tetrad frame and the original Eulerian observer frame follows and Note that the upper and lower spatial tetrad components are the same Therefore, we can define F(μ) as the tetrad transformation of F µ in the form Here for momentum components ( FS ( ĵ) ) (μ) of F(μ) we need to perform one more tetrad transformation due to its tensorial nature.Since we only focus on the flux along r direction, the components F( t) and F(r) are written as where The inverse transformation is given by which gives In addition, we can reformulate the conservation form Eqs. ( 24) and ( 25) with the tetrad basis.Note that the indexes (t, r, θ, ϕ) and ( t, r, θ, φ) are interchangeable with (x 0 , x 1 , x 2 , x 3 ) and (x ( 0) , x ( 1) , x ( 2) , x ( 3) ) respectively.Making use of the following invariance property and transformation rule we can get (see also [84]) where ĝ := det(η ij ) = −1.This gives the volume integral of F t (24) and integrated spatial flux of F r (25) in local tetrad basis as with nonzero interface velocity from a nonzero drift in the direction of interest, in agreement with [69,88].
With tetrad basis formulation, the procedure to obtain the numerical flux across the first spatial direction involves the following steps: 1. Obtain the values of the primitive variables (ρ, P, u i ) and tetrad basis e (ν) 2. Construct the conserved variable F( 0) and flux F( 1) for the left and right state in the tetrad frame.
3. Solve the Riemann problem in the tetrad frame with a nonzero interface velocity V (r) interface .
4. Once we have the updated solution of F( 0) and F( 1) , we can obtain the numerical flux across the first spatial direction in the Eulerian observer frame according to Eq. ( 40).

B. HLLC Riemann Solver in the tetrad frame
We solve the Riemann problem in the tetrad frame by adopting a special relativity form.We calculate the HLLC flux F(r) by solving the one-dimensional conservation law [70]: with Given an initial condition at cell interface rj+ 1  2 described by three characteristic waves and four states will be established inside the Riemann fan as and the corresponding numerical flux across interface rj+ 1 2 is F(r) where λL/R is the characteristic speed of the left/right going nonlinear wave and ( F(r) ) Explicitly, we have the left or the right state as To reduce the number of unknowns and have a wellposed problem, we assume that S * (r) = (τ * + P * + D * ) λ * (see [68]).If one defines E := τ + D and performs the calculation of (53b) − λ * × [(53a) + (53e)], one will get the following expression, giving λ * in terms of P * [68]: By imposing P * L = P * R across the contact discontinuity, we find the following quadratic equation for λ * where Once we obtain the speed of the contact discontinuity λ * , P * can be obtained from Eq. ( 54).The conserved quantities in the intermediate states are given by The left and right characteristic speeds λL/R follow Davis's estimate [68] λL = min( λ− ( ÛL ), λ− ( ÛR )), (58) λR = max( λ+ ( ÛL ), λ+ ( ÛR )), (59) with where v 2 = v( î) v( î) and c s is the speed of sound Equivalent expressions for the (θ, ϕ) directions can be easily obtained.In the Eulerian observer frame, the minimum and maximum characteristic speeds λ ± are given by [84,90,91]: Making use of Eq. ( 36) and Eq. ( 60), we can get the following relation: Note that Eq. ( 60) and Eq. ( 62) are derived in the special relativistic and general relativistic setting, respectively.Equation ( 63) establishes their relationship with the tetrad method consistently.
In addition, we implement the HLLE Riemann solver [67,92] for comparison.We adopt the same tetrad formulation.The HLLE Riemann solver is constructed by assuming an average intermediate state between the fastest and slowest waves in the tetrad frame.The two characteristic waves and three states inside the Riemann fan become: The corresponding numerical flux across interface rj+ 1 2 is: where Û * and ( F(r) ) * are the intermediate state and flux.They can be derived from the jump condition [see Eq. 52] as: and C. HLLC Riemann Solver for the Moving-Mesh GR For the moving mesh in the simulation domain, naturally, we need to solve the Riemann problem on the moving interface with its own coordinate velocity V := dx dt = (V r , V θ , V ϕ ).Let us denote the corresponding four-velocity as u µ gridface .In general, when we consider the 3 + 1 spacetime foliation Σ t , we define a unit normal vector as n µ , and this unit normal vector corresponds by definition to the four-velocity of the Eulerian observer [75].When we define the fluid's four-velocity as u µ , the velocity of the fluid with respect to the Eulerian observer (v µ = (0, v i )) has the following relation: where is the Lorentz factor of the fluid with respect to the Eulerian observer.When we move from a given hypersurface to the next following the normal direction, the change in the spatial coordinates is given as [75]: β µ = (0, β i ) being the shift vector.Then v i is related to the coordinate velocity V by v i = 1 α (V i +β i ).In our case, only the cell interface orthogonal to the radial direction can move with a coordinate velocity denoted as V r .Then the four-velocity of our radially-moving interface is In the above, we illustrate the explicit definition of different velocities for clarity.For our moving-mesh code, the grid moves radially, the integral of the radial flux at a short time interval dt becomes With Note that the above velocity equation relates to Eq. ( 46) and Eq. ( 63).From Eqs. (39), and (40), we have: Compared with the tetrad formulation for the static mesh, we replace the interface velocity γ rr to incorporate the effect of the moving interface into the flux integral.
In principle, the coordinate velocity for the moving interface can be set freely.At each instantaneous time, on the cell interface, the three characteristic waves and four states inside the Riemann fan depends only on the values of the primitive variables on the left and right sides of the interface.The interface velocity will influence which state the numerical flux across the interface will be selected [see Eqs.(50], and ( 51) ).Based on this flexibility, we choose the contact discontinuity velocity as the interface velocity: We find this choice performs well for the simulation of ultrarelativistic jets.
For the derivation of the tetrad formulation and HLLC Riemann solver, we express every metric and fluid variable in the coordinate basis.For the implementation, we utilize those variables in the orthonormal basis instead.For example, in the tetrad basis calculation, we will use In this way, the geometric factors will not directly appear in the tetrad basis calculation.The derivation itself remains the same because of the invariance of spacetime interval under coordinate transformation Making use of this invariance principle, we can handle the moving mesh in another way.First, boost the coordinate basis into the comoving coordinate basis of the interface: Second, boost the primitive velocities into the comoving coordinate basis: Third, making use of the invariance, calculate the corresponding metric g ⟨µ⟩⟨ν⟩ components: Once we have the new lapse, shift and spatial metric in the comoving frame, we can derive the tetrad basis in the comoving coordinate basis, and solve the HLLC Riemann problem accordingly.We lay out this approach for readers' interest as well as for a more complete discussion.

IV. NUMERICAL TECHNIQUES A. Implementation of equations
For the numerical implementation, we discretize the volume averages of Eq. ( 9).Using divergence theorem, the discretized version of equation 9 in the cell (i, j) can be expressed as (Since our code is 2.5D, we will ignore the discretization in the ϕ direction.)[85]: where the cell volume and volume average are defined as while the surface area and surface average is defined as Note that when we perform the volumn average ⟨•⟩ or surface average ⟨•⟩ i , we could strip out the geometric factor from the tensorial expressions in coordinate basis and integrate them together with the volume factor √ γ.In this way, the tensorial variables in orthonormal basis like S {θ} become truly independent of the underlining geometry.For example, in the spherical coordinates, the volume average for the conserved momentum q S θ will be calculated as For our moving-mesh scheme, the cells in the radial direction will continuously merge and divide.When we perform the above integral, the variables in orthonormal basis like MS {θ} will be better conserved.As an example, if we assume MS {θ} is constant across cell (i,j) and cell (i+1,j) , when we merge these two cells, it gives a combined conserved momentum as where cell (inew,j) is the combined cells of cell (i,j) and cell (i+1,j) with ∆V inew,j = ∆V i,j + ∆V i+1,j .From the combined momentum, we can recover the variable MS {θ} accurately.
In code implementation, the contribution of the source term to the conserved variables inside a cell is defined as ⟨s⟩ i,j ∆V i,j .We perform volume integral on the singular factors such as 1/r and cot(θ) that appear in the source term [see Eq. ( 17)].Explicitly, the integral of 1/r factor gives (r 2 + − r 2 − )/2 while the integral of cot(θ) leads to (sin(θ + ) − sin(θ − )).This practice turns out to reduce numerical error for the source term calculation near singular points.
Finally, to work out the cell volume, cell surface, we make the following definition

B. Recovery of primitive variables
There are many possible ways to make the conversion between conserved variables and primitive variables (e.g., [93]).Our current research focuses on relativistic jets propagating in an ambient medium.We need to deal with large variations of density and pressure in the jet simulations.The following cons-to-prim method proves to be robust for such a task.We use ρ, P, W v {i} = u {i} + W β {i} /α as our primitive variables where W v {i} is the projected fluid velocity in orthonormal basis.For the equation of state (EOS), we only consider the case of a single-component perfect gas for now.In this case, the specific enthalpy h ≡ 1 + ε + P/ρ is a function of a temperaturelike variable Θ = P/ρ (see [94]).In the literature, the most widely used EOS is the ideal gas EOS:P = (Γ − 1)ρε, where P is the gas pressure, ε is the specific internal energy density.Which can be expressed as: where Γ is the adiabatic index.The ideal gas EOS has been applied to the gas of either subrelativistic temperature with Γ = 5/3 or ultrarelativistic temperature with Γ = 4/3.For our simulations of relativistic outflow propagating in a cold ambient medium, a variable equivalent adiabatic index Γ eq = (h−1)/(h−1−Θ) is desirable to account for transitions between the nonrelativistic and the relativistic temperature regime.There have been efforts to find EOSs that better describe the thermal dynamics of relativistic gas.Synge and Morse derives the correct EOS for the single-component perfect gas in a relativistic regime using modified Bessel functions.Mignone et al.
proposes an approximate EOS (denoted as TM EOS) that is consistent with the Taub's inequality [96]: for all temperatures.It differs by less than 4% from the theoretical value given in [94].Ryu et al. proposes a new EOS (RC EOS), which better fits the theoretical value.
Let us write the expression of the specific enthalpy for the RC EOS: Following the definition of the general form of polytropic index n and the general form of sound speed c s : their values can be calculated for RC EOS as: For both TM and RC, we have correctly c 2 s → 5Θ/3 in the nonrelativistic temperature limit and c 2 s → 1/3 in the ultrarelativistic temperature limit [97].
We can use these expressions to convert the conservative variables into primitive ones with a standard Newton-Raphson method (NRM) [98], using Θ as our independent variable.We will use the (known) values of the conservative variables.
First, by squaring the momentum equation, we get with h = h(Θ) given by the EOS.Using the relation p = DΘ/W , we get the energy density (excluding rest mass), τ = DhW − DΘ/W − D. We can then derive the following identity [98] : Together with Eq. ( 93), the derivative df /dΘ has the form: where the relation W ′ = −h ′ (W 2 − 1)/(hW ) has been used (derived from Eq. ( 93), see also [98]).The derivative dh/dΘ depends on the particular EOS used.We adopt the RC EOS [see Eq. ( 89)] for the simulations of relativistic jets and the ideal gas EOS for the remaining numerical tests.

C. Reconstruction
We reconstruct the primitive variable (denoted with Q ) to the left and right sides of each cell with the total variation diminishing (TVD) method described in [99]: where ∆ξ i = ξ i+1/2 − ξ i−1/2 is the cell width, and ξ is the cell center.And ∆Q i is a slope-limited gradient function written in terms of a nonlinear limiter function φ(v): We adopt the same modified monotonized central (MC) limiter in [99]. where To reconstruct the left and right state of the cell i, the above stop limiter utilizes the cell average values of ⟨Q⟩ i−1 ,⟨Q⟩ i , and ⟨Q⟩ i+1 , defined at the cell center ξi−1 , ξi , ξi+1 .This algorithm takes into account nonuniform spacing.The cell center position ξi can be taken as the volume-averaged cell center ("centroids of volume") or arithmetic-mean cell center.In this study, we adopt the arithmetic-mean cell center for our simulations.

D. Treatment of numerical conditions
Robust numerical simulations require the treatment of several numerical conditions.One of them is the Courant-Friedrich-Levy (CFL) condition [100], which limits the time step size in explicit numerical methods.The simulation domain of the JET code allocates cells at the same temporal level.A global time step will be used to evolve simulation time.To find the global time step, we first calculate the time step δt of individual cells in the domain according to: where CFL is the CFL number.Its value has been set to 0.4 for simulations performed in this study.(λ r ) ± and (λ θ ) ± are again the minimum and maximum characteristic speeds for the cell in the radial and polar direction, respectively.V r cell is the cell's radial velocity which approximates the cell's upper interface velocity.We then pick the smallest time step as the global one.The subtraction of the cell's radial velocity in Eq. ( 103) leads to a much larger time step, making the long-term simulation of relativistic jets computationally efficient.Another numerical condition that needs to be taken care of is the boundary condition.For our cell-centered grid structure in spherical polar coordinates, we follow the boundary treatment described in [73,80].We first allocate two layers of ghost zones for each of the four boundaries (two in the radial direction, and two in the polar direction), and then fill the boundary ghost zones at the radial origin, and at the θ boundary with values copied from the corresponding points in the interior of the grid, accounting for appropriate parity factors.For the outer boundary in the radial direction, we adopt the Dirichlet boundary condition and use the initial data routine to set their ghost zone values.

E. The adjusted moving-mesh scheme
Since the initial development of the JET code [57], the moving-mesh scheme has kept being updated to improve the accuracy and efficiency of relativistic jet simulations.The adjusted moving-mesh scheme in this study contains the following rules: inside the simulation domain, the radial interface of a grid cell will move at local contact discontinuity velocity of the flow.Each radial track moves independently.The inner and outer radial boundaries of the domain can also move.At each time step, the longest and shortest cell in each radial track will be marked for refinement or derefinement according to the maximum or minimum aspect ratio of grid cell (a := ∆r/r∆θ) allowed in the simulation (see [57] for more information).In ultrarelativistic jet simulations, we find the domain cells can squeeze into an ultrathin shell with the cell's aspect ratio reaching 1/100 or even smaller.In order to resolve the relativistic thin shell, only cells with length ∆r < r/(8Γ 2 ) will be marked for derefinement.In addition, we define an approximate second derivative of a fluid variable as a measurement of error E i to mark the region of interest.At each time step, the cell along each radial track with the maximum measurement of error will be marked for refinement if its aspect ratio is larger than twice the minimum aspect ratio and its measurement error E i > 0.9.The cells with E i < 0.002 will be considered for derefinement.The cell to be derefined is the one that has the smallest time step (see [64]).To reduce load imbalance of CPUs, the number of grids in each radial track will be balanced dynamically during the simulation.

V. FIXED-MESH NUMERICAL SIMULATIONS A. Bondi accretion in maximally sliced trumpet coordinates
We first consider spherically symmetric, radial fluid accretion onto a nonrotating black hole (ingoing Bondi flow) [101,102].Following previous work (e.g., [70,103]), we perform simulations of Bondi flow in maximally sliced trumpet coordinates [104,105].The transformation between Schwarzschild coordinate and maximally slicing trumpet coordinate is illustrated as a reference in Appendix B. We set the fluid parameter according to Table 1 of [103]: the accretion rate Ṁbondi = 10 −4 M , the adiabatic index Γ = 4/3, and the critical radius R s = 10M where M is the mass of the central black hole.For simplicity, M is set to 1 in the simulation.
The simulation domain is in an axisymmetric spherical coordinate, spanning the region r ∈ [0.4 M, 10 M ], θ ∈ [0, π/2].We employ logarithmic grid spacing in the radial direction with a cell's aspect ratio a set to one (i.e. a := ∆r/r∆θ = 1).The finest cell, located closest to the inner boundary, has a spacing ∆r = r min (θ max − θ min )/(Nt), where Nt is the number of cells in the azimuthal direction.
To maintain the unity aspect ratio of the cell, the number of cells in the radial direction Nr is calculated as We conduct simulations with three different resolutions: low resolution with Nt = 128, medium resolution with Nt = 256, and high resolution with Nt = 512.For the benefit of convergence test, we set the number of grids in the radial direction Nr = 2 Nt.In this case, the cell's aspect ratio will deviate from one slightly.In Fig. 1, we show the radial profiles of the fluid restmass density (top) and the fluid velocity (middle) at time T = 0 M (init) and T = 50 M (fin) for the medium resolution simulation.The profile of the Bondi flow has been maintained throughout the simulations.In the bottom panel, we plot the L1-norm of error for the rest-mass density.The L1-norm of error is defined as [69] The Bondi simulations demonstrate second-order convergence for the L1-norm of error with respect to the resolution.The code adopts the second-order RK2 time integrator and the second-order piecewise linear reconstruction method (PLM), described in Sec.IV C. The presented convergence result is as expected and agrees with previous studies (see e.g.[69,70]).For the implementation of a higher-order reconstruction scheme for our unstructured grid in spherical geometry, like the piecewise parabolic method (PPM) [106], weighted essentially nonoscillatory (WENO) [107][108][109][110], or the monotonicity preserving scheme (MP5) [109], we will refer to future work.

B. Tolmann-Oppenheimer-Volkoff star
The next numerical test we consider is the Tolman-Oppenheimer-Volkoff (TOV) star with the structure of a spherically symmetric body of isotropic material in equilibrium [111,112].We conduct two TOV star tests based on [60]: the stationary case and the one with pressure depletion.The initial profile for the TOV star has a central rest-mass density ρ c (0) = 1.28 × 10 −3 .We adopt the polytropic EOS P (ρ) = Kρ Γ , with (K, Γ) = (100, 2) for the initial data.As for the evolution, we adopt the ideal gas law.Additional parameters for the initial profile can be found in Table I in the cgs unit.
In Fig. 2, we plot the central maximum density variation as a function of dynamical time ( ρ c (0) t) for both cases.For the stationary case, we find the central maximum density varies within 0.5% for 14 dynamical times for the Nt = 64 simulation, confirming the stability of the star.When we increase the resolution to Nt = 128, the result gets better.For the pressure depletion simulation, we reduce the TOV initial pressure profile by ten percent.The star falls out of equilibrium and undergoes radial oscillations.We conduct simulations with two different resolutions (Nt = 64 and 128) and find consistent oscillation pattern, as shown in the bottom panel of Fig. 2. The result is equivalent to the test result in [60].

C. Fishbone-Moncrief torus around a Schwarzschild black hole
Our next test concerns a stationary, axisymmetric, isentropic torus around a Schwarzschild black hole [113].We consider a particular instance of the Fishbone-Moncrief solution where the spin of the black hole is set to zero.
We generate the initial data in the Schwarzschild coordinate with its radius denoted by R.However, we will evolve the system in the isotropic coordinate of the Schwarzschild metric with its radius denoted by r (see Appendix B).The initial profile generator follows the implementation in [79,86,114].Table II shows the key variable values for the torus.For the ambient atmosphere, we set ρ = ρ min (R/R g ) −3/2 , P = P min (R/R g ) −5/2 , where ρ min = 10 −8 , P min = 10 −10 , R g = GM/c 2 is the black hole gravitational radius and M is the black hole mass.For the simulation, we employ an ideal gas EOS: P = (Γ − 1)ρε, with Γ = 4/3.In the azimuthal direction, the simulation domain extends from θ min = π/3 to θ max = 2π/3.In the radial direction, the grid covers the region from r min = 4 to r max = 20.At the location of maximum pressure r = 10.98 (R = 12), the orbital period of the torus is around 238.9.We set the final time of the simulation to be 2000, roughly eight orbits.We conduct two simulations with grid resolution Nt = 256 and Nt = 512, and find consistent results.
Figure 3 illustrates the contour plots of the black holetorus system at the beginning (top panel) and at the end of the simulation (middle and bottom panels), taken from the Nt = 512 simulation for better visual effect.The top panel shows the initial contour plot for the logarithmic density ρ.Comparing these two contour snapshots, we first find that throughout the simulation, the torus maintains its density structure.We check that the maximum rest-mass density always keeps the original value within 4% during the simulation, and its radial position varies within 2%.Because the torus stays close to the black hole, the ambient gas falls into the black hole and blows the torus surface in the infalling process.A bow shock appears in front of the torus and a trailing tail fills in the inner region between the torus and the central black hole.The falling gas slows down when it crosses the bow shock as can be seen from the velocity contour plot.The stability of the torus structure near the black hole showcases the code's robustness in the handling of fluid rotation under strong gravity.

D. Rayleigh-Taylor instability for a modified Bondi flow
Previous work [57] with the original JET code has captured the detailed nonlinear features of Rayleigh-Taylor instability in a relativistic fireball.It uses the HLLC Riemann solver described in [55].To test our general relativistic HLLC Riemann solver, we modify the Bondi flow to induce Rayleigh-Taylor instability under strong gravity.The setup is similar to a Strömgren sphere around the central black hole-the low-density hot gas is surrounded by a high-density gas with gravitational acceleration [116].Within a radius of r < 3(1 + 0.1 × (1 + cos(80θ))/2.0),the density and pressure of the Bondi flow have been modified as ρ = 0.1ρ bondi , P = 50P bondi , v r = 0. ρ bondi and P bondi are taken from the Bondi profile in Sec.V A. This setup creates a hot low-density bubble inside the Bondi flow with a curly interface.As the hot low-density gas pushes against the heavier Bondi flow, Rayleigh-Taylor instability (or sometimes referred to as Richtmyer Meshkov instability in this case) develops.We perform this simulation with an azimuthal resolution of Nt=512, covering the azimuthal angle from π/3 to 2π/3. Figure 4 shows its time evolution.Initially, the hot gas pushes outward and compresses the incoming Bondi flow into higher density as shown at t = 5 M. Instability fingers develop and evolve inside the low-density region.Nonlinear features of the instability continuously evolve at t = 15 M. Later on, due to the attraction of the central black hole, the turbulent gas flows into the black hole.The implemented HLLC Riemann solver is able to capture the detailed structure of the instability in the strong field regime.It performs better than the HLLE Riemann solver.

VI. MOVING-MESH NUMERICAL SIMULATIONS A. Spherical shock tube test
One advantage of our moving-mesh code is that the cell face is able to move with the contact velocity of the flow in the radial direction.It has been shown that the contact discontinuity is much better preserved when employing HLLC on the moving mesh (see Fig. 7 of [55]).What is more, the flow naturally adjusts the cell width in the radial direction.Combined with robust refinement and derefinement schemes, the simulation domain will be able to resolve the region of interest [64].To test the accuracy of the moving-mesh scheme, we conduct the identical spherical shock tube test as shown in [61]: within the radius of 0.25 (r < 0.25), the density ρ and pressure P is set to 1. Outside of this region, the value of density and pressure is 0.1.We adopt the Minkowskian frame for the test.Since the tetrad formulation for the HLLC Riemann solver also works for the Minkowskian metric, we do not take any additional steps for the special relativistic simulations.
In azimuthal direction, the simulation domain extends from 0 to π/2 with Nt = 128.In the radial direction, the grid covers the region from r min = 0.01 to r max = 0.5.We adopt logarithmic spacing in the radial direction and set the initial cell's aspect ratio to one.We first conduct the spherical shock tube test with different Riemann solvers in fixed-mesh simulation.Both the HLLE and HLLC Riemann solver handle the test well and give almost the same results (as shown in Fig. 5).In Fig. 6a, we compare the end profile for simulations with the fixed mesh setup and the moving mesh setup.The density plot exhibits a sharp transition at the contact discontinuity for the moving mesh and a relatively smooth one for the fixed mesh.Following the compression of the fluid in the shocked region, the cells squeeze between the contact discontinuity and the forward shock.The P/ρ plot reveals a jump at the contact discontinuity.We find this appears in the moving-mesh simulation here as well as in the literature [55,61].It may come from the physical squeezing of the fluid as the grid moves together with the flow in the moving-mesh simulation or the TVD reconstruction scheme requires some adjustification for the moving mesh.
To investigate the jump's dependence on numerical resolution, we have performed additional moving-mesh simulations with different numerical resolution: the Nt = 32, Nt = 128 and Nt = 2048 resolution.Results from these three simulations (see Fig. 7) demonstrate that the jump feature persists and its magnitude is invariant under different numerical resolution.We also perform additional fixed-mesh simulations with higher resolution and find no presence of the jump feature in these simulations.Since the jump's magnitude does not increase with time and spatial resolution, considering its minimal impact on the fluid dynamics in moving-mesh simulations, we will leave this numerical phenomena to the research community for now.
In the rarefaction region, the cells get elongated, leading to an aspect ratio larger than one.Because of the increase in the aspect ratio (i.e. the reduction of radial resolution), we find the peak of the velocity profile for the moving-mesh simulation becomes less sharp compared to the fixed mesh simulation.
However, since we have full control over the grid refinement, we can specify the maximum aspect ratio in the simulation.We conduct another moving-mesh sodtube simulation which sets the maximum aspect ratio to 1.5.When the elongated cell reaches such a threshold, it will split into two cells.To show the effect of such a refinement scheme on the sod-tube simulation, we compare the profiles for the moving-mesh simulation with or without maximum aspect ratio control in Fig. 6b.With the maximum aspect ratio control, the resolution in the region where the cell's aspect ratio gets to the threshold value increases.The peak of the velocity profile becomes sharper compared to the peak for the moving-mesh simulation without aspect ratio control.Overall, the implemented HLLC Riemann solver on the moving mesh is robust for simulating relativistic outflow.

B. Relativistic jet emerged from a black hole-torus system
The detection of the gravitational wave (GW) signal GW170817, coupled with the observations of its electromagnetic (EM) counterpart signifies the commencement of the multimessenger astronomy era [118].Research has demonstrated that the structure of the emerged relativistic outflows plays a crucial role in shaping the afterglow emission in GRB170817A [51,[119][120][121][122][123][124][125].This event provides an ideal candidate for utilizing the electromagnetic observations of the emerged outflow to infer the BNS merging physics.While the presented moving-mesh code is capable of simulating relativistic jets out of various progenitor systems, in the following, we will use a pseudomodel inspired by the outcome of compact binary merger simulations (see e.g., [28,[126][127][128][129][130]).We set up a black hole-torus system in the isotropic coordinate of Schwarzschild metric, with the mass of the central black hole M having been set to M = 3M ⊙ and the torus mass set to 0.2 M ⊙ .The radius of the inner edge of torus is 6 M , and the radius of its pressure maximum is set to 16 M [131].For the simulation domain, the radius of the inner boundary locates at 12.And we use Nt = 256 grids to cover the half spherical domain with cell's initial aspect ratio been set to 1.We adopt the reflecting boundary condition in the azimuthal direction.Outside the torus, the domain is filled with an ejecta cloud with a total mass of 0.02 M ⊙ .The cloud density structure follows: ρ 0 = (n−3)M ejecta /(4πr 3 0 ) is derived to give a total ejecta mass M ejecta .The pressure is P = Kρ Γ + P floor , K = 0.54, Γ = 4/3.We also add a density floor ρ floor = 6.2 × 10 −5 [g/cm 3 ]( i.e. 10 −22 in the code unit) and a pressure floor P floor = 10 −4 ρ floor to the initial profile to avoid numerical precision error.We set the density slope index n = 3.5 to represent the postmerger ejecta profile.Here, we ignore the ejecta profile velocity for simplicity.The reference radius r 0 is set to 6M ⊙ .A jet engine with a variable luminosity of L 0 = 2×10 51 exp(−t/t decay ) [erg/s] operates for 30 ms, in the polar region just above the black hole-torus plane.The engine decay timescale t decay has been set to 100 ms.This gives a total injected jet engine energy 5.2 × 10 49 [erg].We choose this low-energy jet engine injection to test the code's capability of launching a relativistic jet under constraint.In jet simulations, it becomes easier to successfully launch a relativistic jet given a higher energy injection (see e.g.[51,64]).The profile of the jet engine features a narrow nozzle with an opening angle of 0.1 rad.For the complete jet engine profile, we refer readers to the description in Appendix C as well as in [63,64].
Figure 8 shows the jet launching process during the first 30 ms.At the beginning of the simulation, the cloud flows into the black hole.In the polar direction, at a location centered around 130 [km], a small amount of relativistic gas with a terminal Lorentz factor 100 (i.e.jet engine) gets injected into the cloud.The injected gas has an initial boost velocity in the radial direction (see Appendix C).The addition of the relativistic gas slightly pushes the cloud gas in the polar direction, leading to a non-negative radial velocity (as can be seen from the radial velocity plot at t = 0.1 [ms]).The continuous injection of hot relativistic gas drives shocks and changes the temperature profile in the polar direction.By the time t = 5 [ms], a shocked cocoon develops and reveals a two-layer structure: a high-density layer which results from the forward shock, meanwhile the inner cocoon which heats up by the jet engine and reverse shock gets to a low-density regime [ 41,132,133].Inside the inner cocoon, the shocked gas accelerates to a high velocity with a maximum Lorentz factor around 8 at t = 5 [ms].The moving-mesh scheme dynamically allocates cells to resolve the shocked region.The interfaces of the double-layer structure can be seen in the contour plot for the cell's radial resolution: the first interface lies in the shock front between the cocoon and the unperturbed cloud, the second interface is between the cocoon's inner low-density hot relativistic core and its high-density colder part.At the bottom of the cocoon, the shock front hits the torus.At t = 10 [ms], the shock front starts to move beyond the torus and wrap around it.At the head of the cocoon, the loaded matter diverts part of the shocked gas sideways.Below this region, the inner core of the cocoon accelerates to a higher Lorentz factor of 13.Throughout the acceleration period, the maximum Lorentz factor of the jet reaches 20 (which happens at about t = 13 [ms]), smaller than the terminal Lorentz factor of the injected relativistic gas.This is largely due to the engine's relative low-energy budget (we refer readers to more energetic jet simulations in [51,64]).By the end of the jet engine injection t = 30 [ms], at the base of the grid domain, the frontier of the shocked cocoon has passed the torus region.A relative high-density buffer zone appears between the torus and the cocoon (see the density and temperature contour plots).The torus itself rotates stably during the jet launching process, as illustrated by the inner contour plots in Fig. 8.The head of the shocked cocoon expands beyond the initial grid domain boundary.More cells will be allocated in front of the boundary as the shock front propagates.The radius of the new boundary will make sure that the head of the shock front will stay below 0.8 of this new radius during the simulation.
Figure 9 illustrates the continued evolution of the relativistic cocoon that emerged from the black hole-torus system when the jet engine had been turned off.For computational efficiency, we cut the grid domain within a radius of 300 (around 440 [km]).We let the inner boundary move with a velocity of a fraction of the fluid's local maximum velocity.When the inner and outer boundary expands outward together, the simulation gets into a (weak) scaling region where the simulation time step ∆t increases with the absolute time itself.In this way, the moving-mesh code in spherical coordinates can simulate long-term evolution of relativistic jets over multiple orders of magnitude of time.When the jet engine turns off, it stops accelerating nearby gas.The inner cocoon turns into a shrinking relativistic bubble.The relativistic bubble keeps pushing against the mass-loaded head, exhausting its internal energy and kinetic energy.By the time t = 50 [ms], the top of the bubble starts to decelerate dramatically.The collision with the mass-loaded head converts part of its kinetic energy into thermal energy.The collision drives a wave passing through the relativistic bubble, increasing the temperature of inner cocoon all the way to its bottom (see the simulation video).Eventually, the relativistic bubble turns into a relativistic thin shell (see the velocity contour at t = 300 [ms]).The relativistic shell features a relativistic core with a mildly relativistic sheath, similar to previous special relativistic jet simulations (see, e.g., [51,64,126,134]).The outer layer of the cocoon goes through adiabatic expansion.The density within this layer keeps decreasing while the temperature structure roughly remains the same (see the contour plots for the temperaturelike variable from t = 50 [ms] to t = 1 [s]).The interface between the inner and outer cocoon features a high-density pillar.The relativistic shell keeps sweeping through the medium while depleting its kinetic energy.By the time t = 100 s, the relativistic thin shell is replaced by a mass-loaded slowmoving core.We see a morphological change in the outer shell structure.Finally, by the end of our simulation t = 1000 s, the outflow velocity becomes completely Newtonian.Now we have seen the complete life cycle of a relativistic jet from its birth at a black hole scale to the distance of its dissipation.In the following, we would like to discuss two dynamical features for this specific simulation.The first feature of interest focuses on the base of cocoon.In Fig. 9, we use white circle to indicate this region of interest.We find it appears after the shock front of the cocoon passed over the torus.It originates from the buffer zone or shock zone between the original torus and the remaining cocoon.It propogates subrelativistically and maintains its hump shape before t = 1 [s].Later on, the shock front accumulates enough matter and slows down to Newtonian velocities.When this happens, via hydrodynamical interaction, the morphology of the region changes and the hump shape disappears, leaving behind a broken filament as shown at t = 100 s and t = 1000 s.The second feature of interest relates to the density pillar at the interface between the inner and outer cocoon.We mark it with a red square in the figure.Its formation, to some extent, comes from the shutdown of the central jet engine during the jet launching period.At the beginning, when the central jet engine inflates a cocoon, it drives mass and energy into the cocoon outer layer while creating a low-density hot inner funnel to generate relativistic outflow.When the central engine shuts down, the inner cocoon quickly gets cold and stops pushing the outer layer (see snapshots at t = 50 [ms] and t = 300 [ms]).Then the adiabatic expansion of the outer layer further separates this interface from the shock front as shown at t = 1 [s].The interface pillar also has positive radial velocity and moves with the outer shell.However, the part, connecting to the outer shell, moves faster.To a point, the pillar detaches itself from the outer shell and falls back to the inner region.This is what happens from t = 100 s to t = 1000 s.Because of the long-term simulation of the relativistic jet, we are able to capture such detailed hydrodynamics evolution, which may provide insights for the study of morphologies of astronomical jets.
Throughout the simulation, the maximum grid resolution in the radial direction remains below 60-a value we set initially.We see that the moving-mesh scheme, combined with the dynamical grid refinement and derefinement can capture the detailed dynamical features for the relativistic jet simulation over many orders of magni-FIG.9.The full-time-domain evolution of the relativistic jet emerged from a black hole-torus system.The inner boundary moves outward at a fraction of the local maximum velocity.The simulation video is available from Youtube at [117] with high-definition video quality available.
tude of space and time.Also the adjusted moving-mesh scheme makes the simulation of relativistic jets computationally efficient.The presented simulation has been performed on a single high-performance computing node with 32 Intel Xeon Gold 6148 CPUs.The whole simulation consumes around 6400 core hours.

VII. CONCLUSION
This paper presents an advancement in computational astrophysics: developing and implementing a general relativistic moving-mesh hydrodynamic code featuring an advanced Riemann solver in curvilinear coordinates.We showcase the details of integrating a general relativistic framework into the hydrodynamic simulation code JET, achieved through applying the reference metric method.
With its ability to elegantly handle the intricate spacetime geometries inherent in general relativity, the tetrad formulation is an ideal choice to address the HLLC Riemann problem under strong gravity.The achievement of our work is the successful adaptation of the tetrad formulation to incorporate the HLLC Riemann solver into the general relativistic moving-mesh code.We have conducted a series of numerical simulations to validate and demonstrate the efficacy of our newly developed code.These simulations encompass both fixed-mesh and moving-mesh scenarios, allowing us to test the code's performance under various conditions.The results from these simulations are particularly noteworthy in demonstrating the code's robustness and reliability in simulating fluid flows under the influence of strong gravitational fields.
Compared to the fixed-mesh approach, a moving-mesh scheme can increase the time step for fluid regions with high velocity since it removes the limitation imposed by the bulk velocity (see, e.g., [61]).The moving-mesh approach makes the long-term simulation of relativistic jets computationally feasible (see, e.g., [51,63,64]).By extending the JET code's capability of handling relativistic jets from an astronomical scale to the scale of a black hole, we have opened new avenues for the full-timedomain simulation of relativistic jets, from their genesis to dissipation.To demonstrate this possibility, we design a representative prototype model which features a torus around a central black hole.A jet is manually launched in the polar direction, near the black hole-torus system.For the underlining jet launching mechanism, we refer readers to Blandford-Payne [4] or Blandford-Znajek [3] and related dynamo processes (see e.g.[13,14]).In this work, we prescribe an engine profile to imitate the jet launching process.This setup allows us to explore the dynamics of the jet's journey from its origin near the black holetorus system to its final Newtonian phase.We found multiple new hydrodynamical features from this end-toend simulation.For the first time, we have been able to simulate the complete life cycle of a relativistic jet, providing insights into the detailed structures of the cocoon and emerged jet over the whole journey.Furthermore, these full-time-domain jet simulations enable the joint investigation of various electromagnetic phenomena associated with relativistic outflows.For the case of BNS mergers, the observational phenomena include the kilo-nova emission from the remnant ejecta (see e.g.[24,25,135]), the GRB prompt and afterglow emission (see e.g.[23,136,137]), and other related processes.By combining our simulations with GRMHD models of jet-launching processes, we will be able to extend the evolution of outflows to distances relevant to long-term electromagnetic radiation observations.This integrative approach aligns perfectly with the era of multimessenger astronomy, allowing for an unprecedented understanding of the underlying physics in jet-launching systems.
While this paper sets the foundational steps in this direction, the complete realization of these ambitious goals remains a pursuit for future research.The potential for further advancements and discoveries in the field is vast, and our work may catalyze the next generation of astrophysical jet simulations, potentially revolutionizing our understanding of relativistic jets and their associated physics.

FIG. 1 .
FIG. 1. Top: radial rest-mass density profile of the Bondi flow for the medium resolution Nt = 256.The mass accretion rate is Ṁbondi = 10 −4 M .The blue dotted line indicates the initial analytical solution at T = 0 M , while the yellow dashed line denotes the numerical solution at T = 50 M .The green dashdotted line shows their difference.The black dotted vertical line indicates the radial location of the event horizon of the black hole at r = 0.78 M in the maximally sliced trumpet coordinate.Middle: the radial fluid velocity −W v r profile.W is the Lorentz factor, and v r is the fluid velocity with respect to the Eulerian observer.Bottom: the L1-norm of error for the rest-mass density as a function of the azimuthal grid number Nt.The dash line indicates second-order convergence.

FIG. 2 .
FIG. 2. Normalized variation for the central maximum density as a function of dynamical time for the TOV star at two resolutions (Nt = 64 and Nt = 128).Top: the time evolution for the stationary case.Bottom: the time evolution for the pressure-depleted star whose equilibrium pressure has been reduced by ten percent globally.

a
M is the mass of the central Schwarzschild black hole.R in is the location of the inner edge of the torus.Rmax is the pressure maximum location of the torus in Schwarzschild coordinates.ρmax is the peak density in the torus, Φ Rmax is the angular velocity at Rmax.(calculated from the simulation with Nt = 512 resolution).l = u t u ϕ is the constant specific angular momentum.κ and Γ define the EOS of the initial torus P = κρ Γ .Unless specified, the presented variable value follows the G = c = M ⊙ = 1 unit convention.

FIG. 3 .
FIG. 3. Contour plot of a stationary torus around a static black hole.The simulation adopts an axisymmetric spherical domain.The azimuthal angle extends from π/3 to 2π/3.The radial grid extends from rmin = 4 to rmax = 20.The top panel shows the initial profile of the logarithmic density.The middle panel represents the density profile at the end of the simulation t = 2000.The bottom panel represents the contour plot of the radial velocity W v r .

FIG. 4 .
FIG.4.Time snapshots for the density contour of a hot low-density bubble embedded inside a cold Bondi flow.The left panel shows the simulation results with the HLLE Riemann solver while the right panel presents the comparative results from the simulation with the HLLC Riemann solver.The simulation video for the case with HLLC Riemann solver is available from Youtube at[115].

FIG. 5 .
FIG. 5. Profiles for the spherical shock tube test at t = 0.3 in fixed-mesh simulations with the HLLC (solid line) and HLLE (dotted line) Riemann solver.
FIG. 6.(a) Shock tube test in spherical coordinate for the fixed-mesh (black-dot) and moving-mesh (MM) (red-square) simulations.The profiles are presented at t = 0.3.Initial discontinuity locates at x = 0.25.The inner (outer) state of the shock tube is ρ = 1, P = 1, v = 0 (ρ = 0.1, P = 0.1, v = 0).The EOS follows the ideal gas law with Γ = 4/3.(b) Shock tube test in spherical coordinate for the moving-mesh simulations with different control schemes for the cell's aspect ratio.The black-dot line represents the simulation where the cell's maximum aspect ratio has not been set.The red-square line shows the simulation where the maximum aspect ratio of the cell has been set to 1.5.

FIG. 7 .
FIG. 7. Profiles for the spherical shock tube test at t = 0.3 in the moving-mesh simulations with different grid resolution.The purple diamond shows results from the simulation with resolution Nt = 2048 while the red square and blue circle represent results from the Nt = 128 simulation and the Nt = 32 simulation, respectively.

FIG. 8 .
FIG.8.Jet launching process from a black hole-torus system, visualized at four-time snapshots.From left to right, the contour plots represent the logarithmic density, temperature-like variable Θ, normalized radial velocity W v r /c, and the cell's radial resolution (i.e.inverse aspect ratio).The inner contour plots zoom in on the central region of the domain.The inner boundary is stationary before t = 30[ms].The simulation video is available from Youtube at[117] with high-definition video quality available.

TABLE I .
Parameter values for the initial profile of TOV star.

TABLE II .
Parameters for the stationary torus around a black hole a