\section{Introduction}
The Building Information Model (BIM)\cite{azhar11bim} is an instance of a populated data model of buildings. It constitutes of all the different kinds of data pertaining to a specific building. These data invariably or uniquely define the same building. It is not a temporal representation; however a static projection of the building that defines the building from various different perspectives across disciplines. BIM has replaced traditional onsite construction process in being a significant workload percentage of the AEC industry \cite{darko2020building}.

In the recent years, BIM has been extended to include specialised tools such as energy analysis, structural analysis for energy and space management \cite{volk2014building}. However, creation of accurate BIMs is time consuming and expensive \cite{ahmed2018barriers} if it is done for the existing buildings and the advantages it offers for energy efficiency calculations \cite{eleftheriadis2017life, ghaffarianhoseini2017application} are lost. Hence it is imperative that we look for scalable and low cost approximate BIM modeling for existing infrastructure in a city or a local community.

In this paper, first we create the notion of the prototype buildings which serve as the state of the art implementation of energy efficiency techniques given particular conditions of space, number of residents and energy consumption are considered. Secondly,
we give a method to suggest a better prototype alternatives for the existing building. We propose using Point-GNN for detecting the buildings and determining their bounding boxes and then suggesting the efficient prototypes for the detected buildings.

Contributions of the paper are as follows:
\begin{itemize}
    \item Introducing the novel prototype buildings for capturing the energy profiles.
    \item Scalable framework for energy efficiency analysis and recommendation.
    \item Use Point-GNN for identifying existing buildings and suggesting energy efficient prototypes for them
    \item We perform the experiments on completely synthetic data for demonstrating our proposed approach
\end{itemize}

The remainder of this paper is organized as follows. We review some of the important related studies in section \ref{sec:related}. We discuss the prototype building information model (BIM) in detail in section \ref{sec:prototypeBIM}. We then discuss out proposed method in section \ref{sec:methodology}. We then conclude with section \ref{sec:conclusion} and discuss the future works.

\section{Related Works}
\label{sec:related}
Climate change and impending global energy crisis has moved the scientists along with Architecture, Engineering and Construction (AEC) industry have called for wide use of BIM-based energy analysis for building performance assessment \cite{chen2017green}. Although creating BIMs for exisiting building and retrofitting them to improve the existing building management and maintenance has been done in practice, it is important to note a few caveats. In \cite{di2015bim}, the researchers have tried to update a school building. Even though results are significant, so is the time as well as expense of creating an accurate BIM.

Similar exercises were also carried out by the authors in \cite{bruno2018historic} for historic buildings where they added the diagnosis part for automating acquisition and update of knowledge base for heritage sites. Deconstruction of existing buildings were also taken up by the research community by building 3D models through devices such as kinect for generating low accurate and optimized 3D BIMs in real time \cite{volk2018deconstruction}. However, the problem of scalability still remained an open challenge as none of the methods could give us a quick answer for a larger community such as a city to go through the tedious and time consuming process of the creating individual BIMs for buildings.

There are three broad approaches that are taken to handle the tasks related to point clouds using Deep Neural Networks. The first approach takes considers the point cloud as a set of points. The Multi Layered Perceptron (MLP) is used to extract the features corresponding to each of the nodes. \cite{zaheer2017deep} and \cite{qi2017pointnet} are such networks that aggregate the feature vectors obtained by using neural networks for the set of points.

The other approach considers the point cloud as a grid. The points in the 3D grids units called voxels are considered in \cite{zhou2018voxelnet} and applies 3D convolutions onto it for object detection. Some approaches \cite{chen2017multi}, \cite{yang2018pixor} project the point cloud into 2D space to apply the 2D convolutions onto it. However, projections into grids has drawback of regularity that is added due to the regular nature of the grid, unlike the point cloud.

The third approach is to model the point cloud using an irregular structure of graphs. The graph neural networks in its every layer updates the vertex/edge features by aggregating them across their neighborhood\cite{wu2020comprehensive}, \cite{Joshi2022}. Unlike 2D/3D convolutions, graph convolutions does not reduce the number of points in every iterations. \cite{Point-GNN} first systematically converts the point cloud into a graph for object detection and also gives the mechanism for determining the bounding box of the detected objects.

\section{Prototype Building Information Model (BIM)}
\label{sec:prototypeBIM}

\subsection{3D Model}
3D model gives a visual feedback for the geometry the structure\cite{baz04vbe}. The BIM model is abstract which contains only information but it is not easy to visualize BIM. Since the BIM also contains information about the geometry of the building, they are easily accurately rendered to produce a visual information BIM. In fact, the 3D model is used to accurately find the geometry of the various structure in BIM, and not the other way around, since, obtaining the geometry is difficult humanly without a visual feedback. 3D model with various levels of details are shown in Figure \ref{fig:lod}.

\begin{figure}[h]
    \centering
    \includegraphics[width=\linewidth]{images/levelOfDetail.jpg}
    \caption{3D view of different level of details. (Top left) Original photo. (Bottom left) High level of detail. (Top right) Low level of detail. (Bottom right) Medium level of detail.}
    \label{fig:lod}
\end{figure}
\subsection{Data collection}
The first step for the creation of BIM for an entire city requires careful planning and collection of data. The BIM in itself is a very data intensive model. A assortment of BIMs for each building in a city increases the amount of data collected by many folds. Some of which may be redundant while many are not. But even without taking into factor the huge amount of data to be collected, there are many other factors that come to play.

First, it is equally challenging to figure out which data can be provided by whom and at what time. Even though a overall summary of the data maybe made available, the exact data points are not easy to be maintained and many legacy systems put in place in modern infrastructure do not even collect the data. Since, many systems both modern and legacy infrastructure is in place, there is always a risk of obtaining incomplete data and unbeknownst to the observer, the data might be scaled down and incorrect information creeps into the model.
Second, in a city, many private parties are involved who do not want to share data due to privacy and other concerns.

With the above facts taken into consideration, the data to be collected is then divided into two general category, the geometric and non-geometric data. The geometric data is concerned mostly with the visual look of the buildings and the non geometric data with any other information that the BIM is concerned with.

The geometric data contains the GIS location, structure, space, district/area around the building. The geometric data can be used for the visualization of BIMs and as the whole, the BIM integration of the entire city. The BIM of a given building is different from that of a city, as the later contains more information than the former. Along with the above information, the BIM of the city has information, for instance, the type of area the building is located in industrial, residential or commercial. The geometric data like the roads, bridges, and other public utilities are also information that not the part any given building but of the city. This will be helpful when deciding the prototypes of the building. And the data can be used to estimate the other non-geometric data as well.

The non-geometric data contains, information like mechanical data water requirements, electrical data like electric power consumption, HVAC data, safety or maintenance data. These non-geometric data is specific to the building's BIM, the BIMs for the entire city also needs data like the sewage, traffic, garbage (and waste disposal), drainage system etc, has to be collected.

Another type of data that is required for creating the BIMs of the entire city, which is not required by individual building BIM is, source. For instance these source include, power like the electric grid, the water tanks that provide the water for the city and the traffic load a road can take. Such data required for creating the model for the resource constraint.

Once, all the required data that has to be collected is determined, depending on how many parameters the model has to work on, the data has to be sanitized. The data sanitization step mostly consists of checking for erroneous data, missing or incomplete data, statistically improbable data, but there are other sanitization steps specific to the type of data. The data is also to be normalized, brought to reasonable range with respect to both time and space, convert all the different data units to one uniform unit.

The data collection is a continuous process. As the city grows, the new buildings are created and old buildings are either demolished or renovated, the model evolves with time as the new data enters the system. Optionally, the model should be able to make prediction about the BIM instead of only querying or suggest future development and planning. This adds new data to the system.

\subsection{Prototype Buildings}
There are many buildings in a city and each building is different from the rest. It is true that no two building are similar in all the features of data collected. BIMs generally deal with individual buildings and therefore, making a detailed model is not an issue while, for a city with creating BIMs for all the building is not only a huge task, it doesn't provide any additional benefit to the model. Any perceived advantage is further diminished, since it becomes a hassle and un-maintainable on a large scale, while not providing any significant advantage over a statistical measure.

A alternative approach is, based on the data, to create a different prototypes of buildings and use the prototypes of building to replace all the buildings with their respective BIMs. The prototypes maybe detailed with respect to a certain feature, while being very coarse in an different feature. A feature being coarse means that a prototype can be used to represent many buildings whose feature in question fall in a certain range. While, a detailed feature begs for an exact or near exact match of the feature.

Creating a prototype is not a routine part of the model. The prototype is created when the model is initially created and again when there is a necessity for it's creation. Initially, when BIMs for a city is created from the ground up, a bulk of the work goes into the creation of the prototypes since many of them has to be created. But when the model has been set up, there is not a need to create prototypes anymore. But as the model evolves with more data and new prototype buildings has to be made. But creation of new prototypes is few and far in between and not a regular part of the life cycle as illustrated in the Fig. \ref{fig:cycle}.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{images/cycle.png}
    \caption{Life cycle of BIM for cities}
    \label{fig:cycle}
\end{figure}

The creation of prototypes also depends on the coarse or the detailed feature of the data. For features that need exact or near exact match of the features of the data of the buildings with the prototype, a prototype for every such match is required. Therefore, these features must be prominent features and must broadly divide the buildings into a small limited category. For instance, the type of the building, whether residential, commercial or industrial is such a feature, but the power consumption is not such a feature, since for every value of power consumption, there has to be a prototype, which defeats the entire purpose of making prototypes.

Power consumption, however, is coarse feature for a prototype and therefore it's prototype must span over a large range of values. The ranges of these features ideally may be chosen appropriately so that the buildings in the city may fall as evenly as possible, but this may not be possible as the number of features increases.

There are many different types of buildings at NISER. The buildings are mostly categorized based from their top down view as shown in Table \ref{tab:niserbuild} . What is most important for a building is that it should retain an overall coarse structure for any air or drone traffic. Maintaining the detail structure of the building is not only adds a lot of data points to the storage but also adds a hit to the rendering of the buildings and also for the machine learning algorithm. Making a detailed structure to the building is a ``hard'' task. But making the structure of the building too coarse, implies that the building becomes a cuboid, which is what is to be avoided. A trade-off is needed and hence a top down view of the building is sufficient enough, while maintaining the overall structure. Also, for a machine learning algorithm it is easier to create the top-down view of the buildings.

\begin{table*}
\caption{Array of 3D models and their corresponding satellite images}
\label{tab:niserbuild}
\begin{tabular}{|c|c|c|c|}
		\hline
	\includegraphics[width=.22\linewidth]{./images/soh.png}&\includegraphics[width=.22\linewidth]{./images/SoH-m.png} &
	\includegraphics[width=.22\linewidth]{./images/doh.png}&\includegraphics[width=.22\linewidth]{./images/DoH-m.png}\\\hline
	\includegraphics[width=.22\linewidth]{./images/library.png}&\includegraphics[width=.22\linewidth]{./images/libray-m.png} &
	\includegraphics[width=.22\linewidth]{./images/admin.png}&\includegraphics[width=.22\linewidth]{./images/SoH-m.png}\\\hline
	\includegraphics[width=.22\linewidth]{./images/school-ext.png}&\includegraphics[width=.22\linewidth]{./images/school-ext-m.png} &
	\includegraphics[width=.22\linewidth]{./images/school-w-ext.png}&\includegraphics[width=.22\linewidth]{./images/school-w-ext-m.png}\\\hline
	\includegraphics[width=.22\linewidth]{./images/shopping.png}&\includegraphics[width=.22\linewidth]{./images/shopping-m.png} &
	\includegraphics[width=.22\linewidth]{./images/health.png}&\includegraphics[width=.22\linewidth]{./images/health-m.png}\\\hline
	\hline
\end{tabular}
\end{table*}

Many mapping services like Google Maps, Here Maps, OpenStreetMaps, Bing Maps etc, have top down views of the area where buildings can be easily made out. Looking at a top down view, the buildings at NISER can be categorised as such. The entire map with the protoypes at NISER for demonstration purpose is mapped in Figure \ref{fig:NISERcomplete}.

\begin{figure}
    \centering
    \includegraphics[width=\linewidth]{images/entire.png}
    \caption{Entire Map of NISER}
    \label{fig:NISERcomplete}
\end{figure}

\subsection{Storing building prototypes}

Building prototypes have to be saved in data structures. There were many candidates for the data structures that were considered. One way was to have other data in the building's geometry while the other was to incorporate the geometry in the building data structure, shown below or have a separate data structure for the buildings that includes geometry:
\begin{itemize}
    \item Building - vertices, edges, faces, power, water, traffic
    \item Building - Geometry (Vertices, edges, faces) , power, water, traffic
\end{itemize}

In the long term, the scalability of such a data structure is important. The data structure must be able to add new data categories and integrate well. Finally, keeping in mind that the ultimate goal is to have a learning algorithm to automate many of these tasks, it is imperative that the data structure should also facilitate modification automatically.

\section{eBIM-GNN}
\label{sec:methodology}

We are interested in a system that matches the right prototype to the buildings. Towards this goal,
we propose a two-step system to address this problem, where the first part detects the objects in the scene and classifies into one of the two classes, building or not a building. Along with this classification, the bounding boxes of each of the detected objects are predicted, these are used to localize the buildings. The second part considers the buildings detected by the first part and further assigns the closest prototype to it.

\begin{figure}
    \centering
    \subfloat{
        \includegraphics[width=\linewidth]{images/point-gnn.png}
    }

    \subfloat{
        \includegraphics[width=\linewidth]{images/prototype-pred.png}
    }
    \caption{Predicting prototypes corresponding to the detected objects}
    \label{fig:prototype-pred}
\end{figure}

\subsection{Graph Construction}
\label{subsection:prelim}
A point cloud is defined as a set $P$ of points, given as $P= \{ p_1, ... p_N \}$, where $p_i = ( x_i, s_i )$ is a point with $x_i \in \mathbb{R}^3$ are the co-ordinates of the points and $s_i \in \mathbb{R}^m $ is the feature vector of length $m$ corresponding to the point $i$. We construct a graph $G=(V, E)$ from the given point cloud. We consider each point from $P$ as the vertices in the graph. Therefore, $V = P$, and the edges in the graph $G$ we connect two nodes if they lie within a distance $r$ from each other. This is given in equation \ref{eq:edges}.

\begin{equation}
    \label{eq:edges}
    E = \{(p_i, p_j) \mid \norm{x_i-x_j}_2 \leq r \}
\end{equation}
where $r$ is the radius and $x_j$ and $x_j$ are the nodes in the graph.
We form a graph from the given point cloud as described previously. The graph formed using the entire point could might overburden the computation. Hence, before constructing the graph, the point cloud can be downsampled. Point-GNN makes use of voxel downsampling.

The initial feature vector $s_i$ of node $i$ is obtained by using $Max$ function to aggregate the embeddings generated by a multilayer perceptron (MLP) that embeds the relative coordinates and other features.

\subsection{Graph Neural Networks}
The Graph Neural Networks (GNNs) address the questions related to systems that can be modeled as graphs.
In each layer $l+1$ of a GNN, the features of the node's neighboring nodes are extracted by aggregating them, as shown in the equation \ref{eq:node_update}.

\begin{align}
    v_i^{l+1} &= g^l (\rho(\{ e_{ij}^l \mid (i, j) \in E\}), v_i^l) \label{eq:node_update} \\
    e_{ij}^{l} &= f^t(v_i^l, v_j^l) \label{eq:edge_update}
\end{align}
where $v^l$ and $e^l$ are feature vectors corresponding to vertex $i$ and edge $ij$ respectively. The set function $\rho$ aggregates the edge features for each vertex, function $g^l$ updates the next layer vertex features, and the function $f^l$ aggregates the vertex features to update the edge features.

\subsection{Point-GNN}
\subsubsection{Overview}
The method Point-GNN\cite{Point-GNN} makes use of GNNs for detecting objects given a point cloud, it also localizing the bounding box of the object. Poing-GNN uses an auto-registration mechanism that enables the points to align their coordinates based on their features. The auto-registration mechanism reduced the translation variance. The alignment offset is predicted using the structural features from the previous layers in a GNN. We modify the equation \ref{eq:node_update} to update the feature vector $s_i$ corresponding to each node $i$, and add auto-registration to it.
The auto-registration mechanism is as follows:
\begin{align}
    \Delta x_i^l &= h^l (s_i^l) \\
    s_i^{l+1} &= g^l \left( \rho \left( \left\{ f^l(x_j - x_i, s_i^l) \mid (i, j) \in E \right\} \right), s_i^l \right)
\end{align}
where $\Delta x_i^l$ is coordination offset, calculated by the function $h^l$ using the center node $i$'s feature vector $s_i^l$ from previous layer $l$.

The functions $f^l, g^l, $ and $ h^l$ are modeled using MLPs and $\rho$ is chosen to be $Max$ function. For each layer of the Point-GNN, the MLPs are different that are not shared across the layers. As the Point-GNN is used to handle two tasks, classification and localization, the representations produced by Point-GNN after $L$ layers is passed to two branches to handle each of the two tasks.

\subsubsection{Box Merging and Scoring}
In order to improve the localization of the bounding box, Point-GNN uses the modified Non-maximum suppression(NMS). Multiple bounding boxes can be predicted as there are many points on the same object. Along with the bounding box, the corresponding confidence score is calculated. It is done using the ratio of occupied volumen, called occlusion factor. A box $b_i$ has length $l_i$, width $w_i$ and height $h_i$ and three vectors $v_i^l, v_i^w, v_i^h$ are the unit vectors in their directions respectively. Then the occlusion factor is given as
\begin{equation}
    o_i = \frac{1}{l_iw_ih_i} \prod_{v \in {v_i^l, v_i^w, v_i^h}} \max_{p_j \in b_i}(v^T x_j) - \min_{p_j \in b_i}(v^T x_j)
\end{equation}

\subsection{Loss}

The classification branch outputs a multiclass probability distribution corresponding to each of the vertex. It is given by $\{ p_{c_1}, p_{c_2}, \dots, p_{c_M} \}$ for $M$ classes of objects.
The cross-entropy loss contributed by the first branch of classification, is given as follows:
\begin{equation}
    l_{cls} = - \frac{1}{N} \sum_{i=1}^{N} \sum_{j=1}^{M} y_{c_j}^i log(p_{c_j}^i)
\end{equation}
where $p_{c}^{i}$ and $y_{c}^i$ are the predicted probabilities and one-hot encoding of the predicted label of $i$th class.

The bounding box is predicted for every detected objects as a 7-tuple $b=(x, y, z, l, b, w, \theta)$ with 7 degrees of freedom. It consists of the coordinates of center $x, y, z$, length, width and height $l, h, w$ respectively and the yaw angle $\theta$ of the bounding box. The bounding box is encoded with vertex coordinates $(x_v, y_v, z_v)$ as follows:

\begin{align}
    \begin{split}
        \delta_x &= \frac{x-x_v}{l_m}, \delta_y = \frac{y-y_v}{h_m}, \delta_z = \frac{z-z_v}{w_m} \\
        \delta_l &= \log\frac{l}{l_m}, \delta_h = \log\frac{h}{h_m}, \delta_w = \log\frac{w}{w_m} \\
        \delta_\theta &= \frac{\theta - \theta_)}{\theta_m}
    \end{split}
\end{align}

    The loss from the localization branch that predicts the bounding box, is computed using Huber loss \cite{huber1964robust} on ground truth and prediction. The vertices that need not be localized or lie outside the bounding box contribute nothing to the loss. The average localization loss is given by:
\begin{equation}
    l_{loc} = \frac{1}{N} \sum_{i=1}^N \mathbb{I}(v_i \in b_{interest}) \sum_{\delta \in \delta_{b_i}} l_{huber}(\delta - \delta^{gt})
\end{equation}

\begin{ColorPar}{black}
\section{Prototype Matching}

We need to match the buildings with the prototypes, so that depending on the energy efficiency, a better prototype could be suggested for the given building. We propose two approaches to do this.
\end{ColorPar}

\begin{ColorPar}{black}
\subsection{Approach 1: Parameterized Aggregation}
The bounding box that inscribes the object is predicted for each of the identified object by point-GNN. We predefine the prototypes based on the prior knowledge about ideal buildings, their energy usage and energy efficiency. For each of the building type $b_i$, we consider a set of prototypes $P_{b_i}$, where $1 \leq i \leq B$ with total $B$ different types of buildings. For each of these prototypes, there are associated predefined specifications of each of their features. To $k$th feature $f_{jk}$ of the prototype $j \in P_{b_i}$ we assign a coefficient that denoted the contribution of the particular feature in deciding the better matching of the identified objects with prototypes.

We calculate the mean $\mu_l, \mu_w, \mu_h$ and standard deviation $\sigma_l, \sigma_w, \sigma_h$ of the length, width and height respectively of the identified bounding boxes for each of the classes. The prototypes are chosen with dimensions $(\mu_l, \mu_w, \mu_h), (\mu_l \pm \sigma_l, \mu_w \pm \sigma_w, \mu_h \pm \sigma_h)$ and $(\mu_l \pm 2\sigma_l, \mu_w \pm 2\sigma_w, \mu_h \pm 2\sigma_h)$ of for all different types of buildings.

To match the identified buildings to the prototypes, we consider the closest values of the dimensions of the prototypes to that of the buildings. The range for of the five prototypes corresponding to every building type is $(\mu_l - 2.5 \sigma_l, \mu_w - 2.5 \sigma_w, \mu_h - 2.5 \sigma_h) - (\mu_l - 1.5 \sigma_l, \mu_w - 1.5 \sigma_w, \mu_h - 1.5 \sigma_h)$ , $(\mu_l - 1.5 \sigma_l, \mu_w - 1.5 \sigma_w, \mu_h - 1.5 \sigma_h) - (\mu_l - 0.5 \sigma_l, \mu_w - 0.5 \sigma_w, \mu_h - 0.5 \sigma_h)$,
$(\mu_l - 0.5 \sigma_l, \mu_w - 0.5 \sigma_w, \mu_h - 0.5 \sigma_h) - (\mu_l + 0.5 \sigma_l, \mu_w + 0.5 \sigma_w, \mu_h + 0.5 \sigma_h)$,
$(\mu_l + 0.5 \sigma_l, \mu_w + 0.5 \sigma_w, \mu_h + 0.5 \sigma_h) - (\mu_l + 1.5 \sigma_l, \mu_w + 1.5 \sigma_w, \mu_h + 1.5 \sigma_h)$, and
$(\mu_l + 1.5 \sigma_l, \mu_w + 1.5 \sigma_w, \mu_h + 1.5 \sigma_h) - (\mu_l + 2.5 \sigma_l, \mu_w + 2.5 \sigma_w, \mu_h + 2.5 \sigma_h)$ respectively.

We can then match each of the identified building to its closest prototype.
The loss in matching the prototypes of each of the building type is calculated as follows:
\begin{equation}
    \label{eq:loss-matching}
    l_{pro} = \sum_{j=i}^{N} \sum_{i=1}^{F} \phi_i (f_{a}^{ij} - f_{mp}^{ij})
\end{equation}
where, $F$ is the number of features of the buildings, $\phi_i$ is the coefficient corresponding to $i$th features, $f_{a}^{ij}$ is the actual value of the $i$th feature of $j$th example and $f_{mp}^{ij}$ the value of that feature in the matched prototype. $\phi_i$ can be learnt when the data is annotated.

\end{ColorPar}

\begin{ColorPar}{black}
\subsection{Approach 2: Learning with MLP}
\end{ColorPar}
The second part of the system in figure \ref{fig:prototype-pred} classifies the detected objects into classes (prototypes). This network also gives a multi-class probabilities $\{ p_{p_1}, p_{p_2}, \dots, p_{p_T} \}$ corresponding to the $T$ prototypes for each of the detected objects. We compute the cross entropy loss as follows.

\begin{equation}
    l_{pro} = - \frac{1}{N} \sum_{i=1}^N \sum_{j=1}^T z_{p_j}^i log(p_{p_j}^i)
\end{equation}
where $z_{p_j}$ is the one-hot vector indicating the prototype corresponding to the example.

The total loss is given by:
\begin{equation}
        l_{total} = \alpha l_{cls} + \beta l_{loc} + \gamma l_{pro} + \kappa l_{reg}
\end{equation}
where $\alpha, \beta, \gamma, \kappa$ are the constants to balance losses.

\subsection{Efficiency}

With the use of eBIM-GNN, we get the better and efficient prototypes for the detected objects. To calculate the gain in efficiency, we compare the total energy consumption of the current building with the total energy consumption of the prototype building. The total energy consumption of the detected object is given by $E_o$, and the total energy consumption of suggested prototype is given by $E_p$. We find the percentage gain in total energy efficiency as given in equation \ref{eq:gain}.
\begin{equation}
\label{eq:gain}
    gain_{efficiency} = \frac{E_o - E_p}{E_o} \times 100
\end{equation}

\begin{ColorPar}{black}
\section{Evaluation on Synthetic Data}
\subsection{Dataset}
We consider the KITTI dataset \cite{geiger2012CVPR} that consists of point clouds and camera images. The dataset has 7481 training examples and 7518 testing examples.  Currently due to unavailability of the annotated data consisting of point clouds of buildings, we demonstrate the results on a similar annotated dataset with objects beings cars, pedestrians and cyclists instead of buildings. Since we are considering the KITTI dataset for demonstration, we associate the energy usage to the passenger capacity of each of the three categories. We consider the categories as building types for the sake of convenience. Category 'car' is considered as building type 1, 'pedestrian' as building type 2 and 'cyclist' as building type 3. Once identified, we consider the bounding boxes along with the identified class for to match it with the predefined prototypes.

We get five prototypes corresponding to each class. The ranges of dimensions of bounding boxes for matching the prototypes are as shown in the table \ref{tab:bounding-boxes}. These values are $-0.5 \sigma$ to $+0.5 \sigma$ to each of the dimensions of prototypes. We calculate the total efficiency, as well as the efficiency of each of the building types.

\begin{table*}[ht]
    \caption{Range of dimensions of prototype classes in $(h, w, l)$ format.}
    \label{tab:bounding-boxes}
    \centering
    \begin{tabular}{|c|p{3cm}|p{3cm}|p{3cm}|}
        \hline
        Prototype & Building type 1 & Building type 2 & Building type 3 \\
        \hline
        1 & (1.6971, 1.7051, 4.3538) - (1.8080, 1.7622, 4.6548) & (1.8890, 0.7882, 1.0638) - (1.9667, 0.8725, 1.2325) & (1.8250, 0.6876, 1.9565) - (1.8860, 0.7517, 2.0837)\\
        \hline
        2 & (1.5862, 1.6480, 4.0527) - (1.6971, 1.7051, 4.3538) & (1.8112, 0.7039, 0.8951) - (1.8890, 0.7882, 1.0638) & (1.7640, 0.6235, 1.8293) - (1.8250, 0.6876, 1.9565)\\
        \hline
        3 & (1.4753, 1.5909, 3.7518) - (1.5862, 1.6480, 4.0527) & (1.7334, 0.6196, 0.7264) - (1.8112, 0.7039, 0.8951) & (1.7030, 0.5594, 1.7021) - (1.7640, 0.6235, 1.8293)\\
        \hline
        4 & (1.3644, 1.5338, 3.4508) - (1.4753, 1.5909, 3.7518) & (1.6556, 0.5353, 0.5577) - (1.7334, 0.6196, 0.7264) & (1.6420, 0.4953, 1.5749) - (1.7030, 0.5594, 1.7021)\\
        \hline
        5 & (1.2535, 1.4767, 3.1498) - (1.3644, 1.5338, 3.4508) & (1.5778, 0.4510, 0.3890) - (1.6556, 0.5353, 0.5577) & (1.5810, 0.4312, 1.4477) - (1.6420, 0.4953, 1.5749) \\
        \hline
    \end{tabular}
\end{table*}

We define the energy usage for each of the building types randomly for our experiments. For building type 1 the energy usage is considered as 40 units, for building type 2 it is 20 units, and for building type 3 it is considered as 10 unit of energy. However, each of the prototype has different energy usage values associated with them. These values are as shown in table \ref{tab:ideal-usage}.

\begin{table}[ht]
    \caption{Predefined energy usage values for each of the building prototypes}
    \label{tab:ideal-usage}
    \centering
    \begin{tabular}{|c|r|r|r|r|}
        \hline
        Prototype & Building type 1 & Building type 2 & Building type 3 \\
        \hline
        1 & 40 & 10 & 20 \\
        \hline
        2 & 38 & 8 & 18 \\
        \hline
        3 & 36 & 6 & 16 \\
        \hline
        4 & 34 & 4 & 14 \\
        \hline
        5 & 32 & 2 & 12  \\
        \hline
    \end{tabular}
\end{table}

\subsection{Results}

\begin{table}[ht]
    \caption{Counts of buildings in the range of dimensions of respective prototype classes }
    \label{tab:counts}
    \centering
    \begin{tabular}{|c|r|r|r|r|}
        \hline
        Prototype & Building type 1 & Building type 2 & Building type 3 \\
        \hline
        1 & 76 & 0 & 0 \\
        \hline
        2 & 273 & 269 & 110 \\
        \hline
        3 & 1451 & 2104 & 145 \\
        \hline
        4 & 886 & 109 & 26 \\
        \hline
        5 & 1 & 21 & 1 \\
        \hline
    \end{tabular}
\end{table}

In order to calculate efficiency after matching with the prototype, we count the number of examples belonging to each of the prototypes. This is shown in table \ref{tab:counts} By multiplying the energy usage corresponding to each of the prototypes of every class and summing it up, we get the total energy usage. We show the energy usage in table \ref{tab:usage} along with the reduction in energy usage corresponding to each of the building type. The total reduction in energy usage is 16.22\%.

\begin{table}[ht]
    \caption{Percentage reduction in energy usage}
    \label{tab:usage}
    \centering
    \begin{tabular}{|c|r|r|r|}
        \hline
        Prototype & Actual usage & Prototype usage & Red. in usage (\%)\\
        \hline
        Building type 1 & 107480 & 95806 & 10.86\\
        \hline
        Building type 2 & 25030  & 15254 & 39.05\\
        \hline
        Building type 3 & 5640   & 4676  & 17.09\\
        \hline
    \end{tabular}
\end{table}

\end{ColorPar}

\section{Conclusion}
\label{sec:conclusion}
eBIM-GNN suggests the energy efficient prototypes for the buildings detected in a 3D point cloud. The building detection and determination of its bounding box is done using the method of Point-GNN. After the determination of bounding box, with the use of other features corresponding to the detected objects, the second part of system determines the appropriate prototype corresponding to the object. We demonstrate our approach on synthetic dataset.

In future, we would like to collect real data of the buildings and tune out the possible scalability issues related with the resolution of the point clouds and further optimize the matching against the prototypes. Our system benefits from the number of prototypes it has in the system to more accurately predict the energy efficiency gains.

\bibliographystyle{plain}


\section{Introduction}
\label{sec:intro}
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}\,{\rm erg\,s^{-1}}$ while the GRB jets can reach up to $10^{52}\, {\rm 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 \citep{hh06}. 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 \citep{gfp03,rsfp10, fgr10}, 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_{\rm g}$) of the accretion disk \citep{jbl99, dl12}. 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 \citep{1976ApJ...204..187S, ct95, 1997MNRAS.288..958G}. 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 \citep{wamkp14,m19,skr21}. 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 \citep{fpvj12, kp21,fmgp21}. 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 \citep{mm84,k98,p06,r07} and these equations have been used extensively in steady-state investigations of the radiatively driven outflows around compact objects. \cite{i89} used relativistic equation with near disk approximation of radiation field above a Shakura-Sunyaev disk \citep{ss73} to obtain an upper limit of jet terminal speed. \cite{f99} 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, \cite{fth01} showed that electron-proton jets can reach relativistic terminal speeds
by considering a hybrid disk with inner advection dominated flow (ADAF) \citep{nkh97} and an outer Keplerian disk. The disk model considered by \cite{ct95} 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 \citep{1976ApJ...204..187S,1997MNRAS.288..958G}. Numerical simulations of advective disks show that the extra thermal gradient term in hot post-shock region automatically generates the bipolar outflows \citep{dcnm14,lcksr16} which can be accelerated by the radiation field of the disk.
In the relativistic radiation hydrodynamic regime, investigations by \cite{c05} and \cite{vc17,vc18,vc19}
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 \citep{cc02a,csnr12,rvc21,jcl22}. To study the effect of radiation on winds there are some simulations of line driven winds \citep{psk2000,ybl18}, 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  \cite{vc17,vc18,vc19} and \cite{ftrt85}
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 \cite{ftrt85,vkmc15} 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 \cite{cr09}. 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 \citep{kss2000}, AGNs \citep{wamkp14}, as well as GRB jets \citep{rm94}. Shock very close to the jet base around a stellar mass BH was considered to explain the high-energy power law emission \citep{lrwc11}.
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 \citep{i89,f99,cc02,cdc04,c05}. Therefore if the disk radiation is time-dependent, 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 \ref{sec:equations} we present the assumptions and governing equations used in this paper. In section \ref{sec:method} we outline the methodology to obtain the solutions. Section \ref{sec:results} is devoted to explain and discuss the key findings of this paper. In the end, in section \ref{sec:concl} we conclude our work.

\section{Assumptions and governing equations} \label{sec: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_\phi=0$), axis-symmetric ($\partial/\partial\phi=0$), on-axis ($\theta=0$) jet. We assume a conical geometry of the jet with a narrow opening angle. The energy-momentum tensor for jet material ($T^{\alpha\beta}_M$) and radiation field ($T^{\alpha\beta}_R$) are given as

\begin{equation}
T^{\alpha\beta}_M= \rho h u^\alpha u^\beta +p g^{\alpha\beta};\,\,\, T^{\alpha\beta}_R= \int I_{\nu}l^{\alpha}l^{\beta} d\nu d\Omega
\label{eq:em_tensor}
\end{equation}

where $u^{\alpha}$ represent the components of four velocity, $l^{\alpha}$s are the direction cosines, $I_\nu$ is the specific intensity of the radiation field. $d\Omega$ is the solid angle subtended by the field point on the jet axis to the source point on the accretion disk. $\rho$ 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

\begin{equation}
h=\frac{(e+p)}{\rho}
\label{eq:enthalpy}
\end{equation}

The space-time metric is given by

\begin{equation}
ds^2=-g_{tt}dt^2+dr^2+r^2d\theta^2+r^2sin^2\theta d\phi^2
\label{eq:metric}
\end{equation}
Where $g_{tt}$ is given as;

\begin{equation}
g_{tt}=\left(1+\Phi\right); \,\,\, \Phi=-\frac{1}{2(r-1)}
\label{eq:gtt}
\end{equation}

We have used the geometric unit system where $2G=M_{\rm B}=c=1$, $M_{\rm 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 $\Phi$ is given by the Paczy\'nsky-Wiita potential \citep{pw80}, instead of the Newtonian one.
This way, we retain the conservative nature of equations of motion.
In this unit system, the Schwarzschild radius ($r_{\rm g}=2G M_{\rm B}/c^2$) is the unit of length and $t_{\rm g}=2GM_{\rm B}/c^3$ is the unit of time.

The equations of motion for relativistic radiation hydrodynamics have been derived before \citep{p06, r07}, 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
$(\rho u^\alpha)_{;\alpha}=0$.
The equations of motion are essentially the conservation of energy momentum tensor or, $(T^{\alpha\beta}_M)_{;\beta}=-(T^{\alpha\beta}_R)_{;\beta}=G^\alpha$.

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:

\begin{equation}
\frac{\partial D}{\partial t}+\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\sqrt{g_{tt}}Dv\right)=0
\label{eq:cont_cons}
\end{equation}

\begin{equation}
\frac{\partial M}{\partial t}+\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\sqrt{g_{tt}}(Mv+p)\right)+E\frac{\partial \sqrt{g_{tt}}}{\partial r}-\frac{2\sqrt{g_{tt}}p}{r}=G^r
\label{eq:mom_cons}
\end{equation}

\begin{equation}
\frac{\partial E}{\partial t}+\frac{1}{r^2}\frac{\partial}{\partial r}\left(r^2\sqrt{g_{tt}}(E+p)v\right)+(E+p)v\frac{\partial \sqrt{g_{tt}}}{\partial r}=G^t
\label{eq:eng_cons}
\end{equation}

Here $v$ is the three velocity of the jet. The four velocity ($u^r$) can be written in terms of $v$ as
$u^r=\gamma v$, where $\gamma$ is the Lorentz factor given as $\gamma^2=-u_tu^t=1/(1-v^2)$. The conserved quantities $D,\,M,\,{\rm and}\,E$ are the mass density, momentum density, and total energy density in the laboratory frame given as \citep[also see][]{rcc06,skrhc21}

\begin{eqnarray}
\label{eq:conservq1}
D=\gamma \rho \\
\label{eq:conservq2}
M= \gamma^2 \rho h v \\
E= \gamma^2 \rho h -p
\label{eq:conservq}
\end{eqnarray}

$G^r$ and $G^t$ are the components of the radiation force, which under the present set of assumptions are explicitly given as
\citep{mm84,k98,p06}

\begin{equation}
G^t=\frac{\gamma}{\sqrt{g_{tt}}}vR^r,~\&~ G^r=\gamma^3 R^r
\label{eq:gtfinal}
\end{equation}
where
\begin{equation}
R^r=\rho_e\frac{\sigma_{\rm T}}{m_{{\rm e}} c}\left[(1+v^2)F_{\rm rd}-v(E_{\rm rd}+P_{\rm rd})\right]=\rho_e\left[(1+v^2)\mathcal{F}-v(\mathcal{E}+\mathcal{P})\right]
\label{eq:radterm}
\end{equation}

Here $\rho_e$ is the total lepton density and $\mathcal{F}={\sigma_{\rm T} F_{\rm rd}}/({m_{{\rm e}} c}),\,\mathcal{E}={\sigma_{\rm T} E_{\rm rd}}/({m_{{\rm e}} c}),\,\mathcal{P}={\sigma_{\rm T} P_{\rm rd}}/({m_{{\rm e}} c})$, also, $E_{\rm rd},~F_{\rm rd},~\&~P_{\rm 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 (\ref{eq:cont_cons}-\ref{eq:eng_cons}) we need an additional closure relation which relates the thermodynamic variables $h, p$, and $\rho$. This closure relation is known as the equation of state (EoS). In this study, we consider an EoS for relativistic multispecies fluid proposed by \cite{cr09} (abbreviated as CR EoS). CR EoS is a very close fit to the exact one given by \cite{c39}. The CR has already been used for different astrophysical problems \citep{camrs14,sc19,scl20,jcry21,jcl22}.

The EoS is given as
\begin{equation}
e=\rho f
\label{eq:eos}
\end{equation}

where,
\begin{equation}
f=1+(2-\xi)\Theta\left[\frac{9\Theta+6/\tau}{6\Theta+8/\tau}\right]+\xi\Theta\left[\frac{9\Theta+6/\eta\tau}{6\Theta+8/\eta\tau}\right]
\label{eq:eos2}
\end{equation}
In equations (\ref{eq:eos}, \ref{eq:eos2}) $\rho$ is the mass density of fluid given as $\rho=\Sigma_i n_im_i=n_{{\rm e}^-} m_{{\rm e}} (2-\xi+\xi/\eta)$, where $\xi=n_{\rm p}/n_{{\rm e}^-}$, $\eta=m_{{\rm e}}/m_{\rm p}$ and
$n_{{\rm e}^-}$, $n_{\rm p}$, $m_{{\rm e}}$ and $m_{\rm p}$ are the electron number
density, the proton number density, the electron rest mass, and proton rest mass. $\Theta=p/\rho$ is a measure of temperature and $\tau=2-\xi+\xi/\eta$. The expression for sound speed is given as

\begin{equation}
a^2=\frac{1}{h}\frac{\partial p}{\partial \rho}=-\frac{\rho}{Nh}\frac{\partial h}{\partial \rho}=\frac{\Gamma\Theta}{h}
\label{eq:soundsp}
\end{equation}
Here, $N$ is the polytropic index given as

\begin{eqnarray}
N =  \rho \frac{\partial h}{\partial p}-1=\frac{\partial f}{\partial \Theta}
 =6\left[(2-\xi)\frac{9\Theta^2+24\Theta/\tau+8/\tau^2}{(6\Theta+8/\tau)^2}\right]+6\xi \left[\frac{9\Theta^2+24\Theta/(\eta \tau)+8/(\eta \tau)^2}{\{6\Theta+8/(\tau \eta)\}^2}\right]
\label{eq:poly}
\end{eqnarray}

The polytropic index is a function of temperature. Hence, we do not need to supply it as a free parameter. From equation (\ref{eq:poly}) one can see that for higher values of temperature $\Theta>>1$, $N \rightarrow 3$, and $N \rightarrow 3/2$ for $\Theta<<1$. The adiabatic index is related with the polytropic index as
\begin{equation}
\Gamma=1+\frac{1}{N}
\label{eq:adindx}
\end{equation}

\subsection{Steady state equations of motion}

In the steady-state, we impose ${\partial }/{\partial t}\equiv 0$ on equations of motion (\ref{eq:cont_cons}-\ref{eq:eng_cons}).
Integrating equation (\ref{eq:cont_cons}) then leads to the mass-outflow rate.

\begin{equation}
\dot{M}=\rho\gamma v\mathcal{A}
\label{eq:mdot}
\end{equation}
where $\mathcal{A}(\propto r^2)$ is the cross-section of the jet.
The energy balance equation or 1$^{\rm st}$ law of thermodynamics
is given by $u_\mu(T^{\mu \nu}_M)_{;\nu}=-u_\mu(T^{\mu \nu}_R)_{;\nu}$ which gives the temperature gradient
\begin{equation}
\frac{d \Theta}{d r}=-\frac{\Theta}{N}\left[\frac{2}{r}+\frac{\gamma^2}{v}\frac{dv}{dr}\right]
\label{eq:dthdr}
\end{equation}
Integrating the above and combining with equation \ref{eq:mdot} gives the adiabatic relation
\citep[also, see][]{kscc13,jcl22}
\begin{equation}
\dot{\mathcal{M}}=\gamma v r^2 \Theta^{3/2}(3\Theta+4/\tau)^{k_1}(3\Theta+4/\eta\tau)^{k_2}{\rm exp}(k_3)
\label{eq:scriptm}
\end{equation}

where
\begin{equation}
k_1=\frac{3}{4(2-\xi)},\, k_2=\frac{3\xi}{4},\, k3=-\frac{3}{\tau}\left[\frac{2-\xi}{3\Theta+4/\tau}+\frac{\xi}{\eta(3\Theta+4/\eta\tau)}\right]
\end{equation}
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 [$(\delta^r_\mu+u^ru_\mu)(T^{\mu \nu}_M+T^{\mu \nu}_R)_\nu=0$] can be simplified to the following form,
\begin{equation}
\gamma^2\frac{dv}{dr}=\frac{{2a^2}/{r}-\frac{1}{2(1+\Phi)}\frac{d\Phi}{dr}+{\gamma R^r}/{\rho h}}{v-\frac{a^2}{v}}
\label{eq:std_eul}
\end{equation}
where, the last term in the numerator $\frac{\gamma R^r}{\rho h}$ is the radiation momentum deposition term, which can be simplified to the form
\begin{equation}
\frac{\gamma R^r}{\rho h}=\frac{\gamma(2-\xi)}{\tau h}\left[(1+v^2)\mathcal{F}-v(\mathcal{E}+\mathcal{P})\right]
\label{eq:rad_mom_deps}
\end{equation}
We can integrate equations (\ref{eq:std_eul}) and (\ref{eq:dthdr}) to obtain another constant of motion for steady flow known as the generalized relativistic Bernoulli parameter (${\cal B}_{\rm e}$).

\begin{equation}
{\cal B}_{\rm e}=-hu_t{\rm exp}(-X); \,\,\, X=\int \frac{\gamma R^r}{\rho h} dr
\label{eq:genbp}
\end{equation}

\begin{figure}
\gridline{\fig{fig1a.pdf}{0.4\textwidth}{(a)}
          \fig{fig1b.pdf}{0.6\textwidth}{(b)}}
\caption{(a) Schematic representation of the disk-jet system. The shock location ($x_{\rm s}$), and various components of disk-jet system are shown. The blue region marks the funnel above the PSD. (b)- Cartoon diagram to show the cross sectional view of the disk. The blue portion is inner compact corona. SKD is shown with solid red lines. N is the field point on jet axis where radiative moments are calculated. M and ${\rm M^\prime}$ are the source point on the corona and SKD.}
\label{fig:cartoon}
\end{figure}

We consider an advective type accretion disk \citep{f87,c89,lrc11,2017MNRAS.469.4221K} for this study. Figure \ref{fig:cartoon}(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 \citep{scl20}. 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_{\rm 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 \citep{lcksr16} 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 $\theta_{\rm sk}=85^{\rm o}$ for SKD and intercept of SKD on jet axis is taken to be $d_0=0.4h_{\rm s}$, where $h_{\rm s}$ is the height of corona taken as $h_{\rm s}=2.5x_{\rm s}$. The outer edge of SKD is taken at $x_{\rm o}=3500$. The details of obtaining the radiative moments are presented in \cite{cc00,c05,vkmc15}, here we present only a brief description. \\

In Fig.\ref{fig:cartoon}(b) we have shown a cross sectional view of disk and jet sytem. $M$ and $M^\prime$ are the source points on PSD and SKD, respectively. $TM$ and $T^\prime M^\prime$ 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_{\rm s}$ but due to the shadow effect \citep{c05}, the observer at $N$ sees the SKD starting from $x_{\rm i}$ given as

\begin{equation}
x_i (r)=\frac{r-d_0}{(r-h_{\rm s})/x_{\rm s}+{\rm cot}\theta_{\rm sk}}
\label{eq:skd_in_edge}
\end{equation}

We need the intensities of various disk elements to obtain the radiative moments of respective disk components. In SKD, we assume the synchrotron emission \citep{st83} to be the dominant emission mechanism.

\begin{equation}
I_{\rm sk}=\left[\frac{16}{3}\frac{e^2}{c}\left(\frac{eB_{\rm sk}}{m_{{\rm e}} c}\right)^2\Theta_{\rm sk}^2n_{\rm sk}\right]\frac{\left(d_0\,{\rm sin\theta_{sk}}+x\,{\rm cos\theta_{sk}}\right)}{3}\,{\rm erg\, cm^{-2} s^{-1}}
\end{equation}

where $B_{\rm sk},\,\Theta_{\rm sk},\,n_{\rm 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_{\rm mag}=\beta{p_{\rm gas}}$.
The general definitions of radiative moments are
\begin{equation}
 E_{\rm rd}=\frac{1}{c}\int I d\Omega,~~ F^i_{
 \rm rd}=\int I l^i d\Omega, \mbox{ and } P^{ij}_{\rm rd}=\frac{1}{c}\int I l^i l^j d\Omega.
 \label{eq:moments}
\end{equation}
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 \ref{eq:num_mom}. The details of computing these numerically are given in \citet{vkmc15} and also in Appendix \ref{sec:app1} 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 \ref{eq:skd_ez}-\ref{eq:psd_pz}) with approximate analytical functions (equations \ref{eq:eng_alb}-\ref{eq:pres_alb} for SKD and \ref{eq:eng_alb_ps}-\ref{eq:pres_alb_ps} for PSD in Appendix \ref{sec:app2}) which depend on disk parameters like $x_{\rm s}, ~h_{\rm s} ~\& ~ x_{\rm in}$. The sub-Keplerain accretion rate controls the position of $x_{\rm s}$ (equation \ref{eq:xs_mdot}) and luminosities (equation \ref{eq:skd_lum_func}-\ref{eq:psd_lum_func}) 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_{\rm s}=15$ is shown in Fig.(\ref{fig:Rad_moms}). The fitting is good, and the advantage is, any time variation of $\dot{m}_{\rm 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,
${\cal E},~{\cal F}~\mbox{and}~{\cal P}$
obtained from equation \ref{eq:eng_alb}--\ref{eq:pres_alb}, and \ref{eq:eng_alb_ps}--\ref{eq:pres_alb_ps} for a given $\dot{m}_{\rm sk}$ are used in equation \ref{eq:radterm}.

\begin{figure}[!h]
 \centerline {\includegraphics[width=14cm,height=14cm, keepaspectratio]{fig2.pdf}}
  \caption{The radiative moments from PSD and SKD (excluding the constant factors). Lines (solid, dashed, dashed-dotted) represent the correct numerically integrated moments and open circles (red, green and cyan) represent
  moments obtained by fitting functions.}
  \label{fig:Rad_moms}
\end{figure}

\section{Method of obtaining solutions}
\label{sec:method}

\subsection{Steady state solutions}
The steady-state solutions are obtained by simultaneously integrating equations (\ref{eq:std_eul}) and (\ref{eq:dthdr}). 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\rightarrow0/0$ gives the sonic point conditions

\begin{equation}
v_c=a_c
\label{eq:spcond1}
\end{equation}

\begin{equation}
\frac{2a_c^2}{r_c}-\left(\frac{1}{2(1+\Phi)}\frac{d\Phi}{dr}\right)_c+\left(\frac{\gamma R^r}{\rho h}\right)_c=0
\label{eq:spcond2}
\end{equation}

From equation (\ref{eq:spcond2}) we obtain the temperature ($\Theta_c$) at the sonic point and $dv/dr$ is obtained by using ${\rm L'H\hat{o}pital's}$ rule. We integrate the equations (\ref{eq:std_eul}) and (\ref{eq:dthdr}) inwards and outwards starting from sonic points to obtain the solutions. The steady-state solutions provides the injection values for the simulation code.

\subsection{Time dependent solutions and numerical setup}
We perform the time-dependent calculation using a relativistic total variation diminishing (TVD) routine.
TVD scheme introduced by \cite{h83} 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 \citep[see][for details]{mm03, mm15}. A relativistic TVD code to incorporate CR EoS has been built before \citep{rcc06, crj13}. The detailed description of the code is presented in \citep{rcc06}, here we present only brief details. \\

The set of equations (\ref{eq:cont_cons})-(\ref{eq:eng_cons}) can be written in the form

\begin{equation}
\frac{\partial \vb*{q} }{\partial t}+\frac{1}{r^2}\frac{\partial (r^2 \vb*{F})}{\partial r}=\vb*{S}
\label{eq:rhd_conservative}
\end{equation}
with the state and flux vectors given as

\begin{equation}
\vb*{q}=\begin{pmatrix} D \\ M \\ E	 \end{pmatrix},\,\,\,\,\, \vb*{F}=\begin{pmatrix} \sqrt{g_{tt}}D v \\ \sqrt{g_{tt}} (Mv+p) \\ \sqrt{g_{tt}}(E+p)v \end{pmatrix}
\label{eq:state_vect}
\end{equation}

and the source term $\vb*{S}$ as

\begin{equation}
\vb*{S}=\begin{pmatrix} 0 \\ G^r-E\frac{\partial \sqrt{g_{tt}}}{\partial r}+\frac{2\sqrt{g_{tt}}p}{r} \\ G^t-(E+p)v\frac{\partial \sqrt{g_{tt}}}{\partial r}	 \end{pmatrix}
\label{eq:source}
\end{equation}

The set of equations \ref{eq:rhd_conservative} is solved in two steps. In the first step, the state vectors $\vb*{q}$ at the cell center $i$ at a given time step $n$ are updated using the modified fluxes $\bar{f}_{i+1/2}$ calculated at cell boundaries.

\begin{equation}
\vb*{L}\vb*{q}^{n}=\vb*{q}^{n-1}-\frac{\Delta t^n}{\Delta x}\left(\bar{f}_{i+1/2}-\bar{f}_{i-1/2}\right)
\end{equation}
where,
\begin{equation}
 \bar{f}_{i+1/2} = \frac{1}{2}\left[\vb*{F}^n_i
+\vb*{F}^n_{i+1}\right]-\frac{\Delta x}{2\Delta t^n}
\sum_{k=1}^3\beta_{k,i+1/2}\vb*{R}_{k,i+1/2}^n,
\end{equation}
where,
\begin{equation}
 \beta_{k,i+1/2} = Q_k\left(\frac{\Delta t^n}{\Delta x}\lambda_{k,i+1/2}^n
+\gamma_{k,i+1/2}\right)\alpha_{k,i+1/2}-\left(g_{k,i}+g_{k,i+1}\right).
\end{equation}
Here,
$$
\alpha_{k,i+1/2} = \vb*{L}_{k,i+1/2}^n\cdot\left(\vb*{q}_{i+1}^n
-\vb*{q}_i^n\right)
$$

Explicit forms of eigenvalues $\lambda_{k,i}$, eigenvectors $\vb*{L}~\&~\vb*{R}$, flux limiters $g_{k,i}$, $\gamma_{k,i+1/2}$, $Q_k$ are given in \cite{rcc06}. 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 ($\rho, v, p$) from the updated conserved variables $D,~M,~E$ which is done by inverting the equations \ref{eq:conservq1}-\ref{eq:conservq}. These equations along with  EoS can be expressed as a polynomial of $\gamma$, given as

\begin{equation}
\sqrt{\gamma^2-1}D+g_1\gamma\left[\frac{(2-\eta\tau)}{2(1-\eta)}\left(\frac{9\gamma g_1+6D\sqrt{\gamma^2-1}/\tau}{3\gamma g_1+4D\sqrt{\gamma^2-1}/\tau}\right)+\eta\frac{(\tau-2)}{2(1-\eta)}\left(\frac{9\eta\gamma g_1+6D\sqrt{\gamma^2-1}/\tau}{3\eta\gamma g_1+4D\sqrt{\gamma^2-1}/\tau}\right)+1\right]-M=0
\label{eq:lofac_polyn}
\end{equation}
where $g_1=\gamma M-E\sqrt{\gamma^2-1}$. We solve equation \ref{eq:lofac_polyn} using Newton-Raphson method to obtain $\gamma$, and once $\gamma$ is known the primitive variables can be obtained as
\begin{equation}
v=\frac{\sqrt{\gamma^2-1}}{\gamma}
\label{eq:vcal}
\end{equation}

\begin{equation}
\rho=\frac{D}{\gamma}
\label{eq:rhocal}
\end{equation}

Since $v$, $\gamma$, $E$ and $M$ are known, then using equations \ref{eq:conservq2} and \ref{eq:conservq}, the pressure $p$ is obtained as the root of a cubic equation

\begin{eqnarray}
a_1p^3+a_2p^2+a_3p+a_4=0, \\ \mbox{where, }
a_1=27\eta, \,\, a_2=9\left[\frac{4}{\tau}(1+\eta)\rho-\eta(E-Mv)\right] \nonumber \\
a_3=\frac{12\rho}{\tau}\left[\frac{2\rho}{\tau}+(1+\eta)\rho-\eta(E-Mv)\right], \,\,\,
a_4=\frac{16\rho^2}{\tau^2}[\rho-E+Mv] \nonumber
\label{eq:pres_recov}
\end{eqnarray}

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.

\subsubsection{Time dependent disk as observed by an inertial observer on the axis}
\label{sec:disc_inertial}
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 \citep{ndmc12}. If there is a fluctuation of the accretion rate $\dot{m}_{\rm sk}$ at the outer edge of the disk, the fluctuation would travel inwards. The time variability of the accretion rate $\dot{m}_{\rm sk}$
is given in equation \ref{eq:accretvar}, which induces a sinusoidal
variation in the size of PSD or $x_{\rm s}$, is of the form,
\begin{equation}
x_{\rm s}^\prime=x_{\rm s0}+A\,{\rm sin} (2\pi f_{\rm s} t)
\label{eq:shockt}
\end{equation}
where the oscillation is about a mean $x_{s0}$ (distance OG fo Fig. \ref{fig:retdpos}) with an amplitude of $A$ (NG=GN$^\prime$), as shown in Fig.(\ref{fig:retdpos}) and frequency $f_{\rm 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 \citep[see][]{jcl22}.
In other words, the radiation emitted from $x_{\rm s}^{\prime\prime}$ at $t$ reaches $P$ after a time $t+\delta t$ ($\delta t=l^{\prime\prime}/c$, see Fig. \ref{fig:retdpos}), and in the same time-interval the PSD boundary will move from $x_{\rm s}^{\prime\prime}$ to $x_{\rm s}^\prime$. Therefore,

\begin{equation}
\delta t=l^{\prime\prime}/c=\sqrt{{x_{\rm s}^{\prime\prime}}^2+(r-h_{\rm s}^{\prime\prime})^2}  \mbox{;  in units of } c=1
\label{eq:delta_t}
\end{equation}

The position $x_{\rm s}$ is updated after each time interval $dt$ which is determined by the TVD code.
So we can write
\begin{equation}
x_{\rm s}^{\prime \prime}-x_{\rm s}^{\prime}=s dt
\label{eq:pos_diff}
\end{equation}
Where $s=0.5(v_{{\rm s},\,t}+v_{{\rm s},\,t+dt})$ is the average velocity in the time interval between $t$ to $t+dt$. The instantaneous velocities can be obtained from equation \ref{eq:shockt},

\begin{equation}
 v_{\rm s}=A2\pi f_{\rm s} {\rm cos}(2\pi f_s t)
 \label{eq:instvelsok}
\end{equation}

From equations (\ref{eq:delta_t}) and (\ref{eq:pos_diff})
we can write
\begin{equation}
(x_{\rm s}^{\prime\prime}-x_{\rm s}^{\prime})^2=s^2\left[{x_{\rm s}^{\prime\prime}}^2+(r-h_s^{\prime\prime})^2\right]
\label{eq:ret_pos}
\end{equation}

We solve equation (\ref{eq:ret_pos}) at each $r$ and at every time step to obtain the $x_{\rm s}$ as seen by the observer there, and the corresponding radiative moments are calculated using the analytical estimates of the radiative moments (equations \ref{eq:eng_alb}-\ref{eq:pres_alb_ps}).

\begin{figure}[!h]
 \centerline {\includegraphics[width=15cm,height=6cm]{fig3.pdf}}
  \caption{Schematic representation of time-dependent disk configuration as seen by an observer P on jet axis. $x_{s0}$ is the mean position of shock and $x_{\rm s}^\prime$ is the actual position of shock at time t.}
  \label{fig:retdpos}
\end{figure}

\section{Results}
\label{sec:results}
\subsection{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 \citep[see][]{jcry21}.
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. \ref{fig:steady_comp} we have shown the solutions obtained from the TVD code marked by open red circles and analytical solutions shown by solid black lines. The evolution of various variables like flow velocity $v$, Mach number $M~(\equiv v/a)$, $\Theta$ which is the measure of temperature, and adiabatic index $\Gamma$ along the jet is shown in panels (a)-(d), respectively. The SKD accretion rate $\dot{m}_{\rm sk}$ is taken to be 10, which produces a PSD with $x_{\rm s}=8.212$ and disk luminosity is $l=1.0$. The SKD accretion rate $\dot{m}_{\rm sk}$ and disk luminosity $l$ are in the units of the Eddington accretion rate ($\dot{M}_{\rm Edd}=1.44\times10^{17}M_{\rm{B}}/\rm{M}_\odot \,{\rm g/s}$) and the Eddington luminosity ($\dot{L}_{\rm Edd}=1.26\times10^{38}M_{\rm{B}}/\rm{M}_\odot \, {\rm erg /s}$) respectively. The injection parameters with which the jet is launched are chosen from the analytical solution and are given by $v_{\rm in}=0.065,\,\,\Theta_{\rm in}=0.037,\,\, r_{\rm 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 $1000 r_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.

\begin{figure}
\centerline {\includegraphics[width=12cm,height=12cm, keepaspectratio]{fig4.pdf}}
\caption{Comparison between the steady state solutions obtained from TVD code (red open circles) and analytical method (solid black lines), for $v_{\rm in}=0.065,\,\,\Theta_{\rm in}=0.037,\,\, r_{\rm in}=5.0$. The jet is electron-proton or $\xi=1.0$ fluid.}
\label{fig:steady_comp}
\end{figure}

\subsection{Effect of disk luminosity on jet solution}

In Fig. \ref{fig:lum_comp}a, we compared solutions of electron-proton ($\xi=1$) jets corresponding to different disk luminosities, all starting with the same injection parameters
$v_{\rm in}=0.018,\,\Theta_{\rm in}=0.063,\,r_{\rm 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_{\rm T}=0.13c$. It is expected that higher disk luminosity would
accelerate the jets to higher $v_{\rm T}$. And indeed Fig. \ref{fig:lum_comp}b shows, $v_{\rm 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_{\rm T}\sim 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. \ref{fig:Rad_moms}, it is clear that the radiative
flux ${\cal F}_{\rm ps}<0$ in the funnel like region between the
inner walls of the PSD. So with higher $l$, the ${\cal F}_{\rm ps}$ will actually suppress the jet speed within the funnel region, although at higher altitude would push the jet. Therefore $v_{\rm T}$ increases with $l$.

\begin{figure}
\centerline {\includegraphics[width=13cm,height=13cm, keepaspectratio]{fig5.pdf}}
\caption{The effect of radiation on the jet solutions. Velocity profiles for various solutions with different disk luminosities are plotted in panel(a) and panel (b) shows the terminal velocities.}
\label{fig:lum_comp}
\end{figure}

\subsection{Effect of jet composition}
\label{sec:comp_effect}
\begin{figure}[!h]
\centerline {\includegraphics[width=13cm,height=13cm, keepaspectratio]{fig6.pdf}}
\caption{(a)-The velocity profiles for different plasma compositions of the jet. Panel (b) shows the variation of terminal Lorentz factor with respect to composition parameter for various disk luminosities.}
\label{fig:comp_effect}
\end{figure}

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 (\ref{eq:radterm}). To study the effect of composition we fix injection and disk parameters and the solutions are obtained for various values of $\xi$. In Fig. (\ref{fig:comp_effect} a) we have plotted the velocity profiles of various jet solutions, with the injection parameters
$v_{\rm in}=0.11,\, \Theta_{\rm in}=0.22, \, r_{\rm in}=2.0$. The disk luminosity is $l=0.82$. The radiation term $R^r$ for solutions with higher fraction of protons (higher $\xi$) 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 turn out to be faster at lower $r$. So the electron-proton jet ($\xi=1.0$), plotted with dash-dotted blue line is fastest near the jet base. The composition parameter significantly affects the terminal velocities which can be clearly seen from Fig.(\ref{fig:comp_effect})(a). The terminal velocity for electron-positron flow ($\xi=0$), plotted with dashed-green line is $v_{\rm T}=0.999876$ and the corresponding terminal Lorentz factor is $\gamma_{\rm T}=63.51$. In the right panel we have plotted the variation of terminal Lorentz factor $\gamma_T$ with respect to the jet composition for luminosities $l=0.82$ (dashed-red), $l=0.66$ (dash-dotted blue), and $l=0.54$ (solid-blue).
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.

\subsection{Jet solutions under the time dependent radiation field}
\begin{figure}[!h]
\centerline {\includegraphics[width=24cm, height=11cm]{fig7.pdf}}
\caption{The evolution of jet solution under the influence of disk oscillations. The variation of velocity (${\rm a_1}-{\rm a_5}$), Mach number (${\rm b_1}-{\rm b_5}$), radiative term (${\rm c_1}-{\rm c_5}$) is plotted at different epochs. The steady state solution (blue, solid) is obtained for steady luminosity $l=0.597$ and steady PSD configurtation at $x_{{\rm s}0}=12$. The PSD oscillation is for amplitude $A=7$, time period of $T=2000t_{\rm g}$ about an equilibrium value $x_{{\rm s}0}$ of the size of PSD.}
\label{fig:Shockcolsn}
\end{figure}

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. \ref{fig:Shockcolsn}
\& \ref{fig:collision} the injection parameters are
$v_{\rm in}=0.076,\,\,\,\Theta_{\rm in}=0.064,\,\,\,r_{\rm 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_{\rm s}=12.0$ (obtained from equation \ref{eq:xs_mdot}). 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. \ref{fig:Shockcolsn} ($\rm a_1$), ($\rm b_1$), ($\rm 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
$\dot{m}_{\rm sk}$ (equation \ref{eq:accretvar}) such that it imposes an oscillation on $x_{\rm s}$ with the amplitude $A=7.0$, about the equilibrium value $x_{{\rm s}0}=12$ and the time period being $T=2000 t_{\rm g}$. In this study we have taken the frequency and ampltiude of oscillation as free parameters. As $\dot{m}_{\rm sk}$ decreases, $x_{\rm s}$ moves outward and the net disk luminosity decreases resulting in a weaker radiation field (compare blue and red curves in Figs. \ref{fig:Shockcolsn}c$_1$-c$_5$). The  magnitude of net radiative force term $|R^r|$ for $x_{\rm s}>x_{{\rm s}0}$ is smaller than the steady-state values (blue) and is larger for $x_{\rm s} < x_{{\rm s}0}$. So in the funnel region, $v$ is higher when $x_{\rm s} > x_{{\rm s}0}$ but eventually produces a jet whose $v_{\rm T}$ is weaker and vice versa. Figures \ref{fig:Shockcolsn} 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. \ref{fig:Shockcolsn}($\rm 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 ${\rm S_1,\,S_2,\,S_3}$ as shown in Fig.  \ref{fig:Shockcolsn}($\rm a_4$). The shock $\rm S_1$ being faster, catches up and collides with the shock in front of it ($\rm S_2$), shown in Fig. \ref{fig:Shockcolsn}($\rm a_5$). In Fig. \ref{fig:collision} we have plotted the zoomed region of jet domain $300<r<800$ to highlight the changes in mass density ($\rho$) distribution due to collision of the shocks. Panels
 \ref{fig:collision}$\rm(e_1)-(e_{5})$ show two shocks S$_1$ and S$_2$ collide and are plotted in the epoch between Fig.\ref{fig:Shockcolsn}$\rm (a_4)-(a_5)$. Figures \ref{fig:collision}$\rm(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.

\begin{figure}[!h]
\centerline {\includegraphics[width=20cm, height=10cm]{fig8.pdf}}
\caption{Density distribution $\rho$ of the jet at vaious time steps are plotted to the collision between shock $\rm S_1$ and $ \rm S_2$ in panels $\rm (e_1-e_5)$. Panels $\rm (e_6-e_{10})$ show the $\rho$ after the collision. The PSD oscillation is same as in Fig. \ref{fig:Shockcolsn}.}
\label{fig:collision}
\end{figure}

In Fig. \ref{fig:pos_ratio} (a), (b) we plot the variation of positions of jet shocks $\rm S_1$ and $\rm S_2$ and their compression ratios $R$ as a function of time $t$. Figure \ref{fig:pos_ratio} (a) shows that the slope of $r_{\rm sh}$ vs $t$ graph is steeper for the shock $\rm S_1$, which clearly suggests that the speed of $\rm S_1$ is faster than that of $\rm S_2$ and ultimately catches up with $\rm S_2$ at $r_{\rm sh} \sim 535$. The compression ratio is the ratio of post-shock and pre-shock densities ($R=\rho^+/\rho^-$) 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 $\rm S_1$ becomes stronger as it moves outwards, reaching a maximum value of $R\sim4.4$ and gradually starts to slow down. The compression ratio for second shock $\rm S_2$ reduces initially but starts to increase to reach a maximum value of $R\sim4.55$ and then it becomes weaker.  \\

\begin{figure}[!h]
\centerline {\includegraphics[width=14cm, height=7cm]{fig9.pdf}}
\caption{The variation of shock positions and the compression ratio $R$ is plotted in panels (a) and (b), respectively.
Jets as described in Figs. \ref{fig:Shockcolsn} \& \ref{fig:collision}.}
\label{fig:pos_ratio}
\end{figure}

\begin{figure}[!h]
\centerline {\includegraphics[width=22cm, height=9cm]{fig10.pdf}}
\caption{The time dependent jet solution corresponding to injection parameters $v_{\rm in}=0.14,\,\,\,\Theta_{\rm in}=0.226,\,\,\,r_{\rm in}=2.0$. }
\label{fig:int_shock}
\end{figure}

A closer look at equation \ref{eq:rad_mom_deps} 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_{\rm in}=0.14,\,\,\,\Theta_{\rm in}=0.226,\,\,\,r_{\rm in}=2.0$, considered from the analytical jet solution for a steady disk configuration $x_{\rm s}=11.02$ and $l=0.664$. The velocity ($v$) and Mach number ($M$) variation for the steady-state solution is plotted in Figs. \ref{fig:int_shock}($\rm a_1$) and ($\rm b_1$) with the solid blue line. The amplitude of oscillation for $x_{\rm s}$ is $A=6.0$ and the time period of oscillation is $T=500 t_{\rm g}$. In Figs. \ref{fig:int_shock}a$_1$---b$_5$ we have plotted jet solution for different values of $x_{\rm s}$ as $\dot{m}_{\rm sk}$ is varied with time. So the variation of $\dot{m}_{\rm sk}$ determines the $l$ and $x_{\rm 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_{\rm sh} \sim10$) as shown in Fig. \ref{fig:int_shock}(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 above the PSD.
It may be noted that such shocks in jets have been invoked to explain high energy emissions \citep{jfk08,lrwc11}.

\subsection{Effect of disk oscillation frequency on jet solutions}
\begin{figure}[!h]
\centerline {\includegraphics[width=22cm, height=12cm]{fig11.pdf}}
\caption{The evolution of jet solutions ($v$ vs $r$) under the influence of disk oscillations of frequencies $f_s=2000t_{\rm g}$ (solid, red)
and $10^5t_{\rm g}$ (dashed, green). The injection and disk parameters are similar to Fig.(\ref{fig:Shockcolsn}). The two jet solutions in each panel corresponds to same $x_{\rm s}$ and $l$, and the time in which the disk configuration is similar for the two cases are mentioned in the panels}.
\label{fig:freq_effect}
\end{figure}

We have discussed earlier in section \ref{sec:disc_inertial} 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.(\ref{fig:freq_effect}) we have compared the jet solutions corresponding to two different frequencies of the disk oscillation. The injection parameters are $v_{\rm in}=0.076,\,\,\,\Theta_{\rm in}=0.064,\,\,\,r_{\rm in}=3.0$ (same as Fig. \ref{fig:Shockcolsn}). 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_{\rm s}$). The time at which the disk reaches same $l$ and $x_{\rm s}$, will be different for the two frequencies, and has been shown on each panel. Assuming an $M_{\rm B}=10\,M_{\odot}$ 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 $\sim 7 M_{\odot})$ itself shows the QPOs at various frequencies extending from $0.1-10 {\rm Hz}$ \citep{cdp08,rmmbo98}. Hence, as a representative case for two different frequencies of disk oscillations we have taken the values $100\,{\rm mHz}$ and $5\,{\rm 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.

\section{Discussion and Conclusion}
\label{sec:concl}
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$\lower.5ex\hbox{$\; \buildrel < \over \sim \;$}$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 $\gamma_{\rm T}\lower.5ex\hbox{$\; \buildrel > \over \sim \;$} 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 ultra-relativistic speeds ($\gamma_{\rm 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 \cite{kcm14} 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_{\rm in}$ (jet base), $v_{\rm in}$ is higher for higher luminosity, but difference in sound speed or $\Theta_{\rm 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_{\rm in}$ with a frequency similar to that of the oscillating $x_{\rm 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_{\rm 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\sim 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_{\rm s}$ or $\dot{m}_{\rm 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.

\section*{acknowledgments}
 We thank the anonymous referee for valuable comments which improved the quality of the paper.