# A Numerical Study of Relativistic Oblique Shock Reflection

Prasanta Bera1 , Jonathan Granot1,2,3, Michael Rabinovich1,2, Paz
Beniamini1,2,3
1Astrophysics Research Center of the Open University (ARCO), The Open
University of Israel, P.O Box 808, Ra’anana 4353701, Israel
2Department of Natural Sciences, The Open University of Israel, P.O Box 808,
Ra’anana 4353701, Israel
3Department of Physics, The George Washington University, 725 21st Street NW,
Washington, DC 20052, USA<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Shocks are ubiquitous in astrophysical sources, many of which involve
relativistic bulk motions, leading to the formation of relativistic shocks.
Such relativistic shocks have so far been studied mainly in one dimension, for
simplicity, but the complex nature of the relevant astrophysical flows often
requires higher dimensional studies. Here we study the two-dimensional problem
of the reflection of a planer shock off of a wall for a general incidence
angle and a cold unshocked medium. We use primarily relativistic hydrodynamic
numerical simulations, and elaborately compare the results to an analytic
treatment. The simulations are performed both in the rest frame S of the
unshocked fluid, where the dimensionless proper speed of the singly shocked
fluid is $u_{1}=\Gamma_{1}\beta_{1}$ and the shock incidence angle is
$\alpha_{1}$, and in the rest frame S′ of the point P of intersection of the
incident shock and the wall for regular reflection (RR). Good agreement is
obtained between the simulations in these two frames and with the analytic
solution. The establishment of a steady flow in frame S′ is explored, along
with the transition between the strong and weak shock RR solutions. The
transition line between RR and Mach reflection (MR) is studied numerically in
the $u_{1}$ \- $\alpha_{1}$ plane and found to coincide with the analytic
detachment/sonic line. The flow properties along the sonic line are
investigated in detail focusing on how they vary between the Newtonian and
relativistic limits.

###### keywords:

shock waves – relativistic processes – methods: numerical – hydrodynamics

††pagerange: A Numerical Study of Relativistic Oblique Shock
Reflection–References††pubyear: 0000

##  1 Introduction

A steady single-phase subsonic inviscid flow maintains a smooth variation over
different locations, excluding interfaces or boundaries. However, supersonic
fluid velocities (with relative speeds between different parts of the fluid
that exceed the sound speed) may form a discontinuity in matter density,
pressure and velocity, which is termed a shock.

The location of the discontinuity (i.e., the shock) generally moves in space.
The fluid crosses the shock from the upstream region to the downstream region
and in the process, its density, pressure and specific entropy increase (see
e.g., Landau & Lifshitz, 1987; Thorne & Blandford, 2017). Rankine–Hugoniot
conditions specify the relationship between the fluid variables across the
discontinuity (Rankine, 1870; Hugoniot, 1887). In the rest frame of the
upstream fluid, the shock front moves supersonically, and the downstream
shocked fluid carries nonzero momentum, kinetic energy and thermal energy.
Shocks are very abundant in terrestrial and astrophysical fluids in supersonic
motion.

In terrestrial phenomenon, the motion of a fluid (e.g. air, water) can attain
a speed larger than the respective sound speed in the medium and this can form
a shock. The head-on interaction of this discontinuity with a rigid wall,
produces a reflected shock with a subsonic downstream region. In the case of
an oblique incidence, the strength of the incident shock and the angle of the
incidence determine the characteristics of the reflection. When the reflected
shock and the incident shock intersect at a reflection point P on the wall, it
is said to be regular reflection (RR). Otherwise, it is considered to be
irregular reflection (IR), the most common configuration of which is called a
Mach reflection (MR). In the case of MR, there exists a triple point ahead of
the wall, where three lines intersect: the incident shock, the reflected
shock, and a Mach stem (Von Neumann, 1963; Courant & Friedrichs, 1948;
Chester, 1954; Hornung, 1986; Olim & Dewey, 1992; De Rosa et al., 1992; Tabak
& Rosales, 1994). For large values of the shock incidence angle (defined as
the angle between the shock and the wall), only IR/MR is possible, whereas for
small incidence angle values only RR is possible.

Shock reflection of non-relativistic oblique shocks was investigated in
different experimental setups (Heilig, 1969; Itoh et al., 1981; Henderson &
Gray, 1981) and numerical studies (Mignone et al., 2007; Gvozdeva & Chulyunin,
2015; Wu et al., 2019). One of the main purposes of these studies was to
pursue the transition criteria from RR to MR, and vice versa. Some theoretical
criteria for this transition are known in the literature (see, e.g. Von
Neumann, 1963; Hornung et al., 1979; Ben-Dor, 1987). For our purposes in this
paper, the important criterion is the sonic criterion, which is discussed in
detail below.

In an astrophysical environment, a fluid element can achieve a speed close to
the speed of light, $c$, and form relativistic shocks, capable of generating
significant radiation (Blandford & McKee, 1976). The microscopic properties of
the fluid are affected by relativistic thermal particle motions as reflected
in the equation of state (Taub, 1948; Thorne, 1973). Relativistic bulk motions
have observational implications such as relativistic beaming effects (Rees,
1966; Gold, 1969). The strength (i.e. the Lorentz factor) of relativistic
shocks may be inferred from the modeling of astrophysical objects, such as
gamma-ray burst afterglows (Sari, 1997; Rees & Mészáros, 1998).

Shocks play an important role in various astrophysical scenarios, such as: i)
accretion by compact object (Salpeter, 1964) ii) free-falling accretion onto
the surface of a star, iii) interaction of stellar wind with the interstellar
medium, iv) high-velocity ejecta from explosive transients, e.g. a nova,
supernova or magnetar giant flare, v) in relativistic jets or outflows, such
as gamma-ray bursts (GRBs), micro-quasars, active galactic nuclei (AGN), tidal
disruption events, fast radio bursts or pulsar wind nebulae (PWNe), where
shocks can form either due to collisions between different parts of the
outflow (internal shocks) or due to its interaction with the ambient medium
(external shocks). These astrophysical sources form some regions with very
high internal energy density and the shocks accelerate both thermal and non-
thermal electrons that produce bright radiation. Therefore the shock dynamics
play a significant role in generating the observable radiation from many
astrophysical sources.

Such astrophysical shocks may experience reflection by an obstacle. Some
possible examples are: i) reflection of a supernova shock by the companion
star in a binary stellar system (Istomin & Soloviev, 2008), ii) reflection of
a GRB afterglow shock (Lamberts & Daigne, 2018), iii) reflection of shock
formed at the magnetosphere by the stellar surface of a neutron star or the
Sun, iv) reflection of a collimation shock at the jet-cocoon interface with a
cocoon in the cylindrical phase (Adamson & Nicholls, 1958; Norman et al.,
1982). To understand the underlying physics we can build a theoretical model
relating the flow dynamics to the expected observable emission signatures. We
follow a simplified approach to study the fluid dynamics of shock reflection
in relativistic and non-relativistic regimes. We consider a perfectly
reflecting wall as the reflector of the incident shocks.

In particular, in this work we numerically study the reflection of an incident
oblique shock having Newtonian up to relativistic speeds, at different
incidence angles. From direct relativistic hydrodynamic numerical simulation,
we identify the characteristics of the reflected shock and find the criteria
of RR. We compare our numerical results to analytic results derived in a
companion paper (Granot & Rabinovich, 2023, hereafter GR23)

Initially, in § 2, we describe the physical setup of the numerical
experiments. The underlying basic mathematical formulation is presented in §
3. Our results are presented in § 4 and the conclusions are discussed in § 5.

##  2 Physical setup

###  2.1 Lab frame S & steady-sate frame S′

Figure 1: Schematic diagram of our setup for the shock reflection problem for
RR, showing the location of the discontinuities. Left: In the lab frame S the
unshocked cold fluid (region 0) is at rest and a piston moving at speed
$v_{1}$ at an angle of $\alpha_{1}$ relative to a wall drives a shock ($s1$)
into it (the shock front moving at speed $v_{s1}$) creating a singly shocked
fluid region 1. The shock $s1$ hits the wall producing a reflected shock
($s2$) with a shock front moving at speed $v_{s2}$, and a doubly-shocked fluid
region 2, whose velocity $v_{2}$ is parallel to the wall. The point $P$ where
the two shocks intersect at the wall moves along the wall at a speed
$v_{p}=v_{s1}/\sin\alpha_{1}=v_{s2}/\sin\alpha_{2}$. Right: In the rest frame
S′ of point P the flow is steady, and the fluid velocity in regions 0 and 2 is
parallel to the wall. This rest frame exists only in the sub-luminal regime
where $v_{p}<c\Leftrightarrow u_{s1}<\tan\alpha_{1}$ ($u_{s1}$ being the
proper speed of shock s1).

The shock reflection is studied in two different frames of reference: i) the
lab-frame S, where the unshocked region 0 is at rest, and ii) the moving frame
S′, where the flow is steady for RR (Figure 1).

Initially, we set up the problem in the lab frame S. There are two shocks
labeled 1 (incident shock) and 2 (reflected shock), which divide the flow into
three regions, labeled 0, 1 and 2, corresponding to the number of times the
fluid in each region was shocked. The unshocked region 0 is adjacent to a
perfectly reflecting static wall and considered to be at rest in frame S
(velocity $v_{0}=0$) and cold (pressure
$p_{0}/\rho_{0}c^{2}\ll\min(1,u_{1}^{2})$, where $\rho_{0}$ is its proper
rest-mass density and $u_{1}=\Gamma_{1}\beta_{1}$ is the dimensionless proper
speed of region 1). In frame S, the incident shock (‘shock 1’) moves with a
velocity $v_{s1}$ along its normal and makes an angle $\alpha_{1}$ with
respect to the wall. It can be thought of as generated by a piston moving at
velocity $v_{1}<v_{s1}$ and driving a shock with a velocity $v_{s1}$ (see Fig.
1). The proper rest-mass density $\rho_{1}$ and pressure $p_{1}$ in region 1
are determined by the jump conditions of shock 1. The collision of the
incident shock 1 with the wall creates a reflected ‘shock 2’, with a velocity
$v_{s2}$ along its normal, and making an angle $\alpha_{2}$ with respect to
the wall. A post-shock region 2 forms between the wall and shock 2, with
proper rest-mass density $\rho_{2}$, pressure $p_{2}$ and velocity $v_{2}$.

As the incident shock 1 is oblique, $\alpha_{1}>0$, it intersects the wall at
a point P, which moves along the wall at a velocity $\boldsymbol{v}_{p}$,whose
magnitude is given by

$v_{p}=\frac{v_{s1}}{\sin\alpha_{1}}=\frac{v_{s2}}{\sin\alpha_{2}}\ .$ (1)

In the case of a semi-infinite steady oblique incident shock, the incident and
reflected shocks move in a self-similar pattern with respect to point P. In
the sub-luminal regime that corresponds to $\beta_{p}=v_{p}/c<1$, one can
transform (through a boost at $-\boldsymbol{v}_{p}$) to a rest frame S’, in
which for RR, point P is at rest and the flow is steady. In frame S′, region 0
moves with a velocity $\boldsymbol{v}^{\prime}_{0}=-\boldsymbol{v}_{p}$.
Similarly, the velocity of region 1 and the incident shock are Lorentz-boosted
by $-\boldsymbol{v}_{p}$ from the S-frame values. The proper rest-mass
densities ($\rho_{0},\rho_{1}$) and pressures ($p_{0},p_{1}$) are invariant.
The angles formed by the incident and reflected shocks with the wall in S′ are
also Lorentz boosted, i.e.,

$\tan\alpha^{\prime}_{i}=\frac{\tan\alpha_{i}}{\Gamma_{p}}\quad\quad(i=1,\,2)\
,$ (2)

where $\Gamma_{p}=(1-\beta_{p}^{2})^{-1/2}$.

###  2.2 The General Structure of the Parameter Space

Figure 2: The different regions and critical lines in the $u_{1}$ –
$\alpha_{1}$ parameter space, shown in terms of $\log_{10}(u_{1})$ versus
$\log_{10}(\tan\alpha_{1})$. The luminal line (in black;
$\beta_{p}=1\Leftrightarrow u_{s1}=\tan\alpha_{1}$) separates the super-
luminal region (cyan shading) and the sub-luminal attachment region (green
shading), which is in turn separated from the detached region (in white, where
there is no regular reflection – RR) by the detachment line (in dashed red),
which almost coincides with the sonic line for the weak RR solution (in dashed
blue; $\beta^{\prime}_{2,w}=\beta_{c_{s},2,w}$, see GR23 for details).

We study oblique reflected shocks with different incidence angles,
$\alpha_{1}$, and different proper velocities of the incident fluid,
$u_{1}=\Gamma_{1}\beta_{1}$. Figure 2 shows the analytic expectation (as
derived in GR23) for the different regions in the $u_{1}$ – $\alpha_{1}$
parameter space, and the critical lines that separate between them. This is
displayed by showing $\log_{10}(u_{1})$ in the $y$-axis versus
$\log_{10}(\tan\alpha_{1})$ in the $x$-axis. The luminal line (in black;
defined by the condition $v_{p}=c$ or, equivalently, $u_{s1}=\tan\alpha_{1}$)
separates the super-luminal region (in cyan shading) and the sub-luminal
regions. The sonic line for the weak shock RR solution (in dashed blue;
defined by $\beta^{\prime}_{2,w}=\beta_{c_{s},2,w}$ where the subscript ‘w’
stands for the weak shock RR solution) is found (GR23) to almost coincide with
the detachment line (in dashed red), which bounds the region with RR
solutions. We shall therefore not make the distinction between them here, and
refer mainly to the sonic line. The sonic line always lies in the sub-luminal
region111This is since it corresponds to
$\beta_{p}=(\beta_{2,w}+\beta_{c_{s},2,w})/(1+\beta_{2,w}\beta_{c_{s},2,w})<1$,
i.e. the sonic condition implies that $v_{p}$ is equal to the lab-frame speed
of a sound wave moving in region 2 parallel to the wall, which must therefore
be less than $c$. and separates between the sub-sonic (or detachment) region
(in white), where there is no RR solution (but instead, there is IR – a more
complicated type of shock reflection, such as MR), and the super-sonic (or
attachment) regions. The region between these two critical lines – the sub-
luminal super-sonic (or attachment) region is marked in green shading.

In the following, we find numerically that, as expected analytically, the
sonic line bounds the region of RR. There can in principle also be a dual
region where both MR and RR are possible for the same $(u_{1},\alpha_{1})$
values. Such a dual region borders the sonic line on the super-sonic side but
is located well within the sub-luminal region. The fact that we do not find
such a dual region might be since the RR weak shock solution is a more stable
attractor solution, such that the MR solution is not found in the simulations,
similar to the RR strong shock solution that is discussed in § 4.1.1.

Therefore, the sonic line is of particular physical importance. While it was
extensively studied in the Newtonian regime, it was not studied before in the
relativistic regime. We study it here in detail, stressing the differences
between the Newtonian and relativistic regimes, and how the system transitions
between these two limits.

Figure 3: Scatter marks represent the parameter space coverage of the
numerical calculations presented in section 4. The blue, magenta and red
circles correspond to the snapshots of RR in S and S′ and MR in S′ (section
4.1) respectively. For a small inclination angle, the S-frame captures RR
effectively. The frame S′ is applicable for the incident angle higher than the
luminal boundary.

We study shock reflection in the ($u_{1}-\alpha_{1}$) parameter space. We use
the above-mentioned reference frames S and S′. Figure 3 shows the points for
which we performed special relativistic hydrodynamic numerical simulations
(section 3) to obtain the outcome of the shock reflection by a wall.

##  3 Numerical Method

The conservation equations for total mass, momentum and energy in the special
theory of relativity may be written as:

$\displaystyle\partial_{\mu}(\rho u^{\mu})$
$\displaystyle=\frac{\partial(\rho\Gamma)}{\partial
t}+\nabla\cdot(\rho\textbf{{u}})=0\ ,$ (3)
$\displaystyle\partial_{\nu}T^{i\nu}$
$\displaystyle=\frac{\partial(w\Gamma\textbf{{u}})}{\partial
t}+\nabla\cdot(w\textbf{{u}}\textbf{{u}}+p\boldsymbol{I})=0\ ,$ (4)
$\displaystyle\partial_{\nu}T^{0\nu}$
$\displaystyle=\frac{\partial(w\Gamma^{2}-p)}{\partial
t}+\nabla\cdot(w\Gamma\textbf{{u}})=0\ ,$ (5)

where $T^{\mu\nu}=wu^{\mu}u^{\nu}+p\eta^{\mu\nu}$ is the stress-energy tensor,
$\eta^{\mu\nu}$ is the Minkowski metric, $u^{\mu}$ is the 4-velocity,
$\mathbf{u}=\Gamma\boldsymbol{\beta}=u\hat{\boldsymbol{u}}$ is the proper
velocity of the fluid,
$\boldsymbol{\beta}=\boldsymbol{v}/c=\beta\hat{\boldsymbol{\beta}}$,
$\Gamma=(1-\beta^{2})^{-1/2}$ is the Lorentz factor, $\rho$ is the proper rest
mass density, $p$ is the pressure, $w=e+p=\rho c^{2}+e_{\rm int}+p$ is the
proper enthalpy density, and $e$ ($e_{\rm int}$) is the proper (internal)
energy density. Here $\frac{\partial}{\partial t}$, $\nabla$ and
$\boldsymbol{I}$ are the time derivative, spatial derivative and the unit
$3\times 3$ matrix, respectively.

In the presence of a 1D shock, the fluid variables on its two sides (0,1: pre-
& post-shock regions) satisfy the following (Rankine–Hugoniot) jump conditions
conditions,

$\displaystyle\rho_{0}\Gamma_{0,s1}\beta_{0,s1}$
$\displaystyle=\rho_{1}\Gamma_{1,s1}\beta_{1,s1}\ ,$ (6) $\displaystyle
w_{0}\Gamma_{0,s1}^{2}\beta_{0,s1}^{2}+p_{0}$
$\displaystyle=w_{1}\Gamma_{1,s1}^{2}\beta_{1,s1}^{2}+p_{1}\ ,$ (7)
$\displaystyle w_{0}\Gamma_{0,s1}^{2}\beta_{0,s1}$
$\displaystyle=w_{1}\Gamma_{1,s1}^{2}\beta_{1,s1}\ ,$ (8)

where quantities relating to the upstream (pre-shock) region and the
downstream (post-shock) region are denoted by subscripts 0 and 1,
respectively. These jump conditions may be obtained by equating the fluxes of
matter, momentum and energy on the two sides of the shock, in the frame where
the shock front is at rest and the fluid velocities are normal to it. The
velocities $\beta_{0,s1}$ and $\beta_{1,s1}$ are those of regions 0 and 1 in
the shock 1 rest-frame, while $\Gamma_{i,s1}=(1-\beta_{i,s1}^{2})^{-1/2}$ are
the corresponding Lorentz factors. The pressure, density and normal component
of velocity are discontinuous across the shock.

To solve the above set of fluid equations and to obtain the values of
downstream fluid for the given upstream values we need to provide the equation
of state (EoS). To capture the relativistic and non-relativistic regimes we
consider the following equation (Mignone & McKinney, 2007):

$\displaystyle\left(h-\Theta\right)\left(h-4\Theta\right)=1\ ,$ (9)

where $\Theta=p/\rho c^{2}$ and the enthalpy per unit rest energy $h$ and the
effective adiabatic index $\hat{\gamma}$ are give by

$\displaystyle h$ $\displaystyle=$ $\displaystyle
1+\frac{\hat{\gamma}\Theta}{\hat{\gamma}-1}=\frac{5}{2}\Theta+\sqrt{1+\frac{9}{4}\Theta^{2}}\
,$ (10) $\displaystyle\hat{\gamma}$ $\displaystyle=$
$\displaystyle\frac{\frac{\partial h}{\partial\Theta}}{\frac{\partial
h}{\partial\Theta}-1}=\frac{1}{6}\left(8-3\Theta+\sqrt{4+9\Theta^{2}}\right)\
.\quad\quad$ (11)

This EoS satisfies the Taub (1948) inequality of relativistic matter. The
corresponding dimensionless sound speed, $\beta_{c_{s}}=c_{s}/c$, is given by
(Ryu et al., 2006)

$\displaystyle\beta_{c_{s}}^{2}=\frac{\partial p}{\partial
e}=\frac{\Theta}{h}\frac{\frac{\partial h}{\partial\Theta}}{\frac{\partial
h}{\partial\Theta}-1}=\frac{3\Theta^{2}+5\Theta\sqrt{\Theta^{2}+4/9}}{12\Theta^{2}+2+12\Theta\sqrt{\Theta^{2}+4/9}}\
.$ (12)

Here we aim to find the values of the downstream quantities from the direct
hydrodynamic simulations. To this end we study the shock reflection in frames
S & S’ described above.

###  3.1 Numerical setup in frame S

In order to study the shock reflection problem described in § 2 in frame S, we
set up the incident shock s1 in this frame by prescribing regions 0 and 1 in
the computation domain. We then numerically solve the time evolution of the
computation domain, identifying the different regions and the critical lines
that separate them. In particular, we track the formation of the reflected
shock s2 and the doubly shocked region 2 as it is described in § 2.1. To
reduce the artifacts from the numerical scheme, we avoid ultra-low values of
pressure in region 0 and choose a moderately low value of
$\Theta_{0}=p_{0}/\rho_{0}c^{2}\sim 10^{-9}$ to represent a cold medium of
region 0. The fluid is at rest in region 0 ($\boldsymbol{v}_{0}=0$) while the
velocity of region 1 is $\boldsymbol{v}_{1}$. The pressure ($p_{1}$) and
proper rest-mass density ($\rho_{1}$) in region 1, as well as the velocity of
shock s1 along its normal in frame S ($\boldsymbol{v}_{s1}$) are obtained by
solving the shock jump conditions (equations (6)-(8)).

We use the pluto code (Mignone et al., 2007) to solve the hydrodynamic
equations (3)-(5) in a fixed linear spaced grid. We use the piece-wise
parabolic reconstruction scheme, with the second-order Runge-Kutta time
integration and HLLC Reimann solver. We choose the initial shock location
along the diagonal of the computational domain connecting the top-left and the
bottom-right. We consider a few hundreds to thousands of grid points in each
side of the two-dimensional computational domain maintaining near near-equal
aspect ratio of grid-spacing. The presence of the wall at the r.h.s. boundary
is obtained by implementing reflecting boundary conditions. The region 1
inflow boundary conditions are implemented at the left and at the bottom of
the computation domain. The top of the computation domain is maintained with a
free outflow condition.

In the lab frame S, region 1 increases its area as the point of contact P
moves along the wall. The shape of the post-shock region i.e. the angle
$\alpha_{2}$ is not known a priori. To start the simulation we use the inflow
condition from the left and lower boundaries. As the doubly-shocked region 2
develops and it interferes with this fixed inflow condition along the lower
boundary. To prevent the impact of the boundary conditions on the results we
focus our analysis on a region sufficiently close to the point of contact P,
such that it is not affected by the lower boundary condition.

###  3.2 Numerical setup in frame S′

By construction, the shock s1 remains static in frame S′ and the region 0 is
confined between the shock s1 and the wall which makes an angle
$\alpha_{1}^{\prime}$ at P, given by
$\tan\alpha_{1}^{\prime}=\Gamma_{p}^{-1}\tan\alpha_{1}$. Region 0 has proper
rest mass density $\rho_{0}$, pressure $p_{0}$ and velocity
$\boldsymbol{v}^{\prime}_{0}=-\boldsymbol{v}_{p}$. Region 1 has proper rest
mass density $\rho_{1}$, pressure $p_{1}$ and the velocity is given by a
Lorentz boost by $-\boldsymbol{v}_{p}$ from the frame S value
$\textbf{{v}}_{1}=\beta_{1}c(-\cos\alpha_{1},\,\sin\alpha_{1})$,

$\boldsymbol{v}^{\prime}_{1}=\frac{\left[-v_{1}\cos\alpha_{1},\,\Gamma_{p}(v_{1}\sin\alpha_{1}-v_{p})\right]}{\Gamma_{p}(1-\beta_{p}\beta_{1}\sin\alpha_{1})}\
.$ (13)

In the pluto setup we start the numerical simulation with the fluids in
regions 0 and 1. The top and left edges of the computational domain maintain
inflow boundary conditions of regions 0 and 1, respectively, while the bottom
edge maintains a free outflow boundary condition. As the incident shock s1
impacts the wall, it forms the reflected shock s2 and the doubly-shocked
region 2 develops.

We also test the dynamical stability of the RR strong shock solution, by
adding the corresponding algebraic solution for region 2 to the initial
conditions of the simulation (as this solution is unstable and does not
otherwise develop naturally in the numerical simulation). In this case the
point of transition between the inflow boundary conditions of regions 0 and 1
is no longer at the top left corner, but is instead located at a fixed point
along the top edge of the computational domain.

We consider the evolution in frame S′ in the vicinity of point P, as much as
possible. The advantage of this frame is that the incident and reflected
shocks are static, and the flow is steady for RR. However, since
$\beta_{p}=\beta_{1s}/\sin\alpha_{1}$, for high incident shock speeds
$\beta_{1s}$ and/or small incidence angles $\alpha_{1}$, the velocity of point
P might rise above the speed of light, and in this super-luminal regime frame
S’ does not exist. We start a numerical evolution of fluid in region 1 and
gradually it develops the region 2. As the post-shock region (region 2)
develops, we select a region away from the boundary and find the location of
the discontinuity.

When performing simulations in frame S’ we inject fluid into region 0 and 1
with velocities $\boldsymbol{v}^{\prime}_{0}$ and
$\boldsymbol{v}^{\prime}_{1}$, from the upper and left boundaries of the
simulation box, respectively.

We have the freedom to choose arbitrary time units, $t_{\rm unit}$ in S and
$t^{\prime}_{\rm unit}$ in S′ (corresponding to length units $l_{\rm
unit}=ct_{\rm unit}$ and $l^{\prime}_{\rm unit}=ct^{\prime}_{\rm unit}$), to
design the frame for the direct numerical study. In frame S we measure the
simulation time in units of shock crossing time, $t_{p}=L_{y}/v_{p}$, where
$L_{y}$ is the simulation box size along the wall. In frame S′ we use as our
time unit the sound crossing time of the doubly shocked region 2,
$t^{\prime}_{sc}=L^{\prime}/c_{s,2}$, where $L^{\prime}$ is its length along
the wall and $c_{s,2}$ is the sound speed in region 2.

##  4 Results

###  4.1 The $u_{1}$ – $\alpha_{1}$ parameter space

Here we summarize the results obtained from our numerical simulations of shock
reflection for different proper speeds $u_{1}$ of the singly shocked region 1
and different incidence angles $\alpha_{1}$. For a given $u_{1}$, RR is
expected for small enough values $\alpha_{1}$ (see Figure 2).

####  4.1.1 Regular Reflection (RR)

Figure 4: Snapshots from a numerical simulation of shock reflection in the lab
frame S for RR, at different times: (a) $t/t_{p}=0.0$, (b) $t/t_{p}=0.57$, (c)
$t/t_{p}=0.86$, where $t_{p}=L_{y}/v_{p}$ is the box crossing time of point P
(using equal aspect ratio). The computation domain of lengths ratio
$L_{y}:L_{x}=3.23:1$, where the reflecting wall is along its right boundary
while the incident shock s1 is initially along its top-left to bottom-right
diagonal. The unit of length is arbitrary. The incident shock s1 leaves the
computation domain at $t=t_{p}$. The color scale indicates the fluid’s proper
rest-mass density while the red arrows show its velocity vector. This
simulation is initialized with $\alpha_{1}=0.3$ and $u_{1}=1$, where the
latter implies $u_{1s}=1.37$ and $\rho_{1}/\rho_{0}=4\Gamma_{1}=4\sqrt{2}$ for
$p_{0}\ll\rho_{0}c^{2}u_{1}^{2}$. Panel (a) shows the initial conditions.
Other snapshots (panels (b), (c)) indicate the gradual growth of the doubly
shocked region 2 with a higher density. The reflected shock forms an angle
$\alpha_{2}=0.137$ with the wall. Region 2 remains uniform unless it is
affected by the lower boundary condition. Figure 5: Snapshots from a numerical
simulation of shock reflection in the frame S’ for RR, at different times: (a)
$t^{\prime}/t_{sc}^{\prime}=0$, (b) $t^{\prime}/t_{sc}^{\prime}=0.72$, (c)
$t^{\prime}/t_{sc}^{\prime}=1.43$ (where $t_{sc}^{\prime}$ is defined in §
2.1). The height-to-width ratio of the computation domain is
$L^{\prime}_{y}/L^{\prime}_{x}=4.13$. The doubly-shocked region 2 develops
between the wall and the reflected shock s2, which forms an angle
$\alpha^{\prime}_{2}$ (with $\tan\alpha_{2}^{\prime}=0.33$) relative to the
wall. Figure 7 shows zoomed-in snapshots (of the region below the dashed white
line in panel (b)), more densely sampled in time to better illustrate the
formation of region 2.

We have captured the time evolution of incident shock s1 and the development
of reflected shock s2 by performing numerical simulations in two different
rest frames. In the lab frame S, both the incident and the reflected shocks
move (left panel of Figure 1 and snapshots from the numerical study in Figure
4). The same incident and reflected shocks remain steady at the moving frame
$S^{\prime}$ (right panel of Figure 1 and snapshots from the numerical study
in Figure 5).

In frame S we start the numerical simulation of an incident shock s1 with
$u_{1}=1,\beta_{1}=1/\sqrt{2}\approx 0.7071,\beta_{1s}=4\sqrt{2}/7\approx
0.8081$, $\rho_{1}/\rho_{0}=4\sqrt{2}\approx 5.657$, $p_{1}=4/3$) and an
incidence angle $\alpha_{1}=0.3$ relative to the reflecting wall (at the
bottom right corner of Figure 4a). Therefore, initially there is no doubly-
shocked region 2. The intersection point P of the incident shock s1 and the
wall moves along the wall at a speed $v_{p}=\beta_{1s}/\sin\alpha_{1}$. We
consider a computation box (a $568\times 918$ grid in the $x$-$y$ plane) with
its vertical length ($L_{y}$) along the wall being 3.23 times larger than the
horizontal length ($L_{x}$). We measure the time in units of the shock
crossing time i.e. $t_{p}=L_{y}/v_{p}$. As time evolves, a high-density
doubly-shocked region 2 develops, between the wall and the reflected shock s2,
which makes an angle $\alpha_{2}=0.137$ relative to the wall (Figure 4b & 4c).
The similarity between the Figures 4b & 4c indicates the self-similar nature
of the flow with respect to the point P. The fluid in region 2 moves along the
wall, relatively slowly, at a proper speed $u_{2}=0.247$ for
$(u_{1},\alpha_{1})=(1,0.3)$. The inflow boundary condition at the lower
boundary is unphysical within region 2, and its effects become more
significant for a higher value of $\alpha_{1}$.

In frame S’ we initialized the numerical simulations as described in § 3.2.
Figure 5 shows snapshots from such a simulation with
$(u_{1},\alpha_{1})=(1,1)$ (or $\tan\alpha_{1}=1.557$), such that
($\rho_{0}=1$, $p_{0}=10^{-9}\rho_{0}c^{2}$) imply
($\beta_{0y}^{\prime}=-\beta_{p}=-0.907$, $\rho_{1}=5.657$, $p_{1}=1.333$,
$\tan\alpha^{\prime}_{1}=0.434$). We start the numerical simulation with
regions 0 and 1 (Figure 5a) in the computational box (of $252\times 1038$ grid
points) with its vertical length along the wall ($L^{\prime}_{y}$) being 4.13
times the horizontal length ($L^{\prime}_{x}$). The doubly-shocked region 2
develops with time (Figure 5b & 5c) as the incident shock s1 remain static in
this frame. The reflected shock s2 settles down at an angle
$\alpha_{2}^{\prime}$ (with $\tan\alpha_{2}^{\prime}=0.33$) relative to the
wall. The matter in region 2 ($\rho_{2}/\rho_{0}=16.171$,
$p_{2}/\rho_{0}c^{2}=7.039$) moves along the wall with a proper speed
$u^{\prime}_{2}=1.207$.

Figure 6 shows that the doubly-shocked region 2 forms and settles down over a
timescale close to its sound crossing time,
$t^{\prime}_{sc}=L^{\prime}/c_{s,2}$, where $L^{\prime}$ is its length along
the wall. We identify region 2 through its higher density relative to region
1. Region 2 contains proper-density fluctuations of a few percent ($\lesssim
5\%$) relative to the mean value. The non-uniformity in region 2 is due to the
gradual transition at the boundary and the fluctuations.

Figure 7 shows the snapshots displaying the density and velocity of the
doubly-shocked region 2 at different times. Frame S′ is suitable to study the
shock interaction for high enough $\alpha_{1}$ values, corresponding to the
sub-luminal region.

For RR, our simulations (in frame S′) that initially did not contain the
doubly-shocked region 2, evolved to and settled at the ‘weak’ shock RR
solution. From the algebraic solution, one may also obtain a ‘strong’ shock RR
solution (in the super-sonic sub-luminal region – see GR23), corresponding to
higher values of $\alpha^{\prime}_{2}$ (and therefore $\alpha_{2}$),
$\rho_{2}$ and $p_{2}$. To study the stability of this strong shock solution
we numerically evolve the fluid variables starting from an initial
configuration that also includes region 2 with properties corresponding to the
algebraic strong shock solution (see Table 1, Figures 8 and 9).

Frame S | Frame S′
---|---
$\rho_{0}$ | 1 | $\rho_{0}$ | 1
$\beta_{0y}$ | 0 | $\beta_{0y}^{\prime}$ | -$v_{p}$=-0.9555363
$\alpha_{1}$ | 1.0079245 | $\alpha_{1}^{\prime}$ | 0.43718
$u_{1}$ | 1 | $u_{1}^{\prime}$ | 1.795769
$\rho_{1}$ | 5.6568 | $\rho_{1}$ | 5.6568
$p_{1}$ | 1.33333 | $p_{1}$ | 1.33333
| Weak | Strong |  | Weak | Strong
$\alpha_{2}$ | 0.8942 | 1.40331 | $\alpha_{2}^{\prime}$ | 0.35191 | 1.05019
$\beta_{2}$ | 0.737 | 0.9033 | $\beta_{2}^{\prime}$ | -0.73887 | -0.38151
$\frac{p_{2}}{\rho_{0}c^{2}}$ | 7.127 | 18.8817 | $\frac{p_{2}}{\rho_{0}c^{2}}$ | 7.127 | 18.8817
$\rho_{2}/\rho_{0}$ | 16.3425 | 27.6949 | $\rho_{2}/\rho_{0}$ | 16.3425 | 27.6949
Table 1: Weak and strong solutions of a shock collision. Figure 6: Average
density ($\rho_{2}$), pressure ($p_{2}$) and velocity component along the wall
($\beta_{2y}$) of the doubly-shocked region 2, as it forms and settles down
(to $\rho_{\rm 2,set}$, $p_{\rm 2,set}$, and $\beta_{2y{\rm,set}}$,
respectively) over about a sound crossing time ($t_{sc}^{\prime}$). The
vertical dotted lines indicate the time stamps of the snapshots shown in
Figure 7. Figure 7: The gradual development of the doubly-shocked region 2 is
shown in this sequence of snapshots in frame S’, depicting the lower half of
the computational domain from Figure 5 at the times indicated by the vertical
dotted line in Figure 6.

Figure 8 shows the resulting evolution of the mean density and pressure of the
doubly-shocked region, in terms of its fractional deviation from the weak and
strong shock solutions. The system quickly transitions from the algebraic
strong shock RR solution to its numerical counterpart, in which the density
and pressure in region 2 differ by $\sim\\!1.5\\!-\\!2\%$. The system then
starts to linearly deviate from this solution with a growth rate of about
$0.5~{}{t_{sc}^{\prime}}^{-1}$ or $e$-folding time about $2t_{sc}^{\prime}$
(two sound crossing times). Subsequently, the transition between the strong
and weak shock solutions enters a non-linear phase. Finally, the weak shock RR
solution is approached at about 10 sound crossing time ($t_{sc}^{\prime}$),
and the deviation from this solution appears to decrease exponentially with
time.

Figure 9 shows snapshots from the corresponding simulation (in frame S’) in
Figure 9 displaying the fluid variables at different times (indicated by the
dashed vertical lines in Figure 8). Panels b & $\bar{b}$ show a small (linear
order) change in density, pressure and shape of the doubly-shocked region 2
from its initial state (which is shown in panels a & $\bar{\rm a}$). Panels c
& $\bar{c}$ indicate the significant (non-linear) changes with the transient
appearance of a new (third) shock and a contact discontinuity, which bound a
triply-shocked region at the bottom-right corner of the snapshot. Panels d &
$\tilde{d}$ show the system reached the weak shock solution and it stays
there.

Figure 8: The evolution of the mean proper rest-mass density ($\rho_{2}$) and
pressure ($p_{2}$) in the doubly-shocked region 2, shown in terms of their
fractional deviations from their values for the ‘strong’ ($\rho_{\rm strong}$,
$p_{\rm strong}$) and ‘weak’ ($\rho_{\rm weak}$, $p_{\rm weak}$) shock RR
solutions. The simulation starts from the algebraic ‘strong’ shock solution
and the system moves to the ‘weak’ shock solution within several sound-
crossing times of region 2 ($t_{\rm sc}^{\prime}$). The vertical dashed lines
denote the times of the snapshots shown in Figure 9. Figure 9: Snapshots of
proper rest-mass density (top panels a-d) and pressure (bottom panels
$\bar{\rm a}$-$\bar{\rm d}$) from a simulation in frame S’ starting from the
algebraic strong shock RR solution (panels a, $\bar{\rm a}$). The red arrows
are velocity vectors, whose size indicates the fluid speed at their starting
point. This time sequence captures the evolution between the strong and weak
shock solutions, in the initial linear phase (panels b, $\bar{\rm b}$) and
subsequent non-linear phase (panels c, $\bar{\rm c}$). The system finally
settles in the weak shock solution (panels d, $\bar{\rm d}$).

####  4.1.2 Irregular Reflection (IR) or Mach Reflection (MR)

In the sub-sonic region, if there was RR then the dense, high-pressure fluid
in the doubly-shocked region 2 would be in causal contact with point P, and
cause it to detach from the wall, thereby leading to IR. For this reason there
is no RR in the sub-sonic region, and instead only IR types of shock
reflection such as MR. As mentioned in § 2.2, there may be a dual region
within the super-sonic region where both RR and IR/MR are possible, but the
weak shock RR solution appears to be the most stable attractor solution that
generically appears in our simulations.

Therefore, we generally expect the formation of IR/MR in our simulations in
the sub-sonic regime. For a given $u_{1}$, this corresponds to sufficiently
large incidence angles $\alpha_{1}$. For such incidence angles, the post-shock
region develops multiple zones separated by discontinuities.

Figure 10 shows an example of a simulation for such a case, where MR develops.
For this numerical simulation we considered an incident shock characterized by
$(\alpha_{1},\,u_{1})=(1.1,\,1.0)$, with an unshocked region 0 of
($\rho_{0}=1$, $p_{0}/\rho_{0}c^{2}=10^{-9}$, $u_{0}=0$). To calculate the
fluid variables in frame S′ we consider the corresponding boost of
$\beta_{p}=0.9067722$, which implies an S’ frame incidence angle of
$\tan\alpha_{1}^{\prime}=0.8283838$. We evolve the fluid maintaining the
boundary conditions mentioned in section 3.2.

Figure 10 shows snapshots from this simulation. The meeting point of the
incident and the reflected shocks, P, detaches from the reflecting wall
(panels b, c) and a Mach stem develops behind which there is singly shocked
fluid at pressure equilibrium across a contact discontinuity with doubly-
shocked fluid behind the reflected shock s2. One can clearly see the
development of Kelvin–Helmholtz instability along this contact discontinuity
due to the velocity shear (discontinuous parallel velocity component) across
it. The reflected shock assumes an irregular non-triangular shape. The Mach
stem slowly moves upward in frame S’ where this simulation is performed, at a
constant speed, such that its length (or the distance of point P from the
wall) increases linearly with time. In this paper, we do not explore in detail
the characteristics of this IR/MR and instead leave this for a future work.

Figure 10: Sanpshots from a simulation with $(\alpha_{1},\,u_{1})=(1.1,\,1.0)$
in which Mach reflection (MR) develops, at times
$t^{\prime}/t_{sc}^{\prime}=0$, $t^{\prime}/t_{sc}^{\prime}=0.2$, and
$t^{\prime}/t_{sc}^{\prime}=0.4$.Each panel shows a colormap of the proper
rest-mass density and red velocity vectors. The side ratio of the computation
domain is $L^{\prime}_{y}/L^{\prime}_{x}=2.46$.

####  4.1.3 consistency of numerical results in frames S and S’

Figure 11: Snapshots (proper rest-mass density colormap and red velocity
vectors) from simulations of RR for $(\alpha_{1},\,u_{1})=(0.465,\,0.316)$
performed in: (a) the lab frame S ($\tan\alpha_{2}=0.281$) at $t/t_{p}=0.65$,
and (b) the moving frame S′ ($\tan\alpha_{1}^{\prime}=0.246$,
$\tan\alpha_{2}^{\prime}=0.139$) at $t^{\prime}/t^{\prime}_{sc}=1.91$

.

Here we show the consistency of the numerical results obtained through
relativistic hydrodynamic simulations performed in rest frames S and S′.
Figure 11 shows the snapshots from frames S and S′ for the same physical shock
reflection case of $(\alpha_{1},\,u_{1})=(0.465,\,0.316)$. In the lab frame S
the frame S′ moves with velocity $v_{p}=0.87c$ upward along the wall, such
that the incidence angle $\tan\alpha_{1}=0.501$ in S transforms to
$\tan\alpha_{1}^{\prime}=0.246$ in S′. From the numerical evolution studies,
we obtain the development of the doubly-shocked region 2. In frame S we derive
a value of $\tan\alpha_{2}=0.281$ for the angle of the reflected shock s2
relative to the wall, and in frame S′ we derive a corresponding value of
$\tan\alpha^{\prime}_{2}=0.139$. Using a Lorentz boost from frame S′ to S we
obtain $\tan\alpha_{2}=0.284$ from the evolution in frame S′. The algebraic
solution corresponding to the same input parameters gives
$\tan\alpha_{2}=0.282355$. Therefore, the angle of reflection is consistent in
both frames S and S′ to within 1%, and both are consistent with the analytic
value.

####  4.1.4 consistency of the numerical and analytic results

Figure 12: Comparing the matter proper rest-mass density in the doubly-shocked
region 2 ($\rho_{2}$, normalized by $\rho_{0}$) from the analytic algebraic
solution (solid green lines) to our hydrodynamic simulation results, in the
lab frame S (blue x symbols) and in the moving frame S′ (red + symbols), for
$u_{1}=10,1,0.1$ (from top to bottom). The vertical dotted black line (dashed
cyan line) corresponds to the luminal (sonic) line. Frame S is well-suited for
low $\alpha_{1}$ values, while frame S′ is a favorable option near the sonic
line.  Figure 13: Similar to Fig. 12 but comparing the region 2 pressure.
Figure 14:  Similar to Fig. 12 but comparing $\tan\alpha_{2}/\tan\alpha_{1}$.

We compare the results of our numerical simulations with the exact algebraic
solution for the relatively simpler case of RR, for which such an analytic
solution can be obtained (GR23). The critical incidence angle along the sonic
line, $\alpha_{1,{\rm sonic}}(u_{1})$, below which RR is possible, increases
as the incident shock velocity $\beta_{s1}$ increases. A detailed comparison
is shown for the proper rest-mass density $\rho_{2}$ (Figure 12), and pressure
$p_{2}$ (Figure 13) of the doubly-shocked region 2, as well as the ratio of
the tangens of the angles relative to the wall of the reflected and incident
shock fronts, $\tan\alpha_{2}/\tan\alpha_{1}$ (Figure 14).

We note that lab frame S simulation is more accurate for a smaller incidence
angle $\alpha_{1}$. As the value of $\alpha_{1}$ increases the doubly-shocked
region 2 in frame S is more strongly affected by the imposed inflow lower
boundary condition. Hence, the numerical results deviate from the expected
value as $\alpha_{1}$ approaches the critical value for a RR (i.e. the sonic
line). Frame S′ is more suitable for simulations in this regime, and can
follow the shock reflection for longer times, as the flow becomes steady in S′
for RR.

###  4.2 The Sonic Line

####  4.2.1 The Significance of the Sonic Line

The sonic line corresponds to the condition

$\beta^{\prime}_{2,w}=\beta_{c_{s},2,w}\ \ \Longleftrightarrow\ \
\beta_{p}=\frac{\beta_{2,w}+\beta_{c_{s},2,w}}{1+\beta_{2,w}\beta_{c_{s},2,w}}\
.$ (14)

The first condition is that in the rest frame S′ where the flow is steady the
velocity in region 2 (of the doubly-shocked fluid) for the weak shock RR
solution, $\beta^{\prime}_{2,w}$, is equal to its sound speed,
$\beta_{c_{s},2,w}$. Once $\beta^{\prime}_{2,w}$ drops below
$\beta_{c_{s},2,w}$, i.e. in the subsonic regime, region 2 comes into causal
contact with point $P$, and can then potentially cause it to separate from the
wall resulting in MR. On the other hand, in the super-sonic regime
($\beta^{\prime}_{2}>\beta_{c_{s},2}$) region 2 is not in causal contact with
point $P$ (for an initial unperturbed weak shock RR solution) so it cannot
affect it and therefore point $P$ cannot separate from the wall and allow a
transition to IR/MR. This may potentially suppress a transition between the
weak shock RR solution and MR (in the dual region between the sonic line and
the mechanical equilibrium line; see e.g. GR23), and require a sufficiently
large perturbation for it to occur. The fact that the strong shock RR solution
is always subsonic ($\beta^{\prime}_{2,s}<\beta_{c_{s},2,s}$) may potentially
account for its instability, e.g. as found in § 4.1.1. Since the sonic
condition is that for causality between region 2 and point $P$, it can also be
expressed in the lab frame S such that the speed of a sound wave propagating
in region 2 along the wall towards point $P$,
$(\beta_{2}+\beta_{c_{s},2})/(1+\beta_{2}\beta_{c_{s},2})$, equals that of
point $P$, $\beta_{p}$.

####  4.2.2 The Flow Properties along the Sonic Line

The sonic line’s physical significance makes it interesting to study in detail
the flow properties along it. The analytic solution for the flow properties
along the sonic line is derived in Appendix A, along with analytic expressions
for all of the flow quantities in the Newtonian and relativistic limits.

Figure 15 shows the values of different hydrodynamic variables along the sonic
line ($\beta^{\prime}_{2,w}=\beta_{c_{s},2,w}$), conveniently parameterised
according to the value of $u_{1}$ along this line. These values are found by
numerically solving the set of algebraic equations for RR, together with the
sonic condition, in the frame S′ (Appendix A) or equivalently in the lab frame
S (as is done in GR23). The weak and strong shock RR solutions exactly
coincide at the detachment line, which almost coincides with the sonic line,
such that both solutions are extremely close along the sonic line. There is
excellent agreement with both the Newtonian and the relativistic limits that
are found analytically in Appendix A. Moreover, it can be seen from Fig. 15
that in the relativistic limit $u_{p}>u_{2}>u_{1}\gg 1$ while $u_{12}$ is of
order unity, such that the first (incident) shock $s1$ is ultra-relativistic,
while the second (reflected) shock $s2$ is mildly relativistic.

Figure 15: The values of different hydrodynamic variables along the sonic line
($\beta^{\prime}_{2}=\beta_{c_{s},2}$) are shown as a function of the value of
$u_{1}$ along this line (thick solid lines). The Newtonian and relativistic
limits from Equations (A) and (23), respectively, are indicated by thin dashed
lines (which fall on top of the thick solid lines). Figure 16: The asymptotic
Newtonian (right panel) and relativistic (left panel) flow configurations
along the sonic lines, for which the flow parameters are given in Eqs. (A) and
(23), respectively. Figure 17: The location of the transition from RR to MR
(red + symbols; $\alpha_{1,{\rm crit}}(u_{1})$) found from our numerical
simulations in frame S′, match the analytically calculated location of the
sonic line (solid cyan line; $\alpha_{1,{\rm sonic}}(u_{1})$). The dotted
black line indicates the luminal line, which is shown for reference.  Figure
18: The critical transition angle in frame S′, $\tan\alpha^{\prime}_{1,{\rm
crit}}(u_{1})$, from our numerical simulations (red + symbols), compared to
the analytically derived corresponding angle for the sonic line,
$\tan\alpha^{\prime}_{1,{\rm sonic}}(u_{1})$, whose asymptotic Newtonian and
relativistic limits are shown by horizontal black dotted lines.

Figure 16 shows the flow configurations for the asymptotic Newtonian and
relativistic limits along the sonic line, in the rest frame S′ where the flow
is steady and point $P$ is at rest. These limits are particularly simple and
can be fully solved analytically (Appendix A). In frame S′ the angles
$\alpha^{\prime}_{1}$, $\alpha^{\prime}_{2}$ and $\chi^{\prime}$ do not vary
that drastically between these two limits (see also Fig. 15).

Figure 17 shows the location of the transition from RR to MR (red + symbols).
For each $u_{1}$ values the incidence angle $\alpha_{1}$ is gradually
increased between different simulations until we find the critical value
$\alpha_{1,{\rm crit}}(u_{1})$ at which point P detaches from the wall,
signaling the transition from RR to MR. It is determined more accurately by
performing several iterations for each $u_{1}$ value. These numerical values
are in excellent agreement with the analytically calculated location of the
sonic line (solid cyan line in Fig. 17), $\alpha_{1,{\rm
crit}}(u_{1})=\alpha_{1,{\rm sonic}}(u_{1})$. These simulations clearly shows
that: (i) there is indeed no RR in the sub-sonic region, and (ii) in the
super-sonic region the simulations reach the weak shock RR solution and not
the strong shock RR solution (in the sub-luminal region) or MR (in the dual
region).

Figure 18 shows the critical transition angle in frame S′,
$\tan\alpha^{\prime}_{1,{\rm crit}}(u_{1})$, from our numerical simulations
(red + symbols), compared to the analytically derived corresponding angle for
the sonic line, $\tan\alpha^{\prime}_{1,{\rm sonic}}(u_{1})$, whose asymptotic
Newtonian and relativistic limits are shown by horizontal black dotted lines.
They agree to within about $1\%$. This critical angle in frame S′ decreases by
about $10\%$ between the Newtonian and relativistic limits (see also Fig. 15).

##  5 Conclusions

We have studied relativistic shock reflection, mainly numerically using two
dimensional relativistic hydrodynamic simulations, and with detailed
comparisons to analytic results. Our simulations were performed in two
different rest frames: the lab frame S where the cold unshocked fluid (region
0) is at rest, and the rest frame S′ where for RR the point P of intersection
of the incident and reflected shocks with the reflecting wall is at rest, and
the flow is in a steady state.

Our numerical simulations have validated the analytic results (derived mainly
in GR23, but also in Appendix A for the sonic line). We have also pointed out
the importance of using a suitable reference frame in different simulations of
shock interactions. In particular, for small incidence angles $\alpha_{1}$ the
lab frame S is more suitable, while the frame S′ that exists only in the sub-
luminal region is more suitable closer to the sonic line.

We have also studied the transition between RR and IR (namely MR). The
transition from RR to IR/MR maintains similar characteristics in the Newtonian
and relativistic regimes. In our numerical study, this transition occurred
exactly at the sonic/detachment line.

Moreover, in the super-sonic region the simulations always reached the weak
shock RR solution. We have found the alternative strong shock RR solution,
which exists in the super-sonic sub-luminal region, to be unstable. Moreover,
we numerically studied how it transitions to the weak shock RR solution, which
appears to be a stable attractor solution. While a dual region where either RR
or MR are possible should exist from analytic considerations (on the super-
sonic side of the sonic line but well within the sub-luminal region), it was
never reached in our simulations, suggesting that it is not an attractor
solution (and may also be unstable).

## acknowledgement

P. Bera is supported by the Israel Academy of Sciences and Humanities &
Council for Higher Education Excellence Fellowship Program for International
Postdoctoral Researchers. This research was funded in part by the ISF-NSFC
joint research program under grant no. 3296/19 (J.G.) and by the United
States-Israel Binational Science Foundation (BSF) under grant no. 2020747 (P.
Beniamini).

## Data availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

##  Appendix A The Solution along the Sonic Line

Here we derive an analytic solution for RR along the sonic line. The
conditions in region 1, for a cold region 0 ($p_{0}=0$,
$e_{0}=w_{0}=\rho_{0}c^{2}$ and $h_{0}=1$), can be conveniently calculated in
the lab frame S and are given by

$\displaystyle\rho_{1}=4\Gamma_{1}\rho_{0}\ ,\quad
p_{1}=\frac{4}{3}u_{1}^{2}\rho_{0}c^{2}\ ,\quad
e_{1}=4\Gamma_{1}^{2}\rho_{0}c^{2}\ ,\quad$ $\displaystyle e_{\rm
int,1}=\frac{4\Gamma_{1}u_{1}^{2}}{\Gamma_{1}+1}\rho_{0}c^{2}\ ,\quad
w_{1}=4\Gamma_{1}^{2}\left(1+\frac{\beta_{1}^{2}}{3}\right)\rho_{0}c^{2}\ ,$
(15) $\displaystyle\beta_{s1}=\frac{4\Gamma_{1}u_{1}}{4\Gamma_{1}^{2}-1}\
,\quad u_{1}=\frac{1}{2}\sqrt{u_{s1}^{2}-2+\sqrt{4+5u_{s1}^{2}+u_{s1}^{4}}}\
,$
$\displaystyle\beta_{1,s1}=\frac{\beta_{s1}-\beta_{1}}{1-\beta_{1}\beta_{s1}}=\frac{\beta_{1}}{3}\
,$

(GR23) where the last equation means that for our equation of state, in the
rest frame of the downstream fluid (region 1), the speed at which the shock is
receding is a third of the incoming upstream speed.

Since the sonic line is always in the sub-luminal regime, it can conveniently
be analyzed in frame S′. In this frame the flow is steady and there are two
oblique shocks: $s1$ and $s2$, at angles $\alpha^{\prime}_{1}$ and
$\alpha^{\prime}_{2}$, respectively, relative to the wall. The velocity of
region 1 in frame S′ can be expressed through

$\displaystyle\textbf{{u}}^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\left[-u_{1}\cos\alpha_{1},\;\Gamma_{p}\Gamma_{1}(\beta_{1}\sin\alpha_{1}-\beta_{p})\right]\
,\quad\quad$ $\displaystyle\Gamma^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\sqrt{1+u^{\prime\,2}_{1}}=\Gamma_{1}\Gamma_{p}(1-\beta_{1}\beta_{p}\sin\alpha_{1})$
(16) $\displaystyle=$
$\displaystyle\frac{3\,\Gamma_{1}\sin\alpha_{1}}{\sqrt{(4\Gamma_{1}^{2}-1)^{2}\sin^{2}\alpha_{1}-16\Gamma_{1}^{2}u_{1}^{2}}}\
,\quad\quad$ $\displaystyle\tan\chi^{\prime}$ $\displaystyle=$
$\displaystyle\frac{u^{\prime}_{1x}}{u^{\prime}_{1y}}=\frac{\beta_{1}\cos\alpha_{1}}{\Gamma_{p}(\beta_{p}-\beta_{1}\sin\alpha_{1})}\
,$ $\displaystyle\tan\alpha^{\prime}_{1}$ $\displaystyle=$
$\displaystyle\frac{\tan\alpha_{1}}{\Gamma_{p}}\
,\quad\quad\quad\tan\alpha^{\prime}_{2}=\frac{\tan\alpha_{2}}{\Gamma_{p}}\ .$

The remaining conditions are the oblique shock jump conditions in frame S′ and
the sonic condition, which read

$\displaystyle\rho_{1}u^{\prime}_{1}\sin\alpha^{\prime}_{+}$ $\displaystyle=$
$\displaystyle\rho_{2}u^{\prime}_{2}\sin\alpha^{\prime}_{2}\ ,$ $\displaystyle
w_{1}u^{\prime\,2}_{1}\sin^{2}\alpha^{\prime}_{+}+p_{1}$ $\displaystyle=$
$\displaystyle w_{2}u^{\prime\,2}_{2}\sin^{2}\alpha^{\prime}_{2}+p_{2}\ ,$
$\displaystyle w_{1}\Gamma^{\prime}_{1}u^{\prime}_{1}\sin\alpha^{\prime}_{+}$
$\displaystyle=$ $\displaystyle
w_{2}\Gamma^{\prime}_{2}u^{\prime}_{2}\sin\alpha^{\prime}_{2}\ ,$ (17)
$\displaystyle\beta^{\prime}_{1}\cos\alpha^{\prime}_{+}$ $\displaystyle=$
$\displaystyle\beta^{\prime}_{2}\cos\alpha^{\prime}_{2}\ ,$
$\displaystyle\beta^{\prime}_{2}$ $\displaystyle=$
$\displaystyle\beta_{c_{s},2}\ ,$

where we denote $\alpha^{\prime}_{+}=\chi^{\prime}+\alpha^{\prime}_{2}$ and

$\beta_{c_{s},2}^{2}=\frac{3\Theta_{2}^{2}+5\Theta_{2}\sqrt{\Theta_{2}^{2}+4/9}}{12\Theta_{2}^{2}+2+12\Theta_{2}\sqrt{\Theta_{2}^{2}+4/9}}\
,\quad\quad\Theta_{2}=\frac{p_{2}}{\rho_{2}c^{2}}\ .$ (18)

In the Newtonian limit ($\beta_{1}<\beta_{p}\ll 1$) the adiabatic index is
$\hat{\gamma}=5/3$ and the equations reduce to

$\displaystyle\frac{\rho_{1}}{\rho_{0}}=4\
,\quad\frac{p_{1}}{\rho_{0}}=\frac{4}{3}v_{1}^{2}\ ,\quad\frac{e_{\rm
int,1}}{\rho_{0}}=2v_{1}^{2}\ ,\quad\beta_{s1}=\frac{4}{3}\beta_{1}\ ,$
$\displaystyle\textbf{{v}}^{\prime}_{1}=v_{1}\left(-\cos\alpha_{1},\;\;\sin\alpha_{1}-\frac{4}{3\sin\alpha_{1}}\right)\
,$
$\displaystyle\tan\chi^{\prime}=\frac{v^{\prime}_{1x}}{v^{\prime}_{1y}}=\frac{3\cos\alpha_{1}\sin\alpha_{1}}{1+3\cos^{2}\alpha_{1}}\
,$ (19)
$\displaystyle\rho_{1}v^{\prime}_{1}\sin\alpha^{\prime}_{+}=\rho_{2}v^{\prime}_{2}\sin\alpha^{\prime}_{2}\
,$
$\displaystyle\rho_{1}v^{\prime\,2}_{1}\sin^{2}\alpha^{\prime}_{+}+p_{1}=\rho_{2}v^{\prime\,2}_{2}\sin^{2}\alpha^{\prime}_{2}+p_{2}\
,$ $\displaystyle
5\frac{p_{1}}{\rho_{1}}+v^{\prime\,2}_{1}\sin^{2}\alpha^{\prime}_{+}=5\frac{p_{2}}{\rho_{2}}+v^{\prime\,2}_{2}\sin^{2}\alpha^{\prime}_{2}\
,$ $\displaystyle
v^{\prime}_{1}\cos\alpha^{\prime}_{+}=v^{\prime}_{2}\cos\alpha^{\prime}_{2}\
,$ $\displaystyle v^{\prime\,2}_{2}=\frac{5}{3}\frac{p_{2}}{\rho_{2}}\ ,$

which have the following simple solution:

$\displaystyle
v_{p}=v^{\prime}_{0}=\frac{4}{\sqrt{3}}v_{1}=\sqrt{3}\,v_{s1}=\frac{4}{\sqrt{11}}v^{\prime}_{1}=2v^{\prime}_{2}=2v_{2}\
,$ $\displaystyle\rho_{2}=\frac{5}{2}\rho_{1}=10\rho_{0}\ ,\quad\quad
p_{2}=6p_{1}=8\rho_{0}v_{1}^{2}\ ,$
$\displaystyle\frac{u_{2}}{u_{1}}\to\frac{\beta_{2}}{\beta_{1}}=\frac{2}{\sqrt{3}}\
,\quad\quad\frac{u^{\prime}_{2}}{u^{\prime}_{1}}\to\frac{\beta^{\prime}_{2}}{\beta^{\prime}_{1}}=\frac{2}{\sqrt{11}}\
,$ $\displaystyle\frac{u_{12}}{u_{1}}\to\frac{\beta_{12}}{\beta_{1}}=1\
,\quad\quad\frac{u_{p}}{u_{1}}\to\frac{\beta_{p}}{\beta_{1}}=\frac{4}{\sqrt{3}}\
,$
$\displaystyle\tan\alpha^{\prime}_{1}=\tan\alpha_{1}=\tan\alpha^{\prime}_{2}=\tan\alpha_{2}=\frac{1}{\sqrt{2}}\
,$
$\displaystyle\sin\alpha^{\prime}_{1}=\sin\alpha_{1}=\sin\alpha^{\prime}_{2}=\sin\alpha_{2}=\frac{1}{\sqrt{3}}\
,$
$\displaystyle\cos\alpha^{\prime}_{1}=\cos\alpha_{1}=\cos\alpha^{\prime}_{2}=\cos\alpha_{2}=\sqrt{\frac{2}{3}}\
,$ $\displaystyle\tan\chi^{\prime}=\frac{\sqrt{2}}{3}\
,\quad\quad\sin\chi^{\prime}=\sqrt{\frac{2}{11}}\
,\quad\quad\cos\chi^{\prime}=\frac{3}{\sqrt{11}}\ ,$

From Fig. 15 it can be seen that in the relativistic limit
($u_{p}>u_{2}>u_{1}\gg 1$) $u_{12}$ is of order unity, such that while the
first (incident) shock $s1$ is ultra-relativistic (with a relative upstream to
downstream proper speed of $u_{1}\gg 1$), the second (reflected) shock $s2$ is
only mildly relativistic. Therefore, while region 0 is cold, both regions 1
and 2 are relativistically hot, with an adiabatic index of $\hat{\gamma}=4/3$
and $p=e_{\rm int}/3\gg\rho c^{2}$ while $w\approx 4p$ and
$\beta_{c_{s},2}=1/\sqrt{3}$, such that the sonic condition implies
$\beta^{\prime}_{2}=1/\sqrt{3}$, $u^{\prime}_{2}=1/\sqrt{2}$,
$\Gamma^{\prime}_{2}=\sqrt{3/2}$. This also implies
$\frac{p_{2}}{p_{1}}\approx\frac{e_{2,{\rm int}}}{e_{1,{\rm
int}}}\approx\frac{e_{2}}{e_{1}}\approx\frac{w_{2}}{w_{1}}$. Therefore, in the
relativistic limit along the sonic line the equations reduce to:

$\displaystyle\frac{\rho_{1}}{\rho_{0}}=4\Gamma_{1}\
,\quad\frac{4w}{3\rho_{0}c^{2}}\approx\frac{e_{1,{\rm
int}}}{\rho_{0}c^{2}}\approx\frac{e_{1}}{\rho_{0}c^{2}}=4\Gamma_{1}^{2}\ ,$
$\displaystyle\frac{p_{1}}{\rho_{0}c^{2}}=\frac{4}{3}u_{1}^{2}\approx\frac{4}{3}\Gamma_{1}^{2}\
,\quad u_{s1}\approx\sqrt{2}u_{1}\
,\quad\beta_{1,s1}=\frac{\beta_{1}}{3}\approx\frac{1}{3}\ ,$
$\displaystyle\rho_{1}u^{\prime}_{1}\sin\alpha^{\prime}_{+}=\rho_{2}\frac{1}{\sqrt{2}}\sin\alpha^{\prime}_{2}\
,$ $\displaystyle
4p_{1}u^{\prime\,2}_{1}\sin^{2}\alpha^{\prime}_{+}+p_{1}=2p_{2}\sin^{2}\alpha^{\prime}_{2}+p_{2}\
,$ $\displaystyle
4p_{1}\Gamma^{\prime}_{1}u^{\prime}_{1}\sin\alpha^{\prime}_{+}=2\sqrt{3}p_{2}\sin\alpha^{\prime}_{2}\
,$ (21)
$\displaystyle\beta^{\prime}_{1}\cos\alpha^{\prime}_{+}=\frac{1}{\sqrt{3}}\cos\alpha^{\prime}_{2}\
.$

In the relativistic limit $\alpha_{1},\,\alpha_{2}\approx\frac{\pi}{2}$ along
the sonic line, so it is convenient to use the angle
$\bar{\alpha}_{1}=\frac{\pi}{2}-\alpha_{1}$ to express the solution to the
above equations in terms of

$\displaystyle a\equiv\Gamma_{1}\bar{\alpha}_{1}$ $\displaystyle\to$
$\displaystyle\sqrt{\frac{49\sqrt{3}+\sqrt{6195-3576\sqrt{3}-84}}{16(9-5\sqrt{3})}}$
(22) $\displaystyle\approx$ $\displaystyle 0.6004187327198\ ,$

where $a=\Gamma_{1}\bar{\alpha}_{1}\approx u_{1}\bar{\alpha}_{1}\approx
u_{1}\cos\alpha_{1}\approx u_{1}/\tan\alpha_{1}$ approaches a constant values
in this relativistic limit,

$\displaystyle\frac{\Gamma_{p}}{\Gamma_{1}}$ $\displaystyle\to$
$\displaystyle\frac{u_{p}}{u_{1}}\to\sqrt{\frac{2}{1-2a^{2}}}\approx
2.677423238004\ ,$ $\displaystyle\frac{\Gamma_{p}}{\Gamma_{2}}$
$\displaystyle\to$
$\displaystyle\frac{u_{p}}{u_{2}}\to\frac{\sqrt{2}}{\sqrt{3}-1}\approx
1.931851652578\ ,$ $\displaystyle\frac{\Gamma_{2}}{\Gamma_{1}}$
$\displaystyle\to$
$\displaystyle\frac{u_{2}}{u_{1}}\to\frac{\sqrt{3}-1}{\sqrt{1-2a^{2}}}\approx
1.385936251591\ ,$ $\displaystyle\Gamma_{12}$ $\displaystyle\to$
$\displaystyle\frac{3\sqrt{3}-1-4a^{2}}{4\sqrt{1-2a^{2}}}\approx
1.303551928752944\ ,$ $\displaystyle\beta^{\prime}_{1}$ $\displaystyle\to$
$\displaystyle\frac{1}{3}\sqrt{1+16a^{2}}\approx 0.867182056601\ ,$
$\displaystyle u^{\prime}_{1}$ $\displaystyle\to$
$\displaystyle\sqrt{\frac{1+16a^{2}}{8(1-2a^{2})}}\approx 1.741360042442\ ,$
$\displaystyle\Gamma^{\prime}_{1}$ $\displaystyle\to$ $\displaystyle
3/\sqrt{8(1-2a^{2})}\approx 2.008067428503\ ,$
$\displaystyle\sin\alpha^{\prime}_{1}$ $\displaystyle\to$
$\displaystyle\sqrt{1-2a^{2}}\approx 0.528199480119\ ,$
$\displaystyle\cos\alpha^{\prime}_{1}$ $\displaystyle\to$
$\displaystyle\sqrt{2}a\approx 0.849120314915\ ,$
$\displaystyle\tan\alpha^{\prime}_{1}$ $\displaystyle\to$
$\displaystyle\sqrt{1-2a^{2}}/\sqrt{2}a\approx 0.622054932430\ ,$ (23)
$\displaystyle\sin\alpha^{\prime}_{2}$ $\displaystyle\to$
$\displaystyle\frac{1-\sqrt{3}+4a^{2}}{\sqrt{2(2-\sqrt{3})(1+4a^{2})}}\approx
0.620609996819\ ,$ $\displaystyle\cos\alpha^{\prime}_{2}$ $\displaystyle\to$
$\displaystyle\frac{a\sqrt{8(1-2a^{2})}}{\sqrt{2(2-\sqrt{3})(1+4a^{2})}}\approx
0.784119398975\ ,$ $\displaystyle\tan\alpha^{\prime}_{2}$ $\displaystyle\to$
$\displaystyle\frac{1-\sqrt{3}+4a^{2}}{a\sqrt{8(1-2a^{2})}}\approx
0.791473846497\ ,$ $\displaystyle\sin\chi^{\prime}$ $\displaystyle\to$
$\displaystyle\frac{a\sqrt{8(1-2a^{2})}}{\sqrt{1+16a^{2}}}\approx
0.344798730926\ ,$ $\displaystyle\cos\chi^{\prime}$ $\displaystyle\to$
$\displaystyle\frac{1+4a^{2}}{\sqrt{1+16a^{2}}}\approx 0.938676640357\ ,$
$\displaystyle\tan\chi^{\prime}$ $\displaystyle\to$
$\displaystyle\frac{a\sqrt{8(1-2a^{2})}}{1+4a^{2}}\approx 0.367324290498\ ,$
$\displaystyle\frac{\rho_{2}}{\rho_{1}}$ $\displaystyle\to$
$\displaystyle\frac{1-\sqrt{3}+4(4-\sqrt{3})a^{2}}{2(1-\sqrt{3}+4a^{2})\sqrt{1-2a^{2}}}\approx
3.384471042815\ ,$ $\displaystyle\frac{p_{2}}{p_{1}}$ $\displaystyle\to$
$\displaystyle\frac{\sqrt{3}-3+4(4\sqrt{3}-3)a^{2}}{4(1-\sqrt{3}+4a^{2})(1-2a^{2})}\approx
5.549111674227\ .$

## References

  * Adamson & Nicholls (1958) Adamson T. C., Nicholls J. A., 1958, Journal of the Aerospace Sciences, 26, 16
  * Ben-Dor (1987) Ben-Dor G., 1987, Journal of Fluid Mechanics, 181, 467
  * Blandford & McKee (1976) Blandford R. D., McKee C. F., 1976, Physics of Fluids, 19, 1130
  * Chester (1954) Chester W., 1954, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 45, 1293
  * Courant & Friedrichs (1948) Courant R., Friedrichs K. O., 1948, Supersonic flow and shock waves
  * De Rosa et al. (1992) De Rosa M., Famà F., Palleschi V., Singh D. P., Vaselli M., 1992, Phys. Rev. A, 45, 6130
  * Gold (1969) Gold T., 1969, Nature, 221, 25
  * Granot & Rabinovich (2023) Granot J., Rabinovich M., 2023, arXiv e-prints, p. arXiv:2309.11293
  * Gvozdeva & Chulyunin (2015) Gvozdeva L., Chulyunin A. Y., 2015, in Proceedings of the 6th European Conference for Aeronautics and Space Sciences (EUCASS), Krakov.
  * Heilig (1969) Heilig W., 1969, The Physics of Fluids, 12, I
  * Henderson & Gray (1981) Henderson L., Gray P., 1981, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 377, 363
  * Hornung (1986) Hornung H., 1986, Annual review of fluid mechanics, 18, 33
  * Hornung et al. (1979) Hornung H., Oertel H., Sandeman R., 1979, Journal of Fluid Mechanics, 90, 541
  * Hugoniot (1887) Hugoniot H., 1887, Journal de l’École Polytechnique., 57, 3
  * Istomin & Soloviev (2008) Istomin Y., Soloviev F., 2008, arXiv e-prints, p. arXiv:0809.1244
  * Itoh et al. (1981) Itoh S., Okazaki N., Itaya M., 1981, Journal of Fluid Mechanics, 108, 383
  * Lamberts & Daigne (2018) Lamberts A., Daigne F., 2018, MNRAS, 474, 2813
  * Landau & Lifshitz (1987) Landau L. D., Lifshitz E., 1987, Fluid Mechanics: Volume 6. Elsevier
  * Mignone & McKinney (2007) Mignone A., McKinney J. C., 2007, Monthly Notices of the Royal Astronomical Society, 378, 1118
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Norman et al. (1982) Norman M. L., Winkler K. H. A., Smarr L., Smith M. D., 1982, A&A, 113, 285
  * Olim & Dewey (1992) Olim M., Dewey J., 1992, Shock Waves, 2, 167
  * Rankine (1870) Rankine W. J. M., 1870, Phil. Trans. R. Soc.160277–288
  * Rees (1966) Rees M. J., 1966, Nature, 211, 468
  * Rees & Mészáros (1998) Rees M. J., Mészáros P., 1998, ApJ, 496, L1
  * Ryu et al. (2006) Ryu D., Chattopadhyay I., Choi E., 2006, ApJS, 166, 410
  * Salpeter (1964) Salpeter E. E., 1964, ApJ, 140, 796
  * Sari (1997) Sari R., 1997, ApJ, 489, L37
  * Tabak & Rosales (1994) Tabak E. G., Rosales R. R., 1994, Physics of Fluids, 6, 1874
  * Taub (1948) Taub A. H., 1948, Physical Review, 74, 328
  * Thorne (1973) Thorne K. S., 1973, ApJ, 179, 897
  * Thorne & Blandford (2017) Thorne K. S., Blandford R. D., 2017, Modern Classical Physics: Optics, Fluids, Plasmas, Elasticity, Relativity, and Statistical Physics
  * Von Neumann (1963) Von Neumann J., 1963, Theory of games, astrophysics, hydrodynamics and meteorology. Vol. 6, Pergamon Press
  * Wu et al. (2019) Wu K., Zhang G., Kim H. D., 2019, Journal of Visualization, 22, 283


11institutetext: Software Engineering Research Lab (SERL)
Blekinge Institute of Technology
11email<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>22institutetext:
fortiss GmbH, Germany

# Continuous Software Engineering in the Wild

Eriks Klotins 11 0000-0002-1987-2234 Tony Gorschek 1122 0000-0002-3646-235X

###### Abstract

Software is becoming a critical component of most products and organizational
functions. The ability to continuously improve software determines how well
the organization can respond to market opportunities. Continuous software
engineering promises numerous advantages over sprint-based or plan-driven
development. However, implementing a continuous software engineering pipeline
in an existing organization is challenging.

In this invited position paper, we discuss the adoption challenges and argue
for a more systematic methodology to drive the adoption of continuous
engineering. Our discussion is based on ongoing work with several industrial
partners as well as experience reported in both state-of-practice and state-
of-the-art.

We conclude that the adoption of continuous software engineering primarily
requires analysis of the organization, its goals, and constraints. One size
does not fit all purposes, meaning that many of the principles behind
continuous engineering are relevant for most organizations, but the level of
realization and the benefits may still vary. The main hindrances to continuous
flow of software arise from sub-optimal organizational structures and the lack
of alignment. Once those are removed, the organization can implement
automation to further improve the software delivery.

###### Keywords:

Continuous software engineering process improvement continuous integration and
delivery

##  1 Introduction

Software is a critical component of most products, services, manufacturing
processes, and back-office functions. The ability to continuously improve
software is crucial for organizations to respond to market opportunities
swiftly and remain competitive. Software is also becoming increasingly more
complex. Organizations seek to improve both the effectiveness and the
efficiency of software engineering to enable further growth without increasing
overhead and losing flexibility.

Continuous software engineering is a paradigm aiming to streamline software
engineering by delivering software frequently and in small increments, and by
doing so reaping different benefits from fast customer feedback, to continuous
value delivery to said customers. Smaller increments are potentially easier to
plan, develop, integrate, and verify. On the customers’ side, more minor
updates ought to create less disruption and are easier to adopt in contrast to
big-bang software updates requiring downtime and catching up with the new
features. Software vendors can collect more focused telemetry and customer
feedback to steer further product development [2]. More frequent and smaller
software updates enable customers to provide more focused feedback empowering
them to participate in steering the product more actively. In turn, this
allows collaborative and experience-based business models [2].

Continuous integration and delivery (CI/CD) is part of continuous software
engineering. We consider how organizations can apply continuous principles
throughout the whole software engineering process throughout inception,
development, integration, verification, delivery, operation, use, collection
of feedback, and planning the next software iteration steps [1]

However, few organizations have adopted continuous engineering beyond
automating tests and other repetitive development tasks. To the best of our
knowledge, an industrial-scale end-to-end pipeline from inception of a feature
to collecting and analyzing customer feedback has yet to be demonstrated in
peer-reviewed literature [7].

The potential benefits of continuous software engineering have gained much
attention. Benefits like flexibility, efficiency, and improved time-to-market
are appealing to most companies. However, the ability to retrofit an existing
organization, what parameters determine the suitability, and what trade-offs
are associated with adopting continuous software engineering remain largely
unexplored.

This invited position paper discusses the potential, and highlights challenges
hindering the widespread adoption of continuous software engineering. This is
a position paper based on ongoing research with a dozen companies from
multiple domains, featuring different market positions and customers, as part
of the KKS research profile project Software Engineering Rethought (see
http://rethought.se).

Overall, we observe that end-to-end continuous software engineering, as per
many recommendations, may not be applicable everywhere [2, 1]. However,
organizations can adopt parts of the pipeline to streamline software
engineering [4]. The challenge lies in the critical evaluation of the current
situation and the goals of the company to maximize the gains from adopting
continuous software engineering principles.

This paper is structured as follows: Section 2 presents the potential of an
end-to-end continuous software engineering. In section 3, we discuss
challenges and the need for further research, before concluding our paper in
section 4.

##  2 Continuous Software Engineering in a Nutshell

The idea of continuous software engineering originates from lean principles in
manufacturing. One of the key principles in lean is to reduce waste and
maximize customer value by implementing the flow. That is, linking all
relevant production steps together and minimizing the lead time of each step.
In software engineering, this principle can be implemented by delivering small
increments of software [6].

To enable the development and delivery of small software increments, software
vendors need to implement a software delivery pipeline. The pipeline picks up
the latest changes in the source code and automatically performs testing,
integration, delivery, and other required steps to make the latest changes
available to end-users. Once end-users start using the software, feedback and
telemetry is relayed back to the software vendor for analysis and decision
support to steer further product decisions.

In an idealized scenario, end-users gain access to the latest features and
start generating feedback minutes after developers have finished the
development [1, 2].

The state-of-the-art view on the end-to-end pipeline, along with the key
steps, is shown in Fig, 1. The figure shows an idealized scenario that
requires adaptations to fit any real-life scenario.

In the figure, we show the key steps of the pipeline denoted with rectangular
boxes. The arrows represent the flow of software and related artifacts through
the pipeline. Dashed lines represent levels of stakeholders involved in
product development.

Figure 1: An overview of state-of-the-art continuous software engineering
pipeline

We differentiate between four levels of stakeholders, present how continuous
software engineering affects each level, and present several potential
opportunities (PT) associated with adopting continuous engineering. The PT’s
are idealized positive outcomes that can be a result of realizing good
continuous practices, used here to highlight potential benefits only for
illustration purposes.

### Level I: Engineering organization

receives plans from the product planning organization and turns these plans
into working software for deployment (Step 3). The engineering organization
implements build, test, integration, and deployment automation; ensures the
organization and software architecture supports incremental and parallel work.
Note that the exact engineering activities are case-specific. In the figure,
we illustrate the most common activities.

PT 1: Automation and parallel work on small increments are shown to improve
efficiency quality, reduce stress and improve developer satisfaction.
Independent, cross-functional teams taking responsibility for specific
features and modular architectures allow scaling of development organization
with minimal need for additional overhead.

### Level II: Product organization

uses various inputs to devise plans for further product development (Step 2).
These plans guide the engineering organization. Once the engineering is
complete, the engineering organization returns working software (Step 3).

PT 2: Working with lightweight plans and quick turnaround time allows product
organization to rapidly adjust and explore new market opportunities. In
dynamic markets, extensive market research and analysis could be
counterproductive as the results are already outdated when they arrive.
Instead, organizations could adopt a more experiment-driven approach and try
out new ideas. It is possible if testing new ideas and recovering from
unsuccessful experiments are reasonable effort and risk.

### Level III: Operations

takes the working software from the product organization and makes it
available for the end-users (Step 4). The working software generates feedback
and telemetry, enabling monitoring of customer behavior (Step 6) and input for
further planning (Step 1).

PT 3: Automating deliveries and frequently delivering to customers helps
improve time-to-market and reduces release pains. Release pain is associated
with developers required to put in extra effort to make sure the upcoming
release is ready for customers. For large releases, it is associated with
additional stress and workload. However, when releasing small incremental
updates, the stress and workload associated with each release are minimized.

PT 4: Implementing telemetry and automated collection of feedback allows
monitoring the software in use and making adjustments as needed. Quick
turnaround time allows rapid release of patches and bug fixes, thus reducing
the adverse effects of slipped defects. From the customers perspective,
automated, seamless updates remove the need to update the software manually.

### Level IV: End-user/customer organization

receives and uses the software (Step 5). Frequent upgrades and continuous
access to new features enable service and experience-oriented business models.
Such models encourage trust (Step 5.1), and the end-users are incentivised to
actively participate in the product development by providing feedback (Step
5.2).

PT 5:  Both customers and the vendor can benefit from more closer experience
and service based collaboration model. It facilitated trust and created
incentives for customers not to switch to another vendor.

### Cross-cutting concerns

transcends organizational levels. For instance, continuous improvement
(Concern 7.6) attempts to measure and fine-tune the whole cycle continuously.
Continuous experimentation allows a product organization to set up quick
experiments to test market responses to new ideas (Concern 7.5).

PT 6:  Continuous collection of data and frequent execution of the delivery
process allows to systematically and continuously improve the process with
every iteration.

PT 7:  Transparency in planning and road maps of the development company can
achieve a number of benefits. (i) Internal coordination in the development
organization of what is in the pipe short- and medium term. (ii) Customers can
be prepared of what is coming to prepare and plan for changes and benefits.
(iii) Customers and the development organization can rather via transparency
act as partners as interested customers can be active in the planning and
release work.

In summary, adopting an end-to-end continuous software delivery pipeline as
shown in Fig. 1 allows the software vendor to increase internal efficiency and
streamline value creation.

##  3 Challenges and Future Needs

Despite many potential advantages, the utilization of continuous software
engineering remains scarce relative to plan-driven or sprint-based engineering
models. Several of our partner companies have started adopting initiatives
moving towards continuous software engineering. However, the adoption is far
from straightforward and hindered by many challenges. In this section, we
present challenges that, in our view and experience, should be solved to
support an industry-wide adoption of continuous engineering principles.

### Challenge 0: Should you do it?

The tendency to go for one-size-fits-all solutions and the power of “the next
thing everyone is doing” can cause as much damage as benefit. A key is to
start with an analysis of both the capacity and goals of the development
organization, as well as the capacity and needs of the customer’s
organizations.

The development organization needs to specify goals and to break down these
goals to a level where success (or failure) to achieve these goals can be
”measured” by the organization, continuously as changes are realized. This is
paramount to a risk analysis, and also a cost benefit analysis anchored in
reality and objective observations. A significant bonus is that this activity
also enables creating a common understanding of the potential, implications
and direction of the changes from items impacting individual teams, to entire
departments and the whole organization, as well as how they are dependent and
tie together. In this, it is important to separate between the treatment
(change/tool/action/new way of working) and the goals and attainment of the
goals to enable an evidence-based change continuously.

At the same time, knowing your customer and the impact of continuous on them
and the relationship with them is critical. What is the benefit for the
customer? Is the domain and agreements between the customer and the
development company suitable for a continuous environment? One specific
example could be if the customer wants continuous changes in the product, and
is there a telemetrics/data feedback system in place to report on new features
and changes enabling feedback? In some domains and for some development
organizations a continuous model is non-controversial, but in some other
domains most of the benefits attained via try-and-learn-improve are impossible
to achieve. For example, an organization and customer operating in (even
partly) safety-critical products, or where down-time of products incur
substantial consequences, the benefit/risk/cost calculation is significantly
different.

Significant research needs to be conducted, and usable and useful models need
to be developed, to establishing cost-effective ways to continuously evaluate
and course correct the continuous changes needed to achieve a continuous
product development environment; maybe inspired by a hierarchically connected
Goal-Question-Metric model [3] \- back-filled by metrics and data collected
via the continuous feedback cycle. As such, a goal and measurement program is
created (and scaled to be usable and useful for the organization in question)
and it is important to realize that it might lead to the conclusion that
continuous engineering in its entirety and complete idealized form might not
fit your organization. However, there are probably many parts of continuous
engineering that are beneficial for most organizations. Start there.

### Challenge 1: Determining adoption goals and constraints.

We observe that organizations often put forward aims like improving speed and
efficiency to drive the adoption of continuous engineering. However, such aims
are too vague to be measured and drive systematic improvements. When
interpreted by different parts of the organization, vague goals may come at
odds, or even worse result in different interpretations of said goals
resulting in different direction of work and sub-optimizations.

For example, an R&D unit could interpret the efficiency as maximizing the
delivery of new experimental features. For operations, efficiency could mean
minimizing resources to ensure services availability. Without a joint view of
what efficiency means and how it is measured in the given organizational
context, internal deadlocks may arise, hindering the company’s adoption effort
and optimal operation.

Attaining specific goals often imply trade-offs. As in the earlier example, it
could be challenging to launch many innovative features and reach high
stability of services simultaneously. Such trade-offs need to be identified
and analyzed to understand the associated constraints and degrees of freedom.

The analysis of organizational goals, constraints, and trade-offs should drive
the organizational change towards continuous software engineering. One
powerful tool in coordinating an organization, and making goals clear is to
break-down goals into measurable effects or metrics. This requires a number of
steps. The terms ”value”, ”efficiency”, or ”effectiveness” need to be defined
for each part of the organization [3]. These definitions need to be
coordinated and streamlined - in essence shared. Then, metrics on how to
ascertain level of success have to be detailed very early. How to measure if
you are improving towards a goal is a prerequisite before adopting any
treatment or change. For example, if you say a specific practice or set of
practices should be realized to improve customer value, what type of customer
value, and how do you measure it? This analysis would also pinpoint
bottlenecks, inefficiencies, and areas of improvement, in addition to acting
as a coordinating force as an organization realizes changes. You need to be
able to measure benefits of a change as well as you measure the cost of said
change.

### Challenge 2: Considering the return-of-investment perspective

Retrofitting an organization with a new continuous engineering pipeline and
new ways of working is a substantial investment. The investments should be
justified with potential benefits and be aligned with organizational
objectives (Challenge 1). Notably, the organization must be prepared to
realize the potential and materialize the benefits.

For example, an organization may invest in data collection (PT 5) and gain the
potential of improved data-driven decisions (PT 2). However, if the rest of
the organization is not ready to use the data in decision support, the
potential is not realized, and the investment is wasted. Furthermore,
customers may be slow in adopting new features, thus delaying the feedback and
nullifying its value to the organization.

The organization should perform a cost-benefit analysis to gauge the viability
of any goals, new practices, and working methods. Adoption champions should do
such analysis in parallel with determining adoption goals, trade-offs, and
constraints, see Challenge 1.

To address Challenges 1-2, we propose a model supporting inventory of goals,
trade-offs, and return-of-investment calculation to support the systematic
adoption of continuous practices. As a part of this analysis, the break-down
of goals into defining the terms used to establish a common vocabulary, then
establishing how to ”measure” goal attainment by breaking down goals into
measurable items. This also allows for course corrections during the changes
associated with retrofitting.

### Challenge 3: Focus on cost savings instead of value and potential
creation.

Our industry partners often mention the need to reduce cost and improve
efficiency (PT 1) as reasons for considering continuous software engineering.
Potential benefits like new value streams, improved time-to-market, and new
business models are rarely, if at all, mentioned.

The question if adopting continuous practices would provide some cost savings
is flawed. We compare such a question to asking whether eating healthier will
be cheaper.Both eating healthier and adopting continuous engineering will
likely cost more. Developing tests, maintaining test suites and automation
infrastructure, collecting and maintaining test data, refactoring software
architectures, and driving organizational adjustments will push the overall
cost of software engineering upwards. At the same time, speed and flexibility
will open up opportunities for new offerings enabling the organization to
become faster at responding to market opportunities, among other benefits.
Eating healthier may cost more, but may enable you a better and longer life.

We propose emphasizing the business value arising from streamlining software
delivery over the exclusive focus on potential cost savings. More cross-
disciplinary studies are needed to explore the business value perspective. We
observe that continuous shares many characteristics with so called ”digital
transformation”. Considering how continuous engineering fits into the broader
organizational transformation could help to shift the focus from cost savings
to unlocking the organizational potential [5].

### Challenge 4: Remembering Conway’s law.

In 1967, M. Conway formulated an adage that organizations develop systems that
mimic their communication structure. In software engineering terms, software
architecture should follow organizational structures. To change the prior, one
needs to change the latter first.

Best practices of continuous software engineering dictate that development
work should be done by real cross-functional teams that own the development of
whole features. A team should have full responsibility from generating
improvement ideas to development, testing, delivery, and telemetry analysis.
This practice allows maintaining a modular architecture and to minimize the
gap between organizational structures and software architectures.

From discussions with our partner companies, we learned that organizations
often find their software architectures monolithic and poorly suited for
parallel development, automation, and modular deployments. The software
travels through various organizational silos, each performing a specific
function without the complete picture. Handover from one silo to another
creates friction and bottlenecks.

Attempts to break the software monolith often lead to failure as the
surrounding organization remains the same. Changing the organizational
structures is often extremely challenging due to internal inertia and
resistance to changes.

We propose to explore the adoption of continuous engineering practices from
the organizational view first. That is, analyze what organization silos and
bottlenecks currently limit the continuous flow of software. Once
organizational inefficiencies are addressed, the software delivery process can
be further improved with tools such as automation. In addition, the detailed
goals and following follow-up (metrics) should be owned by the parts of the
organization (the teams) that own the module or part of the architecture. This
way you can measure the fit of the organization to the goals of development.
For example, a change in organizational structure might result in (metric)
less waiting times between teams.

### Challenge 5: Ever-increasing complexity

Contemporary market-driven software engineering exists in a dynamic
environment. It faces mercurial market influences, changing organizational
goals, technologies, and the growing size and complexity of software and the
surrounding organization.

A single person, or a small group of people, can no longer grasp the
complexity and make optimal and timely decisions using their expertise alone.
This increase in complexity has far-reaching implications for how
organizations make decisions.

One viable way forward could be to consider data collection and analysis as an
integral part of the product and the engineering process. The organization can
use data to support decisions both on what features to develop (What to
build?) and how to improve the engineering process (How to build?).

Furthermore, retrospective analysis, that is, analyzing past events, has
limited use in an increasingly changing environment. There is a potential to
explore the applicability of inferential statistical methods, simulations, and
machine learning techniques to make relevant predictions. This is however
subject to significant research initiative to develop transparent tools and
methods that can be trusted by practitioners to predict and simulate
developments, at least in the short-term.

##  4 Conclusions

This invited position paper analyzes continuous software engineering and
pinpoints several challenges of adopting ongoing software engineering from our
collaboration with industrial partners. The challenges emphasize the need for
systematic methods to analyze organizational goals, context, structures,
constraints, among other contextual factors, to remove inefficiencies and
realize the full potential of software-intensive products and services.

We wish to highlight that most inefficiencies and obstacles to streamlined
value delivery can be traced to the lack of organizational alignment and
coordination. The spirit of continuous engineering is to identify and remove
such hindrances systematically. In the end, continuous software engineering is
more about being pragmatic and disciplined in engineering than it is about
automation.

## References

  * [1] Fitzgerald, B., Stol, K.J.: Continuous software engineering: A roadmap and agenda. Journal of Systems and Software 123, 176–189 (2017)
  * [2] Humble, J., Kim, G.: Accelerate: The science of lean software and devops: Building and scaling high performing technology organizations. IT Revolution (2018)
  * [3] Khurum, M., Gorschek, T., Wilson, M.: The software value map—an exhaustive collection of value aspects for the development of software intensive products. Journal of software: Evolution and Process 25(7), 711–741 (2013)
  * [4] Klotins, E., Gorschek, T.: Towards cost-benefit evaluation for continuous* software engineering activities. In review by Empirical Software Engineering Journal xx(x), xx–xx (2021)
  * [5] Klotins, Eriks, P.A.E.: The unified perspective of digital transformation and continuous software engineering. In: Proceedings of 5th International Workshop on Software-intensive Business (IWSIB) (2022)
  * [6] Poppendieck, M., et al.: Principles of lean thinking. IT Management Select 18(2011), 1–7 (2011)
  * [7] Shahin, M., Babar, M.A., Zhu, L.: Continuous integration, delivery and deployment: a systematic review on approaches, tools, challenges and practices. IEEE Access 5, 3909–3943 (2017)