We use version 4.6.2 of the adaptive-mesh refinement (AMR) hydrodynamical FLASH code [1]} in three dimension (3D). We turn off the radiative cooling since the regions we simulate are optically thick.
We include radiation pressure, electrons pressure, and ions pressure as we assume an adiabatic index of \(\gamma =5/3\)  and \(\mu = 0.5\) , as we assume a fully ionised pure hydrogen.
We set outflow conditions on all the boundary surfaces of the 3D grid. We use resolution with 7 refinement levels, i.e., the ratio of the length of the sides of the grid to the shortest cell in the grid is \(2^{9}\) . The finest cell in the grid has a side of \(1.05 AU\) . We simulate the whole space with a total size of the Cartesian numerical grid of \((534 AU)^3\) , i.e., \((L_x,L_y,L_z) = \pm 267 AU\) .
We simulate the entire space, rather than only one half of the symmetry plane \(z=0\) , to avoid numerical artifacts at the symmetry plane had we simulated only one half of the interaction region. We will present the flow structure in the entire computational grid to emphasize the bipolar structure that the jets shape and the flow in the equatorial plane.
