Shocks in radiatively driven time dependent, relativistic jets around black holes

We study time-dependent relativistic jets under the influence of radiation field of the accretion disk. The accretion disk consists of an inner compact corona and an outer sub-Keplerian disk. The thermodynamics of the fluid is governed by a relativistic equation of state (EoS) for multispecies fluid which enables to study the effect of composition on jet-dynamics. Jets originate from the vicinity of the central black hole where the effect of gravity is significant and traverses large distances where only special relativistic treatment is sufficient. So we have modified the flat metric to include the effect of gravity. In this modified relativistic framework we have developed a new total variation diminishing (TVD) routine along with multispecies EoS for the purpose. We show that the acceleration of jets crucially depends on flow composition. All the results presented are transonic in nature, starting from very low injection velocities, the jets can achieve high Lorentz factors. For sub-Eddington luminosities, lepton dominated jets can be accelerated to Lorentz factors>50. The change in radiation field due to variation in the accretion disk dynamics will be propagated to the jet in a finite amount of time. Hence any change in radiation field due to a change in disk configuration will affect the lower part of the jet before it affects the outer part. This can drive shock transition in the jet flow. Depending upon the disk oscillation frequency, amplitude and jet parameters these shocks can collide with each other and may trigger shock cascades.


INTRODUCTION
Astrophysical jets are collimated outflows of plasma associated with various Galatic and extragalactic sources, such as neutron stars, young stellar objects (YSOs), gamma-ray bursts (GRBs), microquasars, and active galactic nuclei (AGNs). These jets span over a wide range of energy and length scales. For example, the jets of YSO travel a distance of several parsecs while AGN jets are observed to remain collimated for hundreds of kiloparsec range. In terms of jet kinetic luminosity, YSO jets have the luminosity of the order 10 32 erg s −1 while the GRB jets can reach up to 10 52 erg s −1 . Since the advent of radio telescopes in the era of modern astronomy, the understanding of these ubiquitous objects has improved significantly, but there is no consensus about the formation, composition, and collimation of these objects (Harris & Krawczynski 2006). Black holes (BH), which sit in the heart of microquasars and AGNs, do not have hard surface and therefore the origin of the jet should be the accreting matter itself. Interestingly, simultaneous radio and X-ray observations show a strong correlation between different spectral states of accretion disks and jets states in microquasars (Gallo et al. 2003;Rushton et al. 2010;Fender et al. 2010), which indicates that the jet originates from the accretion disk. Moreover, the jet is launched from the inner part (less than 100 Schwarzschild radius or r g ) of the accretion disk (Junor et al. 1999;Doeleman et al. 2012). While the illusive corona i.e., the source of hard power law radiation, has been identified as the inner hot part of the accretion disk in many different models (Shapiro et al. 1976;Chakrabarti & Titarchuk 1995;Gierlinski et al. 1997). Therefore, the compact corona may be responsible for the jet activity. Additionally, the jet which originates from a region very close to the black hole (BH) horizon, will travel through the radiation field of the disk, and should certainly be influenced by it.
In recent years, numerical simulations of relativistic jets have played a crucial role in explaining the non-linear dynamics of the system. There are many studies to investigate the propagation of supersonic jets (Walg et al. 2014;Martí 2019;Seo et al. 2021a). These simulations present a global picture of the jet with a forward bow shock, backflow of the fluid, and multiple shocks within the jet beam. Moreover, the role of magnetic field in the jet formation and collimation has also been studied (Sheikhnezami et al. 2012;Komissarov & Porth 2021;Fuentes et al. 2021). But the numerical simulations of radiatively driven outflows and jets are very limited, in spite of the fact that the interaction between the radiation field and plasma is not a new subject. The equations of radiation hydrodynamics have been developed by many authors (Mihalas & Mihalas 1984;Kato et al. 1998;Park 2006;Takahashi 2007) and these equations have been used extensively in steady-state investigations of the radiatively driven outflows around compact objects. Icke (1989) used relativistic equation with near disk approximation of radiation field above a Shakura-Sunyaev disk (Shakura & Sunyaev 1973) to obtain an upper limit of jet terminal speed. Fukue (1999) showed that under the influence of radiation arising from disk corona, the jet can achieve significant collimation. So, the radiation not only accelerates the jet, it also helps in collimation. Later, Fukue et al. (2001) showed that electron-proton jets can reach relativistic terminal speeds by considering a hybrid disk with inner advection dominated flow (ADAF) (Narayan et al. 1997) and an outer Keplerian disk. The disk model considered by Chakrabarti & Titarchuk (1995) with a mixture of matter with sub-Keplerian and Keplerian angular momentum, showed that sub-Keplerian disk can go through a shock transition and create a hot and compact post-shock disk, which can act as an illusive corona. It may be noted however, inner, puffed hot disk acting like a corona has been suggested by other authors too (Shapiro et al. 1976;Gierlinski et al. 1997). Numerical simulations of advective disks show that the extra thermal gradient term in hot post-shock region automatically generates the bipolar outflows (Das et al. 2014;Lee et al. 2016) which can be accelerated by the radiation field of the disk. In the relativistic radiation hydrodynamic regime, investigations by Chattopadhyay (2005) and Vyas & Chattopadhyay (2017, 2019 clearly highlight that the radiation field of the accretion disk plays a significant role in the acceleration of the jets. In addition to these steady-state investigations there are limited numerical simulations in the non-relativistic regime (Chattopadhyay & Chakrabarti 2002a;Chattopadhyay et al. 2012;Raychaudhuri et al. 2021;Joshi et al. 2022). To study the effect of radiation on winds there are some simulations of line driven winds (Proga et al. 2000;Yang et al. 2018), but the line driven force is only effective when the wind is cold, so the line driven force is not applicable in hot plasma which is being considered in this paper.
In this paper, we have made an attempt to understand the dynamics of the jet in the time-dependent radiation field of the accretion disk. The steady-state investigations by Vyas & Chattopadhyay (2017, 2019 and Ferrari et al. (1985) have shown that the radiation field can produce internal shocks in jets but how do these shocks evolve under a time-dependent radiation field is not known yet. So, the central theme of this paper is to investigate the formation and evolution of these shocks in relativistic jets due to its interaction with the radiation field of the accretion disk. We study jets which are launched with subsonic speeds close to the central BH. Since the jets are launched close to the BH so one has to consider the effect of gravity at the jet base. However, due to the huge length scales jets traverse, in most parts gravity is negligible and special relativity suffices. Previous studies like Ferrari et al. (1985); Vyas et al. (2015) used the Newtonian or pseudo-Newtonian potential as an additive source term for gravity in the special relativistic Euler equation but this destroys the conservative nature of equations of motion which is essential in the Eulerian upwind numerical schemes we are using. Hence, we include the gravity by modifying the special relativistic metric itself and thereby retaining the conservative form of the equations of motion. The presence of gravity makes the jet transonic, while we avoid general relativistic details. The temperature of the jet varies by a large value so we use a relativistic equation of state (EoS) in our analysis given by Chattopadhyay & Ryu (2009). This EoS (abbreviated as CR) handles the variation of the adiabatic index and also allows us to study the effect of plasma composition on the jet solution.
The internal shocks in the jet are identified as bright knots in the jets. These shocks are the sites for particle acceleration in the jets, resulting in non-thermal emissions. The jet models with shock are very popular and have been used to explain the outbursts in microquasars (Kaiser et al. 2000), AGNs (Walg et al. 2014), as well as GRB jets (Rees & Meszaros 1994). Shock very close to the jet base around a stellar mass BH was considered to explain the high-energy power law emission (Laurent et al. 2011). Numerical simulations of jets show that the shocks form in the jet beam whenever a fast moving supersonic plasma catches up with a slower moving flow ahead. It is well known that the radiation field accelerates by pushing the jet material. However, due to relativistic effects, radiation would also decelerate the jet close to the disk via a process called radiation drag (Icke 1989;Fukue 1999;Chattopadhyay & Chakrabarti 2002b;Chattopadhyay et al. 2004;Chattopadhyay 2005). Therefore if the disk radiation is timedependent, can it accelerate the jet at lower altitudes to trigger shocks by crashing the lower altitude faster blobs onto the slower portion of the jet above it? As these jet shocks are time-dependent, how does the shock strength evolve in time? If the disk is oscillating, can the resulting time-dependent radiation field produce multiple shocks? Further, would these shocks collide with each other and generate shock cascades? How would the jet solutions depend on the baryon fraction of the jet? In this paper, we address these questions.
In section 2 we present the assumptions and governing equations used in this paper. In section 3 we outline the methodology to obtain the solutions. Section 4 is devoted to explain and discuss the key findings of this paper. In the end, in section 5 we conclude our work.

ASSUMPTIONS AND GOVERNING EQUATIONS
We study the jets driven by the radiation field of an accretion disk around a non-rotating black hole. We do not include the accretion disk directly in our study, the disk plays a supportive role by supplying the radiation field. All the disk parameters used in our study are required only for the calculation of the radiative moments. We have not explored the direct connection between disk and jet generation, instead, we assume the jet injected with a finite subsonic initial velocity and some initial temperature. In this study we restrict our calculations for a relativistic non-rotating (u φ = 0), axis-symmetric (∂/∂φ = 0), on-axis (θ = 0) jet. We assume a conical geometry of the jet with a narrow opening angle. The energy-momentum tensor for jet material (T αβ M ) and radiation field (T αβ R ) are given as where u α represent the components of four velocity, l α s are the direction cosines, I ν is the specific intensity of the radiation field. dΩ is the solid angle subtended by the field point on the jet axis to the source point on the accretion disk. ρ is the rest mass density of fluid and h is the specific enthalpy of the fluid which is related to the internal energy density e and pressure p as The space-time metric is given by ds 2 = −g tt dt 2 + dr 2 + r 2 dθ 2 + r 2 sin 2 θdφ 2 Where g tt is given as; We have used the geometric unit system where 2G = M B = c = 1, M B is the mass of the BH, G is the gravitational constant, and c is the speed of light. The information of gravity is supplied through g tt , somewhat in the spirit of weak-field approximation, although the Φ is given by the Paczyńsky-Wiita potential (Paczyńsky & Wiita 1980), instead of the Newtonian one. This way, we retain the conservative nature of equations of motion. In this unit system, the Schwarzschild radius (r g = 2GM B /c 2 ) is the unit of length and t g = 2GM B /c 3 is the unit of time.
The equations of motion for relativistic radiation hydrodynamics have been derived before (Park 2006;Takahashi 2007), here we present only a brief description to obtain the equations of motion. The continuity equation which represents the conservation of mass flux is given as (ρu α ) ;α = 0. The equations of motion are essentially the conservation of energy momentum tensor or, (T αβ M ) ;β = −(T αβ R ) ;β = G α . For a conical jet along the axis of symmetry implies that u r is the only significant component of the four-velocity and conservations of mass flux and energy-momentum tensor can be written as a set of three conservation laws: Here v is the three velocity of the jet. The four velocity (u r ) can be written in terms of v as u r = γv, where γ is the Lorentz factor given as γ 2 = −u t u t = 1/(1 − v 2 ). The conserved quantities D, M, and E are the mass density, momentum density, and total energy density in the laboratory frame given as (also see Ryu et al. 2006;Seo et al. 2021b) G r and G t are the components of the radiation force, which under the present set of assumptions are explicitly given as (Mihalas & Mihalas 1984;Kato et al. 1998;Park 2006) where Here ρ e is the total lepton density and F = σ T F rd /(m e c), E = σ T E rd /(m e c), P = σ T P rd /(m e c), also, E rd , F rd , & P rd are the zeroth moment or energy density, first moment or the flux, and second moment or the pressure of the radiation field, respectively.
To solve the equations of motion (5-7) we need an additional closure relation which relates the thermodynamic variables h, p, and ρ. This closure relation is known as the equation of state (EoS). In this study, we consider an EoS for relativistic multispecies fluid proposed by Chattopadhyay & Ryu (2009) (abbreviated as CR EoS). CR EoS is a very close fit to the exact one given by Chandrasekhar (1939). The CR has already been used for different astrophysical problems (Cielo et al. 2014;Singh & Chattopadhyay 2019;Sarkar et al. 2020;Joshi et al. 2021Joshi et al. , 2022.
The EoS is given as where, f = 1 + (2 − ξ)Θ 9Θ + 6/τ 6Θ + 8/τ + ξΘ 9Θ + 6/ητ 6Θ + 8/ητ In equations (13, 14) ρ is the mass density of fluid given as ρ = Σ i n i m i = n e − m e (2 − ξ + ξ/η), where ξ = n p /n e − , η = m e /m p and n e − , n p , m e and m p are the electron number density, the proton number density, the electron rest mass, and proton rest mass. Θ = p/ρ is a measure of temperature and τ = 2 − ξ + ξ/η. The expression for sound speed is given as Here, N is the polytropic index given as The polytropic index is a function of temperature. Hence, we do not need to supply it as a free parameter. From equation (16) one can see that for higher values of temperature Θ >> 1, N → 3, and N → 3/2 for Θ << 1. The adiabatic index is related with the polytropic index as

Steady state equations of motion
In the steady-state, we impose ∂/∂t ≡ 0 on equations of motion (5-7). Integrating equation (5) then leads to the mass-outflow rate.Ṁ = ργvA (18) where A(∝ r 2 ) is the cross-section of the jet. The energy balance equation or 1 st law of thermodynamics is given by Integrating the above and combining with equation 18 gives the adiabatic relation (also, see Kumar et al. 2013;Joshi et al. 2022 where Since interaction of radiation with the jet in the Thompson scattering regime is isentropic, so the entropy-outflow rate is a constant along a streamline in the steady-state, but discontinuously jumps up in case there is a shock transition in the jet. The relativistic Euler equation can be simplified to the following form, where, the last term in the numerator γR r ρh is the radiation momentum deposition term, which can be simplified to the form We can integrate equations (22) and (19) to obtain another constant of motion for steady flow known as the generalized relativistic Bernoulli parameter (B e ).
We consider an advective type accretion disk (Fukue 1987;Chakrabarti 1989;Lee et al. 2011;Kumar & Chattopadhyay 2017) for this study. Figure 1(a) shows a typical cartoon diagram for the disk and jet system. As mentioned earlier, the disk plays only a supportive role in our calculations by supplying the radiation field. Depending upon the accretion rate advective disks have moderate to low radiative efficiency (Sarkar et al. 2020). The disk has two components, an inner puffed up post shock disk called PSD and a sub-Keplerian disk (SKD). The sub-Keplerain part of the disk supplies the soft photons and the geometrically thick PSD also known as corona supplies the hard photons. This type of disk structure can mimic hard to hard-intermediate spectral states in microquasars. The PSD terminates at the horizon but we take the inner edge of the corona at 1.5 r g as the disk is expected to emit negligible radiation from a region within it. Numerical simulations of advective accretion flows suggest that the SKD is flatter than the PSD (Lee et al. 2016) and the ratio of height to the horizontal extent of corona can vary from 1.5 to 10, so we take a semi-vertical angle θ sk = 85 o for SKD and intercept of SKD on jet axis is taken to be d 0 = 0.4h s , where h s is the height of corona taken as h s = 2.5x s . The outer edge of SKD is taken at x o = 3500. The details of obtaining the radiative moments are presented in Chattopadhyay & Chakrabarti (2000); Chattopadhyay (2005); Vyas et al. (2015), here we present only a brief description.
In Fig.1(b) we have shown a cross sectional view of disk and jet sytem. M and M are the source points on PSD and SKD, respectively. T M and T M are the local normals on disk surface. We calculate the radiative moments on the field point N . The inner edge of SKD is at x s but due to the shadow effect (Chattopadhyay 2005), the observer at N sees the SKD starting from x i given as We need the intensities of various disk elements to obtain the radiative moments of respective disk components. In SKD, we assume the synchrotron emission (Shapiro & Teukolsky 1983) to be the dominant emission mechanism.
where B sk , Θ sk , n sk , x represent magnetic field, local dimensionless temperature, electron number density, and horizontal distance from central object respectively. We assume a stochastic magnetic field in SKD with a constant magnetic to gas pressure ratio p mag = βp gas . The general definitions of radiative moments are and at each field point (e. g. N) these quantities have to be computed from both the sources SKD and PSD. All the components of the radiation moments should ideally be numerically computed in equations C7. The details of computing these numerically are given in Vyas et al. (2015) and also in Appendix C of this paper. However, numerically computing the radiative moments in the simulation code at every cell and at every time step would slow down the code significantly. So we fitted the numerically computed radiative moments (equations C8-C14) with approximate analytical functions (equations D15-D17 for SKD and D20-D22 for PSD in Appendix D) which depend on disk parameters like x s , h s & x in . The sub-Keplerain accretion rate controls the position of x s (equation A1) and luminosities (equation B5-B6) of the disk components. The comparison between the moments obtained from numerical integration (solid, dashed, dashed-dotted lines) and the analytical fitted functions (red, green and cyan open circles) for x s = 15 is shown in Fig.(2). The fitting is good, and the advantage is, any time variation ofṁ sk will produce a time dependent radiation field. But the moments are given as analytical functions and need not be numerically computed at every cell and at every time step. It may be noted that, E, F and P obtained from equation D15-D17, and D20-D22 for a givenṁ sk are used in equation 12. The steady-state solutions are obtained by simultaneously integrating equations (22) and (19). The jet at its base is hot and slow and therefore subsonic. The jet is powered by the thermal gradient term and the radiative term against gravity and becomes supersonic at some value of r. In other words, jets are transonic in nature and at some point r c called sonic point, the flow velocity v c crosses the local sound speed a c . At r c , dv/dr → 0/0 gives the sonic point conditions From equation (29) we obtain the temperature (Θ c ) at the sonic point and dv/dr is obtained by using L Hôpital s rule. We integrate the equations (22) and (19) inwards and outwards starting from sonic points to obtain the solutions. The steady-state solutions provides the injection values for the simulation code.

Time dependent solutions and numerical setup
We perform the time-dependent calculation using a relativistic total variation diminishing (TVD) routine. TVD scheme introduced by Harten (1983) is a second-order, Eulerian, finite difference scheme, which was initially proposed to solve the set of hyperbolic, non-relativistic hydrodynamic equations. However, many relativistic numerical codes have been built based on the same philosophy (see Martí & Müller 2003. A relativistic TVD code to incorporate CR EoS has been built before (Ryu et al. 2006;Chattopadhya et al. 2013). The detailed description of the code is presented in (Ryu et al. 2006), here we present only brief details.
The set of equations (5)-(7) can be written in the form with the state and flux vectors given as and the source term S as The set of equations 30 is solved in two steps. In the first step, the state vectors q at the cell center i at a given time step n are updated using the modified fluxesf i+1/2 calculated at cell boundaries. where, Here, α k,i+1/2 = L n k,i+1/2 · q n i+1 − q n i Explicit forms of eigenvalues λ k,i , eigenvectors L & R, flux limiters g k,i , γ k,i+1/2 , Q k are given in Ryu et al. (2006). In the second step, the radiation and curvature terms are added as source term in the updated values of state vectors. The next step consists of recovering the updated primitive variables (ρ, v, p) from the updated conserved variables D, M, E which is done by inverting the equations 8-10. These equations along with EoS can be expressed as a polynomial of γ, given as where g 1 = γM − E γ 2 − 1. We solve equation 36 using Newton-Raphson method to obtain γ, and once γ is known the primitive variables can be obtained as Since v, γ, E and M are known, then using equations 9 and 10, the pressure p is obtained as the root of a cubic equation a 1 p 3 + a 2 p 2 + a 3 p + a 4 = 0, where, a 1 = 27η, a 2 = 9 4 We employ a continuous inflow boundary condition at the injection cell and outflow boundary conditions at the outer boundary using ghost cells. In this time-dependent study, we want first to obtain the steady state solutions as is predicted by theoretical analytical investigations. And then we want to study the influence of accretion disk oscillation on the jet dynamics, while keeping the injection parameters same as the steady state values.
3.2.1. Time dependent disk as observed by an inertial observer on the axis In this paper, the time dependence of the jet is not imposed through the injection parameters, but by making the radiation field from the accretion disk time-dependent. In microquasars, quasi-periodic oscillations (QPO) are observed in hard X-rays associated with the oscillation of the inner part of the accretion disk (Nandi et al. 2012). If there is a fluctuation of the accretion rateṁ sk at the outer edge of the disk, the fluctuation would travel inwards. The time variability of the accretion rateṁ sk is given in equation A2, which induces a sinusoidal variation in the size of PSD or x s , is of the form, where the oscillation is about a mean x s0 (distance OG fo Fig. 3) with an amplitude of A (NG=GN ), as shown in Fig.(3) and frequency f s . The information about the disk-fluctuation propagates at the speed of light to the observer P, i. .e, it takes a finite amount of time to reach the jet axis. Hence the different points on the jet axis see a different disk configuration (see Joshi et al. 2022). In other words, the radiation emitted from x s at t reaches P after a time t + δt (δt = l /c, see Fig. 3), and in the same time-interval the PSD boundary will move from x s to x s . Therefore, The position x s is updated after each time interval dt which is determined by the TVD code. So we can write Where s = 0.5(v s, t + v s, t+dt ) is the average velocity in the time interval between t to t + dt. The instantaneous velocities can be obtained from equation 40, v s = A2πf s cos(2πf s t) From equations (41) and (42) we can write We solve equation (44) at each r and at every time step to obtain the x s as seen by the observer there, and the corresponding radiative moments are calculated using the analytical estimates of the radiative moments (equations D15-D22).

Code verification
The relativistic Riemann problem is a great tool to test the time-dependent hydrodynamic code, and the special relativistic TVD code with CR EoS was indeed tested against the analytical relativistic Riemann problem (see Joshi et al. 2021). To check the correctness and accuracy of the code in the relativistic radiation hydrodynamic regime in presence of gravity, presently we compare the numerical result with the steady-state semi-analytical radiatively driven relativistic jet solution. In Fig. 4 we have shown the solutions obtained from the TVD code marked by open red The injection parameters with which the jet is launched are chosen from the analytical solution and are given by v in = 0.065, Θ in = 0.037, r in = 5.0. The sonic point of the steady-state solution corresponding to these parameters is at r c = 14.0. We divide the computational domain of total length 1000r g into uniform 6000 cells. It may be noted that the simulation of outflows requires higher resolution compared to accretion disk simulations, because of the presence of sharp gradients at outflow injection boundaries. The comparison shows that the code regenerates the solution very accurately. The adiabatic index of the jet at the base is moderately relativistic, but becomes non-relativistic at large r, which suggests that the use of CR EoS to describe the thermodynamics of the jet is appropriate.

Effect of disk luminosity on jet solution
In Fig. 5a, we compared solutions of electron-proton (ξ = 1) jets corresponding to different disk luminosities, all starting with the same injection parameters v in = 0.018, Θ in = 0.063, r in = 2.5. These injection parameters are taken from a thermally driven steady-state solution characterized by the sonic point location r c = 15.0. Each curve represents thermally driven jet (solid, black curve), jets driven by disk luminosities l = 0.53 (solid, red), l = 1.0 (dash-dot, blue) and l = 2.4 (dashed, green), all in units of Eddington luminosity. The terminal speed of the thermally driven jet is v T = 0.13c. It is expected that higher disk luminosity would accelerate the jets to higher v T . And indeed Fig.  5b shows, v T increases with l, to the extent that an electron-proton jet for sub-Eddington disk luminosity (l < 1) may reach mildly relativistic speeds (v T ∼ 0.4 − 0.5c), while for the super-Eddington luminosities the jet can achieve relativistic terminal velocities. Interestingly, higher l suppresses the jet velocity v in the funnel like region above the PSD. From Fig. 2, it is clear that the radiative flux F ps < 0 in the funnel like region between the inner walls of the PSD. So with higher l, the F ps will actually suppress the jet speed within the funnel region, although at higher altitude would push the jet. Therefore v T increases with l.

Effect of jet composition
The interaction between the jet plasma and radiation is dominated by Thompson scattering, so the amount of momentum transferred to the jet will be higher for the lepton dominated jet which can be clearly seen from the equation (12). To study the effect of composition we fix injection and disk parameters and the solutions are obtained for various values of ξ. In Fig. (6 a) we have plotted the velocity profiles of various jet solutions, with the injection parameters v in = 0.11, Θ in = 0.22, r in = 2.0. The disk luminosity is l = 0.82. The radiation term R r for solutions with higher fraction of protons (higher ξ) will be weaker resulting in lower terminal speeds. But the radiation term inside the funnel is negative so the lower magnitude of R r means the low value of deceleration, so the proton dominated flows    These results clearly suggest that the interaction of disk radiation with the jet plasma can play a crucial role in the acceleration of the jet. In contrast to the electron-proton jets, the lepton dominated jets can reach ultra-relativistic speeds even for the sub-Eddington luminosities.

Jet solutions under the time dependent radiation field
To study the effect of disk oscillation on jet dynamics, we first obtain a steady-state solution corresponding to injection parameters of steady-state analytical solution. The steady-state time is taken to be t = 30000 t g . In Figs. 7 & 8 the injection parameters are v in = 0.076, Θ in = 0.064, r in = 3.0. The disk luminosity for the steady-state configuration is l = 0.597 and the outer edge of the PSD is at x s = 12.0 (obtained from equation A1). The velocity (v), Mach number (M ) and net radiation force term (R r ) for the steady-state solutions are plotted (blue, solid curve) in Figs. 7 (a 1 ), (b 1 ), (c 1 ), respectively. The sonic point for steady-state solution is at r c = 13.0. Once the steady-state solution is obtained we induce a variability inṁ sk (equation A2) such that it imposes an oscillation on x s with the amplitude A = 7.0, about the equilibrium value x s0 = 12 and the time period being T = 2000t g . In this study we have taken the frequency and ampltiude of oscillation as free parameters. Asṁ sk decreases, x s moves outward and the net disk luminosity decreases resulting in a weaker radiation field (compare blue and red curves in Figs. 7c 1 -c 5 ). The magnitude of net radiative force term |R r | for x s > x s0 is smaller than the steady-state values (blue) and is larger for x s < x s0 . So in the funnel region, v is higher when x s > x s0 but eventually produces a jet whose v T is weaker and vice versa. Figures 7 a 2 , b 2 and c 2 , show the solution for the configuration when the accretion disk is brightest and the flow is suppressed to remain subsonic within the funnel. However, the flow velocity is boosted in the range 20 < r < 40. The accelerated jet material encounters the slowly moving outer domain (r > 40). As the supersonic fluid is resisted by relatively slower fluid ahead, it drives a shock in the jet, shown in Fig. 7(a 3 ). The formation of shock further reduces the downstream velocity of the fluid, as a result we can see the formation of multiple shocks in the flow namely S 1 , S 2 , S 3 as shown in Fig. 7(a 4 ). The shock S 1 being faster, catches up and collides with the shock in front of it (S 2 ), shown in Fig. 7(a 5 ). In Fig. 8 we have plotted the zoomed region of jet domain 300 < r < 800 to highlight the changes in mass density (ρ) distribution due to collision of the shocks. Panels 8(e 1 ) − (e 5 ) show two shocks S 1 and S 2 collide and are plotted in the epoch between Fig.7(a 4 ) − (a 5 ). Figures 8(e 6 ) − (e 10 ) clearly shows that the collision results in a cascade effect and induces additional shocks in the flow. Of course, if the oscillation is stopped, the jet eventually returns back to its steady-state. In Fig. 9 (a), (b) we plot the variation of positions of jet shocks S 1 and S 2 and their compression ratios R as a function of time t. Figure 9 (a) shows that the slope of r sh vs t graph is steeper for the shock S 1 , which clearly suggests that the speed of S 1 is faster than that of S 2 and ultimately catches up with S 2 at r sh ∼ 535. The compression ratio is the ratio of post-shock and pre-shock densities (R = ρ + /ρ − ) which is an indicator for the strength of the shock. Panel(b) shows that both the shocks are very strong but the variation of R with respect to time has a complicated profile for both shocks. The shock S 1 becomes stronger as it moves outwards, reaching a maximum value of R ∼ 4.4 and gradually starts to slow down. The compression ratio for second shock S 2 reduces initially but starts to increase to reach a maximum value of R ∼ 4.55 and then it becomes weaker.
A closer look at equation 23 indicates that the radiative effects are strong for the colder flows and relatively weaker for the hotter jets. Hence, we also study a hotter jet solution with injection parameters v in = 0.14, Θ in = 0.226, r in = 2.0, considered from the analytical jet solution for a steady disk configuration x s = 11.02 and l = 0.664. The velocity (v) and Mach number (M ) variation for the steady-state solution is plotted in Figs. 10(a 1 ) and (b 1 ) with the solid blue line. The amplitude of oscillation for x s is A = 6.0 and the time period of oscillation is T = 500t g . In Figs. 10a 1 -b 5 we have plotted jet solution for different values of x s asṁ sk is varied with time. So the variation ofṁ sk determines the l and x s at various t. One of the distinct feature of this hotter solution is that it harbors an internal shock very close to the jet base (r sh ∼ 10) as shown in Fig. 10(a 3 , b 3 ) and this shock is transported outwards as the jet solution evolves in time. Since this is a hot jet near the base, radiative force is less effective, and it would be thermally accelerated to supersonic speeds within the first few Schwarzschild radii and more importantly within the funnel. A supersonic jet has low temperature and radiative terms become significant. Inside the funnel R r < 0, so resistance to the supersonic part of the jet drives a shock. This type of shock was not present in the previously shown solution because the jet base was relatively colder and so the jet was strongly decelerated to remain subsonic throughout the funnel region just  r (e 10 ) x s =12.0, t =32999.53 Figure 8. Density distribution ρ of the jet at vaious time steps are plotted to the collision between shock S1 and S2 in panels (e1 − e5). Panels (e6 − e10) show the ρ after the collision. The PSD oscillation is same as in Fig. 7. above the PSD. It may be noted that such shocks in jets have been invoked to explain high energy emissions (Jamil et al. 2008;Laurent et al. 2011).

Effect of disk oscillation frequency on jet solutions
We have discussed earlier in section 3.2.1 that the information regarding the time variability of the disk (or portions of it) does not reach the entire jet simultaneously. The photons take some finite amount of time to reach different parts of the jet. Hence, the different points on the jet axis at a time t would see the outer edge of the PSD at different locations. If the oscillation frequency of the disk is very low then the emitted photons can reach to a larger distance of the jet axis before the outer edge of the PSD reaches a new position. In Fig.(11) we have compared the jet solutions corresponding to two different frequencies of the disk oscillation. The injection parameters are v in = 0.076, Θ in = 0.064, r in = 3.0 (same as Fig. 7). The time period of the disk oscillation is T = 2000 t g for solution plotted with solid-red line and T = 10 5 t g for the solution plotted with dashed-green line. Here the jets are compared at the same disk configuration (same l and x s ). The time at which the disk reaches same l and x s , will be different for the two frequencies, and has been shown on each panel. Assuming an M B = 10 M these time periods correspond to oscillation frequencies 5 Hz and 100 mHz, respectively. The microquasars show the quasi periodic oscillations (QPOs) at various frequencies.
For exmaple GRO J1655-40 (with central mass ∼ 7M ) itself shows the QPOs at various frequencies extending from 0.1 − 10Hz (Chakrabarti et al. 2008;Remillard et al. 1998). Hence, as a representative case for two different frequencies of disk oscillations we have taken the values 100 mHz and 5 Hz. These results clearly indicate that for a low frequency oscillations in disk the transient features (shocks) are not present in jet solutions. The jet velocity gets boosted up or slowed down depending on the disk configuration but the solutions remain smooth throughout the domain.

DISCUSSION AND CONCLUSION
In this paper, we studied the relativistic jets driven by the radiation field of advective type disks. The jet gains momentum supplied by the radiation field. We use relativistic radiation hydrodynamic equations along with a relativistic EoS to study the jet properties. We use the steady-state analytical solutions to verify our time-dependent code. The steady-state solutions also provide us the injection parameters at the jet base to be used in time-dependent cases. The accretion disk plays an auxiliary role in our study by supplying only radiation. It may be noted that we have considered sub as well as super-critical accretion rates to compute the radiation field above the disk. In the advective disk regime, super critical accretion rates may make the inner part of accretion disk optically slim (1 <optical depth < ∼ few), although higher infall velocity keeps the outer part of the disk optically thin. In such cases, radiative transfer equation need to be solved to compute the radiative intensity emitted by the disk, especially from the inner part. Since the disk plays an auxiliary role and is not part of the computational domain, we considered simplifying assumptions to compute the radiation field. We show the effect of radiative acceleration on the terminal speeds of the jets. The super-critical disks accelerate the electron-proton jets to terminal Lorentz factor γ T > ∼ 4. And for the jets driven by sub-critical disks also achieve reasonable relativistic speeds. The EoS allowed us to study the solutions corresponding to various plasma compositions. We found out that the electron-positron jets can be accelerated to  Fig.(7). The two jet solutions in each panel corresponds to same xs and l, and the time in which the disk configuration is similar for the two cases are mentioned in the panels . ultra-relativistic speeds (γ T > 50). The radiative acceleration is quite impressive, however, we have not considered the effect of back scattered photons onto the jet base, which may inhibit the jets. That might slightly lower such impressive terminal speeds we have obtained. In a nutshell, it may be summarized that the net radiation force (R r ) is negative inside the funnel region above the PSD and therefore opposes the jet expansion, but above the funnel region the net radiative force is positive and therefore the jet is accelerated. The time dependence of the jet is imposed by the time dependence of the radiation field. It may be recalled that Kumar et al. (2014) studied the radiatively and thermally driven outflows from advective accretion disks and they have shown that the radiative momentum deposition has a more significant effect on jet acceleration than conversion of thermal energy to the kinetic one. Moreover, they also found out that for a given r in (jet base), v in is higher for higher luminosity, but difference in sound speed or Θ in is imperceptible. Hence we have fixed the jet base for this study. It may be noted that, we have checked that by imposing a 20% fluctuation of v in with a frequency similar to that of the oscillating x s , it produces imperceptible change in the jet solution, therefore keeping the jet base values fixed was reasonable. From its steady-state value as x s decreases, l increases, the funnel becomes smaller, so there is larger acceleration at a smaller r, this faster region of the jet eventually is opposed by the slower moving outer jet, which eventually develops into a shock. Therefore an oscillating radiation field in the funnel like region above PSD produces shocks in the jet. The cold jet solution does not harbor a shock very close to the jet base, but jets with higher initial temperatures can have an internal shock very close to the jet base (r s ∼ 10). We have shown that the shocks can collide with each other and create multiple shocks after collision giving rise to shock cascade. We also showed that the production of shocks in jet can also depend on the frequency of oscillation of x s orṁ sk . Oscillations with higher frequency are likely to produce jet shocks compared to oscillations with lower frequency. In this work, we showed how the radiation field of the accretion disk can affect a jet in a complicated manner. The radiation field of the disk can produce mildly relativistic to highly relativistic jets and we have also shown that disk oscillations give rise to the production and collision of shocks in jets.

ACKNOWLEDGMENTS
We thank the anonymous referee for valuable comments which improved the quality of the paper.

A. SUB-KEPLERIAN ACCRETION RATE AND THE SIZE OF PSD
The size of the PSD is related to the sub-Keplerian accretion rate (Vyas et al. 2015) x s = 64.8735 − 14.1476ṁ sk + 1.242ṁ 2 sk − 0.0394ṁ 3 sk (A1) We give a perturbation inṁ sk given aṡ which results in a sinusoidal oscillation in shock position with a frequency f s .ṁ sk0 is the mean value of accretion rate and A 1 and A 2 control the amplitude of oscillation inṁ sk . For example,ṁ sk0 = 7.5, A 1 = 3.6872 and A 2 = 1.313, results in a sinusoidal oscillation in shock position with mean shock position x s = 12 and amplitude of oscillation is A = 7 .

B. LUMINOSITY OF DISK
The luminosities of disk components are obtained by numerically integrating the emissivity over the disk volume.
Where i represents the disk component (i=SKD or PSD), and l i is the luminosity in units of Eddington luminosity. I i is the Synchrotron emissivity (erg cm −3 s −1 ), X (x s ) is the Comptonization parameter fitting function. X (x s ) = 1 for SKD and for PSD it given as (Kumar et al. 2014) X (x s ) = 0.659234 + 0.127851x s − 0.00043x 2 s − 1.13 × 10 −6 x 3 s (B4) The disk luminosities are approximated by the algebraic functions given as The radiative moments with relativistic transformations on the jet axis due to SKD are given as (Vyas et al. 2015). It may be noted that, In the following we present the expressions of radiative moments computed from the SKD (e. g., Q sk ) and PSD (e. g., Q ps ). One need to follow Fig. (1b) to get a proper reference of the following equations, where where S sk is a constant given as S sk = (9.22 × 10 33 e 4 Θ 3 0 βσ Tṁ 2 sk )/(πm 2 e m 2 p c 2 G 2 M 2 ). Here, γ sk is the Lorentz factor, ϑ i is the i-th contra variant component of three-velocity of accreting matter and l i s are direction cosines. The temperature at outer edge of SKD Θ 0 is taken to be Θ 0 = 0.2.ṁ sk represents the accretion rate of SKD in units of Eddington accretion rate Ṁ Edd = 1.44 × 10 17 M B /M gs −1 . ϑ sk is the radial velocity of the accretion disk estimated from free fall velocity (see Vyas et al. 2015). The specific angular momentum λ 0 at outer edge is taken to be λ 0 = 1.7.
The specific intensity of PSD (Chattopadhyay & Chakrabarti 2002b) is given by I ps = l ps L Edd /πA ps , where A ps and L Edd are surface area of PSD and Eddington luminosity, respectively and luminosity of PSD l ps is in units of L Edd . The exact form of radiative moments from PSD are Where S is a constant given as S = (1.3 × 10 38 l ps σ T )/(2πm e A ps GM ). Equations (C8-C14) are computed numerically.

D. ANALYTICAL ESTIMATION OF RADIATIVE MOMENTS
In simulation code, we use the analytical estimates of radiative moments (mentioned in C7), the explicit expressions for these analytical functions are given below.